Three-dimensional localization spectroscopy of individual nuclear spins with sub-Angstrom resolution

Nuclear magnetic resonance (NMR) spectroscopy is a powerful method for analyzing the chemical composition and molecular structure of materials. At the nanometer scale, NMR has the prospect of mapping the atomic-scale structure of individual molecules, provided a method that can sensitively detect single nuclei and measure inter-atomic distances. Here, we report on precise localization spectroscopy experiments of individual 13C nuclear spins near the central electronic sensor spin of a nitrogen-vacancy (NV) center in a diamond chip. By detecting the nuclear free precession signals in rapidly switchable external magnetic fields, we retrieve the three-dimensional spatial coordinates of the nuclear spins with sub-Angstrom resolution and for distances beyond 10 Å. We further show that the Fermi contact contribution can be constrained by measuring the nuclear g-factor enhancement. The presented method will be useful for mapping atomic positions in single molecules, an ambitious yet important goal of nanoscale nuclear magnetic resonance spectroscopy.

O ne of the visionary goals of nanoscale quantum metrology with NV centers is the structural imaging of individual molecules, for example proteins, that are attached to the surface of a diamond chip 1 . By adapting and extending measurement techniques from nuclear magnetic resonance (NMR) spectroscopy, the long-term perspective is to reconstruct the chemical species and three-dimensional location of the constituent atoms with sub-Angstrom resolution 2,3 . In contrast to established structural imaging techniques like X-ray crystallography, cryo-electron tomography, or conventional NMR, which average over large numbers of target molecules, only a single copy of a molecule is required. Conformational differences between individual molecules could thus be directly obtained, possibly bringing new insights about their structure and function.
In recent years, first experiments that address the spatial mapping of nuclear and electron spins with NV-based quantum sensors have been devised. One possibility is to map the position into a spectrum, as it is done in magnetic resonance imaging. For nanometer-scale imaging, this requires introducing a nanomagnet [4][5][6] . Another approach is to exploit the magnetic gradient of the NV center's electron spin itself, whose dipole field shifts the resonances of nearby nuclear spins as a function of distance and internuclear angle. Refinements in quantum spectroscopic techniques have allowed the detection of up to 8 individual nuclear spins 7,8 as well as of spin pairs [9][10][11] for distances of up to~30 Å 12,13 . Due to the azimuthal symmetry of the dipolar interaction, however, these measurements can only reveal the radial distance r and polar angle θ of the inter-spin vector r = (r, θ, ϕ), but are unable to provide the azimuth ϕ required for reconstructing three-dimensional nuclear coordinates. One possibility for retrieving ϕ is to change the direction of the static external field 12 , however, this method leads to a mixing of the NV center's spin levels which suppresses the ODMR signal 14 and shortens the coherence time 15 . Other proposed methods include positiondependent polarization transfer 16 or combinations of microwave and radio-frequency fields [17][18][19] .
Here we demonstrate three-dimensional localization of individual, distant nuclear spins with sub-Angstrom resolution. To retrieve the missing angle ϕ, we combine a dynamic tilt of the quantization axes using a high-bandwidth microcoil with high resolution correlation spectroscopy 20,21 . Our method provides the advantage that manipulation and optical readout of the electronic spin can be carried out in an aligned external bias field. This ensures best performance of the optical readout and the highest magnetic field sensitivity and spectral resolution of the sensor.

Results
Theory. We consider a nuclear spin I = 1/2 located in the vicinity of a central electronic spin S = 1 with two isolated spin projections m S = {0, −1}. The nuclear spin experiences two types of magnetic field, a homogeneous external bias field B 0 (aligned with the quantization axis e z of the electronic spin), and the local dipole field of the electronic spin. Because, the electronic spin precesses at a much higher frequency than the nuclear spin, the latter only feels the static component of the electronic field, and we can use the secular approximation to obtain the nuclear free precession frequencies, Here, γ n is the nuclear gyromagnetic ratio and is the secular hyperfine vector of the hyperfine tensor AðrÞ that gives rise to the hyperfine magnetic field m S A z (r)/γ n (see Fig. 1b).
To obtain information about the distance vector r, a standard approach is to measure the parallel and transverse components of the hyperfine vector, a || = A zz and a ? ¼ , and to relate them to the field of a point dipole, where μ 0 = 4π × 10 −7 TmA −1 is the vacuum permeability, ħ = 1.054 × 10 −34 Js is the reduced Planck constant, |γ e | = 2π × 28 GHz T −1 is the electron gyromagnetic ratio, and where we have included a Fermi contact term a iso (set to zero for now) for later discussion. Experimentally, the parallel projection a || can be inferred from the precession frequencies f m S using Eq. (1), and the transverse projection a ⊥ can be determined by driving a nuclear Rabi rotation via the hyperfine field of the central spin and measuring the rotation frequency 21 . Once a || and a ⊥ are known, Eqs. (3) and (4) can be used to extract the distance r and polar angle θ of the distance vector r = (r, θ, ϕ). Due to the rotational symmetry of the hyperfine interaction, however, knowledge of a || and a ⊥ is insufficient for determining the azimuth ϕ.
To break the rotational symmetry and recover ϕ, we apply a small transverse magnetic field ΔB during the free precession of the nuclear spin. Application of a transverse field tilts the quantization axes of the nuclear and electronic spins. The tilting modifies the hyperfine coupling parameters a || and a ⊥ depending on the angle between ΔB and A z , which in turn shifts the nuclear precession frequencies f m S . To second-order in perturbation theory, the m S -dependent precession frequencies are given by: 22 where α(m S ) is a small enhancement of the nuclear g-factor. The enhancement results from non-secular terms in the Hamiltonian that arise due to the tilting of the electronic quantization axis, and is given by 22 Here D = 2π × 2.87 GHz is the ground-state zero-field splitting of the NV center. By measuring the shifted frequencies f m S and comparing them to the theoretical model of Eqs. (5) and (6), we can then determine the relative ϕ angle between the hyperfine vector and ΔB.
Experimental setup. We experimentally demonstrate threedimensional localization spectroscopy of four 13 C 1-4 nuclei adjacent to three distinct NV centers. NV 1 is coupled to two 13 C spins, while NV 2 and NV 3 are each coupled to a single 13 C spin. For readout and control of the NV center spin, we use a custombuilt confocal microscope that includes a coplanar waveguide and a cylindrical permanent magnet for providing an external bias field of B 0~1 0 mT applied along the NV center axis e z . Precise alignment of the bias field is crucial for our experiments and is better than 0.3°(Methods section).
To dynamically tilt the external field we implement a multiturn solenoid above the diamond surface (Fig. 1d). The coil produces~2.5 mT field for 600 mA of applied current and has a rise time of~2 μs. We calibrate the vector magnetic field of the coil with an absolute uncertainty of less than 15 μT in all three spatial components using two other nearby NV centers with different crystallographic orientations (ref. 23 and Methods section).
Mapping of r and θ. We begin our 3D mapping procedure by measuring the parallel and perpendicular hyperfine coupling constants using conventional correlation spectroscopy 21 with no coil field applied, ΔB = 0 (Fig. 2). The parallel coupling a || is determined from a free precession experiment (sequence 1 in Fig. 2) yielding the frequencies f 0 and f −1 (Fig. 2b). The coupling constant is then approximately given by a || /(2π) ≈ f −1 − f 0 . The transverse coupling a ⊥ is obtained by driving a nuclear Rabi oscillation via the NV spin, using sequence 2, and recording the oscillation frequency f R , where a ⊥ /(2π) ≈ πf R (Fig. 2c). Because, the Zeeman and hyperfine couplings are of similar magnitude, these relations are not exact and proper transformation must be applied to retrieve the exact coupling constants a || and a ⊥ (ref. 21 and Methods section). Once the hyperfine parameters are known, we can calculate the radial distance r = 8.58(1) Å and the polar angle θ = 52.8(1)°of the nuclear spin by inverting the point-dipole formulas (Eqs. 3,4). The measurement uncertainties in r and θ are very small because correlation spectroscopy provides high precision estimates of both a || and a ⊥ .
Mapping of ϕ. In a second step, we repeat the free precession measurement with the coil field turned on (sequence 3), yielding a new pair of frequency values f ′ 0 , f ′ À1 (Fig. 2d). We then retrieve ϕ by computing theoretical values for f ðthÞ 0 , f ðthÞ À1 based on Eq. (5) and the calibrated fields in Table 1, and minimizing the cost function with respect to ϕ. To cancel residual shifts in the static magnetic field and improve the precision of the estimates, we compare the frequency difference between m S states rather than the absolute precession frequencies.
In Fig. 3a, we plot |ξ(ϕ)| for three different coil positions and opposite coil currents for 13 C 1 . We use several coil positions because a single measurement has two symmetric solutions for ϕ, and also because several measurements improve the overall accuracy of the method. The best estimate ϕ = 239(2)°is then given by the least squares minimum of the cost functions (dashdotted line in Fig. 3a). To obtain a confidence interval for ϕ, we calculate a statistical uncertainty for each measurement by Monte Carlo error propagation taking the calibration uncertainties in B 0 and ΔB, as well as the measurement uncertainties in the observed precession frequencies into account (Methods section). Values for all investigated 13 C nuclei are collected in Supplementary Tables 2-9.
Fermi contact interaction. Thus far we have assumed that the central electronic spin generates the field of a perfect point dipole. Previous experimental work 22,24 and density functional theory (DFT) simulations 25,26 , however, suggest that the electronic wave function extends several Angstrom into the diamond host lattice. Fig. 1 Coordinate systems for spins and magnetic fields. a Reference frame of the central nitrogen-vacancy (NV) sensor spin (red) with a nuclear spin (blue) located at the three-dimensional position r = (r, θ, ϕ). The quantization axis of the NV center defines the z-axis. The hyperfine field of the NV spin (red field lines) provides the magnetic field gradient for imaging. b Sketch of two nuclear spins I 1 and I 2 experiencing the same hyperfine interaction (red) [Eq. (2)]. Application of a transverse field ΔB (purple) reduces (I 1 ) or increases (I 2 ) the total magnetic field B ′ tot (blue) experienced by the nuclear spins depending on the ϕ angle, allowing us to discriminate the nuclear locations. B 0 is the static external field (green). c Geometry of the experimental setup in the laboratory frame of reference. A small solenoid on top of the diamond chip provides a rapidly switchable magnetic field ΔB. To change the vector orientation of ΔB, we translate the coil over the diamond The finite extent of the spin density leads to two deviations from the point dipole model: modified hyperfine coupling constants A ij , and a non-zero Fermi contact term a iso . In the remainder of this study we estimate the systematic uncertainty to the localization of the nuclear spins due to deviations from the point dipole model.
We first consider the influence of the Fermi contact interaction, which arises from a non-vanishing NV spin density at the location of the nuclear spin. The Fermi contact interaction adds an isotropic term to the hyperfine coupling tensor, A + a iso 1, which modifies the diagonal elements A xx , A yy , and A zz . DFT simulations 25,26 indicate that a iso can exceed 100 kHz even for nuclear spins beyond 7 Å. It is therefore important to experimentally constrain the size of a iso .
To determine a iso , one might consider measuring the contact contribution to the parallel hyperfine parameter a || , which is equal to A zz . This approach, however, fails because a measurement of a || cannot distinguish between dipolar and contact contributions. Instead, we here exploit the fact that the gyromagnetic ratio enhancement α depends on A xx and A yy , and hence a iso . To quantify the Fermi contact coupling we include a iso as an additional free parameter in the cost function (7). By minimizing ξ(ϕ, a iso ) as a joint function of ϕ and a iso and generating a scatter density using Monte Carlo error propagation, we obtain maximum likelihood estimates and confidence intervals for both parameters (Fig. 3b). The resulting contact coupling and azimuth for nuclear spin 13 C 1 are a iso /(2π) = 9(8) kHz and ϕ = 238(2)°,  Supplementary Tables 2-9. b Free precession signal of the nuclear spin as a function of time t, using sequence 1. Right panel shows the corresponding power spectrum. The two frequencies f 0 and f −1 are approximately equal to γ n B 0 /(2π) and (γ n B 0 + a || )/(2π), respectively, see text. c Application of periodic π pulses on the NV center during t (sequence 2) causes a Rabi nutation of the nuclear spin, whose oscillation frequency f R is approximately equal to (a ⊥ /π)/(2π). d Activation of a transverse microcoil field ΔB during the nuclear precession (sequence 3) leads to shifted frequencies f ′ 0 and f ′ À1 . All measurements were conducted on 13 C 1 . Extracted frequencies are listed in Table 1 Table 1 Data base of measured precession frequencies and calibrated external magnetic fields used to determine the 3D position of 13 C 1

Quantity Value
Reference   Table 2 and 3 Five further measurements of f ′ 0 ; f ′ À1 À Á were made to improve the localization accuracy (data given in Supplementary Tables 2, 3). Vector magnetic fields refer to the NV coordinate system defined in Fig. 1 respectively; data for 13 C 2-4 are collected in Table 2. Because the gyromagnetic ratio enhancement α is only a second-order effect, our estimate is poor, but it still allows us constraining the size of a iso . By subtracting the Fermi contact contribution from a || , we further obtain refined values for the radial distance and polar angle, r = 8.3(2) Å and θ = 58(4)°. Note that introducing a iso as a free parameter increases the uncertainties in the refined r and θ, because the error in a iso is large. This leads to disproportionate errors for distant nuclei where a iso is small. Once nuclei are beyond a certain threshold distance, which we set to r = 10 Å in Table 2, it therefore becomes more accurate to constrain a iso = 0 and apply the simple point dipole model.
Extended electronic wave function. The second systematic error in the position estimate results from the finite size of the NV center's electronic wave function. Once the extent of the wave function becomes comparable to r, the anisotropic hyperfine coupling constants A ij are no longer described by a point dipole, but require integrating a geometric factor over the sensor spin density 25 . While we cannot capture this effect experimentally, we can estimate the localization uncertainty from DFT simulations of the NV electron spin density. Following ref. 26 , we convert the calculated DFT hyperfine parameters of 510 individual lattice sites to (r, θ) positions using the point-dipole formula (Eqs. 3, 4), and compute the difference to the DFT input parameters (r DFT , θ DFT ). The result is plotted in Fig. 4a. We find that the difference 〈Δr〉 = r − r DFT decreases roughly exponentially with distance, and falls below 0.2 Å when r > 10 Å (gray dots and curve). Figure 4b summarizes our study by plotting the reconstructed locations for all four carbon atoms in a combined 3D chart. The shaded regions represent the confidence areas of the localization, according to Table 2, projected onto the Cartesian coordinate planes. We note that the DFT simulations are in good agreement with our experimental results.

Discussion
The accuracy of our present experiments is limited by deviations from the point-dipole model, which dominate for small r (Fig. 4a). For larger r ≳ 1 nm, this systematic uncertainty becomes negligible, and the localization imprecision is eventually dictated by the NMR frequency measurement. In the present study, which probed isolated 13 C nuclei with a narrow intrinsic line-width, the frequency precision was limited by the accuracy of our detection protocol to~100 Hz. This corresponds to a radial localization error of~0.75 Å at a distance of r = 3 nm (solid black line in Fig. 4a and Supplementary Fig. 6). Improving the frequency precision to 3 Hz 11,27 would extend this distance to r~7 nm (dashed black line).
Our work demonstrates a basic strategy for mapping spatial positions of single nuclei in 3D with high precision. Extending these experiments to single molecules outside a diamond chip poses a number of additional challenges, and overcoming them will require the combination of several strategies. To isolate single molecules, they can be embedded in a spin-free matrix layer deposited on the diamond surface 28 or immobilized by a linker   29 . Nuclear dipole interactions can be suppressed using homo-and heteronuclear decoupling 13 , taking advantage of the existing microcoil. Line-widths and spectral complexity can be further reduced by polarizing the constituent nuclei, and by spin dilution and isotope labeling of molecules 30 . Alternatively, measurements of inter-spin couplings will allow constraining the structure and size of molecules 31 . To sensitively detect very weakly coupled nuclei in distant molecules, spin precession can be recorded by repetitive weak measurements 32,33 . Further improvement of sensitivity is possible by optimizing the optical detection efficiency compared to our present setup 34 , and possibly by cryogenic operation 11 . How well these strategies will work we do not know at present, but we believe that the prospect of a general single-molecule MRI technique, which will have many applications in structural biology and chemical analytics, provides sufficient motivation to warrant these efforts.

Methods
Diamond sample. Experiments were performed on a bulk, electronic-grade diamond crystal from ElementSix with dimensions 2 × 2 × 0.5 mm with 110 h i edges and a 100 h i front facet. The diamond was overgrown with a layer structure of 20 nm enriched 12 C (99.99%), 1 nm enriched 13 C (estimated in-grown concentratioñ 5-10%) and a 5 nm cap layer of again enriched 12 C (99.99%). Nitrogen-vacancy (NV) centers were generated by ion-implantation of 15 N with an energy of 5 keV, corresponding to a depth of~5-10 nm. After annealing the sample for NV formation, we had to slightly etch the surface (at 580°C in pure O 2 ) to remove persistent surface fluorescence. The intrinsic nuclear spin of the three NV centers studied in our experiments were confirmed to be of the 15 N isotope. Further characterizations and details on the sample can be found in a recent study (sample B in ref. 35 ).
Coordinate systems. In Supplementary Fig. 2a both laboratory and NV coordinate system are shown in a combined schematic. The laboratory coordinate system (x Lab , y Lab , z Lab ) is defined by the normal vectors to the diamond faces, which lie along 110 h i, 110 h i and 001 h i, respectively. The reference coordinate system of the NV center is defined by its quantization direction, which is labeled z NV and lies along 111 h i. The x NV -and y NV -axis are pointing along the 11 2 h i and 110 h i direction, respectively.
Experimental apparatus. A schematic of the central part of the experimental apparatus is shown in Supplementary Fig. 1. The diamond sample is glued to a 200 μm thick glass piece and thereby held above a quartz slide with incorporated microwave transmission line for electron spin control. Below the quartz slide we placed a high numerical aperture (NA = 0.95) microscope objective for NV excitation with a 532 nm laser and detection using a single photon counting module (SPCM). We applied static, external magnetic bias fields with a cylindrical NdFeB permanent magnet (not shown in Supplementary Fig. 1). The magnet is attached to a motorized, three-axis translation stage. The NV control pulses were generated by an arbitrary waveform generator (Tektronix, AWG5002C) and upconverted by I/Q mixing with a local oscillator to the desired~2.6 GHz.
Planar, high-bandwidth coil. The planar coil is positioned directly above the diamond sample and attached to a metallic holder, which can be laterally shifted to translate the coil. Due to the thickness of the diamond (500 μm) and the glass slide the minimal vertical stand-off of the coil to the NV centers is~700 μm. Design parameters of the planar coil, used in our experiments, are listed in Supplementary Table 1. These were found by numerically maximizing the magnetic field at the position of the NV center, located at a planned vertical stand-off of~700-1000 μm ( Supplementary Fig. 1). The coil had an inductance of ≤2.5 μH and a resistance of ≤0.5 Ω. The coil was manufactured by Sibatron AG (Switzerland) and it is mounted onto a copper plate that acts as a heat-sink, using thermally conducting glue. For the coil control, a National Instruments NI PCI 5421 arbitrary waveform generator was used, to generate voltage signals that controlled a waveform amplifier (Accel Instruments TS-250) which drives the coil current.
Calibration of the coil field ΔB. We calibrated the vector field generated by the coil ΔB using the target NV, coupled to nuclear spins of interest, and two auxiliary NV centers with different crystallographic orientation. All three NV centers were located in close proximity to each other, with a distance of typically ≤5 μm ( Supplementary Fig. 2c). Over this separation the magnetic field of the coil can be assumed to be homogeneous. We determined the orientation of the symmetry axis of many NV centers by moving the permanent magnet over the sample and observing the ODMR splitting. The azimuthal orientation of the target NVs defines the x-axis in the laboratory and NV frame (ϕ = 0). This orientation was the same for all target NV centers investigated in this work. The auxiliary NV centers were selected to be oriented along ϕ a,1 = 90°and ϕ a,2 = 270°( Supplementary Fig. 2b). To calibrate the coil field, we removed the permanent magnet and recorded ODMR spectra for the target NV center and both auxiliary NV centers with the field of the coil activated. In this way, we record in total 6 ODMR lines, with 2 lines per NV center.
A numerical, nonlinear optimization method was used to determine the magnetic field ΔB from these ODMR resonances. For each of the three NV centers we simultaneously minimized the difference between the measured ODMR lines and the eigenvalues of the ground-state Hamiltonian: Here, the magnetic field (ΔB) i acting onto the specific NV center is obtained by a proper rotation of ΔB into the respective reference frame.
Precise alignment of the bias field B 0 . Precise alignment of the external bias field to the quantization axis of the NV center (z-axis) is critical for azimuthal localization measurements, because residual transverse fields of B 0 modify the precession frequencies in the same way as the field of the coil. The coarse alignment of the magnet and a rough adjustment of the magnitude of the field, to~10 mT, was achieved by recording ODMR spectra of the target NV center for different b Reconstructed locations of the four distant nuclear spins 13 C 1-4 . Shaded regions mark the 2σ-confidence area of the localization projected onto (xy, yz, xz)-planes of the coordinate system. Gray points represent carbon lattice positions projected onto the same planes. The origin is set to the expected center of gravity of the spin density at 2.29 Å from the nitrogen nucleus on the NV symmetry axis 25,26 . Due to the inversion symmetry of the hyperfine interaction, our method cannot distinguish between sites in the upper and lower hemisphere; all 13 C are therefore plotted in the upper hemisphere (x, y, z)-positions of the magnet. Afterwards, we iteratively optimized the alignment of the magnet. In each iteration, we reconstructed the vector field B 0 acting on the target NV centers using the method used for the calibration of ΔB. Subsequently, we moved the magnet in the lateral (x, y)-plane of the laboratory frame. The direction and step size was determined from a field map of the permanent magnet and the residual transverse components of the field B 0 . We terminated this iterative process when the residual transverse field components were smaller than 50 μT.
Determination of hyperfine couplings (a || , a ⊥ ) from (f 0 , f −1 , f R ). The hyperfine couplings a || and a ⊥ in the limit 2πf 0 ) a jj ; a ? are given by: In our experiments the hyperfine couplings and the nuclear Larmor frequency f 0 were of similar magnitude, and we used the following transformations 21 to obtain the hyperfine couplings.
Monte Carlo error propagation. Confidence intervals for ϕ and a iso were obtained using standard Monte Carlo error propagation 36 . The Monte Carlo simulation took calibration uncertainties in the external fields B 0 , ΔB and in the observed precession frequencies f m s into account. All parameters subject to uncertainty were assumed to follow a normal distribution. Precession frequencies were determined using a nonlinear, least-squares fitting algorithm and their measurement uncertainties were obtained from the fit error 21 . The uncertainty in the magnetic field components was estimated from the residuals between calculated and measured ODMR lines in the calibration method for B 0 , ΔB.
Nuclear g-factor enhancement. The nuclear g-factor enhancement factor α(m S ) given in Eq. (6) which is also valid in the limit of large magnetic fields γ e B 0 ) D and provides, in principle, more accurate theory values for small B 0 . We have analyzed our experimental data using this expression and found deviations to Eq. (6) that are smaller than the frequency resolution in our experiments.
Code availability. Custom code was programmed to perform the Monte Carlo simulations. The code is available from the corresponding author upon request.

Data availability
The data that support the findings of this study are available from the corresponding author upon request.