Coherence and entanglement of inherently long-lived spin pairs in diamond

Understanding and protecting the coherence of individual quantum systems is a central challenge in quantum science and technology. Over the last decades, a rich variety of methods to extend coherence have been developed. A complementary approach is to look for naturally occurring systems that are inherently protected against decoherence. Here, we show that pairs of identical nuclear spins in solids form intrinsically long-lived quantum systems. We study three carbon-13 pairs in diamond and realize high-fidelity measurements of their quantum states using a single NV center in their vicinity. We then reveal that the spin pairs are robust to external perturbations due to a unique combination of three phenomena: a clock transition, a decoherence-free subspace, and a variant on motional narrowing. The resulting inhomogeneous dephasing time is $T_2^* = 1.9(3)$ minutes, the longest reported for individually controlled qubits. Finally, we develop complete control and realize an entangled state between two spin-pair qubits through projective parity measurements. These long-lived qubits are abundantly present in diamond and other solids, and provide new opportunities for quantum sensing, quantum information processing, and quantum networks.

Understanding and protecting the coherence of individual quantum systems is a central challenge in quantum science and technology. Over the last decades, a rich variety of methods to extend coherence have been developed [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. A complementary approach is to look for naturally occurring systems that are inherently protected against decoherence. Here, we show that pairs of identical nuclear spins in solids form intrinsically long-lived quantum systems. We study three carbon-13 pairs in diamond [14][15][16] and realize high-fidelity measurements of their quantum states using a single NV center in their vicinity. We then reveal that the spin pairs are robust to external perturbations due to a unique combination of three phenomena: a clock transition, a decoherence-free subspace, and a variant on motional narrowing. The resulting inhomogeneous dephasing time is T * 2 = 1.9(3) minutes, the longest reported for individually controlled qubits [8]. Finally, we develop complete control and realize an entangled state between two spin-pair qubits through projective parity measurements. These long-lived qubits are abundantly present in diamond and other solids, and provide new opportunities for quantum sensing [17], quantum information processing [4,18,19], and quantum networks [20].
To realize solid-state quantum systems that are inher- * T.H.Taminiau@TUDelft.nl ently long-lived, we investigate pairs of identical interacting nuclear spins. Such spin pairs are naturally and abundantly present in solids like diamond, silicon, siliconcarbide, germanium, graphene and MoS 2 [14-16, 22, 23]. Traditionally, the dynamics of such spin pairs have been regarded as a primary noise source for solid-state spin qubits [10,14,22,23,28]. In contrast, we show that spin pairs themselves provide individually controllable and decoherence-protected quantum systems.
The system that we investigate is illustrated in Fig.  1a. We consider three pairs of coupled 13 C nuclear spins in the vicinity of an NV center in a diamond at 3.7 K. The NV center provides a controllable electron spin with long coherence times that can be initialised and measured optically (Methods) [2,4,14,16,[18][19][20]. Because the NV spin creates a switchable local magnetic-field gradient over each pair, it can be used to sense and manipulate the spin pairs [14][15][16], despite their excellent protection from external influences.
A spin-1/2 pair is described by four states: |↑↑ , |↑↓ , |↓↑ and |↓↓ . We focus on the dynamics in the antiparallel subspace and define a pseudo-spin spanned by |⇑ = |↑↓ and |⇓ = |↓↑ (Methods) [14][15][16]. The pseudo-spin Hamiltonian is: in whichÎ z andÎ x are spin-1/2 operators. X is the dipolar coupling between the 13 C spins, which creates the evolution |⇑ ↔ |⇓ (flip-flops). m s = {−1, 0, +1} is the NV spin projection and Z is the difference between the two NV-13 C hyperfine couplings (Methods). Pair A and pair B are nearest-neighbour pairs oriented along the external magnetic field with X A = X B = 2π · 2080.9900(3) Hz, Z A = 2π · 130(1) Hz and Z B = 2π · 91(2) Hz (see measurements below). Pair C has a larger spatial separation between the spins resulting in X C = 2π · 188.33(2) Hz, and Z C = 2π · 2802(2) Hz. In the following, we develop initialisation, control and measurement for pairs A and B, for which X Z (see Methods for pair C control, for which Z X). Previous work has demonstrated that the pseudo-spin of pairs can be detected through decoupling sequences that toggle the m s ZÎ z term by periodically inverting the NV electron spin (Fig. 1b) [14][15][16]. For an interpulse delay of 2τ = π/ω r , with ω r = X 2 + (Z/2) 2 , the se- System and basic spectroscopy. a. We study three 13 C spin pairs (A, B and C) in a diamond. The pairs are detected and controlled using a nearby NV center. The insets show the spatial configuration of the pairs. Pair A and B are nearest-neighbour pairs oriented along the external magnetic field Bz. For pair C we show one of the three possible orientations (Supplementary Section IV). The main source of decoherence is the surrounding bath of 13 C spins (1.1% abundance). b. Sensing the pair pseudo-spins [14][15][16]. The NV electron spin is prepared in a superposition and a periodic sequence of π pulses is applied. If the interpulse delay is resonant with the dynamics of a pair, a loss of electron coherence is observed. We set τ = m2π/ωL with m an integer and ωL the 13 C Larmor frequency to avoid interactions with individual 13 C spins [14,16]. The vertical lines mark the values for τ used in this work for the three pairs (τA = τB = 120.330 μs and τC = 177.026 μs). The NV spin is prepared (RS) and read out (RO) optically (see Methods).
quence is resonant with the pseudo-spin dynamics and the effective NV-pair interaction is of the formŜ zÎz , witĥ S z the spin operator for the NV electron spin [14][15][16]. The NV center thus accumulates a phase that depends on the z-projection of the pseudo-spin. We use the NV center as a sensor to detect the spin pairs in its environment by sweeping τ (Fig. 1b) [14][15][16] and find the resonances for pair A and B (τ = 120.330 μs) and pair C (τ = 177.026 μs).
We start by developing projective single-shot measurements. Unlike all previous work, which was limited to manipulating mixed states of the parallel and antiparallel subspaces [14], these measurements enable us to initialise and measure the complete state of the spin pairs with high contrast.
Our method is based on repeated non-destructive measurements and illustrated in Fig. 2. Each repetition comprises an interaction period between the NV and the pair pseudo-spin before optical readout. During the interaction the NV electron spin accumulates a positive (negative) phase if a pair is in |⇑ (|⇓ ). For a pair in the parallel subspace (|↑↑ or |↓↓ ), the NV spin does not accumulate any phase. We choose τ such that pairs A and B interact with the NV spin simultaneously. Therefore, the NV spin accumulates a phase that depends on which of the 16 possible states the two pairs are in (Fig. 2c). By repeatedly applying this sequence, we realize a projective measurement that can distinguish multiple states in a single shot and with high contrast.
We combine these sequences to realize high-fidelity initialisation and measurement (Fig. 2e). We first apply the parity measurement sequence (20 repetitions) to herald preparation in an even parity state, and to exclude the cases for which one or both pairs are in their parallel subspace. Then, we apply a spin measurement (30 repetitions) to herald either |⇑⇑ or |⇓⇓ . Finally, we measure the pseudo-spin state with another spin measurement (30 repetitions). The resulting conditional histograms show well-isolated distributions (Fig. 2f) and an optimization of the measurement decision threshold (Fig. 2g) yields a combined initialisation and readout fidelity of 98.2(7)%. We refer to the Supplementary Information for the full optimization procedure.
We now use the developed high-contrast measurements to investigate the coherence of pair A and B. First, we perform a free-evolution experiment with the NV spin in m s = −1 (Fig. 3a), for which the NV-pair coupling is on. Because the pseudo-spin precession frequency √ X 2 + Z 2 is different for the pairs (Z A = Z B ), this measurement reveals the presence of the two pairs and characterizes their couplings Z to the NV. The two frequencies observed give Z A = 2π · 130(1) Hz and Z B = 2π · 91(2) Hz (Methods). We obtain the dephasing times from a Fourier transform The NV electron spin starts in ms = 0. TheŜzÎz interaction sequence (τ = 120.330 μs and Ns = 14) maps the state of the two pairs onto the NV spin. The NV spin is subsequently read out (RO) and reset (RS) to ms = 0. We synchronize subsequent repetitions by calibrating a waiting time τ s φ = 323.5 μs to compensate for theÎx evolution during NV readout. This ensures that the full sequence duration is a multiple of 1/X. b. Sequence to measure the pseudo-spin parity (Np = 20, Supplementary Section VII). We set τ p φ = 81 μs to synchronize subsequent measurements (sequence duration a multiple of 1/(2X)). c. XYplane of the NV Bloch sphere showing the possible phases accumulated in the spin measurement. The NV spin starts along x and picks up a positive (negative) phase for a pair in |⇑ (|⇓ ) and no phase for a pair in a parallel state (|↑↑ or |↓↓ ). Reading out along the y-axis distinguishes the 4 pseudo-spin states (blue). Note that |↓↓ (not shown) gives the same result as |↑↑ . d. XY-plane of the NV Bloch sphere under parity readout. The initial state (x-axis) and the readout axis (x-axis) are identical so that the parity of pair A and B is measured. e. Measurement sequence to calibrate single-shot readout and initialisation. The top right of each block indicates the number of repetitions. The optimal number of spin readouts is 30 (Supplementary Section VI). f. Conditional histograms for 30 spin readouts after initialisation in |⇑ |⇑ (green) and |⇓ |⇓ (blue). The initialisation conditions for the 30 preceding spin readouts are indicated in red. g. Combined initialisation and readout fidelity for |⇑ |⇑ (green) and |⇓ |⇓ (blue) for 30 spin readouts. We find an optimum of F = 98.2(7)% for a decision threshold of 14 out of 30. . These values are one to two orders of magnitude larger than for individual 13 C spins in the same sample [4].
Second, we perform the same experiment with the NV spin in m s = 0, so that the coupling to the NV is effectively turned off. Now both pairs precess with frequency X A = X B = 2π · 2080.9900(3) Hz (Fig. 3b) and a co-herent oscillation that extends past 70 s is observed. To extract the dephasing time, we measure the oscillation amplitude at various times (Fig. 3c). The resulting decay yields T * 2 = 1.9(3) minutes, a four order-of-magnitude improvement over an individual spin [4] and the longest inhomogeneous dephasing time reported for any individually controllable quantum system [8].
We now elucidate the mechanisms which lead to these remarkable coherence properties. We add a magneticfield noise term ∆Z(t) to the pseudo-spin Hamiltonian: (2) The first mechanism which enhances the coherence is the decoherence-free subspace [7] formed by the pseudo-spin states. Because the spins are identical, ∆Z(t) is given by the fluctuations of the magnetic field difference between the two spins. The atomic distance between the spins ensures near-complete immunity to noise from distant sources, such as the external magnetic field and the control signals. The main source of noise is the surrounding 13 C spin bath. Hence, ∆Z(t) can be approximated as a Gaussian distribution with a correlation time τ c [29, 30] and variance b 2 = 1 k ) is the dipolar coupling of bath spin k to the first (second) spin of the pair. We calculate the typical effective noise strength b ∼ 10 Hz by numerically simulating many spin-bath configurations. This is a noise reduction by a factor ∼ 2 due to the decoherence-free subspace (Supplementary Section II).
We first analyze the case of the NV electron spin in m s = −1 (Fig. 3a), which enables us to extract the strength of the noise due to the spin bath. Because X Z ∆Z(t), the Hamiltonian can be approximated as (Supplementary Section I) with ω −1 = √ X 2 + Z 2 . Additionally, the NV spin creates a field gradient that suppresses spin flip-flops in the bath (a frozen core [4,12]). Therefore, the noise can be treated as quasi-static and the signal decay is Gaussian [29], as experimentally observed (Fig. 3a). The dephasing time is given by [29] In this case, the coherence is enhanced by a factor ω−1 Z ≈ 20 in addition to the enhancement by the decoherencefree subspace. Finally, inserting the measured dephasing times into equation (4) yields noise strengths b A = 2π · 13.9(2) Hz and b B = 2π · 12.5(4) Hz. These values are consistent with the inter-pair distance and 13 C concentration (Supplementary Section II).
Second, we analyze the case with the NV electron spin in m s = 0 (Fig. 3b) The eigenenergies are now first-order insensitive to ∆Z(t) as the spin pair forms a clock transition due to the coupling X, a second mechanism that enhances coherence. Note that the clock transition in this system does not require a specific magnetic field value, as the simultaneous decoherence-free subspace removes the dependence on global magnetic fields.
The decoherence-free subspace and clock transition alone cannot yet explain the observed m s = 0 dephasing time. In particular, for quasi-static or slow noise the coherence would be limited to ∼ 10 s (Supplementary Section I). However, the increased coherence, in combination with the lack of a frozen core for m s = 0, unlocks a new regime in which the nuclear-spin bath fluctuations become relatively fast (τ c X/b 2 ). A mathematically equivalent Hamiltonian was analysed theoretically by Dobrovitski et al. [30]. The resulting time constant is The dependence on the correlation time τ c reveals a third mechanism, similar to motional narrowing [30], that further enhances the coherence. Inserting the parameters obtained from the m s = −1 measurement and a typical value for τ c ∼ 0.1 s [18], inhomogeneous dephasing times of ∼ 100 s are predicted. Together these three mecha-nisms thus provide an explanation for the observed dephasing times.
To further analyse the different physical regimes that play a role, we investigate pair C (Fig. 1). Because Z X, the dynamics are different and the clock transition can be switched on (m s = 0) and off (m s = −1) (Supplementary Section I). We develop complete control, initialisation and single-shot readout of such pairs in the Methods.
For evolution under m s = 0, the situation is similar to pairs A and B. We find T * 2 = 0.6(1) s, which is reduced compared to pairs A and B because the smaller coupling X makes the clock transition less effective (Extended Data Fig. 3). Additionally, similar values obtained for spin echo (T 2 = 0.7(1) s) and relaxation measurements (T 1 = 0.9(2) s) indicate that relaxation plays a role in limiting the coherence (Supplementary Section I). For m s = −1, a frozen core is formed and the clock transition is turned off, so that the noise ∆Z(t) affects the eigenfrequencies linearly. We find T * 2 = 18(1) ms with Gaussian decay, indicating quasi-static noise [29], which is consistent with spin echo (T 2 = 0.3(2) s T * 2 ) and relaxation measurements (T 1 1 s) (Extended Data Fig.  3). In this case, there is no significant coherence protection and the results are similar to individual 13 C spins [4]. Finally, we demonstrate the potential of the spin pairs as qubits by demonstrating an entangled state of pair A and pair B. We create entanglement through subsequent projective measurements of the σ y σ y and σ z σ z pseudo-spin parity (Fig. 4a). We herald on outcomes σ y σ y = +1 and σ z σ z = −1, so that the resulting state is 1 √ 2 (|⇑⇓ + |⇓⇑ ). This state is a 4-spin entangled state 1 √ 2 (|↑↓↓↑ + |↓↑↑↓ ) that encodes two qubits of information in two long-lived pseudo-spin states.
To characterize the resulting state we first measure parity oscillations by varying the evolution time t (Fig.  4a). The observed frequency is 4.20(4) kHz, which equals 2X/2π, as expected (Fig. 4b). To determine the state fidelity, we measure the pseudo-spin parity operators σ x σ x , σ y σ y and σ z σ z . We realize the required single-qubit rotations through waiting times (for x-rotations) and dynamical decoupling sequences with the NV spin in an eigenstate (for z-rotations) (Fig. 4a). Figure 4c shows the resulting expectation values, which yield a fidelity F = 0.74(2), confirming entanglement (F > 0.5). This result highlights the high-fidelity initialisation, control, and non-destructive measurements realized.
In conclusion, we have developed complete control over multiple nuclear-spin pairs.
These spin pairs provide inherently long-lived quantum states due to a unique combination of three physical phenomena: a decoherence-free subspace, a clock transition and a variant of motional narrowing. This inherent coherence protection makes spin pairs promising qubits for a variety of applications, including quantum sensing [17] and quantum information processing [4]. In particular, due to the small effective coupling to the NV spin, they might provide robust quantum memories for optically connected quantum networks [20,21]. Such spin pairs are available for most NV centers (Supplementary Section IV) and are present in a variety of other materials. Therefore, our results reveal a new, promising, and abundantly available resource for quantum science and technology.

METHODS
Sample and NV center. The experiments are performed on a naturally occurring NV center in a cryogenic confocal microscope (3.7 K). The diamond was homoepitaxially grown using chemical vapor deposition and cleaved along the 111 axis (Element Six). There is a natural abundance (1.1%) of 13 C. The NV centre was selected on the absence of strongly coupled 13 C spins exceeding ∼ 500 kHz hyperfine coupling, but without any other criteria on the spin environment or spin pairs. The NV electron spin has a dephasing time of T * 2 = 4.9(2) μs and a spin echo time of T 2 = 1.182(5) ms. The electron relaxation (T 1 > 1 h) at this temperature is negligible [14]. We measure the NV spin state in a single shot using spin-selective optical readout [18]. The readout fidelities are 0.905(2) (0.986(2)) for the m s = 0 (m s = −1) state with an average fidelity of F = 0.946 (1). The dynamical decoupling sequences follow the XY8-scheme to mitigate pulse errors [31]. Pseudo-spin Hamiltonian. The Hamiltonian for two 13 C spins in the vicinity of an NV center in an appropriate rotating frame and under the secular approximation can be written as: where ω L = γ c B is the 13 C spin Larmor frequency, with γ c the 13 C gyromagnetic ratio. B is the magnetic field along the NV-axis. I (i) are the spin-1 2 operators acting on spin i, m s = {−1, 0, +1} are the NV electron spin states and A (i) = [A x , A y , A z ] is the NV-13 C hyperfine interaction vector of spin i. H D is the dipolar interaction between two 13 C spins. For a large magnetic field compared to the dipolar (X) and hyperfine couplings (A (1) , A (2) ) H D can be written as: where µ 0 is the vacuum permeability, r 12 is the vector between the two 13 C atoms and θ 12 the angle between the magnetic field axis and the pair axis. Since ω L = γ c B = 2π · 432.140 kHz is large compared to the dipolar (X) and hyperfine couplings (A (1) , A (2) ), the antiparallel states |↑↓ and |↓↑ form an isolated subspace in which we define a pseudo-spin 1 2 as |⇑ = |↑↓ and |⇓ = |↓↑ Z originates from the difference of the hyperfine couplings of the two spins to the NV electron spin, and is given by [28] where We obtain T = 0.53(4) s, n = 2.1(4), f A = 9.07(6) Hz and f B = 7.0(1) Hz (measured with a 10 Hz detuning with respect to 2086 Hz). Using f = √ X 2 + Z 2 and X = 2π · 2080.9900(3) Hz, the values for f A and f B yield Z A = 2π·130(1) Hz and Z B = 2π·91(2) Hz. The observed shape of the decay (n = 2.1(4)) is in agreement with the predicted Gaussian (n = 2) decay for quasi-static noise (Supplementary Section I).
Note that the precise value obtained for X deviates from simple theoretical estimates and is analyzed in Supplementary Section X. The Fourier transform is fitted to We obtain f 0 = 0.2402(3) Hz and σ 0 = 0.0074(3) Hz. The data in Fig. 3c is fitted to exp(−(t/T ) n ) obtaining T = 1.9(3) min and n = 0.23 (2). Note that n deviates from the simple exponential decay (n = 1) associated to equation 6, indicating that other effects beyond pure dephasing contribute to the decoherence (Supplementary Section I).
Pseudo-spin control methods. The control methods for pairs A and B are given in the main text and summarized in Extended Data Fig. 1.
Pair C has a dipolar coupling X = 2π · 188.33(2) Hz and a hyperfine difference Z = 2π·2802(2) Hz. Therefore, Z X, in contrast to pairs A and B for which X Z. This changes the dynamics in two ways. First, for m s = −1, Z is the dominant term in the pair frequency ω −1 = √ X 2 + Z 2 and thus sets the location of the resonance in Fig. 1b. Second, the effective NV-pair interaction during the dynamical decoupling sequence becomesŜ zÎx [14] (Extended Data Fig. 1).
For spin pairs with Z X, the timing of repetitions is complicated by the fact that the m s = 0 and m s = −1 evolution frequencies and eigenstates differ significantly. Here, we mitigate this by minimizing the NV readout time (RO, ∼ 5 μs) and applying a fast reset of the NV spin (RS), so that the potential time spent in m s = −1 is small. Because the states that the measurement projects onto ( 1 √ 2 (|⇑ ± |⇓ )) are eigenstates of the m s = 0 evolution, there is no timing requirement after resetting the NV and we simply concatenate subsequent measurements.
For pairs A and B (X Z), we use free evolution and dynamical decoupling sequences to realize universal single-qubit control for the pseudo-spins ( Fig. 4; Extended Data Fig. 1). For pair C (Z X) all single-qubit operations can be obtained by letting the system evolve freely. Evolution with the NV electron spin in m s = 0 implements a rotation around the x-axis, and evolution under m s = −1 realizes a rotation around the z-axis (Extended Data Fig. 1). Note that the z-axis rotation is approximate as Z is finite. In principle, this can be corrected for but this is not done here. We use the pair C control to measure the pseudo-spin dephasing time T * 2 , the spin echo time T 2 and the relaxation time T 1 in both NV electron spin states, see Extended Data Fig. 3. Spectroscopy and control of the full Hilbert space. Most of the work presented is focused on initialising, controlling and measuring the states in the antiparallel subspace of spin pairs, i.e. |↑↓ and |↓↑ . In Extended Data Fig. 4, we demonstrate that the entire Hilbert space of the spin pairs can be controlled by RF driving the singlespin-flip transitions of pair C.
The single-spin transition frequencies are ω 1 = 2π · 429.314(5) kHz and ω 2 = 2π · 432.122(7) kHz (Extended Data Fig. 4a). Since the frequency of a single-spin (1) = 2π · 2826(5) Hz and A (2) = 2π · 18 (7) Hz. Note that these values assume that A ⊥ is of similar magnitude, so that it can be neglected. The frequencies observed are consistent with the characteristic 13 C frequencies (ω L = 2π · 432.140 kHz), further corroborating our assignment of 13 C-13 C pairs as the source of the signals. These results also demonstrate selective initialisation, control and measurement of an individual carbon spin with negligible coupling to the NV by using its coupling to neighbouring spins. Spin 2 couples negligibly to the NV (18(7) Hz), so that it overlaps in precession frequency with most of the spin bath. Nevertheless, it can be initialised and controlled selectively by using the NV to directly detect its flip-flops with spin 1 (i.e. pseudo-spin dynamics).   A, B). The top row indicates the sequence performed on the NV electron spin, the middle row the corresponding pair dynamics in the XZ-plane of the pseudo-spin Bloch sphere and the bottom row indicates the effective pseudo-spin Hamiltonian term under that sequence. For the left two columns the rotation frequencies are given inside the Bloch spheres. From left to right the sequences are free evolution in ms = 0 and ms = −1, and a dynamical decoupling sequence with τ resonant, i.e. 2τ = π/ X 2 + (Z/2) 2 . Rotations that are unconditional on the NV electron spin state can be obtained by setting τ = π/ X 2 + (Z/2) 2 (unconditional z-rotation) and by setting τ far off-resonant (unconditional x-rotation), but these are not shown or used here. Note that the z-rotation frequency depends on the hyperfine field difference Z, so that pair A and B can, in principle, be controlled individually. b. Pseudo-spin dynamics of pairs with Z X (pair C). Like above, z-and x-rotations that are unconditional on the NV electron-spin state can be obtained by setting different values for τ (not shown). Together these operations enable universal control of the system consisting of the three pseudo-spins and the NV center.   Fig. 2. Projective spin and parity measurements for Pair C. a. Sequence to measure the spin state of pair C. The NV starts in ms = 0. Because Z X, the effective interaction between the NV spin and pseudo-spin isŜzÎx. During the NV readout (RO) and reset (RS), the NV can spend an unknown time in ms = −1 which causes dephasing of the pair spin. Additionally pair C undergoes a deterministic z-rotation for any known time spent in ms = −1. To minimize these effects, we use a fast readout and reset. Ns = 8. b. Sequence to measure the parity of the two spins that make up pair C. Note that in this case the timing of the sequence is unimportant, as evolution in ms = −1 or ms = 0 does not change the parity. Np = 14. c. XY-plane of the NV Bloch sphere during the NV-pair interaction in a. The NV picks up a positive or negative phase depending on the x-projection of the pair pseudo-spin and no phase when the pair is in the parallel subspace. d. XY-plane of the NV Bloch sphere during the NV-pair interaction in b.   Sequence to reveal the transitions between the subspaces. First, the pair is initialised in the antiparallel subspace through a parity measurement, then an RF pulse with variable frequency with the NV in ms = −1 is applied, and finally the subspace population is measured using another parity measurement. If the frequency of the RF pulse is resonant with a single-spin transition, the spin pair changes its subspace. c. Measurement result. Four transitions are observed corresponding to the marked transitions in a. The green dashed line corresponds to the bare Larmor frequency ωL = 2π · 432.140 kHz. We fit the data to four Lorentzians and extract ω1 = 2π · 429.314(5) kHz, ω2 = 2π · 432.122(7) kHz. For the left (right) dips we also obtain X = 2π · 184(3) (2π · 194(4)) Hz. These results corroborate the assignment of the signals to 13 C pairs and enable complete control over the full pair state. In this section we provide the derivations for the spin-pair coherence alongside a more detailed discussion. The key idea is that different physical regimes are accessed depending on the NV electron-spin state and the spin-pair parameters. In particular, the analysis shows that the long coherence times observed for the m s = 0 state are made possible by a unique combination of a decoherence-free subspace, a clock transition and a variant of motional narrowing.

Supplementary information for "Coherence and entanglement of inherently long-lived spin pairs in diamond"
We first focus on interactions and noise that cause dephasing or relaxation within the pseudo-spin subspace. Potential relaxation out of the pseudo-spin subspace, for example due to flip-flops with the surrounding bath spins, is discussed separately below. The Hamiltonian for a pair in the pseudo-spin subspace (spanned by |⇑ = |↑↓ and |⇓ = |↓↑ ), including a magnetic-field noise term ∆Z(t), is where X is the dipolar coupling between the two 13 C spins, m s = {0, −1} the electron state (m s = +1 is not populated during the experiments) and Z is due to the hyperfine field gradient induced by the electron spin. Because the two spins are identical, the pseudo-spin forms a decoherence-free subspace. It is therefore only sensitive to field gradients. Noise from distant sources, such as fluctuations of the applied magnetic field (distance to the closest magnet: ∼ 1 cm) and noise from the control electronics on the on-chip stripline (distance: ∼ 10 µm) has negligible influence on the pseudo-spin dynamics, see section III. Therefore, the surrounding 13 C bath is the main source of noise and the origin of ∆Z(t).
The interaction between two central spins and a bath can be described as where A Only the latter term I whereÎ z is the pseudo-spin operator and the sum is over all k bath spins. The pair forms a decoherence-free subspace: only the difference in coupling of a bath spin k causes noise on the pseudo-spin of the pair. For a single spin A (2) k = 0 andÎ z becomes the single-spin operator I z .
Following analyses by Anderson et al. [1] and Klauder et al. [2], we now take a classical approach to the noise in the pseudo-spin subspace. We model ∆Z(t) as an Ornstein-Uhlenbeck process with correlation function ∆Z(t)∆Z(0) = b 2 exp(−t/τ c ) where τ c is the correlation time of the bath due to the intra-bath dynamics H B and b 2 = 1 k ) 2 is the variance of the noise. We calculate the distributions for b numerically (see section II). For nearest-neighbour pairs (like pairs A and B), we find a typical b ∼ 10 Hz. For the parameters of pair C we find a larger value b ∼ 15 Hz, consistent with the larger distance between the spins leading to less correlated noise and a less effective decoherence-free subspace. As a reference, for an individual 13 C spin one has b ∼ 20 Hz. For the correlation time when the NV spin is in m s = 0, a typical value τ c = 0.1 s can be estimated from previous experiments [3].
A. Pair C and ms = −1: no clock transition, a frozen core, quasi-static noise For pair C we have Z X. For the NV in m s = −1, the effect of the coupling X becomes negligible, so that equation (S1) can be approximated as In this case, ∆Z(t) directly and linearly modifies the eigenenergies. There is no coherence protection related to a clock transition. Additionally, in the m s = −1 state, the NV magnetic field gradient creates a frozen core, in which nuclear spin flip-flops are suppressed [4,5]. As a result, the dephasing time (T * 2 ) is short compared to the correlation time of the noise τ c , and ∆Z(t) can be described as quasi-static.
This regime leads to a Gaussian decay of exp −b 2 τ 2 /2 with T * 2 = √ 2 b [1]. Therefore, we can extract b from the experimental data. The measured T * 2 time of 18(1) ms yields b C = 2π · 12.5(7) Hz, consistent with the expected distribution for the inter-pair distance of pair C (Supplementary Fig. 2). The hypothesis that the noise can be treated as quasi-static is further corroborated by the fact that a large increase in coherence is observed for a spin echo (T 2 = 0.3(2) s, Extended Data Fig. 3d).
In summary, for pair C and m s = −1, we probe a regime where the clock transition is turned off, the decoherencefree subspace has a reduced influence (larger inter-pair distance) and the bath noise can be treated as quasi-static (frozen core and T * 2 τ c ). This regime and the resulting T * 2 is similar as for an individual nuclear spin in the same environment [4]. No significant enhancement of coherence is obtained for the spin pair.
B. Pair A, B and ms = −1: a detuned clock transition, a frozen core, quasi-static noise For pair A, B we have X Z. Additionally we typically have a situation in which X Z ∆Z(t). Taking only terms that contribute to dephasing, we can therefore approximate the Hamiltonian with the NV in m s = −1 using a Taylor series expansion as where Similarly to the case of pair C, we expect a Gaussian decay shape, but now with . While the coupling X now creates a clock transition, the system is detuned from the ideal point, because Z ∆Z(t). As a result, the effect of the noise is reduced by a factor ω−1 Z ≈ 20, consistent with the increase of T * 2 of pair A and B compared to pair C.
The experimental data contains the decay of both pair A and B that are generally not equal. We extract two decay times from the Gaussian fit of the Fourier transform (Fig. 3a), obtaining T * 2,A = 0.26(2) s and T * 2,B = 0.39(6) s. For pair A with Z A = 2π · 130(1) this corresponds to b A = 2π · 13.9(2) Hz and for pair B with Z B = 2π · 91(2) Hz this corresponds to b B = 2π · 12.5(4) Hz. The values agree with typical values for b (Supplementary Fig. 2).
In conclusion, in this case we probe a regime where the pair interaction X creates a clock transition (as X Z), but the system is detuned because Z ∆Z(t). Additionally, the NV electron spin creates a frozen core and the noise can thus be treated as quasi-static. The resulting dephasing time is enhanced by a factor √ X 2 + Z 2 /Z.

C. Pair A, B and ms = 0: clock transition, motional narrowing
In m s = 0, the noise cannot be treated as quasi-static anymore, because flip-flops between 13 C spins occur more frequently. We therefore have to take into account the correlation time τ c of the bath as well as its strength b. The Hamiltonian for a pair in the pseudo-spin subspace in m s = 0 is For pair A, B it holds that X ∆Z(t) for typical values of ∆Z(t) (Supplementary Fig. 2). We initially assume that the bath has no significant frequency components leading to direct transitions between the pair eigenstates (X 1/τ c ), but will come back to this effect at the end of the section. Then, following the analysis by Dobrovitski et al. [6], we can approximate the Hamiltonian as The system now forms an effective clock transition and the noise term ∆Z(t) enters quadratically in the Hamiltonian.
This model can be solved analytically for the expectation value S z of the pair pseudo-spin [6]: where P = R 2 − 2ib 2 R/X and R = 1/τ c where τ c is the correlation time of the bath.
Equation (S8) holds generally, but it is possible to consider three different regimes separately. These regimes are defined by the rate of the noise bath fluctuations R compared to the effective coupling to the noise bath of b 2 /2X. First, we consider a quasi-static bath. This leads to a non-exponential decay of the form [6,7] Second, the noise bath can be 'slow' compared to the coupling strength to the noise: R b 2 /X. For short times, during which the bath is static, the decay follows equation (S9). For longer times, during which slow bath dynamics have to be taken into account, the decay is of the form [6] M (t) = exp −bt R/4X . (S10) Third, the bath dynamics can be fast compared to the magnitude of the noise: R b 2 /X. In this regime we can approximate the solution in equation (S8) as [6] The expected decay time of the Ramsey experiment is T * 2 = 4X 2 R b 4 . In this regime the pair coherence benefits from an effect that is similar to motional narrowing [6]: T * 2 linearly increases with the rate of fluctuations R. The result differs from the standard case of motional narrowing in that there is an additional frequency shift of b 2 /2X [6].
We now discuss which of these regimes governs the coherence of pair A and B. As shown in the previous section b A = 2π · 13.9(2) Hz and b B = 2π · 12.5(4) Hz. From the main text we know X = 2π · 2080.9900(3) Hz (Fig. 3b). First we consider the expected dephasing times for slow baths. Equation (S8) is plotted in Supplementary Fig. 1 for pair A (b A = 2π · 13.9(2) Hz) with a correlation time of τ c = 10 s. From Supplementary Fig. 1 it is clear that a slow bath cannot explain the long inhomogeneous dephasing time observed in the main text (T * 2 = 1.9(3) min, Fig. 3c). Because typical measured relaxation times T 1 for individual 13 C spins indicate τ c = 0.1 − 0.3 seconds [3], one would indeed expect that the bath noise enters a fast regime: the fast bath condition R b 2 /X is satisfied for bath correlation times of τ c 0.5 s.
The envelopes for a fast bath are shown in Supplementary Fig. 1. These results show that significantly longer dephasing times are expected in this regime compared to the slow bath regime. We conclude that the measured T * 2 = 1.9(3) min can only be explained in this fast bath regime where an effect similar to motional narrowing further enhances the inhomogeneous dephasing time.
Since the coupling to the noise of pair A (b A ) and B (b B ) were determined above, the bath fluctuation rate can be estimated under the assumption that dephasing is limiting. Using T * 2 = 1.9(3) min and X = 2π · 2080.9900(3) Hz we obtain R A = 10(2) Hz and R B = 6(2) Hz. These values are consistent with previously measured values for T 1 of individual 13 C spins with the NV electron spin in m s = 0 [3].
Finally, we consider relaxation as a potential mechanism limiting the dephasing time. Spectral components around frequency X could lead to direct transitions between the antiparallel pair eigenstates . In the regime of a fast bath this can be treated analytically and introduces a factor that multiplies equation (S8) by exp −b 2 Rt/2X 2 [6]. If this type of relaxation is dominant we would obtain T * 2 = 2X 2 b 2 R . The identified values for b A , b B , X and T * 2 give R A = 400(75) Hz and R B = 500(100) Hz, which is inconsistent with previously measured single 13 C spin T 1 values [3] and with the typical 13 C-13 C couplings for this 13 C concentration. We conclude that such relaxation within the pseudo-spin subspace is unlikely to contribute to the observed coherence curves.
In summary, the long observed dephasing times for pair A, B are the result of three different physical effects working together. First, the effective noise b is reduced because correlated noise does not affect the spin pair, i.e. the pair forms a decoherence-free subspace (DFS). Second, since X ∆Z(t), the pair pseudo-spin forms a clock transition (equation (S7)). It is therefore only second-order sensitive to noise following ∆Z 2 (t)/2X. Third, the DFS and clock transition alone are not sufficient to explain a dephasing time of 1.9(3) min. Only in the regime of a fast bath, in which the pair benefits from an effect similar to motional narrowing, can such dephasing times be realized.

D. Pair A, B and ms = 0: other decoherence mechanisms
Above, we show how a combination of three effects strongly suppresses dephasing for the spin pairs. This strong suppression of dephasing is a necessary condition to obtain the long T * 2 times observed. However, it does not imply that the obtained T * 2 times and decay curves are explained by and purely limited by dephasing. Indeed, the decay envelope observed for pairs A and B (Fig. 3c) deviates from the simple exponential decay obtained in equation (S11). In particular, the fit yields a decay curve following e −(t/T * 2 ) n with n = 0.23(2) (Methods) and the data suggests additional features in the decay shape that are not captured by the fit curve. These observations indicate that other mechanisms, like coherent interactions with the bath or relaxation of the pair spins due to flip-flops with the bath spins are contributing to decoherence. Such effects strongly depend on the microscopic environment of the pairs and are challenging to treat generally. Future research could aim to understand the microscopic environment of the pairs and determine the mechanisms limiting the observed coherence times.

E. Pair C and ms = 0
The analysis for pair C is analogous to the above analysis. However, there is an important difference between pair A, B and pair C. Namely, the dipolar coupling X is an order of magnitude smaller for pair C. Since the effective noise strength for pair C (b C = 2π · 12.5(7) Hz) is similar to the noise for pair A (b A = 2π · 13.9(2) Hz) and pair B (b B = 2π · 12.5(4)), the approximation of a fast bath cannot be easily made. Therefore we are in an intermediate regime where the full expression in equation (S8) describes the dephasing. The pair C coherence in m s = 0 therefore still benefits from the decoherence-free subspace and clock transition, but to a lesser degree from motional narrowing. For pair C, spin echo (T 2 ) and relaxation (T 1 ) measurements in both electron states were taken (Extended Data Fig. 3). When the NV electron spin is in m s = −1, the NV spin creates a hyperfine field gradient that slows down spin flip-flops in the bath (a frozen core [4,5]). The noise a 13 C pair experiences is therefore expected to be quasistatic. Given quasi-static noise, a spin echo is expected to increase coherence. This is in agreement with the marked increase in coherence time: T 2 = 0.3(2) s (Extended Data Fig. 3d). Similarly, the frozen core also suppresses flip-flops involving one of the 13 C spins of the pair, leading to long relaxation times. For pair C we do indeed find a relaxation time T 1 1 s, comparable to that of single 13 C spins [4]. When the NV electron spin is in m s = 0, 13 C spins are no longer detuned and flip-flops can become limiting for spin coherence. For pair C we find a spin echo time of T 2 = 0.7(1) s. The relaxation time is T 1 = 0.9(2) s or T 1 = 3.6(7) s depending on the initial state. This suggests that the relaxation mechanism is dependent on initialisation in the singlet ((|⇑ − |⇓ )/ √ 2) or triplet state ((|⇑ + |⇓ )/ √ 2). Furthermore these results indicate that the coherence of pair C may be limited by relaxation.

II. DECOHERENCE-FREE SUBSPACE: DISTRIBUTIONS FOR b
The noise ∆Z(t) on a spin pair originates from the surrounding 13 C spins. As a pair is only sensitive to field gradients (a decoherence-free subspace), distant external noise sources can generally be neglected. There are k bath spins that each create a field difference A k on the pair (Supplementary Fig. 2a). As outlined above, we model ∆Z(t) as an Ornstein-Uhlenbeck process with a variance b 2 = 1 The question that we address in this section is what b is for typical spin baths. We numerically generate 10 5 different baths (1.1% 13 C abundance) surrounding a pair (Supplementary Fig. 2b,c) or single spin (Supplementary Fig. 2d) in a volume of 15 × 15 × 15 unit cells. For each bath we calculate b 2 = 1 4 > 50 Hz, i.e. we exclude strongly coupled spins for which the system would not be a well-defined spin pair anymore. The expectation is that the closer the spins of the pair are, the more correlated the noise and the smaller b is.
The result for a nearest neighbour pair oriented along the magnetic field (like pair A or B) is shown in Supplementary  Fig. 2b. We find a mean of 10 Hz and a standard deviation of 4 Hz. For the parameters of pair C ( Supplementary  Fig. 2c) we find a mean of 14 Hz and a standard deviation of 5 Hz. Lastly for an individual 13 C spin ( Supplementary  Fig. 2d) we find a mean of 20 Hz and a standard deviation of 6 Hz. A decrease in the effective noise is observed for the pairs compared to an individual spin. Furthermore, the closer the pair spins are, the smaller the effective noise is.

III. DECOHERENCE-FREE SUBSPACE AND EXTERNAL NOISE
The decoherence-free subspace makes pairs nearly immune to noise from distant sources. In this section we consider the effect of two such external sources. For single 13 C spins in the same sample typical inhomogeneous dephasing times of 10 ms are observed [4], which sets a bound on the noise strength. If we take the extreme case that all noise comes from an external source, i.e. not from the spin bath, this gives an upper limit of the noise magnitude b = √ 2/T * 2 = 2π · 22.5 Hz. Since we are interested in an order of magnitude estimate, we take b ∼ 2π · 10 2 Hz, corresponding to magnetic field fluctuations of δB ∼ 10 −5 T. Now we consider that these fluctuations originate from either the on-chip MW line that we use to apply microwaves or from the external magnets that we use to apply a magnetic field.

A. Microwave line
We approximate the microwave line as an infinite wire that generates a field at a distance r from the wire of magnitude B(r) = µ 0 I/2πr where µ 0 is the vacuum permeability and I the current through the wire. Given r ∼ 10 μm and δB ∼ 10 −5 T, we obtain I = 2πrδB/µ 0 ∼ 10 −4 A. Now we turn to the effect of this noise on a decoherencefree subspace formed by a 13 C pair. The positions of the pair spins are r a = 10 μm and r b = 10 μm + a 0 where for a 0 we take a conservative value of ∼ 10 −9 m. Given I = 10 −4 A, the MW line would add a field difference to the decoherence-free subspace of ∆B = B(r a ) − B(r b ) = µ0I 2π ( 1 ra − 1 r b ) ∼ 10 −10 T. That corresponds to ∼ 10 −3 Hz which has a negligible effect on the coherence.

B. Magnet
The external magnetic field comes from a cylindrical permanent magnet. We calculate the effect of that field (∼ 0.04 T) on the decoherence-free subspace of a pair. From the above we know that the maximum field fluctuations are on the order of δB ∼ 10 −5 T. We consider a magnet with an NV-magnet distance r ∼ 10 −2 m, the radius of the magnet R = 5 mm, the length L = 5 mm and the remanence field B r = 1.5 T. To calculate the magnetic field at r we use At r ∼ 10 −2 m the magnetic field is ∼ 0.04 T. The expected effect of the permanent field on the decoherence-free subspace is then B(r a ) − B(r b ) ∼ 10 −8 T or 10 −1 Hz. We have used r a = 1 cm and r b = 1 cm + a 0 where for a 0 we take a conservative value of 10 −9 m. This is a constant field difference (c.f. Z) added to the pair. However, the field fluctuations are δB ∼ 10 −5 T, more than three orders of magnitude less than the permanent field of ∼ 0.04 T. The influence of δB on the decoherence-free subspace is therefore < 1 mHz which has a negligible effect on the coherence.

IV. EXPECTED NUMBER OF NEAREST-NEIGHBOUR PAIRS PER NV
In this section we address how many nearest neighbour pairs with similar Z as pair A and B one would expect surrounding a typical NV center. To that end we generate 10 4 different 13 C baths with 1.1% abundance surrounding an NV center in a volume of 15 × 15 × 15 diamond unit cells. For every generated bath, we look for the nearest neighbour pairs along the magnetic field axis and calculate the hyperfine field difference Z due to the NV, assuming a dipolar NV-13 C interaction. Then we estimate a controllable region of 2π · 50 < Z < 2π · 500 Hz. The upper bound comes from the required condition X Z and the (approximate) lower bound is a limit due to the detrimental effect of electron dephasing for a large number of dynamical decoupling pulses. Additionally for smaller Z resolving a pair from the background bath of pairs with small Z is expected to be challenging. For every generated bath, we determine how many pairs satisfy this condition and plot the result in Supplementary Fig. 3. The expected number of such nearest neighbour pairs per NV is 1 ± 1. Moreover, more than 70% of simulated NVs host at least one nearest-neighbour pair, indicating that such pairs can commonly be found.
In the above we only consider nearest-neighbour pairs along the magnetic field axis. Other pairs with smaller X and larger Z can also be detected and controlled (pair C, see Methods). In Supplementary Table I we show the ten largest values of X with their corresponding occurrence and vector r between the two 13 C spins. The expected number of pairs with 2π · 50 < Z < 2π · 500 Hz for a volume of 15 × 15 × 15 unit cells surrounding an NV center. We estimate that pairs with a Z within the indicated region to be controllable with high fidelity. When Z is larger, it becomes comparable to X and the control methods presented in this paper do not hold anymore. When Z is much smaller, more pulses are required in the dynamical decoupling sequence resulting in more electron dephasing, and resolving the pair from the background bath of pairs becomes more challenging. Supplementary Table I. The occurrence of a given coupling X when fixing one of the 13 C spins of the pair and moving the other around the lattice. The vector between the two 13 C spins of the pair is given for each X coupling in units of a0/4 where a0 is the lattice constant of diamond. All permutations of the entries of r give the same coupling X.

V. PAIR INITIALISATION AND READOUT CALIBRATION
In this section we describe the optimization of the parameters used for initialisation and single-shot measurements. The most important trade-off lies in the number of repetitions of the measurement sequences. On the one hand, increasing the number of repetitions improves the fidelity because the different states can be distinguished better and the effect of the NV electron spin dephasing is diminished. On the other hand, the pair spin state decoheres during the measurement, limiting the maximum number of repetitions. Therefore, there is an optimum in the number of repetitions and the corresponding decision thresholds used.
We first describe our approach to optimize the parameters in general. We will call the two states that we want to optimally distinguish |a and |b . In an initialisation step, we generally use k repetitions and we record the number of counts (m s = 0 outcomes) as N (k). This initialisation step is then defined by N (k) > N a and N (k) < N b where N a and N b are the thresholds set (red lines in Fig. 2f,g). In case there is a two-step initialisation process (see e.g. Fig.  2e) we denote the number of counts in the first step as N (p) with condition N (p) > N 0 where N 0 is the threshold set. In the readout step, we use m repetitions and obtain N (m) counts. Two histograms are obtained (see e.g. Fig.  2f), one corresponding to each initialised state. To optimally distinguish these states, we sweep a threshold T (see e.g. Fig. 2g) and obtain the combined initialisation and readout fidelity as We then optimize the fidelity for the number of repetitions m, the measurement decision threshold T , and the initialisation thresholds N a , N b and if applicable N 0 . For the initialisation, the number of repetitions k is not as important. The main trade-off now lies in the initialisation thresholds: we pick values that balance the resulting fidelity and the success probability (experimental rate). Namely, the stricter the threshold is, the higher the fidelity but the lower the experimental rate.

VI. CALIBRATION OF THE SPIN MEASUREMENT FOR PAIR A AND B
For pair A and B we calibrate the spin measurement to optimally distinguish |a = |⇑ |⇑ and |b = |⇓ |⇓ . The sequence is shown in Supplementary Fig. 4a. First, we initialise the parity of pair A and B with p = 20 repetitions and a threshold of N 0 = 12. To initialise the two states of interest we then use k = 30 spin readouts and set thresholds of N a = 25 and N b = 3 counts. For a varying number of readouts m we now determine the optimal threshold and corresponding fidelity as in Fig. 2g in the main text. In Supplementary Fig. 4b, we plot the average fidelity and corresponding optimal threshold for a varying number of readouts. The optimum is at m = 35 readouts with a threshold of T = 16. We obtain a combined initialisation and readout fidelity of F = 98.4(7)%. For the results in the main text we use m = 30 and T = 14 which gives, within error, the same fidelity (F = 98.2(7)%). The sequence used to calibrate the spin readout. First we select on > 14/20 counts in 20 parity readouts to initialise the pairs in 1 2 |⇑ |⇑ ⇑| ⇑| + 1 2 |⇓ |⇓ ⇓| ⇓|. Then we either initialise the pairs in |⇑ |⇑ (> 25/30) or |⇓ |⇓ (< 3/30) with 30 spin readouts. Finally the spin is read out using m spin readouts. b. We vary the number of readouts m and find the optimal threshold for each m. The average fidelity (equation (S13)) is plotted against the number of readouts m and the corresponding optimal threshold is indicated for that number of readouts. The initial increase of fidelity with the number of repetitions originates from the limiting effect of dephasing of the NV electron spin during the interaction sequence. For a small number of repetitions, this dephasing limits the measurement fidelity. But with more repetitions the conditional histograms (Fig. 2f) separate and NV electron dephasing plays a limited role. For a large number of repetitions a slow decay in the fidelity is observed, which is due to decoherence of the pair spin during the measurement sequence.

VII. CALIBRATION OF THE PARITY MEASUREMENT FOR PAIR A AND B
For pair A and B the parity measurement is optimized to distinguish the pair parity. To calibrate the measurement, we use the states |a = 1 2 |⇑ |⇑ ⇑| ⇑| + 1 2 |⇓ |⇓ ⇓| ⇓| and |b = 1 2 |⇑ |⇓ ⇑| ⇓| + 1 2 |⇓ |⇑ ⇓| ⇑|. To prepare these two states, we use a two-step initialisation (Supplementary Fig. 5a). First, we use p = 20 parity initialisations with a threshold N 0 = 14 to obtain |a . Then we apply a π/2 pulse around x to obtain a fully mixed state in the antiparallel subspace ({|⇑⇑ , |⇑⇓ , |⇓⇑ , |⇓⇓ }). This first step of the initialisation has reduced the Bloch sphere in Fig. 2d to just contain antiparallel states (the pseudo-spin states), improving subsequent initialisation steps. The second step of the initialisation uses k = 20 readouts with initialisation thresholds of N a = 15 and N b = 2. We set the number of decoupling elements for the parity sequence to N p = 20 in Fig. 2b. This is less than twice the decoupling elements for the spin sequence (N s = 14 in Fig. 2a). The optimal interaction time is smaller than a π-rotation due to the loss of NV electron spin coherence with increasing N p .
In Supplementary Fig. 5b the histograms for |a (green) and |b (blue) are shown for m = 18 readouts. We use equation (S13) to obtain the fidelity of the individual states ( Supplementary Fig. 5c) and the average fidelity while sweeping the threshold T . This process is repeated for a varying number of readouts m. In Supplementary Fig. 5d the average fidelity is shown as a function of the number of readouts m. The optimum is found for m = 16 readouts with a threshold T = 7. We obtain a combined initialisation and readout fidelity of F = 94(1)%. For the results in the main text we use m = 18 and T = 7 which gives, within error, the same fidelity. The sequence used to calibrate the pair A and B parity measurement. We select on > 14/20 counts in the first 20 readouts and apply a π/2 around x to obtain a mixed state in the antiparallel subspace. Then another 20 parity readouts are used to initialise pair A and B in even parity (> 15/20) or odd parity (< 2/20) states. Finally m parity readouts are performed that we aim to calibrate. b. Conditional histograms for |a (green) and |b (blue) for m = 18 readouts. The optimal threshold of T = 7 is indicated. c. Combined initialisaton and readout fidelity of the individual states |a (green) and |b (blue) for m = 18 readouts. d. We vary the number of readouts m and find the optimal threshold for each m. The average fidelity (equation (S13)) is plotted against the number of readouts m and the corresponding optimal threshold is indicated for that number of readouts.

VIII. CALIBRATION OF THE SPIN MEASUREMENT FOR PAIR C
For pair C we optimize the spin readout to optimally distinguish |a = 1 √ 2 (|⇑ + |⇓ ) and |b = 1 √ 2 (|⇑ − |⇓ ) (Extended Data Fig. 2). To initialise these states, we first initialise the pair in the antiparallel subspace and then initialise the spin state ( Supplementary Fig. 6a). We set N 0 = 9 in the p = 10 parity readouts and set N a = 6 and N b = 1 for the k = 7 spin readouts.
In Supplementary Fig. 6b the histograms for |a (green) and |b (blue) are shown for m = 6 readouts. We calculate the combined initialised and readout fidelity of |a and |b using equation (S13) for varying thresholds T (Supplementary Fig. 6c). In Supplementary Fig. 6d we vary the number of readouts m and plot it against the average fidelity. The indicated threshold is the one that gives the maximum average fidelity. The optimum is found for m = 6 with a threshold T = 3. We obtain a combined initialisation and readout fidelity of 90 (2) Figure 6. Spin readout calibration of pair C. a. The sequence used to calibrate the pair C spin readout. We select on > 9/10 in the 10 parity readouts to obtain a mixed state in the antiparallel subspace. Then we initialise 1 √ 2 (|⇑ + |⇓ ) (> 6/7) or 1 √ 2 (|⇑ − |⇓ ) (< 1/7). Finally we do m spin readouts that we calibrate to optimally distinguish these states. b. Conditional histograms for 1 √ 2 (|⇑ + |⇓ ) (green) and 1 √ 2 (|⇑ − |⇓ ) (blue). The optimal threshold of T = 3 is indicated. c. Combined initialisation and readout fidelity of 1 √ 2 (|⇑ + |⇓ ) (green) and 1 √ 2 (|⇑ − |⇓ ) (blue) for m = 6 readouts. d. The average fidelity (equation (S13)) is plotted as a function of the number of readouts m and the corresponding optimal threshold for the given number of readouts.
In Supplementary Fig. 7b the histograms for |a (green) and |b (blue) are shown for m = 16 readouts. We calculate the combined initialisation and readout fidelity of |a and |b using equation (S13) for varying thresholds T (Supplementary Fig. 7c). In Supplementary Fig. 7d we vary the number of readouts m and plot it against the average fidelity. The indicated threshold is the one that gives the maximum average fidelity. The optimum is found for m = 16 with a threshold T = 7. We obtain a combined initialisation and readout fidelity of 95.9(6)%.  Figure 7. Parity readout calibration of pair C a. The sequence used to calibrate the pair C parity readout. We initialise the pair in the antiparallel (> 9/10) or parallel (< 1/10) subspace. Then we do m parity readouts that are calibrated to optimally distinguish the subspaces. b. Conditional histograms for the antiparallel (green) and parallel (blue) subspace for m = 16 readouts. The optimal threshold of T = 7 is indicated. c. Combined initialisation and readout fidelity of the antiparallel (green) and parallel (blue) subspace for m = 16 readouts. d. The average fidelity (equation (S13)) is plotted as a function of the number of readouts m and the corresponding optimal threshold for the given number of readouts.

X. VALUE OF THE DIPOLAR COUPLING X
The dipolar coupling X between two carbon spins is given by [8] X = µ 0 γ 2 c 8πr 3 (1 − 3 cos 2 θ) (S14) as defined in Methods, with γ c = 67.2828 · 10 6 rad s −1 T −1 and µ 0 = 4π · 10 −7 H/m. A nearest neighbour pair along the field, such as pair A and B, has r = a0 4 [1, 1, 1] where a 0 = 3.5668Å is the lattice constant of diamond at 3.7 K [9]. Consequently θ = 0 and we obtain X = 2π · 2062.37 Hz. This is significantly different from the observed value of X = 2π · 2080.9900(3) Hz. For pair C the theoretically predicted value is X = 2π · 186.92 Hz whereas the observed value is X = 2π · 188.33(2) Hz. Notably pair A and B have the same or very similar X values, suggesting a mechanism that is not dependent on the local environment. Furthermore, the observed change in X for pair A, B and C is consistent with a reduction of the lattice constant by ∼ 0.01Å. In this section, we analyze several mechanisms that could affect the value of X.

A. Strain
The E x and E y optical transitions for the NV considered in this work are split by ∼ 4 GHz, implying a strain of δ = 2 GHz [10]. This is transversal strain with respect to the NV-axis. Axial strain is aligned with the pair A, B axis. We take a typical strain-induced splitting on the order of ∼ 10 GHz and perform an order of magnitude estimate. 1 GPa of external stress corresponds to a 10 3 GHz splitting [11,12], so ∼ 10 GHz corresponds to 0.01 GPa. Taking a Young's modulus of ∼ 10 3 GPa we obtain a deformation of ∼ 10 −5 . If the axial strain is comparable to the transversal strain, strain cannot explain the observed increase in X, since the required reduction in the lattice constant is on the order of ∼ 10 −3 . However, a precise value for the axial strain is not known and we can therefore not rule out this mechanism.
B. Effect of 13 C isotope Pair A and B each consist of two 13 C spins surrounded by mostly 12 C spins, because of the natural 13 C abundance of 1.1%. The value for a 0 used to get to X = 2π · 2062 Hz does not take into account variations in the 13 C abundance. It has been shown that the diamond lattice constant decreases upon increasing the 13 C abundance [13]. This suggests that local variations in the diamond lattice constant due to 13 C might play a role in the observed value of X. Additional research is required to be able to quantitatively compare this microscopic effect to the measured values.
C. Frequency shifts due to the 13 C bath From equation (S11) we find that the observed frequency in m s = 0 is actually X + b 2 2X . The noise contribution δf = 1 2π b 2 2X is relatively small and therefore we quote the observed value as X in the main text. However, we can estimate the noise term using the values of b A = 2π · 13.9(2) Hz and b B = 2π · 12.5(4), obtained from the T * 2 measurement with the NV electron spin in m s = −1. Thus δf A = 0.046(1) Hz and δf B = 0.038(2) Hz. Notably, the difference is 0.009(3) Hz. It is therefore possible that the long T * 2 measurements presented in the main text are affected by this frequency difference. However the frequency changes due to the noise b 2 2X are too small to explain the difference between expected and measured value of X.

D. Electron-mediated coupling
The coupling strength measured can also be modified due to the presence of the electron spin especially in the presence of a misaligned magnetic field [14]. To study the electron-mediated coupling effects on the measured coupling strength, we consider the Hamiltonian describing a system of an NV center and a pair of coupled 13 C spins with dipolar coupling X = 2π · 2062.37 Hz. We numerically solve the full system Hamiltonian (i.e. not considering the pseudo-spin model) to obtain the eigenenergies of the system, from which we calculate the effective coupling strength. Supplementary Figure 8