Patrick Harvey-Collard
[Correspondence to<EMAIL_ADDRESS>QuTech and Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands

Jurgen Dijkema
QuTech and Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands

Guoji Zheng
QuTech and Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands

Amir Sammak
QuTech and Netherlands Organization for Applied Scientific Research (TNO), 2628 CJ Delft, The Netherlands

Giordano Scappucci
QuTech and Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands

Lieven M. K. Vandersypen
[Correspondence to<EMAIL_ADDRESS>QuTech and Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands

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.

§ 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 high-quality-factor superconducting resonators at cryogenic temperatures to enable the coupling and readout of superconducting qubits <cit.>.
Meanwhile, spin qubits in gate-defined semiconductor quantum dots (QDs) are also promising for quantum computing <cit.>, 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.

Significant research efforts have been dedicated to bringing the benefits of circuit QED to the platform of spin qubits <cit.>, in order to leverage different characteristics of each system. Experiments have first incoherently coupled photons to QD charges or spins in various materials <cit.>.
Following multiple theory proposals <cit.>,
experiments using spin-charge hybridization have reached the strong coupling regime of circuit QED with single electron spins in silicon <cit.>, with multispin qubits in GaAs <cit.>, and with carbon nanotubes <cit.> (albeit not resonantly).
Concerning distant interactions, many works have demonstrated the coupling of combinations of gallium-arsenide qubits and transmons qubits <cit.>.
Recently, resonant spin-spin-resonator coupling has been demonstrated in silicon <cit.>. However, the resonant spin-spin-resonator regime does not allow a straightforward two-qubit gate <cit.>.

Despite these groundbreaking realizations, some important hallmarks of circuit QED experiments have remained elusive for spins, in part due to the insufficient spin-resonator interaction strength, in combination with decoherence and fabrication challenges. These hallmarks include dispersive interaction between two spins mediated by virtual resonator photons <cit.> and photon-number-dependent spin dispersive shifts <cit.>, both requiring a higher level of interaction-to-decoherence ratio than previously achieved <cit.>. 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 <cit.> and universal control <cit.>.

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 high-impedance superconducting resonator. Each single spin is trapped in a double quantum dot (DQD) formed in a  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 <cit.>. Finally, we resolve the photon-number states from the discrete shifts they induce in the spin transition frequency.

§ METHODS

The device is shown in fig:device.

Spin-spin coupling device.    a Optical image of the resonator and QD areas with colored overlays. The ground plane, resonator, microwave (MW) in and out ports, and DC tap are thin superconducting . The substantial kinetic inductance and resonator impedance result in a $6.916$-GHz resonator that is only $250 \um$ long.    b Colorized angled scanning electron microscope image of a (nominally identical) DQD gate structure and micromagnets. The  resonator end is contacted with an Al plunger electrode of the DQD, which is biased through the DC tap. (Image rotated $90 \deg$.)    c Schematic of the DQD heterostructure and gate stack using the same color scheme as in panel (b). A photon couples to an electron spin through spin-orbit (micromagnet) and orbital-photon (voltage zero-point fluctuations $V_\text{zpf}$) interactions.       d Micromagnets are tilted $\pm 15 \deg$, allowing one to fine-tune each spin's Zeeman energy through the external magnetic field angle $\phi$. Here, $\vct B_\text{ext} = (B_r, \theta=-90\deg, \phi)$.

A 5-to-7-nm-thin film of  is deposited on the surface of a  heterostructure and patterned to form the superconducting resonator, ground planes, gate filters <cit.>, and gate fan-out lines (). The sheet kinetic inductance (around $140 \pHpsq$), the narrow width of the resonator center conductor ($170 \nm$), and the retracted ground planes combine into a high effective resonator impedance, $Z_r \approx 3.0 \kohm$ <cit.>; this results in a fixed coupling $g_c/2\pi = 192 \MHz$ to the DQD charge degree of freedom (“charge qubit”) at a resonator frequency of $\omega_r/2\pi = 6.916 \GHz$ for the half-wave mode (fig:devicea). Achieving a large $g_c \propto \alpha_c V_\text{zpf}$ (approximately 5 times larger than borjans2020) is a result mostly of the increased voltage zero-point fluctuations $V_\text{zpf} \propto \omega_r\sqrt{Z_r}$ of the high-impedance resonator, and of the $(0,1){-}(1,0)$ interdot transition lever arm $\alpha_c$ <cit.>.
A combination of loss mitigation strategies () results in an undercoupled resonator with a linewidth $\kappa_r/2\pi = 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 \um$ long, which is a consequence of the substantial kinetic inductance that translates to a high effective magnetic permeability <cit.>.
The DQD potential is shaped by applying suitable voltages to surface gate electrodes, as shown in fig:deviceb-c.
Cobalt micromagnets provide a transverse magnetic field difference $\Delta B_\perp = 42 \mT$ between the two dots while minimizing magnitude differences $\Delta B_\parallel$ ($\Delta B_\perp$ is taken as constant, see  for details).
Unlike $g_c$, which is, to a large extent, fixed by the device structure, the spin-photon coupling $g_s/2\pi$ at zero charge detuning approximatively scales as

g_s ≈g_c g_eμ_BΔB_⊥/4(2t_c-ħω_r)

for $2t_c \gg \{\hbar\omega_s, \hbar\omega_r\}$, and it is therefore tunable via the DQD tunnel splitting $2t_c$ <cit.>. Here, $g_\text{e} = 2$ is the electron Landé $g$-factor, $\mu_\text{B}$ is the Bohr magneton, and $\omega_s/2\pi$ is the spin transition frequency. In any case, $g_s \leq g_c$.
These device characteristics combine to enable values of $g_s/2\pi$ up to around $40 \MHz$ in this work.
The micromagnets are tilted $\pm 15 \deg$ relative to the vertical direction (fig:deviced), allowing one to fine-tune each spin's Zeeman energy by rotating the in-plane external magnetic field of magnitude $B_r = \abs{\vct B_\text{ext}}$ by an angle $\phi$ <cit.> using a vector magnet.
Device fabrication and experimental setup details are given in  and , respectively.

§ RESULTS

First, the system is tuned to the regime where both spins are resonant with the photons <cit.>, which happens at a field angle of $\phi = 10.5 \deg$ due to a $6.5$-mT difference in the micromagnet fields at $\phi = 0 \deg$.

Resonant spin-photon-spin coupling.     a-b Spin-photon vacuum Rabi splitting of DQD1 (a) and DQD2 (b) for tunnel splittings $(2t_c^1, 2t_c^2)/h = (13.2, 13.7) \GHz$ and a field angle of $\phi = 10.5 \deg$. When the gate voltage is set on the zero-charge-detuning point, which corresponds to the middle of the $(0,1){-}(1,0)$ interdot charge transition, the charge-photon interaction $g_c/2\pi = 192 \MHz$ between the DQD charge and the resonator dispersively shifts the resonator frequency (insets) and enables the spin-photon interaction through the artificial spin-orbit interaction. The dashed lines are the Hamiltonian model transitions. The probe power is $-117.5 \dBm$.     c When both dots interact simultaneously with the resonator, the dispersive shifts are additive, and the vacuum Rabi splitting is enhanced by a factor of approximately $\sqrt{2}$, an effect of the coherent collective spin-photon interaction. Model transitions are adjusted to the individual interaction data and used to predict the simultaneous interaction data.

In the main panels of Figs. <ref>a and <ref>b, DQD1 and DQD2 are set separately to zero charge detuning, $\eps = 0 \ueV$, 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 spin-photon coupling) can be switched off simply by biasing the DQD charge detuning ($\epsilon \gg 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 $\phi$), an automated fitting procedure is used to extract a data cut along $\eps = 0 \ueV$ and reconstruct the two-dimensional data. This procedure ensures that the measurements are protected against long-term drift, and it is further detailed in .
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 $\omega_r/2\pi$, the charge-photon coupling $g_c^i/2\pi$, and the DQD charge detuning and tunnel coupling $\eps^i/h$ and $t_c^i/h$ (respectively). Here, $h$ is the Planck constant, and $i$ is an index identifying the DQD. Full details about the model are given in . 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:resonanta and fig:resonantb, we extract a spin-photon coupling of $(g_s^1, g_s^2)/2\pi = (11.8, 11.0) \MHz$ for DQD1 and DQD2, respectively, with tunnel splittings $(2t_c^1, 2t_c^2)/h = (13.2, 13.7) \GHz$. Both DQDs achieve the strong spin-photon coupling regime, i.e., $g_s > \{\kappa'_r, \varGamma_s\}$, where $\varGamma_s \leq 6 \MHz$ is the spin linewidth.
When both DQDs are set to zero charge detuning simultaneously (fig:resonantc), 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_s^{12}/2\pi = 32.3 \MHz$ from the two-spin ensemble. The enlarged splitting matches well the predicted $2\sqrt{(g_{s}^1)^2 + (g_{s}^2)^2}/2\pi$ value, very close to a factor $\sqrt{2}$ larger than for single spins in this case. The Hamiltonian model transitions, which are calibrated solely on the one-at-a-time interaction data, predict very well the outcome of the simultaneous interaction. The $\sqrt{2}$ enhancement indicates simultaneous, coherent, and resonant interaction of both spins with the resonator, as demonstrated in prior work <cit.>.

Spin-spin coupling via virtual photons.     a-b Dependence of the independently interacting spin transition frequencies as a function of the external magnetic field angle $\phi$ at $B_r = 52 \mT$. See  for the reconstruction procedure responsible for the vertical stripes.     c When both spins interact simultaneously with the resonator, an avoided crossing is observed. The spin states hybridize with a splitting $2J/2\pi = 19.0 \MHz$, larger than any of the transition widths $2\varGamma_s/2\pi \leq (13.2\pm4.2) \MHz$. A dark state is observed along the upper branch, an effect caused by the symmetry of the coherently hybridized spin states.     d Plot showing how the exchange interaction is reduced when the spin-photon detuning is increased.

To resolve the spin-spin exchange splitting $2J/2\pi$ mediated by virtual photons, a larger spin-photon coupling is needed than in the resonant case (see, e.g., ). Given a fixed ratio $g_s/\varDelta_s$, where $\varDelta_s = \omega_s - \omega_r$ is the spin-photon detuning, and given that $J \approx (g_s)^2/\varDelta_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 \approx 8.8 \GHz$, yielding $g_s/2\pi \approx 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 $\phi$ using two-tone spectroscopy. A pump tone at frequency $f_\text{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_\text{probe}$ set to the dispersively shifted resonator frequency for each DQD at zero charge detuning (e.g., as in the insets of fig:resonant).
Line cuts along the charge zero-detuning axis of data like those of  are assembled in a two-dimensional diagram, resulting in fig:dispersivephi.
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 $\phi$ 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 $\phi \in [9, 13] \deg$ is almost aligned with the DQD2 micromagnet angle $15 \deg$ and farther from the one of DQD1. When both spins interact simultaneously (fig:dispersivephic), an avoided crossing is observed, while the upper transition becomes dark close to spin-spin resonance <cit.>. With $(2t_c^1, 2t_c^2)/h = (8.82, 8.80) \GHz$ and $B_r = 52 \mT$, we extract $(g_s^1, g_s^2)/2\pi = (32.4, 32.7) \MHz$, $2J/2\pi = 19.0 \MHz$, and $\varDelta_{2s}/2\pi = -79 \MHz$ (here, $\varDelta_{2s}$ is for two spins since each DQD contributes an additive charge dispersive shift $\chi_c$ to $\w_r$, as explained in ). In fig:dispersivephid and , additional spin-spin hybridization results show that the exchange interaction is reduced when either $\abs{\varDelta_s}$ is increased or $g_s$ is decreased, as expected. The full width at half minimum of the dip, $2\varGamma_s/2\pi$, is $2\varGamma_s^-/2\pi = (11.7 \pm 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\varGamma_s^1, 2\varGamma_s^2)/2\pi = (10.4\pm3.0, 13.0\pm3.2) \MHz$ for the individual spin-interaction case on resonance ($\phi = 10.5 \deg$), and $(2\varGamma_s^1, 2\varGamma_s^2)/2\pi = (8.8\pm4.1, 13.2\pm4.2) \MHz$ for the simultaneous interaction away from resonance ($\phi = 13 \deg$). Since $2J > (2\varGamma_s^1 + 2\varGamma_s^2)/2$, the two spins are coherently hybridized through virtual photons, achieving a long-standing goal for the field.
Arguably, the data presented in this figure are not very deep into the dispersive regime. In , we present additional data at larger $\varDelta_s$ where $2J$ is still larger than $2\varGamma_s$.
The ratio of interaction strength to decoherence should be sufficient to enable two-qubit gates <cit.> 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 <cit.>. The shift $2\chi_s/2\pi$ is expected to scale as $\chi_s \approx (g_s)^2/\varDelta_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 $\varGamma(n)$ formula below).
The photon-number-dependent dispersive shift of DQD1 is shown in fig:dispersivephotonnumber.

Photon-number-resolved spin dispersive shift.     a Spin transition frequency of DQD1 showing an extra dip below the main transition for larger probe power. The dip matches well the prediction from the Hamiltonian model (dashed lines) for $\ket{\downarrow,n} \leftrightarrow \ket{\uparrow,n}$ transitions, with the first two transitions visible in this plot.     b Line cut of the data in panel (a) for $\phi = 10.55 \deg$, and fit to four Lorentzian dips (first three shown with dashed lines). A dispersive shift $2\chi_s/2\pi = -13.1 \pm 2.2 \MHz$ is extracted.     c Relative area of the dips for $\phi = 10.55 \deg$, and comparison with thermal [$P_\text{thermal}(n) = \bar n_r^n/(\bar n_r+1)^{n+1}$] and coherent [$P_\text{coherent}(n) = \E{-\bar n_r}\bar n_r^n/n!$] photon-number distributions. The error bar represents the $95 \pc$ confidence interval.     d Relative area for $\phi = 11.6 \deg$. At this angle, the resonator has a smaller $\kappa'_r$, resulting in a larger steady-state photon number. The coherent-state distribution has better agreement.

A larger probe power is used ($-119 \dBm$, $0 \MHz$ detuning) than for the data of fig:dispersivephi ($-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 $\ket{\downarrow,1} \leftrightarrow \ket{\uparrow,1}$ spinlike transition.
In fig:dispersivephotonnumberb, a line cut from fig:dispersivephotonnumbera is extracted, and the dip areas and separations are fit to a sum of Lorentzian dips with widths $2\varGamma(n)/2\pi = 2\gamma_s/2\pi + (n + \bar n_r)\kappa'_r/2\pi$. A value of $2\chi_s/2\pi = -13.1 \pm 2.2 \MHz$ is extracted from the fit (with parameters $\kappa'_r/2\pi = 3.0 \pm 0.2 \MHz$, $\bar n_r = 0.62$), slightly larger than the linewidths $(2\varGamma(0), 2\varGamma(1))/2\pi = (8.6, 11.6) \MHz$. The Hamiltonian model yields $g_s/2\pi = 33.4 \MHz$ and $\varDelta_s/2\pi = -102 \MHz$. The dressed linewidth $\kappa'_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$ <cit.>. The relative areas are plotted in fig:dispersivephotonnumberc and compared with thermal- and coherent-state distributions. The coherent state with $\bar n_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 .

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, $\chi_s > \{\varGamma_s, \kappa'_r\}$. It enables the preparation and detection of quantum photon states such as number states or cat states <cit.>, and therefore paves the way towards bosonic codes <cit.> 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.

§ DISCUSSION

The strong backaction of the probe photons on the spin observed here highlights the limits of continuous-wave measurements, and the necessity for future work to include time-domain control or dedicated readout and coupling resonators, as is now standard with superconducting qubits.
Since the probe power needs to be much below $\bar n_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 $\varGamma(n)$ suggests a fundamental spin linewidth $\gamma_s/2\pi = 3.4 \pm 1.0 \MHz$ when subtracting photon losses ($\kappa'_r/2\pi = 3.0 \pm 0.2 \MHz$) and for $g_s/2\pi = 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 <cit.>, 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 \leq g_c$ <cit.>. Other contributing factors include the engineered spin-orbit interaction through the micromagnet's $\Delta B_\perp$, and state-of-the-art resonator losses despite the increased resonator coupling to its environment <cit.>. Interestingly, the charge-qubit linewidth $\varGamma_c/2\pi \lesssim 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 <cit.>, and could be improved further <cit.>.
The cooperativity reaches a demonstrated value of $(g_s)^2/\kappa'_r\varGamma_s^1 = 72$ (taking a conservative $\varGamma_s^1/2\pi = 5.2 \MHz$) or a projected $(g_s)^2/\kappa'_r\gamma_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\pc$ <cit.> could potentially be achieved with the current device. To fully explore the optimal parameter space, a device with improved long-term stability, larger resonator coupling rates for readout, and tailored gate filters to allow driving signals without attenuation would all be beneficial.

§ 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-to-decoherence 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 <cit.>, single-shot dispersive spin-qubit readout (without spin-to-charge conversion) <cit.>, bosonic codes through preparation and detection of quantum photon states with spins <cit.>, coherent links between dense spin-qubit networks <cit.>, or quantum simulation with spin QED networks <cit.>.
For future improvements, we believe that decoupling the spin from the charge noise is a promising path. An increased charge-photon coupling would allow to further detune the charge qubit from the spin and photon while preserving the spin-photon coupling, effectively suppressing the effects of the charge linewidth. This is because $g_s \propto (2t_c-\hbar\omega_r)^{-1}$ while $\varGamma_s \propto (2t_c-\hbar\omega_r)^{-2}$ is suppressed more quickly <cit.>. Furthermore, improved longitudinal magnetic gradient symmetry (we measured $\Delta B_\parallel \sim 1 \mT$) could help reduce the spin's noise sensitivity. Finally, improvements to materials and fabrication could help reduce the charge noise itself.

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  film deposition, F. Alanis Carrasco for assistance with sample fabrication, and other members of the spin-qubit team at QuTech for useful discussions.

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 contributions

L.M.K.V.,  and G.Z. conceived and planned the experiments.
, J.D. and G.Z. performed the electrical cryogenic measurements.
fabricated the device.
A.S. contributed to sample fabrication.
A.S. grew the heterostructure with G.S.’s supervision.
, J.D. and L.M.K.V. analyzed the results.
wrote the manuscript with input from all co-authors.
L.M.K.V. supervised the project.

*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>.

§ DEVICE FABRICATION

The 10-nm-thick  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 ^31P implantation masked with photolithography and activated with a $700 \degC$ rapid thermal anneal.
The $5{-}7$-nm superconducting  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 \pHpsq$ (measured at $140 \pHpsq$ for this device).
The $10$-nm Al2O3 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  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  film is etched via SF6/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  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:extdevice) to be wire bonded onto a printed circuit board for cryogenic measurements.

Composite optical image of the device chip without wire bonds.

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 \tesla$, sufficient for our spin-qubit experiments, while the  structures maintain low losses up to several tesla <cit.>. To mitigate the microwave losses through the gate fan-out lines <cit.>, microwave low-pass filters are patterned on the gate fan-out lines using a combination of nanowire inductors and thin-film capacitors <cit.>.

§ 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:setup.

Schematic of the measurement setup.

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.

§ 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 $\phi$), an automated fitting procedure is used to extract a data cut along $\eps = 0 \ueV$ 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:fitcentera. 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.

Fitting of zero charge detuning.     a Example of a single frame of data used in a vacuum Rabi splitting reconstruction.     b Cut along the black line from panel (a) corresponding to the bare resonator frequency to extract the zero-detuning gate voltage $x_0$. An empiric fitting function is used to find the center of the signal dip due to the dispersive shifts.     c Cut along the fitted zero-detuning value $x_0$ (white line) and fit to a single Fano resonance.

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:fitcenterb. 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:fitcenter 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 fig:fitcentera one by one).

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 $(\chi_c^1+\chi_c^2)/2\pi$; 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 spin-photon 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.

§ SPIN-PHOTON HAMILTONIAN MODEL

§.§ Model definition

The spin-photon system is modeled with the following Hamiltonian to numerically calculate the exact transition energies between eigenstates:

\begin{align}
H &= \hbar\w_r{a^\dag}a + \sum_{i=1}^2 \ofb{H_i + \hbar g_c^i(a+a^\dag)\tau_z^i} , \label{eq:suph}\\
H_i &= \frac{1}{2}\ofb{ \eps^i\tau_z^i + 2t_c^i\tau_x^i + \of{\vct h^i + \vct{\Delta h}^i\tau_z^i/2}\cdot\vct\sigma^i } .
\label{eq:supmodel}
\end{align}

Here, $h$ ($\hbar$) is the (reduced) Planck constant, $a$ is the photon annihilation operator, $\omega_r/2\pi$ is the resonator frequency, $i$ is a DQD index, $g_c^i/2\pi$ is the charge-photon coupling, $\tau_z^i = \ketbra{L^i}{L^i}-\ketbra{R^i}{R^i}$ and $\vct\sigma^i$ are charge and spin Pauli operators (respectively), with $\vct\sigma^i = \sigma_x^i \uvct{x} + \sigma_y^i \uvct{y} + \sigma_z^i \uvct{z}$ a vector of Pauli operators, and $\eps^i/h$ and $2t_c^i/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 $\vct h^i$ and $\vct{\Delta h}^i$ at the left ($L$) and right ($R$) dot positions,

\begin{align}
\vct h^i/g_\text{e}\mu_\text{B} &= (\vct B^i_L + \vct B^i_R)/2 \\
&= \vct B_\text{ext} + \vct B^i_\text{um} , \\
\vct{\Delta h}^i/g_\text{e}\mu_\text{B} &= \vct B^i_L - \vct B^i_R \\
&= \vct{\Delta B}_\text{um}^i ,
\end{align}

with $g_\text{e} = 2$ the electron Landé $g$-factor and $\mu_\text{B}$ the Bohr magneton. Vectors like the external magnetic field $\vct B_\text{ext} = (B_r, \theta=-90\deg, \phi)$ are conveniently expressed in spherical coordinates, with $\phi$ the polar angle. The micromagnet average field is modeled with the empirical formula

\begin{align}
\vct B_\text{um}^i =
(B_\text{um0}^i + \chi_\text{um}^i(\vct B_\text{ext}-\vct B_\text{ext0}^i) \cdot \uvct u_\text{um}^i) \uvct u_\text{um}^i
\end{align}

to account for the susceptibility ($\chi_\text{um}^1 \approx 0.67$, $\chi_\text{um}^2 \approx 0.63$), and with $\uvct u_\text{um}^1 = (1, 90\deg, 15\deg)$ and $\uvct u_\text{um}^2 = (1, 270\deg, 15\deg)$. Empirical parameters are $B_\text{um0}^1 = 152.6 \mT$, $B_\text{um0}^2 = 146.9 \mT$, and $\vct B_\text{ext0}^i = (B_\text{r0}=-20\mT, -90\deg, \phi) \parallel \vct B_\text{ext}$. However, the micromagnet difference field is taken as constant because of practical difficulties in measuring its response,

\begin{align}
\vct{\Delta B}_\text{um}^i = (\Delta B_\perp, 180\deg, 90\deg) = -\Delta B_\perp \uvct x ,
\end{align}

with $\Delta B_\perp = 42 \mT$. We have not noticed discrepancies in $g_s$ (through $\chi_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 $\Delta B_\parallel \sim 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:fitcenter; 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.

Dispersive regime transitions, and spin-spin interaction versus DQD detuning.     a Schematic representation of the transition energy shifts in the dispersive regime. See sec:model for details.     b-c Measurement of spin transitions for DQD1 (b) and DQD2 (c) as a function of the charge detuning $\epsilon$ with dispersive readout, with $2t_c/h \approx 8.7 \GHz$ and $B_r = 52 \mT$. A pump tone is sent down the gate lines to excite the spins, resulting in a signal dip when the pump and spin frequencies match. Here, $(f_\text{probe}^1, f_\text{probe}^2, f_\text{probe}^{12}) \approx (6.9038, 6.9061, 6.8981) \GHz$ for $\phi = 10.8 \deg$. The spin transition frequency at $\epsilon = 0 \ueV$ 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, $\phi = 13.0 \deg$, their energies are minimally perturbed, as can be seen by the model lines and the white guides. In the resonant case, $\phi = 10.8 \deg$, the two spins hybridize due to the exchange interaction $2J/2\pi = 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.

Some derived, effective, spin-photon quantities (such as $g_s$, $\varDelta_s = \w_s-\w_r$, $\chi_s$, $J$) are used by analogy to their idealized versions, without the charge degree of freedom, with, for example, a Tavis-Cummings Hamiltonian <cit.>

\begin{align}
\begin{split}
H_\text{TC} &=  \hbar\w_r{a^\dag}a \\
&\phantom{= \ }
+ \sum_{i=1}^2 \ofb{ \frac{\hbar}{2}\w_s^i \sigma_z^i + \hbar g_s^i\of{a^\dagger \sigma_-^i + a \sigma_+^i}}
\end{split} \label{eq:hres}
\end{align}

or a dispersive Hamiltonian <cit.>

\begin{align}
\begin{split}
H_\text{disp} &=  \hbar\w_r'{a^\dag}a + \sum_{i=1}^2 \ofb{ \frac{\hbar}{2}\of{\w_s'^i + 2\chi_s^ia^\dagger a}\sigma_z^i } \\
&\phantom{= \ }
+ \hbar J\of{\sigma_+^1 \sigma_-^2 + \sigma_-^1 \sigma_+^2} .
\end{split} \label{eq:hdisp}
\end{align}

Here, $\w_r' = \w_r - \chi_s$ and $\w_s' = \w_s + \chi_s$.
The spin frequency $\w_s/2\pi$ and photon frequency $\w_r/2\pi$ depend strongly on whether each dot interacts with the resonator ($\eps = 0$) or not ($\eps \rightarrow \infty$). The spin frequency is lowered when the interaction is on because of the artificial spin-orbit coupling <cit.>, $\w_s \rightarrow \w_s + \chi_\text{SO}$,

ħχ_SO    = -Δh^2/8h
- 1/16 Δh^2/2t_c-h
+ 1/16 Δh^2/2t_c+h
< 0 ,

for $\vct{\Delta h} \perp \vct{h}$, $2t_c-\abs{\vct{h}} \gg \abs{\vct{\Delta h}}$ and $\abs{\vct{h}} \gg \abs{\vct{\Delta h}}$.
The photon frequency is lowered by the dispersive interaction with the charge, $\w_r \rightarrow \w_r - \chi_c^i$, for $2t_c > \hbar \w_r$. We find that the charge dispersive shift is well predicted by taking into account the counter-rotating terms (Bloch-Siegert shift),

χ_c = (g_c)^2/2t_c/ħ-_r + (g_c)^2/2t_c/ħ+_r > 0 .

Failing to do so can lead, for example, to an overestimated $g_c/2\pi = 220 \MHz$ instead of $192 \MHz$. A schematic representation of various shifts in the dispersive regime is shown in fig:dipersivedetuninga. While the analytical forms are insightful, we use exact numerical values calculated from eq:suph instead in this work.
To avoid ambiguities, we define $\varDelta_s^i/2\pi$ as the bare (i.e., theoretical but including $\chi_c$ and $\chi_\text{SO}$) spin-photon detuning for spin $i$ individually interacting with the resonator (see fig:dipersivedetuninga), and $\varDelta_{2s}^i/2\pi$ 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 $\phi = 13 \deg$ in fig:dipersivedetuning).
The dispersive approximations $\chi_s \approx (g_s)^2/\varDelta_s$ and $J \approx g_s^1g_s^2(1/\varDelta_{2s}^1+1/\varDelta_{2s}^2)/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 $\varDelta_s \gg g_s$ or $\bar n_r \ll n_\text{crit} = (\varDelta_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 <cit.>.
Instead, and in spite of Eqs. (<ref>) and (<ref>), we define $2g_s/2\pi$ as the vacuum Rabi splitting gap, $2\chi_s/2\pi$ as the shift in the spin frequency induced by one photon, and $2J/2\pi$ as the spin-spin splitting; these quantities are extracted from the data or from Eqs. (<ref>) and (<ref>). This avoids the approximation pitfalls mentioned previously.

§.§ Parameter extraction from the data

Example of data used to extract model parameters.     For this set, $2t_c/h = 13.5 \GHz$, $B_r = 52 \mT$, and $\phi = 10.5 \deg$. 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:supmodel. For these tunnel splittings and spin-photon detunings, the model predicts $g_s/2\pi = 11.4 \MHz$ and $2J/2\pi = 3.3 \MHz$. Here, $2J/2\pi$ is too small to be resolved. Note that transition linewidths are power broadened.

Summary of model parameters. The resonator linewidth is a relaxation rate ($\kappa_r=\kappa_1$), while the charge and spin (the “qubits”) linewidths are dephasing rates, i.e., $\gamma_2=\gamma_1/2+\gamma_\phi$, as is standard in the field <cit.>. In general, the frequencies and shifts depend strongly on whether one, two, or zero of the DQDs are interacting with the resonator ($\epsilon^i = 0 \ueV$).

Quantity, symbol     units     DQD1     DQD2     resonator

noninteracting resonator frequency, $\omega_r/2\pi$			   GHz         6.916

noninteracting resonator linewidth, $\kappa_r/2\pi$				   MHz         1.8

charge-photon coupling, $g_c/2\pi$					   MHz   192   192

charge-qubit linewidth ($2t_c/h=12\GHz$), $\varGamma_c/2\pi$		   MHz   $\lesssim60$   $\lesssim60$

transverse magnetic field difference, $\Delta B_\perp$	   mT   42   42

5lResonant coupling,

tunnel splitting, $2t_c/h$							   GHz   13.2   13.7

charge dispersive shift, $\chi_c/2\pi$					   MHz   7.689   7.214

spin-orbit shift, $\chi_\text{SO}/2\pi$					   MHz   -34.4   -33.6

spin-photon coupling, $g_s/2\pi$					   MHz   11.8   11.0

half vacuum Rabi splitting (two spins), $g_s^{12}/2\pi$	   MHz         16.2

spin linewidth, $\varGamma_s/2\pi$					   MHz   $\leq6$   $\leq6$

interacting resonator linewidth, $\kappa'_r/2\pi$			   MHz   2.6   2.3

5lDispersive coupling,

tunnel splitting, $2t_c/h$							   GHz   8.82   8.80

charge dispersive shift, $\chi_c/2\pi$					   MHz   21.46   21.66

spin-orbit shift, $\chi_\text{SO}/2\pi$					   MHz   -62.3   -62.6

spin-photon coupling, $g_s/2\pi$					   MHz   32.4   32.7

spin-photon detuning, $\Delta_s/2\pi$				   MHz   -99   -101

spin-photon detuning (two-spins), $\Delta_{2s}/2\pi$		   MHz   -79   -79

spin dispersive shift, $\chi_s/2\pi$					   MHz   -7.8   -7.8

spin-spin coupling, $J/2\pi$						   MHz         9.5

spin linewidth, $\varGamma_s/2\pi$					   MHz   5.2   6.6

mean photon number, $\bar n_r$					            <0.5

5lPhoton number splitting, , $\phi = 10.55 \deg$

tunnel splitting, $2t_c/h$							   GHz   8.75

charge dispersive shift, $\chi_c/2\pi$					   MHz   22.18

spin-orbit shift, $\chi_\text{SO}/2\pi$					   MHz   -63.7

spin-photon coupling, $g_s/2\pi$					   MHz   33.4

spin-photon detuning, $\Delta_s/2\pi$				   MHz   -102

spin dispersive shift (predicted), $\chi_s/2\pi$			   MHz   -8.0

spin dispersive shift (measured), $\chi_s/2\pi$			   MHz   -6.6

zero-photon-peak linewidth, $\varGamma(0)/2\pi$			   MHz   4.3

one-photon-peak linewidth, $\varGamma(1)/2\pi$			   MHz   5.8

spin linewidth without photon losses, $\gamma_s/2\pi$	   MHz   3.4

interacting resonator linewidth, $\kappa'_r/2\pi$			   MHz         3.0

mean photon number, $\bar n_r$					            0.62

Here, we describe the procedure to extract the model parameters from the experiment.
Parameters are successively 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 spin-photon 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:extvrs1 (vacuum Rabi splitting, and two-tone spin measurement versus $B_r$ and $\phi$) for each dot independently.
Notably, the size of the spin-photon gap ($2g_s/2\pi$) in a vacuum Rabi splitting measurement is most affected by $\Delta B_\perp$. The micromagnet susceptibility $\chi_\text{um}$ and offset field $B_\text{um0}$ are mainly fixed by the slope of the spin transition frequency versus $B_r$ and the spin-photon resonance condition, respectively. Then, $B_\text{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:params.

§ 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\pi = 34.2 \MHz$. A pump tone at frequency $f_\text{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_\text{probe}$ set to the dispersively shifted resonator frequency for each DQD at zero charge detuning (e.g., as in the insets of ).
The results are plotted in fig:dipersivedetuning for two values of the magnetic field angle $\phi$, corresponding to off-resonant and resonant spin transition energies. The spin transition frequency is reduced at zero charge detuning compared with the localized dot states because of the spin-charge hybridization shift $\chi_\text{SO}/2\pi$. The signal shows a peak centered on $\epsilon = 0 \ueV$ when the resonator frequency gets pushed down by the charge dispersive shift $\chi_c/2\pi$. 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\chi_s/2\pi$, visible as a dip in the transmission that follows the spin transition energy.
When the two spins are not resonant ($\phi = 13.0 \deg$), 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 ($\phi = 10.8 \deg$ 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 from the symmetry of the hybridized spin states and is expected in the case of coherent spin-spin interaction <cit.>.
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\pi = 20.2 \MHz$, with $B_r = 52 \mT$ and $\varDelta_{2s} = -80 \MHz$.

§ 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:dispersivebr.

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 ; 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 $\phi$ is caused by jitter of the DQD2 tunnel coupling over the time it takes to reconstruct all angles. The $2t_c/h \approx 8.7 \GHz$ is retuned between plots, and yields $g_s/2\pi \approx 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 $t_c$.     f Splittings $2J/2\pi$ extracted from the model as a function of spin-photon detuning (through $B_r$). Comparison with the linewidth is impaired by the photon-number broadening of the transitions. As expected, $2J/2\pi$ decreases as $|\varDelta_{2s}/2\pi|$ increases.

To speed up data acquisition, a larger probe power is used than in the main text data of , which, as a consequence, also significantly broadens the width of the spin transitions due to the photon-number-dependent spin dispersive shift. The apparent jitter of the transition energy as a function of $\phi$ 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 $|\varDelta_{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:dispersivetc.

Extended spin-spin hybridization data versus spin-photon coupling strength.     a-d Dispersive spin-spin hybridization measured for decreasing spin-photon coupling values. The $g_s$ is adjusted through its dependence on $t_c$ and is approximately the same for both spins, while $B_r = 52 \mT$ is fixed. The $\varDelta_{2s}/2\pi \approx 78 \MHz$ is approximately unchanged between plots because changes in $\chi_\text{SO}$ are approximately compensated by changes in $\chi_c$. The model fully captures the observed transition frequencies within plots and between the plots, and the only adjustment is the experimentally measured $t_c$. As expected, $2J/2\pi$ decreases as $g_s/2\pi$ decreases.

In this work, we have measured spin-spin interactions ($J$) with values of spin-photon interactions $g_s/2\pi$ up to $40 \MHz$. A practical limit on how large $g_s$ can be is the broadening of $\kappa'_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  is chosen empirically based on linewidth, exchange coupling and readout signal.

In the experiment of the main text , the critical photon number $n_\text{crit} = (\varDelta_{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 number $\bar n_r$ to be $\bar n_r \ll n_\text{crit}$ <cit.>. We use the photon-number-dependent spin dispersive shift to establish an upper bound of $\bar n_r < 0.5$ for the data of the main text  ($\bar n_r \approx 0.4 \pm 0.1$). We can no longer distinguish between the coherent and thermal distributions since they converge and the data are too noisy. The $n_\text{crit}$ could be optimized for future qubit experiments. For instance, in fig:dispersivebrb, the $n_\text{crit} = 2$ is already much higher while $2J/2\pi$ is still $17 \MHz$. Time-domain control can also help reduce $\bar n_r$ by probing only during readout. The exact requirements will depend on the target two-qubit gate and desired fidelity (amongst other things).

§ EXTENDED PHOTON-NUMBER-DEPENDENT SPIN DISPERSIVE SHIFT DATA

In this section, we extend the analysis of the photon-number-dependent spin dispersive shift data by looking at data for $\phi = 11.6 \deg$.

Extended photon-number-dependent spin dispersive shift data.     a Photon-number-dependent spin dispersive shift of DQD1 (same as in the main text).     b Line cut of the data in panel (a) for a different angle $\phi = 11.6 \deg$, and fit to four Lorentzian dips (first three shown with dashed lines). A dispersive shift $2\chi_s/2\pi = -12.5 \MHz$ is extracted.      The coherent-state distribution has better agreement than the thermal distribution; see the main text for number distribution results.

In fig:photonnumberextb, a line cut from fig:photonnumberexta is extracted, and the dip areas and separations are fit to a sum of Lorentzian dips with $2\varGamma(n) = 2\gamma_s + (n + \bar n_r)\kappa'_r$. A value of $2\chi_s/2\pi = -12.5 \MHz$ is extracted from the fit ($g_s/2\pi = 33.4 \MHz$, $\varDelta_s/2\pi = -122 \MHz$), slightly larger than the average of the linewidths $(2\varGamma(0), 2\varGamma(1))/2\pi = (10.1, 13.1) \MHz$.
The coherent state has a higher photon number, $\bar n_r = 0.82$, than at $\phi = 10.55 \deg$ (see main text, $\bar n_r = 0.62$). This can be attributed to the resonator's higher dressed quality factor for $\phi > 11.2 \deg$ (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 $|\text{S21}|$ value at fixed probe power. The difference between the coherent-state and thermal-state distributions is more pronounced than in the main text.

In tab:params, a small discrepancy is observed between the predicted and measured values of $\chi_s$. This could be due to driven system dynamics, which are known to modify the splitting <cit.>.


# BayesDLMfMRI: Bayesian Matrix-Variate Dynamic Linear Models for Task-based
fRMI Modeling in R

Johnatan Cardona Jiménez
Facultad de Ingeniería
Institución Universitaria Pascual Bravo
Cl 73 No. 73A - 226, Medellín, ZIP 050034, Colombia
<EMAIL_ADDRESS>
Use footnote for providing further information about author (webpage,
alternative address)— _not_ for acknowledging funding agencies.

###### Abstract

This article introduces an R package to perform statistical analysis for task-
based fMRI data at both individual and group levels. The analysis to detect
brain activation at the individual level is based on modeling the fMRI signal
using Matrix-Variate Dynamic Linear Models (MDLM). Therefore, the analysis for
the group stage is based on posterior distributions of the state parameter
obtained from the modeling at the individual level. In this way, this package
offers several R functions with different algorithms to perform inference on
the state parameter to assess brain activation for both individual and group
stages. Those functions allow for parallel computation when the analysis is
performed for the entire brain as well as analysis at specific voxels when it
is required.

_Keywords_ First keyword $\cdot$ Second keyword $\cdot$ More

##  1 Introduction: fMRI packages with options for Bayesian Analysis

Statistical modeling of processed fMRI data is a challenging problem, which
has caught the attention of the statistical community in the past two decades.
Large amount of observations are usually obtained with only one subject in an
fMRI session which makes the implementation of more appropriate and
sophisticated statistical models that can account for the spatiotemporal
structures that are usually present in this type of data really challenging.
From a Bayesian perspective, there have been published and proposed different
types of models to model fMRI data (see for instance Zhang et al. (2016),
Eklund et al. (2017), Bezener et al. (2018), Yu et al. (2018) and a complete
review by Zhang et al. (2015)). Despite that, just a few software
implementations of Bayesian methods are available for fMRI data analysis. For
instance, the FSL software (Jenkinson et al., 2012) offers a Bayesian
procedure for the group stage analysis, which depends on a Frequentist or
Classical output from the individual stage. Another popular software among
practioners in the neurosicience community with implementations of Bayesian
methods is the SPM package (Penny et al., 2011). It has alternatives for
Bayesian analysis to both individual and group stages, but as in the case of
the FSL package, the group analysis depends on a Frequentist individaul stage
output. There is also available a MATLAB GUI called NPBayes-fMRI (Kook et al.,
2019), which is an implementation of the work proposed by Zhang et al. (2016),
a fully Bayesian modeling for individual and group stages. In R (R Core Team,
2018), tools for different types of fMRI data analysis are provided in
packages like fmri (Tabelow and Polzehl, 2011), oro.nifti (Whitcher et al.,
2011), neurobase (Muschelli, 2018), and neuRosim (Welvaert et al., 2011). For
fMRI data analysis under a Bayesian approach there was a package called
cudaBayesreg (da Silva et al., 2011), which was removed recently from the CRAN
repository. Thus, the package BHMSMAfMRI (Sanyal and Ferreira, 2019) is to our
knowledge the only Bayesian option for fMRI data analysis in R to this date.

In this work, we are presenting an R package called BayesDLMfMRI, which is an
implementation of the method proposed in Jiménez et al. (2019) and Jiménez
(2019). The BayesDLMfMRI permits to perform fMRI individual and group analysis
based on the Matrix-Variate Dynamic Linear Model (MVDLM) proposed by Quintana
(1985). Given this type of analysis usually involved large amounts of data, we
take advantage of the packages Rcpp (Eddelbuettel et al., 2011) and
RcppArmadillo (Eddelbuettel and Sanderson, 2014) in order to speed up the
computation time. BayesDLMfMRI also depends internally on the package pbapply
(Solymos and Zawadzki, 2019), which allows the user to visualize a progress
bar while the process is run in either sequence or parallel. In order to run
an analysis using our package the user must provide a 4D array (or a 4D$\times
N$ array for a group analysis with $N$ subjects) containing the sequence of
processed images and a design matrix whose columns are related to the so-
called expected blood-oxygen-Level dependent (BOLD) response and/or some other
covariates related to particular subjects’ characteristics. To process the raw
images, we recommend the use of packages such as FSL or SPM. And to build the
expected BOLD response the user has options like the fmri package, which
allows defining different types of models to represent the haemodynamic
response function (HRF). BayesDLMfMRI package is intended to be just another
well tested and assessed tool for the practitioners who are interested to
perform statistical analysis of processed fMRI images from a Bayesian
perspective. In the next section, we give a brief description of the model and
methods on which the BayesDLMfMRI package is based and present its R functions
to perform fMRI data analysis. In section three, we present some examples to
illustrate the use of the package and in the final section, we give some
concluding remarks.

##  2 Methods and software

Table 1:  Functions implemented in the BayesDLMfMRI package for fMRI (group)
individual analysis. LLT stands for Transformation Linear Test, which is the
same as the average cluster effect (ACE) distribution defined in ArchiveRef.

Function |  Description
---|---
(ffdGroupEvidenceFEST) ffdEvidenceFEST |  It returns 3D arrays to build (group) individual activation evidence maps (based on outputs from MVDLM) fitting an MVDLM at individual level and using the FEST algorithm to assess voxel activation. There are two options related to the posterior distribution of $\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}$: LLT and Joint. Each one runs independently and must be set by the user as an input parameter.
(ffdGroupEvidenceFFBS) ffdEvidenceFFBS |  It returns 3D arrays to build (group) individual activation evidence maps (based on outputs from MVDLM) fitting an MVDLM and using the FFBS algorithm to assess voxel activation. The options LLT and Joint related to the posterior distribution of $\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}$ are simultaneously executed in the same run.
(ffdGroupEvidenceFSTS) ffdEvidenceFSTS |  Same features as (ffdGroupEvidenceFFBS) ffdEvidenceFFBS, though using the FSTS algorithm.
(GroupSingleVoxelFEST) SingleVoxelFEST |  Produces some usefull outputs from a single voxel analysis related to the FEST algorithm
(GroupSingleVoxelFFBS) SingleVoxelFFBS |  Same features as (GroupSingleVoxelFEST) SingleVoxelFEST, though using the FFBS algorithm
(GroupSingleVoxelFSTS) SingleVoxelFSTS |  Same features as (GroupSingleVoxelFEST) SingleVoxelFEST, though using the FSTS algorithm

The type of MVDLM which this package relies on is the version originally
developed by Quintana (1985) and Quintana (1987). Here, we just give a brief
description of the model and methods. For a better understanding of the method
implemented in this package for the individual stage, see Jiménez et al.
(2019). Let $\mathbf{Y}_{\scaleto{[i,j,k]t\mathstrut}{6pt}}^{(z)}$ be a
$q\times 1$ random vector representing the cluster of observed BOLD responses
at position $(i,j,k)$ in the brain image, time $t$ and subject $z$, for
$i=1,\ldots,d_{1}$, $j=1,\ldots,d_{2}$, $k=1,\ldots,d_{3}$, $t=1,\ldots,T$ and
$z=1\ldots,N$. Thus, the cluster of BOLD signals is modeled as

$\begin{array}[]{lccl}\text{Observation:}&\mathbf{Y}_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}&=&\mathbf{F}^{{}^{\prime}}_{t}\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}+\bm{\nu}_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{{}^{\prime}(z)\mathstrut}{4pt}}\\\
\text{Evolution:}&\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}&=&\mathbf{G}_{t}\Theta_{\scaleto{[i,j,k]t-1\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}+\Omega_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}},\\\
\end{array}$ (1)

where, for each $t$ we have a $q\times 1$ vector
$\bm{\nu}_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{{}^{\prime}(z)\mathstrut}{4pt}}$
of observational errors, a $p\times q$ matrix
$\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}$ of state parameters, a $p\times
q$ matrix
$\Omega_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}$ of
evolution errors. The $1\times p$ and $p\times p$ matrices
$\mathbf{F}^{{}^{\prime}}_{t}$ and $\mathbf{G}_{t}$ respectivelly are common
to each of the $q$ univariate DLMs. The covariates related to the design being
used, either a block or an event-related design as well as other
characteristics of the subjects, can be included in the columns of
$\mathbf{F}^{{}^{\prime}}_{t}$. For individual analysis the model (1) is
fitted to every cluster of voxels related to position $[i,j,k]$ in the brain
image and the cluster size depends on the $r$ distance, which is a parameter
defined by the user. To identify wether or not there is significant evidence
of cluster activation or brain reaction at region $[i,j,k]$, three different
algorithms (FEST, FSTS and FFBS) proposed in Jiménez et al. (2019) are
implemented. Those algorithms are sampling schemes that allow to draw on-line
trajectory curves related to the state parameter
$\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}$,
and with those resulting samples a Monte Carlo evidence related to the event
of cluster activation is computed. The information obtained from the first
stage can be combined in different ways to produce several measures of
evidence for the group activation. For a better understanding of the method
implemented in this package for the group stage, see Jiménez (2019). In table
2, the package’s functions for (group) individual analysis and a brief
description of them are presented.

##  3 Illustrations

In order to show some practical illustrations about the use of the
BayesDLMfMRI package, we use data related to an fMRI experiment where a sound
stimilus is presented. That experiment is intended to offer a "voice
localizer" scan, which allows rapid and reliable localization of the voice-
sensitive "temporal voice areas" (TVA) of human auditory cortex (Pernet et
al., 2015). The data of this "voice localizer" scan is freely available on the
online platform OpenNEURO (Gorgolewski et al., 2017). In the original
experiment, the voice and non-voice sounds are separately analyzed, but here
we merge both sounds in one block as it were just one stimulus (see figure 1).
For the individual analysis, we select one from the 217 subjects whose data
are available on OpenNEURO, specifically, we take the data from sub-007. To
illustrate the use of the group analysis functions, we take 20 subjects
(sub-001:sub-021). The raw fMRI data is preprocessed using the standard
processing pipelines implemented on the FSL software for motion correction,
spatial smoothing and other necessary procedures to perform statistical
analysis for both individual and group stages.

Figure 1:  Merged stimuli function for voice and non-voice sounds with its
respective expected BOLD response function for the "voice localizer" example.

### Individual analysis

In order to run any of the functions related to individual analysis, the user
must provide two inputs: an array of four dimensions containing the MRI images
and a matrix whose columns contain the covariates to model the observed BOLD
response. Thus, we read the data.nii.qz file, which contains the MRI images,
using the function readNIFTI from the package oro.nifti and the covariates.cvs
file which contains the expected BOLD response (show in figure 1) and its
derivative respectively.

R> library("oro.nifti")R> fMRI.data <\- readNIfTI("./fMRIData/sub-007.nii.gz",
reorient=FALSE)R> fMRI.data <\-<EMAIL_ADDRESS>dim(fMRI.data)

[1] 91 109 91 310

R> Covariates <\- read.csv("./covariates.csv", header=FALSE, sep="")R>
dim(Covariates)

[1] 310 2

To perform a voxel-wise anaylsis 111The BayesDLMfMRI package depends on Rcpp,
RcppArmadillo, RcppDist and pbapply packages. to obtain a 3D array with
measurements of activation evidence for every voxel, the user has three
options (ffdEvidenceFEST, ffdEvidenceFFBS and ffdEvidenceFSTS), which at the
same time can yield three different types of evidence measurements for voxel
activation. Just to ilustrate their use, we run them using the "voice
localizer" example.

R> library(devtools)R> install_github("JohnatanLAB/BayesDLMfMRI")R>
library(BayesDLMfMRI)R> res <\- ffdEvidenceFEST(ffdc = fMRI.data, covariates =
Covariates,\+ m0 = 0, Cova = 100, delta = 0.95, S0 = 1, n0 = 1, Nsimu1 = 100,
Cutpos1 = 30,\+ r1 = 1, Test = "LTT", Ncores = 15)

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 05m 36s

The arguments m0, Cova, S0 and n0 are the hyper-parameters related to the
joint prior distribution of
$(\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}},\Sigma_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}})$.
For this example, we are setting a "vague" prior distribution acording to
Quintana and West (1987), where m0 = 0 define a null matrix with zero values
in all its entries and both Cova = 100 and S0 = 1 define diagonal matrices. r1
is the euclidean distance, which defines the size of the cluster of voxels
jointly modeled. Test is the parameter related to the test selected by the
user, for which there are two options: "LTT" and "Joint". "Ncores" is the
argument related to the number of cores when the process is executed in
parallel. Nsimu1 is the number of simulated on-line trajectories related to
the state parameter
$\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}$.
From our own experience dealing with different sets of fMRI data, we recommend
Nsimu1 = 100 as a good number of draws to obtain reliable results. Cutpos1 is
the time up from where the on-line trajectories are considered in order to
compute the activation evidence and delta is the value of the discount factor.
For a better understanding about the setting of these two last arguments, see
(Archive1).

R> str(res)

List of 2 $ : num [1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91,
1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ...

R> dim(res[[1]])

[1] 91 109 91

The output for the ffdEvidenceFEST function depends on the type of Test set by
the user. For Test = "LTT" the function returns a list of the type res[[p]][x,
y, z], where [[p]] represents the column position in the covariates matrix and
[x, y, z] represent the voxel position in the brain image. Thus, for the
"voice localizer" example res[[1]] and res[[2]] are the 3D arrays related to
the evidence for brain activation related to the BOLD response for the
auditory stimuli and its derivative respectively. When Test = Joint the output
returned is an array

R> library(neurobase)R> res.auxi <\- res[[1]]R> ffd <\-
readNIfTI("./standard.nii.gz")R> Z.visual.c <\- nifti(res.auxi, datatype=16)R>
ortho2(ffd, ifelse(Z.visual.c > 0.95, Z.visual.c, NA),\+ col.y =
heat.colors(50), ycolorbar = TRUE, ybreaks = seq(0.95, 1, by = 0.001))

The neurobase package is one amongst several options available in R to
visualize MRI images. In this example, we use its ortho2() function in order
to plot the evidence activation map. The standar.nii.gz contains the MNI brain
atlas, which is used in this work as a reference space for individual and
group analysis. For a better understanding of the use of brain atlas, see
Brett et al. (2002).

|
---|---
Figure 2: Activations maps obtained with the functions ffdEvidenceFEST (left)
and ffdEvidenceFFBS (right) using the "Joint" as sampler distribution.

R> res <\- ffdEvidenceFEST(ffdc = fMRI.data, covariates = Covariates,\+ m0 =
0, Cova = 100, delta = 0.95, S0 = 1, n0 = 1, Nsimu1 = 100, Cutpos1 = 30,\+ r1
= 2, Test = "JointTest", Ncores = 15)

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 28m 57s

R> str(res)

List of 4 $ : num [1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91,
1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91, 1:109, 1:91] 0 0 0 0 0 0 0
0 0 0 ...

For both ffdEvidenceFFBS and ffdEvidenceFSTS the input arguments and output
structures are the same. Below we run the "voice localizer" example using the
ffdEvidenceFFBS function. It returns a list of the form res[[T]][p,x,y,z],
where T defines the type of test (T = 1 for "Marginal", T = 2 for "JointTest",
and T = 3 for "LTT"), p represents the column position in the covariates
matrix and x, y, z represent the voxel position in the brain image.

#Change JointTest for JointR> res <\- ffdEvidenceFFBS(ffdc = fMRI.data,
covariates = covariables, m0 = 0,\+ Cova = 100, delta = 0.95, S0 = 1,n0 = 1,
Nsimu1 = 100, Cutpos1 = 30,\+ r1 = 1, Ncores = 15)

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 19m 08s

R> str(res)

List of 3 $ : num [1:2, 1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num
[1:2, 1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:2, 1:91, 1:109,
1:91] 0 0 0 0 0 0 0 0 0 0 ...$

R> library(neurobase)R> res.auxi <\- res[[3]][1,,,]R> ffd <\-
readNIfTI("standard.nii.gz")R> Z.visual.c <\- nifti(res.auxi, datatype=16)R>
ortho2(ffd, ifelse(Z.visual.c > 0.95, Z.visual.c, NA),\+ col.y =
heat.colors(50), ycolorbar = TRUE, ybreaks = seq(0.95, 1, by = 0.001))

The BayesDLMfMRI package also has functions that allow taking a closer look
for specific voxels defined by the user. For instance, let’s suppose we are
interested to see some of the output elements related to the FEST algorithm
for an active voxel when using the LTT test.

R> #Identifying active voxels for a probability threshold of 0.99R>
active.voxels = which(res[[1]] > 0.99, arr.ind = TRUE)R> head(active.voxels)

dim1 dim2 dim3[1,] 28 75 22[2,] 22 71 24[3,] 22 72 24[4,] 23 72 24[5,] 23 73
24[6,] 24 73 24

R> N1 <\- dim(covariables)[1]R> res.indi <\- SingleVoxelFEST(posi.ffd = c(14,
56, 40), covariates\+ = Covariates, ffdc = fMRI.data, m0 = 0, Cova = 100,
delta = 0.95, S0 = 1,\+ n0 = 1, Nsimu1 = 100, N1 = N1, Cutpos1 = 30, Min.vol =
0.10, r1 = 1,\+ Test = "LTT")R> str(res.indi)

List of 4$ Eviden : num [1, 1:2] 0.98 0$ Online_theta: num [1:280, 1:2, 1:100]
2.57 2.52 2.54 2.62 2.46 ...$ Y_simu : num [1:280, 1:100] 1.785 0.521 0.624
1.028 -0.165 ...$ FitnessV : num 0.762

The function IndividualVoxelFEST() requires just a few additional input
arguments: the position of the voxel in the brain image (posi.ffd), the last
period of time of the temporal series
($\mathbf{Y}_{\scaleto{[i,j,k]t=1\mathstrut}{5pt}}:\mathbf{Y}_{\scaleto{[i,j,k]t=N1\mathstrut}{5pt}}$)
that are used in the analysis (N1) and Min.vol, which helps to define a
threshold for the voxels that are considered in the analysis. For example,
Min.vol = 0.10 means that all the voxels with values below to max(fMRI.data) *
Min.vol are discarted from the analysis. The output is a list containing a
vector (Eviden) with the evidence measure of activation for each one of the p
covariates considered in the model, the simulated online trajectories of
$\Theta_{\scaleto{[i,j,k]t\mathstrut}{4pt}}^{\scaleto{(z)\mathstrut}{4pt}}$
(Online_theta), the simulated BOLD responses (Y_simu) and measure to examine
the goodnes of fit of the model
($100*|\mathbf{Y}_{\scaleto{[i,j,k]t\mathstrut}{5pt}}-\mathbf{\hat{Y}}_{\scaleto{[i,j,k]t\mathstrut}{5pt}}|/\mathbf{\hat{Y}}_{\scaleto{[i,j,k]t\mathstrut}{5pt}}$)
for that particular voxel (FitnessV).

R> res.indi2 <\- singleVoxelFEST(posi.ffd = c(14, 56, 40),\+ covariates =
covariates, ffdc = fMRI.data, m0 = 0, Cova = 100,\+ delta = 0.95, S0 = 1,n0 =
1, Nsimu1 = 100, N1=N1, Cutpos1 = 30, Min.vol=0.10,\+ r1 = 1, Test =
"JointTest")R> str(res.indi2)

List of 5$ EvidenMultivariate: num [1, 1:2] 0.81 0$ EvidenMarginal : num [1,
1:2] 0.98 0$ Online_theta : num [1:2, 1:280, 1:100] 2.38 -1.7 2.42 -1.77 2.35
...$ Y_simu : num [1:100, 1:7, 1:280] 0.597 0.883 1.716 0.11 0.344 ...$
FitnessV : num 0.805$

When Test = "Joint", the function IndividualVoxelFEST returns two measures of
voxel activation related to the joint and marginal tests along with the
already explained remain list elements.

R> res.indi3 <\- SingleVoxelFFBS(posi.ffd = c(14, 56, 40), covariates =
covariates,\+ ffdc = fMRI.data, m0 = 0, Cova = 100, delta = 0.95, S0 = 1, n0 =
1,\+ Nsimu1 = 100, N1 = N1, Cutpos1 = 30, Min.vol=0.10, r1 = 1)R>
str(res.indi3)

List of 5 $ Eviden_joint : num [1, 1:2] 0.94 0.03 $ Eviden_margin : num [1,
1:2] 1 0.08 $ eviden_lt : num [1, 1:2] 1 0.08 $ Online_theta : num [1:310,
1:2, 1:100] 0 0 0 0 0 0 0 0 0 0 ... $ Online_theta_mean: num [1:310, 1:2,
1:100] 0 0 0 0 0 0 0 0 0 0 ...$

For both functions IndividualVoxelFFBS and IndividualVoxelFSTS the input
arguments and output estructures are exactly the same. Three measures of
evidence related to the joint, marginal and linear transformation (LTT) tests
are generated in the same run. There are also two elements from the list
output related to the online simulated trajectories of the state parameter.
Online_theta are the draws for the marginal model and Online_theta_mean are
the mean draws obtained from the joint model.

R> frame()R> plot.window(xlim=c(30, 311), ylim=c(-5, 5))axis(1, at=seq(30,
310, by = 20), lwd = 2, xlab = "Time")R> box(lwd = 2)R> for(i in
1:dim(res.indi$Y_simu)[2]){lines(31:dim(covariates)[1],\+ res.indi$Y_simu[,
i])}R> lines(31:dim(covariates)[1], covariates[31:dim(covariates)[1], 1],\+
col = "red", lwd = 2)

R> frame()R> plot.window(xlim=c(30, 311), ylim=c(-0.5, 4.5), ylab =
expression(theta))R> axis(1, at=seq(30, 310, by = 20), lwd = 2, xlab =
"Time")R> axis(2, at=-1:4, lwd=2)R> box(lwd = 2)R> for(i in
1:dim(res.indi$Online_theta)[2]){lines(31:dim(covariates)[1],\+
res.indi$Online_theta[, i])}R> abline(h = 0, col = "red" , lwd = 2)

posi.ffd = c(14, 56, 40) | posi.ffd = c(28, 67, 15)
---|---
|
Figure 3: Simulated Bold response (black lines) and expected Bold response (red line) for an active (left) and non-active (right) voxel respectively. posi.ffd = c(14, 56, 40) | posi.ffd = c(28, 67, 15)
---|---
|
Figure 4:  Simulated online trajectories of the state parameter for an active
(left) and nonactive (right) volxel respectively.

### Group analysis

Now we ilustrate how to run an fMRI group analysis as it is described in
(RefARchive2). First, we read the fMRI images of 21 subjects from the "voice-
localizer" example:

R> names <\- list.files("~/fMRIData/")R> names

using the "Joint" as sampler distribution.

[1] "sub-001.nii.gz" "sub-002.nii.gz" "sub-003.nii.gz" [4] "sub-004.nii.gz"
"sub-005.nii.gz" "sub-006.nii.gz" [7] "sub-007.nii.gz" "sub-008.nii.gz"
"sub-009.nii.gz"[10] "sub-010.nii.gz" "sub-011.nii.gz" "sub-012.nii.gz"[13]
"sub-013.nii.gz" "sub-014.nii.gz" "sub-015.nii.gz"[16] "sub-016.nii.gz"
"sub-017.nii.gz" "sub-018.nii.gz"[19] "sub-019.nii.gz" "sub-020.nii.gz"
"sub-021.nii.gz"

R> DataGroups <\- function(x){\+ ffd.c2 <\-
readNIfTI(paste("~/fMRIParietal/",x, sep=""), reorient = FALSE)\+ ffd.c <\-
<EMAIL_ADDRESS>return(ffd.c)\+ }R> system.time(DatabaseGroup <\-
parLapply(names, DataGroups, cl = 7)) user system elapsed637.016 923.008
400.528R> str(DatabaseGroup)

List of 21 $ : num [1:91, 1:109, 1:91, 1:310] 0 0 0 0 0 0 0 0 0 0 ... $ : num
[1:91, 1:109, 1:91, 1:310] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91, 1:109, 1:91,
1:310] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91, 1:109, 1:91, 1:310] 0 0 0 0 0 0
0 0 0 0 ... $ : num [1:91, 1:109, 1:91, 1:310] 0 0 0 0 0 0 0 0 0 0 ... $ : num
[1:91, 1:109, 1:91, 1:310] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91, 1:109, 1:91,
1:310] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91, 1:109, 1:91, 1:310] 0 0 0 0 0 0
0 0 0 0 ...

In order to run any of the functions available in this package to perform fMRI
group analysis, the datasets or set of images from each subject must be stored
on a list object as it is shown above. To deal with this huge amount of
information the user has to have a big RAM memory capacity available on the
machine where this process is going to be run. It is also recommended to have
a multi-core processor available in order to speed up computation time. The
arguments or input parameters for any of the functions offered in this package
to run group analysis are almost the same as those required for individual
analysis. There is only an additional argument needed (mask), which adds a 3D
array that works as a brain of reference (MNI atlas) for the group analysis.

R> MASK <\- readNIfTI("~/mask.nii.gz")

R> res <\- ffdGroupEvidenceFEST(ffdGroup = DatabaseGroup,\+ covariates =
Covariates, m0 = 0, Cova = 100, delta = 0.95, S0 = 1,\+ n0 = 1, N1 = FALSE,
Nsimu1 = 100, Cutpos=30, r1 = 1, Test = "Joint",\+ mask = MASK, Ncores = 7)

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 53m 45s

R> str(res)

List of 4 $ : num [1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91,
1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:91, 1:109, 1:91] 0 0 0 0 0 0 0
0 0 0 ... $ : num [1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ...

ffdGroupEvidenceFEST returns an array of dimension $2\times p$ elements, where
$p$ is the number of covariates and $2$ is the number of options evaluated as
sampler distributions: "Joint" and "Marginal". The first $p$ elements are the
3D arrays related to each column of the covariates matrix respectively when
computing the activation evidence using the "Join" distribution. The ramaining
arrays are those related to the marginal distribution.

R> res2 <\- ffdGroupEvidenceFFBS(ffdGroup = DatabaseGroup, covariates =
Covariates,\+ m0=0, Cova=100, delta = 0.95, S0 = 1, n0 = 1, N1 = FALSE, Nsimu1
= 100,\+ Cutpos = 30, r1 = 1, mask = MASK, Ncores = 7)

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 44m 49s

R> str(res2)

List of 3 $ : num [1:2, 1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num
[1:2, 1:91, 1:109, 1:91] 0 0 0 0 0 0 0 0 0 0 ... $ : num [1:2, 1:91, 1:109,
1:91] 0 0 0 0 0 0 0 0 0 0$ ...

ffdGroupEvidenceFFBS returns an 3D array with the same structure and
characteristics as its individual counterpart.

R> library(neurobase)R> res.auxi <\- res2[[1]][1,,,]R> ffd <\-
readNIfTI("standard.nii.gz")R> Z.visual.c <\- nifti(res.auxi, datatype=16)R>
ortho2(ffd, ifelse(Z.visual.c > 0.95, Z.visual.c, NA),\+ col.y =
heat.colors(50), ycolorbar = TRUE, ybreaks = seq(0.95, 1, by = 0.001))

|
---|---
Figure 5: Activations maps obtained with the functions ffdGroupEvidenceFEST
(left) and ffdGroupEvidenceFFBS (right) using the "Joint" as sampler
distribution.

The functioning of the functions for single-voxel analysis at the group stage
is the same as their counterparts at the individual stage.

R> resSingle <\- GroupSingleVoxelFEST(posi.ffd = c(14, 56, 40),
DatabaseGroup,\+ covariates = Covariates, m0 = 0, Cova = 100, delta = 0.95, S0
= 1, n0 = 1,\+ N1 = FALSE, Nsimu1 = 100, r1 = 1, Test = "Joint", Cutpos = 30)

R> frame()R> plot.window(xlim=c(30, 311), ylim = c(-5, 5))axis(1, at=seq(30,
310, by = 20), lwd = 2, xlab = "Time")R> box(lwd = 2)R> for(j in 2:7){for(i in
1:dim(resSingle[[4]])[1]){lines(31:dim(covariates)[1],\+ resSingle[[4]][i,
j,], col = "green")}}R> for(i in
1:dim(resSingle[[4]])[1]){lines(31:dim(covariates)[1],\+ resSingle[[4]][i,
1,], col = "black")}lines(31:dim(covariates)[1],
covariates[31:dim(covariates)[1], 1],\+ col = "red", lwd = 2)

posi.ffd = c(14, 56, 40) | posi.ffd = c(28, 67, 15)
---|---
|
Figure 6: The black curves are simulated BOLD responses for voxel at position posi.ffd and the green curves are the BOLD responses related to its neighbors. The red curve is the expected BOLD response for the "voice localizer" experiment. posi.ffd = c(14, 56, 40) | posi.ffd = c(28, 67, 15)
---|---
|
Figure 7:  Simulated online trajectories of the state parameter for an active
(left) and nonactive (right) volxel respectively.

R> frame()R> plot.window(xlim=c(30, 311), ylim=c(-0.5, 4.5), ylab =
expression(theta))R> axis(1, at=seq(30, 310, by = 20), lwd = 2, xlab =
"Time")R> axis(2, at=-1:4, lwd=2)R> box(lwd = 2)R> for(i in
1:dim(resSingle[[3]])[3]){lines(31:dim(covariates)[1],resSingle[[3]][1, ,
i])}R> abline(h = 0, col = "red" , lwd = 2)

##  4 Conclusions and future work

In this work, we present the BayesDLMfMRI package, which allows performing
statistical analysis for fMRI data at individual and group stages. It offers
different options to assess brain activation for single voxels as well as the
entire brain volume and/or more specific brain regions with the help of a
mask. The low-level functions are written in C++ and options for parallel
computation are available in some of the functions. Currently, some extensions
for this package related to comparisons between groups and comparisons between
tasks as well as other relevant features are being developed.

## Computational details

The results in this paper were obtained using R 3.4.4 with the Rcpp 1.0.2,
RcppArmadillo 0.9.200.5.0 and pbapply 1.4.2 packages on a computer with Linux-
Ubuntu, 32 CPUs and 188GB of RAM. R itself and all packages used are available
from the Comprehensive R Archive Network (CRAN) at
https://CRAN.R-project.org/.

## References

  * Zhang et al. [2016] Linlin Zhang, Michele Guindani, Francesco Versace, Jeffrey M Engelmann, Marina Vannucci, et al.  A spatiotemporal nonparametric bayesian model of multi-subject fmri data.  _The Annals of Applied Statistics_ , 10(2):638–666, 2016.
  * Eklund et al. [2017] Anders Eklund, Martin A Lindquist, and Mattias Villani.  A bayesian heteroscedastic glm with application to fmri data with motion spikes.  _NeuroImage_ , 155:354–369, 2017.
  * Bezener et al. [2018] Martin Bezener, John Hughes, Galin Jones, et al.  Bayesian spatiotemporal modeling using hierarchical spatial priors, with applications to functional magnetic resonance imaging.  _Bayesian Analysis_ , 2018.
  * Yu et al. [2018] Cheng-Han Yu, Raquel Prado, Hernando Ombao, and Daniel Rowe.  A bayesian variable selection approach yields improved detection of brain activation from complex-valued fmri.  _Journal of the American Statistical Association_ , (just-accepted):1–61, 2018.
  * Zhang et al. [2015] Linlin Zhang, Michele Guindani, and Marina Vannucci.  Bayesian models for fmri data analysis.  _Wiley interdisciplinary reviews. Computational statistics_ , 7(1):21–41, 2015.
  * Jenkinson et al. [2012] Mark Jenkinson, Christian F Beckmann, Timothy EJ Behrens, Mark W Woolrich, and Stephen M Smith.  Fsl.  _Neuroimage_ , 62(2):782–790, 2012.
  * Penny et al. [2011] William D Penny, Karl J Friston, John T Ashburner, Stefan J Kiebel, and Thomas E Nichols.  _Statistical parametric mapping: the analysis of functional brain images_.  Elsevier, 2011.
  * Kook et al. [2019] Jeong Hwan Kook, Michele Guindani, Linlin Zhang, and Marina Vannucci.  Npbayes-fmri: Non-parametric bayesian general linear models for single-and multi-subject fmri data.  _Statistics in Biosciences_ , 11(1):3–21, 2019\.
  * R Core Team [2018] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2018.  URL https://www.R-project.org/.
  * Tabelow and Polzehl [2011] Karsten Tabelow and Jörg Polzehl.  Statistical parametric maps for functional mri experiments in R: The package fmri.  _Journal of Statistical Software_ , 44(11):1–21, 2011.  URL http://www.jstatsoft.org/v44/i11/.
  * Whitcher et al. [2011] Brandon Whitcher, Volker J. Schmid, and Andrew Thornton.  Working with the DICOM and NIfTI data standards in R.  _Journal of Statistical Software_ , 44(6):1–28, 2011.  URL http://www.jstatsoft.org/v44/i06/.
  * Muschelli [2018] John Muschelli.  _neurobase: ’Neuroconductor’ Base Package with Helper Functions for ’nifti’ Objects_ , 2018.  URL https://CRAN.R-project.org/package=neurobase.  R package version 1.27.6.
  * Welvaert et al. [2011] Marijke Welvaert, Joke Durnez, Beatrijs Moerkerke, Geert Verdoolaege, and Yves Rosseel.  neurosim: An r package for generating fmri data.  _Journal of Statistical Software_ , 44(10):1–18, 2011.
  * da Silva et al. [2011] AR Ferreira da Silva et al.  cudabayesreg: parallel implementation of a bayesian multilevel model for fmri data analysis.  _Journal of Statistical Software_ , 44(4):1–24, 2011.
  * Sanyal and Ferreira [2019] Nilotpal Sanyal and Marco A.R. Ferreira.  _BHMSMAfMRI: Bayesian Hierarchical Multi-Subject Multiscale Analysis of Functional MRI Data_ , 2019.  URL https://CRAN.R-project.org/package=BHMSMAfMRI.  R package version 1.3.
  * Jiménez et al. [2019] Johnatan Cardona Jiménez, Carlos A. de B. Pereira, and Victor Fossaluza.  Assessing dynamic effects on a bayesian matrix-variate dynamic linear model: an application to fmri data analysis.  _arXiv:1910.12058_ , 2019.
  * Jiménez [2019] Johnatan Cardona Jiménez.  fmri group analysis based on outputs from matrix-variate dynamic linear models.  _arXiv:1911.00708_ , 2019.
  * Quintana [1985] Jose Mario Quintana.  A dynamic linear matrix–variate regression model, 1985.
  * Eddelbuettel et al. [2011] Dirk Eddelbuettel, Romain François, J Allaire, Kevin Ushey, Qiang Kou, N Russel, John Chambers, and D Bates.  Rcpp: Seamless r and c++ integration.  _Journal of Statistical Software_ , 40(8):1–18, 2011.
  * Eddelbuettel and Sanderson [2014] Dirk Eddelbuettel and Conrad Sanderson.  Rcpparmadillo: Accelerating r with high-performance c++ linear algebra.  _Computational Statistics and Data Analysis_ , 71:1054–1063, March 2014.  URL http://dx.doi.org/10.1016/j.csda.2013.02.005.
  * Solymos and Zawadzki [2019] Peter Solymos and Zygmunt Zawadzki.  _pbapply: Adding Progress Bar to ’*apply’ Functions_ , 2019.  URL https://CRAN.R-project.org/package=pbapply.  R package version 1.4-1.
  * Quintana [1987] Jose Mario Quintana.  _Multivariate Bayesian forecasting models_.  PhD thesis, University of Warwick, 1987.
  * Pernet et al. [2015] Cyril R Pernet, Phil McAleer, Marianne Latinus, Krzysztof J Gorgolewski, Ian Charest, Patricia EG Bestelmeyer, Rebecca H Watson, David Fleming, Frances Crabbe, Mitchell Valdes-Sosa, et al.  The human voice areas: Spatial organization and inter-individual variability in temporal and extra-temporal cortices.  _Neuroimage_ , 119:164–174, 2015.
  * Gorgolewski et al. [2017] Krzysztof Gorgolewski, Oscar Esteban, Gunnar Schaefer, Brian Wandell, and Russell Poldrack.  Openneuro - a free online platform for sharing and analysis of neuroimaging data.  _Organization for Human Brain Mapping. Vancouver, Canada_ , 1677, 2017\.
  * Quintana and West [1987] Jose Mario Quintana and Mike West.  An analysis of international exchange rates using multivariate dlm’s.  _The Statistician_ , pages 275–281, 1987.
  * Brett et al. [2002] Matthew Brett, Ingrid S Johnsrude, and Adrian M Owen.  The problem of functional localization in the human brain.  _Nature reviews neuroscience_ , 3(3):243, 2002\.