Three-dimensional magnetic resonance tomography with sub-10 nanometer resolution

We demonstrate three-dimensional magnetic resonance tomography with a resolution down to 5.99 +- 0.07 nm. Our measurements use lithographically fabricated microwires as a source of three-dimensional magnetic field gradients, which we use to image NV centers in a densely doped diamond by Fourier-accelerated magnetic resonance tomography. We also present a compressed sensing scheme for imaging of a spatially localized ensemble from undersampled data, which allows for a direct visual interpretation without numerical optimization. The resolution achieved in our work approaches the positioning accuracy of site-directed spin labeling, paving the way to three-dimensional structure analysis by magnetic-gradient based tomography.

scanners, these comprise thousands of voxels, and similar data volumes are expected for particle detectors or imaging of a densely spin-labeled protein.Yet, threedimensional Fourier-accelerated imaging at the nanoscale has remained elusive.
Imaging by less scalable techniques has been demonstrated multiple times.Three-dimensional imaging of nuclear spins with atomic resolution has been achieved using the intrinsic field gradient emerging from the magnetic dipole field of a color center 10,11 .However, only the closest few nanometers around a defect can be imaged by this approach so that it is limited to intrinsic spins in the diamond so far.One-dimensional and two-dimensional 12 imaging in static gradients has been demonstrated, including resolving two adjacent centers by the gradient field of a hard-drive write head 13 .Three-dimensional images of intrinsic electron spins in a diamond have been obtained using a static gradient positioned by a scanning probe 14 .Fourier acceleration has remained out of reach of these approaches where gradients cannot be switched within a spectroscopy sequence.Fourier acceleration of imaging has been demonstrated 7,9 using quickly switchable conductors as gradient sources, but experiments with nanoscale resolution have remained limited to one-and twodimensional proofs of the concept.One-dimensional imaging in MRFM by current-driven gradients has very recently even achieved sub-Ångstrom resolution 15 .
Here we demonstrate Fourier-accelerated threedimensional imaging with nanometer-scale resolution.
The key is a device to produce three linearly independent magnetic field gradients from a two-dimensional layout of conductors (Fig. 1a).Three microfabricated wires, arranged in a U-shape structure, create linearly inde- pendent gradient fields in a plane few microns beneath the structure.This device is fabricated via lift-off photolithography on a diamond substrate hosting a dense ([NV] ≈ 0.13 ppb) ensemble of NV centers Fig. 1(a).The U-microstructure consists of a 200 nm gold film on top of a 10 nm thick titanium layer.Each of its three arms is 5 µm long and 500 nm wide.The top arm also serves as a microwave antenna to implement single-qubit gates.We generate switchable magnetic field gradients by sending currents, labeled I 1 , I 2 , and I 3 in Fig. 1(a), into the three arms of the U microstructure.All currents are terminated with a 50 Ω resistance at the same vertex of the U structure.In all measurements a homogeneous bias magnetic field of B 0 ≈ 76 G is applied along one of the four NV axes.This device is used to implement gradient echo pulse sequences, like the one-dimensional example shown in Fig. 1(b-d).A Hahn echo sequences decouples the NV centers from static and slowly fluctuating background fields, to enable T 2 -limited sensing.A magnetic field gradient pulse, created by the current I 1 , is applied during one half of the echo sequence.This phase-encodes the position, because an NV center at point ⃗ x acquires a position-dependent phase shift where ω( ⃗ B I (⃗ x, τ )) denotes the shift in the Larmor frequency induced by the current I and the approximation holds for pulses close to a rectangular shape.At the end of the Hahn echo sequence, this phase shift translates into an oscillatory spin signal For a distribution of NV centers, the oscillatory signals of all centers will linearly superpose to a characteristic beating pattern (assuming a trailing π x /2 pulse).The various ω(⃗ x NV ) can be recovered from this signal by an inverse Fourier transform, creating a one-dimensional image of the NV centers Fig. 1(d).
One major challenge of this experiment consists in creating sufficiently rectangular pulses to satisfy the approximation of eq. ( 1).This requires a stable current supply that moreover has to be controlled with a fast (100 MHz) bandwidth to ensure that the rising and falling edges are quasi-instantaneous, i.e. much shorter than one period of the current-induced Larmor frequency ω( ⃗ B I ).Stability within every current pulse is required, because any variation of ω( ⃗ B(⃗ x, τ )) over the pulse will introduce a chirp in the time domain signal (Fig. 1(c)), which will blur the image in the frequency domain (Fig. 1(d)).Stability between successive experimental repetitions is required, because shot-to-shot fluctuations of the magnetic field induce decoherence (see below).
We experimentally address these constraints by two means (Fig. 1(e)).First, the current pulses are generated by switching a stable voltage source (Keithley 2230G-30-6) using fast switches (ic-Haus HGP), ensuring nearly rectangular pulses.Second, we correct for residual nonlinearities and fluctuations by measurement and online post-processing.We acquire the current integral t 0 I(τ )dτ for every pulse of every experimental repetition by hardware-integration on a fast A/D converter (Spectrum M4i.4451-x8), and use this value to define a new time axis for all photonic measurements that removes chirps.Specifically, we make the following approximation (valid for nearly rectangular pulses and a nearly linear Zeeman shift) Here, ω I,ref (⃗ x) is the shift in Larmor frequency induced by some reference current I 0 and t eff denotes an "effective pulse duration".We thus absorb minor fluctuations of the current over the pulse into this redefined time coordinate t eff .All time-domain plots in this paper will use t eff as time axis unless noted otherwise.This correction also suppresses shot-to-shot fluctuations, improving coherence 16 .We now extend this one-dimensional magnetic resonance tomography to three-dimensions, employing the three magnetic field gradients provided by the currents I 1 , I 2 , I 3 of our device.Note that these gradients are linearly independent if the focal spot of the microscope is placed a few micrometers below the plane of the Ustructure.
During the Hahn Echo sequence the pulses of these three currents are applied consecutively (see Fig. 2 (a)).
Since the accumulated phase will just add up linearly, the resulting spin signal is given by: (assuming a trailing π x /2 pulse).
In analogy to one-dimensional tomography the set of (ω I1 (⃗ x NV ), ω I2 (⃗ x NV ), ω I3 (⃗ x NV )) can be recovered from the three-dimensional time domain data (Fig. 2(b)) by a 3D inverse Fourier transform, forming a threedimensional image.We note that this resulting image is distorted because the gradients, while linearly independent, are not fully orthogonal.While this distortion can in principle be corrected by computing and inverting the exact spatial distribution of the frequency shift ω I (⃗ x), the raw result of the Fourier transform is still a true three-dimensional image.This resulting image reveals individual NV centers in the diamond.Only NV centers within the confocal volume of the microscope can be imaged.For the given NV density 5 − 15 centers are expected to be seen.
One challenge of such multidimensional measurements is the large number of required data points in Fourier space, e.g. 10 6 points in Fig. 2.This number can be reduced by compressed sensing, where only a subset of points is acquired, and the image is reconstructed by numerical techniques like L 1 minimization, exploiting the a priori knowledge that the signal is a sparse set of discrete points 9 .Interestingly, our experimental setting allows for another compressed sensing approach.It does not require elaborate numerical reconstruction and exploits a different kind of a priori knowledge: that the signal is restricted to a narrow region of interest, i.e. a narrow band in frequency space.In this special case, we can implement an effective "zoom" into this region of interest by undersampling the signal in the time domain, lowering the amount of data points.Undersampling leads to aliasing of the signal in frequency space.For suitable parameters, this will shift the signal frequency band to a contiguous low frequency window where it can still be recovered by the inverse Fourier transform, effectively implementing a zoom.We demonstrate a proof of concept of this idea in Fig. 3.The simulated one-dimensional time and frequency domain plots (left part Fig. 3 (a)) display a limited frequency band.When the time-domain signal is undersampled, i.e. sampled at a rate that the Nyquist frequency f Nyq is smaller than the highest signal frequency, any signal at f > f Nyq will be aliased to a frequency where N is the integer minimizing |f −2N f Nyq |.In (Fig. 3 (a), middle plot), the signal band (around f ≈ 30 MHz) is  ,c) The Nyquist frequencies for each gradient direction (x and y axes) were set to be larger than the highest frequency in that direction (no aliasing).For (d,e) and (f,g) the measurement was done using an aliased grid, the aliased factor for each direction can be read in orange above FFT panel of each measurement.In (d,e) the undersampling parameters are such that the signal is flipped in both axes.For (f,g) the signal is flipped in the horizontal axis, but remains unflipped in the vertical axis.The aliased measurement (f,g) reduces the acquisition time by a factor of 10.
close to twice the Nyquist frequency (2 • f Nyq = 35 MHz) and hence aliased to a region close to f = 0 MHz.Note that the aliasing involves mirroring of the signal band, because the signal is at a lower frequency than the closest even multiple of f Nyq .
We show this concept acquiring three separate twodimensional measurements (Fig. 3 (b-g)), which display NV centers in a limited region of interest (upper right corner in Fig. 3 (b,c)), defined by the confocal volume of the microscope .This process of undersampling and aliasing implements a zoom into the region of interest (Fig. 3 (d-e)).Note that this process requires the signal to be confined to a limited window of frequencies.Since frequencies equal to an integer multiple of 2f N yq will appear at the same f obs (Equation .2),signals outside the zoom window will fold back into the signal of interest.To prevent contamination of the resulting image, the signal should be bandpass limited, i.e. values outside the frequency range of the signal should be zero and the Nyquist frequency should not fall below the bandwidth of the NV signal spectrum.For a suitable parameter choice of the undersampling, a zoom can be achieved that exactly covers the region of interest (Fig. 3 (f,g)), allowing for the acquisition of a full image with a greatly reduced number of data points.In the specific example (Fig. 3 (f,g)), the two dimensions are undersampled by a factor of 6 and 3, reducing the number of data points by a factor of 18, i.e. more than an order of magnitude.Note that reconstruction and visualization are still feasible by an inverse Fourier transform (Fig. 3 (f)).L 1 minimization is not required for reconstruction, but can still be implemented to improve the quality of the image and/or further reduce the number of data points required (Fig. 3 (c,e,g)).
We finally note a constraint of the technique.The undersampled data points have to be placed equidistantly in time, since any jitter or chirp will lead to a spectral broadening in the frequency-space image.Since a variation of the gradient current over the duration of a pulse is indistinguishable from a variation in timing, this also places higher demands on the constancy of currents, i.e. the requirement of rectangular current pulses discussed above.
We finally analyze the spatial resolution that is achieved in our measurement.This is defined by the magnitude of the magnetic field gradient, and the spectral resolution of the spectroscopy.The frequency resolution of Fourier-transformed data is given as the inverse of the length of the time domain signal.Analogous to that, the frequency (and thus spatial) resolution of our magnetic resonance tomography depends on how long we can make the gradient pulse length and still observe an oscillatory spin signal.The longest usable pulse is limited by the fact that the spin signal decays over time on a timescale of ≈ 10 µs (see e.g.Fig. 4 (a-b)) because of shot-to-shot fluctuations of the gradient currents, which result in the decoherence of the NV centers.Denoting the timescale of this decay by T 2,I (i.e. the coherence time in the presence of the gradient current), the frequency resolution is where ∆f denotes full width at half-maximum (FWHM) of the peak in frequency space.See supplementary information for explanation of the √ 2 π factor.We extract T 2,I from a long 1D tomography data set, extending to several multiples of T 2,I .We calculate the Fourier transform for short time window and "slide" this window over the whole range of the time domain signal.The resulting spectrogram (Fig. 4 (b)) shows the evolution of the NV spectrum with increasing gradient pulse lengths.The signal from a single spin appears as a horizontal line in a specific frequency band.The decaying power of the signal with increasing time in this band defines the SNR over the measurement (see supplementary).We Fit the SNR curve (Fig. 4 (c)) with a Gaussian (e −t 2 eff /2T 2 2 ) to obtain T 2,I .For the data of Fig. 4 we thus arrive at a coherence time of T 2,I2 = 8.64 ± 0.1 µs.Combined with a gradient of ||∇ω(⃗ x)||/2π = 6.34 kHz/nm, obtained from a numerical simulation of the gradient field (see supplementary), this corresponds to a spatial resolution of σ x,I2 = 8.22 ± 0.10 nm.Similarly, we obtain σ x,I1 = 5.99 ± 0.07 nm and σ x,I3 = 14.47 ± 0.50 nm for the other two gradient currents.This resolution could be limited by several effects.First, shot-to-shot fluctuations of the current could shorten T 2,I .We try to suppress this by hardware integrating every single current pulse (see above) and applying post-processing corrections, but this process is equally limited by electronic noise at a lower level.Second, a spatial drift of the current path between successive experimental repetitions can equally lead to a decrease of T 2 .A spatial drift could arise from heat expansion of the diamond and the conductors, but an expansion on the level of 10 −3 would require a temperature difference of ≈ 1000 K. Which seems unlikely.A drift of the current path within the conductor, due to local heating appears more reasonable.Intriguingly the product of ω NV T 2 , i.e. the relative stability of the gradient field differs between the three wires.This tentatively suggests that spatial drifts of the current in the wires are the limiting factor rather than electrical fluctuations, which would be expected to be the same in all wires.
In summary, we have demonstrated Fourieraccelerated 3D imaging of single spins with nanoscale resolution.We have also presented a compressed sensing scheme, which exploits a limited field of view, rather than sparseness of the data.Our experiments demonstrate that resolution in the sub-10 nm range can be achieved by switchable magnetic field gradients.
While our experiment has been performed on NV centers inside the diamonds, the device and measurement technique could equally be applied to dark spins outside of the diamond.Here a single NV center would merely serve as a detector to enable electron/nuclear spin spectroscopy on spins, while the entire process of imaging could be performed by the device presented here.Our compressed sensing technique of "Fourier zooming" will be especially advantageous in this setting where all the spins are confined to the nanoscale detection volume of a shallow NV center.Such a direct 3D imaging technique could image an arbitrary number of spins and constrain inter-spin distances larger than 80 Å, which is not possible by current electron spin resonance spectroscopy.Shrinking the structure by one order of magnitude would even push the resolution into the range of Å. Notably the T 2 of established spin labels is sufficiently long for the spectroscopy presented here 17 .

Figure 1 .
Figure 1.Experimental setup and one-dimensional magnetic resonance tomography of NV centers.(a) Electron micrograph of a device as used in the present study.Currents in the three gold wires of a microfabricated U-Structure create three linearly independent magnetic field gradients in the densely doped diamond below the structure.(b) Pulse sequence for one-dimensional imaging.A magnetic gradient pulse (length t eff ) inserted into a Hahn echo sequence phase-encodes position.The NV spin state is initialized and the spin projection ⟨ Ŝz⟩ is read out optically.(c) Measurement result of (b).π/2x and π/2y denote the phase of the trailing π/2 pulse in (b).(d) Fourier transform of a dataset like (c) extending to t = 60 µs.Every NV center gives rise to one peak at the Larmor frequency set by the magnetic field BI 1 (d) of the wire.d denotes the distance from wire I1.(e) Experimental setup.A single microwave generator, a 90 • splitter and two microwave switches are used to implement the Hahn Echo sequence.A confocal microscope with an avalanche photodiode (APD) as a detector is used for NV center polarization and readout.The gradient currents I1, I2, I3 are created from a constant voltage source and can be pulsed by a fast switch.The voltage drop across the resistor is recorded by an A/D-converter and the pulse integral Idt is saved for every single current pulse.

Figure 2 .
Figure 2. Three-dimensional magnetic resonance tomography.(a) Pulse sequence.The sequence of Fig. 1 is extended to contain three magnetic gradient pulses from different wires.(b) Time domain data recorded from the sequence of (a) ending with π/2x.(c) Three-dimensional Fourier transform of the data in (b).The plot shows the square of the absolute value (spectral power) of the Fourier transform.d1, d2, d3 denote the distance to the respective wire (see supplementary).The bottom, left, and back faces show projections of the 3D data.

Figure 3 .
Figure 3. Aliasing magnification and speed-up of Fourier magnetic imaging.The plots in (a) show a simulated signal and its Fourier transform.Going through the columns from left to right the sampling rate is reduced resulting in a slower oscillation (red trace) and a lower Nyquist frequency (black-dashed line).Aliasing around the closest even multiple of the undersampled Nyquist frequency shifts the signal to a window close to f = 0, but does not change its shape.The Nyquist frequency of the undersampled signal has to be chosen large enough to cover the entire signal bandwidth.The signal is shifted to negative frequencies, and hence appears flipped in a frequency axis using |f |, if it sits on the left of the closest even multiple of the Nyquist frequency.Panels b-g display measured 2D images of NV centers acquired by taking the FFT of the time domain signal (b, d, f) or by doing an L1 minimization of the time domain signal (c, e, g).(b,c) The Nyquist frequencies for each gradient direction (x and y axes) were set to be larger than the highest frequency in that direction (no aliasing).For (d,e) and (f,g) the measurement was done using an aliased grid, the aliased factor for each direction can be read in orange above FFT panel of each measurement.In (d,e) the undersampling parameters are such that the signal is flipped in both axes.For (f,g) the signal is flipped in the horizontal axis, but remains unflipped in the vertical axis.The aliased measurement (f,g) reduces the acquisition time by a factor of 10.

Figure 4 .
Figure 4. Benchmarking of the spatial resolution for I2.(a) time-domain signal of a one-dimensional tomography (sequence of Fig. 1(b)).Excerpts at different time windows are shown.(b) Spectrogram (windowed Fourier transform) of (a).The signal produced by the NV centers decays over a timescale of ≈ 10 µs.(c) Signal-to-noise ratio of (b), computed by integrating the power in the signal window marked in (b) and referencing it to the noise observed outside this window (see supplementary).(c) A Gaussian fit to the data yields a decay timescale T2,I 2 = 8.64 ± 0.1 µs.(d) Fourier transform (absolute value) of the time domain signal in (a).