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

Mohammad T Amawi  Institute for Physics, Quantum technology, University of
Rostock, Germany  Physics Department, Technical University of Munich, Germany
Department of Life, Light and Matter, Quantum Technology, University of
Rostock,Germany  Munich Center for Quantum Science and Technology (MCQST),
Germany  Andrii Trelin  Institute for Physics, Quantum technology, University
of Rostock, Germany  Department of Life, Light and Matter, Quantum Technology,
University of Rostock,Germany  You Huang  Institute for Physics, Quantum
technology, University of Rostock, Germany  Department of Life, Light and
Matter, Quantum Technology, University of Rostock,Germany  CAS Key Laboratory
of Microscale Magnetic Resonance and School of Physical Sciences, University
of Science and Technology of China, Hefei 230026, China  CAS Center for
Excellence in Quantum Information and Quantum Physics, University of Science
and Technology of China, Hefei 230026, China  Paul Weinbrenner  Institute for
Physics, Quantum technology, University of Rostock, Germany  Department of
Life, Light and Matter, Quantum Technology, University of Rostock,Germany
Francesco Poggiali  Institute for Physics, Quantum technology, University of
Rostock, Germany  Department of Life, Light and Matter, Quantum Technology,
University of Rostock,Germany  Munich Center for Quantum Science and
Technology (MCQST), Germany  Joachim Leibold  Physics Department, Technical
University of Munich, Germany  Munich Center for Quantum Science and
Technology (MCQST), Germany  Martin Schalk  Physics Department, Technical
University of Munich, Germany  Munich Center for Quantum Science and
Technology (MCQST), Germany  Friedemann Reinhard  friedemann.reinhard@uni-
rostock.de Institute for Physics, Quantum technology, University of Rostock,
Germany  Department of Life, Light and Matter, Quantum Technology, University
of Rostock,Germany  Munich Center for Quantum Science and Technology (MCQST),
Germany

(August 28, 2024)

###### Abstract

We demonstrate three-dimensional magnetic resonance tomography with a
resolution down to $5.99\pm 0.07\leavevmode\nobreak\ $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.

In recent years, various nano-sensors, most prominently magnetic resonance
force microscopy (MRFM) and nitrogen-vacancy (NV) centers, have enabled the
detection of small ensembles of electron[1, 2] and nuclear[3, 4, 5] spins,
partially down to the level of single spins. Translating this power from a
mere detection to a three-dimensional imaging technique promises
transformative applications. Three-dimensional imaging of color centers would
enable selective addressing and readout of networks of coherently coupled
color centers in densely doped samples[6, 7], or the detection of elementary
particles by high-resolution mapping of the crystal strain they induce upon
impact[8]. Applied to electron spins, it would enable three-dimensional
imaging of spin-labeled proteins. Such a technique would in particular provide
distance constraints for label distances of $>80\leavevmode\nobreak\
\textup{\AA}$ and for proteins labeled with arbitrary many electron spins,
filling two blind spots of present electron spin resonance spectroscopy.
Applied to nuclear spins, it would provide an ultimate microscope, able to
image within opaque samples with label-free chemical contrast.

Conceptually, the step from detection to imaging is straightforward. Applying
a magnetic field gradient is all it takes to turn a magnetic resonance
spectrum into a one-dimensional image. Multiple gradients along linearly
independent directions can encode multi-dimensional images, most beautifully
illustrated in the output of clinical magnetic resonance imaging scanners. If
the gradients can be switched faster than the duration of the spectroscopy
sequence, Fourier-accelerated techniques can acquire extended volumes in
reasonable time[9]. In clinical scanners, these comprise thousands of voxels,
and similar data volumes are expected for particle detectors or imaging of a
densely spin-labeled protein. Yet, three-dimensional 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[9, 7] using quickly
switchable conductors as gradient sources, but experiments with nanoscale
resolution have remained limited to one- and two-dimensional proofs of the
concept. One-dimensional imaging in MRFM by current-driven gradients has very
recently even achieved sub-Ångstrom resolution[15].

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_{\textit{eff}}$)
inserted into a Hahn echo sequence phase-encodes position. The NV spin state
is initialized and the spin projection $\langle\hat{S}_{z}\rangle$ is read out
optically. (c) Measurement result of (b). $\pi/2_{x}$ and $\pi/2_{y}$ denote
the phase of the trailing $\pi/2$ pulse in (b). (d) Fourier transform of a
dataset like (c) extending to $t=60\leavevmode\nobreak\ $µs. Every NV center
gives rise to one peak at the Larmor frequency set by the magnetic field
$B_{I_{1}}(d)$ of the wire. $d$ denotes the distance from wire $I_{1}$. (e)
Experimental setup. A single microwave generator, a $90^{\circ}$ 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
$I_{1},I_{2},I_{3}$ 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 $\int Idt$ is saved for every single
current pulse.

Here we demonstrate Fourier-accelerated three-dimensional 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 independent 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 ($[\textit{NV}]\approx
0.13\leavevmode\nobreak\ $ppb) ensemble of NV centers Fig. 1(a). The
U-microstructure consists of a $200\leavevmode\nobreak\ $nm gold film on top
of a $10\leavevmode\nobreak\ $nm thick titanium layer. Each of its three arms
is 5 µm long and $500\leavevmode\nobreak\ $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 I1, I2, and I3 in Fig.
1(a), into the three arms of the U microstructure. All currents are terminated
with a $50\leavevmode\nobreak\ \Omega$ resistance at the same vertex of the U
structure. In all measurements a homogeneous bias magnetic field of
$B_{0}\approx 76\leavevmode\nobreak\ $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 $\vec{x}$ acquires a
position-dependent phase shift

$\phi(\vec{x},t)=\int_{0}^{t}\omega(\vec{B}_{I}(\vec{x},\tau))d\tau\approx\omega(\vec{B}_{I}(\vec{x}))t$
(1)

where $\omega(\vec{B}_{I}(\vec{x},\tau))$ 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

$\langle\hat{S}_{z}\rangle(t)=\begin{cases}(1+\cos(\omega(\vec{x})t)/2&\text{
trailing $\pi_{x}/2$ pulse}\\\ (1+\sin(\omega(\vec{x})t)/2&\text{ trailing
$\pi_{y}/2$ pulse}\\\ \end{cases}$

For a distribution of NV centers, the oscillatory signals of all centers will
linearly superpose to a characteristic beating pattern

$\langle\hat{S}_{z}\rangle(t)=\Bigg{(}1+\sum_{\vec{x}_{\textit{NV}}}\cos(\omega(\vec{x}_{\textit{NV}})t)\Bigg{)}/2$

(assuming a trailing $\pi_{x}/2$ pulse). The various
$\omega(\vec{x}_{\textit{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\leavevmode\nobreak\ $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 $\omega(\vec{B}_{I})$. Stability within
every current pulse is required, because any variation of
$\omega(\vec{B}(\vec{x},\tau))$ 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 $\int_{0}^{t}I(\tau)d\tau$ 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)

$\int_{0}^{t}\omega(\vec{B}_{I}(\vec{x},\tau))d\tau\approx\omega_{I,\textit{ref}}\frac{\int_{0}^{t}I(\tau)d\tau}{I_{0}}=\omega_{I,\textit{ref}}t_{\textit{eff}}$

Here, $\omega_{I,\textit{ref}}(\vec{x})$ is the shift in Larmor frequency
induced by some reference current $I_{0}$ and $t_{\textit{eff}}$ denotes an
"effective pulse duration". We thus absorb minor fluctuations of the current
over the pulse into this redefined time coordinate $t_{\textit{eff}}$. All
time-domain plots in this paper will use $t_{\textit{eff}}$ as time axis
unless noted otherwise. This correction also suppresses shot-to-shot
fluctuations, improving coherence[16].

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 $\pi/2_{x}$. (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. $d_{1},d_{2},d_{3}$ denote the
distance to the respective wire (see supplementary). The bottom, left, and
back faces show projections of the 3D data.

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 U-structure.

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:

$\langle\hat{S}_{z}\rangle(t)=\Bigg{(}1+\sum_{\vec{x}_{\textit{NV}}}\cos(\omega_{I_{1}}(\vec{x}_{\textit{NV}})t_{1}+\\\
\omega_{I_{2}}(\vec{x}_{\textit{NV}})t_{2}+\omega_{I_{3}}(\vec{x}_{\textit{NV}})t_{3})\Bigg{)}/2$

(assuming a trailing $\pi_{x}/2$ pulse).

In analogy to one-dimensional tomography the set of
($\omega_{I_{1}}(\vec{x}_{\textit{NV}})$,
$\omega_{I_{2}}(\vec{x}_{\textit{NV}})$,
$\omega_{I_{3}}(\vec{x}_{\textit{NV}})$) can be recovered from the three-
dimensional time domain data (Fig. 2(b)) by a 3D inverse Fourier transform,
forming a three-dimensional 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
$\omega_{I}(\vec{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.

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.

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_{\textit{Nyq}}$ is smaller
than the highest signal frequency, any signal at $f>f_{\textit{Nyq}}$ will be
aliased to a frequency

$f_{\textit{obs}}=|f-2N\cdot f_{\textit{Nyq}}|\leavevmode\nobreak\ ,$ (2)

where $N$ is the integer minimizing $|f-2Nf_{\textit{Nyq}}|$. In (Fig. 3 (a),
middle plot), the signal band (around $f\approx 30\leavevmode\nobreak\ $MHz)
is close to twice the Nyquist frequency ($2\cdot
f_{\textit{Nyq}}=35\leavevmode\nobreak\ $MHz) and hence aliased to a region
close to $f=0\leavevmode\nobreak\ $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_{\textit{Nyq}}$.

We show this concept acquiring three separate two-dimensional 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_{Nyq}$ 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.

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
$\approx 10\leavevmode\nobreak\ $µ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 $T_{2,I_{2}}=8.64\pm
0.1\leavevmode\nobreak\ $µs. (d) Fourier transform (absolute value) of the
time domain signal in (a).

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 $\approx 10\leavevmode\nobreak\ $µ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

$\Delta f=\frac{\sqrt{2}}{\pi T_{2,I}}$

where $\Delta f$ denotes full width at half-maximum (FWHM) of the peak in
frequency space. See supplementary information for explanation of the
$\frac{\sqrt{2}}{\pi}$ 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_{\textit{eff}}^{2}/2T_{2}^{2}}$) to
obtain $T_{2,I}$. For the data of Fig. 4 we thus arrive at a coherence time of
$T_{2,I_{2}}=8.64\pm 0.1\leavevmode\nobreak\ \mu$s. Combined with a gradient
of $||\nabla\omega(\vec{x})||/2\pi=6.34\leavevmode\nobreak\ $kHz/nm, obtained
from a numerical simulation of the gradient field (see supplementary), this
corresponds to a spatial resolution of $\sigma_{x,I_{2}}=8.22\pm
0.10\leavevmode\nobreak\ $nm. Similarly, we obtain $\sigma_{x,I_{1}}=5.99\pm
0.07\leavevmode\nobreak\ $nm and $\sigma_{x,I_{3}}=14.47\pm
0.50\leavevmode\nobreak\ $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 $\approx 1000\leavevmode\nobreak\ $K. Which seems
unlikely. A drift of the current path within the conductor, due to local
heating appears more reasonable. Intriguingly the product of
$\omega_{\textit{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 Fourier-accelerated 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\leavevmode\nobreak\ \textup{\AA}$, 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].

This work has been supported by the Deutsche Forschungsgemeinschaft (DFG,
grants RE3606/1-2, RE3606/3-1 and excellence cluster MCQST EXC-2111-390814868,
SFB 1477 “Light–Matter Interactions at Interfaces” (Project No. 441234705))
and the European Union (ASTERIQS, Grant Agreement No. 820394). Y.H.
acknowledges financial support from the China Scholarship Council. The authors
acknowledge the help of Regina Lange and Anja Clasen with taking the SEM
picture and helpful technical discussions with John Marohn.

## References

  * Rugar  _et al._ [2004] D. Rugar, R. Budakian, H. J. Mamin, and B. W. Chui, “Single spin detection by magnetic resonance force microscopy,” Nature 430, 329–332 (2004).
  * Shi  _et al._ [2015] F. Shi, Q. Zhang, P. Wang, H. Sun, J. Wang, X. Rong, M. Chen, C. Ju, F. Reinhard, H. Chen, J. Wrachtrup, J. Wang, and J. Du, “Single-protein spin resonance spectroscopy under ambient conditions,” Science 347, 1135–1138 (2015).
  * Degen  _et al._ [2009] C. L. Degen, M. Poggio, H. J. Mamin, C. T. Rettner, and D. Rugar, “Nanoscale magnetic resonance imaging,” Proceedings of the National Academy of Sciences 106, 1313–1317 (2009).
  * Staudacher  _et al._ [2013] T. Staudacher, F. Shi, S. Pezzagna, J. Meijer, J. Du, C. A. Meriles, F. Reinhard, and J. Wrachtrup, “Nuclear Magnetic Resonance Spectroscopy on a (5-Nanometer)3 Sample Volume,” Science 339, 561–563 (2013).
  * Mamin  _et al._ [2013] H. J. Mamin, M. Kim, M. H. Sherwood, C. T. Rettner, K. Ohno, D. D. Awschalom, and D. Rugar, “Nanoscale Nuclear Magnetic Resonance with a Nitrogen-Vacancy Spin Sensor,” Science 339, 557–560 (2013).
  * Zhang  _et al._ [2017] H. Zhang, K. Arai, C. Belthangady, J. C. Jaskula, and R. L. Walsworth, “Selective addressing of solid-state spins at the nanoscale via magnetic resonance frequency encoding,” npj Quantum Information 3, 31 (2017).
  * Artzi  _et al._ [2023] Y. Artzi, O. Zgadzai, B. Solomon, and A. Blank, “Three-dimensional fourier imaging of thousands of individual solid-state quantum bits – a tool for spin-based quantum technology,” Physica Scripta 98, 035815 (2023), publisher: IOP Publishing.
  * Rajendran  _et al._ [2017] S. Rajendran, N. Zobrist, A. O. Sushkov, R. Walsworth, and M. Lukin, “A method for directional detection of dark matter using spectroscopy of crystal defects,” Physical Review D 96, 035009 (2017), publisher: American Physical Society.
  * Arai  _et al._ [2015] K. Arai, C. Belthangady, H. Zhang, N. Bar-Gill, S. J. DeVience, P. Cappellaro, A. Yacoby, and R. L. Walsworth, “Fourier magnetic imaging with nanoscale resolution and compressed sensing speed-up using electronic spins in diamond,” Nature Nanotechnology 10, 859–864 (2015).
  * Zopes  _et al._ [2018] J. Zopes, K. S. Cujia, K. Sasaki, J. M. Boss, K. M. Itoh, and C. L. Degen, “Three-dimensional localization spectroscopy of individual nuclear spins with sub-Angstrom resolution,” Nature Communications 9, 4678 (2018), number: 1 Publisher: Nature Publishing Group.
  * Abobeih  _et al._ [2019] M. H. Abobeih, J. Randall, C. E. Bradley, H. P. Bartling, M. A. Bakker, M. J. Degen, M. Markham, D. J. Twitchen, and T. H. Taminiau, “Atomic-scale imaging of a 27-nuclear-spin cluster using a single-spin quantum sensor,” arXiv:1905.02095 [cond-mat, physics:quant-ph] (2019), arXiv: 1905.02095.
  * da Silva Barbosa  _et al._ [2021] J. F. da Silva Barbosa, M. Lee, P. Campagne-Ibarcq, P. Jamonneau, Y. Kubo, S. Pezzagna, J. Meijer, T. Teraji, D. Vion, D. Esteve, R. W. Heeres, and P. Bertet, “Determining the position of a single spin relative to a metallic nanowire,” Journal of Applied Physics 129, 144301 (2021), publisher: American Institute of Physics.
  * Bodenstedt  _et al._ [2018] S. Bodenstedt, I. Jakobi, J. Michl, I. Gerhardt, P. Neumann, and J. Wrachtrup, “Nanoscale spin manipulation with pulsed magnetic gradient fields from a hard disc drive writer,” Nano letters 18, 5389–5395 (2018).
  * Grinolds  _et al._ [2014] M. S. Grinolds, M. Warner, K. D. Greve, Y. Dovzhenko, L. Thiel, R. L. Walsworth, S. Hong, P. Maletinsky, and A. Yacoby, “Subnanometre resolution in three-dimensional magnetic resonance imaging of individual dark spins,” Nature Nanotechnology 9, 279–284 (2014).
  * Haas  _et al._ [2022] H. Haas, S. Tabatabaei, W. Rose, P. Sahafi, M. Piscitelli, A. Jordan, P. Priyadarsi, N. Singh, B. Yager, P. J. Poole, D. Dalacu, and R. Budakian, “Nuclear magnetic resonance diffraction with subangstrom precision,” Proceedings of the National Academy of Sciences 119, e2209213119 (2022), publisher: Proceedings of the National Academy of Sciences.
  * Braunbeck  _et al._ [2021] G. Braunbeck, M. Kaindl, A. M. Waeber, and F. Reinhard, “Decoherence mitigation by real-time noise acquisition,” Journal of Applied Physics 130, 054302 (2021).
  * Soetbeer  _et al._ [2018] J. Soetbeer, M. Hülsmann, A. Godt, Y. Polyhach, and G. Jeschke, “Dynamical decoupling of nitroxides in o-terphenyl: a study of temperature, deuteration and concentration effects,” Physical Chemistry Chemical Physics 20, 1615–1628 (2018).