Quantum probe hyperpolarisation of molecular nuclear spins

Hyperpolarisation of nuclear spins is important in overcoming sensitivity and resolution limitations of magnetic resonance imaging and nuclear magnetic resonance spectroscopy. Current hyperpolarisation techniques require high magnetic fields, low temperatures, or catalysts. Alternatively, the emergence of room temperature spin qubits has opened new pathways to achieve direct nuclear spin hyperpolarisation. Employing a microwave-free cross-relaxation induced polarisation protocol applied to a nitrogen vacancy qubit, we demonstrate quantum probe hyperpolarisation of external molecular nuclear spins to ~50% under ambient conditions, showing a single qubit increasing the polarisation of ~106 nuclear spins by six orders of magnitude over the thermal background. Results are verified against a detailed theoretical treatment, which also describes how the system can be scaled up to a universal quantum hyperpolarisation platform for macroscopic samples. Our results demonstrate the prospects for this approach to nuclear spin hyperpolarisation for molecular imaging and spectroscopy and its potential to extend beyond into other scientific areas.

where S x,y,z are the Pauli spin matrices of the spin-1 system of the NV, D = 2.87 GHz is the corresponding zero-field splitting, B 0 is the external field strength, γ NV and γ n are the gyromagnetic ratios of the NV and target spins, I (j) x,y,z are the Pauli spin matrices of nuclear spin j, A j is the hyperfine tensor describing the spin-spin interaction between the NV and spin-j, B jk is the tensor describing the magnetic dipole interaction between spins j and k; and summation over j, k refer to all spins in the environment. In Eq. 2, γ NV is defined positive, while γ n can be positive or negative depending on the species considered.

Cross-relaxation resonances
Transforming to the interaction picture and applying the rotating wave approximation to the NV-target resonances, we may restrict ourselves to the NV subspace spanned by {|0 , |−1 }, since the |0 ↔ |+1 transition is detuned from the nuclear spin transitions by approximately 5.7 GHz. We thus define the |0 ↔ |−1 NV transition frequency as ω NV ≡ 2πD − γ NV B 0 , and the target transition frequency as ω n ≡ γ n B 0 . The cross-relaxation resonances occur when |ω NV | = |ω n |. Since ω NV changes sign at the NV ground state level crossing at B 0 = 2πD/γ NV ≈ 1024 G, there are two resonance points before and after the crossing, at magnetic fields B < (γ n ) ∼ 2D/(γ NV + γ n ) and B > (γ n ) ∼ 2D/(γ NV − γ n ), respectively. Here and throughout the paper, we assume γ n > 0. A negative sign would simply swap the two resonances and their associated NV-target interaction strengths, relative to the crossing.
Near the NV level crossing, hyperfine interaction of the NV electron spin with its own intrinsic nuclear spin (associated with the nitrogen atom) causes the NV levels to mix and anti-cross, adding small corrections to the above Resonance position of NV and target nuclear spins. a Simplified depiction of the energy levels of the NV spin at the GSLAC as a function of the applied magnetic field, and the polarisation transfer mechanism when on resonance with a target nuclear spin. b Transition frequencies as a function of the applied magnetic field, showing the NV-nuclear spin resonance conditions for multiple different nuclear spin species. The black lines represent the two allowed NV transitions which act as the probe frequency. resonance conditions. A detailed description of this NV ground state level anti-crossing (GSLAC) may be found in Supplementary Ref. [1]. The relevant energy levels of the NV electron spin and target nuclear spin are depicted in Supplementary Fig. 1a. As the polarisation transfer is mediated by the NV-target hyperfine interaction, the resulting target spin state (|↑ or |↓ ) depends on which nuclear spin transition the NV spin is tuned to. Before the GSLAC (e.g. at B < ( 13 C) = 1023.9 G), CRIP will polarise nuclear spins into the |↑ state ( Supplementary Fig. 1a, left side). On the other hand, when applied at the resonance condition after the GSLAC (e.g. B > ( 13 C) = 1024.9 G, or B > ( 1 H) = 1026.2 G) the target is polarised into the |↓ state ( Supplementary Fig. 1a, right side). The NV-nuclear resonance conditions for various species of nuclear spins is shown are Supplementary Fig. 1b. One can see that some target nuclear spin species, including 1 H for example, do not have a resonance condition before the GSLAC, as the requisite NV spin transition is forbidden. After the GSLAC, however, there is no significant state mixing, which allows all spin species to be addressed via the CRIP protocol. In the following, we will consider only cases where state mixing is negligible, so that we do not have to include the NV nuclear spin in the model.

Definitions
Under this simplification the hyperfine component of the Hamiltonian becomes S · A j · I = A (j) xx S x I (j) x + A (j) xy S x I (j) y + S y I (j) x + A (j) yy S y I (j) which may be written more compactly for near-resonant cases before and after the GSLAC as respectively, where |0 , | − 1 refer to the NV spins states, and As these two transitions are spectrally distinct, from this point forward, we will simply refer to the transverse hyperfine coupling between the NV spin and spin j as A j , where it is understood that its particular functional definition (A < k or A > k ) depends on the choice of resonant field strength, B < or B > .
For a radial separation distance R j and polar angle Θ j , the effective transverse coupling rates between the NV spin and environmental spin j are given by for magnetic fields of B < , and B > , respectively, corresponding to resonant transitions either side of the NV GSLAC. For the B 0 ∼ 1000 G magnetic field strengths studied in this work, the Zeeman energies of the nuclear spins far exceed their mutual nuclear dipole couplings, hence any spin flip-flop dynamics they exhibit must be magnetisation conserving, and we may apply the secular approximation to the Hamiltonian describing their mutual dipole-dipole interaction where and r jk and θ jk are the separation distance and polar angle between spins j and k.

C Environment
For the purposes of discussion of the relative strengths of the two transitions, it is useful to average over the angular components of Eq. 6, giving for before or after the GSLAC, respectively where a = µ 0 4π γ NV γ n .
Determination of the total hyperfine field strength then proceeds via summation over R j , where we now index R j according to the expected distance from the NV to its j th nearest 13 C spin; the distribution of which is given by where n t = 1.95 nm −3 is the density of 13 C nuclei in diamond. Using this result in Eq. 9, we find A < ≈ 6.24an t = 1.43 × 10 6 rad s −1 = 228 kHz (11) A > ≈ 15.28an t = 3.51 × 10 6 rad s −1 = 559 kHz.  Figure 2. Coupling before and after the GSLAC comparison. Contours of constant coupling for after (a) and before (b) the GSLAC with 111 and 100 surfaces shown in purple. c Polarisation profile at the diamond surface for all four possible scenarios, considering an NV located 10 nm below the surface and polarising a bath of PMMA spins.

H Environment
For the 1 H bath in PMMA, the NV-environment separation distance (approximately 10 nm) is much larger than the separation distances between adjacent 1 H spins. As such, the total hyperfine field strength may be evaluated via a standard integral over the effectively continuous, semi-infinite PMMA slab.
Due to the avoided crossing in NV spin level structure before the GSLAC (see Supplementary Fig. 1), the only resonance point between the NV spin and 1 H spins occurs after the GSLAC at B > = 1026.2 G. In addition, the angular dependence of the interaction (see Eq. 6) must be considered when choosing the NV orientation relative to the surface. For instance, an NV normal to the surface of a 111 cut crystal, which is optimal as far as photon collection is concerned, would place the nearest 1 H spins in the node of the hyperfine coupling (see Supplementary  Fig. 2). Here we use an NV in a standard 100 cut diamond crystal, which gives a relatively good coupling. As such, we transform the original coordinate system, R = (x, y, z) to the new coordinate system, R = (X, Y, Z), defined by where α = arccos 1/ √ 3 ≈ 54.7 • . The hyperfine integral then proceeds as where n t = 56 nm −3 is the density of 1 H spins in PMMA, and h is the depth of the NV below the diamond-PMMA interface (measured to be approximately 10 nm in this work).
Theoretical description of CRIP protocol in the absence of NV spin dephasing In this section, we solve for the evolution of a number of special cases of Eq. ?? in order to motivate and develop a general continuum description for an environment containing an arbitrary number of spins.

Single hyperfine coupled spin
Assuming a starting density matrix of where p denotes the probability of the target spin being in its |↑ state, the probability of finding the target spin in its |↑ state at some later time, t, is given by where δ = ω NV − ω n is the detuning between the transition frequencies of the NV and the target spin. By tuning the magnetic field such that δ = 0, and choosing the evolution time to facilitate the maximum possible magnetisation exchange, t = √ 2π A ≡ τ , the probability of finding the spin in its up state becomes p(τ ) = 0, indicating that the polarisation of the NV has been transferred to the target spin. In the following, we will take δ = 0 for simplicity. The effect of quasistatic magnetic noise, which amounts to a fluctuating detuning, will be be incorporated in the continuum description by introducting the NV dephasing rate, Γ 2 .

Two independent hyperfine coupled spins
In the case of two target spins coupled to the NV, complete polarisation cannot be achieved within a single interaction. Assuming an initial state of the time-dependent |0 state population of the NV is Therefore, the maximum possible contrast corresponds to the minimum of p NV , given by when δ = 0 and Substituting τ into the expressions for p 1 and p 2 , we find the change in populations to be The terms proportional to p 1 − p 2 are mixing terms that facilitate faster polarisation of the farthest spin due to the back action of the NV spin in response to the closest spin.

Two dipole coupled spins
We now add an additional magnetic dipole-dipole interaction between the two spins. In general, this system does not exhibit a closed-form solution, however, the timescale of the dipole coupling between the two environmental spins, τ B ∼ 1/B 12 , is much longer than that of the hyperfine dynamics, τ A ∼ 1/A. Hence, we may expand the dipole-dipole dynamics for times t 1/B 12 to give For cases where both spins exist outside the strongly hyperfine coupled core surrounding the NV, we have that the two spins will thermalise on the timescale of their dipole coupling, τ B . Hence, for outside the core on timescales of τ , we have General case and continuum description

General case and continuum description
For the case of N S spins, each having hyperfine coupling A j to the NV, mutual dipole couplings B jk , and initial state, the population of the NV |0 state is at a minimum when t = τ = √ 2π A0 , and is given by where A 2 0 ≡ j A 2 j is the total hyperfine coupling field. The change in probability of spin j is given by where B jk is the magnetic dipole coupling between spins j and k, and C 2 jk ≡ A 2 j A 2 k /A 2 is the effective hyperfinemediated diffusion strength. For approximately N S ∼ 10 4 spins, the evolution of discrete spins in this system may be solved numerically, an example of which is given in Fig. 3.
Spins residing in the hyperfine core are characterised by having a stronger coupling to the NV than to all other spins in the environment, i.e. A 2 j > k B 2 jk . Typically we are interested in total times t τ (hours vs µs, respectively) for which the polarisation region has far exceeded the hyperfine core, as a result of dipole-mediated flip-flops between the polarised core and the unpolarised outer region. As such, we are not interested in the relatively fast dynamics associated with intra-core diffusion, and instead only concern ourselves with the rate at which the NV polarises spin in the core (A k ), and the rate at which this effect is communicated to the rest of the bath B jk . We thus take C jk ∼ 0.
In modelling these systems, we typically want to consider regions of at least 200 nm in size, which means of order N S ∼ 10 10 spins and ∼ 10 20 couplings for the case of PMMA, thereby making modelling of discrete spin states computationally unfeasible. We therefore map Eq. 26 to a temporally and spatially continuous description via the following: a. Time-dependent probability field -As the regions of polarisation considered in this work are larger than the hyperfine core of the environment, which itself is comprised of the order of 10 4 spins, instead of considering the timedependent populations of discrete spins, we instead consider the time-dependent probability field of the environmental spins having an average density, n t , In essence, p (R, t) is the average probability of spins in an infinitesimal volume at position R from the NV being in state |↑ at time t. b. Probability field dynamics -In the discrete description, the probability of finding spin j in its |↑ state is monitored at discrete multiples of the optimal dark time, τ . Ignoring probability diffusion, this would behave as a geometric series. As we are interested in times t τ , we map changes in p j over time to c. Hyperfine coupling field -As with the population field, we also map the discrete hyperfine couplings to a continuous field whose strength is determined by its position relative to the NV, R.
d. Total hyperfine coupling strength -The summation over all hyperfine couplings is mapped to an integral, e. Probability diffusion -We put p j → p (R, t) and p k → p (R + r, t), and discretise the Laplacian operator to get (31) f. Additional sources of relaxation -To account for additional sources of environmental spin-lattice relaxation, where applicable, we add an additional probability sink in which any probability outside equilibrium decays at a rate Γ SL = 1 TSL , whose rate is given by The differential equation describing the evolution of this system is thus given by where is the effective source, or cooling, coefficient at position R relative to the NV, and is the effective probability diffusion coefficient related to the intra-bath interactions B jk . Recasting Eq. 33 in terms of the polarisation P (R, t) = 1 − 2p(R, t), we obtain Equation (1) of the main text.  Table 1. Summary of parameters used to describe the different scenarios examined in this work: the intrinsic 13 C bath in diamond (first two column, for before and after GSLAC resonances, respectively), and the 1 H bath formed by a PMMA slab on diamond with a 10-nm-deep NV (third column). Before GSLAC After GSLAC For cases in which the dephasing rate of the NV spin, Γ 2 , is greater than the hyperfine coupling to the desired target, i.e. Γ 2 A 0 , the population of |0 state of the NV spin will exhibit an exponential time dependence. As such, for this regime we instead make the replacement where the details of this process are discussed in the following section. Experimentally, Γ 2 is obtained via the T 1 -relaxometry spectrum (with the interleaved sequence), since Γ 2 defines the width of the cross-relaxation resonances [2]. We note that Γ 2 is expected to decrease upon polarisation of the 13 C bath. However, the effect is relatively small (see section Supplementary Note 2), and as such in the calculations shown in main text Fig. 2 we kept Γ 2 constant and equal to the off-resonance value (i.e. unpolarised case). For the 1 H case, the dephasing is dominated by surface effects, and is therefore unaffected by the 1 H polarisation.
A summary the parameters associated with the three systems studied is given in Table 1, detailing: the diamond 13 C bath before the GSLAC ( Supplementary Fig. 8) and after the GSLAC (main text Fig. 2), and the PMMA 1 H bath after the GSLAC (main text Fig. 3).
The reconstructed radial profiles for the 13 C polarisation in the main text are shown in Supplementary Fig. 4a for both before and after the GSLAC. In a similar fashion the PMMA polarisation can be reconstructed. 2D slices of the polarisation extent in the PMMA for different polarisation times are shown in Supplementary Fig. 4b.
Theoretical description of CRIP protocol in the presence of NV spin dephasing The above development outlines the behaviour of the spin bath under CRIP for cases in which NV spin dephasing is not appreciable over the interaction time, τ . In what follows, we discuss how this behaviour is modified in the opposite regime. To motivate the discussion, we recall the semi-classical result from Supplementary Ref. [2], which gives the response of an NV spin prepared in the |0 state to an effective hyperfine field of magnitude A, and transverse spin relaxation (dephasing) rate, Γ 2 , where In the semi-classical sense, A refers to the coupling of the effective magnetic field felt to the transverse components of the NV spin due to all target spins the environment. The magnitude of this field depends on the orientation of the individual spins, and is related to the population of finding them in the | ↑ state, p j , and the hyperfine coupling of the NV spin to spin j, A j , via Naturally a bath in which every spin is polarised in its | ↓ state will not be able to flip the NV spin, leading to no relaxation effect; however a bath in which all spins are in the | ↑ state will have the greatest effect. For cases in which the baths are sufficiently large that the spins effectively comprise a continuum, this becomes which is the case for both the native 13 C and PMMA 1 H baths considered in this work. In the case of the latter, the distance between the NV and the nearest target spin is 10 nm, whereas their average separation is only a few angstrom. Whilst this is not the case in general for the former, the target spins close to the NV polarise so quickly under CRIP that the effective standoff between the NV and the unpolarised spins in the target bath very quickly becomes larger than the separation between target spins. In what follows, we derive a quantum mechanical description of the NV spin relaxation, and show that it agrees with the classical result for the cases considered in this work. As with the no dephasing case, we consider examples of only a few spins, and show how this generalises to large target baths.

Single target spin
For the case of a single spin, we follow Supplementary Ref. [2] and consider the evolution of the density matrix as described by where σ z = |0 0| − | − 1 −1| operates on the NV spin; subject to initial state ρ 11 = p 1 ; ρ 22 = 1 − p 1 , with all other ρ jk = 0. This gives rise to the following differential equations, which, in the limit of Γ 2 2 A 2 1 , results in an NV |0 population of where, as in the classical case, Two target spins For the case of a two spin target bath, there are now three states that couple to the |0 state of the NV spin: |↑ |↑ , |↑ |↓ , and |↓ |↑ . Following the same process as for the single spin case, but for these three transitions, we have where p 1,2 is the probability of spin 1,2 being in its | ↑ state. From this, we see that the populations of the four possible basis states decay with rates of respectively. Thus we see that the decay rate of the NV spin depends on the initial state of the target bath, and consequently, on the individual polarisations of the bath constituents.

General case
Extending this to the many spin case, we see that the NV |0 state population in response to the target being initialised in the product basis state |σ 1 σ 2 . . . , where σ k = 1, 0 corresponds to | ↑ , | ↓ respectively, is given by where, is the decay rate out of product state |σ 1 σ 2 . . . , and q(σ 1 , σ 2 , . . .) is the initial probability of being in that state. In order to relate this process to the individual spin probabilities, p j , we note that for a random initial state, each spin has an associated initial density state given by where the initial full target density state is ρ(0) = k ρ k (0). From this, we see that The full NV |0 state population is just the sum over all possible initial states associated with Eq. 50, For cases in which the NV-target standoff is much larger than the average distance between the target spins, many spins will couple to the NV spin simultaneously and with similar strength. The contribution from each individual spin ∼ A 2 j /2Γ 2 will therefore be necessarily small as compared with the total NV spin relaxation rate ∼ j A 2 j /2Γ 2 , meaning that we may formally expand both the exponential term and the logarithm above to first order in A 2 j / A 2 j , giving in agreement with the classical result (Eqs. 37-40).

Scale-up strategy
Having analysed the CRIP protocol above for the case of a single NV spin, we now discuss how this technique may be scaled up for the purpose of polarising macroscopic quantities of arbitrary nuclear spin labels for clinical applications such as MRI imaging.

Polarisation rates
From the analysis above, the maximum rate (spins/s) at which polarisation is delivered to the target spin ensemble from a single NV is equal to the effective hyperfine field, where n t is the numerical density of the target spins, h NV is the distance of the NV from the diamond/target interface; and a = µ0 4π γ NV g T µ N , where g T is the g-factor of the target species, and µ N is the nuclear magneton. For a planar NV density within the diamond substrate of σ, the total polarisation delivery rate is (57)

Time-dependent polarisation
Assuming that the spatial diffusion of target spins is sufficiently fast to ensure perfect mixing (i.e. uniform polarisation throughout the cell on timescales of ∼ 1/A), the polarisation of the target spins is governed by where P is the average polarisation of the target spins, and T SL is their spin-lattice relaxation time. Solution of this equation yields where N ≡ n t × Area × h cell is the total number of target spins in the cell, and h cell is the height of the cell.

Optimal geometry
From Eq. 59, we see that the steady-state polarisation and characteristic polarisation time are given by respectively. It is clear that these expressions define a critical effective number of spins, N crit ≡ RT SL , that may be polarised before that polarisation is lost to spin-lattice relaxation effects. As such, any delivery system of spinpolarised contrast agents necessitates that N crit N . As both of these quantities are proportional to the area of the NV array, we may rewrite this constraint in terms of the coupling parameters h N V and cell height h cell to give where where N m is the number of target spins in the target molecule, and we have assumed a target concentration of 1 M. This allows us to define an effective length scale by which to characterise the suitability of candidate systems to this scheme. From Eq. 59, the time required to achieve a desired polarisation, P , is and the volume delivery rate of polarised contrast agent is where d is the dilution factor, typically of order 10 3 to reach mM concentrations. In the limit of RT SL N , this becomes Notice that T SL does not enter into Eq. 66, as we have assumed RT SL N , and hence h cell h crit . Values of these quantities are given in Table 2 for various candidate contrast agents. Supplementary Figure 5. Cross-relaxation between a single NV and the 13 C bath. a T1-relaxometry spectrum performed around the GSLAC (middle feature) with two signals (outer features) corresponding to the 13 C resonances, with a wait time τ = 5µs. The black and red lines (|0 and |−1 ) are obtained with a CRIP and CRIP −1 interleaved pulse sequence. The pulse sequence acts to prevent polarising the bath and allows the NV-13 C interaction to be measured. The solid line is a guide to the eye. b The difference between the |−1 and |0 state (∆PL = PL(|−1 )− PL(|0 )) against the transition frequency measured in ODMR for before the GSLAC (black) and after (blue). c Short time spin dynamics using the same interleaved pulse sequence and varying the interaction time, τ . The measurement is taken on resonance with the 13 C bath (after GSLAC), revealing flip-flop interaction between the NV and the 13 C bath with a total hyperfine coupling of 250 ± 40 kHz.

C detection and polarisation
Before proceeding with the polarisation the 13 C spin bath, the T 1 -relaxometry spectrum of the 13 C bath was first obtained using the interleaved sequence, while scanning the magnetic field across the GSLAC. The full spectrum is shown in Supplementary Fig. 5a for the NV used in main text Fig. 2. The spectrum resolves three peaks, the two outer peaks correspond to the signal from the 13 C (before and after the GSLAC), the inner peak is an intrinsic feature of the GSLAC Supplementary Ref. [1]. The GSLAC feature of the |0 state (black, readout from the CRIP sequence) is narrower than that of the |−1 state (red, readout from the CRIP −1 sequence) due to experimental errors in implementing a π-flip in this region. This is caused by the rotating wave approximation breaking in this regime. This is because the driving field term in the Hamiltonian becomes similar in size to the quantisation axis itself. This needs to be considered for pulse experiments near the GSLAC as it can lead to artefacts in the measurements.
To elucidate the origin of the signal the PL difference (∆PL) between the |0 and |−1 states versus the transition frequency (ω NV ) is shown in Supplementary Fig. 5b. The expected transition frequency, in the weakly coupled regime, is ω C = γ C B z ≈ 1.09 MHz. There is a slight shift in the frequency measured which is possibly due to a hyperfine coupling between the NV and the nearest 13 C spins that remain slightly polarised. Measuring the short time dynamics of the NV spin on resonance shows a coherent evolution caused by the hyperfine coupling to the surrounding 13 C spins. The frequency of the oscillation corresponds to the sum of the total hyperfine coupling of the bath to the NV (A 2 0 = k A 2 k , see section S1) but will be dominated by the nearest 13 C. The damping is given by fluctuations in the spin configuration of the remaining bath spins, which causes the detuning, δ, to randomly vary. This is shown in Supplementary Fig. 5c and has an approximate coupling strength of 250 ±40 kHz.
The NV-bath coupling strength depends on the particular configuration of the bath around a given NV. We therefore expect a variability in the spectrum and degree of polarisation depending on the specific location of the 13 C spins in relation to the central NV spin. The 13 C polarisation has been tested over a variety of diamonds including deep NVs near surface NVs in bulk diamond, and NVs in micro-pillars. All of the different samples and NVs have shown the capability to polarise the 13 C spin bath provided that the transition is observable around the GSLAC, which requires the spin dephasing time of the NV, T * 2 , to be sufficiently long Supplementary Ref. [1]. The full CRIP and interleaved spectra (including the GSLAC feature) for the NV used in main text Fig. 2 are shown in Supplementary Fig. 6a. Also shown in Supplementary Fig. 6b are the spectra from another NV. The time dynamics for this NV on resonance with the 13 C bath is shown in Supplementary Fig. 6c, revealing a similar hyperfine coupling. Supplementary Figure 6. Polarising the 13 C bath. a The 13 C spectrum (orange) obtained using the depolarisation sequence and the polarised spectrum (blue) using CRIP. This is the same NV as the one in main text Fig. 2. b Spectrum of 13 C around the GSLAC with polarisation (blue) and without polarisation (orange) on a different NV. With the interleaved sequence, both NV initialisation states are shown, |0 (upper orange) and |−1 (lower orange). c Time dynamics of the NV-bath interaction at the after-GSLAC 13 C resonance.

C polarisation dynamics
The polarisation dynamics was investigated for the NV used in main text Fig. 2, when on resonance with the 13 C bath (after the GSLAC). To this end, a series of N = 30 pulses was performed with a π pulse (i.e. CRIP −1 ), which polarises the 13 C in the |↑ state. This is followed by N pulses without the π pulse (i.e. CRIP) polarising the 13 C bath in the |↓ state. The pulse sequence is shown in Supplementary Fig. 7a. The PL readout for each laser pulse as a function pulse number for three interaction times τ = 3, 5, 10 µs is shown in Supplementary Fig. 7b. The resulting polarisation curves were fitted with exponentials (PL = a exp (N/T pol ) + y 0 ) in order to determine the characteristic polarisation time. The result of the fits is shown in Supplementary Fig. 7c and shows a clear trend of decreasing the  number of polarisation pulses required as the interaction time is increased. It is important to note that due to the relatively short interaction time used in these measurements, the contrast obtained is dominated by the few nearest spins and is therefore an indication of the polarisation of nearby 13 C spins only, which appear to be efficiently polarised by the CRIP protocol after a few tens of µs. More remote 13 C spins can be probed using longer τ times, as done in main text Fig. 2b, but the significantly longer acquisition time for these measurements prevents a systematic study of the polarisation dynamics as done in Supplementary Fig. 7.
Finally in order to test the robustness of the sequence to the resonance position we performed a magnetic sweep across both resonances and the GSLAC, shown in Supplementary Fig. 7d. Four features are present, two related to 13 C (outside peaks) and two related to the GSLAC (inside peaks). The width of the polarisation region is roughly 0.1 -0.2 MHz. The width of this polarisation region is governed by the T * 2 of the NV spin. A comparison of the number of pulses required to polarise the nearest 13 C spins on both sides of the GSLAC shows a distinct difference, which is explained by the difference in the angular dependence of the dipole-dipole coupling before and after the GSLAC.

Spatial extent of 13 C polarisation
To monitor drifts in the magnetic field (caused by thermal fluctuations), the full-length T 1 measurements used to estimate the extent of polarisation (main text Fig. 2b) were accumulated while monitoring the NV frequency via an ODMR spectrum taken every 30 minutes. In general, the NV frequency remained on resonance with the 13 C transition frequency (1.09 MHz) for several hours, before it drifted away by more than the NV intrinsic linewidth (≈ 200 kHz). For the data in main text Fig. 2b, the total acquisition time was 5 hours, during which the NV was on resonance within the NV linewidth, i.e. ω NV in the range 1.0-1.2 MHz. This ensures nearly optimal interaction between NV and 13 C bath, which is important not only to polarise the 13 C bath but also to accurately probe the remaining NV-13 C interaction which we used to estimate the polarisation extent. For the reference measurements off resonance, we repeated this procedure by maintaining the NV at lower frequency (< 0.9 MHz), and then at higher frequency (> 1.3 MHz), for more than 10 hours each. The final curves were then fitted with exponentials to determine the total relaxation rate Γ tot , which is expressed as Γ tot = Γ bg + Γ CR on resonance, and simply Γ tot = Γ bg off resonance. The full-length T 1 curves are shown in Supplementary Fig. 8a for the NV used in main text Fig. 2 near the after-GSLAC resonance, under the three conditions explained above. The extracted total decay times are shown in Supplementary  Fig. 8b for both before and after GSLAC cases. The data before the GSLAC (green) shows a slight change in the decay rate measured on resonance, i.e. Γ CR ≈ 62 Hz. After the GSLAC however, there is no discernible difference outside of error, Γ CR Γ error = 19 Hz, given by the uncertainty of the fit parameter.
Effect of 13 C polarisation on NV spin coherence As the limiting factor of the dephasing time of the NV spin (T * 2 ) is noise from the 13 C bath (for the deep NVs considered here), it is expected that polarisation of the 13 C spins will result in an increase in the T * 2 of the NV. This effect was observed in previous works that employed a Hartmann-Hahn resonance to transfer polarisation Supplementary Ref. [3,4]. To measure T * 2 , one usually performs a free induction decay (FID, or Ramsey) experiment. However, such a measurement cannot be done at the cross-relaxation resonance because it would interfere with the flip-flop dynamics Supplementary Figure 9. Effect of 13 C polarisation on the free induction decay. a Pulse sequence for CRIP polarisation at the resonant frequency followed by an FID measurement at a different frequency, using a magnetic field offset pulse. b Schematic of the experimental set up with the coil for the field offset. c Example ODMR of the NV during the CRIP phase (left) and during FID measurements (right). d FID signal as a function of NV frequency as measured during the CRIP phase (top axis) and during the FID phase (bottom axis). The FID wait time is τ = 4 µs. A detuning of 0.2 MHz was applied when driving the NV in the FID measurement. e FID curves obtained from a similar sweep as in d. The two curves show the case where the CRIP phase is on resonance (blue, polarisation) and off resonance (black, no polarisation). f Long FID measurement taken at a fixed frequency with the CRIP phase on resonance (blue, polarisation) and off resonance (black, no polarisation). The applied driving detuning was different in the two case in order to obtain the same oscillation frequency.
of the NV-13 C interaction. A solution is to polarise the bath at the resonant field then quickly shift the field to off the resonance and perform the T * 2 measurement, following the sequence shown in Supplementary Fig. 9a. To implement this, a coil was attached to the board holding the diamond (schematic shown in Supplementary Fig. 9b). A TTL pulse was used to generate a DC field from the coil to act as a quick field switching mechanism, allowing the polarisation to occur at the resonance frequency and the measurement to be performed at another frequency, as illustrated by the ODMR spectra shown in Supplementary Fig. 9c. The field offset thus applied was B offset ≈ 0.7 G, turned on while a series of polarisation pulses was implemented. The offset was then turned off and the FID measurement performed.
Using the same NV as in main text Fig. 2, we first performed a single-τ FID measurement as a function of the magnetic field, across the 13 C resonance (after GSLAC). The resulting spectrum is shown in Supplementary Fig. 9  d, where ω NV (measured by an ODMR spectrum at each magnet position) was scanned from -1.0 MHz to 1.5 MHz during the CRIP polarisation phase (with the field offset on), corresponding to 1 MHz to 3.5 MHz during the FID measurement (field offset off). For the polarisation phase, we used N = 10 CRIP pulses, and for the FID we used τ = 4 µs, and a driving frequency detuned by 0.2 MHz from the NV frequency so as to induce an oscillation in the FID response. The CRIP+FID sequence is repeated many times at each magnet position, so that the 13 C bath should be efficiently polarised when on resonance. The data reveals a feature centred on the NV-13 C resonance, i.e. when ω NV = 1.1 MHz during the polarisation phase. To interpret this feature, we measured full-length FID curves while scanning the magnetic field. Supplementary Fig. 9e shows the resulting curves with the CRIP at the resonance or off the resonance. The main difference between the two FID curves is in the frequency of the oscillation, which differ from each other by 58 (27) kHz. This difference is attributed to a change in the DC magnetic field seen by the NV, induced by polarisation of the nearest 13 C spins. However, the envelope of the oscillation shows a similar decay time (T * 2 ) in both cases, which means the 13 C polarisation was not sufficient to significantly reduce the magnetic noise seen by the NV. To increase the extent of the 13 C polarisation, we performed a longer measurement at just two fixed frequencies, on and off resonance. In addition, we adjusted the driving frequency detuning such that the FID oscillation showed the same frequency with and without polarisation, for ease of comparison. The results are shown in Supplementary  Fig. 9d, and reveal a ∼50% increase of T * 2 from 3.6(5) µs to 5.5(5) µs upon polarisation, similar to previous results with Hartmann-Hahn based polarisation Supplementary Ref. [3,4]. The improvement is limited by the stability of the resonance condition during CRIP, as the switching of the magnetic field to bring the NV into resonance is non ideal with our current setup (in particular, showing significant overshoot).
Another way to probe a change in T * 2 is by examining the width of the resonance feature in the T 1 -relaxometry spectrum, since it is directly given by the dephasing rate 1/T * 2 Supplementary Ref. [2]. To test this, we implemented the sequence of Supplementary Fig. 7a, with N = 10 and a varying number of CRIP −1 pulses, N . The results are shown in Supplementary Fig. 10. When N = 0, there is no signal as the 13 C bath is polarised. As N is increased, the width of the feature increase from 55(12) kHz for N = 2 to 90(13) kHz for N = 6. No difference in the width of the feature was observed for larger values of N . The FWHM was extracted from a Lorentzian fit with a slight linear background to account for the tail of the GSLAC feature.

H polarisation
The polarisation of the 1 H spin bath in PMMA was performed on multiple shallow NVs. Supplementary Fig. 11a shows a spectrum showing the polarisation effect for a different NV (NV2) than that used in main text Fig. 3. Similarly, no 1 H signal was observed for the polarised sequence (blue) using an interaction time τ = 20µs, whereas a signal is observed when the depolarisation sequence is used. We note that no measurable 1 H signal was observed in the absence of PMMA (removed with dichloromethane) even with the depolarisation sequence, and that the polarisation effect was observed again after reapplying PMMA. This suggests that the 1 H spins we detect and polarise are mainly from the PMMA, although we cannot exclude some contribution from contaminants trapped under the PMMA, e.g. water.
Like for 13 C, the full T 1 measurements at the 1 H resonance were obtained by regularly monitoring the NV frequency Supplementary Figure 11. Additional 1 H polarisation measurements. a Spectrum of 1 H unpolarised (orange) and polarised (blue) for NV2. b Full-length T1 measurements for NV3 on resonance unpolarised (orange), polarised (blue) and a background off resonance measurement (green).
to make sure it remains on resonance within the NV linewidth (≈ 500 kHz here). Because the T 1 time of shallow NVs is shorter than for bulk NVs, we could limit the total acquisition time to about 30 minutes per T 1 curve. To compare with the background T 1 , the measurements were done on resonance (ω NV ≈ 4.4 MHz) and off resonance (ω NV ≈ 3.0 MHz). Supplementary Fig. 11b shows full-length T 1 curves from an additional NV (NV3). The data fit resulted in a ratio of relaxation rates of Γ unpol CR /Γ pol CR = 2.4(3) which is consistent with the measurement from NV1 presented in the main text of 2.8(3). Error is given by the uncertainty in the decay rate fit parameter.