### Set-up and sample preparation

Our measurements were performed under ultrahigh vacuum (base pressure, *p* < 10^{−10} mbar) with a home-built conductive-tip atomic force microscope equipped with a qPlus sensor^{34} (resonance frequency, *f*_{0} = 30.0 kHz; spring constant, *k* ≈ 1.8 kNm^{−1}; quality factor, *Q* *≈* 1.9 × 10^{4} and 2.8 × 10^{4}) and a conductive Pt-Ir tip. The microscope was operated in frequency-modulation mode, in which the frequency shift Δ*f* of the cantilever resonance is measured. The cantilever amplitude was 0.55 Å (1.1 Å peak to peak), except if specified otherwise. Constant-height AFM images were taken at tip-height changes Δ*z* with respect to the set point, as indicated. Positive Δ*z* values indicate being further away from the surface.

As a sample substrate, we used a cleaved mica disc, on which we deposited gold in a loop structure (diameter, *d* = 10.5 mm; thickness, *t* = 300 nm) by means of electron-beam physical vapour deposition. This gold structure contained a 100-μm-wide constriction, on which the measurements were performed. A non-conducting spacer material was introduced below the mica disc to prevent eddy-current screening of the RF magnetic field. The sample was prepared by short sputtering and annealing cycles (annealing temperature, *T* ≈ 550 °C) to obtain a clean Au(111) surface. On half of the sample, a thick NaCl film (>20 monolayers) was grown at a sample temperature of approximately 50 °C; the other half of the sample was used for tip preparation, presumably resulting in the tip apex being covered with gold. Part of the data was measured with a CO-functionalized tip apex. To this end, a sub-monolayer coverage of NaCl was also deposited on the whole surface at a sample temperature of approximately 35 °C, to grow two monolayer NaCl islands also on the half of the sample used for tip preparation. After preparing a tip by indenting into the remaining gold surface, a CO molecule was picked up from the two monolayer NaCl islands, after which the tip was transferred to the thick NaCl film^{35}. The NaCl film inhibits any electrons to tunnel to or from the gold structure. The voltage that is applied to the gold structure with respect to the tip represents a gate voltage (*V*_{G}), gating the molecular electronic states against the chemical potential of the conductive tip. The measured molecules (pentacene-h_{14} and PTCDA-h_{8}, Sigma-Aldrich; pentacene-d_{14}, Toronto Research Chemicals) and CO for tip functionalization were deposited in situ onto the sample inside the scan head at a temperature of approximately 8 K. Pentacene was reported to adsorb centred above a Cl^{−} anion with the long molecular axis aligned with the polar direction of NaCl, resulting in two equivalent azimuthal orientations^{36}.

The AC voltage pulses were generated by an arbitrary waveform generator (TGA12104, Aim-TTi), combined with the DC voltage, fed to the microscope head by a semi-rigid coaxial high-frequency cable (Coax Japan Co. Ltd.) and applied to the gold structure as *V*_{G}. The high-frequency components of the pulses of *V*_{G} lead to spikes in the AFM signal because of the capacitive coupling between the sample and the sensor. To suppress these spikes, we applied the same pulses with opposite polarity and adjustable magnitude to an electrode that capacitively couples to the sensor.

The RF signal was produced by a software-defined radio (bladeRF 2.0 micro xA4, Nuand), low-pass filtered to eliminate higher-frequency components and amplified in two steps (ZX60-P103LN+, Mini Circuits; KU PA BB 005250-2 A, Kuhne electronic). The RF was pulsed using RF switches (HMC190BMS8, Analog Devices), which were triggered by the arbitrary waveform generator, allowing synchronization with *V*_{G} and control over the pulse duration. The pulsed RF signal was fed into the microscope head by a semi-rigid coaxial high-frequency cable (Coax Japan Co. Ltd.) ending in a loop, inductively coupling the RF signal to the gold loop on the sample. These two loops are in the surface plane of the sample, such that the inductive coupling adds a vertical *z* component to the magnetic field. The field generated by the microstrip is associated to field lines looping around the microstrip (Ampère’s law). At the position of a molecule placed above the strip, the local magnetic field resulting from the microstrip is expected to be homogeneous, in the surface plane and perpendicular to the direction of the microstrip.

The RF signal transmission of the cables including the loop for inductive coupling was detected by a magnetic field probe and can be well approximated to be constant over intervals of tens of megahertz around the T_{X}–T_{Z} transition, that is, wider than the spectral features observed in the experiments. Although the microstrip will contribute to the overall transmission of the signal to the local magnetic field, it is expected to not introduce any resonances in the frequency range of interest. Note that the RF signal at a frequency of 1,500 MHz has a wavelength roughly three times the entire circumference of the loop of the microstrip.

To excite the entire broadened ESR resonance for the lifetime measurements of the triplet state T_{1} with RF, shown in Fig. 1b, we used IQ modulation to generate a broadband RF pulse. We created a chirped pulse with a width of 12 MHz, a repetition time of 5 μs and a centre frequency of 1,544 MHz. Thereby, the RF signal spans the range 1,538–1,550 MHz in frequency space.

### ESR-AFM pulse sequence and data acquisition

The description of the measurement of the triplet-state lifetime can be found in ref. ^{15}. The ESR-AFM experiments were performed with a similar voltage-pulse sequence, which is shown in Extended Data Fig. 1. Between each individual voltage-pulse sequence, the voltage is set to *V*_{deg}, the bias voltage, at which the respective ground states of the positively charged (D_{0}) and the neutral (S_{0}) molecules are degenerate. This way, the spin states of the molecule are converted to different charge states and detected^{15} by charge-resolving AFM (ref. ^{17}). The dwell voltage pulse duration *t*_{D} was fixed to 100.2 μs for pentacene and, simultaneously, an RF pulse with a variable frequency was applied. In the case of PTCDA, the triplet lifetimes were determined to be 350 ± 43 µs, 170 ± 13 µs and 671 ± 62 µs, and *t*_{D} was set to 501 µs for the ESR-AFM spectra shown. To reduce the statistical uncertainty for a given data-acquisition time, we repeated the pump–probe pulse sequence 160 or 320 times per second (instead of eight times per second for the lifetime measurements of the triplet state T_{1}). Note that, to prevent the excitation of the cantilever, the durations of the voltage pulses were set to an integer multiple of the cantilever period (33.4 μs). At this high repetition rate of the voltage-pulse sequence, the charge states cannot be read out individually. Instead, the AFM signal, that is, the frequency shift Δ*f*, was averaged over an interval of 20 s. This average frequency shift ⟨Δ*f*⟩ reflects the ratio of the charged and neutral states and thus the triplet and singlet states, but because the change in Δ*f* is very small, it is also sensitive to minor fluctuations in the tip–sample distance (see Extended Data Figs. 2c and 5b). To minimize the fluctuations in tip–sample distance, the tip–sample distance was reset by shortly turning on the Δ*f* feedback either after every sweep of the RF or after a fixed time (15 to 60 min). To minimize the dependence of ⟨Δ*f*⟩ on the remaining fluctuations in tip height, ⟨Δ*f*⟩ was normalized using the frequency shifts of the charged Δ*f*^{+} and neutral Δ*f*^{0} molecule, as \({\Delta f}_{{\rm{norm}}}=\frac{\langle \Delta f\rangle -{\Delta f}^{0}}{{\Delta f}^{+}-\Delta {f}^{0}}\). These frequency shifts were determined at the beginning and end of every 20-s data trace (see Extended Data Fig. 2a); the charge state was changed by applying small voltage pulses (\({V}_{{\rm{set}}}^{0}={V}_{{\rm{\deg }}}+0.3\,{\rm{V}}\), \({V}_{{\rm{set}}}^{+}={V}_{{\rm{\deg }}}-0.3\,{\rm{V}}\)). Tunnelling events during the readout of these frequency shifts were minimized by using a tip–sample distance at which the decay constant for the decay of the D_{0} state into the S_{0} state during a pulse of *V*_{deg} + 1.2 V was around 4 μs (note that this requirement restricts the possible distances to a small range, as the tip height should also be small enough such that the tunnelling processes are considerably faster than the triplet decay; see also discussion of the spatial resolution in the main text).

If still a charging event happened, the data trace was discarded. To maximize the rate of the tunnelling processes during the voltage-pulse sequence, the beginning and end of the voltage pulses were synchronized with the closest turnaround point of the cantilever movement. The data-acquisition and renormalization scheme to derive Δ*f*_{norm} is shown in Extended Data Fig. 2. As can be seen in Extended Data Figs. 2c and 5b, also the raw ⟨Δ*f*⟩ signal exhibits the ESR features but with stronger baseline drift.

Note that Δ*f*_{norm} typically deviates from the triplet population, but that, for a given measurement, a linear relation between them exists. This deviation arises from the voltage pulses that are for pentacene turned on for 4.3% of the time, during which the frequency shift corresponds to the applied voltages and thus crucially depends on the exact shape of the Kelvin probe force parabola^{17}. This explains the differences in the baseline of the Δ*f*_{norm} signal (without RF or RF off-resonance) for different measurements—even for those above the same molecule—owing to differences in the position above the molecule. Quantitative results (right axis of Fig. 1c,d) can be obtained from a calibration measurement in which the population was determined by counting the individual outcomes after each pulse sequence at a repetition rate of eight per second. This calibration was performed for an RF corresponding to the maximum of the ESR signal, as well as an RF that was off-resonance. For both cases, 7,680 pump–probe cycles were recorded.

Note that we do not observe any appreciable change in the damping of the cantilever during an RF sweep.

### Experimental uncertainties and statistical information

To determine the uncertainty on the ESR-AFM data points, the 20-s data traces were repeated several times and the error bars were extracted as the s.d. of the mean of these repetitions. This way, any type of non-systematic uncertainty will be accounted for, irrespective of its source (see next paragraph). Note that the hydrogen spins can have a different configuration for every individual readout^{19}. Given our large number of sampling events, we acquire an average over the possible nuclear spin configurations.

The three main sources of uncertainty of Δ*f*_{norm} are the statistical uncertainty from the finite number of repeats^{15}, the remaining drift of the tip height and the noise on the frequency shift Δ*f*. We choose the number of repeats per data point such that the statistical uncertainty becomes comparable with the other two sources of uncertainties; depending on the exact experimental conditions, any of these three sources can dominate. Probe tips that give a strong response to charging (a large charging step in the Kelvin parabola) provide a better signal-to-noise ratio and, therefore, a smaller relative uncertainty. To minimize this contribution to the uncertainty, we only used tips for which the charging step was large compared with the noise in Δ*f* (size of charging step 0.2–0.4 Hz for tips with a Δ*f* setpoint around −1.5 Hz at zero bias; Δ*f* averaged over 1 s exhibits a typical uncertainty of 1 mHz). In the case of the data in Fig. 4c, top and Extended Data Fig. 8c, bottom, the drift was clearly dominating the error margins. Therefore, the noise resulting from drift was, for these datasets, further minimized by setting the average of every repeat equal to the average over all data points.

In the case of pentacene, the T_{X}–T_{Z} transition was measured for 19 individual pentacene-h_{14} molecules, 20 pentacene-d_{14} and 1 pentacene-h-d_{13}; for 18 of these molecules, we also measured the T_{X}–T_{Y} transition. The T_{Y}–T_{Z} transition of PTCDA was measured in 20 individual spectra for 2 molecules, whereas T_{X}–T_{Y} and T_{X}–T_{Z} transitions (not shown) were measured 3 and 5 times, respectively. In total, 16 different tips were used for these measurements. Molecule-to-molecule variations of the resonance frequencies for two different tips are shown in Extended Data Table 1.

### Spin Hamiltonian and eigenstates

The pulse sequence prepares the molecule into an electronically excited state with one unpaired electron in both the HOMO and the LUMO. These electrons couple through exchange interaction, leading to a large energy difference of Δ*E* ≈ 1.1 eV (refs. ^{37,38}) between the excited singlet S_{1} and the excited triplet T_{1} states. We note in passing that this energy difference allows us to selectively occupy T_{1} instead of S_{1}. With respect to exchange interaction, all three triplet sub-states of T_{1} are degenerate. These are typically represented in the basis of magnetic quantum numbers *m*_{S} = −1, 0 and +1, which—in the representation of the two coupled spins—is T_{−1} = |↓↓⟩, T_{0} = (|↑↓⟩ + |↓↑⟩)/√2 and T_{+1} = |↑↑⟩, respectively.

As explained in the following, the magnetic dipole–dipole interaction between the two electron spins, which is orders of magnitude weaker than the exchange interaction, lifts this degeneracy, leading to a splitting of the triplet states, called the zero-field splitting. We note that the zero-field splitting may also have contributions from spin–orbit interaction. The dipole–dipole interaction is described by the Hamiltonian^{39}

$${\mathscr{H}}=-\frac{{\mu }_{0}{\gamma }_{{\rm{e}}}^{2}}{4\pi {r}^{3}}(3({{\bf{S}}}_{1}\cdot \widehat{{\bf{r}}})({{\bf{S}}}_{2}\cdot \widehat{{\bf{r}}})-{{\bf{S}}}_{1}\cdot {{\bf{S}}}_{2}){\hbar }^{2},$$

with the two spins **S**_{1} and **S**_{2} at a distance *r* in a relative direction \(\widehat{{\bf{r}}}={\bf{r}}/r\). *μ*_{0} is the magnetic constant, *γ*_{e} the gyromagnetic ratio, ħ the reduced Planck constant and **r** the vector connecting the two spins. Notably, the magnetic dipole–dipole interaction is highly anisotropic, that is, for given spin orientations, it strongly differs and even changes sign for different relative positions of the two spins (see Extended Data Fig. 3a). The spatial positions of the electron spins are given by the orbital densities of the two electrons, the confinement of which is very different along the three molecular axes (see Extended Data Fig. 3b). Note that, for pentacene and PTCDA, the *z* direction is perpendicular to the molecular plane and, thereby, perpendicular to the surface plane; *x* points along the long molecular axis^{21}. The anisotropy of the dipole–dipole interaction together with the non-uniformity of orbital densities gives rise to an energy difference in the range of microelectronvolts for the spins pointing in different real-space dimensions. This zero-field splitting is thus a fingerprint of the orbital densities and thereby the molecular species (see Fig. 2).

The corresponding eigenstates are no longer T_{−1}, T_{0} and T_{+1} but T_{X}, T_{Y} and T_{Z.} The latter eigenstates expressed in the basis of the former read T_{X} = (T_{−1} − T_{+1})/√2, T_{Y} = (T_{−1} + T_{+1})i/√2 and T_{Z} = T_{0}, whereas expressed as the states of the two individual spins |m_{s1} m_{s2}〉, they are T_{X} = (|↓↓⟩ − |↑↑⟩)/√2, T_{Y} = (|↓↓⟩ + |↑↑⟩)i/√2 and T_{Z} = (|↑↓⟩ + |↓↑⟩)/√2. Further, they have the property that the expectation value of the total spin ⟨T_{i}|**S**|T_{i}⟩ vanishes for all three states T_{i=X,Y,Z}, whereas \(\langle {{\rm{T}}}_{i}| {{\bf{S}}}_{j}^{2}| {{\rm{T}}}_{i}\rangle =\left(1-{\delta }_{ij}\right)\). Here *δ*_{ij} is 0 for *i* ≠ *j* and 1 for *i* = *j*. ⟨**S**⟩ = 0 renders these triplet states relatively insensitive to external perturbations; an external magnetic field affects the system and energies only to the second order.

The spin Hamiltonian \({\mathscr{H}}\) for the zero-field splitting and an external magnetic field **B** (excluding hyperfine terms) is \({\mathscr{H}}={\bf{S}}\widehat{D}{\bf{S}}+{g}_{{\rm{e}}}{\mu }_{{\rm{B}}}{\bf{SB}}\), with the dipole–dipole-interaction tensor \(\widehat{D}\). Explicitly expressed in the basis of the zero-field split states T_{X}, T_{Y} and T_{Z}, it reads^{39}

$${\mathscr{H}}=\left[\begin{array}{ccc}{{\epsilon }}_{{\rm{X}}} & {-{\rm{i}}g}_{{\rm{e}}}\,{\mu }_{{\rm{B}}}{B}_{{\rm{Z}}} & {{\rm{i}}g}_{{\rm{e}}}\,{\mu }_{{\rm{B}}}{B}_{{\rm{Y}}}\\ {\rm{i}}{g}_{{\rm{e}}}\,{\mu }_{{\rm{B}}}{B}_{{\rm{Z}}} & {{\epsilon }}_{{\rm{Y}}} & {-{\rm{i}}g}_{{\rm{e}}}\,{\mu }_{{\rm{B}}}{B}_{{\rm{X}}}\\ -{\rm{i}}{g}_{{\rm{e}}}\,{\mu }_{{\rm{B}}}{B}_{{\rm{Y}}} & {\rm{i}}{g}_{{\rm{e}}}\,{\mu }_{{\rm{B}}}{B}_{{\rm{X}}} & {{\epsilon }}_{{\rm{Z}}}\end{array}\right]$$

Here *μ*_{B} is the Bohr magneton, *g*_{e} is the electron *g*-factor and *ϵ*_{X}, *ϵ*_{Y} and *ϵ*_{Z} are the zero-field energies of T_{X}, T_{Y} and T_{Z}, respectively. With increasing external magnetic field, the eigenstates will gradually change and asymptotically become the states T_{−1}, T_{0} and T_{+1} in the limit of large magnetic fields (for example, see Extended Data Fig. 3c).

### Selection rules

It follows from the above Hamiltonian that any two of the three zero-field split states are coupled by means of the magnetic field component pointing in the remaining third real-space dimension^{32}. For example, T_{X} and T_{Z} are only coupled through *B*_{Y}, such that only the latter can drive the T_{X}–T_{Z} transition. Because *x*, *y* and *z* are defined with respect to the molecular axes, they might not coincide for different individual molecules (for example, see Fig. 4c).

### Origin of asymmetric lineshape

The hyperfine interaction in protonated pentacene can be described as an effective magnetic field *B*_{HFI} created by the nuclei with non-zero spins acting on the electron spins. Assuming a random orientation of the 14 proton nuclear spins at a given point in time, *B*_{HFI} will fluctuate around zero-field (see Extended Data Fig. 3c) and point in a random direction. Because of the many fluctuating nuclear spins acting together at random, the probability distribution of *B*_{HFI} has its maximum around zero and falls off towards larger absolute values. The influence of *B*_{HFI} is always small compared with the zero-field splitting such that *B*_{HFI} shifts the energies of the triplet states only to the second order, that is, to \(\propto {B}_{{\rm{HFI}}}^{2}\) (ref. ^{21}). This is depicted in Extended Data Fig. 3c for the *B*_{HFI} component in the *z* direction, *B*_{HFI,Z}, in which it becomes clear that the broadening is single-sided in this case. The *x* and *y* components of *B*_{HFI} contribute much less to the broadening. From Extended Data Fig. 3c, it becomes clear that the curvature around *B*_{HFI} = 0 of the hyperbolic avoided crossing is responsible for the asymmetric broadening. This curvature is inversely proportional to the energy difference of the respective pair of states. As the T_{X}–T_{Y} transition has the smallest energy splitting of all possible pairs, the broadening is dominated by their avoided crossing occurring along the *z* component of *B*_{HFI} (ref. ^{21}). Specifically, for the case of pentacene, this effect is smaller by roughly one order of magnitude in the other two directions.

Different individual isotopes contribute differently to *B*_{HFI}, such that different isotopologues give rise to a different probability distribution of *B*_{HFI} when considering all possible nuclear spin configurations. The assignment to the isotopologue is done in comparison with previous work^{29}, based on the line profile, which not only includes the width but also its shape. Because the hyperfine interaction enters as a second-order term, the mere presence of one nucleus (for example, a proton) with strong hyperfine interaction also influences how strongly all the other nuclei (for example, deuterons) affect the line, thereby changing its overall shape.

Analogously, the hyperfine interaction of the eight proton nuclear spins in PTCDA gives rise to its asymmetric lineshape (see Fig. 2b). Note that the asymmetric shoulder appears here at the low-frequency side. Such a lineshape is expected for the T_{Y}–T_{Z} signal^{21}, as becomes clear from Extended Data Fig. 3c (when considering the T_{Y}–T_{Z} transition instead of the T_{X}–T_{Z} transition that is explicitly illustrated). The T_{Y}–T_{Z} signal is the largest for PTCDA; small signals were also observed for the T_{X}–T_{Y} transition (at 252 MHz) and the T_{X}–T_{Z} transition (at 1,501 MHz). Note that this is in contrast to pentacene, for which the T_{X}–T_{Z} signal is the largest and the T_{Y}–T_{Z} transition was not detected (because of the similarity in the lifetimes of its T_{Y} and T_{Z} states).

### Fitting of the lineshapes

As explained in the previous section, the hyperfine interaction is the origin of the asymmetric lineshape of the ESR signals. The lineshape of the T_{X}–T_{Z} transition can be well approximated by a sudden onset at the frequency *f*_{onset} followed by an exponential decay of width *f*_{decay} (ref. ^{20}) as

$$\Theta (f-{f}_{{\rm{o}}{\rm{n}}{\rm{s}}{\rm{e}}{\rm{t}}}\,)\exp (\,-\,(f-{f}_{{\rm{o}}{\rm{n}}{\rm{s}}{\rm{e}}{\rm{t}}}\,)/{f}_{{\rm{d}}{\rm{e}}{\rm{c}}{\rm{a}}{\rm{y}}}\,),$$

in which Θ(*x*) denotes the Heaviside function.

A second contribution to the overall lineshape results from the finite lifetimes of the involved states. This leads to a lifetime broadening, resulting in a Lorentzian of the form

$${\pi }^{-1}\Gamma /({(f-{f}_{{\rm{res}}})}^{2}+{\Gamma }^{2})$$

centred around each resonance frequency *f*_{res} with a full width at half maximum of 2Γ. Accordingly, the experimental resonances are fit to a convolution of the above two functions, allowing to extract the broadening owing to the hyperfine interaction and the finite lifetimes separately. We note that non-Markovian processes^{20,24} may lead to a deviation from the idealized Lorentzian and that power broadening was avoided in the measurements of the ESR-AFM signals. The effect of power broadening is illustrated in Extended Data Fig. 4.

### Rabi oscillations simulations and delay time

The Rabi oscillations were measured using an RF pulse applied around the middle of the dwell voltage pulse with a varying duration and a frequency corresponding to the maximum of the T_{X}–T_{Z} or T_{X}–T_{Y} ESR signal. To illustrate the effect of such an RF pulse for the T_{X}–T_{Z} transition, the evolution of the populations of the three triplet states and the singlet state during the dwell voltage pulse were simulated, as shown in Extended Data Fig. 6.

These simulations were performed using the Maxwell–Bloch equations^{40}, analogous to the model used for ODMR^{41}. The Rabi oscillation data are a temporal average of a single molecule, which—according to the ergodic assumption—is the same as an ensemble average. Therefore, we can use the density-matrix formalism^{42} to simulate our data. Note that, on driving the T_{X}–T_{Z} transition, T_{Y} is decoupled from the T_{X} and T_{Z} dynamics and simply decays independently. The Bloch equations in the density-matrix formalism can therefore be restricted to the two coupled states^{43}, here T_{X} and T_{Z}, whereas the occupation of the third triplet state is treated separately as a simple exponential decay function. With respect to the two coupled states, the system is described by the density matrix^{42}

$$\rho =\left[\begin{array}{cc}{\rho }_{{\rm{ZZ}}} & {\rho }_{{\rm{ZX}}}\\ {\rho }_{{\rm{XZ}}} & {\rho }_{{\rm{XX}}}\end{array}\right]$$

and evolves according to the Liouville equation^{42}

$$\frac{{\rm{d}}\rho }{{\rm{d}}t}=-\frac{i}{\hbar }\left[{\mathscr{H}},\rho \right].$$

The Hamiltonian of the molecule interacting with the RF field (with Rabi rate Ω) at resonance with the T_{X}–T_{Z} transition (with resonance frequency *ω*_{Z} − *ω*_{X}) can be written as^{41}

$${\mathscr{H}}=\left[\begin{array}{cc}\hbar {\omega }_{{\rm{X}}} & -\hbar \Omega \cos \left({{\omega }_{{\rm{Z}}}-\omega }_{{\rm{X}}}\right)\\ -\hbar \Omega \cos \left({{\omega }_{{\rm{Z}}}-\omega }_{{\rm{X}}}\right) & \hbar {\omega }_{{\rm{Z}}}\end{array}\right]$$

The time evolution of the density operator in the rotating-frame approximation, with phenomenologically added relaxation and dephasing terms, can be described as^{41}

$$\frac{{\rm{d}}{\rho }_{{\rm{XX}}}}{{\rm{d}}t}=\frac{i\Omega }{2}\left({\rho }_{{\rm{ZX}}}-{\rho }_{{\rm{XZ}}}\right)-\frac{{\rho }_{{\rm{XX}}}}{{\tau }_{{\rm{X}}}}$$

$$\frac{{\rm{d}}{\rho }_{{\rm{ZZ}}}}{{\rm{d}}t}=\frac{i\Omega }{2}\left({\rho }_{{\rm{XZ}}}-{\rho }_{{\rm{ZX}}}\right)-\frac{{\rho }_{{\rm{ZZ}}}}{{\tau }_{{\rm{Z}}}}$$

$$\frac{{\rm{d}}{\rho }_{{\rm{XZ}}}}{{\rm{d}}t}={-\rho }_{{\rm{XZ}}}\left(\frac{1}{{T}_{2}}+\frac{1}{2{\tau }_{{\rm{X}}}}+\frac{1}{2{\tau }_{{\rm{Z}}}}\right)+\frac{i\Omega }{2}\left({\rho }_{{\rm{ZZ}}}-{\rho }_{{\rm{XX}}}\right)$$

$$\frac{{\rm{d}}{\rho }_{{\rm{ZX}}}}{{\rm{d}}t}={-\rho }_{{\rm{ZX}}}\left(\frac{1}{{T}_{2}}+\frac{1}{2{\tau }_{{\rm{X}}}}+\frac{1}{2{\tau }_{{\rm{Z}}}}\right)+\frac{i\Omega }{2}\left({\rho }_{{\rm{XX}}}-{\rho }_{{\rm{ZZ}}}\right)$$

The time evolution of T_{Y} is simply given by

$$\frac{{\rm{d}}{\rho }_{{\rm{YY}}}}{{\rm{d}}t}=-\frac{{\rho }_{{\rm{YY}}}}{{\tau }_{{\rm{Y}}}}$$

These last five equations were used for the simulation for Extended Data Fig. 6. As input parameters for the simulation, we used the parameters that were experimentally derived for pentacene-h_{14}: the decay constants of the triplet states: *τ*_{X} = 20.8 μs, *τ*_{Y} = 66.6 μs and *τ*_{Z} = 136.1 μs; the decay constant of the Rabi oscillations: *T*_{2} = 2.2 μs; the initial populations *ρ*_{XX} = *ρ*_{YY} = *ρ*_{ZZ} = 0.8/3 and coherences *ρ*_{XZ} = *ρ*_{ZX} = 0; the starting time of the RF pulse *t*_{S} = 45.1 μs and its duration of 4 and 4.5 Rabi-oscillation periods, respectively. Note that interconversion between T_{X} and T_{Z} resulting from spin-lattice relaxation is assumed to be negligible compared with *τ*_{X} and *τ*_{Z}.

Here the initial occupation of the T_{X}, T_{Y} and T_{Z} states are assumed to be all equal to 0.8/3. Simulations and data in ref. ^{15} show that the triplet state is initially approximately 80% occupied (this value depends on the exact tip position, as the two competing tunnelling rates to form the T_{1} and S_{0} states depend on the wave-function overlap between tip and the LUMO and the HOMO, respectively). We assume that the probability to tunnel in the three states is equal (same spatial distribution and tunnelling barrier, as their energy differences are negligibly small). The Maxwell–Bloch simulations were performed to guide the understanding of our Rabi-oscillation measurements. For this purpose, we disregarded non-Markovian effects^{20,24} and modelled the relaxation with a single phenomenological time constant *T*_{2}.

The delay time *t*_{S}, at which the RF pulses started, was fixed for one Rabi-oscillation sweep. The optimal *t*_{S} was experimentally determined by sweeping the timing of a π RF pulse over the range of the dwell pulse. A *t*_{S} > 0 is needed to initiate an imbalance between the T_{X} and T_{Z} states. Similarly, a decay time after the RF pulses is required such that the final triplet population is dominated by only one of these two triplet states. The optimal delay time is, therefore, shortly before the middle of the dwell voltage pulse. Furthermore, it is important that, on increasing the duration of the RF pulse, the sensitivity for differentiating T_{X} and T_{Z} does not greatly reduce, otherwise a further decay of the Rabi oscillations is induced by the readout. Therefore, we chose 30 μs as a delay time for the Rabi oscillations of pentacene-d_{14}, which were probed up to an RF pulse duration of 30 μs.

### Rabi oscillations baseline fit

The baseline of the Rabi-oscillation experiment represents the situation of equal populations in the coupled states T_{X} and T_{Z} during the pulse; even if the Rabi signal is not yet decayed, it is oscillating around the baseline. The decay of the baseline arises from the decay of the (on average) equally populated T_{X} and T_{Z} states into the singlet state during the RF pulse. As the final population of T_{Y} is independent of the RF signal, it will only give rise to a constant background and will be disregarded in the following.

Hence, the baseline is defined by the following: in the initial phase 0 < *t* < *t*_{S}, all three triplet states decay independently from each other. At the beginning of the RF pulse, that is, at *t* = *t*_{S}, the sum of populations in T_{X} and T_{Z} is

$${P}_{{\rm{XZ}}}({t}_{{\rm{S}}})={P}_{0}/3\left(\exp \left(-{k}_{{\rm{X}}}{t}_{{\rm{S}}}\right)+\exp \left(-{k}_{{\rm{Z}}}{t}_{{\rm{S}}}\right)\right),$$

in which \({k}_{{\rm{X}}}={\tau }_{{\rm{X}}}^{-1}\) and \({k}_{{\rm{Z}}}={\tau }_{{\rm{Z}}}^{-1}\) are the decay rates of T_{X} and T_{Z}, respectively, and *P*_{0} is the initial total population in the triplet state, such that *P*_{0}/3 is the initial population in each T_{X}, T_{Y} and T_{Z}. During the RF pulse, that is, for *t*_{S} < *t* < *t*_{E} (with *t*_{E} being the end of the RF pulse), the RF signal equilibrates (on average) the populations of two of the states, thus at the end of the RF pulse

$${P}_{{\rm{XZ}}}\left({t}_{{\rm{E}}}\right)={P}_{{\rm{XZ}}}\left({t}_{{\rm{S}}}\right)\exp \left(-\left({k}_{{\rm{X}}}+{k}_{{\rm{Z}}}\right)\left({t}_{{\rm{E}}}-{t}_{{\rm{S}}}\right)/2\right).$$

Finally, for *t*_{E} < *t* < *t*_{D}, the states decay again independently, giving at the end of the dwell time

$${P}_{{\rm{XZ}}}\left({t}_{{\rm{D}}}\right)={P}_{{\rm{XZ}}}\left({t}_{{\rm{S}}}\right)\exp \left(-\left({k}_{{\rm{X}}}+{k}_{{\rm{Z}}}\right)\left({t}_{{\rm{E}}}-{t}_{{\rm{S}}}\right)/2\right)\{\exp \left(-{k}_{{\rm{X}}}\left({t}_{{\rm{D}}}-{t}_{{\rm{E}}}\right)\right)+\exp \left(-{k}_{{\rm{Z}}}\left({t}_{{\rm{D}}}-{t}_{{\rm{E}}}\right)\right)\}/2,$$

which can be rearranged to

$$\begin{array}{l}{P}_{{\rm{XZ}}}({t}_{{\rm{D}}})={P}_{{\rm{XZ}}}({t}_{{\rm{S}}})\{\exp (-{k}_{{\rm{X}}}({t}_{{\rm{D}}}-{t}_{{\rm{S}}}))\exp ({t}_{{\rm{RF}}}({k}_{{\rm{X}}}-{k}_{{\rm{Z}}})/2)\\ \,\,\,\,+\exp (-{k}_{{\rm{Z}}}({t}_{{\rm{D}}}-{t}_{{\rm{S}}}))\exp (-{t}_{{\rm{RF}}}({k}_{{\rm{X}}}-{k}_{{\rm{Z}}})/2)\}/2.\end{array}$$

Note that *P*_{XZ}(*t*_{S}) does not depend on *t*_{RF} = *t*_{E} − *t*_{S} and therefore just represents a constant prefactor. The two terms provide contributions to the baseline that rise and fall exponentially with *t*_{RF}, respectively. For the specific case and parameters considered here, the prefactor of the rising term is much smaller than that of the falling term and is therefore neglected. Because the decay rates were determined (Fig. 1b) for the pentacene-h_{14} molecule, for which the Rabi oscillations were measured, these rates were used for the fitting of the Rabi oscillations of the pentacene-h_{14} molecule (Fig. 3a). In case of pentacene-d_{14} (Fig. 3c), we set (*k*_{X} − *k*_{Z})/2 = 0.012 μs^{−1} based on the measured decay rates of another individual pentacene-d_{14} molecule. In the experiment, other effects (for example, a thermal expansion owing to RF-induced heating) may also add to a temporal evolution of the baseline. These contributions were not separately accounted for but they are fitted as part of the falling term described above.

Source link