Coherent spin-spin coupling mediated by virtual microwave photons

We report the coherent coupling of two electron spins at a distance via virtual microwave photons. Each spin is trapped in a silicon double quantum dot at either end of a superconducting resonator, achieving spin-photon couplings up to around $g_s/2\pi = 40 \ \text{MHz}$. As the two spins are brought into resonance with each other, but detuned from the photons, an avoided crossing larger than the spin linewidths is observed with an exchange splitting around $2J/2\pi = 20 \ \text{MHz}$. In addition, photon-number states are resolved from the shift $2\chi_s/2\pi = -13 \ \text{MHz}$ that they induce on the spin frequency. These observations demonstrate that we reach the strong dispersive regime of circuit quantum electrodynamics with spins. Achieving spin-spin coupling without real photons is essential to long-range two-qubit gates between spin qubits and scalable networks of spin qubits on a chip.


I. INTRODUCTION
There is tremendous interest in the realization of quantum computers, and architectures based on solid-state devices offer significant advantages to achieve this goal. Circuit quantum electrodynamics (QED) leverages highquality-factor superconducting resonators at cryogenic temperatures to enable the coupling and readout of superconducting qubits [1][2][3]. Meanwhile, spin qubits in gate-defined semiconductor quantum dots (QDs) are also promising for quantum computing [4,5], having achieved high-fidelity quantum operations, long coherence and relaxation times, and operation above 1 kelvin. Spin qubits in silicon could eventually leverage the advanced manufacturing capabilities of the microelectronics industry, which is a compelling argument towards their development.
Despite these groundbreaking realizations, some important hallmarks of circuit QED experiments have remained elusive for spins, in part due to the insufficient * Correspondence to: P.Collard@USherbrooke.ca † Correspondence to: L.M.K.Vandersypen@tudelft.nl spin-resonator interaction strength, in combination with decoherence and fabrication challenges. These hallmarks include dispersive interaction between two spins mediated by virtual resonator photons [3] and photon-numberdependent spin dispersive shifts [32], both requiring a higher level of interaction-to-decoherence ratio than previously achieved [30]. The former is required for most two-qubit gate schemes and arguably represents the next frontier of the field, while the latter enables higher signal with dispersive readout, and photon state measurement [33] and universal control [31,34].
In this work, we overcome previous challenges and demonstrate both spin-spin interaction mediated by virtual photons and photon-number-dependent spin dispersive shifts using single spins at either end of a highimpedance superconducting resonator. Each single spin is trapped in a double quantum dot (DQD) formed in a 28 Si/SiGe heterostructure, and tunable spin-charge hybridization is enabled by a micromagnet. We first reach the resonant strong spin-photon-spin coupling regime; then, we bring both spins in resonance with each other but detuned from the resonator photons, and we observe a spin-spin avoided crossing showing coherent remote interaction. This differs from previously reported work where the virtual coupling could not be achieved [30]. Finally, we resolve the photon-number states from the discrete shifts they induce in the spin transition frequency.

II. METHODS
The device is shown in Fig. 1. A 5-to-7-nm-thin film of NbTiN is deposited on the surface of a 28 Si/SiGe heterostructure and patterned to form the superconducting resonator, ground planes, gate filters [35,36], and gate fanout lines (Fig. 5). The sheet kinetic inductance (around 140 pH/ ), the narrow width of the resonator center conductor (170 nm), and the retracted ground planes combine into a high effective resonator impedance, Z r ≈ 3.0 kΩ [36,37]; this results in a fixed coupling g c /2π = 192 MHz to the DQD charge degree of freedom ("charge qubit") at a resonator frequency of ω r /2π = 6.916 GHz for the half-wave mode (Fig. 1a). Achieving a large g c ∝ α c V zpf (approximately 5 times larger than Ref. [30]) is a result mostly of the increased voltage zero-point fluctuations V zpf ∝ ω r √ Z r of the high-impedance resonator, and of the (0, 1)−(1, 0) interdot transition lever arm α c [17]. A combination of loss mitigation strategies (Appendix A) results in an undercoupled resonator with a linewidth κ r /2π = 1.8 MHz, limited by resistive or dielectric losses near the DQDs. Each end of the resonator terminates as one of the dot plunger gates and is biased through a DC tap. The resonator is only 250 µm long, which is a consequence of the substantial kinetic inductance that translates to a high effective magnetic permeability [36]. The DQD potential is shaped by applying suitable voltages to surface gate electrodes, as shown in Fig. 1b-c. Cobalt micromagnets provide a transverse magnetic field difference ∆B ⊥ = 42 mT between the two dots while minimizing magnitude differences ∆B (∆B ⊥ is taken as constant, see Appendix D for details). Unlike g c , which is, to a large extent, fixed by the device structure, the spin-photon coupling g s /2π at zero charge detuning approximatively scales as for 2t c { ω s , ω r }, and it is therefore tunable via the DQD tunnel splitting 2t c [22,23]. Here, g e = 2 is the electron Landé g-factor, µ B is the Bohr magneton, and ω s /2π is the spin transition frequency. In any case, g s ≤ g c . These device characteristics combine to enable values of g s /2π up to around 40 MHz in this work. The micromagnets are tilted ±15 • relative to the vertical direction (Fig. 1d), allowing one to fine-tune each spin's Zeeman energy by rotating the in-plane external magnetic field of magnitude B r = |B ext | by an angle φ [30] using a vector magnet. Device fabrication and experimental setup details are given in Appendix A and Appendix B, respectively.

III. RESULTS
First, the system is tuned to the regime where both spins are resonant with the photons [30], which happens at a field angle of φ = 10.5 • due to a 6.5-mT difference in the micromagnet fields at φ = 0 • . In the main panels of Figs. 2a and 2b, DQD1 and DQD2 are set separately to zero charge detuning, = 0 µeV, allowing each spin to interact with the photons while the other is decoupled. For what follows, it is useful to keep in mind that the effective charge-photon coupling (and therefore the effective spinphoton coupling) can be switched off simply by biasing the DQD charge detuning ( t c ) such that the electron is confined to one dot. Sweeping B r results in a vacuum Rabi splitting measurement for DQD1 and DQD2. The device used in this experiment has slow drift of the (0, 1)−(1, 0) interdot charge transitions. Hence, for every magnetic field setting (B r or φ), an automated fitting procedure is used to extract a data cut along = 0 µeV and reconstruct the two-dimensional data. This procedure ensures that the measurements are protected against long-term drift, and it is further detailed in Appendix C. A series of measurements is used to successively constrain and extract system parameters, and the calculated transition frequencies are then plotted over the measured data. The model parameters include the resonator frequency ω r /2π, the charge-photon coupling g i c /2π, and the DQD charge detuning and tunnel coupling i /h and t i c /h (respectively). Here, h is the Planck constant, and i is an index identifying the DQD. Full details about the model are given in Appendix D. Experimentally, we observe that there is a high level of symmetry between the two DQDs, and unless mentioned explicitly, we omit the DQD index and take the dot parameters to be the same. From the vacuum Rabi splitting measurements of Fig. 2a and Fig. 2b, we extract a spin-photon coupling of (g 1 s , g 2 s )/2π = (11.8, 11.0) MHz for DQD1 and DQD2, respectively, with tunnel splittings (2t 1 c , 2t 2 c )/h = (13.2, 13.7) GHz. Both DQDs achieve the strong spin-photon coupling regime, i.e., g s > {κ r , Γ s }, where Γ s ≤ 6 MHz is the spin linewidth. When both DQDs are set to zero charge detuning simultaneously (Fig. 2c), both spins interact with the resonator, yielding a larger resonator dispersive shift due to the additive effects of the two DQDs, as well as an enhanced vacuum Rabi splitting 2g 12 s /2π = 32.3 MHz from the two-spin ensemble. The enlarged splitting matches well the predicted 2 (g 1 s ) 2 + (g 2 s ) 2 /2π value, very close to a factor √ 2 larger than for single spins in this case. The Hamiltonian model transitions, which are calibrated solely on the one-at-atime interaction data, predict very well the outcome of the simultaneous interaction. The √ 2 enhancement indicates simultaneous, coherent, and resonant interaction of both spins with the resonator, as demonstrated in prior work [30].
To resolve the spin-spin exchange splitting 2J/2π mediated by virtual photons, a larger spin-photon coupling is needed than in the resonant case (see, e.g., Fig. 9). Given a fixed ratio g s /∆ s , where ∆ s = ω s −ω r is the spin-photon detuning, and given that J ≈ (g s ) 2 /∆ s in the dispersive regime, increasing g s should allow J to become larger than the spin linewidth. The spin-photon coupling strength is increased by reducing the DQD tunnel splitting to 2t c /h ≈ 8.8 GHz, yielding g s /2π ≈ 33 MHz. Insight into how the two spin states hybridize is gained by measuring the spin transition frequencies as a function of the external magnetic field angle φ using two-tone spectroscopy. A pump tone at frequency f pump is sent down a gate line to each DQD to generate an excited spin-up population, while the transmission coefficient S21 is probed at a fixed frequency f probe set to the dispersively shifted resonator frequency for each DQD at zero charge detuning (e.g., as in the insets of Fig. 2). Line cuts along the charge zerodetuning axis of data like those of Fig. 8 are assembled in a two-dimensional diagram, resulting in Fig. 3. The spin transition of the independently interacting DQD1 and DQD2 is visible as a dip in the S21 signal magnitude. The slope in the spin transition frequency as a function of φ is caused by the relative angle of the field and each micromagnet, which allows one to tune each spin's transition energy, as explained earlier. The DQD2 slope is smaller than the one of DQD1 because φ ∈ [9,13] • is almost aligned with the DQD2 micromagnet angle 15 • and farther from the one of DQD1. When both spins interact simultaneously (Fig. 3c), an avoided crossing is observed, while the upper transition becomes dark close to spin-spin resonance [3]. With (2t 1 c , 2t 2 c )/h = (8.82, 8.80) GHz and B r = 52 mT, we extract (g 1 s , g 2 s )/2π = (32.4, 32.7) MHz, 2J/2π = 19.0 MHz, and ∆ 2s /2π = −79 MHz (here, ∆ 2s is for two spins since each DQD contributes an additive charge dispersive shift χ c to ω r , as explained in Appendix D). In Fig. 3d and Appendix F, additional spin-spin hybridization results show that the exchange interaction is reduced when either |∆ s | is increased or g s is decreased, as expected. The full width at half minimum of the dip, 2Γ s /2π, is 2Γ − s /2π = (11.7 ± 3.6) MHz for the two-spin lower branch on resonance (the upper branch's visibility is too low for a reliable fit). While this is the most relevant linewidth, for completeness, we also report (2Γ 1 s , 2Γ 2 s )/2π = (10.4 ± 3.0, 13.0 ± 3.2) MHz for the individual spin-interaction case on resonance (φ = 10.5 • ), and (2Γ 1 s , 2Γ 2 s )/2π = (8.8 ± 4.1, 13.2 ± 4.2) MHz for the simultaneous interaction away from resonance (φ = 13 • ). Since 2J > (2Γ 1 s + 2Γ 2 s )/2, the two spins are coherently hybridized through virtual photons, achieving a longstanding goal for the field. Arguably, the data presented in this figure are not very deep into the dispersive regime. In Appendix F, we present additional data at larger ∆ s where 2J is still larger than 2Γ s . The ratio of interaction strength to decoherence should be sufficient to enable two-qubit gates [18][19][20][21] in future experiments.
The photon-number-dependent dispersive shift, where a single (probe) photon will shift the qubit frequency by more than its linewidth, is also a hallmark of circuit QED [32]. The shift 2χ s /2π is expected to scale as χ s ≈ (g s ) 2 /∆ s in the dispersive regime, a dependence reminiscent of the qubit-qubit exchange interaction J. It therefore seems reasonable to observe both effects if the linewidths are sufficiently narrow. However, while virtual spin-spin coupling requires zero photons and therefore has lower sensitivity to photon losses, this effect requires finite photon population, which exposes the spin to measurement broadening (see the Γ (n) formula below). The photon-number-dependent dispersive shift of DQD1 is shown in Fig. 4. A larger probe power is used (−119 dBm, 0 MHz detuning) than for the data of Fig. 3 (−123 dBm, −1.2 MHz detuning), which populates the resonator with more photons. An extra dip appears below the main transition, and its frequency shows good agreement with the prediction of the Hamiltonian model for the |↓, 1 ↔ |↑, 1 spinlike transition. In Fig. 4b, a line cut from Fig. 4a is extracted, and the dip areas and separations are fit to a sum of Lorentzian dips with widths 2Γ (n)/2π = 2γ s /2π+(n+n r )κ r /2π. A value of 2χ s /2π = −13.1 ± 2.2 MHz is extracted from the fit (with parameters κ r /2π = 3.0 ± 0.2 MHz,n r = 0.62), slightly larger (c) Relative area of the dips for φ = 10.55 • , and comparison with thermal [P thermal (n) =n n r /(nr + 1) n+1 ] and coherent [P coherent (n) = e −nrnn r /n!] photon-number distributions. The error bar represents the 95% confidence interval. (d) Relative area for φ = 11.6 • . At this angle, the resonator has a smaller κ r , resulting in a larger steady-state photon number. The coherent-state distribution has better agreement.
than the linewidths (2Γ (0), 2Γ (1))/2π = (8.6, 11.6) MHz. The Hamiltonian model yields g s /2π = 33.4 MHz and ∆ s /2π = −102 MHz. The dressed linewidth κ r accounts for exact experimental conditions at the time of the measurement, including Purcell decay caused by the charge qubit that could otherwise change slightly over time and conditions. The area under the photon-number dips should be proportional to the probability of each photon number n [32,38]. The relative areas are plotted in Fig. 4c and compared with thermal-and coherent-state distributions. The coherent state withn r = 0.62 shows better agreement. This is consistent with the observation that the areas of the extra photon-number dips are reduced when using lower probe powers. Additional analysis can be found in Appendix G.
The ability to resolve quantized photon-number shifts in the qubit spectrum is a feature of the so-called strong dispersive regime of circuit QED, χ s > {Γ s , κ r }. It enables the preparation and detection of quantum photon states such as number states or cat states [33], and therefore paves the way towards bosonic codes [34,39] with spin qubits. Reciprocally, this also entails a shift of the resonator frequency larger than its linewidth, enabling fast and strong qubit readout. Finally, it dramatically highlights the consequences of residual photons on the qubits' dephasing.

IV. DISCUSSION
The strong backaction of the probe photons on the spin observed here highlights the limits of continuouswave measurements, and the necessity for future work to include time-domain control or dedicated readout and coupling resonators, as is now standard with supercon-ducting qubits. Since the probe power needs to be much belown r = 1 photons, this entails a low signal-to-noise ratio, especially without a parametric amplifier. Because of slow drift in the DQD interdot transition specific to this device (of the scale of minutes), long averaging times become problematic. However, this is not a fundamental problem for the platform in the future, as spin-qubit devices can be much more stable.
The fit to the formula for Γ (n) suggests a fundamental spin linewidth γ s /2π = 3.4 ± 1.0 MHz when subtracting photon losses (κ r /2π = 3.0 ± 0.2 MHz) and for g s /2π = 33.4 MHz. This linewidth value should be cited with care given that there are many assumptions involved; however, it compares favorably with literature values [22][23][24], and especially when considering the large g s achieved here. As with previous work, the spin linewidth seems limited by charge noise coupling in through the artificial spin-orbit interaction.
The large g s in this experiment is key to achieving coherent spin-photon interaction in the dispersive regime. This would not be possible without the large g c enabled by the high-impedance resonator, considering that g s ≤ g c [20,21]. Other contributing factors include the engineered spin-orbit interaction through the micromagnet's ∆B ⊥ , and state-of-the-art resonator losses despite the increased resonator coupling to its environment [36]. Interestingly, the charge-qubit linewidth Γ c /2π 60 MHz at 2t c /h = 12 GHz is not particularly small, suggesting that the good overall performance of the device comes from other factors, like the low g s /g c ratio [40], and could be improved further [41]. The cooperativity reaches a demonstrated value of (g s ) 2 /κ r Γ 1 s = 72 (taking a conservative Γ 1 s /2π = 5.2 MHz) or a projected (g s ) 2 /κ r γ s = 109 (assuming photon-induced broadening could be eliminated). Given that the exchange splitting is larger than the spin decoherence rates, a two-qubit gate of modest fidelity greater than or close to 75% [21] could potentially be achieved with the current device. To fully explore the optimal parameter space, a device with improved longterm stability, larger resonator coupling rates for readout, and tailored gate filters to allow driving signals without attenuation would all be beneficial.

V. CONCLUSION
In summary, we have demonstrated coherent hybridization of two spins mediated by virtual photons, as well as spin dispersive shifts by single photons, both larger than the spin linewidth. These experiments are more challenging than previous demonstrations of circuit QED with spins because they require a larger coupling-todecoherence ratio (i.e., cooperativity). Admittedly, the cooperativity in this platform is not yet on par with contemporary superconducting qubits. The regime of circuit QED achieved here is quite promising for the platform; it could enable two-qubit gates between spin qubits mediated by resonators [20,21], single-shot dispersive spin-qubit readout (without spin-to-charge conversion) [40], bosonic codes through preparation and detection of quantum photon states with spins [33,34,42], coherent links between dense spin-qubit networks [5], or quantum simulation with spin QED networks [43]. For future improvements, we believe that decoupling the spin from the charge noise is a promising path. An increased chargephoton coupling would allow to further detune the charge qubit from the spin and photon while preserving the spinphoton coupling, effectively suppressing the effects of the charge linewidth. This is because g s ∝ (2t c − ω r ) −1 while Γ s ∝ (2t c − ω r ) −2 is suppressed more quickly [20,21]. Furthermore, improved longitudinal magnetic gradient symmetry (we measured ∆B ∼ 1 mT) could help reduce the spin's noise sensitivity. Finally, improvements to materials and fabrication could help reduce the charge noise itself.

ACKNOWLEDGMENTS
The authors thank T. Bonsen and M. Russ for helpful insights involving input-output simulations, L. P. Kouwenhoven and his team for access to the NbTiN film deposition, F. Alanis Carrasco for assistance with sample fabrication, and other members of the spin-qubit team at QuTech for useful discussions.
Funding This research was funded in part by the European Research Council (ERC Synergy Quantum Computer Lab), the Dutch Ministry for Economic Affairs through the allowance for Topconsortia for Knowledge and Innovation (TKI), and the Netherlands Organization for Scientific Research (NWO/OCW) as part of the Frontiers of Nanoscience (NanoFront) program.
Author Competing interests The authors declare no competing interests.
Data availability The data reported in this paper are archived online at https://dx.doi.org/10.4121/ 15015453.

Appendix A: Device fabrication
The 10-nm-thick 28 Si/SiGe quantum-well heterostructure is grown on a 100-mm Si wafer via reduced-pressure chemical vapor deposition. The SiGe barrier thickness is 30 nm. Photolithography alignment markers are plasma etched into the surface with a Cl/HBr chemistry. Doped contacts to the quantum well are formed by 31 P implantation masked with photolithography and activated with a 700 • C rapid thermal anneal. The 5−7-nm superconducting NbTiN film is deposited via magnetron sputtering, preceded by a hydrofluoric acid dip and Marangoni drying, and followed by liftoff of the resist-covered quantum dot areas. The sheet inductance is targeted to be around 115 pH/ (measured at 140 pH/ for this device). The 10-nm Al 2 O 3 gate oxide is grown by atomic-layer deposition, followed by wet etching with buffered hydrofluoric acid everywhere except for the resist-covered quantum-dot areas. Contacts to implants, contacts to the NbTiN film, and electron-beam-lithography alignment markers are patterned with Ti/Pt evaporation preceded with buffered hydrofluoric acid dip and followed by liftoff. The wafer is diced into pieces for further electron-beam-lithography steps. The 25-nm Al gates are deposited via evaporation followed by liftoff. The NbTiN film is etched via SF 6 /He reactive ion etching to define the resonator, inductors, capacitors, and gate lines in a single electron beam lithography step, leaving a 40-nm step after the etch. The thin-film capacitor is patterned by first sputtering 30 nm of silicon nitride in a conformal deposition, and then evaporating 5 nm of Ti and 100 nm of Au in a directional deposition, allowing for a single patterning and liftoff step. The SiN z conformal deposition covers the 40-nm steps created during the etch of the NbTiN film. The micromagnets are patterned by first sputtering 30 nm of silicon nitride in a conformal deposition, and then evaporating 5 nm of Cr and 200 nm of Co in a directional deposition, allowing for a single patterning and liftoff step. Pieces are diced into individual 4-mm-by-2.8-mm device chips (Fig. 5) to be wire bonded onto a printed circuit board for cryogenic measurements.
To reduce resonator losses due to resistive currents in the gate structure, the gates are made of Al, and they maintain superconductivity up to in-plane magnetic fields of around 0.5 to 0.6 T, sufficient for our spin-qubit exper-  iments, while the NbTiN structures maintain low losses up to several tesla [37]. To mitigate the microwave losses through the gate fan-out lines [35], microwave low-pass filters are patterned on the gate fan-out lines using a combination of nanowire inductors and thin-film capacitors [36].

Appendix B: Experimental setup
The device is cooled using an Oxford Instruments Triton 400 dilution refrigerator with a base temperature of approximately 8 mK. The refrigerator is equipped with a (6, 1, 1)-T vector magnet. The equipment setup is shown in Fig. 6.
The resonator probe tone is generated using a mixer to allow for rapid sweeping of the probe frequency. The heterodyne detection is performed with an IQ mixer and at a variable intermediate frequency (IF) in the range [10,110] MHz. Combined with voltage ramps on the plunger gates or IQ modulation of the pump tone with an arbitrary waveform generator (AWG), this configuration allows for rapid acquisitions of two-dimensional data.
Appendix C: Data reconstruction along zero charge detuning The device used in this experiment has slow drift of the (0, 1)−(1, 0) interdot charge transitions. Hence, for every magnetic field setting (B r or φ), an automated fitting procedure is used to extract a data cut along = 0 µeV of a data frame and to reconstruct the two-dimensional data. This procedure ensures that the measurements are protected against long-term drift.
Each data frame is acquired in a single digitizer call, as in Fig. 7a. The fast axis is usually gate voltage, swept with an AWG ramp. The slow axis is usually a pump or a probe frequency, also controlled with AWG I or IQ modulation. The waveforms are repeated with a period of 4 ms (typical) and averaged into a frame consisting of 1000 to 3000 repetitions (4 to 12 s of cumulative integration time). Data transfer overheads mean that the frames take between 20 and 60 s to acquire and process.
The empirical fitting functions can be found in the source code files used for the data processing. The algorithm works by extracting a line cut at or near the bare resonator frequency, as in Fig. 7b. The maximum signal is always away from zero charge detuning, and it is symmetrically centered on zero detuning in most cases. For certain values of B r , the spin dispersive shift (positive) compensates the charge dispersive shift (negative) and can lead to features inside the dip, which can look asymmetric. Because these features are always lower in amplitude than the main edges that set the symmetry point, the algorithm is usually robust to this. In some cases, we also improve the fit by accounting for these artifacts, also empirically. The case presented in Fig. 7 is one of those more difficult situations, which works nonetheless (near the spin-photon resonance). The consistency of the whole procedure is validated by visual inspection of the zero-detuning fit results (e.g., by looking at plots like For data where both spins interact simultaneously with the resonator, one of the dots is fixed at zero charge detuning while the other dot detuning is rapidly swept, and a similar fitting procedure as described above is applied. This makes the acquisition tolerant of drift for the swept detuning, but the other dot must remain fixed near zero detuning long enough for the automated acquisition to be completed. Certain validation procedures can be applied to verify that the fixed dot did not drift after acquisition. For example, the dispersive shifts must add to (χ 1 c + χ 2 c )/2π; away from zero detuning, the shift is less. Conditions are sometimes set to recenter the detuning and reacquire the frame if some conditions are not met. These postvalidation procedures are used best with dispersive spin sensing, but they are harder to implement with vacuum Rabi splitting reconstructions (near spinphoton resonance). They are not always necessary, and are only applied when too much drift is observed.
Occasional bad data lines are caused by missed zero charge detunings and can easily be seen in the reconstructed data as vertical lines of poor or anomalous signal. Although they are admittedly not pretty, we are confident that they do not impact the validity or interpretation of the results.
The full source data and analysis code is available in the online repository.  energies between eigenstates: Here, h ( ) is the (reduced) Planck constant, a is the photon annihilation operator, ω r /2π is the resonator frequency, i is a DQD index, g i c /2π is the charge-photon coupling, τ i z = L i L i − R i R i and σ i are charge and spin Pauli operators (respectively), with σ i = σ i xx +σ i yŷ +σ i zẑ a vector of Pauli operators, and i /h and 2t i c /h are the DQD charge detuning and tunnel splitting (respectively). In cases where the DQD parameters are the same, or in the context of formulas that describe only one DQD, we omit the DQD index. The external and micromagnet magnetic fields are parametrized using the average and difference field energies h i and ∆h i at the left (L) and right (R) dot positions, with g e = 2 the electron Landé g-factor and µ B the Bohr magneton. Vectors like the external magnetic field B ext = (B r , θ = −90 • , φ) are conveniently expressed in spherical coordinates, with φ the polar angle. The micromagnet average field is modeled with the empirical formula to account for the susceptibility (χ difficulties in measuring its response, with ∆B ⊥ = 42 mT. We have not noticed discrepancies in g s (through χ s or J) that could be specifically attributed to this approximation. An effect that is left out is the fact that the Zeeman energy is not the same in the left and right dots. Experimentally, we observe ∆B ∼ 1 mT by tracking the spin-photon resonance condition versus DQD charge detuning. This is responsible for some asymmetry in plots like the one of Fig. 7; however, we find that taking this into account is not necessary to model the energy levels at zero detuning. The parametrization described in this paragraph is sufficient to capture the magnitude and angular dependence of the spin Zeeman energies over the range of interest. Some derived, effective, spin-photon quantities (such as g s , ∆ s = ω s − ω r , χ s , J) are used by analogy to their idealized versions, without the charge degree of freedom, with, for example, a Tavis-Cummings Hamiltonian [31] H TC = ω r a † a or a dispersive Hamiltonian [31] Here, ω r = ω r − χ s and ω s = ω s + χ s . The spin frequency ω s /2π and photon frequency ω r /2π depend strongly on whether each dot interacts with the resonator ( = 0) or not ( → ∞). The spin frequency is lowered when the interaction is on because of the artificial spin-orbit coupling [44], ω s → ω s + χ SO , for ∆h ⊥ h, 2t c − |h| |∆h| and |h| |∆h|. The photon frequency is lowered by the dispersive interaction with the charge, ω r → ω r − χ i c , for 2t c > ω r . We find that the charge dispersive shift is well predicted by taking into account the counter-rotating terms (Bloch-Siegert shift), Failing to do so can lead, for example, to an overestimated g c /2π = 220 MHz instead of 192 MHz. A schematic representation of various shifts in the dispersive regime is shown in Fig. 8a. While the analytical forms are insightful, we use exact numerical values calculated from Eq. (D1) instead in this work. To avoid ambiguities, we define ∆ i s /2π as the bare (i.e., theoretical but including χ c and χ SO ) spin-photon detuning for spin i individually interacting with the resonator (see Fig. 8a), and ∆ i 2s /2π by taking the individually interacting bare spin frequency and the simultaneously interacting bare photon frequency (which we find to be a quite accurate proxy; see φ = 13 • in Fig. 8). The dispersive approximations χ s ≈ (g s ) 2 /∆ s and J ≈ g 1 s g 2 s (1/∆ 1 2s + 1/∆ 2 2s )/2 are insightful but not quantitatively accurate because of various shifts in the frequencies, rotating wave approximations, dressing by the DQD charge degree of freedom, or violations of the dispersive approximation [requiring ∆ s g s orn r n crit = (∆ s /2g s ) 2 ]. It is not within the scope of this work to derive these quantities from system parameters, as this has been tackled elsewhere [16,17,20,21]. Instead, and in spite of Eqs. (D9) and (D10), we define 2g s /2π as the vacuum Rabi splitting gap, 2χ s /2π as the shift in the spin frequency induced by one photon, and 2J/2π as the spin-spin splitting; these quantities are extracted from the data or from Eqs. (D1) and (D2). This avoids the approximation pitfalls mentioned previously.

Parameter extraction from the data
Here, we describe the procedure to extract the model parameters from the experiment. Parameters are suc- q R e f A X D 2 t V 7 e s H R y 6 r f H 0 3 C s u z 0 + O z C F H d i Z u p X W f H 6 9 z B 6 / z M / E q 5 8 X 5 8 T p 2 8 D o + E 6 9 q R q 3 L r e 9 N r 1 K 8 T 7 B a j w q L H I q G O 2 r s W 6 N G K I q K G 2 F 4 U 8 w p V 7 O U u 0 H + S 6 1 5 I V + C e r R 7 3 F U r e 2 4 b I q G O 5 J 4 P d 6 w z 4 U 4 g j o o G Y 0 Y 7 y O 6 a 4 d i K N w 7 G 0 2 N B G G L f i t g P R q R 1 T B t P 7 e P Q u z r y F E L t P d c x w r x 5 z + r z n z O q b y Z 4 H k G e + Z Y t 9 p g 6 J U l U C q 2 q O c N + t c e 2 q N + X F v Q v m f 6 u g E l / i v Q X W C F z z 0 w t 5 f 3 e B l R 7 2 6 v 9 1 S q r q n m r n Z w 3 K E U l 9 v 4 i 9 L 9 K u r f e n a l 2 m t j k L T 6 + K p 3 3 8 r f z 0 6 j P b 3 w o y P i / P n d 9 t / / W 1 1 Y s n e 9 u 7 v 9 3 e / e b 6 + S 9 + L f 4 T p X n 6 T 5 w r s 0 s 8 n I 0 q m 9 4 j K y C 0 A 6 t e K + s e P s G 3 g       Here, (f 1 probe , f 2 probe , f 12 probe ) ≈ (6.9038, 6.9061, 6.8981) GHz for φ = 10.8 • . The spin transition frequency at = 0 µeV is marked with a white dashed line for comparison between plots. (d) When the spins are simultaneously interacting with the resonator, the spin states hybridize depending on their energy difference. In the nonresonant case, φ = 13.0 • , their energies are minimally perturbed, as can be seen by the model lines and the white guides. In the resonant case, φ = 10.8 • , the two spins hybridize due to the exchange interaction 2J/2π = 20.2 MHz mediated by virtual resonator photons. The upper state is hardly visible because its symmetry makes it dark to the resonator probe, which is expected for a coherent spin-spin interaction.

W Q t j B s X k d h b D Z n 7 S Z B Z 6 t J T 4 D a t p E u e U z N b m l Y 7 n f H 9 D f H z M o q 3 Y 8 6 E e V 7 Y 3 T K f R 3 M s h O 1 d I 5 Z y e R O R h c u t M n R 1 k a C C 9 I 8 V I F p 7 8 O X L w w t X P U 4 D d e O m j x R p N 6 g 6 8 4 B 9 T d f 7 G c y 6 J w B V E 7 0 Z + 4 d B p 2 q O w L A a W T 7 G Q T 1 p F v / 1 A J 0 D h G H P P Q M K 0 5 X v Z E H e s J 7 W 3 C O r e s E B F j 7 d f E O b 5 e 2 9 e r S T Q / q 2 O p 5 x 8 K 3 d u 3 E 7 o u a F 8 h 5 m / 2 J f e m x Z + h h 1 6 m i d Z 5 4 t c C z V + x q D d N V s 1 P t X u B p y D B s q v P Y p a U I p m p F 8 O 1 a K 6 q 5 d S G r S N v Y z j 2 2 U H + h N 8 x t c x a e 3 h B P +
cessively constrained using specific experiments for each DQD one at a time. First, the resonator bare frequency and bare linewidth are easily extracted from a probe frequency sweep while the dots are in Coulomb blockade. These quantities are later dressed by the interaction with charge and spin. For example, the resonator linewidth is significantly affected by Purcell decay from the charge qubit, and this practically limits the achievable spinphoton coupling and spin measurement sensitivity. The DQD lever arm is extracted from bias triangles. The charge-qubit transition frequency 2t c /h is measured with two-tone spectroscopy, and the corresponding resonator dispersive shift (away from spin-photon resonance) then allows one to uniquely calibrate the value of g c . Next, the micromagnet parameters are extracted by simultaneously adjusting the spin-photon transition frequencies to the experimental values of Fig. 9 (vacuum Rabi splitting, and two-tone spin measurement versus B r and φ) for each dot independently. Notably, the size of the spin-photon gap (2g s /2π) in a vacuum Rabi splitting measurement is most affected by ∆B ⊥ . The micromagnet susceptibility χ um and offset field B um0 are mainly fixed by the slope of the spin transition frequency versus B r and the spin-photon resonance condition, respectively. Then, B r0 is tweaked to get the best simultaneous agreement with both dots and for various angles.

Parameter table
A summary of the main model parameters for the key results is given in Tab. I.

Appendix E: Spin transitions as a function of DQD detuning
The DQD tunnel splitting is set to 2t c /h = 8.7 GHz, yielding g s /2π = 34.2 MHz. A pump tone at frequency f pump is sent down a gate line to each DQD to generate an excited spin-up population, while the transmission coefficient S21 is probed at a fixed frequency f probe set to the dispersively shifted resonator frequency for each DQD at zero charge detuning (e.g., as in the insets of Fig. 2). The results are plotted in Fig. 8 for two values of the magnetic field angle φ, corresponding to off-resonant and resonant spin transition energies. The spin transition frequency is reduced at zero charge detuning compared  Figure 9. Example of data used to extract model parameters. For this set, 2tc/h = 13.5 GHz, Br = 52 mT, and φ = 10.5 • . The first column shows vacuum Rabi splitting reconstructions for individually and simultaneously interacting spins. The second column shows dispersive spin sensing (the data near spin-photon resonance are scrambled because of the close proximity of the levels). The third column shows the angular dependence of the spin levels with dispersive spin sensing. The two spinlike and the photonlike transitions are shown together. These data, and more, are used to calibrate the micromagnet parameters for use in the Hamiltonian model Eq. (D2). For these tunnel splittings and spin-photon detunings, the model predicts gs/2π = 11.4 MHz and 2J/2π = 3.3 MHz. Here, 2J/2π is too small to be resolved. Note that transition linewidths are power broadened.
with the localized dot states because of the spin-charge hybridization shift χ SO /2π. The signal shows a peak centered on = 0 µeV when the resonator frequency gets pushed down by the charge dispersive shift χ c /2π. When the pump frequency matches the spin transition frequency, the excited spin-state population is increased, and the resonator frequency is pushed down further by the spin dispersive shift 2χ s /2π, visible as a dip in the transmission that follows the spin transition energy. When the two spins are not resonant (φ = 13.0 • ), each spin can be independently measured while interacting with the resonator one at a time. When the two spins are simultaneously interacting, they can both still be sensed, albeit with an adjusted probe frequency to account for the dual charge dispersive shift, and the spin transitions are mostly unperturbed from their independent values, as seen from the white dashed line serving as a guide. When the two spins are set to resonance (φ = 10.8 • at this magnetic field), the two states avoid each other when both spins are simultaneously interacting. Comparing the simultaneous spin-interaction results in the two cases, we see that the lower state has enhanced visibility while the upper state has a reduced one. This is consistent with the formation of a dark state, an effect that results Table I. Summary of model parameters. The resonator linewidth is a relaxation rate (κr = κ1), while the charge and spin (the "qubits") linewidths are dephasing rates, i.e., γ2 = γ1/2 + γ φ , as is standard in the field [31]. In general, the frequencies and shifts depend strongly on whether one, two, or zero of the DQDs are interacting with the resonator ( i = 0 µeV). from the symmetry of the hybridized spin states and is expected in the case of coherent spin-spin interaction [3]. As in the other cases, model transitions are adjusted to the individual interaction data and used to predict the simultaneous interaction data. The simultaneous interaction is again well predicted by the Hamiltonian model, as can be seen from the orange dashed lines. From the model, we extract a minimum separation between the spin states of 2J/2π = 20.2 MHz, with B r = 52 mT and ∆ 2s = −80 MHz.
Appendix F: Extended spin-spin coupling data In this section, we present additional data that demonstrate the hybridization of the two spin states. Hybridization as a function of spin-photon detuning (through B r ) is shown in Fig. 10. To speed up data acquisition, a larger probe power is used than in the main text data of Fig. 3, which, as a consequence, also significantly broadens the width of the spin transitions due to the photon-numberdependent spin dispersive shift. The apparent jitter of the transition energy as a function of φ is caused by jitter of the DQD2 tunnel coupling over the time it takes to reconstruct all angles. As expected, the spin splitting is reduced at larger |∆ 2s |. The upper branch is also dark in all plots near spin-spin resonance. These observations are consistent with spin-spin hybridization mediated by virtual resonator photons.
Next, hybridization of the two spin states is shown as a function of the spin-photon coupling (through t c ) in Fig. 11. In this work, we have measured spin-spin interactions (J) with values of spin-photon interactions g s /2π up to 40 MHz. A practical limit on how large g s can be  Figure 10. Extended spin-spin hybridization data versus spin-photon detuning. (a-e) Dispersive spin-spin hybridization measured for different spin-photon detuning values. The probe power is larger for these data than for the data in Fig. 3; the shadow visible below the main transition is caused by the photon-induced dispersive shift of the spin. The apparent jitter of the transition energy as a function of φ is caused by jitter of the DQD2 tunnel coupling over the time it takes to reconstruct all angles. The 2tc/h ≈ 8.7 GHz is retuned between plots, and yields gs/2π ≈ 34 MHz for both spins. The model fully captures the observed transition frequencies within plots and between the plots, and the only adjustment is the experimentally measured small variations in tc. (f) Splittings 2J/2π extracted from the model as a function of spin-photon detuning (through Br).
Comparison with the linewidth is impaired by the photon-number broadening of the transitions. As expected, 2J/2π decreases as |∆2s/2π| increases.
is the broadening of κ r at small charge-photon detunings, mainly caused by Purcell decay from the charge qubit. At some point, the resonator becomes too undercoupled and the signal too small. This could therefore be improved by using a dedicated readout resonator with an optimized coupling, by adding a near-quantum-limited amplifier at the mixing chamber, or by reducing the charge linewidth. The setting used in the main text data Fig. 3 is chosen empirically based on linewidth, exchange coupling and readout signal.
In the experiment of the main text Fig. 3, the critical photon number n crit = (∆ 2s /2g s ) 2 = 1.5 is quite low. This number is often used to quantify the validity of the dispersive regime, which requires the expectation value of the photon numbern r to ben r n crit [31]. We use the photon-number-dependent spin dispersive shift to establish an upper bound ofn r < 0.5 for the data of the main text Fig. 3 (n r ≈ 0.4 ± 0.1). We can no longer distinguish between the coherent and thermal distributions since they converge and the data are too noisy. The n crit could be optimized for future qubit experiments. For instance, in Fig. 10b, the n crit = 2 is already much higher while 2J/2π is still 17 MHz. Timedomain control can also help reducen r by probing only during readout. The exact requirements will depend on the target two-qubit gate and desired fidelity (amongst other things).
Appendix G: Extended photon-number-dependent spin dispersive shift data In this section, we extend the analysis of the photonnumber-dependent spin dispersive shift data by looking at data for φ = 11.6 • . In Fig. 12b, a line cut from Fig. 12a is extracted, and the dip areas and separations are fit to a sum of Lorentzian dips with 2Γ (n) = 2γ s + (n +n r )κ r . A value of 2χ s /2π = −12.5 MHz is extracted from the fit (g s /2π = 33.4 MHz, ∆ s /2π = −122 MHz), slightly larger than the average of the linewidths (2Γ (0), 2Γ (1))/2π = (10.1, 13.1) MHz. The coherent state has a higher photon number,n r = 0.82, than at φ = 10.55 • (see main text, n r = 0.62). This can be attributed to the resonator's higher dressed quality factor for φ > 11.2 • (which is independently verified), and it leads to a larger probe photon population at a steady state. This is also consistent with the background showing a larger relative |S21| value at fixed probe power. The difference between the coherentstate and thermal-state distributions is more pronounced than in the main text. Figure 11. Extended spin-spin hybridization data versus spin-photon coupling strength. (a-d) Dispersive spin-spin hybridization measured for decreasing spin-photon coupling values. The gs is adjusted through its dependence on tc and is approximately the same for both spins, while Br = 52 mT is fixed. The ∆2s/2π ≈ 78 MHz is approximately unchanged between plots because changes in χSO are approximately compensated by changes in χc. The model fully captures the observed transition frequencies within plots and between the plots, and the only adjustment is the experimentally measured tc. As expected, 2J/2π decreases as gs/2π decreases. In Tab. I, a small discrepancy is observed between the predicted and measured values of χ s . This could be due to driven system dynamics, which are known to modify the splitting [32].