All-optical nuclear quantum sensing using nitrogen-vacancy centers in diamond

Solid state spins have demonstrated significant potential in quantum sensing with applications including fundamental science, medical diagnostics and navigation. The quantum sensing schemes showing best performance under ambient conditions all utilize microwave or radio-frequency driving, which poses a significant limitation for miniaturization, energy efficiency, and non-invasiveness of quantum sensors. We overcome this limitation by demonstrating a purely optical approach to coherent quantum sensing. Our scheme involves the 15N nuclear spin of the Nitrogen-Vacancy (NV) center in diamond as a sensing resource, and exploits NV spin dynamics in oblique magnetic fields near the NV’s excited state level anti-crossing to optically pump the nuclear spin into a quantum superposition state. We demonstrate all-optical free-induction decay measurements—the key protocol for low-frequency quantum sensing—both on single spins and spin ensembles. Our results pave the way for highly compact quantum sensors to be employed for magnetometry or gyroscopy applications in challenging environments.


INTRODUCTION
Spin based quantum sensors can be employed to measure a wide range of relevant physical quantities, including magnetic [1] or electric fields [2], temperature [3], or rotary motion [4].This abundance of potential observables, combined with their high sensitivity at the nanoscale makes quantum sensors highly interesting for many fields of application, such as life sciences [5], geological sciences [6,7], navigation [8] and material sciences [9].
Nitrogen-Vacancy (NV) centers in diamond (Fig. 1a) are a particularly promising platform for such spin based quantum sensing applications, because they host a single electron spin [10] with long coherence times [11,12] even at room temperature [13].Upon optical excitation with green light, the NV center emits spin-dependent red photoluminescence (PL) [14], which enables all-optical electron spin readout.At the same time, such optical excitation pumps the NV electron spin [15,16] into a specific spin eigenstate, enabling all-optical spin initialization.Time-varying (AC) driving fields, mostly in the microwave (MW) or radio-frequency (RF) domain can then be used to coherently control the spin, and create superposition states for sensing.This combination of optical initialization, readout, and coherent spin manipulation by AC driving fields form the basis of almost all established spin based approaches to sensing [17].
The NV electron spin is inherently coupled to the nuclear spin of its Nitrogen atom -a spin which exhibits significantly longer coherence times compared to the NV electron spin [18,19] and therefore provides another interesting resource for quantum technology applications.Specifically, nuclear spins have been exploited as a quan-tum register for quantum communication [20] and enhanced spin readout techniques [21], but they also offer interesting opportunities for sensing, be it for magnetometry [18,[22][23][24] or for gyroscopy [25,26].
Many Nitrogen spin based quantum sensing schemes in diamond rely on the resonant coupling of the NV spin and the nuclear Nitrogen spin at a magnetic field of about 500 G [27,28], where spin flip-flop processes occur in the NV's orbital excited state at the excited state level anticrossing (ESLAC).It has been shown that optical pumping in the vicinity of the ESLAC results in nuclear spin hyperpolarization [29,30], and that -by virtue of the same mechanism -the NV center shows a nuclear spinstate dependent rate of transient PL [28].As a result, optical pumping close to the ESLAC enables both alloptical readout of the nuclear spin state and initialization into a nuclear spin eigenstate, which, together with RF driving, forms the basis for nuclear spin based sensing schemes [28].
However, the ubiquitous need for AC coherent driving in spin based quantum sensing is a severe limitation for many applications.Specifically, such AC driving fields can adversely affect investigated samples and for integrated or portable sensing devices, their delivery increases power-consumption and overall complexity, and thereby size and cost of the system.Recent experiments have demonstrated microwave-free NV magnetometry schemes that are based on sharp changes in NV PL at level anti-crossings, that occur at specific magnetic fields [31,32].While avoiding the need for MW or RF delivery, these approaches do not exploit quantum coherence (and are therefore limited in sensitivity) and are furthermore highly vulnerable to background drifts in NV PL.
Here, we present a novel method for coherent, microwave-free quantum sensing using the 15 N nuclear spin of the NV center in diamond.Our approach is based on optical driving of the NV center near the ESLAC in the presence of a small, static magnetic field transverse to the NV symmetry axis denoted by the unit vector e z (see Fig. 1a).As we will show, such a small transverse magnetic field component has a striking effect, in that optical pumping prepares the 15 N nuclear spin in a coherent superposition state within the NV's ground state spin manifold.Following optical pumping, this initialization leads to Larmor precession of the nuclear spin about an effective magnetic field; a precession we directly monitor via nuclear spin-state dependent PL [28].Figure 1b shows an example of such an all-optical nuclear free induction (FID) measurement, obtained using the pulse sequence depicted in Fig. 1c.Data for this work were recorded on a home-built confocal optical microscope (see methods) with magnetic field control; here at a magnetic field of strength The negatively charged NV center possesses an electron spin S = 1 quantized along the NV symmetry axis e z .For NV centers formed by 15 N (denoted as " 15 NV" in the following), the Nitrogen nucleus exhibits a spin I = 1/2.The Hamiltonian for the orbital ground (gs) and excited state (es) of such a 15 NV can be expressed as Ĥgs,es ) where Ŝ and Î are the NV electron and nuclear spin operators, γ S = 2.8 MHz/G and γ I = 431.7 Hz/G are the respective gyromagnetic ratios, D gs 0 = 2.87 GHz and D es 0 = 1.42 GHz are the zero-field splittings, and B ext is the applied magnetic field.The hyperfine coupling tensor A gs,es has two independent components A gs = 3.03 MHz and A gs ⊥ = 3.65 MHz for the ground state; and A es = −57.8MHz, A es ⊥ = −39.2MHz for the excited state [33,34].This Hamiltonian is conveniently expressed in a basis of spin eigenstates {|m S , m I }, where m S and m I are the magnetic quantum numbers associated with Ŝz and Îz .The coherent FID dynamics that are studied in this work (Fig. 1b) can be completely encompassed by a reduced subspace spanned by {|0, −1/2 , |0, +1/2 }.

Derivation of an effective Hamiltonian.
In the following we explain the nuclear precession data presented in Fig. 1b by discussing how the presence of the transverse magnetic field component B ⊥ significantly affects the system's optical pumping and subsequent FID dynamics.We start by calculating an effective Hamilto- nian for the 15 N spin in the m S = 0 ground state subspace, {|0, −1/2 , |0, +1/2 }, using Van Vleck perturbation theory [35] (see the supplementary information, Section I, for a detailed derivation).Without loss of generality, we set the transverse magnetic field to point along the unit vector e x .We then obtain the effective Hamiltonian where denotes the correction to the diagonal elements caused by mixing between states of different m S , and is the corresponding correction to the off-diagonal elements.We note that such an effective hyperfine Hamil-tonian for 15 NV's has been discussed earlier as a perturbation in B ⊥ in the limit B z D gs 0 [36], or as a perturbation in A gs ⊥ [28], but never as a perturbation in both simultaneously as we present it here.Additionally, exact analytic expressions for ν ⊥ have previously been derived in [37].
Hamiltonian H m S =0 eff yields that the 15 N nuclear spin in the NV ground state is quantized along an effective magnetic field B eff .Diagonalization of where Interestingly, ϑ is significantly larger than the misalignment angle Φ between B ext and e z , and ϑ has a sign opposite to Φ due to the negative sign of γ I (c.f.Fig. 1d).Finally, note that according Eq. ( 4), ν ⊥ = 0 when B ⊥ = 0, which in turns causes ϑ = 0.In this case, both B eff and B ext are aligned with the NV symmetry axis.Analysis of 15 N spin dynamics.We now discuss how the presence of B ⊥ affects the 15 N nuclear spin dynamics and enables all-optical initialization into a nuclear spin superposition state.The use of 15 N is key to this, since it does not have a quadrupolar spin splitting which would prevent any nuclear Larmor precession.It is clear from Eq. (4) that when B ⊥ = 0, the 15 N nuclear quantization axis depends sensitively on the hyperfine coupling parameter A ⊥ , and the splitting of the involved spinlevels.Therefore, the ground and excited state nuclear spin quantization axes are in general different, because the hyperfine parameters differ in both magnitude and sign between the two cases.This difference in 15 N quantization axes results in optical pumping of the nuclear spin into a state that does not correspond to an eigenstate of the effective ground state Hamiltonian H m S =0 eff .To be specific, optical pumping accumulates NV excited state population in the eigenstate with the largest m S = 0 character (i.e. the state |ψ for which ψ| Ŝz / |ψ is closest to zero), since this state has the lowest probability of shelving into the NV's singlet state -we denote this state as 0es .By the same argument, 0es is also the "brightest" state in that it yields the largest rate of emission of NV fluorescence photons.For state 0es , the expectation value of the nuclear spin lies along the vector e init = 0es Î/ 0es .Vector e init therefore defines both the direction along which the nuclear spin is initialized under green illumination, as well as the measurement axis for optical readout of the nuclear spin.For B ⊥ = 0, e init is not collinear with B eff , and thus optical pumping will initialize the 15 N nuclear spin in a state that is not an eigenstate of H m S =0 eff .Disengaging green laser excitation after optical pumping will therefore result in precession of the 15 N nuclear spin around B eff .Finally, note that for NV centers formed with 14 N, no Larmor precession occurs because its quadrupolar splitting locks B eff onto the NV axis.
Comparison with numerical results.To further verify this picture, we developed a numerical model that simulates the dynamics of the 15 NV system during and after optical pumping.The model is based on classical rate equations for the optical pumping process [16], coupled with master equations describing the quantum-mechanical evolution of the system's density matrix within each relevant orbital manifold: the orbital ground and excited states as well as the singlet state (see further details in the supplementary information, Section II.1).
In Fig. 1d we summarize the numerical and theoretical results in a Bloch sphere representation of the 15 N nuclear spin dynamics for the same magnetic field that was used to obtain the experimental results in Fig. 1b.The effective field B eff is calculated numerically through exact diagonalisation of the ground state Hamiltonian Ĥgs , and e init is calculated via diagonalization of the excited state Hamiltonian Ĥes .The nuclear spin density matrix ρnuc init following optical pumping is obtained by propagating the system density matrix ρ for 3 µs of laser excitation, followed by a 50 ns dark time (to let the system relax fully to the ground state), and by subsequently taking the trace over the NV electron spin degrees of freedom.The nuclear spin precession dynamics is then described by propagating ρ under the influence of Ĥgs .
We make two observations that underline the excellent agreement of this numerical model with our analytical discussion.First, the orientation of B eff obtained from numeric diagonalization of the full ground state Hamiltonian Ĥgs shows perfect agreement with the analytical prediction from Eq. (2) (see further details in the supplementary information, Section II.3).Second, the initial nuclear spin direction Tr( Î ρnuc init ) obtained from our numerical model is perfectly collinear with e init , as long as we set the intersystem crossing rate for the m S = 0 states to zero (see the supplementary information, Section II.3).Both observations strongly support the validity of our model and its applicability to quantitatively describe our data.
All-optical nuclear 15 N precession The presented theory framework allows us to further analyze the data presented in Fig. 1b.The observed FID oscillation frequency was determined by least-squares fitting to f nuc = 251.18± 0.12 kHz, in good agreement with Eq. ( 5), which yields f nuc = 252.71kHz for the experimental conditions |B ext | = 540 G and B ⊥ = 10.6 G.The small remaining discrepancy can be assigned to uncertainties in controlling the tilt angle Φ, and determining the exact field components B ⊥ and B .To demonstrate that the observed oscillations indeed originate from nuclear spin precession, we repeated the same experiment at fixed angle Φ, while varying |B ext |. Figure 1e shows the resulting, near-linear dependance of f nuc on |B ext | and the excellent agreement with the predictions of both Eq. ( 5) and the numeric model.The enhancement of f nuc over the bare Larmor frequency results from the terms ν z and ν ⊥ in Eq. (5).
Angle (Φ) and field (|B ext |) dependance of the all-optical 15 N FID signal.The requirement of applying a transverse magnetic field B ⊥ to obtain an observable, all-optical 15 N FID signal motivates the question of how the FID readout contrast C depends on both Φ and |B ext |. Figure 2a shows single NV data, where we determined C as a function of Φ for a fixed field |B ext | = 533 G.We determined C from the Fourier space amplitude of individual FID curves, for varying field misalignment angles Φ x and Φ y , applied in the x-z and yz-planes, respectively.For the small angles we investigated, Φ ≈ (Φ 2 x + Φ 2 y ) 1/2 .As expected, when Φ = 0, no nuclear FID contrast is observed, because in this case B eff and e init are both collinear with e z such that the nuclear spin is optically pumped into the non-precessing eigenstate |0, +1/2 of Ĥgs .Upon increasing Φ, B eff and e init both tilt away from e z in different directions, resulting in nuclear precession of increasing contrast C. At the same time, increasing Φ (i.e.B ⊥ ) tends to reduce the nuclear hyperpolarization efficiency [28,29]  optical spin readout contrast [16], both of which reduce C. Overall, these counteracting effects imply that there is an optimal tilt angle Φ opt which maximizes C by balancing the magnitude of the nuclear spin coherences with nuclear spin readout efficiency.We call this maximized contrast C max .
To determine Φ opt and C max , we show in Fig. 2b the data from Fig. 2a as a function of transverse field B ⊥ , where for each data point, we determined B ⊥ from the NV's full optically detected magnetic resonance spectrum.Figure 2b reveals a clear maximum in C at B ⊥ ≈ 8.6 G, which corresponds to Φ opt ≈ 0.86 • (see further details in the supplementary information, Section III).These results are in good agreement with the predictions of our numerical model (black curve in in Fig. 2b).The quality of the simulation depends sensitively on the NV intersystem crossing rates, all of which were kept fixed to literature values [16] in our calculations (see further details in the supplementary information, Section II.3).We assign remaining discrepancies between data and simulations to uncertainties on optical transition rates employed in the model.
Interestingly, we find that our all-optical 15 N nuclear FID protocol is relatively resilient to deviations of B ext away from ideal ESLAC conditions.For this, we investigated the dependance of the contrast C on the applied magnetic field |B ext | and tilt angle Φ x , where for each data point, we ensured that Φ y = 0 is maintained to within experimental accuracy.The resulting data show a nontrivial dependance of C on |B ext | and Φ x (Fig. 2c), and in particular reveal that both Φ opt and C max change with |B ext | (Fig. 2d).These dependencies are qualitatively captured by our numerical model.We find a global maximum of C max ≈ 4.2 % for |B ext | = 533 G, and a drop of C over a full-width at half maximum (FWHM) range of ∼ 50 G.
Nuclear spin precession frequency.Further, we investigate the dependance of f nuc on the magnetic field tilt angles Φ x and Φ y .For this, we determine f nuc by Fourier analysis of the FID data for each data point sampled in Fig. 2a.The result is shown in Fig. 3a, along with the corresponding plot of the same data as a function of |B ⊥ | in Fig. 3b.The precession frequency f nuc increases with B ⊥ in a way that is excellently described by both Eq. ( 5) and our numerical model.We again assign small discrepancies between measured and predicted values of f nuc to experimental uncertainties in determining B ext .
Nuclear coherence time.The ease of use of our alloptical 15 N FID experiments enables facile assessment of the 15 N inhomogeneous nuclear spin coherence time T * 2 .To determine T * 2 , we extend the measurement pulse sequence shown in Fig. 1c to longer FID evolution times τ , to resolve the full decay of the signal (Fig. 4a).Fitting an exponentially decaying harmonic function to the data yields C = 3.09 ± 0.11 % and T * 2 = 248.1 ± 12.4 µs for the single NV center under investigation.This decoherence time is somewhat shorter than previously reported values [38], but consistent with the rather short NV electron spin relaxation time of T 1 = 315 ± 16 µs for our shallow NV -a timescale which is known to limit the NV's nuclear spin decoherence time [39,40].
Scaling to NV ensembles An interesting question with particular relevance for potential applications in quantum sensing is whether our all-optical scheme also scales to ensembles of NV centers.To address this question, we repeated our experiments on an ensemble of NV centers in an unstructured, CVD grown diamond sample with a preferential orientation of NV centers along one of the four possible crystal directions [41] (see methods).We maximize the contrast C in the same fashion as done in the single NV case, and determine T * 2 through the full FID trace shown in Fig. 4b.Using a least-squares fit as before, we find C = 2.08 ± 0.04%, and T * 2 = 508.5 ± 17.4 µs for this NV ensemble.While C is comparable to the single NV case (with a slight deterioration due to the minority of non-aligned NVs), T * 2 almost doubles.This value of T * 2 , however, is still short of the best reported values of ≈ 2.2 ms for 15 NV nuclear spin coherence times [26,38].We exclude electron spin T 1 relaxation as a source for this fast nuclear spin decoherence, as we measured T 1 = 5.8 ± 0.5 ms in this sample.Possible other sources for nuclear spin dephasing in our ensemble experiment include fluctuations in external magnetic field or temperature [38,40].

DISCUSSION
Our all-optical 15 N FID scheme lends itself to applications in quantum sensing, e.g. in magnetometry and gyroscopy (rotational sensing).In the following we discuss the predicted performance of such all-optical coherent quantum sensing schemes.
The shot noise limited FID sensitivity for spin based, low-frequency magnetometry is given by [17,42,43] Here, N is the average number of detected photons per readout pulse, C is the readout contrast, and γ is the gyromagnetic ration of the spins employed for sensing.Further sensitivity reductions due to overhead in preparation and measurement of quantum states are not included in this expression, but of little relevance to our conclusion, given the long T * 2 times at hand.Evaluating Eq. ( 6) for our single NV data (T * 2 = 250 µs, C = 4%, N = 0.1) using the effective nuclear gyromagnetic ratio γ = 1.2γI determined from the slope of the data in Fig. 1e (or from Eq. ( 5)), we obtain a photon shot noise limited magnetometry sensitivity of η mag = 154 µT/ √ Hz.Further, given that our approach scales to NV ensembles, we make predictions on future ensemble NV magnetometry sensitivity.For this, we assume a laser power of 100 mW, a 350 ns readout window, and a conversion ratio of excitation photons to detected PL photons of 3.4 % [44], to obtain N = 3.2 • 10 9 .Together with the measured ensemble values T * 2 = 500 µs and C = 2 %, we obtain η ensemble mag = 1.22 nT/ √ Hz.For spin based gyroscopy, the sensitivity is determined in analogy to Eq. ( 6), but with omission of the gyromagnetic ratio, i.e. η gyro = γ • η mag [28].Nuclear spins are therefore particularly attractive for gyroscopy, since their long T * 2 times generally offer them better sensitivities compared to electron spins, while they are less susceptible to magnetic fields and their fluctuations.Employing the same procedure as before, we obtain a projected ensemble gyroscope sensitivity of η ensemble gyro ≈ 135 • / √ hour.To place these estimates in context, we note that best reported magnetometry sensitivities using electron spin ensembles in diamond were ηensemble √ hour [26].While for magnetometry, our projected sensitivity alone is not competitive with the stateof-the-art, the microwave-free modality we present still lends itself to specific applications, e.g.remote sensing through optical fibres [46], or for cases where the MW drive would critically affect the sample of interest.Conversely, for gyroscopy, we project numbers competitive with previous approaches.The added feature of all-optical NV gyroscopy is hereby a key asset, which may enable future integrated and power-efficient NV gyroscopes.
Looking forward, we note that our all-optical nuclear spin sensing scheme is also amenable to alternate high fidelity readout schemes to increase measurement contrast, based on spin-to-charge conversion [47,48].Another potential path to improving contrast C is to dynamically pulse the field misalignment angle between spin initiali- sation and readout, to separately optimise the two processes.
In conclusion, we have presented an all-optical scheme for observing FID dynamics of 15 N nuclear spins in diamond NV centers.Our technique is based on optical pumping of the 15 N nuclear spin into a quantum superposition state -a novel optical pumping process that occurs near the NV's ESLAC in presence of a small transverse magnetic field.These results may find applications in various fields of quantum sensing, most notably alloptical magnetometry and gyroscopy, for which we give benchmark comparisons which compare favourably with the state-of-the-art.Our results also suggest possible extensions to a range of other, relevant scenarios, including analogous dynamics near the NV's ground-state level anti-crossing, or all-optical addressing of nearby 13 C nuclear spins.The nuclear spin dynamics we discussed should generally be observable in any color center exhibiting suitable level anti-crossing dynamics and coupling to nuclear spins, and might as such offer interesting opportunities for sensing with and characterization of novel color centers in a variety of solid state hosts.

Single NV diamond sample
The majority of our experimental results (Figs. [1][2][3] were obtained on a single NV center that was created in an "electronic grade" diamond sample (Element Six) by ion implantation and subsequent sample annealing [49].For implantation, we employed singly charged 15 N ions at a flux of 10 11 cm −2 and an energy of 6 keV, corresponding to a nominal implantation depth of about 9 nm [50].To increase PL collection efficiency, parabolic diamond pillars were fabricated in the diamond surface [51] subsequent to NV creation -a single pillar containing an individual NV center was studied in this work.
NV ensemble diamond sample The NV ensemble sample used to obtain data shown in Fig. 4b was grown on a CVD diamond substrate along the (113) crystal orientation to facilitate Nitrogen incorporation and create NVs preferentially oriented along the NV-axis lying closest to the growth plane [41].A 15 µm thick layer containing NVs was obtained using 12 C and 15 N enriched gas mixture, which led to an estimated NV density of ∼ 300 ppb [52].
Experimental setup A home-built confocal microscope (Olympus LMPLFLN-100 objective, NA = 0.8) was used to focus a green laser (Cobolt 06-MLD; emission wavelength 515 nm) on the sample and to simultaneously collect the emitted red PL.All data shown in this paper were taken by optically exciting the NV(s) near saturation, which for our setup corresponded to a laser power of about 50 µW for the single NV in a nano-pillar, and 2.2 mW for the ensemble of NVs in unstructured diamond.A static magnetic field was applied using a permanent neodymium disk magnet (supermagnete, S-60-05-N) mounted on a linear translation stage, to tune the magnetic field strength at the NV location.For precise magnetic field alignment near the ESLAC, the magnet is mounted on a goniometric stage (SmarAct SGO-60.5 and SGO-77.5).Finally, the laser and photon-detectors were gated with pulses which were created and synchronized using a high-frequency signal generator (Zurich Instruments SHFSG), which also served as a source for microwave pulses used for the characterization of the magnetic field via optically detected magnetic resonance experiment.
Here, latin indices denote states within a given subspace, and greek indicies count over the subspaces.Equation ( 7) is valid if the energy difference between states in different blocks is much larger than the coupling between them, Hamiltonian Ĥgs from the main text is such a block diagonal Hamiltonian.Written in the basis {|m S , m I }, it reads Ĥgs where the blocks are defined by states of equal values of m S , and where, without loss of generality, we define the direction of the transverse magnetic field as the x-direction, e.g. that B y = 0 and B x =: B ⊥ .We now evaluate Eq. ( 7) for Ĥgs = Ĥ0 + V with Ĥ0 = D gs 0 Ŝ2 z + γ S b z Ŝz + γ I b z Îz and V = B ⊥ (γ I Îx + γ S Ŝx ) + Ŝ • A • Î, and α corresponding to the m S = 0 subspace.This leads to the following matrix elements, Since γ I γ S , we simplify in all denominators the terms (γ S ± γ I ) ≈ γ S .In similar fashion, we use the fact that A gs D gs 0 and thus set (A gs ± 2D gs 0 ) ≈ ±2D gs 0 in all denominators.We obtain the following effective Hamiltonian for the 15 N spin in the m S = 0 manifold: Next, we add an energy-offset of D gs 0 (A gs ⊥ ) 2 − 2D gs 0 (γ S B ⊥ ) 2 / (D gs 0 ) 2 − (γ S B z ) 2 to place the energy levels symmetrically around zero, and thereby obtain the Hamiltonian given in the main text, where denotes the correction to the diagonal elements caused by mixing between states of different m S , and is the corresponding correction to the off-diagonal elements.Note that these expressions for ν z and ν ⊥ are diverging for D gs 0 = γ S B z , e.g.near the ground state level anti-crossing.However, Van Vleck formalism is not applicable to that regime since the corresponding electronic subspaces of H 0 are not sufficiently spaced in energy once this condition is approached.

II. Numerical Model for Optical Pumping under ESLAC Conditions
In this section, we present in detail our numerical model for simulating the spin dynamics of the NV center with and without green illumination for a given magnetic field B. To calculate the trajectory of both the electron and nuclear spin under optical pumping it is necessary to consider not only classical rate equations coupling the orbital states, but also to incorporate the quantum mechanical evolution of the spins within each orbital state.Our model follows an approach previously taken to model the effect of chemical reaction kinetics on NMR spectra [53].

II.1 Mathematical Description of the Model
We model the room temperature NV − center as a system made up of three distinct electronic states: The ground state (gs), the excited state (es), and a meta stable singlet state (s).We neglect the distinct E x and E y orbital branches in the excited state as they are efficiently averaged at room temperature, as well as the existence of two singlet states -we assume that there is only one such singlet state.We also neglect laser induced ionization to the NV 0 state.
First we define a separate spin density operator for each orbital state (labelled with α) which evolves coherently per the Liouville-von Neumann equation of motion where the carets denote operators, double carets denote superoperators, and α indexes over the different orbital states.The commutation superoperator Lα in Equation ( 16) can be calculated from the corresponding Hamiltonian Ĥα as where α ∈ {es, gs, s} denotes the orbital, E α is the identity matrix of the same dimensionality as H α and T denotes matrix transposition.The Hamiltonians for each orbital state are given by: Ĥgs,es /h = D gs,es 0 where Ŝ and Î are angular momentum operators for the electron and 15 N nucleus acting on the joint electron/nuclear Hilbert space, and σ are the 2x2 spin-1/2 matrices.The constants and coupling tensors are defined in table I. Next, we couple these differential equations with additional (real valued) superoperators corresponding to the incoherent optical pumping process.These superoperators act to reduce or increase the population of a given spin state and thus take the role of spin-selective relaxation superoperators where P±1,0 are projection superoperators that project ρα onto the NV-electron spin state with m S = ±1, 0 while leaving the dimensionality of ρα unchanged.This superoperator ensures that the rate of inter system crossing (ISC) out of the excited state depends on the instantaneous spin state population of ρes .Next, T e is a partial trace superoperator that acts on a 36-dimensional joint electron/nuclear density operator and traces out the NV-electron degrees of freedom, leaving a 4-dimensional density operator corresponding only to 15 N. Finally, Ŝ⊗ 0,±1 is a direct product superoperator that acts on a 4-dimensional 15 N density operator and turns it into a joint electron/nuclear density operator with the NV-electron in the state ms = 0, ±1.Note that because their effect changes dimensionality of ρα , the matrix representations of T e and Ŝ⊗ 0,±1 are not square.Further, since none of these superoperators act on the 15 N degrees of freedom this model assumes that the nuclear spin state is preserved throughout the optical pumping process.Equation ( 20) is normalized such that the sum of the traces of the spin density operators ρ α is equal to 1, meaning Trace{ρ α } is the fractional population of the total system in orbital state α.We will give procedures for generating the corresponding matrix representations for the superoperators in Eq. ( 20) later in this section.
Figure 5 shows a pictorial representation of the processes modeled by Eq. ( 20).The ground and excited states are coupled via optical excitation with rate k green and radiative decay with rate k red respectively.The spin selective ISC causes non-radiative transitions from the excited state to the singlet, described by the rates k ms=±1,0 ISC .Relaxation from the electronic singlet into the ground state is also electron spin selective, with rates given by k ms=1 and k ms=0 respectively.The precise values for the ISC and relaxation rates have been subject to considerable debate.At the end of these section we present modeling results for different sets of these parameters.
where E 36( 4) is a 36x36(4x4) identity matrix.In the following, we will call the resulting 76x76 dimensional matrix Â.Note that because of how we formed ρ by concatenation it cannot represent coherences between different orbital states, however such coherences are unlikely to be significant and are not necessary to explain the physics of interest in this work.We also note that since we use identity superoperators to represent the radiative processes, spin-spin coherences will be preserved under optical excitation and radiative decay within this model.Since the matrix Â commutes with itself for all values of time t, Eq. ( 21) can easily be integrated and thus the time evolution of the system can be calculated as ρ(t) = e 2π Ât ρ(0), (22) for any time t.At any given time the first 36 entries of ρ(t) corresponds to ρgs (t) written in vector form, the next 4 entries are ρs (t) and the last 36 are ρes (t).We take the predicted instantaneous photoluminescence for ρ(t) as the fractional population in the excited electronic state, We note that in this work we use Eq. ( 22) to calculate the time evolution of the system both in the presence and and absence of green illumination by choosing different values for the parameter k green and calculating the propagation piecewise with different Â matrices.

II.2 Matrix Representations of Superoperators
We define the matrix representations of the 15 N spin operators on the joint electron/nuclear space as Îx/y/z = σx/y/z ⊗ Ê3 and the NV-electron spin operators as Ŝx/y/z = Ê2 ⊗ λx/y/z , where λ are the 3x3 spin-1 matrices.The action of the partial trace superoperator is defined by T e Îx/y/z = Trace{ Ê3 }σ x/y/z .In the numerical implementation of our model we vectorize operators as 36x1 or 4x1 dimensional column vectors.The matrix representation of T e is determined following the procedure layed out in [53] and reproduced here for completeness: where α = (i − 1) where i, n = 1 ...d I count through the degrees of freedom of the first subspace (nuclear in our case), and j, m = 1 ...d S count through the degrees of freedom of the second subspace (NV-electron in our case).The resulting T e matrix is 4x36 dimensional.The matrix representation of the direct product superoperator Ŝ⊗ 0,±1 Ê2 = Ŝ0,±1 is where The numbers i, j, n, m, d S and d I are the same as defined above.The resulting dimensionality of the matrix representation of Ŝ⊗ 0,±1 is 36x4.
Finally, we define the matrix representation for the projection superoperator, P0,±1 = where ρ k 0,±1 = κk nuc ⊗ κ0,±1 is the 36x1 column vector representation of the joint 6x6 density matrix operator and While T e and Ŝ⊗ 0,±1 both change the dimensionality of the state they operate on, P0,±1 is a square 36x36 matrix and thus preserves dimensionality.The sum in equation (29) ensures that P0,±1 projects onto a particular electronic spin state, while leaving the nuclear spin unchanged.

II.3 Details on Simulation Evaluation
In order to simulate the PL of the all-optical pulse sequence introduced in the main text using the model described above, one needs to choose values for the various NV parameters and transition rates k.Here, apart from D gs 0 which we determine experimentally as described later in the SI, we use literature values for these simulation parameters.The exact numbers employed are listed in table I and originate from [33,34].For the optical transition rates, we consider five different sets of numbers, labeled parameter sets 1 to 5, as shown in the table II Note that parameters sets 1 and 2 assume k m S =0 ISC = 0, prohibiting inter system crossing for m S = 0 population in the excited state.Using the values of either of these sets, we fully simulate the all-optical sequence described in the main text, that consists of a green initialization pulse, followed by a interval of free evolution in the dark and subsequent PL readout during a short green laser pulse.First, we simulate optical pumping of ρ(0) at NV saturation for 3000 ns, e.g.setting k green = s • k red with s = 35% as determined experimentally on the single NV sample.To that end, we choose ρ(0) such that all six states in the ground electronic orbital are occupied with each 1/6, and the excited orbital states and the single states are empty.Afterwards, we run the model with k green = 0, propagating the optically pumped state in the dark for 290 ns.This empties out the singlet into the ground state.Subsequently, still with k green = 0, we propagate the system for a duration τ , followed by a readout made up of 350 ns with k green = s • k red .This readout pulse is simulated in steps of 1 ns, and the total PL(τ ) is the sum of the PL of each step.Plotting PL(τ ) yields an oscillation whose frequency f we determine by taking the position of maximum of the corresponding Fourier spectrum.The contrast C of the oscillation in PL(τ ) is also determined via the Fourier spectrum: We compute it as C = 4V /V 0 where V is the amplitude of the detected Fourier space peak, and V 0 is the Fourier space peak at zero.With this procedure, we simulate both the nuclear precession frequency f as well as the precession contrast C, for B-field values corresponding to every pixel shown in Fig 2b in the main text.The numeric result for the precession contrast C is shown in Fig. 6a, where for ease of comparison, the experimental data has been normalized to 1 and the simulation to 0.75.The plot demonstrates that the numeric result depends strongly on the chosen optical transition rates.Specifially, depending on the employed optical set of parameters, we find maximal contrast for B ⊥ between 7.9 and 10.2 G.
In Overall, we find that parameter set 4 has the best agreement with our data, and thus we use set 4 for all simulations shown in the main text's figures.Note that since PL(τ ) is not true to scale, we normalize all numeric results for C shown in the main text to the corresponding data.
Finally, we note that the numeric model agrees very well with the analytic predictions for the nuclear precession frequency.This excellent agreement is demonstrated in the main text's figure 3b.We want to emphasize, however, that not only the precession frequency f nuc = γ I |B eff | of numerical and analytical predictions agree well, but also the tilt angle ϑ of B eff , as shown in Fig. 7. Here, B eff and its tilt angle ϑ are calculated numerically by diagonalization of Ĥgs , or analytically via Eq.(2) in the main text.Combined, this means that both approaches yield the same B eff .
Simulation Analytic FIG. 7. Prediction of the tilt angle ϑ of B eff of both our the numeric simulation and the analytic approach based on Van Vleck perturbation theory.These predictions have excellent agreement.

III. Magnetic Field Map
In the main text's Fig. 2b (and 3b), we have shown the same data as in Fig 1a (and 3a), but plotted against the transverse magnetic field B ⊥ .Here, we elaborate on how this magnetic field is determined experimentally.Since we generate the external magnetic field with a permanent magnet that is mounted on a goniometer, and since this goniometer does not rotate about the NV perfectly, it is in general not accurate to simply measure |B ext | and set B ⊥ = sin(Φ)|B ext |.Instead, we measure the magnetic field vector components B ⊥ and B individually for every other pixel of the main text's Fig. 2a, using the same experimental conditions, by taking an optically detected magnetic resonance (ODMR) spectrum of each individual pixel.This way, we experimentally determine all four transition frequencies for a given magnetic field; two corresponding to the hyperfine transitions from the electron spin manifold |0 to the electron spin manifold |+1 , and two corresponding to the hyperfine transitions from |0 to |−1 .Then, we determine D gs 0 by fitting the ground state Hamiltonian Ĥgs as defined in the main text to the four measured transition frequencies of the center pixel, where we enforce B ⊥ = 0 since Φ x = Φ y = 0. We find D gs 0 = 2870.760402MHz.Using this calibration of D gs 0 , we then determine the magnetic components B (Φ x , Φ y ) and B ⊥ (Φ x , Φ y ) of all other pixels by fitting Ĥgs to the transition frequencies of these pixels, using D gs 0 as determined before.In the end, we interpolate both B (Φ x , Φ y ) and B ⊥ (Φ x , Φ y ) to twice the density in Φ x and Φ y to match the pixel density of Fig. 2a (respectively 3a).This interpolated result is shown in Fig. 8a  Further, these magnetic field data can be plotted against total angle, Φ = (Φ 2 x + Φ 2 y ) 1/2 , as shown in Fig. 8c, revealing a near-linear dependance of B ⊥ on Φ.This is to be expected, since for |B| = 533 G and small angles Φ as investigated here, The linear fit shown in Fig. 8c gives that in our experiment the slope is 10.01 G/deg.This slightly larger value compared to above prediction (and its not perfeclty linear shape) is explained by the fact that the goniometer employed in our experiments is not rotating the permanent magnet perfectly about the NV, resulting a slight change in total magnetic field whenever the magnet is rotated.
Finally, note that we use the slope of 10.01 G/deg to convert B ⊥ to Φ -but only for data taken at |B| = 533 G because this is the only field where we acquired ODMR data for the entire range of the goniometer.Specifically, as discussed in the main text, at |B| = 533 G the best contrast is achieved at B ⊥ = 9.4 G which thus corresponds to Φ opt = 0.94 • .

IV. Ensemble Data
In the main text's Fig. 2 and 3, we discuss the dependance of the of the observed nuclear precession on the applied magnetic field in the context of a single NV.However, all that data can also be reliably reproduced on a full ensemble of NVs.In Fig. 9, we show the results of these ensemble measurements which were taken on the NV ensemble described in the main text: an unstructured, grown diamond with 15 NVs that are preferentially aligned with a specific crystal orientation.These ensemble measurements reproduce the single-NV data in both qualitative and quantitative way, demonstrating that the all-optical nuclear quantum sensing scheme proposed in this work does scale well to NV ensembles.We note, however, that the ensemble shows a contrast C of about half the single-NV contrast, which we assign to the minority of non-aligned NVs.

FIG. 1 .
FIG.1.a Crystal structure of the Nitrogen-Vacancy center with illustration of its associated spins and coordinate axes.b All-optical nuclear spin precession of the 15 N Nuclear spin observed at a magnetic field of |B| = 540 G tilted away from the NV symmetry axis by Φ = 1 • .Fitting of a harmonic function (black) yields a precession frequency 251.18 ± 0.12 kHz.c Pulse sequence employed for b, consisting of a 3 µs green laser pulse separated by a variable delay τ .The first 350 ns of each green pulse are utilized for optical nuclear spin readout, while the remainder of the pulse reinitializes the spin system.d Quantitative Bloch-Sphere representation of the15 N spin in the |ms = 0 ground state manifold.For a magnetic field Bext tilted from the NV symmetry axis by the angle Φ = 1 • , optical pumping initializes the nuclear spin into ρnuc init .The nuclear spin subsequently precesses around an effective magnetic field B eff .The measurement axis for all-optical readout of this precession is given by einit.e Experimentally observed precession frequency (blue crosses) at different values of |Bext| and fixed Φ = 1 • , together with analytic (solid orange line) and numerical predictions (black dots).

FIG. 2 .
FIG. 2. a Measured free induction decay (FID) contrast of a single NV as a function of magnetic field orientation, for a fixed total field of |Bext| = 533 G. b Same data as in a, but plotted against total transverse magnetic field B ⊥ .The black line is the prediction of our numerical model, which is normalized to the mean of the data at Φopt ± 0.03 • .c Nuclear FID contrast as a function of magnetic tilt angle Φ and total magnetic field |Bext|.d Maximal observed contrast Cmax and corresponding tilt angle Φopt for each value of |Bext|.The prediction of our numerical model is shown in black, which is normalized to the maximal data point.

FIG. 3 .
FIG.3.a15N nuclear precession frequency fnuc corresponding to the data shown in Fig.2a.Only pixels for which C > 1% are shown.b Same data as in a, plotted against total transverse magnetic field B ⊥ , together with the numerical model prediction (black), the prediction from Eq. (5) (orange), and the bare Larmor frequency γI Bz (green).

FIG. 4 .
FIG. 4. All-Optical nuclear spin precession of a, a single NV in a diamond nanopillar, and b, an ensemble of NVs in bulk diamond, both measured at |Bext| = 533 G with tilt angles Φ = 0.65 • and Φ = 0.8 • , for a and b respectively.Each data set is fitted with a damped harmonic function to determine the nuclear spin coherence time T * 2 , yielding T * 2 = 248.1 ± 12.4 µs and T * 2 = 508.5 ± 17.4 µs, for a and b respectively.

FIG. 5 .
FIG. 5. Level Structure for the optical pumping model.Each of the three orbital states is governed by its own Hamiltonian.

FIG. 6 .
FIG.6.The data shown in Fig2b and din the main text, together with the results of the numerical simulation, run with each of the five optical transition rate parameter sets.Overall, set 4 has the best agreement with the data, which is why model 4 is used for all numerical results shown in the main text.
Fig 6b we show the results of running the model for B-field values that correspond to the measurement shown in Fig 2d in the main text.Again, it is evident that the choice of optical transition rates has a significant effect on the numeric result.All set of parameters predict a global maximum in C between 530 and 550 G. Interestingly, sets 1 and 2 which dictate k m S =0 ISC = 0 differ drastically from the other parameter sets with k m S =0 ISC = 0 in that they show a second local maximum for C at about 470 G.This is not consistent with our data.

FIG. 8 .
FIG. 8. a Parallel magnetic field component B and b transverse magnetic field component B ⊥ as a function of tilt angles Φx and Φy, measured via optically detected magnetic resonance of all microwave-driving 15 N transitions.Each of the shown pixels corresponds to one pixel in the main text's Fig. 2a and 3a.The total magnetic field is about |B| = 533 G. c Measured transverse magnetic field as a function of total tilt angle Φ, revealing a linear dependance with slope 10.01 G/deg, allowing for a simple conversion of B ⊥ to corresponding Φ and vice versa.

FIG. 9 .
FIG.9.Same measurement as shown in the main text in Fig.2a, c, dand Fig.3a, but here the data were taken on an ensemble of NVs rather than on a single NV.We find the same qualitative and quantitative result, with the only difference being a slightly lower contrast C.
nuc is the expected nuclear spin precession frequency.Importantly, B eff is neither aligned with B ext nor with the NV symmetry axis.Instead, it lies in the plane spanned by B eff and B ext , and is tilted away from e z by an angle ϑ

TABLE I .
. Numeric values for the employed NV parameters.

TABLE II .
Numeric values for of the optical transition rates used in our model.We consider five different parameter sets.In the main text, we use set 4.