Title : Double Electron Spin Resonance of Engineered Atomic Structures on a Surface

Atomic-scale control of multiple spins with individual addressability enables the bottom-up design of functional quantum devices. Tailored nanostructures can be built with atomic precision using scanning tunneling microscopes, but quantum-coherent driving has thus far been limited to a spin in the tunnel junction. Here we show the ability to drive and detect the spin resonance of a remote spin using the electric field from the tip and a single-atom magnet placed nearby. Read-out was achieved via a weakly coupled second spin in the tunnel junction that acted as a quantum sensor. We simultaneously and independently drove the sensor and remote spins by two radio frequency voltages in double resonance experiments, which provides a path to quantumcoherent multi-spin manipulation in customized spin structures on surfaces.

One-Sentence Summary: Using a scanning tunneling microscope, we simultaneously control two spins using one tip, paving the way for multi-spin-qubit operations on surfaces.

Main Text:
Fabricating and controlling functional quantum devices with atomic precision is one of the long-term goals of quantum-coherent nanoscience (1). Quantum-coherent control of individual spins has been achieved in nanoscale systems, including dopants in semiconductors (2), color centers in insulators (3), semiconductor quantum dots (4), and molecular devices (5). However, building quantum systems with atomic scale control of their architecture and the microscopic interactions therein presents a major challenge (6). Individual spins, localized in atoms or molecules on surfaces, are a new candidate for realizing qubits, where atomically-precise fabrication can be achieved using a scanning tunneling microscope (STM). Quantum control of single atoms on surfaces driven by a magnetic tip has been demonstrated through continuous-wave (7)(8)(9)(10)(11) and pulsed (12) electron spin resonance (ESR) in STM. The STM can readily assemble multi-spin quantum systems (13,14) but use of ESR to coherently control a spin that is remote from the tunnel junction has not yet been considered (15)(16)(17)(18).
Here we demonstrate the ability to drive the resonance and sense the resulting spin state of a remote spin that is not directly in the tunnel junction. The key ingredient for the remote driving is a single-atom magnet placed in atomic proximity to the remote spin whose spatially inhomogeneous magnetic field converts the applied radio-frequency (RF) electric field into a driving magnetic field. To detect the quantum state of the remote spin, we positioned a nearby spin-resonant atom as a quantum sensor in the tunnel junction and performed double resonance experiments.
Experiments were performed in an ultra-high vacuum STM systems operating at 0.4 K (19). Two atomic spins were positioned on a surface so that they weakly coupled to each other with energy !,# (Figs. 1A and B). Both spins were a spin-1/2 hydrogenated Ti atom on bilayer MgO/Ag substrate (10,20,21), hereafter referred to as Ti-1 or "sensor spin", and Ti-2 or "remote spin". An Fe atom was placed in close proximity to Ti-2 to give a coupling strength of #,$% . The Fe atom's spin ( = 2) remained aligned perpendicular to the MgO surface due to its large easyaxis magnetic anisotropy and flipped up and down infrequently compared to the Rabi times used here (7,22). The resonance of the sensor spin was driven by the inhomogeneous magnetic field generated from the coupling with the tip (16,18), whereas the remote spin was driven by its interaction with the nearby Fe. Ti-1 Ti-2 To characterize the interactions in the three-atom structure, we first acquired ESR-STM spectra of the Ti spins in the structure by applying a single-frequency RF voltage to the tunnel junction ( Fig. 1C; figs. S1 and S2). The spectrum of Ti-1 (red curve in Fig. 1C) shows two ESR peaks that are separated by the coupling between Ti-1 and Ti-2 of !,# = 0.112 GHz, in good agreement with previous studies (20,21). Coupling between Ti-2 and Fe created additional To demonstrate the remote driving of Ti-2, we performed double-resonance ESR measurements by positioning the tip above Ti-1 and monitoring one of its two resonances while sweeping a second RF frequency across the Ti-2 resonances ( Fig. 2 and fig. S3C). This scheme is a nanoscale implementation of the double-resonance technique used in conventional ensemble ESR (23,24). At )$! = ↕' , the spectra showed two peaks, separated by the coupling energy !,# ( Fig. 2A). Note that a single-frequency ESR signal is negative in our measurements (Fig. 1C), therefore peaks in a double-resonance spectrum correspond to a reduction of the ESR signal of Ti-1. This reduction in the Ti-1 ESR signal stems from a transfer of state population from |00⟩ to |01⟩ due to resonant driving of the remote spin, . For comparison, the spectra obtained at )$! = ↕! (Fig. 2B) showed resonances at the same frequencies, but they appear as dips instead of peaks because of an increase in the ESR signal of Ti-1. We found that the resonance frequencies in these double-resonance spectra were ~400 MHz higher than those in the single-frequency ESR with the tip positioned on Ti-2 ( Fig. 1C), indicating that the tip's magnetic field (fig. S1) is absent when the Ti-2 spin is being remotely driven.
To show the origin of the double-resonance contrast more directly, we acquired a series of ESR spectra of Ti-1 with increasing )$# at one of the Ti-2 resonances (Fig. 1D). For )$# = 0, this experiment is identical to single-frequency ESR of Ti-1 (Fig. 1C), where a larger peak height is seen at ↕' ( ↕' ) than at ↕! ( ↕! ), with the ratio ↕' ↕! ⁄ given by the thermal population of the Ti-2 spin states. With increasing )$# , we observed a pronounced, simultaneous decrease of ↕' and increase of ↕! (Fig. 1D and fig. S4). Therefore, the double-resonance peaks and dips observed in Fig. 2 can be understood as the change of the ESR peak height with )$# in Fig. 1D (figs. S3 and S4).
When )$! was increased beyond the level used in Fig. 2, a splitting of the doubleresonance peaks was observed, whose spacing increased linearly with )$! (Fig. 3). This is a manifestation of the AC Stark effect (25), where the states |00⟩ and |10⟩ were split into two doublets by the strong interaction between the Ti-1 spin and its resonant RF field (19,26,27). The splitting of the |00⟩ state was spectroscopically probed by the '↕ transition (Fig. 3 inset). The splitting in the doublets directly corresponds to the Rabi rate Ω (!) of Ti-1, which we found to be proportional to )$! , allowing us to measure the Rabi rate of Ti-1, Ω (!) /(2π )$! ) = 0.160 ± 0.015 MHz/mV (Fig. 3B). By reversing the roles of the two RF voltages, the AC Stark splitting was also observed when sweeping )$! ( Fig. 1D and fig. S8 S14). We found a much longer lifetime for Ti-2, the remote spin, obtaining ! (#) = 150 ± 10 ns by fitting the double-resonance spectra (fig. S14) and the peak height ratio Fig. 1E; fig. S4D). In the Lindblad simulations we did not need to include a pure dephasing channel, indicating that ! processes were the dominant source of dephasing ( fig. S15). We note that the average interval between successive tunneling electrons at a tunnel current of 20 pA is 8 ns, suggesting the dominant role of tunneling electrons in the energy relaxation process of the sensor spin, Ti-1 (10). On the other hand, the large relaxation time of Ti-2 strongly indicates the absence of tunnel current-induced relaxation for the remote spin. In addition, an easily accessible RF voltage, )$# ≈ 40 mV (Fig. 1E), was enough to remotely drive the populations of Ti-2 states to saturation (Fig. 1E). An intuitive understanding of these population transfers can be developed using a rate equation approach (19).
To further investigate the underlying mechanism of the remote driving, we formed an isolated Ti-Fe pair and performed single-frequency ESR spectroscopy on Ti (Fig. 4A). We measured ESR spectra on Ti as a function of tip-sample distance and found two resonance peaks, corresponding to the two Fe spin states |⇑⟩ and |⇓⟩, respectively. We observed a strongly nonmonotonic change of the height for one of the peaks, which vanished at an intermediate tip-atom distance (Fig. 4, figs. S18 and S19). This non-monotonic behavior was not observed for isolated contributed to data analysis and manuscript preparation.

Competing interests:
The authors declare no competing interests.

Data and materials availability:
The data that support the findings of this study are available from the corresponding authors upon reasonable request.

Materials and Methods Supplementary Text
Figs. S1 to S20 Measurements were performed in three ultrahigh-vacuum (< 10 !" Torr) scanning tunneling microscopes (STMs). The data presented in the main figures and Figs. S1, S2, S4-S8, S18, and S20 were measured in a commercial 3 He-cooled STM (Unisoku, USM1300) at T = 0.4 K equipped with two-axis superconducting magnets. The data presented in Figs. S9-S13 were measured in a home-built STM based on a commercial cryostat (Janis Research, JDR-250) equipped with two-axis superconducting magnets, which was kept at 0.9 K during the measurement by Joule-Thomson cooling with a small amount of 3 He-4 He gas. The data presented in Fig. S19 were measured in a home-built 3 He-cooled STM at T = 1.1 K equipped with a singleaxis superconducting magnet. In this work, we chose field directions mostly along the sample plane with details labelled in the corresponding figures.
High-frequency transmission cables were installed on the STM systems as described in detail elsewhere (11). The single electron spin resonance (ESR) spectra were acquired by sweeping the frequency of an RF voltage #$ generated by an RF generator (Agilent E8257D) across the tunneling junction and monitoring changes in the tunneling current. To drive double resonance, the output signals of two RF generators (Agilent E8257D and E8267D) were combined using a power combiner/splitter (Mini-circuits, ZC2PD-K0244+) before combining with a DC bias voltage through a bias tee (SigaTek, SB15D2). Different lock-in detection schemes were designed to emphasize the desired signals in different measurements as described in Fig. S3.
Bilayer MgO was epitaxially grown on atomically clean Ag(001) single crystals by thermal evaporation of Mg in an O2 pressure of 1 × 10 !% Torr (7,11). In one STM system (Unisoku), Ti and Fe atoms were deposited onto a sample at approximately 100K, in the other two systems at ~10 K using standard e-beam evaporators. The Ti atoms on bilayer MgO are spin-1 2 ⁄ and are most likely hydrogenated (20). All measurements were performed on Ti atoms bound to a bridge site of the MgO (that is, in the middle of two oxygen sites). Before measurements, the STM tip made by Pt/It wire was poked into the Ag(001) surface until satisfactory topographic and spectroscopic features were observed on atoms on MgO, after which Fe atoms were picked up by the STM tip from MgO (by applying a DC voltage pulse of 0.6 V) to create a spin-polarized tip. The tip's spin polarization was calibrated with the asymmetry around zero bias in the dI/dV spectra of Ti on MgO (20).

Supplementary Text
Section 1. Single ESR spectra of an isolated Ti atom and Ti in the Ti-Ti-Fe structure in Fig. 1 To characterize effects of the tip used in collecting the main figures, we used the same tip to obtain a set of ESR spectra of an isolated Ti at different tunnel conductance (Fig. S1A). This allows us to extrapolate the ESR resonance frequency using a linear fit (Fig. S1B) considering the exchange interaction between tip and Ti, and to estimate the Ti resonance frequency in the absence of the tip magnetic field ( &'(,* ). The obtained &'(,* is 0.28 GHz higher than that under the tunneling conditions (20 pA, 50 mV) used for collecting the data in the main text. We also noticed that the ESR peak heights showed a monotonic decrease with deceasing tunnel conductance which vanished at zero conductance (Fig. S1C). This supports the theoretical claim (16) that the magnetic field gradient from the tip magnetic moment drives the ESR of an isolated Ti spin.   . 1C). Black curves are fits to the spectra using two Lorentzian curves on Ti-1 (red) and 4-Lorentzian curves on Ti-2 (blue), which yield the resonance frequencies labelled in the figure ( +, = 20 pA, +, = 50 mV, #$/ = 30 mV, 0.4 < < 0.5 K).

Section 2. Lock-in detection schemes used in double electron resonance experiments
We carefully designed the lock-in detection schemes in double electron resonance measurements to reveal the desired spectral information. We first consider the lock-in detection scheme used in the traditional single ESR-STM measurement of Ti-1 ( Fig. 1C and Fig. S3A).
Here the spin polarization signal was the difference between a driven state (lock-in A cycle) and the thermal state (lock-in B cycle). The measured ESR signal then revealed the effect of RF driving from the thermal equilibrium state (35).
The measurements performed in Fig. 1D can be similarly understood as the effect of RF driving of Ti-1 but from a non-thermal-equilibrium state. This non-equilibrium state was created by continuous driving of the remote spin Ti-2 with a second RF generator in both the lock-in A and B cycles (Fig. S3B). At zero Ti-2 driving ( #$/ = 0 mV), the measurement was identical to the single ESR measurement of Ti-1. With increasing #$/ , a non-thermal-equilibrium state was realized by driving Ti-2 transitions, and, on top of which, the Ti-1 ESR peak strengths were altered as shown in Fig. 1D and Fig. S4. In Fig. 1D, for example, we set #$/ at #$/ = *↕ , which induces population transfers from state |00⟩ to |01⟩. As a result, the Ti-1 transition at #$1 = ↕* (i.e., between states |00⟩ and |10⟩) becomes weaker (because the population of state |00⟩ is reduced by #$/ ), whereas the Ti-1 transition at #$1 = ↕1 (i.e., between states |01⟩ and |11⟩) becomes stronger (because the population of state |01⟩ is increased by #$/ ).
The measurements performed in Figs. 2 and 3 are different from the above measurement in that we perform a frequency sweep across the range of the Ti-2 resonances (i.e., the remote spin) rather than those of Ti-1 (Fig. S3C). These measurements can be thought of as spectroscopically revealing the frequencies at which #$/ reaches Ti-2 resonances. Due to the absence of the tip magnetic field, these resonance frequencies were significantly different from ESR spectra collected with the tip above Ti-2. To detect the remote spin resonance, we kept monitoring the Ti-1 ESR strength (by driving it in both lock-in A and B cycles) but drove Ti-2 in the A cycle only. When the Ti-2 driving was off resonance, little change was induced in the Ti-1 ESR signal between the A and B cycles, which resulted in a null signal in the double resonance spectrum.
When the Ti-2 driving reached its resonances in Fig. 2, population was transferred and the Ti-1 ESR strength became different between the A and B cycles, thus resulting in a nonzero signal as shown in Fig. 2. More explicitly, consider the case where we monitor the Ti-1 transition at #$1 = ↕* (i.e., between states |00⟩ and |10⟩). When we additionally drive the Ti-2 resonance at *↕ (lower frequency peak in Fig. 2A), population is further transferred from state |00⟩ to |01⟩, which results in less population at state |00⟩ and reduces the Ti-1 ESR signal in the lock-in A cycle (compared to the B cycle). Since the Ti-1 ESR signals in the single frequency sweep are negative in our setup (Fig. 1), a reduction of the negative single-frequency ESR signal leads to a positive double resonance peak as seen in Fig. 2A. Similarly, when we monitored a different Ti-1 transition at #$1 = ↕1 whilst still driving the Ti-2 resonance at *↕ (lower frequency peak in Fig. 2B), the population transfer from state |00⟩ to |01⟩ (from Ti-2 driving) resulted in more population at state |01⟩ and therefore enhanced the Ti-1 ESR signal in the lock-in A cycle (compared to the B cycle). This leads to an enhancement of the negative single Ti-1 ESR signal and thus a negative double resonance peak (i.e., a dip) as seen in Fig. 2B.
The Ti-2 driving at the higher resonance frequency 1↕ (higher frequency peaks in Fig.  2A and B) connects two excited states and involves more complicated interplay between driving and relaxation, which can yield a peak, a dip, or more intricate lineshapes depending on the parameters. This, however, is captured by our Lindblad simulations explained below (Sections 5.1-5.3) and shown in the main figures. We note that the Lindblad simulations can be simplified into a simpler rate-equation model (Sections 5.4 and 5.5), which serves as the basis for the population-transfer-based explanations provided above and in the main text.  Fig. 1C. The spin polarization signal was contrasted between a driven state (lock-in A cycle) and the thermal state (lock-in B cycle) while the RF frequency #$1 was swept. (B) Lock-in detection scheme of the double resonance experiment in Fig. 1D. This scheme is similar to A but based on a non-thermal-equilibrium state created by a continuous second driving at a fixed frequency #$/ at a Ti-2 resonance. (C) Lock-in detection scheme of the double resonance experiment in Figs. 2 and 3. In this scheme, #$/ was swept across the range of Ti-2 resonances. Through monitoring the amplitude change of a Ti-1 spin resonance, we could remotely detect the frequency at which #$/ reaches Ti-2 resonances. Section 3. Supplementary data to the Ti-Ti-Fe structure presented in the main text Section 3.1. ESR peak height in double resonances with fixed RF2 frequency (Fig. 1D) In addition to the double resonance spectra at #$/ = *↕ shown in Fig. 1D, we also acquired the double resonance spectra at #$/ = 1↕ and an off-resonance frequency of Ti-2, as shown in Figs. S4B and C. To extract the ESR peak heights in these spectra, we performed a 2-Lorentzian curve fit to each spectrum (black curve). The extracted peak height ratios ( ↕* ↕1 ⁄ ) were plotted as a function of #$/ (the RF driving voltage of Ti-2) in Figs. S4A-C and summarized in Fig. S4D.
The change of ↕* ↕1 ⁄ is noticeably smaller at #$/ = 1↕ than at #$/ = *↕ for the same driving amplitude due to the more significant relaxation effects when driven at 1↕ (this complicated interplay between driving and relaxation is included in our Lindblad simulations, see below).
When driving the Ti-2 spin at an off-resonance frequency (which does not alter the populations), the double resonance spectrum should be equivalent to single-frequency ESR of Ti-1, and one expects a constant ↕* ↕1 ⁄ ratio independent of #$/ . However, in this data set, we observed a slight monotonic decrease of ↕* ↕1 ⁄ shown in Fig. S4C. This stems from the RF power heating of the spin system. To capture these RF heating effects in Figs. 1D and S4, we recorded the sample temperature at each RF power and frequency and used this temperature reading for the simulation of each ESR curve (Fig. S5). A minor overall up-shift of temperature (0.04 K) was included in the simulations to capture the effective spin temperature (which is reflected in the correct intercept for the peak ratio at VRF2 = 0 in Figs. 1E and S4D).
When performing the same experiment as Fig. 1D using another STM at 0.9 K (Fig. S10), the temperature variations were much reduced due to the higher thermal stability at this temperature of this STM as evidenced by the near-constant peak ratio at off resonance. In that case, only a fixed temperature was used in the theoretical simulations to match with the experiments, highlighting the consistence between the experiment and the simulation. Section 3.2. AC Stark effect in double resonance spectra in Fig. 3 In this work, we utilized the AC start effect to directly obtain the Rabi rates of both Ti-1 and Ti-2 spins. An elegant approach to understand the AC Stark effect is the "dressed atom" theory as follows. Consider a quantum system of two atom levels, | ⟩ and | ⟩ ( |3⟩ < |5⟩ ), in a near-resonant electromagnetic field with photon energy ℏ , so that |5⟩ − |3⟩ ≈ ℏ . If photons do not interact with the atom states, the total Hamiltonian is A * = A 6.78 + A 9:';< , and the eigenstates of A * are | , ⟩ and | , ⟩ where labels the photon occupation. We notice that the state of the atom level with the photon number − 1 (| , − 1⟩) is almost degenerate with the state | , ⟩, and so are | , ⟩ and | , + 1⟩. We define the detuning as Δ = |5⟩ − |3⟩ , and the two states | , ⟩ and | , − 1⟩ become fully degenerate at Δ = 0. Now if we turn on an interaction A :=. between atom and field, the states | , ⟩ and | , − 1⟩ become hybridized and their degeneracy lifted. Then, the eigenstates of the system can be described by a linear combination of the two nearly degenerate states of the non-interacting Hamiltonian | ⟩ = 3 | , ⟩ + 5 | , − 1⟩, where the coefficients are obtained from the Schrodinger equation, ( A 6.78 + A 9:';< + A :=. )| ⟩ = + | ⟩.
In our double resonance experiment, four spin states are involved. However, in each experiment, only one ESR transition was strongly driven. For example, only the transition ↕* between states |00⟩ and |10⟩ as shown in Fig. 3A. As a result, the ESR transition between the states |00⟩ and |01⟩ occurs at two distinct frequencies *↕ (±) = *↕ ± Ω (1) 2 ⁄ when probed by another weak electromagnetic field that was swept across the Ti-2 resonance *↕ . The splitting allows us to determine the Rabi rate of Ti-1 as Ω (1) /(2 #$1 ) = 0.160 ± 0.015 MHz/mV. The ESR peak heights at *↕ (±) are sensitive to the detuning Δ due to its influence on the coefficients sin and cos in Eqs. S4 and S5, which is also seen in our simulations in Fig. S16. We omit the quantum numbers of the photon states for convenience in the main figures and text. Section 3.3. AC Stark effect in double resonance with the measurement scheme in Fig. 1D To extract the Rabi rate of Ti-2 (the remote spin), we need to probe a Ti-1 ESR transition while strongly driving a Ti-2 transition, similarly to the scheme in Fig. 1D but with higher #$/ . These data are presented Fig. S8 where we observed ESR peak splitting at #$/ above 50 mV. Such splitting is due to the dressing by photons at frequency *↕ , which split the spin states |00⟩ and |01⟩, as depicted in Fig. S8C. The peak splitting Δ showed a linear dependence on #$/ (Fig. S8B). A linear fit to Δ ( #$/ ) for the lower (higher) frequency peak gave the Rabi rate of the Ti-2 spin, Ω (/) (2 #$/ ) ⁄ = 0.262 ± 0.011 MHz mV ⁄ (0.225 ± 0.016 MHz mV ⁄ ). We used the average of the two Ω (/) (2 #$/ ) ⁄ values in the simulations of the double resonance spectra. We also obtained an additional data set with the same tip but at a slightly different magnetic field with a smaller out-of-plane component, resulting in a slightly smaller Rabi rate of Ti-2 (Figs. S8D and F). This observation indicates that the interaction between Ti-2 and Fe decreases with smaller z-component of the magnetic field, leading to a weaker driving of the Ti-2 spin and smaller AC Stark splitting.

Section 4. Double resonance measurements of a second Ti-Ti-Fe structure
To verify the universality of remote spin driving, we performed double resonance measurements of a second Ti-Ti-Fe structure using a different STM kept at 0.9 K as shown in Figs. S10-S13. The magnetic field was kept at 0.85 T with a ~12° angle off the sample plane (i.e., the magnetic field was mostly along the sample plane). Fig. S9 shows the structure and single ESR spectroscopic characterization of the structure. Figures S10 and S11 present the results measured on the second structure using the measurement scheme of Figs. 1D and S4. Here, due to the higher data acquisition temperature (0.9 K) and the thermal stability of this STM, the sample temperature did not vary as much with increasing RF power (as shown by the near-constant off-resonance peak ratio in Fig. S10D), and thus only a fixed spin temperature was needed in the simulations to match the experimental data (see the discussion in Section 3). This further highlights the agreement between the experiments and the simulations. Figure S12 plots the double resonance spectroscopy results (similar to Fig. 2) at different tip locations on Ti-1. Despite significant shifts of the single ESR peaks of Ti-1 due to the variations of the local tip magnetic field, the frequency of the double resonance peak/dip features stay mostly unaltered, indicating the absence of the tip's magnetic field on the remote spin. Figure  S13 illustrates the reverse double resonance spectroscopy measurement, where the STM tip was placed on Ti-2. We attempted to drive the Ti-1 ESR transitions remotely, but failed as expected, which indicates that Ti-1 cannot be remotely driven due to its negligible coupling to Fe (i.e., absence of a single magnet for remote driving).
Lindblad simulations of this data set show 1 (1) = 5 ns, Ω (1) (2 #$1 ) ⁄ = 0.080 MHz mV ⁄ , 1 (/) = 50 ns, Ω (/) (2 #$/ ) ⁄ = 0.0644 MHz mV ⁄ , and spin temperature = 0.96 K. The remote Rabi rate is smaller in this structure due to the slightly larger distance between Fe and Ti-2 (0.72 nm) compared to the structure studied in the main text (0.59 nm). The shorter 1 (/) is likely due to the higher measurement temperature and/or differences in the local magnetic environment.     where is the density matrix and is the Hamiltonian of the closed system. A real-world quantum system, on the other hand, is always coupled to an environment. To include decoherence effects induced by the environment, a commonly used method is to invoke the Lindblad master equation. The Lindblad formalism assumes a Markovian interaction between the system and the environment and introduces a series of Lindblad operators E that describe such interactions and the resultant non-unitary evolution of the system (29). With the unitary and nonunitary evolution combined, the celebrated Lindblad master equation reads As a simple example, consider the Lindblad formalism of a single spin-1/2 coupled to a thermal bath at temperature . Three relaxation paths can be represented by three Lindblad operators: • Pure dephasing is represented by The effect is that off-diagonal terms exponentially decay towards zero at the rate of Γ G in the absence of the drive, as expected.
• Energy relaxation from |0⟩ to |1⟩ is represented by ! = pΓ ! ! , which yields The effect is that ** exponentially decays towards zero at the rate of Γ ! , and a dephasing process also occurs at the rate of Γ ! /2.
• Energy relaxation from |1⟩ to |0⟩ is represented by > = pΓ > > , which yields (Eq. S10) The effect is that 11 exponentially decays towards zero at the rate of Γ > , and a dephasing process also occurs at the rate of Γ > /2.
Combining all three terms together, the Lindblad equation for a single spin becomes Here we grouped all the off-diagonal decay terms together, which defines a total dephasing rate The ratio between Γ > and Γ ! is not arbitrary but fixed by the thermal Boltzmann distribution. To see this, consider the steady-state value of ** in the absence of drive: ** = Γ > 11 − Γ ! ** = Γ > − (Γ > + Γ ! ) ** .
(Eq. S12) This means that the steady state in the absence of the drive is given by which should be equal to the thermal population where ℏ * is the energy difference between states |1⟩ and |0⟩. We are thus motivated to define a total energy relaxation rate 1 (Eq. S16) Using the single-spin Hamiltonian in the rotating frame = ℏ / z −∆ Ω Ω ∆ |, where ∆ is the detuning frequency and Ω is the Rabi frequency on resonance, the master equation becomes %.
(Eq. S17) This is the expression of the master equation in the Fock-Liouville space and can be shown to be equivalent to the Bloch equations that governs the single spin dynamics.
Section 5.2. Simulation of two coupled spins using the Lindblad formalism In the lab frame, the Hamiltonian describing a coupled spin-1/2 system is given by where ⨂ represents the tensor product, is the identity operator, *1 = − 1 * and */ = − / * are the Larmor frequencies of the spin 1 and 2, Ω (1) = − 1 1 1 and Ω (/) = − / 1 / are the on-resonance Rabi frequencies that represent the strengths of the RF waves, and LM1 and LM/ are the applied RF frequencies. In our weakly-coupled spin pair, is small compared to the energy scale of *1 − */ . The prefactor of ensures that the resulting ESR splitting is . The two RF waves are assumed to be applied in direction.
To construct the Lindblad formalism for two coupled spins, we need to consider different relaxation paths. The processes that originate from single-spin relaxation can be formulated similar to the single-spin case shown in the last section: • Pure dephasing of the first spin is described by G1 = k The four rows and columns are labelled according to |00⟩, |01⟩, |10⟩, and |11⟩.
• Pure dephasing of the second spin is described by G/ = k (Eq. S20) • Energy relaxation by the first spin flipping from |0⟩ to |1⟩ is represented by 1! = pΓ 1! 1! ⨂ , which yields Here we ignored the dependence of the first spin relaxation rate on the second spin state, which should not be important unless becomes the dominant energy scale.
• Energy relaxation by the first spin flipping from |1⟩ to |0⟩ is represented by 1> = pΓ 1> 1> ⨂ , which yields (Eq. S22) • Energy relaxation by the second spin flipping from |0⟩ to |1⟩ is represented by (Eq. S23) • Energy relaxation by the second spin flipping from |1⟩ to |0⟩ is represented by (Eq. S24) • There are four more energy relaxation paths that correspond to the cross-relaxation processes of the two spins. However, our main experimental results are well captured without considering these cross-relaxation processes. The unimportant role of the flip-flop processes can be justified by the very different Larmor frequencies of the two spins in the experiment. Similar to the single-spin case, Γ 1! /Γ 1> and Γ /! /Γ /> are determined by the thermal population. In the absence of the drive, the steady thermal populations of |00⟩, |01⟩, |10⟩, and |11⟩ are given by and the relaxation rate ratios are thus determined to be Γ 1> Γ 1> + Γ 1! = 1 1 + exp(ℏ *1 / I ) , Here we ignored the change of the Zeeman splitting induced by the exchange interaction , which would only induce a negligible difference (less than 1%) of the Boltzmann factors in our parameter regime. We then define a total energy relaxation rate for each spin as which are used as tuning parameters in the simulations. Section 5.3. Numerical solution to the Lindblad master equation in the lab frame Inserting Eq. S19-S25 into the Lindblad equation yields the master equation for the two-coupled-spin density matrix in the lab frame. The time evolution of the coupled two-spin density matrix was simulated using the Lindblad master equation solver as implemented in QuTip (28). In our setup, the time-evolution is performed by integration of the time-dependent Hamiltonian (Eq. S18) for 2000 ns using a step-size of 0.02 ns.
To obtain the steady-state populations, the density matrix was averaged for the last 10% of the simulation time to remove residual oscillations from the calculation and the ESR signal was calculated according to Eq. S29 below. The frequency resolution was 2 MHz corresponding to 150 discrete frequency steps for all spectra shown in the main text.
To reproduce the double resonance spectra, we notice that the observable in our double resonance experiment is the steady-state spin polarization of the spin 1 (the sensor spin under the STM tip), which is given by The double resonance spectra shown in the main text were then obtained by calculating the difference of SP1 between the lock-in A and B cycles (as discussed in the Section 2 above). We note that our simulations has 8 parameters to describe the driven two spin-system: the two Larmor frequencies, the coupling strength between spin-1 and spin-2 (J12), the Rabi rates of spin-1 and spin-2 (Ω (1) and Ω (2) ), and T1 times of spin-1 and spin-2 (T1 (1) , T1 (2) ) as well as an effective spin temperature T.
The Larmor frequencies and the coupling strength were determined from the double resonance experiments. For the simulations of the data presented in the main text, the dependence of the Rabi rates on the RF voltages was determined from the AC Stark splitting (Fig. 3, Fig. S7, and Fig. S8). Analytical temperature functions were implemented according to Fig. S5, and an additional, constant temperature offset of ~40 mK was added to match the peakheight ratio of double resonance spectra in thermal equilibrium.
The remaining two parameters were determined as follows: T1 (2) can be determined with a high accuracy from the peak-height ratios of the double resonance experiments when driving the remote spin at *↕ = 15.071 GHz (Figs. 1D and E). By a series of calculations in the intermediate RF voltage regime (0 < #$/ < 40 mV) we found that a value of T1 (2) = 150 ± 10 ns minimizes the deviations between experiment and simulation (Fig. S14B). T1 (1) was determined by comparing simulated and experimental spectra obtained at #$/ = 40 mV (Fig. S14A) to be T1 (1) = 8 ± 1 ns.
The excellent agreement between experiment and calculations confirms that we are indeed observing an electron-electron double resonance experiment with a narrowly determined set of parameters. The simulations also allow us not to consider pure dephasing (Γ % , see Eqs. S19 and S20) as additional relaxation channels, since all experimental observations can be reproduced with T1 relaxation processes only. The effects of additional pure dephasing on Ti-2 and Ti-1 on the spectra are shown in Fig. S15.
For the second data set of double resonance experiments (section 4) the driving strength could not be obtained from the splitting. In this case we used both the T1 (1) , T1 (2) as well as Ω (1) and Ω (2) as free parameter and optimized them in a similar manner as discussed above.
We note that a small detuning (∆ < 10 MHz) can change the relative peak heights of the Autler-Townes doublet as shown in Fig. S16 (compared to Figs. S8A and D). Such detuning can occur in the experiment due to imperfect RF generator or magnet settings. ⁄ : = I sim ) and experiment (I exp ) expressed as the root mean square deviation (RMSD) for the intermediate voltage range (0 < #$/ < 40 mV) as function of T1 (2) giving 150 ± 10 ns as the best estimate.  To obtain an intuitive understanding of the competition between driving and relaxation in our system, in this and the next subsections we derive a rate equation that involves only population terms. Before diving into the two-coupled-spin case, in this subsection we first consider how to derive a rate equation for a single spin from its master equation (Eq. S17).
To obtain a rate equation involving only population terms, we take the adiabatic approximation (36) where the coherence terms are assumed to respond so fast to changes in populations that the coherences remain in the instantaneous steady states set by the populations, or, 1* / = *1 / = 0 (this approximation is fully satisfied at the real steady states). Using the second and third lines of Eq. S17, this approximation allows us to express coherences with populations as where Γ G .7. = 1/ / = Γ > /2 + Γ > /2 + Γ G . Inserting Eq. S29 to the first and fourth lines of Eq. S17 yields the rate equations for a single spin ** ( ) = −( + Γ ! ) ** ( ) + ( + Γ > ) 11 ( ), where the effective driving is given by The driving is strongest at zero detuning ∆ = 0 and weakens when ∆ becomes comparable to the relaxation rate Γ > /2 + Γ > /2 + Γ G . This rate equation has a straightforward interpretation by simply considering all possible paths of population transfers as shown in Fig. S17A. For example, the population of state |0⟩, ** , can be decreased by transferring populations from state |0⟩ to |1⟩ via either driving or relaxation Γ ! (yielding the negative rate −( + Γ ! ) ** in Eq. S30), and ** can be increased by transferring populations from state |1⟩ to |0⟩ via either driving or relaxation Γ > (yielding the positive rate ( + Γ > ) 11 in Eq. S30). The balance of driving and relaxation yields the steady state solution at ** / = 11 / = 0.

Section 5.5. Rate equation of two coupled spins
To obtain a rate equation for our coupled-spin system, we first cast its master equation timeindependent by a rotating frame transformation and then focus on the population terms.
The obtained rate equations provide an intuitive way to further understand our observations. The peak widths of the double resonance peaks, for example, were not much narrower than single frequency STM-ESR peaks in our data. This is because the detuning of the remote spin Ti-2 in Eqs. S46 and S47 needs to be compared to a combined relaxation rate that involves not only the remote spin relaxation rate Γ /± (which is small) but also the sensor spin relaxation rate Γ 1± (which is large and comparable to the single STM-ESR relaxation rates).
In our continuous-wave experiments presented in this manuscript, the steady states rely on an intricate interplay between relaxation rates Γ's and the effective driving rates 's. Consider, for example, the double resonance scheme in Fig. 2 where we monitor the sensor spin Ti-1 transition ↕* between states |00⟩ and |10⟩. If we strongly drive the remote spin Ti-2 transition *↕ , O becomes dominant over the relaxation rates Γ /± , which tends to make the populations of |00⟩ and |01⟩ equal, i.e., causing a population transfer from state |00⟩ to |01⟩. This further causes a transfer of population from state |10⟩ to |00⟩ due to Γ 1! relaxation and drive O . As a result, we expect to have a smaller spin polarization signal on Ti-1 (SP1 = ** PP + 11 PP − // PP − OO PP ), mostly because the population of |10⟩ ( // PP ) is reduced. This results in the experimental contrast we show in Fig. 2. This type of population-based analysis scheme was used throughout the manuscript.