Detection of nanoscale electron spin resonance spectra demonstrated using nitrogen-vacancy centre probes in diamond

Electron spin resonance (ESR) describes a suite of techniques for characterizing electronic systems with applications in physics, chemistry, and biology. However, the requirement for large electron spin ensembles in conventional ESR techniques limits their spatial resolution. Here we present a method for measuring ESR spectra of nanoscale electronic environments by measuring the longitudinal relaxation time of a single-spin probe as it is systematically tuned into resonance with the target electronic system. As a proof of concept, we extracted the spectral distribution for the P1 electronic spin bath in diamond by using an ensemble of nitrogen-vacancy centres, and demonstrated excellent agreement with theoretical expectations. As the response of each nitrogen-vacancy spin in this experiment is dominated by a single P1 spin at a mean distance of 2.7 nm, the application of this technique to the single nitrogen-vacancy case will enable nanoscale ESR spectroscopy of atomic and molecular spin systems.

Here we develop the theory of relaxation based sensing as used in the main text of this work. As we are considering axial magnetic field strengths, B 0 , such that B 0 ∼ 2πD/2γ e ∼ 512 G (where D = 2.87 GHz is the zero-field splitting of the NV spin, and γ e = 17.6 × 10 6 s −1 G −1 is the electronic gyromagnetic ratio), only the |0⟩ ↔ | − 1⟩ transitions of the NV spin will be appreciably excited, meaning we can disregard any population of the | + 1⟩ state. The time evolution of the associated density matrix, ρ T , is described by where ρ T represents the combined density matrix of the entire NV spin + environment system. The full Hamiltonian is given by H T = H NV + H int + H E where H NV and H E are the self Hamiltonians of the NV centre and environment, which, for a general spin bath environment are given by respectively, where we have assumed that the both the z axis and the external magnetic field are aligned with the principle axis of the NV spin, and E ij is the tensor describing the interaction between spins i and j in the environment, which in general may include both exchange and magnetic-dipole interactions depending on the environment in question. Due to the highly localised nature of the NV wavefunctions, the coupling of the environmental spins to the NV may be described by the magnetic dipolar interaction alone, where B i is the symmetric magnetic dipole tensor describing the interaction of the NV spin with the i th environmental spin, and includes both transverse and longitudinal components, proportional to P x,y and P z of the NV spin respectively. The latter have a pure dephasing effect, resulting in an additional contribution to the intrinsic dephasing rate of the NV. As relaxation processes occur on timescales that are much longer that the typical interaction timescales of the environmental constituents, essentially placing the system in a Markovian regime, the resulting dephasing will be purely exponential. These effects may thus be modeled using a master equation approach for the reduced density matrix of the NV spin, ρ NV , as follows where, in the present context, L is the Lindbladian operator corresponding to a pure dephasing process on the NV spin, and is given by L = √ 2Γ NV P z . The total dephasing rate due to both the local crystal environment and the longitudinal coupling to any external environment is given by . The timescale of the intrinsic dephasing process is described using the inhomogeneous linewidth, (T * 2 ) −1 , since the transverse phase accumulation occurs in the absence of any pulsed microwave control. Subtle tuning effects that modify the sensitivity of this technique to various parts of the environmental spectral density may be achieved by changing the intrinsic dephasing rate via dynamic decoupling techniques.
In what follows, owing to the strong intra-environment and comparatively weak NV-environment couplings, we will treat the coupling of the environment to the transverse components of the NV spin as a semiclassical oscillatory field (these simplifications will be later justified in Supplementary Note 2), where B (ω E ) = b x (ω E ) P x +b y (ω E ) P y ; and b x and b y are the x and y components of the magnetic field. The frequency spectrum (ie the distribution of ω E ) is determined by analysing the interaction between environmental constituents, as described by H E (again, see Supplementary Note 2). To make the solution tractable, we change to the interaction picture. The transformed equation of motion is given by with the interaction Hamiltionian given by V I = e iHNVt Ve −iHNVt . We are then interested in determining the rate at which the NV spin relaxes to its equilibrium state under the influence of the environment. We proceed by reducing the 3 × 3 system of first order linear differential equations described by equation 6 to a higher order differential equation for P 0 ≡ ρ 00 . We then wish to solve this equation, together with the initial conditions of ρ ij = 0 unless i = j = 0, in which case we have ρ 00 = 1, representing the initial polarisation of the NV spin in the |0⟩ state.
To gain insight into the expected analytic solution for the spin-1 NV centre, we consider the simplified case in which only one of the transitions of the NV centre is excited by the environment, and the other is assumed to be too far detuned to have any effect on the population of the spin states. This simplifies the analysis dramatically, yet demonstrates the main properties of relaxation based detection, and is applicable to a spin-1 system for cases of significantly strong Zeeman splittings between the | ± 1⟩ states.
The equation of motion for where b ≡ ⟨ B 2 ⟩ 1 2 is the second moment of the strength of the effective magnetic field operator, B, and δ = 2πD − ω NV − ω E is the detuning between the energy levels of the NV (2πD − ω NV ) and environmental (ω E ) frequencies.
A. Response to a monochromatic transverse field

Resonant case
For the case where the frequency of the environment is resonant with the transition frequency between the probe's spin states (δ = 0), the solution of Supplementary Equation 7 is Typically the spin based environments in which we are interested couple weakly to the NV spin as compared with its intrinsic dephasing rate, implying Γ NV ≫ b. In fact, even a strong coupling will also induce additional dephasing, so even in a worst case scenario we are guaranteed Γ NV > b. In this limit, we have Hence the resonant (and therefore maximal) longitudinal relaxation rate is given by

General case
When a finite detuning, δ, exists, we may examine relative importance of the terms within Supplementary Equation 7, subject to rescaling t in terms of the decay time from the resonant solution. That is, if we consider the dimensionless variable T = Γ max 1 t and retain terms up to and including order O , the solution for an arbitrary detuning becomes For zero detuning, we recover the previous result (Supplementary Equation 10). For finite detuning, the relaxation rate is modified by a Lorentzian factor with a FWHM of Γ NV . The complete decay profile is then obtained by integrating this expression over the spectral density of the environment, implying that the δ-dependent relaxation rate acts acts to filter out the environmental spectrum about δ ∼ 0.

B. Response to a transverse field with an arbitrary spectral density
Even without considering the specifics of the spectral density, the response of the NV spin to an arbitrary spin bath can vary remarkably due the geometric proximity and arrangement of the bath relative to the NV centre. The definitions of the NV spin relaxation and the corresponding filter function are given by ,and respectively, where the filter function, G(ω NV , ω E ) acts to filter out regions of the spectral density (as dependent on the external field strength, B 0 = ω NV /γ e ), and depends explicitly on the geometric arrangement of the environmental constituents. Ultimately, given some measurement record and filter function, it is expression 12 that must be deconvolved to reproduce the spectral density, S(ω E ). In this section, we consider the effects of the geometric arrangement of the environment on the filter function, G, for a general spectral density; and specific case of the internal P1 nitrogen donor electron spin bath in type-1b diamond is considered below in Supplementary Note 2.

Response to an internal (bulk) spin bath
We note that the coupling of the NV to a bath spin located at some distance r may be written b ≡ β(θ, ϕ)/r 3 , where the specific details of the angular dependence of the coupling are incorporated into the parameter β(θ, ϕ). We note that is is necessary to omit further discussion of β until Supplementary Note 2, as different environmental processes will be more readily detectable by the NV spin at different relative angles. Unlike the transverse (spin echo) case, where the precession of the NV spin vector in the x − y plane is sensitive to all longitudinal field sources, the effect on the longitudinal projection is dominated by the coupling to the nearest P1 electron spin. As such, we may write P b (b) = P r (r)P θ,ϕ (θ, ϕ), where r, θ, ϕ are the spherical coordinates associated with the distribution of field sources.
From Ref. 1, the distribution of distances from a given NV centre to its nearest spin impurity is given by where n is the average density of impurities in the bath. Substituting this expression into Supplementary Equation 12, we find where G is the Meijer-G function, and ⟨ |β| ⟩ = ∫ |β(θ, ϕ)| P θ,ϕ (θ, ϕ) dθ dϕ. Thus, we identify the filter function associated with environments inside the diamond lattice to be where A in is a constant, associated with the geometry of the bath, that may be renormalised.

Response to an external (surface) spin bath
In contrast to the bulk spin bath case, spins on the surface are unable to exist arbitrarily close to the NV centre. Typically we consider samples in which NV centres exist at some depth h + δh below the surface, with h being the mean depth, and δh a normally distributed variable with variance ⟨ δh 2 ⟩ ≪ h 2 . In this case, an individual NV spin is exposed to many bath spins, meaning that the effective coupling distribution, P b (b), is normally distributed.
In this case, we may expand Supplementary Equation 11 for small t, which, upon substitution into Supplementary Equation 12 and averaging over P b (b) gives Thus, we identify the filter function associated with environments outside the diamond lattice to be Supplementary Note 2.

Theoretical description of the coupled NV-P1 system
In this section, we discuss the features we expect to be evident in the P1 centre spectrum by examining the effect of a P1 centre on the magnetic field dependent relaxation rate of a near-by NV centre. We conclude this section by demonstrating the equivalence of the semi-classical approach used in this work, and a quantum mechanical treatment of the NV-P1 interaction.

A. P1 Hamiltonian
The Hamiltonian of a P1 centre is given by where A P1 is the hyperfine tensor describing the coupling between the P1 electron, ⃗ S, and 14 N nuclear spin, ⃗ I; and Q is the quadrupole splitting of the nuclear spin. For all field strengths at which there is an appreciable overlap of this spectrum with the NV spin filter function, B 0 ∼ 512 G, we find that the eigenstates of the P1 centre electron spin are predominantly dictated by the external magnetic field. In this instance, the Hamiltonian of the P1 centre becomes for cases where the P1 axis is aligned with the external magnetic field, where A z = 114 MHz, and A x = 81.3 MHz. If the P1 axis is aligned along one of the three other bond directions not aligned with the field, the Hamiltonian may be transformed via the rotation operator R = exp (−iI y θ) = 1 − iI y sin(θ) − I 2 y (1 − cos(θ)) (The other two axes may be realised via a trivial rotation about the z axis), and is given by where a z = 1 9 (8A x + A z ) = 85 MHz, and a x = 1 9 (5A x + 4A z ) = 99 MHz.

B. Coupling of the P1 environment to the NV centre
The interaction between the NV spin and the P1 nuclear spin is ignored on account of its comparative weakness. For the interaction between the NV and the P1 electron, the magnetic dipole interaction is given by wherer is the unit separation vector between the NV and P1 centres, and is the effective dipolar coupling strength.
The components of the NV-P1 interaction responsible for the relaxation of the NV spin are those coupling to its transverse components, namely P x and P y . Without loss of generality, we may rewrite this interaction as an effective quantum mechanical magnetic field, ⃗ B, where B x,y,z couples to P x,y,z . To make use of Supplementary Equation12, we make the semi-classical approximation, and assume that NV-P1 interaction plays very little part in determining the dynamics of the environment. The problem of determining the environmental spectrum then reduces to solving for the environmental evolution exclusively under its own influence as follows.
In order to determine the effective strength of the semi-classical magnetic field used to model the effect of the P1 bath on the NV spin, we must determine which components of the NV-P1 interaction, H int , are relevant, and which may be discarded for a given magnetic field strength. To do this, we transform H int to the interaction picture (we point out that only the Hilbert space associated with the NV and P1 electron spins need be considered, as the interaction between the NV electron spin and the P1 nuclear spin is sufficiently weak that it may be ignored). In matrix form, H int is given (in the In order to simplify H int , we switch to the interaction picture to see which terms are important near 512 G.
where k = −1, 0, +1 is the hyperfine projection of the nuclear spin. Using the rotating wave approximation, we can see that only terms of frequency 2πD −2ω 0 −2πkA z need be retained, and all other off-diagonal terms may be ignored. Transforming back to the Schrodinger picture, the simplified interaction Hamiltonian becomes from which we infer that the effective transverse field strength associated with the allowed NV-P1 |0, ↓⟩ ↔ | + 1, ↑⟩ transitions ((a), (b) in Supplementary Figure 1) is given by Determination of the effective field strength associated with the disallowed transitions ((c), (d) in Supplementary  Figure 1) follows the same approach as above. From this, we find the associated field strength to be

C. P1 dynamics
To determine the dynamic behaviour of the P1 environment, we compute the autocorrelation functions associated with the field components above. Interactions between environmental components may be modeled by damping these autocorrelation functions with a decaying exponential, exp (−Γ P1 t) to describe their relaxation due to mutual flip-flop processes with corresponding relaxation rate Γ P1 .
The corresponding spectra may then be found by computing the Fourier transforms of the autocorrelation functions. From this, we find the spectra associated with the allowed transitions to be S on all (ω) = 1 6π and S off all (ω) = 1 6π for cases of on and off axis P1 centres respectively. Taking the relative proportions of on and off-axis P1 centres to be 25% and 75% respectively, we find the overall spectrum associated with the allowed transitions to be Similarly, the spectra associated with the disallowed transitions are given by and respectively, where The spectrum associated with the disallowed transitions is then S dis (ω) = 1 4 S on dis (ω) + 3 4 S off dis (ω).
By employing the full spectrum, S(ω E ) = S all (ω E ) + S dis (ω E ), in equation 12, we find the resulting external fielddependent relaxation rate of the NV centre to be where the effective couplings, , are due to integration over all possible NV-P1 separations. By taking the average P1 density to be 50 ppm, and the FID rate to be Γ 2 = 5.0 MHz, we may plot the resulting field-dependent relaxation rate of the NV spin (Supplementary Figure 1).