# Fileset

[2406.04199v3_High-Fidelity Electron Spin Gates for Scaling Diamond Quantum Registers.pdf](https://mdr.nims.go.jp/filesets/9763de3e-62d2-43dc-8e32-274299b41b6e/download)

## Creator

T. Joas, F. Ferlemann, R. Sailer, P. J. Vetter, J. Zhang, R. S. Said, [T. Teraji](https://orcid.org/0000-0002-7731-0547), S. Onoda, T. Calarco, G. Genov, M. M. Müller, F. Jelezko

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[High-Fidelity Electron Spin Gates for Scaling Diamond Quantum Registers](https://mdr.nims.go.jp/datasets/7d4f0a2b-5344-497e-acd4-eb50be307670)

## Fulltext

1 High-Fidelity Electron Spin Gates for Scaling Diamond Quantum Registers T. Joas1, F. Ferlemann2,5, R. Sailer1, P.J. Vetter1, J. Zhang1, R. S. Said1, T. Teraji3, S. Onoda4, T. Calarco2,5,6, G. Genov1, M. M. Müller2, F. Jelezko1,7 1 Institute for Quantum Optics, Ulm University, 89081 Ulm, Germany 2 Peter Grünberg Institute - Quantum Control (PGI-8), Forschungszentrum Jülich GmbH, 52428 Jülich, Germany 3 Research Center for Electronic and Optical Materials, National Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan 4 Takasaki Advanced Radiation Research Institute, National Institutes for Quantum Science and Technology (QST), Takasaki, Gunma 370-1292, Japan 5 Institute for Theoretical Physics, University of Cologne, 50937 Köln, Germany 6 Dipartimento di Fisica e Astronomia, Università di Bologna, 40127 Bologna, Italy 7 Center for Integrated Quantum Science and Technology (IQST), Ulm University, 89081 Ulm, Germany   Diamond is a promising platform for quantum information processing as it can host highly coherent qubits that could allow for the construction of large quantum registers. A prerequisite for such devices is a coherent interaction between nitrogen vacancy (NV) electron spins enabling scalable entanglement. Entanglement between dipolar-coupled NV spin pairs has been demonstrated, but with a limited fidelity and its error sources have not been characterized. Here, we design and implement a robust two-qubit gate between NV electron spins in diamond and quantify the influence of multiple error sources on the gate performance. Experimentally, we demonstrate a record gate fidelity of 𝐹2𝑞 = (96.0 ±2.5) % under ambient conditions. Our identification of the dominant errors paves the way towards NV-NV gates beyond the error correction threshold.   I. INTRODUCTION Quantum processors have evolved from a scientifically intriguing concept into powerful devices on the verge of solving certain computational problems efficiently1–3. Their success has been enabled by the increased number of controllable qubits on recent devices, which can reach three-digit numbers3,4. Nevertheless, many compelling computational problems require operating quantum error correction protocols on thousands of individually controllable qubits, connected through high-fidelity multi-qubit gates. Each quantum platform faces technical challenges in scaling towards this goal5–8. Actual devices supporting error correction have only modest size9–12 or have not been tasked with computational problems yet4. Nuclear spin qubits in diamond, which are controlled by a single nitrogen-vacancy (NV) electron spin, offer gate fidelities among the leading platforms13,14, reaching 99.9415 % and 99.9315 % for single-qubit and two-qubit gates, respectively, and support register sizes of up to ten nuclear spins16,17. Their capability to operate under ambient conditions can lower requirements on experimental infrastructure and classical control electronics, potentially allowing preservation of high-fidelity gates as the number of qubits increases. However, scaling up the register size in diamond with a single NV is restricted by the limited range of the nuclear dipolar interaction and spectral crowding, making diamond a challenging candidate for building quantum processors, capable of error correction. 2 To overcome these scaling limitations, photon-mediated entangling gates18 have been implemented as a proof of concept to optically link separate nuclear registers19,20, albeit with low entanglement rates. Complementarily, it is appealing to directly connect multiple nuclear spin registers within the same diamond by coupling their central NV spins via the longer-range electron dipolar interaction. A central capability for such a connection, electron spin entanglement, has already been achieved but with a limited fidelity21–23. Improved performance is possible for state-to-state transfers in NV-NV registers with optimal control microwave pulses24. Yet, optimization of a full two-qubit gate for arbitrary input states remains a significant challenge at room temperature as the available error metrics imply increased measurement time and higher computational overhead15,25,26. Furthermore, experimental implementation is often impeded by the individual, non-linear microwave response of each experimental setup.  In this work, we apply an alternative approach to optimal control for fidelity improvement, where we characterize the physical error sources. Starting with a relatively simple gate, we analyze how errors affect its dynamics via a combination of experiments and simulations. Modeling the physics of the full system of four spins (two electron + two nuclear spins) enables us to quantify gate error sources that determine the entangling gate fidelity and to mitigate their effect in a targeted manner. This allows us to apply simple, sine-envelope27,28 microwave pulses while optimizing the entangling gate parameters (see Methods), which directly contribute to gate infidelity.  To demonstrate the effectiveness of our approach, we use randomized benchmarking to measure two-qubit gate fidelities in an electron spin quantum register in ambient conditions. We achieve a two-qubit gate fidelity of 𝐹2𝑞 = (96.0 ± 2.5) %, outperforming the reported Bell state fidelities in the literature24 – even without the previously needed optimal control pulses.  As a result of the optimization, our gate operates close to the 𝑇2 decoherence limit, and we characterize the remaining coherent errors. For diamond quantum computing hardware beyond the noisy intermediate-scale quantum era (NISQ29), the goal will be to surpass the error correction threshold. Our model projects that this objective is attainable at room temperature, and it identifies the physical mechanisms that promise the largest enhancements in fidelity. The prolonged electron coherence time (up to 𝑇2~1s30) at cryogenic temperature could further improve NV-NV gate fidelity.  Figure 1. Diamond quantum register. (a) Left: The register is fabricated by implanting adenine molecules into an epitaxially-grown diamond layer. Right: The coupled NV system (consisting of two NV centers plus their inherent 14N nuclear spins) is accessed optically with a confocal microscope. Spin states under a bias magnetic field 𝐵⃗  are manipulated through a microwave antenna. (b) Optically detected magnetic resonance spectrum of the register in panel (a). We employ the marked transitions 1 and 2 as “target qubit” and “control qubit”. In principle, these labels are interchangeable. Inset: Sketched energy diagram of a single NV without and with (blue background) magnetic field applied. Here, the control qubit is better aligned with the magnetic field and thus shows a higher splitting in the microwave domain. 3 II. RESULTS AND DISCUSSION a) Diamond quantum register The quantum register for our gate is fabricated by implanting C5N4Hn molecules from an adenine ion source31 into a 12C enriched, epitaxially grown diamond layer (see Fig. 1a and details in Methods section).  Two of the four nitrogen atoms in the molecule can end up in close proximity (≈9 ±  4 nm) in the diamond lattice and form two negatively charged nitrogen vacancy (NV) centers. The estimated probability of finding an NV pair in a confocal spot is ~15 %. NV- initialization and readout are achieved optically: We initialize the spin state with a green laser pulse. If required, the charge state is probed by a weak orange laser pulse exciting only the negative charge state of the NV center32. Note that the NV centers are too close to be optically resolvable based on their spatial separation. Hence, room temperature, state-of-the-art readout allows only to monitor the total fluorescence from both qubits upon green laser illumination. The readout signal is proportional to the summed qubit populations, but lacks information about correlations between them (see Methods)21. Each NV in the pair has a different orientation in the crystal lattice, resulting in spin sublevels of different energies when we apply a magnetic bias field, e.g., |𝐵|~100 G (Suppl. Note 1), of appropriate orientation. A continuous-wave optically detected magnetic resonance (ODMR) experiment thus shows four lines; two below and two above the zero-field splitting energy (ZFS, Fig. 1b). We define our qubit subspace by the choice of microwave transitions shown in the inset of Fig. 1b. The difference in transition frequencies allows for single qubit rotations by applying pulses with frequencies that match the resonance of the respective qubit. We can detect the presence of a magnetic dipolar interaction of 𝑣𝑑𝑖𝑝~0.1 MHz between NVs with a Double Electron Electron Resonance33 (DEER) experiment. As the coupling is smaller than the reciprocal dephasing time 1/𝑇2∗~0.5 MHz, the entangling gate needs to be embedded into a decoherence-protecting dynamical decoupling pulse sequence. Extending upon previous work21 to increase the coherence during the operation, we employ a two-qubit gate that applies robust XY8 microwave pulse sequences34 on both NVs (see Fig. 2a). Tuning the spin flip time, 𝜏2 between the center of the π pulses on the second NV and the times when the first NV is refocused, we partially refocus the dipolar interaction, adjusting the effective interaction time and acquired phase. Equipped with single-qubit rotations and an entangling operation, we have a gate set that in principle allows for arbitrary two-qubit quantum computations35.  b) Entangling gate  Our entangling gate, labelled the √𝑍𝑍 gate is:  𝑈√𝑍𝑍 = 𝑑𝑖𝑎𝑔(1, 𝑖, 𝑖, 1) (1) in the basis of the computational basis states {|00⟩, |01⟩, |10⟩, |11⟩} and up to a global phase (see Suppl. Note 4). We obtain additional intuitive understanding from the spin evolution on the Bloch sphere. To this end, we bring the first NV (“target qubit”) into a superposition state and apply our √𝑍𝑍 gate sequence. Choosing such an input state renders the experiment a dynamically decoupled DEER36 and varying the spin flip time 𝜏2 allows us to measure the dipolar coupling (see Fig. 2).  4  Figure 2. Entangling gate. (a) Microwave and laser sequence for probing the √𝑍𝑍 entangling gate. Blue: Envelope-shaped microwave pulses resonant with NV 1 (target qubit) and NV 2 (control qubit) of labeled phase (X or Y) and pulse area (𝜋/2 or 𝜋). For each NV, 𝑁𝜋 = 8 (phase cycled with XY8-1) 𝜋 pulses inside the grey brackets form the √𝑍𝑍 gate. For the sake of a clear figure, only two of the eight pulses per NV are drawn. The 𝜋/2 pulses around the gate are used to initialize the target qubit to the Bloch sphere equator and map the evolution under the gate onto the 𝜎𝑧 readout axis. Green: 3 μs pulses of the 552 nm laser used for spin initialization and readout (blue stroke). Orange: 3.5 ms pulses of the the 594 nm laser for optional charge state initialization. The unitary matrix on the right is realized for a calibrated 𝜏1 = 800 ns, 𝜏2 = 𝑡√𝑍𝑍/𝑁𝜋, 𝑛 = 1. (b)&(c) Measured evolution of the 𝜎𝑥, 𝜎𝑦 component of the target qubit for varying 𝜏2 in the √𝑍𝑍 gate sequence (𝑡𝑒𝑣𝑜𝑙 = 𝑁𝜋𝑛𝜏2, here 𝜏1 = 3000 ns, 𝑛 = 1, 𝑁𝜋 = 8; fitted 5 with 𝑦 = 𝐴𝑠 sin(2𝜋𝑣𝑑𝑖𝑝 +𝜙𝑠) + 𝑦0,𝑠, where 𝐴𝑠 , 𝜙𝑠 , 𝑦0,𝑠 are free parameters). The √𝑍𝑍 gate is realized after 𝑡√𝑍𝑍. The target qubit is initialized and mapped onto 𝜎𝑧 as given by the respective microwave sequence. Charge initialization with threshold parameter 𝑛𝑡ℎ𝑟𝑒𝑠ℎ = 8, 9 is applied for (c). (d) 𝑇2 coherence measurements for both NVs by applying a XY8-n sequence with fixed pulse spacing 𝜏1=800 ns sequence and varying order 𝑛𝑥𝑦 (total sequence duration 𝑡𝑓𝑟𝑒𝑒 = 8𝑛𝑥𝑦𝜏1) to only one of the NVs, respectively. Fixing 𝜏1 while increasing 𝑛𝑥𝑦 naturally compares to repeated application of our gate sequence that uses fixed 𝜏1, 𝜏2. (e) Idealized spin dynamic under the gate with initial 𝜋/2 pulse as in (b)&(c) on the Bloch sphere equator of the target qubit for the cases control qubit state = |0⟩/|1⟩ (blue/ red). A control qubit conditioned rotation from 𝜏2 = 0 with no acquired phase evolves for 𝜏2 > 0. In the 𝜎𝑥 component, we observe (Fig. 2b) a sinusoidal evolution arises as more phase is acquired through the dipolar interaction by a longer 𝜏2 (sketch in Fig. 2e). The controlled nature of the gate that is required for a computationally complete gate set is revealed by the reversed direction of the 𝜎𝑥 DEER oscillation when the control qubit is initialized into |1⟩. We extract the dipolar coupling 𝑣𝑑𝑖𝑝 =(119.8 ± 1.0) kHz between the two NVs from the fitted sine oscillation frequency. Thus, we find a calibrated evolution time 𝑡√𝑍𝑍 =𝑁𝜋𝜏2,√𝑍𝑍 = (2.17 ± 0.02) μs and describe details and optimization of the 𝜋 pulse spacing 𝜏1 = 800 ns in Suppl. Note 4. When calibrating, we minimize unwanted interactions with the inherent 14N nuclear spins37 due to the misaligned magnetic field to the NV axes.  An important aspect of the gate performance is inferred from the sine oscillation amplitude and offset parameters (ideally: 𝐴𝑠 = 1, 𝑦0,𝑠 = 0). The apparent non-zero offset in Fig. 2b can be explained by the imperfect initialization of the NV charge state. If any of the NVs is neutrally charged (NV0), no ZZ phase is collected during the √𝑍𝑍 gate; the output state is independent of 𝜏2 and thus far off the target state. Applying a green (552 nm) laser pulse yields a steady state charge distribution of ~0.7 NV- and ~0.3 NV0 for a single NV38 and thus a probability of 𝑝(𝑁𝑉−, 𝑁𝑉−) = 0.49 for both NVs to be in the desired negative charge state. Overall, the observed oscillation is an average over the different charge cases. This lowers the contrast of the 𝜎𝑥, 𝜎𝑦 oscillation in Fig. 2b. Second, the 𝜎𝑦 component becomes asymmetric (offset parameter of fitted sine 𝑦0,𝑠 < 0), as this component stays 〈𝜎𝑦〉𝑁𝑉1 = −1 for all 𝜏2 given a wrong charge state and the degree of asymmetry (|𝑦0,𝑠|/𝐴𝑠 ) can be a direct measure of the expected charge state of the NV pair (Suppl. Note 5).  In Fig. 2c, we mitigate the charge initialization error by using a weak orange laser  that employs thresholding (with a parameter nthresh) over data sets containing all detected fluorescence photons32. The DEER oscillations in Fig. 2c, show improvements in contrast and asymmetry of the 𝜎𝑦 component. Yet, there remains a significant deviation from full, symmetric contrast. This is in accordance with the limited charge state fidelity 𝐹𝑁𝑉−,𝑁𝑉− = (83 ±6) % measured in a separate measurement in Suppl. Fig. 4. We also show that the charge state fidelity can be improved by increasing the threshold parameter 𝑛𝑡ℎ𝑟𝑒𝑠ℎ of the post-selection. However, stricter thresholding comes at the cost of a worsened signal-to-noise ratio. Nevertheless, higher photon collection efficiency, provided, e.g., by micro structured lenses39, multi-laser pulse excitation40 or doping-induced NV charging41–43, could enhance charge state fidelity significantly.  Given the minimal gate time 1/(2𝑣𝑑𝑖𝑝) = 4.2 μs and the coherence times 𝑇2,𝑋𝑌8 = (454 ±58) μs, (476 ±30) μs for each NV measured in Fig. 2d, we calculate a coherence limit for the entangling gate fidelity of 99.1 %, and 98.6 % for our exact microwave gate with the actual 𝜏1. However, our analysis in Section E reveals that, 6 in the long coherence regime enabled by dynamical decoupling, gate errors related to the 14N nuclear spins become significant. We note that the slight asymmetry of the applied sequence on NV 2 might not be optimal in terms of robustness44 but such sequences have been demonstrated to be favorable for decoupling in the presence of cross-talk errors45. In addition, asymmetry allows us flexibility in terms of the choice of 𝜏1 and the total gate duration.  c) Repetitive benchmarking Imperfections that are not related to the state evolution during a gate are commonly referred to as state preparation and measurement errors (SPAM)46.  We use repetitive benchmarking, i.e. repeated application of the gate to a certain input state47, to separate the entangling gate fidelity from SPAM errors. We demonstrate that this tool quantifies gate errors independently of the charge state. The latter error can be significantly reduced by charge state mitigation strategies, as described in the previous section.  Separation of the gate error is achieved by repeated application of the gate under test. After n repetitions, single qubit rotations are applied to map back to the ground state and the 〈𝜎𝑧〉 expectation value is measured. In the resulting decay curves, errors occurring during the gate operation are collected in the lifetime parameter of a fitted model which we convert to a pseudo error per gate (pEPG) metric (see Methods).  7  Figure 3. Repetitive benchmarking. (a)&(c) Measured surviving population after n applications of the √𝑍𝑍 gate. For each data point, a reverse circuit is added that returns to |00⟩ in the absence of gate and preparation errors. Data for exemplary input states, an exponential fit (grey line) to extract the pEPG for input state (|0⟩ − 𝑖|1⟩)⨂ |0⟩  (simplifying notation, this abbreviates |𝛹⟩𝑁𝑉1⨂|𝛹⟩𝑁𝑉2 =|0⟩−𝑖|1⟩√2⨂|0⟩ with dropped normalization throughout this work), and simulation results including the charge state behavior are shown (colored solid lines, charge state probabilities annotated). (b)&(d) Extracted pEPG from the data in (a)&(c) for all probed input states. The charge-state-induced modulation in the decay data described in the main text leads to the plotted 1𝜎 fit errors that depend on the input state. Blue shading represents the standard deviation of all data points. Charge initialization is applied for (c)&(d). Mean is calculated as average over the mean pEPGs of each input state group to ensure comparability between (b)&(d).   Figure 3a shows decays for repeated √𝑍𝑍 gates for three input states. We find a mean 𝑝𝐸𝑃𝐺  =0.035 ±  0.017 averaged over all tested input states (see Fig. 3b), corresponding to a pseudo gate fidelity of 𝑝𝐹2𝑞 = 1 − 𝑝𝐸𝑃𝐺 = (96.5 ±1.7) %. This is only a factor of 2.5 away from the T2 limit 𝑝𝐸𝑃𝐺𝑇2 = 0.014 that we calculate from an independent coherence measurement (Fig. 2b). Thus, coherent gate errors are qualitatively small and we will later quantify the different contributions.  8 The input states can be grouped into three classes of pEPGs: The computational basis states show the lowest gate errors due to their magnetic noise insusceptibility. Gate errors increase when one qubit is in a superposition input state because of decoherence during the gate operation. Finally, input states with both qubits in superposition yield entangled states for gate repetitions n with mod(n,2)=1. Thus, not only both qubits will be affected by decoherence, but also some entangled states will pick up this noise more efficiently48. A prominent feature in repetitive benchmarking is an additional modulation with a period of n=4 for input states that are not the computational basis states. This behavior is a consequence of the charge state SPAM error. When one of the NVs is initialized as NV0, the coupling is 𝜈𝑑𝑖𝑝 = 0. This SPAM error has no effect when the initial state is a computational basis state. However, imperfect initialization reduces the output state fidelity when any of the NVs is in a superposition state with double superposition states showing the greatest error (see Fig. 3a). We simulate the repetitive benchmarking experiment in Fig. 3a&c considering all charge state configurations (details in Method section). The experimentally observed dependence of the modulation depth on the chosen input state is well reproduced.  When adding charge state initialization to the repetitive benchmarking in Fig. 3c, the modulation is reduced and we extract a slightly elevated 𝑝𝐸𝑃𝐺 = 0.043 ± 0.027. Comparing this mean 𝑝𝐸𝑃𝐺 for the charge-initialized case in Fig. 3d with the steady state initialized data in Fig. 3b, we find agreement within the standard deviation calculated from the data of all input states.  d) Randomized benchmarking Randomized benchmarking46 has emerged as a standard tool to evaluate gate performance, yielding a fidelity metric of an average computation with a priori limited insight into the nature of the gate error. In contrast to repetitive benchmarking, this fidelity can be directly compared to values obtained on different quantum hardware platforms49–51. We perform two-qubit randomized benchmarking in Fig. 4a&b. (see details in Methods section) 9  Figure 4. Randomized benchmarking and gate error analysis. (a) Sketched gate sequence for a single randomized benchmarking data point. Random Clifford gates are generated and reversed to the ground state. (b) Surviving population after randomized benchmarking in our optimized experimental setting (magnetic field setting 2, Ω𝑟𝑎𝑏𝑖/(2𝜋) = 23.7 MHz, 20 random experiments per Clifford length). From a single-exponential fit, the EPC of the gate set discussed in the main text is extracted. (c) Randomized benchmarking at a fixed Clifford length 𝑛𝑐𝑙𝑖𝑓𝑓 = 2 for varying Rabi frequencies (Ω𝑟𝑎𝑏𝑖 =𝜋𝑡𝜋 with approx. same duration 𝑡𝜋 of sine envelope-shaped 𝜋 pulses for both NVs). Each data point is a mean of 20 random experiments. The solid blue (red) line is a simulation at magnetic field setting 1, ~60 G (magnetic field setting 2, ~100 G), which contains a free parameter to account for the experimental SPAM errors (see Methods section). From the dashed lines, which are simulations at magnetic field 2 where the labeled error sources, i.e. crosstalk + leakage, unpolarized 14N spin, magnetic field orthogonal to the orientations of the NVs ez, are turned off, we extract the relative error contributions in (d). (d) Error contributions to the gate set for varying Rabi frequencies, normalized to the error at the experimental Ω𝑟𝑎𝑏𝑖/(2𝜋) = 23.7 MHz, magnetic field 2. The attribution to the labeled error sources is presented in the Methods section. (e) Simulated gate fidelity of the √ZZgate. The dashed lines selectively disable error sources as in (c), with only decoherence left for the blue dashed line.  In our optimized setting (magnetic field 2 in Suppl. Table 1 unless stated otherwise, Ω𝑟𝑎𝑏𝑖/(2𝜋) = 23.7 MHz) , we achieve an entangling gate fidelity of 𝐹2𝑞= (96.0 ± 2.5) % as extracted from an error per two-qubit Clifford gate 𝐸𝑃𝐶 = (14.9 ± 2.7) % and pre-characterized average effective single-qubit errors of 𝐹1𝑞 = (99.23 ± 0.12) % (see Equation 9 and Suppl. Note 3). Our entangling gate fidelity is better than the best reported NV-NV entangled state fidelity (82.4 %24) that was optimized by state-to-state transfer optimal control algorithms. We emphasize that the improved performance takes into account errors for arbitrary input states, in contrast to a state-to-state transfer fidelity. As we perform no 10 charge state initialization, the measured SPAM error, represented in the amplitude parameter of the fitted exponential decay, is non-negligible. The decay’s lifetime is substantially shorter than the decoherence limit (𝐸𝑃𝐶/𝐸𝑃𝐶𝑇2 =5.9 ± 1.1). We attribute this to other, coherent errors.  e) Gate error sources Three error sources (explained in detail in Suppl. Note 7) mainly limit our gate set fidelity: The unpolarized 14N nitrogen spins, the misaligned magnetic field and microwave crosstalk and leakage. We measure their influence by performing two-qubit randomized benchmarking while varying the Rabi frequency of all pulses of the gate set. To keep the measurement time viable, we determine only the surviving population 〈𝜎𝑧〉𝑁𝑉1+2 for a fixed Clifford length 𝑛𝑐𝑙𝑖𝑓𝑓 = 2. This allows us to extract relative error contributions, as described in the Methods section, but yields no absolute gate error metric25. Increasing the Rabi frequency drives the hyperfine lines due to the 14N spin more homogenously, but is only beneficial until microwave crosstalk and leakage become limiting. For our magnetic field setting 1 (blue line in Fig. 4c, exact parameters in Suppl. Table 1) with modest frequency separation, the optimal Rabi frequency is ~(2π) 15 MHz. Our optimized magnetic field setting 2 (red line in Fig. 4c) exhibits reduced gate errors due to a large frequency separation that mitigates microwave crosstalk and leakage. Simultaneously, these errors are minimized by employing pulse envelope shaping (here: sine envelope) for all microwave pulses27,28. However, at higher magnetic field, state mixing for the misaligned NV 1 becomes relevant. We estimate the SPAM error introduced solely by spin mixing at the higher magnetic field setting to ~17 % by Rabi measurements on both spins in Suppl. Fig. 5a. We employ a model (Methods section) that accurately describes the experimental randomized benchmarking results and quantifies the influence of different gate error sources. In Fig. 4c, the experimental data at Clifford length 𝑛𝑐𝑙𝑖𝑓𝑓 = 2 is well reproduced by our simulation for both magnetic field settings. Turning off one error source at a time, we can extract the relative error contributions to our gate set for our optimized setting in Fig. 4d. At the experimentally chosen Ω𝑟𝑎𝑏𝑖/(2𝜋) =23.7 MHz, the largest error contribution of 55 % is due to the unpolarized 14N spin. The misaligned magnetic field makes up for 11 % of the observed error with crosstalk and leakage contributing 17 %, demonstrating the efficacy of our pulse envelope shaping. We note that each individual coherent error source investigated in Fig. 4 could, in principle, be corrected by adapting existing control techniques52–54. Their interplay and robustness within the complex four-spin Hamiltonian require further research. The relative decoherence contribution, estimated from 𝑇2,𝑋𝑌8 = 454 μs, 476 μs per qubit, is 10 %.  Finally, we project the achievable entangling gate fidelity directly from our model under a realistic scenario for future experiments in Fig. 4e (97.5 % simulated at our experimental settings, incl. decoherence). Addressing all discussed coherent errors would place our gate fidelity at the coherence limit of 98.6 %. Two-qubit gate fidelities of 99.6 % that support realistic error correction protocols55,56, are achievable by additionally extending the coherence time of the register by higher order dynamical decoupling (see Methods).   III. CONCLUSION AND OUTLOOK We have experimentally demonstrated the highest entangling gate fidelity of 𝐹2𝑞= (96.0 ±2.5) % between solid state electron spins at room temperature. Our analysis of the gate error 11 sources reveals that 90 % of those errors are coherent and correctable. The primary gate error sources are the unpolarized nitrogen spin, which could be addressed by adapting one of the existing nitrogen spin polarization techniques57–60 to our magnetic field regime, and the misaligned magnetic field that is unavoidable for NV geometries with different orientation in the crystal lattice. A perfectly aligned magnetic field, would be possible for NV quantum registers of the same orientation that are conceivable when applying strong magnetic gradient fields61,62 (~ 1 G/nm) or leveraging different nitrogen isotopes of each NV, and it would allow a magnetic field strength compatible with straightforward nitrogen spin polarization59.  The smaller remaining microwave errors (crosstalk and leakage) could be mitigated by optimal control53,63–65. We found that we incur a non-negligible SPAM error due to spin mixing in the misaligned magnetic field. Operating at cryogenic temperatures (~4 K) would significantly decrease both of the current initialization errors. First, resonant laser excitation of the sharp optical absorption lines66 would enable fast, high fidelity charge state initialization67,68. Additionally, spin mixing errors could be avoided for arbitrary magnetic field geometries, as spin initialization would be possible without cycling through the NV singlet state67 and thus independent of the off-axis magnetic field component. Furthermore, operating at low temperatures of ~4 K improves coherence times by  several orders of magnitude, up to seconds30, effectively mitigating the detrimental effect of decoherence on gate fidelities. Last, narrow optical lines would allow to distinguish the emitted photons spectrally and combination with single-shot readout techniques67,69 could enable to measure spin correlations between the NVs. Integrating up to four70 coherently coupled NV centers each with nuclear spin registers of ~25 qubits15–17 is possible with our current experimental technique. Our error analysis suggests that the gate fidelities on such a diamond quantum processor could support error-corrected operation of logical qubits. Efficient error suppression will additionally  necessitate9,12,71,72 fast, high-fidelity optical readout of nuclear67,69 and electron spins73,74.IV. METHODS A. Gate Fidelity Optimization To obtain the experimental parameters realizing our highest gate fidelities, we utilize the following optimization strategy: 1. 𝐵⃗ ̂: We align the magnetic field orientation, such that it approximately yields equally spaced microwave transition frequencies and thus reduced crosstalk and leakage. As a result, the field is closely aligned with one of the NV centers. This method also suppresses interactions with the nitrogen nuclear spin of this NV during the two-qubit gate, as the hyperfine tensor 𝐴 remains diagonal.  We continue by iterative application of the following steps, each of which optimizes 𝐹2𝑞 with respect to one of the remaining gate parameters: 2. 𝜏1: The undesired interaction with the nitrogen spin is strongly pronounced on the more misaligned NV.  In Suppl. Note 4, we show how to select a 𝜏1 that minimizes this undesired interaction during the √𝑍𝑍 gate. However, 𝜏1 > 2𝜏2 reduces the effective dipolar interaction during the gate, leading to increased decoherence. 3. |𝐵⃗ |: We simulate the two-qubit gate fidelity and find the magnetic field strength that maximizes 𝐹2𝑞. Note the trade-off discussed in Suppl. Note 6: Higher |𝐵| can increase gate fidelity but lowers the electron spin initialization fidelity. Our final setting 𝐵⃗ 𝑆.2 is found in Suppl. Table 1. 12 4. Ω𝑟𝑎𝑏𝑖: The same model is used to find the Rabi drive power that maximizes 𝐹2𝑞 as in Fig. 4e. We choose equal power for all microwave pulses in the gate set. Further improvement is possible by optimizing the power for all pulses independently (see Suppl. Fig. 6b).  5.  𝜏2: We experimentally calibrate the phase picked up by dipolar interaction during the √𝑍𝑍 gate, which weakly depends on Ω𝑟𝑎𝑏𝑖, as described in Suppl. Note 4. B. Sample The sample used in this study is a type IIa (100) single-crystalline diamond film that is homoepitaxially grown on a high-pressure high-temperature type Ib substrate via microwave-plasma-assisted CVD75. The CVD film thickness is 20 µm. To suppress the effects of 13C on the coherence properties of the NV centers, a 12C-enriched (99.95 %) high-purity (nitrogen concentration <1 ppb) diamond is used. The ion implantation process of this sample was already described elsewhere31. Ionized C5N4Hn ions are extracted from adenine powder and accelerated to 65 keV kinetic energy. The implantation fluence of 108 cm-2 is achieved after 10 s. After ion implantation, the sample is annealed at 1000 °C for 2 h in a forming gas (4 % H2 in Ar) to create NV centers and recover the diamond lattice. The sample is annealed in an oxygen environment at 465 °C for 4 h, followed by cleaning in a 1:1:1 mixture of HNO3, H2SO4, HClO4 at 200 °C under a pressure of 6-8 bar for 30 min. For the NV pair used as a quantum register, we measured a dephasing (Ramsey) and decoherence (Hahn echo, XY8) times of 𝑇2∗ = (2.60 ± 0.11) μs; (2.37 ± 0.13) μs, 𝑇2,𝐻𝐸 = (55 ± 11) μs; (150 ± 9) μs , 𝑇2,𝑋𝑌8 =(454 ± 58) μs; (476 ± 30) μs on NV 1 and NV 2, respectively. Our pulse spacing definition of type 𝜋2− 𝜏2− 𝜋 −𝜏2−𝜋2  with pulse spacing 𝜏 for the Hahn echo differs by a factor of two from ref. 76 to make it comparable with higher order decoupling sequences. Based on our prior work31, we estimate the creation yield of NV pairs with similar coherence to ≈ 1.7 %. The ratio of NV pairs to single NVs in the sample is ≈ 1:6.5.       C. Setup We control and read the NV’s spin and charge state with a standard home-built confocal microscope. Spin initialization & readout laser pulses of 3000 ns at a wavelength of 552 nm (green) are generated by a cw. laser and an acousto-optic modulator (AOM). The orange (594 nm) charge initialization pulses from a different, digitally modulated cw. laser are 3.5 ms long and of circular polarization to ensure equal ionization rates for both NVs during the charge initialization. The polarization of the green laser is linear and adjusted for nearly equal readout contrast for data in Fig. 3. For Figs. 2 and 4, the chosen polarization realized contrasts as listed in Suppl. Table 1. All laser pulse shapes and microwave waveforms are sampled on an arbitrary waveform generator (Keysight M8195A) amplified and applied to the NV through a copper wire of 20 μm diameter placed on top of the diamond. Photoluminescence photons of the NV in the > 680 nm band are collected through an oil-immersion objective lens (Olympus, 60x, NA 1.35), counted by an avalanche photodiode (APD Excelitas SPCM) during a 330 ns gating window at the beginning of the spin readout pulse77 and digitized by a counting card (FAST ComTec MCS6). The experiment is controlled by a custom measurement software (Qudi78) that features a software interface to quantum algorithms generated by Qiskit79.  To convert experimentally measured fluorescence to polarization 〈𝜎𝑧〉𝑁𝑉1+2, we perform every experiment twice with additional 𝜋 pulses on both NVs. The difference Δ𝑎𝑙𝑡 of such alternating experiments is then normalized to the optical contrast I|00⟩ − I|11⟩ between the fluorescence in state |00⟩ and |11⟩. Our measurement relates to the readout of a register state 𝜌 as 〈𝜎𝑧〉𝑁𝑉1+2  =Δ𝑎𝑙𝑡/(I|00⟩ − I|11⟩) = 𝑅(𝜌) with the readout outcome R as defined in Suppl. Note 2.21 D. Simulations For the simulation results of Fig. 4, we use the following Hamiltonian to describe the system of two NV centers in their orbital ground-state  3A2 : 𝐻(𝑡) = 𝐻𝑁𝑉10 + 𝐻𝑁𝑉20 +𝐻12int +𝐻mw(𝑡) (2) where 𝐻𝑁𝑉10  and 𝐻𝑁𝑉20  are the single NV-Hamiltonians, 𝐻12int  the dipole-dipole coupling Hamiltonian and 𝐻mw the microwave Hamiltonian, that is in general time dependent. The single NV center Hamiltonians80 𝐻0 = 𝐷𝑆𝑧2 + 𝜔⃗⃗ 𝑒 ⋅ 𝑆 + 𝑄𝐼𝑧2 + 𝜔⃗⃗ 𝑛 ⋅ 𝐼 + 𝑆 𝐴𝐼    (3) 13 with  𝜔⃗⃗ 𝑒/𝑛 ≡ −𝛾𝑒/𝑛𝐵⃗  contain the electron zero-field splitting 𝐷, the nuclear quadrupole moment 𝑄, the Zeeman splitting due to the static magnetic field 𝐵⃗ , the gyromagnetic ratio of the electron/nuclear spin in angular frequency units, 𝛾𝑒/𝑛, and the hyperfine coupling to the 14N nuclear spin via the tensor 𝐴. Both the electron spin (initialized in its ground state) and the 14N nucleus (initialized in a maximally mixed state) are triplets: 𝑆 = 𝐼 = 1. The two NV centers have distinct crystal axes (c.f. Suppl. Fig 1a), causing misalignment between the magnetic field and the NV axes and thus a logical basis that is given by rotated eigenstates. For simulating the dynamics of the system, we use an effective dipole-dipole interaction (along the quantization axes) in the 𝐻12int  term that is described in detail in Suppl. Note 2 and list the employed estimates of the system parameters, such as the magnetic field and interaction strength in Suppl. Note 1. Decoherence is accounted for by applying a factor (1 − 𝐸𝑃𝐶𝑇2) in Fig. 4c&d,  (1 −𝐸𝑃𝐺𝑇2) in Fig. 4e and (1 − 𝐸𝑃𝐺𝑇2,1𝑞) for the single-qubit gates in Suppl. Fig. 6b. These factors, calculated in Methods section F, yield approximately constant decoherence contributions to the fidelities in the simulated Rabi frequency range.  In order to investigate the effect of the charge state initialization SPAM error, we perform a numerical simulation, which takes into account the probabilities that each of the NV centers is initially prepared in either the NV-, or NV0 charge state. This simulation employs the simplified Hamiltonian described in Suppl. Note 4 and does not take into account the effect of decoherence. Thus, we include a phenomenological single exponential decay of the observed signal in the charge state simulations in Fig. 3a&c, which is calibrated using the measured pEPGs of Fig. 3b&d. E. Repetitive Benchmarking To convert the measured repetitive benchmarking decays into (pseudo) gate errors, we assume a single exponential decay model that is known to well describe usual 𝑇2,𝑋𝑌8 measurements on single NVs81. We express the decay in terms of the gate repetition number 𝑛 as:  𝑦 = 𝑦0 + 𝑎 ∗ exp (−𝑛/𝑁𝑑) = 𝑦0 + 𝑎 ∗ 𝑝𝑛 (4) where we used the identity 𝑝 = exp (−1𝑁𝑑) = (1 −𝑝𝐸𝑃𝐺) to convert the gate decay parameters 𝑁𝑑 , 𝑝 into a gate error 𝑝𝐸𝑃𝐺. The remaining free fit parameters are the amplitude a and the offset 𝑦0.  Similarly, we obtain the 𝑇2 limit decay curve in Fig. 3 by calculating the polarization loss per applied two-qubit gate of length 𝑡𝑔𝑎𝑡𝑒 = 𝑁𝜋𝜏1 from the mean of the measured 𝑇2,𝑋𝑌8 of both NVs in Fig. 2b:  𝑝𝐸𝑃𝐺𝑇2 = 1 − exp (−𝑡𝑔𝑎𝑡𝑒/((𝑇2,𝑋𝑌8,𝑁𝑉1 +𝑇2,𝑋𝑌8,𝑁𝑉2)/2))  (5) We note that the pEPG derived here from the lifetime parameter of an exponential fit is a useful tool for easily implementable benchmarking, but the results are difficult to compare among experiments (which motivates our labeling as “pseudo” EPG). We empirically find a single exponential decay, which is also the decay model widely used in T2 measurements82. While repetitive benchmarking is similar to T2 experiments - except for the additional evolution by the dipolar coupling - the exact nature of the decay depends on the environment of the probed spins81,83.  F. Randomized Benchmarking We use the same single exponential model of Equation 4 to obtain a decay parameter p for the observed randomized benchmarking decays. If the gate-dependency of the errors is small enough, the measured EPC becomes a good estimate of the average gate set infidelity84.  𝐸𝑃𝐶 =2n − 12n(1 − 𝑝) ≈∑ (1 − 𝐹𝑎𝑣𝑔(𝐶̃𝑖, 𝐶𝑖))𝑖∑ 1𝑖 (6) where n is the number of qubits, {𝐶̃𝑖} are the Clifford gates and {𝐶𝑖} the respective ideal ones (to account for the possibility of gate-dependent errors84), and  𝐹𝑎𝑣𝑔(𝐶̃𝑖 , 𝐶𝑖)= ∫ 𝑑𝜓 𝑇𝑟(𝐶̃𝑖[|𝜓〉〈𝜓|]𝐶𝑖[|𝜓〉〈𝜓|]) (7) 14 is the overlap between the actual process of the i-th Clifford 𝐶̃𝑖 and the ideal 𝐶𝑖. Here, 𝐹𝑎𝑣𝑔 is averaged over all possible initial pure states 𝜓.  Our experiments are generated from the Qiskit software package using the gate set {𝜋𝑥 , (𝜋/2)𝑥, 𝜋𝑦, (𝜋/2)𝑦, 𝐶1𝑁𝑂𝑇2} and 20 random experiments per Clifford. While the single-qubit rotations are easily realized on our NV hardware by microwave pulses with appropriate phases, we need to express 𝐶1𝑁𝑂𝑇2 in terms of our available √𝑍𝑍 gate and multiplication with single-qubit unitaries: 𝐶1𝑁𝑂𝑇2 = 𝜋𝑥,1 ∗ (𝜋/2)−𝑦,2  ∗  √𝑍𝑍  ∗ (𝜋/2)𝑥,1 ∗ (𝜋/2)𝑦,2 ∗ (𝜋/2)𝑦,1 ∗ (𝜋/2)−𝑥,2 ∗ (𝜋/2)𝑥,1 (8) After running Qiskit’s automatic transpiler optimization (“Optimize1qGatesDecomposition”) to reduce the number of single-qubit rotations in the gate string, we end up with a randomized benchmarking experiment that contains on average 𝐺𝑃𝐶1𝑞 = 10.5 single-qubit (𝜋𝑥 , (𝜋/2)𝑥, 𝜋𝑦 , (𝜋/2)𝑦) and 𝐺𝑃𝐶2𝑞 = 1.8 two-qubit (√𝑍𝑍) native gates per Clifford gate. The measured EPC, but not 𝐹2𝑞, could be further improved by finding a more efficient gate set with smaller 𝐺𝑃𝐶1𝑞 . Two-qubit randomized benchmarking yields an EPC that is an average error over the whole gate set. For estimating the error of our √𝑍𝑍 gate, we use the EPC definition79: 𝐸𝑃𝐶 = 1 −∏ (1 − 𝐸𝑃𝐺𝑖)𝐺𝑃𝐶𝑖  𝑖= 1− (1 − 𝐸𝑃𝐺2𝑞)𝐺𝑃𝐶2𝑞(1− 𝐸𝑃𝐺1𝑞)𝐺𝑃𝐶1𝑞 (9) where we collected all single-qubit gate errors in an average gate error 𝐸𝑃𝐺1𝑞 . Solving for 𝐸𝑃𝐺2𝑞  allows us to estimate the two-qubit gate error from the measured two-qubit randomized benchmarking (yielding EPC) and the pre-characterized average effective single-qubit error per Clifford 𝐸𝑃𝐶1𝑞 = 1 − (1 − 𝐸𝑃𝐺1𝑞)𝐺𝑃𝐶1𝑞 (data in Suppl. Fig. 2).  Our randomized benchmarking uses the single exponential model in Equation 4 as a decay model and yields results with low statistical uncertainty, if the tail of the decay curve is well captured. Because of limitations of the applicable microwave power, capturing the decay tail is not always possible in our experiment. For the two-qubit case, we thus determine the offset parameter 𝑦0 of the decay separately first, by applying a high number of intentionally miscalibrated pulses, and fix the decay offset parameter. This amplifies gate errors and thus gives a value that is equally representing the high error limit. We consistently observe 𝑦0 < 0.001 in our experiments. We observe that the randomized benchmarking decay lifetime in Fig. 4b is shorter than expected from repetitive benchmarking (factor 𝐸𝑃𝐶/𝐸𝑃𝐶𝑇2 = 5.9 ± 1.1 , 𝑝𝐸𝑃𝐺/𝑝𝐸𝑃𝐺𝑇2 = 2.5 ± 1.2, coherence limit 𝐸𝑃𝐶𝑇2 derived in Methods section). We attribute this to the fact that in repetitive benchmarking, coherent errors are well refocused by the symmetric timing of the decoupling 𝜋 pulses in the √𝑍𝑍 gate sequence. In combination with random single-qubit Clifford gates in randomized benchmarking, this decoupling works less efficiently.  From simulations we find that a single exponential decay with offset 𝑦0 = 0 is well describing the behavior for experimentally accessible numbers of two-qubit Clifford gates. For 𝑛𝑐𝑙𝑖𝑓𝑓 ≫ 10, accumulated crosstalk and leakage become more relevant, which can cause a slower, 2nd exponential decay85 and depolarization into equal populations of the three 𝑚𝑠 = 0,±1 sublevels.   1. Coherence limit For the T2 coherence limit in Fig. 4, we first estimate the error caused by decoherence per two-qubit gate assuming the same single exponential decay as in the repetitive benchmarking case, thus 𝐸𝑃𝐺𝑇2,2𝑞 = 𝑝𝐸𝑃𝐺𝑇2 = 1.4 % (Equation 5). The decoherence errors from single-qubit gates are much smaller, as they feature shorter gate lengths. Hence, we roughly estimate 𝐸𝑃𝐺𝑇2,1𝑞 by inserting the average single-qubit gate length into Equation 5. Knowing the average number of gates per Clifford, we obtain the coherence limit per Clifford 𝐸𝑃𝐶𝑇2 from Equation 9. Higher order decoupling should enable prolonging T2 towards the longer 2𝑇1limit. We thus estimate in Suppl. Note 8 an achievable coherence limit for the √𝑍𝑍 gate of 𝐹2𝑞 = 99.6 % at room temperature. 15  Relative error contributions We extract the relative error contribution in randomized benchmarking (Fig. 4d) from the simulation as follows: First, we simulate the 〈𝜎𝑧〉𝑁𝑉1+2 outcome of the experiment at 𝑛𝑐𝑙𝑖𝑓𝑓 = 2 with all (coherent) error effects and multiply this value (𝑧𝑠𝑖𝑚) with the known 𝐸𝑃𝐶𝑇2 to account for the decoherence of a single average Clifford gate: 𝑧𝑠𝑖𝑚,𝑇2 = 𝑧𝑠𝑖𝑚(1 − 𝐸𝑃𝐶𝑇2) (10) Using the single exponential decay in Equation 4, we convert the simulated outcome of the randomized benchmarking experiment to an average fidelity parameter p:   𝑝𝑓𝑢𝑙𝑙,𝑇2 = √(𝑧𝑠𝑖𝑚,𝑇2 − 𝑦0)/𝑎𝑛𝑐𝑙𝑖𝑓𝑓 (11) where we take the SPAM parameters 𝑎, 𝑦𝑜 from the measurement in Fig. 4a.  The deviation of our simulation (in terms of p) including all errors from a perfect, gate error free experiment with p=1 is Δ𝑓𝑢𝑙𝑙,𝑇2 =  1 −  𝑝𝑓𝑢𝑙𝑙,𝑇2   Δ𝐸𝑟𝑖,𝑇2 =  1 −  𝑝𝐸𝑟𝑖 ,𝑇2  (12) Similarly, we obtain average fidelity parameters 𝑝𝐸𝑟𝑖,𝑇2  and deviations Δ𝐸𝑟𝑖,𝑇2  for turning off specific error sources 𝐸𝑟𝑖. Then, the relative contribution 𝑐𝑟(𝐸𝑟𝑖) is given by the difference between full simulation and the improved result without this error source 𝑝𝐸𝑟𝑖,𝑇2  divided by the swing of the full experiment: 𝑐𝑟(𝐸𝑟𝑖) =𝑝𝐸𝑟𝑖,𝑇2 − 𝑝𝑓𝑢𝑙𝑙,𝑇2Δfull,T2= 1 − Δ𝐸𝑟𝑖,𝑇2/Δ𝑓𝑢𝑙𝑙,𝑇2   (13)  For example, when the error 𝐸𝑟𝑖 is turned off and the result 𝑝𝐸𝑟𝑖,𝑇2  does not change, the error contribution is 𝑐𝑟(𝐸𝑟𝑖) = 0.  Ideally, we would quantify the error introduced by an off-axis magnetic field by comparing to the application of parallel magnetic fields of different strength on both NVs in simulation. However, the exact transition frequencies that are asymmetrically distributed around the ZFS can only be recovered with a misaligned magnetic field in the simulation. Instead, we note that the experimental outcome on the electron spin with polarized nitrogen spin is equivalent to the case where the hyperfine interaction is turned off completely, because no population transferring terms 𝐴𝑖𝑗𝑆𝑖⨂𝐼𝑗 , 𝑖 ≠ 𝑗 exist in 𝐴 for a perfectly aligned magnetic field. Hence, we can attribute the additional error that occurs for the polarized nitrogen spin case when 𝐴𝑖𝑗𝑆𝑖⨂𝐼𝑗 , 𝑖 ≠ 𝑗 terms are present, solely to the misaligned magnetic field: 𝑐𝑟(𝐵⃗ ⊥ 𝑒 𝑧) = cr(𝐻𝐹𝑆) − c𝑟(unpol. _14N )  (14)  Finally, we can turn off all known error sources in the simulation simultaneously - crosstalk, leakage and the hyperfine interaction - and calculate the difference to an experiment with only decoherence.  𝑐𝑟(𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙) = (1 − 𝑐𝑟(𝑇2))− cr(𝐻𝐹𝑆, 𝑐𝑡. +𝑙𝑒𝑎𝑘. ) (15) These small residual errors scale with the microwave power in a range of 𝑐𝑟(𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙) = [0.9% − 0.09%] for the studied Rabi frequencies between 5-30 MHz. Thus, we attribute them mainly to unwanted dipolar interaction occurring during the microwave pulses. Some part of the gate infidelity cannot be attributed to a single error. This cumulative effect increases in parallel with crosstalk and leakage but has a small contribution (7 %) at the experimental setting.  The alternating readout described above for conversion from fluorescence to spin state can reject some of the leakage error. Although we verified in simulation that this effect is small for our magnetic field setting 2, we base Fig. 4d on a simulation without the additional 𝜋 pulses in order to not underestimate leakage. Our simulation of randomized benchmarking in Fig. 4c does not account for experimental SPAM errors like the imperfect charge state initialization or mixing of the initialized spin state. Such SPAM errors do not change 16 the decay parameter p in randomized benchmarking, but will be represented in the amplitude fit parameters 𝑦0, 𝑎. Hence, we use a free parameter that is multiplied onto our simulation result 〈𝜎𝑧〉𝑁𝑉1+2 and which is determined by least square minimization of the difference to the experimental data.  G. Error bars Error bars on experimental fluorescence data are estimated from the photon shot noise as 1/ √𝑛𝑝ℎ𝑜𝑡 and are error propagated to the spin state domain. Errors on parameters extracted from fits are 1𝜎 uncertainties.  V. DATA AVAILABILITY STATEMENT The datasets in this study are available from the corresponding author on reasonable request. Processed data for the plots is publicly available86. VI. ACKNOWLEDGMENTS We thank Michael Olney-Fraser, Thomas Unden, Philipp Neumann and Alex Retzker for helpful discussions. T.J. acknowledges support by the Foundation of German Business (sdw). This work was funded by the German Federal Ministry of Research (BMBF) through future cluster QSENS to Project No. 03ZU2110DB and projects DE-Brill (No. 13N16207), SPINNING (No. 13N16210, No. 13N16209, and No. 13N16215), CoGeQ (No. 13N16101), DIAQNOS (No. 13N16463), quNV2.0 (No. 13N16707), Quamapolis (No. 13N15375), Quantenrepeater.Net (QR.N), EXTRASENS (No. 13N16935), and QuantumHiFi (No. 16KIS1593); DLR via Project QUASIMODO (No. 50WM2170); and Deutsche Forschungsgemeinschaft (DFG) via Project No. 386028944, No. 445243414, No. 387073854, No. 491245864, No. 546850640, and No. 554644981. This project has also received funding from the European Union’s HORIZON Europe program via projects QuMicro (No. 101046911), SPINUS (No. 101135699), C-QuENS (No. 101135359), QCIRCLE (No. 14 101059999), and FLORIN (No. 101086142); European Research Council (ERC) via Synergy grant HyperQ (No. 856432); and CarlZeiss-Stiftung via the Center of Integrated Quantum Science and technology (IQST) and project Utrasens-Vir (No. P2022-06-007). Additional support was provided by JSPS KAKENHI (No. 21H04646 and No. 24K01286), JST Moonshot R&D Grant No. JPMJMS2062, MEXT Q LEAP Grant No. JPMXS0118067395, JST Adopting Sustainable Partnerships for Innovative Research Ecosystem (ASPIRE) via Grant No. JPMJAP24C1. VII. AUTHOR CONTRIBUTION T.J., P.J.V., M.M.M. and F.J. conceived the experiment. T.T., S.O. fabricated, T.J., R.S., S.O. prepared and characterized the sample. T.J. wrote the code for the experiment and performed the measurements for the presented experimental data. F.F. and G.G. developed and analyzed the model and performed the simulations. T.J., F.F. and G.G. analyzed the data and wrote the manuscript. R.S.S., G.G., J.Z., S.O., T.C., M.M.M. and F.J. supervised the project. All authors reviewed and contributed to the manuscript. VIII. COMPETING INTERESTS The authors declare no competing interests.  IX. REFERENCES 1. Arute, F., Arya, K., Babbush, R., Bacon, D., Bardin, J. C., et al. Quantum supremacy using a programmable superconducting processor. Nature 574, 505–510 (2019). 2. Zhong, H. Sen, Wang, H., Deng, Y. H., Chen, M. C., Peng, L. C., et al. Quantum computational advantage using photons. Science (80-. ). 370, 1460–1463 (2020). 3. Kim, Y., Eddins, A., Anand, S., Wei, K. X., van den Berg, E., Rosenblatt, S., Nayfeh, H., Wu, Y., Zaletel, M., Temme, K. & Kandala, A. Evidence for the utility of quantum computing before fault tolerance. Nature 618, 500–505 (2023). 4. Bluvstein, D., Evered, S. J., Geim, A. A., Li, S. H., Zhou, H., et al. Logical quantum processor based on reconfigurable atom arrays. Nature 626, (2023). 5. Bruzewicz, C. D., Chiaverini, J., McConnell, R. & Sage, J. M. Trapped-ion quantum computing: Progress and challenges. Appl. Phys. Rev. 6, (2019). 6. Smith, K. N., Ravi, G. S., Baker, J. M. & Chong, F. T. Scaling Superconducting Quantum Computers with Chiplet Architectures. Proc. Annu. Int. Symp. Microarchitecture, MICRO 2022-Octob, 1092–1109 (2022). 17 7. Wintersperger, K., Dommert, F., Ehmer, T., Hoursanov, A., Klepsch, J., Mauerer, W., Reuber, G., Strohm, T., Yin, M. & Luber, S. Neutral atom quantum computing hardware: performance and end-user perspective. EPJ Quantum Technol. 10, 1–26 (2023). 8. Wang, J., Sciarrino, F., Laing, A. & Thompson, M. G. Integrated photonic quantum technologies. Nat. Photonics 14, 273–284 (2020). 9. Waldherr, G., Wang, Y., Zaiser, S., Jamali, M., Schulte-Herbrüggen, T., Abe, H., Ohshima, T., Isoya, J., Du, J. F., Neumann, P. & Wrachtrup, J. Quantum error correction in a solid-state hybrid spin register. Nature 506, 204–207 (2014). 10. Kelly, J., Barends, R., Fowler, A. G., Megrant, A., Jeffrey, E., et al. State preservation by repetitive error detection in a superconducting quantum circuit. Nature 519, 66–69 (2015). 11. Chen, Z., Satzinger, K. J., Atalaya, J., Korotkov, A. N., Dunsworth, A., et al. Exponential suppression of bit or phase errors with cyclic error correction. Nature 595, 383–387 (2021). 12. Abobeih, M. H., Wang, Y., Randall, J., Loenen, S. J. H., Bradley, C. E., Markham, M., Twitchen, D. J., Terhal, B. M. & Taminiau, T. H. Fault-tolerant operation of a logical qubit in a diamond quantum processor. Nature 606, 884–889 (2022). 13. Foxen, B., Neill, C., Dunsworth, A., Roushan, P., Chiaro, B., et al. Demonstrating a Continuous Set of Two-Qubit Gates for Near-Term Quantum Algorithms. Phys. Rev. Lett. 125, 120504 (2020). 14. Ballance, C. J., Harty, T. P., Linke, N. M., Sepiol, M. A. & Lucas, D. M. High-Fidelity Quantum Logic Gates Using Trapped-Ion Hyperfine Qubits. Phys. Rev. Lett. 117, 060504 (2016). 15. Bartling, H. P., Yun, J., Schymik, K. N., van Riggelen, M., Enthoven, L. A., van Ommen, H. B., Babaie, M., Sebastiano, F., Markham, M., Twitchen, D. J. & Taminiau, T. H. Universal high-fidelity quantum gates for spin-qubits in diamond. arXiv 2403.10633 (2024) doi:10.48550/arXiv.2403.10633. 16. Bradley, C. E., Randall, J., Abobeih, M. H., Berrevoets, R. C., Degen, M. J., Bakker, M. A., Markham, M., Twitchen, D. J. & Taminiau, T. H. A Ten-Qubit Solid-State Spin Register with Quantum Memory up to One Minute. Phys. Rev. X 9, 031045 (2019). 17. van de Stolpe, G. L., Kwiatkowski, D. P., Bradley, C. E., Randall, J., Abobeih, M. H., Breitweiser, S. A., Bassett, L. C., Markham, M., Twitchen, D. J. & Taminiau, T. H. Mapping a 50-spin-qubit network through correlated sensing. Nat. Commun. 15, 1–9 (2024). 18. Simmons, S. Scalable Fault-Tolerant Quantum Technologies with Silicon Color Centers. PRX Quantum 5, 1 (2024). 19. Pompili, M., Hermans, S. L. N., Baier, S., Beukers, H. K. C., Humphreys, P. C., Schouten, R. N., Vermeulen, R. F. L., Tiggelman, M. J., dos Santos Martins, L., Dirkse, B., Wehner, S. & Hanson, R. Realization of a multinode quantum network of remote solid-state qubits. Science (80-. ). 372, 259–264 (2021). 20. Knaut, C. M., Suleymanzade, A., Wei, Y. C., Assumpcao, D. R., Stas, P. J., Huan, Y. Q., Machielse, B., Knall, E. N., Sutula, M., Baranes, G., Sinclair, N., De-Eknamkul, C., Levonian, D. S., Bhaskar, M. K., Park, H., Lončar, M. & Lukin, M. D. Entanglement of nanophotonic quantum memory nodes in a telecom network. Nature 629, 573–578 (2024). 21. Dolde, F., Jakobi, I., Naydenov, B., Zhao, N., Pezzagna, S., Trautmann, C., Meijer, J., Neumann, P., Jelezko, F. & Wrachtrup, J. Room-temperature entanglement between single defect spins in diamond. Nat. Phys. 9, 139–143 (2013). 22. Degen, M. J., Loenen, S. J. H., Bartling, H. P., Bradley, C. E., Meinsma, A. L., Markham, M., Twitchen, D. J. & Taminiau, T. H. Entanglement of dark electron-nuclear spin defects in diamond. Nat. Commun. 12, (2021). 23. Ungar, A., Cappellaro, P., Cooper, A. & Sun, W. K. C. Control of an Environmental Spin Defect beyond the Coherence Limit of a Central Spin. PRX Quantum 5, 1 (2024). 24. Dolde, F., Bergholm, V., Wang, Y., Jakobi, I., Naydenov, B., Pezzagna, S., Meijer, J., Jelezko, F., Neumann, P., Schulte-Herbrüggen, T., Biamonte, J. & Wrachtrup, J. High-fidelity spin entanglement using optimal control. Nat. Commun. 5, 3371 (2014). 25. Kelly, J., Barends, R., Campbell, B., Chen, Y., Chen, Z., et al. Optimal quantum control using randomized benchmarking. Phys. Rev. Lett. 112, 240504 (2014). 26. Vetter, P. J., Reisser, T., Hirsch, M. G., Calarco, T., Motzoi, F., Jelezko, F. & Müller, M. M. Gate-set evaluation metrics for closed-loop optimal control on nitrogen-vacancy center ensembles in diamond. npj Quantum Inf. 1–12 (2024) doi:10.1038/s41534-024-00893-y. 27. Boradjiev, I. I. & Vitanov, N. V. Control of qubits by shaped pulses of finite duration. Phys. Rev. A - At. Mol. Opt. Phys. 88, 013402 (2013). 28. Mihov, I. S. & Vitanov, N. V. Qubit dynamics driven by smooth pulses of finite duration. Phys. Rev. A 052609, 1–13 (2024). 29. Preskill, J. Quantum computing in the NISQ era and beyond. Quantum 2, 1–20 (2018). 30. Abobeih, M. H., Cramer, J., Bakker, M. A., Kalb, N., Markham, M., Twitchen, D. J. & Taminiau, T. H. One-second coherence for a single electron spin coupled to a multi-qubit nuclear-spin environment. Nat. Commun. 9, 1–8 (2018). 31. Haruyama, M., Onoda, S., Higuchi, T., Kada, W., Chiba, A., Hirano, Y., Teraji, T., Igarashi, R., Kawai, S., Kawarada, H., Ishii, Y., Fukuda, R., Tanii, T., Isoya, J., Ohshima, T. & Hanaizumi, O. Triple nitrogen-vacancy centre fabrication by C5N4Hn ion implantation. Nat. Commun. 10, 4–9 (2019). 32. Waldherr, G., Neumann, P., Huelga, S. F., Jelezko, F. & Wrachtrup, J. Violation of a temporal Bell inequality for single spins in a diamond defect Center. Phys. Rev. Lett. 107, 090401 (2011). 33. Grotz, B., Beck, J., Neumann, P., Naydenov, B., Reuter, R., Reinhard, F., Jelezko, F., Wrachtrup, J., Schweinfurth, D., Sarkar, B. & Hemmer, P. Sensing external spins with nitrogen-vacancy diamond. New J. Phys. 13, (2011). 34. Souza, A. M., Alvarez, G. A. & Suter, D. Experimental protection of quantum gates against decoherence and control errors. Phys. Rev. A - At. Mol. Opt. Phys. 86, 050301(R) (2012). 35. Nielsen, M. & Chuang, I. Quantum Computation and Quantum Information. (Cambridge University Press). 36. Shi, F., Zhang, Q., Wang, P., Sun, H., Wang, J., Rong, X., Chen, M., Ju, C., Reinhard, F., Chen, H., Wrachtrup, J., Wang, J. & Du, J. Single-protein spin resonance spectroscopy under ambient conditions. Science 347, 1135–1138 (2015). 37. Liu, Y. X., Ajoy, A. & Cappellaro, P. Nanoscale Vector dc Magnetometry via Ancilla-Assisted Frequency Up-Conversion. Phys. Rev. Lett. 122, 100501 (2019). 38. Aslam, N., Waldherr, G., Neumann, P., Jelezko, F. & Wrachtrup, J. Photo-induced ionization dynamics of the 18 nitrogen vacancy defect in diamond investigated by single-shot charge state detection. New J. Phys. 15, (2013). 39. Marseglia, L., Hadden, J. P., Stanley-Clarke, A. C., Harrison, J. P., Patton, B., Ho, Y. L. D., Naydenov, B., Jelezko, F., Meijer, J., Dolan, P. R., Smith, J. M., Rarity, J. G. & O’brien, J. L. Nanofabricated solid immersion lenses registered to single emitters in diamond. Appl. Phys. Lett. 98, (2011). 40. Wirtitsch, D., Wachter, G., Reisenbauer, S., Gulka, M., Ivády, V., Jelezko, F., Gali, A., Nesladek, M. & Trupke, M. Exploiting ionization dynamics in the nitrogen vacancy center for rapid, high-contrast spin, and charge state initialization. Phys. Rev. Res. 5, 013014 (2023). 41. Groot-Berning, K., Raatz, N., Dobrinets, I., Lesik, M., Spinicelli, P., Tallaire, A., Achard, J., Jacques, V., Roch, J. F., Zaitsev, A. M., Meijer, J. & Pezzagna, S. Passive charge state control of nitrogen-vacancy centres in diamond using phosphorous and boron doping. Phys. Status Solidi Appl. Mater. Sci. 211, 2268–2273 (2014). 42. Geng, J., Shalomayeva, T., Gryzlova, M., Mukherjee, A., Santonocito, S., Dzhavadzade, D., Dasari, D. B. R., Kato, H., Stöhr, R., Denisenko, A., Mizuochi, N. & Wrachtrup, J. Dopant-assisted stabilization of negatively charged single nitrogen-vacancy centers in phosphorus-doped diamond at low temperatures. npj Quantum Inf. 9, 1–8 (2023). 43. Doi, Y., Fukui, T., Kato, H., Makino, T., Yamasaki, S., Tashima, T., Morishita, H., Miwa, S., Jelezko, F., Suzuki, Y. & Mizuochi, N. Pure negatively charged state of the NV center in n -type diamond. Phys. Rev. B 93, 081203(R) (2016). 44. Souza, A. M., Álvarez, G. a & Suter, D. Robust dynamical decoupling. Philos. Trans. A. Math. Phys. Eng. Sci. 370, 4748–69 (2012). 45. Coote, P., Dimov, R., Maity, S., Hartnett, G. S., Biercuk, M. J. & Baum, Y. Resource-efficient context-aware dynamical decoupling embedding for arbitrary large-scale quantum algorithms. arXiv (2024) doi:10.48550/arXiv.2409.05962. 46. Magesan, E., Gambetta, J. M. & Emerson, J. Scalable and robust randomized benchmarking of quantum processes. Phys. Rev. Lett. 106, 180504 (2011). 47. Chow, J. M., Gambetta, J. M., Córcoles, A. D., Merkel, S. T., Smolin, J. A., Rigetti, C., Poletto, S., Keefe, G. A., Rothwell, M. B., Rozen, J. R., Ketchen, M. B. & Steffen, M. Universal quantum gate set approaching fault-tolerant thresholds with superconducting qubits. Phys. Rev. Lett. 109, 060501 (2012). 48. Huelga, S. F., Macchiavello, C., Pellizzari, T., Ekert, A. K., Plenio, M. B. & Cirac, J. I. Improvement of frequency standards with quantum entanglement. Phys. Rev. Lett. 79, 3865–3868 (1997). 49. Chen, J.-S., Nielsen, E., Ebert, M., Inlek, V., Wright, K., Chaplin, V., Maksymov, A., Páez, E., Poudel, A., Maunz, P. & Gamble, J. Benchmarking a trapped-ion quantum computer with 29 algorithmic qubits. Quantum (2023). 50. Evered, S. J., Bluvstein, D., Kalinowski, M., Ebadi, S., Manovitz, T., Zhou, H., Li, S. H., Geim, A. A., Wang, T. T., Maskara, N., Levine, H., Semeghini, G., Greiner, M., Vuletić, V. & Lukin, M. D. High-fidelity parallel entangling gates on a neutral-atom quantum computer. Nature 622, 268–272 (2023). 51. Marxer, F., Vepsäläinen, A., Jolin, S. W., Tuorila, J., Landra, A., et al. Long-Distance Transmon Coupler with CZ -Gate Fidelity above 99.8 percent. 010314, 1–23 (2023). 52. Motzoi, F., Gambetta, J. M., Rebentrost, P. & Wilhelm, F. K. Simple Pulses for Elimination of Leakage in Weakly Nonlinear Qubits. Phys. Rev. Lett. 103, 1–4 (2009). 53. Müller, M. M., Said, R. S., Jelezko, F., Calarco, T. & Montangero, S. One decade of quantum optimal control in the chopped random basis. Reports Prog. Phys. 85, (2022). 54. Kazi, Z., Shelby, I. M., Watanabe, H., Itoh, K. M., Shutthanandan, V., Wiggins, P. A. & Fu, K. M. C. Wide-Field Dynamic Magnetic Microscopy Using Double-Double Quantum Driving of a Diamond Defect Ensemble. Phys. Rev. Appl. 15, 1 (2021). 55. Wang, D. S., Fowler, A. G. & Hollenberg, L. C. L. Surface code quantum computing with error rates over 1 percent. Phys. Rev. A - At. Mol. Opt. Phys. 83, 020302(R) (2011). 56. Fowler, A. G., Mariantoni, M., Martinis, J. M. & Cleland, A. N. Surface codes: Towards practical large-scale quantum computation. Phys. Rev. A - At. Mol. Opt. Phys. 86, 032324 (2012). 57. Schwartz, I., Scheuer, J., Tratzmiller, B., Müller, S., Chen, Q., Dhand, I., Wang, Z. Y., Müller, C., Naydenov, B., Jelezko, F. & Plenio, M. B. Robust optical polarization of nuclear spin baths using Hamiltonian engineering of nitrogen-vacancy center quantum dynamics. Sci. Adv. 4, 1–8 (2018). 58. Scheuer, J., Schwartz, I., Müller, S., Chen, Q., Dhand, I., Plenio, M. B., Naydenov, B. & Jelezko, F. Robust techniques for polarization and detection of nuclear spin ensembles. Phys. Rev. B 96, 174436 (2017). 59. Jacques, V., Neumann, P., Beck, J., Markham, M., Twitchen, D., Meijer, J., Kaiser, F., Balasubramanian, G., Jelezko, F. & Wrachtrup, J. Dynamic polarization of single nuclear spins by optical pumping of nitrogen-vacancy color centers in diamond at room temperature. Phys. Rev. Lett. 102, 057403 (2009). 60. Chakraborty, T., Zhang, J. & Suter, D. Optimization of a quantum control sequence for initializing a nitrogen-vacancy spin register. Phys. Rev. A 105, 022622 (2022). 61. Bodenstedt, S., Jakobi, I., Michl, J., Gerhardt, I., Neumann, P. & Wrachtrup, J. Nanoscale Spin Manipulation with Pulsed Magnetic Gradient Fields from a Hard Disc Drive Writer. Nano Lett. 18, 5389–5395 (2018). 62. Zhang, H., Arai, K., Belthangady, C., Jaskula, J.-C. & Walsworth, R. L. Selective addressing of solid-state spins at the nanoscale via magnetic resonance frequency encoding. npj Quantum Inf. 3, 1–8 (2017). 63. Werninghaus, M., Egger, D. J., Roy, F., Machnes, S., Wilhelm, F. K. & Filipp, S. Leakage reduction in fast superconducting qubit gates via optimal control. npj Quantum Inf. 7, (2021). 64. Rembold, P., Oshnik, N., Müller, M. M., Montangero, S., Calarco, T. & Neu, E. Introduction to quantum optimal control for quantum sensing with nitrogen-vacancy centers in diamond. AVS Quantum Sci. 2, (2020). 65. Glaser, S. J., Boscain, U., Calarco, T., Koch, C. P., Köckenberger, W., Kosloff, R., Kuprov, I., Luy, B., Schirmer, S., Schulte-Herbrüggen, T., Sugny, D. & Wilhelm, F. K. Training Schrödinger’s cat: Quantum optimal control: Strategic report on current status, visions and goals for research in Europe. Eur. Phys. J. D 69, (2015). 66. Irber, D. M., Poggiali, F., Kong, F., Kieschnick, M., Lühmann, T., Kwiatkowski, D., Meijer, J., Du, J., Shi, F. & Reinhard, F. Robust all-optical single-shot readout of nitrogen-vacancy centers in diamond. Nat. Commun. 12, (2021). 67. Robledo, L., Childress, L., Bernien, H., Hensen, B., 19 Alkemade, P. F. A. & Hanson, R. High-fidelity projective read-out of a solid-state spin quantum register. Nature 477, 574–578 (2011). 68. Pfaff, W., Hensen, B. J., Bernien, H., Van Dam, S. B., Blok, M. S., Taminiau, T. H., Tiggelman, M. J., Schouten, R. N., Markham, M., Twitchen, D. J. & Hanson, R. Unconditional quantum teleportation between distant solid-state quantum bits. Science 345, 532–535 (2014). 69. Neumann, P., Beck, J., Steiner, M., Rempp, F., Fedder, H., Hemmer, P. R., Wrachtrup, J. & Jelezko, F. Single-shot readout of a single nuclear spin. Science (80-. ). 329, 542–544 (2010). 70. Kimura, K., Onoda, S., Yamada, K., Kada, W., Teraji, T., Isoya, J., Hanaizumi, O. & Ohshima, T. Creation of multiple NV centers by phthalocyanine ion implantation. Appl. Phys. Express 15, (2022). 71. Unden, T., Balasubramanian, P., Louzon, D., Vinkler, Y., Plenio, M. B., Markham, M., Twitchen, D., Stacey, A., Lovchinsky, I., Sushkov, A. O., Lukin, M. D., Retzker, A., Naydenov, B., McGuinness, L. P. & Jelezko, F. Quantum Metrology Enhanced by Repetitive Quantum Error Correction. Phys. Rev. Lett. 116, 1–6 (2016). 72. Cramer, J., Kalb, N., Rol, M. A., Hensen, B., Blok, M. S., Markham, M., Twitchen, D. J., Hanson, R. & Taminiau, T. H. Repeated quantum error correction on a continuously encoded qubit by real-time feedback. Nat. Commun. 7, 1–7 (2016). 73. Zhang, Q., Guo, Y., Ji, W., Wang, M., Yin, J., Kong, F., Lin, Y., Yin, C., Shi, F., Wang, Y. & Du, J. High-fidelity single-shot readout of single electron spin in diamond with spin-to-charge conversion. Nat. Commun. 12, 1–6 (2021). 74. Joliffe, M., Vorobyov, V. & Wrachtrup, J. Readout of strongly coupled NV center-pair spin states with deep neural networks. arXiv 1–8 (2024) doi:https://doi.org/10.48550/arXiv.2412.19581. 75. Teraji, T., Yamamoto, T., Watanabe, K., Koide, Y., Isoya, J., Onoda, S., Ohshima, T., Rogers, L. J., Jelezko, F., Neumann, P., Wrachtrup, J. & Koizumi, S. Homoepitaxial diamond film growth: High purity, high crystalline quality, isotopic enrichment, and single color center formation. Phys. Status Solidi Appl. Mater. Sci. 212, 2365–2384 (2015). 76. Hahn, E. L. Spin echoes. Phys. Rev. 80, 580–594 (1950). 77. Suter, D. & Jelezko, F. Single-spin magnetic resonance in the nitrogen-vacancy center of diamond. Progress in Nuclear Magnetic Resonance Spectroscopy vols 98–99 50–62 at https://doi.org/10.1016/j.pnmrs.2016.12.001 (2017). 78. Binder, J. M., Stark, A., Tomek, N., Scheuer, J., Frank, F., Jahnke, K. D., Müller, C., Schmitt, S., Metsch, M. H., Unden, T., Gehring, T., Huck, A., Andersen, U. L., Rogers, L. J. & Jelezko, F. Qudi: A modular python suite for experiment control and data processing. SoftwareX 6, 85–90 (2017). 79. {Qiskit contributors}. Qiskit: An Open-source Framework for Quantum Computing. at https://doi.org/10.5281/zenodo.2573505 (2023). 80. He, X. F., Manson, N. B. & Fisk, P. T. H. Paramagnetic resonance of photoexcited N-V defects in diamond. I. Level anticrossing in the 3A ground state. Phys. Rev. B 47, 8809–8815 (1993). 81. Bar-Gill, N., Pham, L. M., Jarmola, A., Budker, D. & Walsworth, R. L. Solid-state electronic spin coherence time approaching one second. Nat. Commun. 4, 1743 (2013). 82. De Lange, G., Wang, Z. H., Ristè, D., Dobrovitski, V. V. & Hanson, R. Universal dynamical decoupling of a single solid-state spin from a spin bath. Science (80-. ). 330, 60–63 (2010). 83. Bauch, E., Singh, S., Lee, J., Hart, C. A., Schloss, J. M., Turner, M. J., Barry, J. F., Pham, L. M., Bar-Gill, N., Yelin, S. F. & Walsworth, R. L. Decoherence of ensembles of nitrogen-vacancy centers in diamond. Phys. Rev. B 102, 134210 (2020). 84. Proctor, T., Rudinger, K., Young, K., Sarovar, M. & Blume-Kohout, R. What Randomized Benchmarking Actually Measures. Phys. Rev. Lett. 119, 130502 (2017). 85. Wood, C. J. & Gambetta, J. M. Quantification and characterization of leakage errors. Phys. Rev. A 97, 032306 (2018). 86. Joas, T. Plot data. doi:10.5281/zenodo.14825963. 87. Balasubramanian, G., Chan, I. Y., Kolesov, R., Al-Hmoud, M., Tisler, J., Shin, C., Kim, C., Wojcik, A., Hemmer, P. R., Krueger, A., Hanke, T., Leitenstorfer, A., Bratschitsch, R., Jelezko, F. & Wrachtrup, J. Nanoscale imaging magnetometry with diamond spins under ambient conditions. Nature 455, 648–651 (2008). 88. Dolde, F., Fedder, H., Doherty, M. W., Noebauer, T., Rempp, F., Balasubramanian, G., Wolf, T., Reinhard, F., Hollenberg, L. C. L., Jelezko, F., Wrachtrup, J. & Nöbauer, T. Electric-field sensing using single diamond spins. Nat. Phys. 7, 459–463 (2011). 89. Jamonneau, P., Lesik, M., Tetienne, J. P., Alvizu, I., Mayer, L., Dréau, A., Kosen, S., Roch, J. F., Pezzagna, S., Meijer, J., Teraji, T., Kubo, Y., Bertet, P., Maze, J. R. & Jacques, V. Competition between electric field and magnetic field noise in the decoherence of a single spin in diamond. Phys. Rev. B - Condens. Matter Mater. Phys. 93, 024305 (2016). 90. Kehayias, P., Turner, M. J., Trubko, R., Schloss, J. M., Hart, C. A., Wesson, M., Glenn, D. R. & Walsworth, R. L. Imaging crystal stress in diamond using ensembles of nitrogen-vacancy centers. Phys. Rev. B 100, 174103 (2019). 91. Knauer, S., Hadden, J. P. & Rarity, J. G. In-situ measurements of fabrication induced strain in diamond photonic-structures using intrinsic colour centres. npj Quantum Inf. 6, 2–7 (2020). 92. Smeltzer, B., McIntyre, J. & Childress, L. Robust control of individual nuclear spins in diamond. Phys. Rev. A - At. Mol. Opt. Phys. 80, 050302 (2009). 93. Chen, M., Hirose, M. & Cappellaro, P. Measurement of transverse hyperfine interaction by forbidden transitions. Phys. Rev. B - Condens. Matter Mater. Phys. 92, 020101 (2015). 94. Tiesinga, E., Mohr, P. J., Newell, D. B. & Taylor, B. N. Codata recommended values of the fundamental physical constants: 2018. Rev. Mod. Phys. 93, 25010 (2021). 95. Antušek, A., Jackowski, K., Jaszuński, M., Makulski, W. & Wilczek, M. Nuclear magnetic dipole moments from NMR spectra. Chem. Phys. Lett. 411, 111–116 (2005). 96. Goerz, M. H., Reich, D. M. & Koch, C. P. Optimal control theory for a unitary operation under dissipative evolution. New J. Phys. 16, (2014). 97. Vetter, P. J., Marshall, A., Genov, G. T., Weiss, T. F., Striegler, N., Großmann, E. F., Oviedo-Casado, S., Cerrillo, J., Prior, J., Neumann, P. & Jelezko, F. Zero- and Low-Field Sensing with Nitrogen-Vacancy Centers. Phys. Rev. Appl. 17, 044028 (2022). 98. Wang, Z. H., De Lange, G., Ristè, D., Hanson, R. & Dobrovitski, V. V. Comparison of dynamical decoupling protocols for a nitrogen-vacancy center in diamond. Phys. Rev. B - Condens. Matter Mater. Phys. 85, 57–61 (2012). 20 99. Degen, C. L., Reinhard, F. & Cappellaro, P. Quantum sensing. Rev. Mod. Phys. 89, 1–41 (2017). 100. Genov, G. T., Schraft, D., Vitanov, N. V. & Halfmann, T. Arbitrarily Accurate Pulse Sequences for Robust Dynamical Decoupling. Phys. Rev. Lett. 118, 133202 (2017). 101. Suter, D. & Álvarez, G. A. Colloquium: Protecting quantum information against environmental noise. Rev. Mod. Phys. 88, 041001 (2016). 102. Casanova, J., Wang, Z. Y., Haase, J. F. & Plenio, M. B. Robust dynamical decoupling sequences for individual-nuclear-spin addressing. Phys. Rev. A - At. Mol. Opt. Phys. 92, 042304 (2015). 103. Neumann, P., Kolesov, R., Naydenov, B., Beck, J., Rempp, F., Steiner, M., Jacques, V., Balasubramanian, G., Markham, M. L., Twitchen, D. J., Pezzagna, S., Meijer, J., Twamley, J., Jelezko, F. & Wrachtrup, J. Quantum register based on coupled electron spins in a room-temperature solid. Nat. Phys. 6, 249–253 (2010). 104. Yudilevich, D., Stöhr, R., Denisenko, A. & Finkler, A. Mapping Single Electron Spins with Magnetic Tomography. Phys. Rev. Appl. 18, 054016 (2022). 105. Tetienne, J. P., Rondin, L., Spinicelli, P., Chipaux, M., Debuisschert, T., Roch, J. F. & Jacques, V. Magnetic-field-dependent photodynamics of single NV defects in diamond: An application to qualitative all-optical magnetic imaging. New J. Phys. 14, 0–15 (2012). 106. Gupta, A., Hacquebard, L. & Childress, L. Efficient signal processing for time-resolved fluorescence detection of nitrogen-vacancy spins in diamond. J. Opt. Soc. Am. B 33, B28 (2016). 107. Waldherr, G., Beck, J., Steiner, M., Neumann, P., Gali, A., Frauenheim, T. H., Jelezko, F. & Wrachtrup, J. Dark states of single nitrogen-vacancy centers in diamond unraveled by single shot NMR. Phys. Rev. Lett. 106, 157601 (2011). 108. Geller, M. R. & Sun, M. Toward efficient correction of multiqubit measurement errors: Pair correlation method. Quantum Sci. Technol. 6, (2021). 109. Gillespie, D. T. The mathematics of Brownian motion and Johnson noise. Am. J. Phys. 64, 225–240 (1996). 110. Senkalla, K., Genov, G., Metsch, M. H., Siyushev, P. & Jelezko, F. Germanium Vacancy in Diamond Quantum Memory Exceeding 20 ms. Phys. Rev. Lett. 132, 26901 (2024). 111. Souza, A. M., Ávarez, G. A. & Suter, D. Robust dynamical decoupling for quantum computing and quantum memory. Phys. Rev. Lett. 106, 1–4 (2011).   21  X. SUPPLEMENTARY MATERIAL S.I. Magnetic Field Parameters    Suppl. Figure 1. Geometry of the two-NV quantum register. (a) Physical geometry of coupled NVs A and B in the crystal lattice, not to scale. The order in the gate sequence (labeled NV 1 and NV 2) is given along with geometry parameters in Suppl. Table 1. (b) Effective picture of the coupling. In this case, each NV is subject to a (unphysical) different magnetic field vector and the dipolar coupling is of zz type.  Most data presented is measured at the optimized magnetic field settings 2 of Suppl. Table 1. There, we use an aligned magnetic field on the control qubit (NV 2) for minimal interaction with the nitrogen nuclear spins. On the misaligned target qubit, the interaction during the entangling gate is much pronounced, but we can use the 𝜏1 degree of freedom to mitigate the issue (see Suppl. Note 4). In Fig. 4b, we also use the un-optimized setting 1 to illustrate the strong influence microwave crosstalk and leakage can have.  To obtain the magnetic field geometries, we determine the microwave transition frequencies by an ODMR measurement. For each magnetic field setting, the ODMR of two NV centers yields a set of four transition frequencies, where the lowest and highest frequencies form a pair as well as the two frequencies in between, which are associated to a magnetic field of higher misalignment with the NV center axis. Given these transition frequencies and the zero-field splitting 𝐷, one is able to calculate the absolute magnetic field value 𝜔𝑒,𝑖/(2𝜋) (in frequency units) and misalignment with each NV center using87: (𝜔𝑒,𝑖2π)2=13(𝑓1,𝑖2 + 𝑓2,𝑖2 − 𝑓1,𝑖𝑓2,𝑖 − (𝐷𝑖2𝜋)2) − (Ei2π)2   (S1) H𝑖 =7(𝐷𝑖2𝜋)3+ 2(𝑓1,𝑖 + 𝑓2,𝑖) (2(𝑓1,𝑖2 + 𝑓2,𝑖2 ) − 5𝑓1,𝑖𝑓2,𝑖 − 9(𝐸𝑖2𝜋)2) − 3 (𝐷𝑖2𝜋)(𝑓1,𝑖2 + 𝑓2,𝑖2 − 𝑓1,𝑖𝑓2,𝑖 + 9(𝐸𝑖2𝜋)2)9 (𝑓1,𝑖2 + 𝑓2,𝑖2 − 𝑓1,𝑖𝑓2,𝑖 − (𝐷𝑖2𝜋)2− 3(𝐸𝑖2𝜋)2) (S2) cos (2𝜃𝑖) =H𝑖 − 𝐸𝑖/(2π)cos (2𝜑𝑖)𝐷𝑖/(2π) − 𝐸𝑖/(2π)cos (2𝜑𝑖) (S3) where we defined H𝑖 for a clear notation, 𝑓1,𝑖 and 𝑓2,𝑖 are the mentioned pair of transition frequencies, 𝜃𝑖 the misalignment angle between the axis of the considered NV center and the magnetic field, 𝜑𝑖 the azimuthal angle (in the plane perpendicular to the NV axis) and 𝑖 ∈ 𝐴, 𝐵 the NV center to which the quantities are assigned. Note that the geometry in Suppl. Fig 1a is defined in terms of the physical NVs A and B, that can differ from their order in the gate sequence (NV 1 and NV 2).  22 Crystal strain can be divided into transversal strain 𝐸𝑖, and axial strain, which we model by including it into the zero-field parameters 𝐷𝑖88. Distinguishing both components from ODMR-like experimental data is not straightforward if considering a system of two NVs. Thus, we use the following methodology to determine 𝐷𝑖, 𝐸𝑖, making use of the assumption that the (absolute value of the) external magnetic field is constant for the closely located NVs: As 𝐸𝑖  is typically small (≪ 1 MHz) in bulk samples88–90 we are setting it to zero. For the magnetic field setting 2, we observe that Suppl. Equations S1-S3 do not lead to a physical solution when assuming a typical zero-field splitting at room temperature of 𝐷/(2π) = 2870 MHz77. Thus, we obtain the solution 𝐷𝐴/(2π) = 2865.42 MHz that minimizes the absolute value of the axial strain on NV A and yields a valid magnetic field vector. Under the constraint that the absolute magnetic field is the same for both NV centers, this leads to an axial strain for NV B of 𝐷𝐵/(2𝜋) = 2867.27 MHz. Our results for 𝐷𝐴, 𝐷𝐵 are not significantly altered if a (pessimistic estimate) transversal strain of 𝐸/(2𝜋) ∼ 1 MHz 91 is introduced, justifying or initial assumption of 𝐸 = 0. Using the same 𝐷𝐴, 𝐷𝐵, obtained from the analysis of magnetic field setting 2, in the lower magnetic field setting 1, we observe a small difference in the absolute magnetic field of both NVs, |𝜔𝑒,𝐴 −𝜔𝑒,𝐵|/γe =1.1 G.  We simplify our notation and use 𝜃 ≡ 𝜃𝐴 and 𝜑 ≡ 𝜑𝐴 in the following. Given the previously determined magnetic field misalignment 𝜃𝑖, we can calculate the remaining angle 𝜑 in Suppl. Fig 1a by making use of simple geometry. The angle between the direction of the magnetic field 𝜔⃗⃗ ̂𝑒 and the axis of NV center B, 𝑛⃗ ̂𝐵 = (sin𝛽 , 0, cos 𝛽)𝑇, is the misalignment angle 𝜃𝐵 : 𝜔⃗⃗ ̂𝑒 = 𝐵⃗ ̂ = (cos𝜑 sin 𝜃, sin𝜑 sin 𝜃, cos 𝜃)𝑇  cos𝜃𝐵 = 𝜔⃗⃗ ̂𝑒𝑇 ⋅ 𝑛⃗ ̂𝐵 = cos𝜑 sin𝜃 sin 𝛽 + cos 𝜃 cos𝛽 (S4) Note that we can find 𝜑 up to a relative minus sign that cannot be determined in a system of two NV centers only, due to its symmetry. At the same time, as we cannot differentiate between the signs, we do not expect that the dynamics are affected by the choice of the sign. Similarly, the geometry is not uniquely defined for 𝜃, with two solutions 𝜃 and 𝜃𝑖′ = 𝜋 − 𝜃𝑖 for 𝜃 ∈ [0, 𝜋). Furthermore, the angle between the NV center axes can take two possible values 𝛽′ = 109.47∘ and 𝛽 =180∘ − 𝛽′ = 70.53° due to the crystal structure of diamond.        23 Suppl. Table 1. Parameters as chosen and estimated in the simulation. The computational basis of the qubit system (defined as |0⟩ = |𝑔⟩, |1⟩ = |𝑒1/2⟩) employs the new spin eigenstates (due to the tilted quantization axis) that are superpositions of the 𝑚𝑠 = |0⟩, ±|1⟩ electron spin states, with |𝑔⟩ referring to the spin ground state and |𝑒1/2⟩ to the excited spin eigenstates (lower/higher transition frequency than ZFS).  Magnetic field setting Setting 1 |𝐵⃗ 𝑆.1|~60G Setting 2 |𝐵⃗ 𝑆.2|~100G NV center (gate sequence/physical order) NV 1 / NV A NV 2 / NV B NV 2 / NV A NV 1 / NV B Zero-field splitting & axial strain 𝐷/(2𝜋) 2865.42 MHz 2867.27 MHz 2865.42 MHz 2867.27 MHz Transversal strain 𝐸/(2𝜋) 0 MHz Transition frequencies 𝑓𝑖 (*): unused 2932.5 MHz 2829.4 MHz (*) 2751.8 MHz 2999.4 MHz (*) 2571.0 MHz 3160.2 MHz (*) 2990.8 MHz 2827.3 MHz (*) Frequency separation 𝑀𝑖𝑛(|𝑓𝑖 − 𝑓𝑗|) 66.9 MHz 163.5 MHz Absolute magnetic field |𝐵| 64.23 G 63.10 G 105.33 G Absolute magnetic field (electron/nuclear) |𝜔𝑒,𝑖|/(2𝜋) = 𝛾𝑒|𝐵|/(2𝜋) |𝜔𝑛,𝑖|/(2𝜋) = 𝛾𝑛|𝐵|/(2𝜋) 180.01 MHz / 19.76 kHz 176.85 MHz / 19.42 kHz 295.18 MHz / 32.41 kHz Misalignment 𝜃𝑖 73.42° 45.53° 3.58° 74.08° Quadrupole moment 𝑄/(2π) -4.945 MHz 92 14N Off-Axis Hyperfine coupling 𝐴𝑥𝑥/(2π) = 𝐴𝑦𝑦/(2π) -2.62 MHz 93 14N On-Axis Hyperfine coupling 𝐴𝑧𝑧/(2π) -2.162 MHz 92 Angle between NV centers 𝛽 70.53° Magnetic field orientation phi 𝜑 47.94° 172.73° Effective dipole-dipole coupling strength 𝜈𝑑𝑖𝑝 = 𝑔/(2𝜋) 0.13995 MHz 0.11261 MHz Computational basis of qubit system |0⟩, |1⟩ |𝑔⟩, |𝑒2⟩ |𝑔⟩, |𝑒1⟩ |𝑔⟩, |𝑒1⟩ |𝑔⟩, |𝑒2⟩ Readout contrast 𝛼𝑖 14.6 % 14.6 % 19.4 % 10.7 % Gyromagnetic ratio (electron/nuclear) 𝛾𝑒/(2𝜋)  𝛾𝑛/(2𝜋) -28.02495 GHz/T 94  3.076272 MHz/T 95 24 S.II. Simulation of Two Interacting NV Centers  a) System Hamiltonian We start with the Hamiltonian80 of a single, negatively charged, NV center with an inherent nitrogen 14N nucleus in its orbital ground-state 3A2: 𝐻0 = 𝐷𝑆𝑧2 + 𝜔⃗⃗ 𝑒 ⋅ 𝑆 ⏟        ≡𝐻e+𝑄𝐼𝑧2 + 𝜔⃗⃗ 𝑛 ⋅ 𝐼 ⏟        ≡𝐻n+𝑆 𝐴𝐼 ⏟  ≡𝐻hfc (S5) where 𝐷 is the electronic zero-field splitting, 𝑄 the nuclear quadrupole moment, 𝐴 the hyperfine tensor, which is diagonal, 𝜔⃗⃗ 𝑒/𝑛 ≡ −𝛾𝑒/𝑛𝐵⃗  with 𝛾𝑒/𝑛 being the gyromagnetic ratio of the electron/nuclear spin in angular frequency units (Hz rad/T) and 𝐵⃗  the magnetic field at the positions of the NV center. In the chosen coordinate system, the z-axis is aligned with the NV center axis. Throughout this work, we use Hamiltonians in angular frequency units, i.e. ℏ = 1.  We call 𝐻𝑒/n the electron/nuclear Hamiltonian and 𝐻hfc  the hyperfine coupling Hamiltonian. The respective spin operators 𝑆𝑖(𝐼𝑖), with 𝑖 ∈ {𝑥, 𝑦, 𝑧}, are defined to act on the electron (nuclear) states of the NV center only. E.g., 𝑆𝑧 is equivalent to 𝑆𝑧⊗𝟙, i.e., a Kronecker product of 𝑆𝑧 for NV1 and the identity operator 𝟙 for the nucleus. Both the electron spin, due to the charge state NV−of the NV center, and the nuclear spin, since we consider the 14N isotope, are triplets (𝑆 = 𝐼 = 1). The system under consideration consists of two NV centers with distinct crystal axes (Suppl. Fig 1a sketches magnetic field setting 2 (Suppl. Table 1), in which NV 1 is more strongly misaligned), interacting with each other via dipole-dipole coupling, that are exposed to a static magnetic field and microwave pulses of arbitrary shape. Note that this section is written with respect to magnetic field setting 2. The labels NV 1 and NV 2 that indicate the order in the gate sequence, would need to be switched for magnetic field setting 1. The free Hamiltonian of the system is given by: 𝐻free = 𝐻𝑁𝑉10 +𝐻𝑁𝑉20 +𝐻12𝑖𝑛𝑡 (S6) with the dipole-dipole interaction Hamiltonian 𝐻12int  𝐻12int = 𝑔0(𝑟12)[𝑆 1 ⋅ 𝑆 2 − 3(𝑆 1 ⋅ 𝑟 ̂12)(𝑆 2 ⋅ 𝑟 ̂12)] = 𝑔0(𝑟12)𝑆 1𝐺(𝑟 ̂12)𝑆 2′  (S7) where we defined 𝑆 𝑘 as the operator 𝑆  acting on the k-th NV center (𝑘 = 1,2) and its components 𝑆𝑖,k (𝑖 =𝑥, 𝑦, 𝑧). The rotation of the spin operator 𝑆 2′  of NV2 is detailed below.   We denote the displacement of NV center 1 from NV center 2 as 𝑟 12 and the absolute coupling strength as 𝑔0(𝑟12), depending on the distance 𝑟12 between the NV centers: 𝑔0(𝑟12) ≡𝜇04𝜋ℏ𝛾𝑒2𝑟123  (S8) We call 𝐺(𝑟 ̂12) the geometric tensor. It does not depend on the distance 𝑟12, but on the normalized relative position vector 𝑟 ̂12 (c.f. Suppl. Fig. 1a).  25   The two NV centers have misaligned axes, i.e. the crystal axis of NV 1 is rotated by an angle 𝛽 around the y-axis with respect to the axis NV 2 (c.f. Suppl. Fig. 1a). Resultingly, the Hamiltonian takes an unhandy shape, e.g., the zero-field splitting term reads 𝐷(cos (𝛽)𝑆𝑧 + sin (𝛽)𝑆𝑥)2. We solve this problem, by transforming the frame of NV 1: We rotate the spin operators and fields of NV 1 around the y-axis by an angle −𝛽 into the 𝑥′, 𝑦′, 𝑧′-frame where the axis of NV 1 coincides with the 𝑧′-axis (c.f. Suppl. Fig. 1a): 𝜔⃗⃗ ′ = 𝑅𝑦(−𝛽) 𝜔⃗⃗ , 𝑆 ′ = 𝑅𝑦(−𝛽) 𝑆 , 𝑆 ≡ (𝑆𝑥𝑆𝑦𝑆𝑧) , 𝑅𝑦(−𝛽) = (cos 𝛽 0 −sin 𝛽0 1 0sin 𝛽 0 cos 𝛽) (S9) In the primed frame, the Hamiltonian of NV 1 takes the again the usual form as given in Suppl. Equation S5, i.e. the zero-field splitting term reads 𝐷(𝑆𝑧′)2. In this picture the NV centers experience a different magnetic field and a modified geometric tensor 𝐺(𝑟 ̂12). As the NV centers are exposed to a constant magnetic field that is (generally) misaligned with their crystal axes, their quantization axes are changed (in length and direction), which results in new eigenvalues and eigenstates (|𝑔⟩, |𝑒1⟩,|𝑒2⟩) of the Hamiltonian. As a consequence, the eigenstates of the electron Hamiltonian 𝐻𝑒 do not coincide with the eigenstates of the 𝑆𝑧,2, 𝑆𝑧,1′  operators, but we can define new spin operators 𝑆̃𝑧,2, 𝑆̃𝑧,1′  that have the same eigenstates as 𝐻e and take the usual (matrix) form in the corresponding basis. Similarly, we define 𝑆̃𝑥/𝑦,2, 𝑆̃𝑥/𝑦,1′ . The dipole-dipole interaction is (typically) small compared to the electron transition frequencies (𝑔0 ≪ 𝜔i), so that a secular approximation can be applied, whereby small, non-diagonal elements (in the basis of the 𝑆̃𝑧,𝑖 eigenstates) in the dipole-dipole interaction Hamiltonian are neglected. This simplifies the dipole-dipole coupling to: 𝐻12int = 𝑔(𝑟 12, 𝛽, 𝐵⃗ )𝑆̃𝑧,2𝑆̃𝑧,1′   (S10) where 𝑔(𝑟 12, 𝛽, 𝐵⃗ ) is the effective coupling strength (𝑣𝑑𝑖𝑝 = 𝑔/(2𝜋) in frequency units). By transforming the frame of NV center 1, transforming the spin operators into the basis of the electron eigenstates and performing the secular approximation, we have developed an effective picture of the two NV center system (c.f. Suppl. Fig. 1b). The total (noiseless) Hamiltonian is given by: 𝐻𝑑𝑟𝑖𝑣𝑒𝑛(𝑡) = 𝐻free +𝐻mw(𝑡) (S11)  where 𝐻𝑁𝑉10  and 𝐻𝑁𝑉20  are the single NV center Hamiltonians of NV centers 1 and 2, and 𝐻mw the microwave Hamiltonian 𝐻mw(𝑡) = ∑  𝑖=1,2√2[𝛺𝑖,𝑥(𝑡) cos(𝜔𝑖𝑡 + 𝜉𝑖) + 𝛺𝑖,𝑦(𝑡) sin(𝜔𝑖𝑡 + 𝜉𝑖)] 𝐻control  (S12) 26   where 𝛺𝑖,𝑥(𝑡), 𝛺𝑖,𝑦(𝑡) are time-dependent amplitudes (or control functions), that vary only slowly compared to the angular carrier frequencies 𝜔𝑖. Phases of the carrier, generally set to zero, are denoted 𝜉𝑖 and 𝐻control  is the control Hamiltonian 𝐻control ≡ 𝑆̃𝑥,1′ + 𝑆̃𝑥,2 + 𝛾̃( 𝐼𝑥,1 + 𝐼𝑥,2) (S13)  with 𝛾̃ ≡ 𝛾𝑛/𝛾𝑒 = 𝜔𝑛/𝜔𝑒. b) Simulation For simulating the time dynamics of the system, we start with single NV Hamiltonians 𝐻1/20 , with the rotation by angle −𝛽 applied on NV center 1, as above and compute the matrix 𝑇 that diagonalizes the sum of the electron Hamiltonians 𝐻𝑒 ≡ 𝐻1𝑒 +𝐻2𝑒, such that 𝑇†𝐻𝑒𝑇 is diagonal. We apply the diagonalization matrix to the sum of the single NV Hamiltonians 𝐻10 + 𝐻20 and then add the effective dipole-dipole coupling term in Suppl. Equation S8 to obtain the (effective) free Hamiltonian, 𝐻free = 𝑇†(𝐻10 +𝐻20)𝑇 + 𝐻12int (S14) From here on, we use the basis of the electronic eigenstates, so that 𝐻12int = 𝑔 𝑆̃𝑧,2𝑆̃𝑧,1′  as above. The spin operators are defined as 𝑆̃𝑖,2 = 𝑆𝑖⊗𝟙⊗ 𝟙⊗ 𝟙, 𝑆̃𝑖,1′ = 𝟙⊗ 𝟙⊗ 𝑆𝑖⊗ 𝟙 𝐼𝑖,2 = 𝟙⊗ 𝐼𝑖⊗ 𝟙⊗ 𝟙, 𝐼𝑖,1 = 𝟙⊗ 𝟙⊗ 𝟙⊗ 𝐼𝑖 (S15) With order (NV2, nuclear spin 2, NV1, nuclear spin 1) and 𝑆𝑥 =1√2(0 1 01 0 10 1 0) = 𝐼𝑥, 𝑆𝑦 =𝑖√2(0 −1 01 0 −10 1 0) = 𝐼𝑦, 𝑆𝑧 = (1 0 00 0 00 0 −1) = 𝐼𝑧  (S16) If the magnetic field is static, the Hamiltonian of the system is given by 𝐻free  and therefore time-independent. In this case, the time evolution operator is calculated via the matrix exponential 𝑈(𝑡) =exp (−𝑖𝐻free 𝑡) and then transformed into the rotating frame via: 𝑈rot (𝑡) = 𝑉†(𝑡)𝑈(𝑡)   with  𝑉(𝑡) = exp (−𝑖𝐻trans 𝑡) (S17)  where 𝐻trans ≡ 𝜔1(𝑆̃𝑧,1′ )2+ 𝜔2(𝑆̃𝑧,2)2 is the Hamiltonian that generates the rotating frame transformation, 𝜔1/2 = 2𝜋𝑓1/2 are the angular carrier frequencies, ideally coinciding with the transition lines 2𝜋𝑓i of the chosen qubit system (Suppl. Table 1). 27  If microwave pulses are applied, the Hamiltonian in Suppl. Equation S11 contains time-dependent terms. In order to calculate the time evolution from 𝑡0 to 𝑡, we approximate the time-evolution via the Riemann summation method, whereby we evaluate the Hamiltonian at 𝑁 discrete time points: 𝑈(𝑡, 𝑡0) = 𝒯exp (−i∫  𝑡𝑡0 𝐻rot(𝑡′)d𝑡′) ≈ 𝑒−𝑖𝐻rot(𝑡𝑁)𝜏…𝑒−𝑖𝐻rot(𝑡2)𝜏𝑒−𝑖𝐻rot (𝑡1)𝜏 (S18)  with 𝑡𝑗 = 𝑡0 + 𝑗 𝜏 and 𝜏 = (𝑡 − 𝑡0)/𝑁, where 𝑁 = 100 GHz (𝑡 − 𝑡0) and 𝒯 is the time ordering operator.  To realize faster convergence, i.e. realize the same precision with larger time-steps, the Hamiltonian is transformed into the rotating frame (at every timestep) in which the Hamiltonian varies only slowly: 𝐻rot(𝑡𝑖) = 𝑉𝑖†(𝐻driven (𝑡𝑖) − 𝐻trans )𝑉𝑖  with  𝑉𝑖 = 𝑉(𝑡𝑖) (S19) When dealing with sequences of pulses and free evolution, we multiply the respective unitaries time-orderly.  Readout results 𝑅(𝜌) are obtained via the POVM element 𝐸𝑔({𝛼𝑖}), similar to the definition in 21, as: 𝑅(𝜌) = Tr𝑒 [𝐸𝑔Tr𝑛 (𝑈𝜌𝑈†)] − Tr𝑒 [𝐸𝑔Tr𝑛 (𝑋12𝑈𝜌𝑈†𝑋12† )] 𝐸𝑔({𝛼𝑖}) =1𝛼1 + 𝛼2(𝛼1|𝑔⟩⟨𝑔|⊗ 𝟙 + 𝛼2𝟙⊗|𝑔⟩⟨𝑔|)  (S20) where 𝑋12 is a 𝜋𝑥-pulse on NV center 1, followed by a 𝜋𝑥-pulse on NV center 2, Tr𝑒/𝑛 is the trace of the electron/nuclear spins, and 𝛼1/2 is a contrast parameter, proportional to the spin contrast observed in experiment. As we only read out diagonal elements in the measurement, the results are the same in the rotating frame as in the lab frame (since 𝑉𝑖 is diagonal). Within the simulation, if not stated otherwise, the electron state is perfectly initialized in the ground state, the nuclear state in the maximally mixed state (thermal state approximation). If the nuclear state is stated to be initialized, then it is initialized into its 𝑚𝐼 = 0 state. If not stated otherwise, the pulses are chosen to have a sinusoidal shape. E.g., for a 𝜋𝑥-pulse on NV center 1, the control functions are given by: 𝛺1,𝑥(𝑡) = Ω𝑚𝑎𝑥𝑠𝑖𝑛(𝜋 ⋅ 𝑡/𝑡𝜋),  𝛺1,𝑦(𝑡) = 𝛺2,𝑥(𝑡) = 𝛺2,𝑦(𝑡) = 0 (S21) where Ω𝑚𝑎𝑥 =𝜋2Ω𝑟𝑎𝑏𝑖 is the peak amplitude of the envelope, Ω𝑟𝑎𝑏𝑖 =𝜋𝑡𝜋, and 𝑡𝜋 is the 𝜋 pulse duration of the sine shaped pulse. For evaluations of the gate fidelity, we use the definition of the average fidelity, given in 96: 28 𝐹avg =1𝑑(𝑑 + 1)∑  𝑑𝑖,𝑗=1  (⟨𝑖|𝐺†𝒟(|𝑖⟩⟨𝑗|)𝐺|𝑗⟩ + ⟨𝑖|𝐺†𝒟(|𝑗⟩⟨𝑗|)𝐺|𝑖⟩) (S22)  where 𝐺 is the target quantum gate and 𝒟 is the dynamical map, i.e. the time evolution of the electron spins (under the influence of the nuclear spins): 𝒟(𝜌) = Tr𝑛 (𝑈𝜌𝑈†) (S23) In our case, 𝑑 = 4 and the {|𝑖⟩} states are the basis of the qubit subspace. We evaluate the gate-fidelity in the rotating frame. The randomized benchmarking circuits, used in the simulation, are the same as in the corresponding experiments. When sweeping the Rabi frequency, the set of pulses and the free evolutions are calculated for each Rabi frequency individually. When crosstalk and leakage are stated to be turned off, this means that the control Hamiltonian is modified to only act on the qubit levels of the desired NV center.   29 S.III. Single-Qubit Randomized Benchmarking To estimate the entangling gate fidelity 𝐸𝑃𝐺2𝑞 from the measured EPC of the gate set, we need to separate the error contribution of the single-qubit gates, as described in the Methods section: 𝐸𝑃𝐶 = 1 − (1 − 𝐸𝑃𝐺2𝑞)𝐺𝑃𝐶2𝑞(1 − 𝐸𝑃𝐺1𝑞)𝐺𝑃𝐶1𝑞 (S24)  This is often done by measuring single-qubit randomized benchmarking on the sub systems (Suppl. Fig. 2a). Due to a particularity of our system, we take a slightly different approach here. We first note that due to limitations of our microwave hardware, we never apply pulses on both NVs at the same time. As a result of the unpolarized nitrogen spins, a single qubit rotation on NV1 thus can cause an unwanted phase pickup during the free evolution on NV2 and vice versa. This effect is neglected when carrying our single-qubit randomized benchmarking on a single NV only. To obtain a one-qubit fidelity that represents all errors that are present when using the two-qubit register, we determine an average effective single-qubit error 𝐸𝑃𝐶1𝑞 = 1 − (1 − 𝐸𝑃𝐺1𝑞)𝐺𝑃𝐶1𝑞 from a modified two-qubit randomized benchmarking experiment in Suppl. Fig. 2b: This approach excludes all entangling gates from the generated experiment (𝐺𝑃𝐶2𝑞 = 0) and appends single-qubit correction pulses to revert to the ground state (resulting modified 𝐺𝑃𝐶1𝑞 = 11.5). On every qubit, the modified experiment retains a random sequence of single-qubit Clifford operations – except for the free evolution that we treat as an error. We then extract the average effective single-qubit error  EPC1q from the fitted decay lifetime parameter 𝑝 =exp (−1𝑁𝑑) as: 𝐸𝑃𝐶1𝑞 = (2n − 1)/2𝑛(1 − 𝑝) = 0.085 ± 0.013 (S25)  With n the number of qubits in the Clifford sequence, here 𝑛 = 2. This translates to an effective single-qubit gate fidelity of 𝐹1𝑞=(99.23 ± 0.12) %.   Suppl. Figure 2. Single-qubit randomized benchmarking. (a) Given the gate set GPC1q = 1.97, we extract bare single-qubit gate fidelities of 𝐹1𝑞,𝑁𝑉1 = (99.49 ± 0.56)% and 𝐹1𝑞,𝑁𝑉2 = (99.53 ± 0.49)%. Offset parameters are fixed by a separate measurement, as described in the Methods section. (b) Average single-qubit error 𝐹1𝑞 = (99.23 ± 0.12)% from modified two-qubit randomized benchmarking with stripped entangling gate.   30 S.IV. Entangling Gate Dynamics & Calibration  Suppl. Figure 3. Entangling gate calibration. (a) XY8-1 experiment performed with pulse spacing τ1 on the misaligned target NV1 at magnetic field settings 2. (b) Measured repeated application of the √ZZ gate with varying τ2 (τ1=800 ns, nrep = 4, Ω𝑟𝑎𝑏𝑖/(2𝜋) = 22.97 MHz) for calibration. The projection pulse after the √ZZ gates for both NVs is a π𝑥/2. The solid line is a fit from which we extract the τ2,√ZZ that realizes our entangling gate. In order to calibrate the duration of our entangling gate, we first analyze how it scales with the duration of our decoupling sequence and the choice of 𝜏1 and 𝜏2.  To this end, we simplify the effective model Hamiltonian (Suppl. Equations S6, S11) further by treating the effect of the nitrogen nuclear spin as effective detuning (similarly to 97) and arrive at: 𝐻 = 𝐻𝑁𝑉10 +𝐻𝑁𝑉20 +𝐻mw +𝐻12int  𝐻𝑁𝑉10 = D𝑆𝑧,12 + Δ1𝑆𝑧,1 𝐻𝑁𝑉20 = D𝑆𝑧,22 + Δ2𝑆𝑧,2 𝐻mw = √2 {Ω1,𝑥(𝑡) cos(𝜔1𝑡 + 𝜉1) + Ω2,𝑥(𝑡) cos(𝜔2𝑡 + 𝜉2)}(𝑆𝑥,1 +𝑆𝑥,2) 𝐻12int = 𝑔𝑆𝑧,2𝑆𝑧,1  (S26) Here, we are considering the effective dipolar coupling 𝑔 = 2𝜋𝜈𝑑𝑖𝑝 between the NV centers98, with the zero-field splitting D of the NVs, Δ𝑘 , 𝑘 = 1,2 is the effective detuning due to the Zeeman and hyperfine splittings, whose effect is typically refocused after the gate, Ω1,𝑥(𝑡) and Ω2,𝑥(𝑡) characterize the Rabi frequencies of the microwave fields with angular frequencies 𝜔1 and 𝜔2 and initial relative phases 𝜉1 and 𝜉2 that we use to drive the two NV centers during the refocusing pulses. We use here the operator labels 𝑆𝑧,1 and 𝑆𝑧,2 for 𝑆̃𝑧,1′  and 𝑆̃𝑧,2 from section 2b, respectively, for simplicity of presentation and without loss of generality. We move to the rotating frame, defined by 𝜔1𝑆𝑧,12  and 𝜔2𝑆𝑧,22 , apply the rotating wave-approximation, using that 𝛺1,2(𝑡) ≪  𝜔1,2, and obtain 31 𝐻𝑟𝑜𝑡𝑁𝑉1 = (D − ω1)𝑆𝑧,12 + Δ1𝑆𝑧,1 +𝛺1(𝑡)√2𝑆𝑥,1 =(   D− ω1 + Δ1𝛺1(𝑡)20𝛺1(𝑡)20𝛺1(𝑡)20𝛺1(𝑡)2D − ω1 − Δ1)   =(   𝛿1𝛺1(𝑡)20𝛺1(𝑡)20𝛺1(𝑡)20𝛺1(𝑡)2δδ1 )     (S27)  𝐻𝑟𝑜𝑡𝑁𝑉2 = (D − ω2)𝑆𝑧,22 + Δ2𝑆𝑧,2 +𝛺2(𝑡)√2𝑆𝑥,2 =(   D− ω2 + Δ2𝛺2(𝑡)20𝛺2(𝑡)20𝛺2(𝑡)20𝛺2(𝑡)2D − ω2 − Δ2)   =(   δδ2𝛺2(𝑡)20𝛺2(𝑡)20𝛺2(𝑡)20𝛺2(𝑡)2𝛿2 )    (S28)  𝐻𝑟𝑜𝑡𝑡𝑜𝑡𝑎𝑙 = 𝐻𝑟𝑜𝑡𝑁𝑉1 + 𝐻𝑟𝑜𝑡𝑁𝑉2 + 𝑔 𝑆𝑧,2𝑆𝑧,1  =(                𝛿1 + δδ2 + 𝑔𝛺1(𝑡)20𝛺2(𝑡)20 0 0 0 0𝛺1(𝑡)2δδ2𝛺1(𝑡)20𝛺2(𝑡)20 0 0 00𝛺1(𝑡)2δδ1 + δδ2 − 𝑔 0 0𝛺2(𝑡)20 0 0𝛺2(𝑡)20 0 𝛿1𝛺1(𝑡)20𝛺2(𝑡)20 00𝛺2(𝑡)20𝛺1(𝑡)20𝛺1(𝑡)20𝛺2(𝑡)200 0𝛺2(𝑡)20𝛺1(𝑡)2δδ1 0 0𝛺2(𝑡)20 0 0𝛺2(𝑡)20 0 𝛿1 + 𝛿2 − 𝑔𝛺1(𝑡)200 0 0 0𝛺2(𝑡)20𝛺1(𝑡)2𝛿2𝛺1(𝑡)20 0 0 0 0𝛺2(𝑡)20𝛺1(𝑡)2𝛿2 + δδ1 + 𝑔⁠)                  (S29)  where we assume that the driving field frequency 𝜔1(𝜔2) is highly detuned from the transition frequencies of NV2 (NV1), so there is negligible crosstalk from pulses applied to NV1 on NV2 and vice versa. In addition, we take 𝜉1 = 𝜉2 = 0 for simplicity of presentation and without loss of generality. The matrix representation of 𝐻𝑟𝑜𝑡𝑁𝑉1, 𝐻𝑟𝑜𝑡𝑁𝑉2 is in the single qutrit basis of each of the NV centers, while 𝐻𝑟𝑜𝑡𝑡𝑜𝑡𝑎𝑙 is represented in the two-qutrit basis. In the experiment, the first driving field is approximately resonant with the |0〉 ↔ |+1〉 transition on NV1 (𝛿1 ≡ D −ω1 + Δ1 ≪ Ω1(𝑡), 𝛿𝛿1 ≡ D −ω1 − Δ1 ≫ Ω1(𝑡)) and the       32 second field - with the |0〉 ↔ |−1〉 transition on NV2 (𝛿2 ≡ D −ω2 − Δ2 ≪ Ω2(𝑡), 𝛿𝛿2 ≡ D −ω2 +Δ2 ≫ Ω2(𝑡)). The Hamiltonian elements, which characterize the evolution of the states, which we can effectively address with the microwave fields (they do not have a very high frequency offset in the diagonal element of the Hamiltonian) are highlighted with dashed lines. Thus, the respective, reduced Hamiltonian, which characterizes the evolution that leads to our two-qubit gate is given by 𝐻𝑟𝑜𝑡,𝑟𝑒𝑑𝑢𝑐𝑒𝑑𝑡𝑜𝑡𝑎𝑙 =(      𝛿1𝛺1(𝑡)2𝛺2(𝑡)20𝛺1(𝑡)20 0𝛺2(𝑡)2𝛺2(𝑡)20 𝛿1 + 𝛿2 − 𝑔𝛺1(𝑡)20𝛺2(𝑡)2𝛺1(𝑡)2𝛿2⁠)        (S30)  The detunings 𝛿1 and 𝛿2 typically cause dephasing but can be refocused by applying dynamical decoupling. The z-z coupling 𝑔 allows us to apply the √𝑍𝑍 gate, as evident from the analysis below. The Hamiltonian during free evolution without microwave pulses and the respective propagator are given by 𝐻𝑓𝑟𝑒𝑒 = (𝛿1 0 0 00 0 0 00 0 𝛿1 + 𝛿2 − 𝑔 00 0 0 𝛿2⁠) ,  𝑈𝑓𝑟𝑒𝑒(𝑇𝑓) = exp(−𝑖 𝐻𝑓𝑟𝑒𝑒𝑇𝑓) = (𝑒−𝑖𝑇𝑓𝛿1 0 0 00 1 0 00 0 𝑒𝑖𝑇𝑓(−𝛿1−𝛿2+𝑔) 00 0 0 𝑒−𝑖𝛿2𝑇𝑓⁠), (S31)  where 𝑇𝑓 is the respective free evolution time. We assume that the Rabi frequencies are much stronger than 𝑔, 𝛿1 and 𝛿2, so we can neglect their effect during the refocusing pulses. Then, the Hamiltonians and propagators of the π pulses on NV 1 and NV 2 are respectively 𝐻𝑝𝑢𝑙𝑠𝑒𝑁𝑉1 ≈(      0𝛺1(𝑡)20 0𝛺1(𝑡)20 0 00 0 0𝛺1(𝑡)20 0𝛺1(𝑡)20⁠)      ,  𝑈1(𝜃1) = exp (−𝑖 𝐻𝑝𝑢𝑙𝑠𝑒𝑁𝑉1 𝑇1) ≈(      cos (𝜃12) −𝑖 sin (𝜃12) 0 0−𝑖 sin (𝜃12) cos (𝜃12) 0 00 0 cos (𝜃12) −𝑖 sin (𝜃12)0 0 −𝑖 sin (𝜃12) cos (𝜃12)⁠)        (S32)  𝐻𝑝𝑢𝑙𝑠𝑒𝑁𝑉2 ≈(      0 0𝛺2(𝑡)200 0 0𝛺2(𝑡)2𝛺2(𝑡)20 0 00𝛺2(𝑡)20 0⁠)      ,   𝑈2(𝜃2) = exp(−𝑖 𝐻𝑝𝑢𝑙𝑠𝑒𝑁𝑉2 𝑇2) ≈(      cos (𝜃22) 0 −𝑖 sin (𝜃22) 00 cos (𝜃22) 0 −𝑖 sin (𝜃22)−𝑖 sin (𝜃22) 0 cos (𝜃22) 00 −𝑖 sin (𝜃22) 0 cos (𝜃22)⁠)         (S33) where 𝜃1 = 𝛺1𝑇1 = 𝜋 and 𝜃2 = 𝛺2𝑇2 = 𝜋 (we assumed that the pulses are rectangular for simplicity of presentation and without loss of generality).  In order to characterize the gate, we consider a sequence of two instantaneous 𝜋 pulses on each NV and analyze the evolution in the interaction basis of the pulses. The Hamiltonians in each free evolution period are given by: Period 1: before the first π pulse on NV1: 33 𝐻𝑖𝑛𝑡,1 = 𝐻𝑓𝑟𝑒𝑒 = ( ⁠𝛿1 0 0 00 0 0 00 0 −𝑔 + 𝛿1 + 𝛿2 00 0 0 𝛿2⁠) , 𝑡1 =τ12,  (S34)  Period 2: after one π pulse on NV1, no π pulses on NV2:  𝐻𝑖𝑛𝑡,2 =  𝑈𝑖𝑛𝑡,2 𝐻𝑓𝑟𝑒𝑒 𝑈𝑖𝑛𝑡,2†= (0 0 0 00 𝛿1 0 00 0 𝛿2 00 0 0 −𝑔 + 𝛿1 + 𝛿2⁠) ,where  𝑈𝑖𝑛𝑡,2 =  𝑈1(𝜋), 𝑡2 =τ12− 𝜏2 (S35)  Period 3: after one π pulse on NV1, one π pulse on NV2:  𝐻𝑖𝑛𝑡,3 =  𝑈𝑖𝑛𝑡,3 𝐻𝑓𝑟𝑒𝑒 𝑈𝑖𝑛𝑡,3†= ( ⁠𝛿2 0 0 00 −𝑔 + 𝛿1 + 𝛿2 0 00 0 0 00 0 0 𝛿1⁠) ,where  𝑈𝑖𝑛𝑡,3 =   𝑈2(𝜋)𝑈1(𝜋), 𝑡3 =τ12+ 𝜏2 (S36)  Period 4: after two π pulses on NV1, one π pulse on NV2:  𝐻𝑖𝑛𝑡,4 =  𝑈𝑖𝑛𝑡,4 𝐻𝑓𝑟𝑒𝑒 𝑈𝑖𝑛𝑡,4†= ( ⁠−𝑔 + 𝛿1 + 𝛿2 0 0 00 𝛿2 0 00 0 𝛿1 00 0 0 0⁠) ,where  𝑈𝑖𝑛𝑡,4 =   𝑈1(𝜋) 𝑈2(𝜋)𝑈1(𝜋), 𝑡4 =τ12− 𝜏2 (S37)  Period 5: after two π pulses on NV1, two π pulses on NV2:  𝐻𝑖𝑛𝑡,5 =  𝑈𝑖𝑛𝑡,5 𝐻𝑓𝑟𝑒𝑒 𝑈𝑖𝑛𝑡,5†= ( ⁠𝛿1 0 0 00 0 0 00 0 −𝑔 + 𝛿1 + 𝛿2 00 0 0 𝛿2⁠) ,where  𝑈𝑖𝑛𝑡,5 =   𝑈2(𝜋)𝑈1(𝜋) 𝑈2(𝜋)𝑈1(𝜋), 𝑡5 = 𝜏2 (S38)  The total evolution in the interaction basis is then given by 𝑈𝑖𝑛𝑡 = exp(−𝑖∑𝐻𝑖𝑛𝑡,𝑘𝑡𝑘) = exp(i 𝜉g)( ⁠1 0 0 00 𝑒𝑖𝑁𝜋𝑔𝜏2 0 00 0 𝑒𝑖𝑁𝜋𝑔𝜏2 00 0 0 1⁠) (S39)   where we used that [𝐻𝑖𝑛𝑡,𝑖, 𝐻𝑖𝑛𝑡,𝑗] = 0 to obtain the first equality, 𝜉g = 𝑁𝜋[𝑔2(−𝜏12+ 𝜏2) +𝜏12(𝛿1 + 𝛿2)] is an irrelevant global phase and 𝑁𝜋 is the total number of applied pulses on each NV center, e.g., 𝑁𝜋 = 2 in our example with two pulses or 𝑁𝜋 = 8 for XY8-1 and 𝑁𝜋 = 16 for XY8-2. Thus, in the following, we define the effective evolution time 𝑡𝑒𝑣𝑜𝑙 = 𝑁𝜋𝜏2.  For demonstrating the controlled oscillations of the √𝑍𝑍 gate with XY8-1 decoupling on the Bloch sphere equator in Fig. 2, we use 𝜏1 = 3000 ns. This choice allows to sweep the effective evolution time 𝑡𝑒𝑣𝑜𝑙 =𝑁𝜋𝜏2 with large dynamic range (constrained by −𝜏1/2 ≤ 𝜏2 ≤ 𝜏1/2), but is not an optimal choice in terms of fidelity for two reasons: 34 First, shorter 𝜏1 allow to realize the same gate unitary with shorter total sequence length (= 𝑁𝜋𝜏1, independent of 𝜏2), and thus less decoherence. Consequently, for realizing the √𝑍𝑍 unitary, 𝜏1 should be chosen close to 𝜏1,𝑚𝑖𝑛 =24𝑁𝜋𝑣𝑑𝑖𝑝~520 ns for our 𝑣𝑑𝑖𝑝 = 120 kHz, 𝑁𝜋 = 8 when using XY8-1. We note that using shorter 𝜏1 for our entangling gate is also possible by increasing 𝑁𝜋 in the dynamical decoupling sequence, e. g., using XY8-2 that consists of 16 pulses, as long as 𝑡𝑒𝑣𝑜𝑙 = 𝑁𝜋𝜏2 with 𝜏2 ≤ 𝜏1/2 is kept constant. Decreasing the pulse spacing 𝜏1 should improve fidelity as the characteristic 𝑇2 time of a decoupling sequence typically increases when the pulse separation is shorter99–101.  Second, due to the misaligned magnetic field on the target qubit, population transfer to and entanglement with the nitrogen nuclear spin can occur. Since we can choose 𝜏1 freely (but not 𝜏2 that is set by 𝜈𝑑𝑖𝑝), we use this degree of freedom to minimize such an unwanted interaction with the nitrogen spin. In Suppl. Fig. 3a, we probe the polarization loss due to interaction with the nitrogen by initializing the target qubit into superposition and applying a standard XY8-1 experiment only on the target qubit with varying 𝜋  pulse spacing 𝜏1. We observe that the simulation is well reproducing the experimental data and find a 𝜏1 = 800 ns that reduces interaction with the nitrogen spin. In future, modified gate microwave sequence could minimize this interaction with the nucleus more strictly102. After determining 𝜏1, our gate calibration for 𝜏2 is carried out as follows: We initialize the input state (|0⟩ − 𝑖|1⟩⨂|0⟩ − 𝑖|1⟩) and apply the √𝑍𝑍 gate sequence, with varying 𝜏2, four times (𝑛 = 4). Repeating the gate improves the accuracy of the calibration, as more oscillation periods are recorded for a 𝜏2 sweep. To the experimentally recorded oscillation (eg. in Suppl. Fig. 3b), we fit a sine and find the 𝜏2,√𝑍𝑍 that is realizing minimal fluorescence, as expected from a quarter rotation (repeated 𝑛 = 4 times) on the Bloch sphere equator. In order to capture the gate dynamic accurately, we repeat the same calibration procedure to obtain a coupling parameter for our model. In addition to the experiment in in Suppl. Fig. 3b at magnetic field setting 2, we experimentally perform a calibration also at magnetic field setting 1 (τ1=800 ns, 𝑛 = 4, Ω𝑟𝑎𝑏𝑖/(2𝜋) = 15.51 MHz, data not shown) and obtain a calibrated 𝜏2,√𝑍𝑍,𝐵1 = 212.6 ns. We then simulate both calibration experiments individually and find the effective dipolar coupling parameters 𝑣𝑑𝑖𝑝 in Suppl. Table 1 that reproduce the experimentally observed 𝜏2,√𝑍𝑍,𝐵1 , 𝜏2,√𝑍𝑍,𝐵2. The dipole-dipole coupling of 119.8 kHz sets an upper limit of (9.54 ± 0.03) nm on the distance between the NV centers103. An exact distance estimate would require knowledge the relative orientation of the dipoles and could be obtained by repeated DEER measurements in a varying bias magnetic field104.   35 S.V. Charge Initialization & Readout Our charge initialization is based on a weak orange laser pulse and post-selection of the recorded spin readout photons. After the green (552 nm) spin initialization laser pulse, the charge state of a single NV in bulk diamond reaches a steady state occupation of typically 𝐴(𝑁𝑉−)/(𝐴(𝑁𝑉0) + 𝐴(𝑁𝑉−))~0.7 38, where 𝐴 denotes the probability to find the NV in a certain charge state. The orange laser probes the NV- absorption with minimal overlap to the NV0 absorption spectrum. At the same time, its weak power avoids ionization. Thus, more readout photons are expected during the orange laser pulse if both of the NVs are in the negative charge state. In Suppl. Fig. 4a, we measure readout photons per 3.5 ms long orange laser on our coupled NV quantum register after spin initialization with a 3 μs green laser. Due to the limited collection efficiency, the three expected Poissonian peaks {(𝑁𝑉0, 𝑁𝑉0); (𝑁𝑉0, 𝑁𝑉−); (𝑁𝑉−, 𝑁𝑉−)} overlap in the histogram. For determination of the charge state, we fit a Poissonian model to the photon probability 𝑝(𝑛𝑝ℎ𝑜𝑡) with free weighting parameters 𝐴𝑖 and fitted mean rates 𝜆𝑖: 𝑝(𝑛𝑝ℎ𝑜𝑡) = 𝐴−−𝑃(𝜆−−) + 𝐴−0𝑃(𝜆−0) + 𝐴00𝑃(𝜆00) (S40)   The extracted weights of the charge state slightly differ from the values 𝐴00 = 0.09,  𝐴−0 = 0.42, 𝐴−− =0.49 expected for an independent combination of two bulk NVs. We speculate that a single charge shared between both NVs is more favorable in our system. In Suppl. Fig. 4b we append another 2.9 ms long orange laser pulse to independently probe the charge state fidelity after the 𝑡𝑝𝑐𝑠 = 3.5 ms long charge initialization laser pulse, without applying post selection yet. As expected, we recover charge state weights that are in agreement with the data in Suppl. Fig. 4a.  We then continue in Suppl. Fig. 4c to analyze only charge readout data for experimental shots where the counted photons during the first initialization orange laser pulse satisfy a threshold photon 𝑛𝑝ℎ𝑜𝑡 > 𝑛𝑡ℎ𝑟𝑒𝑠ℎ. After applying post-selection with the experimental settings (𝑡𝑝𝑐𝑠 = 3.5 ms, 𝑛𝑡ℎ𝑟𝑒𝑠ℎ=9) presented in Fig. 2,3, we measure an increase in (𝑁𝑉−, 𝑁𝑉−) charge state initialization fidelity from (39 ± 6) % (no post-selection) to (83 ± 6) %.  In Suppl. Fig. 4d&e, we characterize the trade-off between higher charge state fidelity and readout noise: Improved charge initialization can be reached for shorter 𝑡𝑝𝑐𝑠 or higher 𝑛𝑡ℎ𝑟𝑒𝑠ℎ, but as more photons are discarded by the thresholding, the photon shot noise of the readout estimated from the total number of detected photons 𝑛𝑝ℎ𝑜𝑡𝑠 as SNR𝑠ℎ𝑜𝑡−1 = 1/√𝑛𝑝ℎ𝑜𝑡𝑠 increases. We tune the charge state by varying 𝑛𝑡ℎ𝑟𝑒𝑠ℎ in Suppl. Fig. 4f. Here, we utilize the 𝜎𝑦 component of the DEER oscillation as an alternative probe of the charge state and observe that its asymmetry reduces to (0.12 ± 0.12) when increasing to 𝑛𝑡ℎ𝑟𝑒𝑠ℎ = 10, indicating a clearly improved charge state initialization.   36  Suppl. Figure 4. Charge initialization & readout. (a) Histogram of photons collected during the 𝑡𝑝𝑐𝑠=3.5 ms orange charge initialization pulse. Before, the register is initialized with a 3 μs green laser pulse. Dashed lines are the single Poissonian contributions with fitted mean rates 𝜆𝑖 and weight 𝐴𝑖. The chosen thresholding photon number 𝑛𝑡ℎ𝑟𝑒𝑠ℎ = 9 used for charge initialization is shown as red dotted line. (b) Photon histogram collected during the 2.9 ms charge readout laser pulse after the charge initialization laser. (c) Photon histogram with same data as (b), but after applying post-selection to analyze only charge readout data if more than 𝑛𝑡ℎ𝑟𝑒𝑠ℎ = 9 photons are detected during the charge initialization pulse. The rates 𝜆𝑖 are fixed during fitting to the values extracted in (b). We conservatively estimate the error on the charge state fidelity by giving the uncertainty of the weight parameter 𝐴−− in (b). (d) NV-/NV- charge initialization fidelity measured for different length 𝑡𝑝𝑐𝑠 of the gating window which counts photons during the charge initialization laser (same data set as in (b)&(c)). A shorter 𝑡𝑝𝑐𝑠 shifts the Poissonian in the histogram to the left and thus increases charge state fidelity. Unstable fits are discarded. (e) Photon shot noise as calculated from the number of analyzed photons (1/√𝑛𝑝ℎ𝑜𝑡𝑠) and normalized to the noise without post-selection.  Increasing the charge state fidelity comes at the cost of higher readout noise. (f) Asymmetry of the 𝜎𝑦 component of the DEER oscillation presented in Fig 2e while varying the post selection threshold parameter 𝑛𝑡ℎ𝑟𝑒𝑠ℎ.     37 S.VI. Spam and Entangling Gate Fidelity in a Magnetic Field Misaligning the magnetic field against the NV axis in the diamond lattice lowers the spin initialization fidelity by green laser excitation and yields a partly mixed spin state. This effect originates from the altered effective decay rates into the new magnetic ground state levels that are superpositions of the aligned spin ground states105. Here, we estimate the SPAM error given in the main text caused by a reduced spin initialization. In our (misaligned) magnetic field setting 2, a reduced readout contrast is observed for the misaligned NV 1 in Suppl. Fig. 5a. We avoid the influence of different orientations of both NVs’ optical dipoles by determining the maximum readout contrast for varying orientation of the linearly polarized green (552 nm) excitation laser by a 𝜆/2 wave plate.  In Suppl. Fig. 5b, we verify that transformation of a 7 level rate equation model105 agrees reasonably well with the experimentally observed relative Rabi contrast max(𝑐𝑜𝑛𝑡𝑟𝑁𝑉1) /max(𝑐𝑜𝑛𝑡𝑟𝑁𝑉2) at different absolute magnetic field values. The misalignment angle 𝜃 ≈ 74° of each datapoint is similar as in magnetic field setting 2, up to a precision of ±2°.  We analyze the model with two sets of rates (see Suppl. Table 2): 1. Optical rates as given in 106. To calibrate the laser pump rate as the only free parameter, we minimize the difference to the experimental data. Comparing to our experimental data, we observe a deviation from the model. 2. Since the intersystem crossing rates are not straightforward to determine, we additionally adapt the rates from the 3E level to the singlet and minimize the difference to experimental data in parallel to the laser pump rate. The resulting rates fit better to our experimental data.   The same model allows us to estimate the SPAM error introduced by a reduced spin polarization. The electron spin initialization fidelity 𝐹𝑖𝑛𝑖𝑡 is given by: 𝐹𝑖𝑛𝑖𝑡 =𝑝0𝑝0 + 𝑝−1 + 𝑝+1  (S41)  Note that the experimentally observed readout contrast at misaligned magnetic field in Suppl. Fig. 5a&b includes contributions from both the reduced spin polarization and the less efficient optical readout. To obtain the spin initialization fidelity, we determine the populations 𝑝𝑚𝑠 of the 𝑚𝑠 = 0,±1 levels after a 3 μs green laser and a waiting time of 1 μs.  In Suppl. Fig. 5c, we plot the simulated spin initialization fidelity from Suppl. Equation S41 and the entangling gate fidelity (predicted like in Fig. 4e) for varying magnetic field. At zero field (equivalent to no misalignment), the predicted spin initialization fidelity (defined in the spin 1 manifold, mean of both set of rates) is 77 %. This value is lower than in experimental work (≥ 92 % 107,(88 ± 4) %103), indicating that the available models might underestimate the absolute degree of spin polarization.  We estimate the SPAM error due to spin mixing to 17 % by dividing the simulated value (mean of both set of rates) at our magnetic field setting 2 by the value at zero field:  𝑒𝑟𝑟𝑆𝑃𝐴𝑀,𝑖𝑛𝑖𝑡 =1 − 𝐹𝑖𝑛𝑖𝑡(𝐵)1 − 𝐹𝑖𝑛𝑖𝑡(𝐵 = 0) (S42)  38 For small NV quantum registers, such a SPAM error by spin mixing seems acceptable, especially when a higher gate fidelity can be reached in turn. On the other hand, larger scale quantum processors will require to mitigate state mixing, as state preparation errors are only inefficiently correctable108.  We observe a trade-off between higher gate fidelity but lower spin initialization fidelity when increasing the magnetic field. Depending on the application, diamond quantum registers with NVs of different orientations may be optimized towards higher entangling gate fidelity or reduced SPAM spin initialization error.    Suppl. Figure 5. SPAM error by spin mixing. (a) Optical contrast as measured from fits to Rabi oscillations on both NVs at different 𝜆/2 wave plate angles with magnetic field setting 2. The maximum optical contrasts of both NVs is determined from the maximum value of a fitted empirical model (sum of two sines). (b) Relative Rabi contrast max(𝑐𝑜𝑛𝑡𝑟𝑁𝑉1) /max(𝑐𝑜𝑛𝑡𝑟𝑁𝑉2) at different absolute values of the magnetic field and misalignment angles 𝜃𝑁𝑉1~74° , experimental precision ± 2°.  The solid line is a simulation based on the rate model and transition rates in Suppl. Table 2. Circles are measurement data taken at different magnetic fields as described in (a). (c) Simulated √𝑍𝑍 gate fidelity and spin initialization fidelity for different absolute values of the magnetic field and a misalignment angles 𝜃𝑁𝑉1 = 74.08° , 𝜃𝑁𝑉2 = 3.58°.    Suppl. Table 2: Rates used for modelling (as in 105) the readout contrast and spin initialization fidelity in a misaligned magnetic field. Rates Model Gupta106, NV 1 Model Gupta, adapted ISC k57 k67 k47 k71 k72 k73 k41 k52 k63 laser pump 𝛽 92.9 MHz 92.9 MHz 11.2 MHz 4.9 MHz 2.03 MHz 2.03 MHz 66.08 MHz 66.08 MHz 66.08 MHz 1.215 90.307 MHz 90.307 MHz 3.004 MHz 4.9 MHz 2.03 MHz 2.03 MHz 66.08 MHz 66.08 MHz 66.08 MHz 1.938    39 S.VII. Simulated Gate Fidelities In the main text, we discuss the error sources determining the gate set fidelity. In Suppl. Fig. 6a, these errors are illustrated in the level structure of the coupled NV system. After spin initialization, both NVs are polarized into the 𝑚𝑠 = 0 sublevels. Due to the unpolarized nuclear 14N spins at ambient conditions, the initialized state is a mixed state with nearly equal classical probability of the nitrogen populations. However, our quantum register is not formed from the 14N levels and initially the electron states are pure (after taking the partial trace).  After and during preparation of arbitrary register input states, multiple gate errors occur:  • As the inherent 14N spins of both NVs are unpolarized, the nuclear spin state is random at the beginning of every experimental shot. Consequently, most of the microwave pulses suffer from a detuning Δ = ±𝐴𝑧𝑧/(2𝜋) =  ±2.16 MHz from the correct microwave transition frequency for the electron, if the nitrogen state is 𝑚𝐼 = ±1. This leads to imperfect driving of the three hyperfine lines, as indicated by the orange sketch of the pulse spectrum. Additionally, the unpolarized nitrogen spin causes free evolution on the non-addressed NV in our single qubit gate implementation (that we discuss in Suppl. Note 3). • In a magnetic field setting that is aligned to the NV quantization axis, the diagonal hyperfine interaction tensor 𝐴 = ∑ 𝐴𝑖𝑖𝑆𝑖⨂𝐼𝑖𝑖=𝑥,𝑦,𝑧  in the Hamiltonian (see Methods, Suppl. Note 2) effectively only contains 𝑆𝑧⨂𝐼𝑧 terms. Here, 𝑆𝑖, 𝐼𝑗 denote the 𝑆 = 𝐼 = 1 spin operator of the NV and the 14N spin, respectively. In our geometry however, magnetic field misalignment is unavoidable and in the tilted basis the hyperfine tensor is no longer diagonal. Then, terms of, e.g., the form 𝑆𝑧⨂𝐼𝑥 appear and can cause entanglement with the nitrogen spin or population transfer from the qubit to the nitrogen spin. Especially during gate operations that contain dynamical decoupling pulses, this decreases the electron spin polarization (see Suppl. Note 4).  • All microwave pulses (driving 𝑚𝑠 = 0,-1 on NV 2, 𝑚𝑠 = 0,1 on NV 1) generate a small, unwanted microwave drive on the other NV (“crosstalk”, dashed arrows in Suppl. Fig. 6a) or transitions to states out of the qubit subspaces (“leakage”, solid arrows).   Note that all error sources may not only change populations, but also cause unwanted phase shifts. In Fig. 4e, we predict the entangling gate fidelity of only the √𝑍𝑍 gate from our model. We use the same simulation to give gate fidelities for the single-qubit gates of the gate set in Suppl. Fig. 6b. As the single-qubit pulses in our experiment don’t employ dynamical decoupling, the fidelity improvement expected from polarized 14N nitrogen spins is more pronounced. At high enough Rabi frequencies, we expect nearly no further improvement by aligning the magnetic field. Consequently, we anticipate that polarizing the nitrogen spins and mitigating crosstalk and leakage by optimal control could enable very high-quality single-qubit control in multi-NV diamond quantum registers.  40  Suppl. Figure 6. Gate errors. (a) Level scheme of the 3A2 level of two NVs including Zeeman and 14N hyperfine interaction. The error sources discussed in the main text are sketched. After spin initialization by the green laser, we obtain pure electron states indicated by blue dashed circles after tracing out the 14N nuclear spins. During a gate, interaction with the thermal nitrogen spin causes depolarization of the electron spin (green dashed ellipse). (b) Simulated gate fidelities for all single-qubit gates (of phase X) of the gate set for different Rabi frequencies. Decoherence is considered as T2 process as described in the Methods section and much smaller for the shorter single-qubit gates.    41 S.VIII. Extended Coherence Limits Using Equation 5 with coherence times for different decoupling sequences, we estimate the coherence-limited entangling gate fidelity in realistic experimental scenarios, also cited in the Methods section F, in Suppl. Table 2. Specifically, we note that decoherence in NV centers and other color centers in diamond can usually be modelled with an Ornstein-Uhlenbeck process81,109,110 and one can estimate the expected memory time limit with ideal, instantaneous π pulses for particular pulse spacing 𝜏1 by110 T2(𝜏1) = (𝑇2,𝐻𝐸𝜏1)2𝑇2,𝐻𝐸    (S43) not taking 𝑇1 decay into account. This results in much higher theoretical values 𝑇2,XY8 of the order of several ms for the particular 𝜏1 = 800 ns, which is then limited by 𝑇1 relaxation. The lower actual 𝑇2,XY8 in our experiment suggest that pulse errors might have an effect on the achieved decoherence times, so using robust dynamical decoupling sequences, e.g., KDD4 or the UR family100,111 is likely to improve 𝑇2. The improvement itself can be characterized in terms of the transition probability error ϵ of a single π pulse, so the fidelity of a single XY8 sequence is110 F𝑋𝑌8 ≈ 1 − 𝑂(ϵ3)  (S44) In comparison, the fidelity of KDD4 (KDD sequence, nested in the XY4 sequence)111 and the UR sequences are100,111  F𝐾𝐷𝐷4 ≈ 1 − 𝑂(ϵ5) (S45) F𝑈𝑅𝑛 ≈ 1 − 𝑂(ϵ𝑛/2) (S46)   Thus, using these advanced robust sequences could, in principle, allow to improve our coherence times towards the 2𝑇1 limit, which should be of the order of milliseconds. Following the fidelity estimation in Equation 5, the decoherence limit of XY8 can be estimated to F2𝑞,𝑋𝑌8 = 1 − 𝑝𝐸𝑃𝐺𝑇2 ≈ 0.986, with the decoherence time for our experiment. Improving the coherence time towards the milliseconds regime would in principle push the expected fidelity towards higher values, e.g., F2𝑞,𝑎𝑑𝑣 ≈ exp(−𝑡𝑔𝑎𝑡𝑒/𝑇2,𝑎𝑑𝑣) ≈0.996 for 𝑇2,𝑎𝑑𝑣 ≈ 2 ms with advanced dynamical decoupling sequences. We note that the assumption of simple exponential decay is conservative as the actual decay rate for dynamical decoupling with an Ornstein-Uhlenbeck noise model is typically a stretched exponential109,110, which would result in even higher fidelities on the time scale of the gate. Another approach for extending the coherence limit is operating at low temperatures, which can be beneficial in several aspects. First, 𝑇1 relaxation due to two-phonon Raman and Orbach-type processes is negligible at temperatures of the order of 4 K, boosting 𝑇1 to the order of one hour30. This results in several orders of improvement in the achievable coherence times e.g., 𝑇2,𝐶𝑃,4𝐾 ≈ 1.58 s30. This would in principle make the effect of decoherence negligible on the time scales of realistic NV-NV dipolar coupling, yielding 1 − F2𝑞,𝑇2,4𝐾 ≈ 4 × 10−6. Apart from negligible decoherence, operating at low temperatures in principle would allow the polarization of the nitrogen nuclear spin, reducing the need for high Rabi frequency to cover all hyperfine transitions. Operating with lower Rabi frequencies will in turn reduce the effect of crosstalk and leakage, minimizing the detrimental effect of both of these errors, as reported in Fig. 4d in the main text.    42 Suppl. Table 2. 𝑇2 coherence times in various experimental scenarios. Assuming the coherent error sources discussed in the context of Fig. 4 are mitigated, the listed decoherence-limited entangling gate fidelities are achievable. For the extended coherence times, the same gate duration 𝑡𝑔𝑎𝑡𝑒 as in the current experiment is assumed (†). The duration 𝑡𝑔𝑎𝑡𝑒 of our two-qubit gate depends on the effective coupling, and considers the optimized 𝜏1 = 800 ns.  Scenario Relevant coherence time 𝑡𝑔𝑎𝑡𝑒 Decoherence-limited gate fidelity 𝐹2𝑞,𝑇2 Previous best, Room temperature, Hahn echo 24 avg(𝑇2,HE) = 332 μs 25.4 μs 0.926 Current experiment, Room temperature, Hahn echo avg(𝑇2,HE) = 103 μs  6.4 μs 0.939 Current experiment, Room temperature, XY8 avg(𝑇2,XY8) = 465 μs 6.4 μs 0.986 Advanced DD, Room temperature avg(𝑇2,aDD) =  2 ms 6.4 μs † 0.996 Carr-Purcell, Cryogenic experiment (4 K) avg(𝑇2,CP) = 1.58 s 30 6.4 μs † 0.999996  To compare our gate design with previous work, we also list the decoherence limits of a two-qubit gate under Hahn echo decoupling for both experiments in Suppl. Table 2. Due to its Hahn echo-like structure, 𝑇2,𝐻𝐸  is the relevant coherence time for the gate in ref. 24.  We observe a slight improvement in the decoherence-limited gate fidelity for our sample, with an improvement factor of (1 − 𝐹2𝑞,𝑇2,𝐻𝐸[24])/(1 −𝐹2𝑞,𝑇2,𝐻𝐸) = 1.21. While the dipolar coupling is higher for our sample, the coherence time avg(𝑇2,𝐻𝐸) is shorter and the ratio 𝑡𝑔𝑎𝑡𝑒/avg(𝑇2)  ultimately determines the fidelity limit. Our gate employs higher order decoupling, thereby extending the limiting coherence time scale to the much longer 𝑇2,𝑋𝑌8. Thus our final infidelity is a factor of (1 − 𝐹2𝑞,𝑇2,𝐻𝐸)/(1 − 𝐹2𝑞) = 1.53 better than our sample's apparent 𝑇2,𝐻𝐸  limit. By the increased coherence time, we access a new regime, where the error sources discussed in the main text become relevant.