# Nuclear Spin Metrology with Nitrogen Vacancy Center in Diamond for
Axion Dark Matter Detection

So Chigusa  0000-0001-6005-4447 Theoretical Physics Group, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Berkeley Center for Theoretical
Physics, Department of Physics, University of California, Berkeley, CA 94720,
USA  International Center for Quantum-field Measurement Systems for Studies of
the Universe and Particles (QUP), High Energy Accelerator Research
Organization (KEK), 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan  Masashi Hazumi
International Center for Quantum-field Measurement Systems for Studies of the
Universe and Particles (QUP), High Energy Accelerator Research Organization
(KEK), 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan  Institute of Particle and
Nuclear Studies (IPNS), KEK, Tsukuba, Ibaraki 305-0801, Japan  Japan Aerospace
Exploration Agency (JAXA), Institute of Space and Astronautical Science
(ISAS), Sagamihara, Kanagawa 252-5210, Japan  Kavli Institute for the Physics
and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, The University of
Tokyo, Kashiwa, Chiba 277-8583, Japan  The Graduate University for Advanced
Studies (SOKENDAI), Miura District, Kanagawa 240-0115, Hayama, Japan  Ernst
David Herbschleb  0000-0001-7579-1765 Institute for Chemical Research, Kyoto
University, Gokasho, Uji-city, Kyoto 611-0011, Japan  Yuichiro Matsuzaki
Department of Electrical, Electronic, and Communication Engineering, Faculty
of Science and Engineering, Chuo university, 1-13-27, Kasuga, Bunkyo-ku, Tokyo
112-8551, Japan  Norikazu Mizuochi  International Center for Quantum-field
Measurement Systems for Studies of the Universe and Particles (QUP), High
Energy Accelerator Research Organization (KEK), 1-1 Oho, Tsukuba, Ibaraki
305-0801, Japan  Institute for Chemical Research, Kyoto University, Gokasho,
Uji-city, Kyoto 611-0011, Japan  Center for Spintronics Research Network,
Kyoto University, Uji, Kyoto 611-0011, Japan  Kazunori Nakayama  International
Center for Quantum-field Measurement Systems for Studies of the Universe and
Particles (QUP), High Energy Accelerator Research Organization (KEK), 1-1 Oho,
Tsukuba, Ibaraki 305-0801, Japan  Department of Physics, Tohoku University,
Sendai, Miyagi 980-8578, Japan

###### Abstract

We present a method to directly detect the axion dark matter using nitrogen
vacancy centers in diamonds. In particular, we use metrology leveraging the
nuclear spin of nitrogen to detect axion-nucleus couplings. This is achieved
through protocols designed for dark matter searches, which introduce a novel
approach of quantum sensing techniques based on the nitrogen vacancy center.
Although the coupling strength of the magnetic fields with nuclear spins is
three orders of magnitude smaller than that with electron spins for
conventional magnetometry, the axion interaction strength with nuclear spins
is the same order of magnitude as that with electron spins. Furthermore, we
can take advantage of the long coherence time by using the nuclear spins for
the axion dark matter detection. We show that our method is sensitive to a
broad frequency range $\lesssim$100\text{\,}\mathrm{H}\mathrm{z}$$
corresponding to the axion mass
$m_{a}\lesssim$4\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V}$$. We
present the detection limit of our method for both the axion-neutron and the
axion-proton couplings and discuss its significance in comparison with other
proposed ideas.

††preprint: KEK-QUP-2024-0013, TU-1233

###### Contents

  1. I Introduction
  2. II NV center metrology
    1. II.1 NV center in diamonds
    2. II.2 Dc magnetometry
    3. II.3 Ac magnetometry
  3. III Axion dark matter
    1. III.1 Setup
    2. III.2 Axion interaction with the ${\ce{{}^{1}4N}}$ spin
  4. IV Sensitivity estimation
  5. V Results
  6. VI Discussion and conclusion
  7. A Spin synthesis
  8. B Calculation of the power spectral density and the quantum noise

## I. Introduction

The existence of dark matter (DM) is one of the most important hints of new
physics in particle physics. While the relic abundance of DM in the current
Universe is known through various cosmological and astrophysical observations,
such as the galaxy rotation curve, weak gravitational lensing, and the cosmic
microwave background (see, e.g. [1, 2, 3, 4, 5, 6, 7, 8, 9]), other properties
of DM remain unrevealed. One approach to study these properties is through
direct detection of DM in lab-based experiments. Given the variety of DM
candidates that can explain the relic abundance, numerous approaches have been
taken to investigate different types of DM interactions with standard model
(SM) particles (see [10, 11, 12] for a review).

Among various DM candidates, the axion stands out as a promising candidate,
motivated by several contexts. The term “axion” can refer to the QCD axions,
such as those proposed in [13, 14, 15, 16], which are introduced to solve the
strong CP problem [17, 18, 19], or the axion-like particles, which represent a
broader set of pseudoscalar particles often predicted in low-energy effective
theories emerging from the string theory [20, 21, 22, 23, 24, 25, 26, 27, 28].
Generally, the axion interacts with SM gauge bosons and fermions, each
interaction controlled by a model-dependent coupling constant. Therefore,
developing various strategies to investigate different couplings is essential
to discover axions and differentiate between axion models.

In this paper, we explore the nitrogen-vacancy (NV) center in diamond, a well-
studied multimodal quantum sensing device, as an apparatus for axion DM
searches. Unlike a previous study [29], where some of the authors used NV
center metrology based on electron spins to detect DM signals, we utilize the
nuclear spin of the ${\ce{{}^{1}4N}}$ atom to search for signals induced by
axion-nucleus couplings. This approach provides new constraints on the axion-
neutron coupling $g_{ann}$ and the axion-proton coupling $g_{app}$, which are
independent of the constraint on the axion-electron coupling $g_{aee}$
obtained in [29]. Our method can be viewed as magnetometry based on nuclear
spins. Although this procedure is not well-suited for detecting ordinary
magnetic fields due to their weak coupling to nuclear spins, it is crucial for
axion DM searches because $g_{ann}$, $g_{app}$, and $g_{aee}$ are independent
parameters. We emphasize that our approach, illustrated in Fig. 2 and Fig. 3,
is completely new and specifically aimed at axion searches.

The rest of the paper is organized as follows. In Sec. II, we review NV center
metrology, starting with an overview of the NV center system (Sec. II.1) and
explaining the protocols used for dc (Sec. II.2) and ac (Sec. II.3)
magnetometry. Sec. III reviews axion properties, where we derive the axion
interaction Hamiltonian with elementary particles (Sec. III.1) and the
${\ce{{}^{1}4N}}$ spin (Sec. III.2). We discuss our detection limit estimation
in Sec. IV and present constraints on the axion coupling constants in Sec. V.
Finally, we provide concluding remarks in Sec. VI.

## II. NV center metrology

### II.1. NV center in diamonds

The NV center is a complex composed of a substitutional nitrogen and an
adjacent vacancy. Among the various possible charge states, $\mathrm{NV}^{-}$
is often used for quantum sensing, where two remnant electrons are localized
to the position of the vacancy. These two electrons form the orbital-singlet,
spin-triplet system at the lowest energy levels. The other possible
combinations of angular momenta, which include the orbital triplet and/or the
spin-singlet states, correspond to excited states. The electron system is
excited to an orbital-triplet state by injecting
$532\text{\,}\mathrm{n}\mathrm{m}$ green light, which can relax either
directly with emitting $600$–$800\text{\,}\mathrm{n}\mathrm{m}$ red light or
through spin-singlet states with emitting infrared light. Since the
probability of direct relaxation depends on whether the initial state of the
two-electron spin $\vec{S}$ is $\ket{S^{z}=0}$ or $\ket{S^{z}=\pm}$, we can
read out the spin state information through the fluorescence measurement. When
relaxing through the spin-singlet states, the electron spin usually ends up in
the lowest energy state $\ket{S^{z}=0}$, thus making the whole procedure work
also as laser cooling.

In addition to the electrons at the NV center, the substitutional nitrogen
also possesses a (nuclear) spin degree of freedom, $\vec{I}$. Since $\sim
99.6\%$ of the nitrogens in nature are ${\ce{{}^{1}4N}}$ with spin $I=1$, we
focus on this isotope. Including the hyperfine interaction between the
electron and nuclear spins, the dynamics of the NV center spin system is
governed by the Hamiltonian

$\displaystyle H$ $\displaystyle=H_{\parallel}+H_{\perp},$ (1)

where the first (second) term corresponds to the interactions parallel to
(perpendicular to) the $z$-axis, which is defined by the NV axis.111 According
to this definition of the $z$-axis, four different configurations of the NV
center in the diamond lattice are effectively distinguished by choosing four
different sets of local coordinates, resulting in different effective magnetic
fields. This affects the resonance frequency of the Rabi cycle we will discuss
below, thus the succeeding spin operation is effective only for a part of four
configurations. However, note that the orientation can be aligned with a
specific fabrication process [30, 31, 32, 33].  They are given by [34]

$\displaystyle H_{\parallel}$
$\displaystyle=\Delta_{0}S^{z2}+Q_{0}I^{z2}+B^{z}(\gamma_{e}S^{z}+\gamma_{N}I^{z})+A_{\parallel}S^{z}I^{z},$
(2) $\displaystyle H_{\perp}$
$\displaystyle=\gamma_{e}\vec{B}_{\perp}\cdot\vec{S}_{\perp}+\gamma_{N}\vec{B}_{\perp}\cdot\vec{I}_{\perp}+A_{\perp}\vec{S}_{\perp}\cdot\vec{I}_{\perp},$
(3)

where $\vec{B}$ is an external magnetic field, while the subscript ⟂ of a
vector denotes components perpendicular to the $z$-axis. $\Delta_{0}\simeq
2\pi\times$2.87\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ and
$Q_{0}\simeq-2\pi\times$4.95\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ are the
zero-field splitting of electron spins and the nuclear quadrupole interaction
parameter, respectively. The gyromagnetic ratios for electron and nuclear
spins are given respectively by $\gamma_{e}\simeq
2\pi\times$28\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}\mathrm{/}\mathrm{T}$$,
$\gamma_{N}=2\pi\times$3.08\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}\mathrm{/}\mathrm{T}$$
[35]. The size of the hyperfine interaction is measured as
$A_{\parallel}\simeq-2\pi\times$2.16\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$
and $A_{\perp}=-2\pi\times$2.62\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$.

Figure 1: Energy levels of the two-qubit subsystem, with $\Delta$ and $Q$
described in the text, and $A_{\parallel}$ the size of the hyperfine
splitting.

The quantum state of electron and nuclear spins can be manipulated by the Rabi
cycle. To see this in more detail, let us first pick up two of the electron
spin states, say $\ket{S^{z}=-,0}$, and two of the nuclear spin states, say
$\ket{I^{z}=0,+}$, to constrain ourselves to an effective two-qubit subsystem
spanned by $\\{\ket{S^{z}I^{z}}=\ket{-0},\ket{-+},\ket{00},\ket{0+}\\}$. For
later convenience, we assume the decomposition
$\vec{B}=B_{0}\hat{z}+\vec{B}_{\perp}\cos\omega t$ with $\hat{z}$ the unit
vector along the $z$-axis. If we treat $H_{\perp}$, or $\vec{B}_{\perp}$ and
$A_{\perp}$, as a perturbation, four states labelled by $\ket{S^{z}I^{z}}$ are
energy eigenstates, whose energy levels are shown in Fig. 1 with
$\Delta\equiv\Delta_{0}-\gamma_{e}B_{0}$ and $Q\equiv Q_{0}+\gamma_{N}B_{0}$.
The effect of the oscillating transverse magnetic field is described by the
effective Hamiltonian

$\displaystyle H_{\mathrm{eff}}(t)=$
$\displaystyle\left[\frac{1}{\sqrt{2}}\gamma_{e}\left(\vec{B}_{\perp}\cdot\vec{\sigma}_{e}\right)+\frac{1}{\sqrt{2}}\gamma_{N}\left(\vec{B}_{\perp}\cdot\vec{\sigma}_{n}\right)\right]$
$\displaystyle\times\cos\omega t,$ (4)

up to constant terms, where $\vec{\sigma}_{e}$ and $\vec{\sigma}_{n}$ are
Pauli matrices $\vec{\sigma}$ acting on the corresponding qubit
$\ket{e}\in{\mathrm{span}\\{\ket{S^{z}=-,0}\\}}$ and
$\ket{N}\in{\mathrm{span}\\{\ket{I^{z}=0,+}\\}}$, respectively. In the vector
space with the basis choice of $\\{\ket{-0},\ket{-+},\ket{00},\ket{0+}\\}$, we
obtain the matrix representation

$\displaystyle
H_{\mathrm{eff}}(t)=\frac{1}{\sqrt{2}}\begin{pmatrix}0&\gamma_{N}B^{-}&\gamma_{e}B^{-}&0\\\
\gamma_{N}B^{+}&0&0&\gamma_{e}B^{-}\\\ \gamma_{e}B^{+}&0&0&\gamma_{N}B^{-}\\\
0&\gamma_{e}B^{+}&\gamma_{N}B^{+}&0\end{pmatrix}\cos\omega t,$ (5)

with $B^{\pm}\equiv B^{x}\pm iB^{y}$.

It is convenient to work in the interaction picture with $H_{\parallel}$
treated as the unperturbed Hamiltonian. The effective Hamiltonian is then
given by $\tilde{H}_{\mathrm{eff}}(t)\equiv
e^{iH_{\parallel}t}H_{\mathrm{eff}}(t)e^{-iH_{\parallel}t}$ or

$\displaystyle\tilde{H}_{\mathrm{eff}}(t)=\frac{1}{\sqrt{2}}\begin{pmatrix}0&\gamma_{N}B^{-}e^{i(Q+A_{\parallel})t}&\gamma_{e}B^{-}e^{i\Delta
t}&0\\\
\gamma_{N}B^{+}e^{-i(Q+A_{\parallel})t}&0&0&\gamma_{e}B^{-}e^{i(\Delta-
A_{\parallel})t}\\\ \gamma_{e}B^{+}e^{-i\Delta
t}&0&0&\gamma_{N}B^{-}e^{iQt}\\\ 0&\gamma_{e}B^{+}e^{-i(\Delta-
A_{\parallel})t}&\gamma_{N}B^{+}e^{-iQt}&0\end{pmatrix}\cos\omega t.$ (6)

By noting that the fast oscillation terms in the above expression can be
neglected, it becomes clear the oscillating magnetic field drives
transformation between two energy levels whose energy gap is close to the
oscillation frequency $\omega$. For example, if we start from
$\ket{\psi(t=0)}=\ket{0+}$ and choose $\omega=\Delta-A_{\parallel}$, which is
typically in the microwave frequency range, the dynamics of the quantum state
are given by

$\displaystyle\ket{\psi(t)}=\exp\left(\frac{i}{\sqrt{2}}\gamma_{e}\vec{B}_{\perp}\cdot\vec{\sigma}t\right)\begin{pmatrix}\ket{-+}\\\
\ket{0+}\end{pmatrix}.$ (7)

Thus, a manipulation of the electron spin state that only affects states with
$\ket{I^{z}=+}$ is possible. Similar dynamics controlled on $\ket{I^{z}=0}$
are realized with the choice $\omega=\Delta$. If we choose $\omega=Q$ instead,
which is typically in the radio frequency range, the dynamics are expressed as

$\displaystyle\ket{\psi(t)}=\exp\left(\frac{i}{\sqrt{2}}\gamma_{N}\vec{B}_{\perp}\cdot\vec{\sigma}t\right)\begin{pmatrix}\ket{00}\\\
\ket{0+}\end{pmatrix}.$ (8)

Thus, a manipulation of the nuclear spin state that only affects states with
$\ket{S^{z}=0}$ is possible also. Similar dynamics controlled on the
$\ket{S^{z}=-}$ state are realized with the choice $\omega=Q+A_{\parallel}$.
Note that energy gaps relevant to the neglected five energy eigenstates
composed of $\ket{S^{z}=+}$ and/or $\ket{I^{z}=-}$ take different values under
non-zero $B_{0}$, so we can stick to the effective two-qubit system of the
total Hilbert space by simply restricting ourselves to the four relevant
frequencies, i.e. $\Delta-A_{\parallel}$, $\Delta$, $Q$, and
$Q+A_{\parallel}$.

The dynamics described in Eq. 7 (Eq. 8) with $\vec{B}_{\perp}\propto\hat{x}$
and $\hat{y}$ represent the (controlled-)$R_{x}$ and $R_{y}$ gates acting on
the qubit $\ket{e}$ ($\ket{n}$), respectively, with a tunable rotation angle
$\theta=\sqrt{2}\gamma_{e}B_{\perp}t$ ($\theta=\sqrt{2}\gamma_{n}B_{\perp}t$).
As is well known, by combining $R_{x}$ and $R_{y}$ gate operations one can
construct an arbitrary $SU(2)$ operation acting on the target qubit. By also
noting that the controlled-$R_{x}$ gate with $\theta=\pi$ (or simply
$\pi_{x}$) works as the CNOT gate up to a global phase factor, an arbitrary
$SU(4)$ operation on the two-qubit system is possible [36]. Finally,
projection measurement of the nuclear spin qubit $\ket{n}$ can be performed by
combining the CNOT gate acting on the electron spin and the preceding
fluorescence measurement [37]. Physically, this final CNOT gate operation is
done with a $\pi$-pulse with frequency $\Delta$ so that the electron spin
state is excited only when the nuclear spin state is $\left|0\right>$. In the
qubit picture, signal strength of the fluorescence measurement is
characterized by

$\displaystyle{F}\equiv\frac{1}{2}\Braket{\varphi}{\sigma^{z}}{\varphi},$ (9)

where $\ket{\varphi}$ is the qubit state of the electron or nuclear spin
depending on the setup.

Thanks to the available spin state manipulation and measurement described so
far, the NV center works as a multimodal quantum sensor [34]. Additionally, we
can operate with either a single NV center or an ensemble of NV centers. In
this paper, we focus on the latter choice with which a large number of NV
centers, $N\gg 1$, helps improve the sensitivity by accumulation of large
statistics.

### II.2. Dc magnetometry

Now, we describe the so-called Ramsey sequence [38] used for dc magnetometry.
Throughout this and the next subsections, we focus on the evolution of a qubit
state $\ket{\varphi(t)}$, which can be either the electron or nuclear spin
state. In the matrix representation, we use the basis
$\\{\ket{S^{z}}=\ket{-},\ket{0}\\}$ for the electron spin and
$\\{\ket{I^{z}}=\ket{0},\ket{+}\\}$ for the nuclear spin.

The Ramsey sequence is sensitive to a dc-like magnetic field $B_{s}(t)\hat{z}$
along the $z$-axis.222 The magnetic field in the $xy$-plane can be neglected
as long as its oscillation frequency is far from the energy gap of the qubit
system. See the calculation of the Rabi cycle.  Let
$\tilde{H}_{\mathrm{int}}(t)$ be the corresponding interaction Hamiltonian in
the interaction picture defined as

$\displaystyle\tilde{H}_{\mathrm{int}}(t)=\frac{1}{2}\gamma
B_{s}(t)\sigma^{z},$ (10)

where $\gamma=\gamma_{e}$ or $\gamma_{N}$ is the suitable choice of the
gyromagnetic ratio. Starting from the lower level
$\ket{\varphi(0)}=(0,1)^{\intercal}$, the qubit state evolution under the
Ramsey sequence is given by333 In this expression and the later discussion, we
neglect the time spent on gate operations for simplicity. It is a reasonable
approximation when the amplitude of the magnetic pulse used for spin
operations is large enough as can be seen from Eq. 8.

$\displaystyle\ket{\varphi(\tau)}=R_{x}^{\pi/2}\exp\left(-i\int_{0}^{\tau}dt\,\tilde{H}_{\mathrm{int}}(t)\right)R_{y}^{\pi/2}\begin{pmatrix}0\\\
1\end{pmatrix},$ (11)

where $\tau$ is the time duration of free precession, while
$R_{\alpha}^{\theta}$ $(\alpha=x,y)$ denotes the corresponding $R_{\alpha}$
gate operation represented in matrices as

$\displaystyle R_{x}^{\theta}$
$\displaystyle=\begin{pmatrix}\cos\frac{\theta}{2}&-i\sin\frac{\theta}{2}\\\
-i\sin\frac{\theta}{2}&\cos\frac{\theta}{2}\end{pmatrix},$ (12) $\displaystyle
R_{y}^{\theta}$
$\displaystyle=\begin{pmatrix}\cos\frac{\theta}{2}&-\sin\frac{\theta}{2}\\\
\sin\frac{\theta}{2}&\cos\frac{\theta}{2}\end{pmatrix}.$ (13)

If the signal magnetic field oscillates slowly as
$B_{s}(t)=B_{s}^{0}\cos(\epsilon t+{\phi})$, the signal strength of the
fluorescence measurement for the state $\ket{\varphi(t)}$ is explicitly
calculated as

$\displaystyle{F}\simeq\frac{\gamma
B_{s}^{0}}{2\epsilon}\left[\sin(\epsilon\tau+{\phi})-\sin{\phi}\right],$ (14)

under the assumption of ${F}\ll 1$. It takes a constant value ${F}\to\gamma
B_{s}^{0}\tau\cos{\phi}/2$ at $\epsilon\to 0$, while the cancellation of fast
oscillations leads to an asymptotic behavior ${F}\propto\epsilon^{-1}$ when
$\epsilon\tau\gtrsim 1$. Thus, this approach is effective for a dc-like signal
with an angular frequency $\epsilon\ll 1/\tau$.

Figure 2: The protocol of dc magnetometry using nuclear spins.

When considering an ordinary magnetic field, electron spins are more useful
than nuclear spins to obtain a sizable effect within a fixed time duration
$\tau$ due to the hierarchy $\gamma_{e}\gg\gamma_{N}$. However, this is not
the case for the axion dark matter detection, because, as we will see below,
the axion interaction strength with electron and nuclear spins have completely
different, and model dependent, relationships. Therefore, it is worthwhile to
perform the dc magnetometry experiment based on nuclear spins in addition to
the one based on electron spins [29]. In Fig. 2, we show the protocol for dc
magnetometry using nuclear spins. Both of the qubits should be initialized to
$(0,1)^{\intercal}$ through laser cooling and an appropriate operation of the
CNOT gates before starting the protocol, preparing the $\ket{00}$ state.
$\varphi_{\tau}\simeq{F}$ is the relative phase factor generated during the
free precession.

We have not taken account of the effects of relaxation in the above
expression. There are two different relaxation time scales for each spin
species, the longitudinal relaxation time $T_{1}$ and the transverse
relaxation, or dephasing, time $T_{2}^{*}$. $T_{1}$ characterizes the spin
flip associated with the energy transfer to or from the environment, which
takes $T_{1e}\sim$6\text{\,}\mathrm{m}\mathrm{s}$$ [39, 40] and
$T_{1N}\sim$4\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$$ [41] for the electron
and nuclear spins, respectively, at room temperature. The dominant source of
the transverse relaxation, on the other hand, is dephasing of spins due to the
inhomogeneous dc magnetic field caused by, e.g. nuclear spins or lattice
defects. $T_{2e}^{*}\sim$1\text{\,}\mu\mathrm{s}$$ [34, 42] and
$T_{2N}^{*}\sim$7.25\text{\,}\mathrm{m}\mathrm{s}$$ [43]444 Although the
measured value of $T_{2N}^{*}\sim$7.25\text{\,}\mathrm{m}\mathrm{s}$$ is for a
single NV center, we use this value as a reasonable estimate of $T_{2N}^{*}$
for an ensemble of NV centers, since for a
$\sim$1\text{\,}\mathrm{p}\mathrm{p}\mathrm{m}$$ concentration of NV centers,
$T_{1e}$ is still in the order of milliseconds [44].  are measured at room
temperature. The large hierarchy $T_{2N}^{*}/T_{2e}^{*}\sim 10^{4}$ is
consistent with the large hierarchy of interaction strengths
$\gamma_{e}/\gamma_{N}\sim 10^{4}$. It should be noted that $T_{2N}^{*}$ is
naturally bounded by $T_{1e}$ since random electron spin flips induce
dephasing of nuclear spins through the hyperfine interaction. On the other
hand, $T_{2e}^{*}$ does not necessarily limit $T_{2N}^{*}$. For example, the
dephasing time scale of the quantum states $(\ket{0+}+\ket{00})/\sqrt{2}$
would be $T_{2N}^{*}$. Our protocol for nuclear spin dc magnetometry
corresponds to this case.

Taking into account the effects of relaxation, the signal strength in Eq. 14
is rescaled as ${F}\to{F}e^{-\tau/T_{2N}^{*}}$. Accordingly, the optimistic
choice of $\tau$ to maximize the sensitivity turns out to be $\tau\sim
T_{2N}^{*}/2$ [45]. We will use this choice for later analysis.

### II.3. Ac magnetometry

As we have seen thus far, the Ramsey sequence is not effective for an ac
magnetic field with an angular frequency $\epsilon\gtrsim 1/\tau$. To realize
another approach sensitive to such high-frequency signals, we can make use of
the Hahn-echo sequence [46] or dynamical decoupling sequences [47, 48] in more
general context. The time evolution of a qubit state under the Hahn echo
sequence is described by

$\displaystyle\ket{\varphi(\tau)}=R_{x}^{\pi/2}\exp\left(-i\int_{\tau/2}^{\tau}dt\,\tilde{H}_{\mathrm{int}}(t)\right)R_{y}^{\pi}\exp\left(-i\int_{0}^{\tau/2}dt\,\tilde{H}_{\mathrm{int}}(t)\right)R_{y}^{\pi/2}\begin{pmatrix}0\\\
1\end{pmatrix},$ (15)

where the only difference from the Ramsey sequence is the $\pi_{y}$ operation
at the middle of the free precession. This operation reverses the effect from
the signal magnetic field and achieves constructive interference of the
oscillating signal effect before and after the $\pi_{y}$ operation when the
angular frequency is $\sim 2\pi/\tau$. The signal strength is explicitly
calculated as follows:

$\displaystyle{F}=\frac{2\gamma
B_{s}^{0}}{\epsilon}\sin^{2}\frac{\epsilon\tau}{4}\sin\left(\frac{\epsilon\tau}{2}+{\phi}\right),$
(16)

which is suppressed in both the dc limit $\epsilon\to 0$ as $S\propto\epsilon$
and the high frequency limit $\epsilon\tau\gg 1$ as $S\propto\epsilon^{-1}$.
On the other hand, it peaks at $\epsilon=2\pi/\tau$ with a peak height
$|S|=(\gamma B_{s}^{0}\sin{\phi})/\pi$ and a peak width $\Delta\epsilon\sim
1/\tau$. This confirms that the Hahn echo sequence has a narrow-band
sensitivity to the angular frequency around $\epsilon\sim 2\pi/\tau$. Fig. 3
shows the protocol of the Hahn echo sequence applied to the nuclear spin qubit
$\ket{n}$.

Figure 3:  The protocol of ac magnetometry using nuclear spins. $\varphi_{1}$
and $\varphi_{2}$ represent the phase factors acquired in the first half
($0<t<\tau/2$) and the second half ($\tau/2<t<\tau$) of the free precession
time, respectively.

In the Hahn echo sequence, the relevant transverse relaxation time $T_{2}$
tends to be longer than $T_{2}^{*}$ for the Ramsey sequence because any dc
magnetic noise effect cancels out due to the $\pi_{y}$ operation. The dominant
contribution to $T_{2}$ comes from the decoherence effect caused by dipole-
dipole self-interaction among spins. The observed value of
$T_{2e}\sim$100\text{\,}\mu\mathrm{s}$$ [49] at room temperature shows a two
orders of magnitude enhancement compared with
$T_{2e}^{*}\sim$1\text{\,}\mu\mathrm{s}$$. Conversely, the observed value of
$T_{2N}\sim$10\text{\,}\mathrm{m}\mathrm{s}$$ at room temperature is
comparable to $T_{2N}^{*}$ because both of them are limited by the single
parameter $T_{1e}$. However, both $T_{1}$ and $T_{2}$ can be further extended
with more sophisticated setups. One possibility is to consider a cryogenic
environment; for example, $T_{1e}\sim 100$ s is reported at
$\lesssim$50\text{\,}\mathrm{K}$$ [39], where
$T_{2N}\sim$100\text{\,}\mathrm{s}$,T_{2N}^{*}\sim$1\text{\,}\mathrm{s}$$
seems to be a reasonable assumption. Another possibility is to perform a
dynamic decoupling (DD) sequence with a large number ($N_{\pi}$) of $\pi_{y}$
pulses during the free precession, which also contributes to making $T_{2}$
longer. In this case, we obtain the signal strength as

$\displaystyle{F}=$ $\displaystyle\frac{\gamma
B_{s}^{0}}{\epsilon}\sin\frac{\epsilon\tau}{2}\sin\left(\frac{\epsilon\tau}{2}+{\phi}\right)\tan\frac{\epsilon\tau}{2(N_{\pi}+1)},$
(17)

which recovers Eq. 16 for $N_{\pi}=1$. Under a DD sequence with $N_{\pi}\gg
1$, we have not only a longer $T_{2N}$, but also a sensitivity peak located at
a higher angular frequency $N_{\pi}/\tau$.

Similar to the dc magnetometry case, the relaxation effect rescales the signal
strength Eq. 16 as ${F}\to{F}e^{-\tau/T_{2N}}$, and the optimal choice of
$\tau$ turns out to be $\tau=T_{2N}/2$.

## III. Axion dark matter

### III.1. Setup

The axion can account for the total relic abundance of DM through mechanisms
such as the misalignment mechanism [50, 51, 52] or production from topological
defects (see Refs. [53, 54, 55] for reviews). A wide range of the axion mass
$m_{a}$ could be consistent with the DM relic abundance, as small as
$m_{a}\sim${10}^{-22}\text{\,}\mathrm{e}\mathrm{V}$$, below which the model is
inconsistent with the existence of DM-dominated dwarf galaxies [56], and as
large as $m_{a}\gtrsim O(1)\,$\mathrm{e}\mathrm{V}$$, where cosmological and
astrophysical constraints on the axion DM tend to become severe (see, for
example, plots in [57]). The axion is described by a classical field
experiencing coherent oscillation

$\displaystyle
a(t,\vec{x})=a_{0}\sin(m_{a}t+m_{a}\vec{v}_{a}\cdot\vec{x}+{\phi}),$ (18)

where $\vec{v}_{a}$ is the axion velocity. Considering the energy density
stored in the coherent oscillation, a relationship
$\rho_{a}=(m_{a}a_{0})^{2}/2$ holds, where
$\rho_{a}\sim$0.4\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$
[58] is the local energy density of DM. Note that the velocity $\vec{v}_{a}$
and the oscillation phase ${\phi}$ are constant only within the de Broglie
wavelength $\sim(m_{a}v_{a})^{-1}$. Thus, $\vec{v}_{a}$ and ${\phi}$ observed
at the laboratory vary with the time scale of
$\tau_{a}\sim(m_{a}v_{a}^{2})^{-1}$, which is called the coherence time.
Assuming that the axion DM halo is virialized, its typical velocity is
estimated as $v_{a}\sim 10^{-3}$, leading to the order estimate of

$\displaystyle\tau_{a}\sim$6.6\text{\,}\mathrm{s}$\left(\frac{${10}^{-10}\text{\,}\mathrm{e}\mathrm{V}$}{m_{a}}\right).$
(19)

Due to its pseudoscalar nature, the axion generally interacts with the SM
fermions $\psi_{\chi}$ of the form

$\displaystyle\mathcal{L}_{\mathrm{int}}=\sum_{\chi}\frac{C_{\chi}}{2f_{a}}\left(\partial_{\mu}a\right)\overline{\psi}_{\chi}\gamma^{\mu}\gamma_{5}\psi_{\chi},$
(20)

where $f_{a}$ is the axion decay constant and $C_{\chi}$ are the model-
dependent coefficients. The index ${\chi}$ labels the SM fermions, including
electron $e$, neutron $n$, and proton $p$.555 In the KSVZ axion model [13, 14]
we have $|C_{e}|\ll 1$ while $C_{p}\sim C_{n}\sim O(1)$. In the DFSZ axion
model [15, 16] or the flavorful axion model [59, 60], we have $C_{e}\sim
C_{n}\sim C_{p}\sim O(1)$.  In the non-relativistic limit, this interaction
term describes the axion interaction with fermion spins $\vec{S}_{\chi}$ given
by

$\displaystyle
H_{\mathrm{int}}=\sum_{\chi}\frac{g_{a{\chi}{\chi}}}{m_{\chi}}\vec{\nabla}a\cdot\vec{S}_{\chi},$
(21)

where the dimensionless coupling constants $g_{a{\chi}{\chi}}\equiv
C_{\chi}m_{\chi}/f_{a}$ are used. It is seen that the axion gradient
$\vec{\nabla}a$ may be regarded as an effective magnetic field, and the
interaction can be rewritten as

$\displaystyle
H_{\mathrm{int}}=\sum_{\chi}\gamma_{\chi}\vec{B}_{\chi}(t)\cdot\vec{S}_{\chi},$
(22)

with the gyromagnetic ratio $\gamma_{\chi}$. Substituting Eq. 18 into the
above expression, the fermion-dependent effective magnetic fields are
calculated as

$\displaystyle\vec{B}_{\chi}(t)\simeq\frac{g_{a{\chi}{\chi}}}{m_{\chi}\gamma_{\chi}}\sqrt{2\rho_{a}}\vec{v}_{a}\cos(m_{a}t+{\phi}),$
(23)

where higher-order terms of $v_{a}\ll 1$ are neglected. The amplitude of the
effective magnetic field $B_{\chi}^{0}$ is estimated as

$\displaystyle
B_{\chi}^{0}\sim$4\text{\,}\mathrm{a}\mathrm{T}$\times\left(\frac{g_{a{\chi}{\chi}}}{10^{-10}}\right).$
(24)

Note that, for typical axion models with $C_{e}\sim C_{n}\sim C_{p}$, this
“effective” magnetic field for nucleons $B_{n,p}^{0}$ is larger by a factor
$\sim m_{n,p}/m_{e}$ than that for the electron $B_{e}^{0}$. Taking account of
the gyromagnetic ratio for the nucleon/electron, it results in the same order
of the interaction strengths with the electron and nucleons. This is
qualitatively different from the ordinary magnetic field, which acts on the
electron spin much more strongly than on nucleons due to the difference of the
gyromagnetic ratio. In this sense, the use of nuclear spins in the NV center
is beneficial for the axion DM search.

For isolated fermion spins, the corresponding $\vec{B}_{\chi}(t)$ works just
the same as the ordinary magnetic field aside from that its amplitude and
oscillation phase vary with the time scale of $\tau_{a}$. It is the effective
magnetic field for electrons, $\vec{B}_{e}(t)$, that is sensed by the ordinary
NV center magnetometry as proposed in [29]. On the other hand, if we use
nuclear spins for magnetometry, both $\vec{B}_{n}(t)$ and $\vec{B}_{p}(t)$ can
be relevant as we will see shortly.

### III.2. Axion interaction with the ${\ce{{}^{1}4N}}$ spin

The interaction between the ${\ce{{}^{1}4N}}$ spin and an ordinary magnetic
field is characterized by its gyromagnetic ratio $\gamma_{N}$. As a rare
stable odd-odd nucleus, $\gamma_{N}$ has contributions from both the neutron
and proton spins and the orbital angular momentum of the proton. However,
$\gamma_{N}$ does not reflect the axion interaction strength with the
${\ce{{}^{1}4N}}$ spin, which is determined by the axion interaction with
neutron and proton spins, $g_{ann}$ and $g_{app}$. To accurately describe the
axion-${\ce{{}^{1}4N}}$ interaction, we need to understand the composition of
the ${\ce{{}^{1}4N}}$ spin $I=1$. In this subsection, we discuss this issue
under the assumption that the nuclear shell model well describes internal
structure of the ${\ce{{}^{1}4N}}$ nucleus.

${\ce{{}^{1}4N}}$ has seven neutrons and seven protons. According to the
nuclear shell model, both nucleons occupy the orbitals as
$1s_{1/2}^{2}1p_{3/2}^{4}1p_{1/2}$, where both $1s_{1/2}$ and $1p_{3/2}$
orbitals form closed shells. Thus, the nuclear spin $I=1$ comes from the
synthesis of the total neutron spin $J_{n}=1/2$ of a neutron in the $1p_{1/2}$
orbital and $J_{p}=1/2$ of a proton in the proton counterpart of the orbital.
In the representation theory of $SU(2)$, this corresponds to the decomposition
$\frac{1}{2}\otimes\frac{1}{2}=1\oplus 0$ with the first term in the right-
hand side is selected, while each spin-$\frac{1}{2}$ representation on the
left-hand side comes from the decomposition $\frac{1}{2}\otimes
1=\frac{3}{2}\oplus\frac{1}{2}$.

Let us explicitly write down the $I=1$ states in terms of the eigenstates of
the spin $\vec{S}_{\chi}$ and the orbital angular momentum $\vec{L}_{\chi}$ of
nucleons ${\chi}=n,p$ in the $1p$ orbitals. Let $\ket{\uparrow}_{\chi}$ and
$\ket{\downarrow}_{\chi}$ be the spin up and down states, and $\ket{m}_{\chi}$
($m=-,0,+$) be the eigenstates of $L_{\chi}^{z}$ for each nucleon $f$.
Corresponding to this decomposition $\frac{1}{2}\otimes
1=\frac{3}{2}\oplus\frac{1}{2}$, the $J_{\chi}=\frac{1}{2}$ component (i.e.,
the $1p_{1/2}$ orbital for the nucleon ${\chi}$) is given by

$\displaystyle\begin{pmatrix}\ket{u_{\chi}}\\\
\ket{d_{\chi}}\end{pmatrix}\equiv\frac{1}{\sqrt{3}}\begin{pmatrix}\ket{\uparrow}_{\chi}\ket{0}_{\chi}-\sqrt{2}\ket{\downarrow}_{\chi}\ket{+}_{\chi}\\\
\sqrt{2}\ket{\uparrow}_{\chi}\ket{-}_{\chi}-\ket{\downarrow}_{\chi}\ket{0}_{\chi}\end{pmatrix}.$
(25)

Using these expressions to evaluate the second decomposition
$\frac{1}{2}\otimes\frac{1}{2}=1\oplus 0$, the nuclear spin $I=1$ component is
expressed as

$\displaystyle\vec{\psi}_{I}^{\intercal}\equiv\begin{pmatrix}\ket{u_{p}}\ket{u_{n}}\\\
\dfrac{1}{\sqrt{2}}\pqty{\ket{u_{p}}\ket{d_{n}}+\ket{d_{p}}\ket{u_{n}}}\\\
\ket{d_{p}}\ket{d_{n}}\end{pmatrix}.$ (26)

To go further, we calculate the matrix elements of the spin operators
$\vec{S}_{\chi}$ (${\chi}=n,p$) in the $36$-dimensional space corresponding to
the all possible choices of $S_{\chi}^{z}$ and $L_{\chi}^{z}$ (${\chi}=n,p$).
In particular, since the nuclear spin $I=1$ states corresponds to the three
dimensional subspace spanned by three vectors in $\vec{\psi}_{I}$, matrix
elements of the spin operators $\vec{S}_{\chi}$ (${\chi}=n,p$) for these basis
vectors represent how the axion interacts with the ${\ce{{}^{1}4N}}$ spin.
From a straightforward calculation, we obtain the following

$\displaystyle\vec{\psi}_{I}^{\dagger}S^{\alpha}\vec{\psi}_{I}=-\frac{1}{6}I^{\alpha},$
(27)

where $I^{\alpha}$ $(\alpha=x,y,z)$ are the spin-1 representations of the
$SU(2)$ generators. This result is consistent with the treatment in [61]. Note
that the spin operators also have non-zero matrix elements outside the three
dimensional subspace, a part of which connects different spin states. These
interactions can in principle invoke the transition from the ground state with
$I=1$ to, e.g. an excited state with $I=0$. However, since the relevant energy
scale of $O(1$–$10)\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ [62] is far beyond the
current setup, we can safely neglect these terms and focus on the terms in Eq.
27 that preserve the nuclear spin structure.

Since the right-hand side of Eq. 27 is proportional to the nuclear spin
operators $I^{\alpha}$, $\vec{S}_{\chi}$ (${\chi}=n,p$) effectively acts as
the $I=1$ spin operators with a non-trivial coefficient $-1/6$. In App. A, we
provide proof that the same interpretation is possible whenever the spin
$S=1/2$ and a general orbital angular momentum $L=\ell$ are considered, and
derive a systematic way to calculate the coefficient. By substituting Eq. 27
in Eq. 21, we obtain an effective axion-${\ce{{}^{1}4N}}$ interaction
Hamiltonian

$\displaystyle H_{\mathrm{int}}=\gamma_{N}\vec{B}_{N}(t)\cdot\vec{I},$ (28)

with the effective magnetic field defined as

$\displaystyle\vec{B}_{N}(t)$ $\displaystyle\equiv
B_{N}\hat{v}_{a}\cos(m_{a}t+{\phi}),$ (29) $\displaystyle\gamma_{N}B_{N}$
$\displaystyle\simeq-\frac{1}{6}\left(\frac{g_{app}}{m_{p}}+\frac{g_{ann}}{m_{n}}\right)\sqrt{2\rho_{a}}v_{a},$
(30)

with $\hat{v}_{a}\equiv\vec{v}_{a}/v_{a}$. For convenience, we define

$\displaystyle\tilde{f}_{a}\equiv\left|\frac{g_{app}}{2m_{p}}+\frac{g_{ann}}{2m_{n}}\right|^{-1},$
(31)

with which $B_{N}\propto\tilde{f}_{a}^{-1}$. Since we can rewrite it as
$\tilde{f}_{a}=2f_{a}/(C_{p}+C_{n})$, $\tilde{f}_{a}$ is of the same order as
$f_{a}$ if coefficients $C_{p}$ and $C_{n}$ are of $O(1)$. A fascinating
consequence of the ${\ce{{}^{1}4N}}$ spin as an odd-odd nucleus is that the
axion coupling is proportional to the combination Eq. 31, and that the
coupling strength is sensitive to the relative sign of $g_{app}$ and
$g_{ann}$.

## IV. Sensitivity estimation

When measurements are repeated ${N_{\mathrm{obs}}}$ times, we obtain time-
sequence data labeled by $j=1,\dots,{N_{\mathrm{obs}}}$, representing the
measurement starting at time $t_{j}$. For simplicity, we assume
$t_{j}=(j-1)\tau$ with $\tau=T_{2n}^{*}/2$ ($T_{2n}/2$) for the dc (ac)
effective magnetometry approach,666 In this context, we neglect the
measurement overhead, including the state preparation and fluorescence
measurement. This is a reasonable approximation given that $\tau\sim
O(1)\,$\mathrm{m}\mathrm{s}$$ while the overhead time is typically of
$O(10$–$100)\,$\mu\mathrm{s}$$ [63].  though it is not necessary for the
following discussion that the measurements are repeated with equal time
intervals. As a result, ${t_{\mathrm{obs}}}\equiv{N_{\mathrm{obs}}}\tau$
denotes the total observation time. Let $\rho_{j}$ be the density matrix
representing the quantum state of the NV center ensemble before the $j$-th
fluorescence measurement. Since our observable is defined as an operator

$\displaystyle
M_{j}^{z}\equiv\frac{1}{2N}\sum_{\ell=1}^{N}\sigma_{j\ell}^{z},$ (32)

where $\sigma_{j\ell}^{z}$ acts on the qubit state of the nuclear spin in the
$\ell$-th NV center at the $j$-th measurement, the data obtained at $t_{j}$
can be calculated as
$\Braket{M_{j}^{z}}_{\rho_{j}}\equiv{\mathrm{Tr}\,\left[{\rho_{j}M_{j}^{z}}\right]}$.
It should be noted that $\Braket{M_{j}^{z}}_{\rho_{j}}$ represents the
$N$-qubit average of the signal at time $t_{j}$, which asymptotes to the
signal strength $F$ with an appropriate choice of the phase factor in the
limit of $N\to\infty$.

The expression above is useful for calculating the ensemble average over
distributions of the axion parameters. It should be noted that $\rho_{j}$ is
equivalent to Eq. 11 and to Eq. 15 with the replacement ${\phi}\to
m_{a}t_{j}+{\phi}$ for the Ramsey and the Hahn echo sequences, respectively.
Therefore, it depends on the axion velocity $\vec{v}_{a}$ and the phase factor
${\phi}$ through the expression of the effective magnetic field Eq. 30. If we
neglect the Earth’s relative motion against the Galactic center, ${\phi}$ is
uniformly distributed in the range $[0,2\pi)$, while the axion velocity has a
random direction with typically a value of $v_{a}\sim 10^{-3}$. Using these
distributions, for example, the ensemble average of the observation result
$\mathcal{M}_{j}$ is calculated as

$\displaystyle\Braket{\mathcal{M}_{j}}\equiv\left.\frac{1}{2\pi}\int
d{\phi}\,\frac{1}{4\pi}\int
d\hat{v}_{a}\,{\mathrm{Tr}\,\left[{\rho_{j}M_{j}^{z}}\right]}\right|_{v_{a}=10^{-3}},$
(33)

where we do not take into account the distribution of $v_{a}$, which highly
depends on the model of the DM profile in our galaxy and only results in an
$O(1)$ modification. It should be noted that, here and hereafter, we neglect
the subscript of the ensemble-averaged quantities for notational simplicity.
As anticipated, the randomness of the signal direction and phase causes
cancellation of the averaged signal, $\Braket{M_{j}}=0$.

To derive meaningful insights from the data, we use two-point functions of the
time-sequence data defined as [64]

$\displaystyle
C_{jj^{\prime}}\equiv\begin{cases}{\mathrm{Tr}\,\left[{(\rho_{j}\otimes\rho_{j^{\prime}})(M_{j}^{z}\otimes
M_{j^{\prime}}^{z})}\right]},&(j\neq j^{\prime})\\\
{\mathrm{Tr}\,\left[{\rho_{j}M_{j}^{z}M_{j}^{z}}\right]}.&(j=j^{\prime})\end{cases}$
(34)

Since the coherence of the axion signal is maintained only for the duration
$\tau_{a}$, $C_{jj^{\prime}}$ behaves differently for
$|t_{j}-t_{j^{\prime}}|<\tau_{a}$ and $|t_{j}-t_{j^{\prime}}|>\tau_{a}$. A
combined expression can be given as

$\displaystyle C_{jj^{\prime}}=$ $\displaystyle\left.\frac{1}{2\pi}\int
d{\phi}\,\frac{1}{4\pi}\int d\hat{v}_{a}\,\frac{1}{2\pi}\int
d{\phi}^{\prime}\,\frac{1}{4\pi}\int
d\hat{v}_{a}^{\prime}\,{\mathrm{Tr}\,\left[{\rho_{jj^{\prime}}M_{jj^{\prime}}^{z}}\right]}\right|_{v_{a}=10^{-3}}$
$\displaystyle\times\left[\Theta(|t_{j}-t_{j^{\prime}}|-\tau_{a})+8\pi^{2}\delta({\phi}-{\phi}^{\prime})\delta(\hat{v}_{a}-\hat{v}_{a}^{\prime})\Theta(\tau_{a}-|t_{j}-t_{j^{\prime}}|)\right],$
(35)

where $\rho_{jj^{\prime}}\equiv\rho_{j}\otimes\rho_{j^{\prime}}$ and
$M_{jj^{\prime}}^{z}\equiv M_{j}^{z}\otimes M_{j^{\prime}}^{z}$ for $j\neq
j^{\prime}$ and $\rho_{jj}\equiv\rho_{j}$ and $M_{jj}^{z}\equiv
M_{j}^{z}M_{j}^{z}$. Also, $\Theta$ is the Heaviside step function. The
integral variables $({\phi},\hat{v}_{a})$ and
$({\phi}^{\prime},\hat{v}_{a}^{\prime})$ correspond to the axion parameters at
time $t_{j}$ and $t_{j^{\prime}}$, respectively, and the delta functions in
the second line denote the coherence of the signal for
$|t_{j}-t_{j^{\prime}}|<\tau_{a}$. In addition, we introduce the power
spectral density (PSD), which is defined as the ensemble-averaged expectation
value ${\mathcal{P}_{k}}\equiv\Braket{\mathcal{O}_{k}}$ of the operator

$\displaystyle\mathcal{O}_{k}\equiv\frac{\tau^{2}}{{t_{\mathrm{obs}}}}\sum_{j,j^{\prime}}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}M_{j}^{z}M_{j^{\prime}}^{z},$ (36)

with $k=0,\dots,{N_{\mathrm{obs}}}-1$. Each ${\mathcal{P}_{k}}$ can be
calculated through the Fourier transformation of the two-point functions as

$\displaystyle{\mathcal{P}_{k}}=\frac{\tau^{2}}{{t_{\mathrm{obs}}}}\sum_{j,j^{\prime}}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}C_{jj^{\prime}}.$ (37)

A detailed calculation of the PSD and the relevant quantum noise is given in
App. B. From Eq. 68, the signal PSD can be defined as
${\mathcal{S}_{k}}\equiv{\mathcal{P}_{k}}-\tau/(4N)$, where the constant shift
ensures that ${\mathcal{S}_{k}}$ is proportional to the axion-induced magnetic
field $B_{N}$, and is given by

$\displaystyle{\mathcal{S}_{k}}\simeq\dfrac{2\mathcal{A}}{{t_{\mathrm{obs}}}\Delta\omega_{k}^{2}}\sin^{2}\dfrac{{t_{\mathrm{obs}}}\Delta\omega_{k}}{2},$
(38)

for ${t_{\mathrm{obs}}}<\tau_{a}$, while

$\displaystyle{\mathcal{S}_{k}}\simeq\dfrac{2\mathcal{A}}{{t_{\mathrm{obs}}}\Delta\omega_{k}^{2}}\sin^{2}\dfrac{\tau_{a}\Delta\omega_{k}}{2}+\dfrac{{t_{\mathrm{obs}}}-\tau_{a}}{{t_{\mathrm{obs}}}\Delta\omega_{k}}\mathcal{A}\sin\left[\tau_{a}\Delta\omega_{k}\right],$
(39)

for ${t_{\mathrm{obs}}}>\tau_{a}$, where
$\Delta\omega_{k}\equiv\omega_{k}-m_{a}$ with $\omega_{k}\equiv 2\pi
k/{t_{\mathrm{obs}}}$ and $\mathcal{A}\propto B_{N}^{2}$ is the protocol-
dependent coefficient defined in Eq. 66. Due to the quantum noise, the
measurement result of the PSD fluctuates even without the axion DM. The
standard deviation of the PSD distribution is calculated as

$\displaystyle{\mathcal{B}_{k}}\equiv\left.\sqrt{\Braket{\mathcal{O}_{k}^{2}}-\Braket{\mathcal{O}_{k}}^{2}}\right|_{B_{N}=0}.$
(40)

As shown in Eq. 72, we obtain ${\mathcal{B}_{0}}\simeq\tau/(2\sqrt{2}N)$ and
${\mathcal{B}_{k\neq 0}}\simeq\tau/(4N)$ for our setup.

Focusing on a single bin $k$, the signal estimation uncertainty can be
evaluated through the well-known formula

$\displaystyle\delta
B_{N}^{2}=\left.\sqrt{\Braket{\mathcal{O}_{k}^{2}}-\Braket{\mathcal{O}_{k}}^{2}}\left(\frac{d\Braket{\mathcal{O}_{k}}}{dB_{N}^{2}}\right)^{-1}\right|_{B_{N}^{2}=0},$
(41)

which determines the estimation error of $B_{N}^{2}$ around the specific
choice of $B_{N}^{2}=0$, i.e. the model without the axion DM. We select
$B_{N}^{2}$ as a parameter to be estimated since $\Braket{O_{k}}$ does not
have a linear term in $B_{N}$ as shown in Eq. 68. By deforming the above
expression with the relationship $\delta B_{N}^{2}=2B_{N}\delta B_{N}$, we can
obtain the $X\sigma$-level detection limit to the axion-induced magnetic field
$B_{N}$ to be $X\delta B_{N}$ (where $X$ depends on the required confidence
level). To gather information from all bins and obtain the best achievable
detection limit with these observables, we follow [64]. Based on the Asimov
dataset [65] rather than the Monte Carlo simulation results, we compute the
test statistic

$\displaystyle
q=2\sum_{k=0}^{{N_{\mathrm{obs}}}-1}\left[\left(1-\frac{{\mathcal{B}_{k}}}{{\mathcal{S}_{k}}+{\mathcal{B}_{k}}}\right)-\ln\left(1+\frac{{\mathcal{S}_{k}}}{{\mathcal{B}_{k}}}\right)\right],$
(42)

with which the $95\,\%$ exclusion limit, which we adopt as the definition of
the detection limit of our approach, is determined by the criteria $q=-2.71$.

It is beneficial to consider two extreme setups and evaluate the scaling of
the detection limit as a function of ${t_{\mathrm{obs}}}$ and $N$. For this
purpose, we first observe that the signal PSD ${\mathcal{S}_{k}}$ has a
resonant structure peaked at $\Delta\omega_{k}=0$ or $\omega_{k}=m_{a}$. The
peak height is evaluated as

$\displaystyle{\mathcal{S}_{k}}=\begin{cases}\mathcal{A}\dfrac{{t_{\mathrm{obs}}}}{2},&({t_{\mathrm{obs}}}<\tau_{a})\\\\[10.0pt]
\mathcal{A}\tau_{a}\left(1-\dfrac{\tau_{a}}{2{t_{\mathrm{obs}}}}\right),&({t_{\mathrm{obs}}}>\tau_{a})\end{cases}$
(43)

where the resonance condition for these setups can be described as
${t_{\mathrm{obs}}}\Delta\omega_{k}\ll 1$ and $\tau_{a}\Delta\omega_{k}\ll 1$,
respectively. When ${t_{\mathrm{obs}}}\ll\tau_{a}$, the signal peak height
grows linearly with ${t_{\mathrm{obs}}}$, while only a single bin enjoys
resonance since ${t_{\mathrm{obs}}}(\omega_{k+1}-\omega_{k})=2\pi$. Then, we
can approximate the test statistic as

$\displaystyle q\simeq
2\left(1-\frac{{\mathcal{B}_{k_{0}}}}{{\mathcal{S}_{k_{0}}}+{\mathcal{B}_{k_{0}}}}\right)-\ln\left(1+\frac{{\mathcal{S}_{k_{0}}}}{{\mathcal{B}_{k_{0}}}}\right),$
(44)

with $k_{0}$ being the label of the resonance bin. Since the above expression
only depends on the ratio ${\mathcal{S}_{k_{0}}}/{\mathcal{B}_{k_{0}}}$, the
detection limit is solely determined by solving $q=-2.71$ for this ratio,
resulting in ${\mathcal{S}_{k_{0}}}/{\mathcal{B}_{k_{0}}}\simeq 8.48$. Since
${\mathcal{S}_{k_{0}}}/{\mathcal{B}_{k_{0}}}$ is proportional to
$N{t_{\mathrm{obs}}}B_{N}^{2}$, the detection limit to $B_{N}$ grows as
$N^{1/2}{t_{\mathrm{obs}}^{1/2}}$ as expected for a coherently oscillating
signal. On the other hand, when ${t_{\mathrm{obs}}}\gg\tau_{a}$, the peak
height is saturated to $\sim\mathcal{A}\tau_{a}$, but the number of bins
involved in the peak grows as ${t_{\mathrm{obs}}}/\tau_{a}$. In this limit,
our setup can be sensitive to small signals with
${\mathcal{S}_{k}}\ll{\mathcal{B}_{k}}$, where we can expand the expression of
the test statistic as

$\displaystyle
q\simeq-\sum_{k}\frac{{\mathcal{S}^{2}_{k}}}{{\mathcal{B}^{2}_{k}}}.$ (45)

Since the number of terms with dominant contributions grow as
${t_{\mathrm{obs}}}$ and the fraction in the summation is proportional to
$N^{2}B_{N}^{4}$, we obtain the detection limit scaling $\propto
N^{1/2}{t_{\mathrm{obs}}^{1/4}}$. Again, this scaling behavior is common for
the signal with randomized direction and phase.

To summarize, the sensitivity to the axion coupling is roughly estimated from

$\displaystyle\sqrt{\mathcal{A}}\sim\begin{cases}\left(\dfrac{\tau}{N{t_{\mathrm{obs}}}}\right)^{1/2},&({t_{\mathrm{obs}}}<\tau_{a})\\\
\left(\dfrac{\tau}{N\tau_{a}}\right)^{1/2}\left(\dfrac{\tau_{a}}{{t_{\mathrm{obs}}}}\right)^{1/4},&({t_{\mathrm{obs}}}>\tau_{a})\end{cases}$
(46)

where $\mathcal{A}$ is defined in Eq. 66.

The analysis explained so far uses the full data set with
$j=1,\dots,{N_{\mathrm{obs}}}$ and their Fourier transformation to look for a
signal. However, this should be interpreted as a way to estimate the best
achievable detection limit curves. In realistic experimental setups, on the
other hand, there are several challenges to performing such an analysis
including memory constraints and limitations on computational power. Given
these constraints and limitations, an alternative analysis procedure is the
one based on the standard deviation [66]. We can show that, by setting an
appropriate data collection time duration, the sensitivities of this procedure
to the target frequencies have the same scaling behavior with
${t_{\mathrm{obs}}}$ and $N$ as shown above.

## V. Results

Figure 4:  The $95\,\%$ exclusion limits on $\tilde{f}_{a}$ as a function of
$m_{a}$ for $T_{2}^{*}=2\tau=$7.25\text{\,}\mathrm{m}\mathrm{s}$$. The total
detector volumes of $(N,{t_{\mathrm{obs}}})=(10^{12},$1\text{\,}\mathrm{s}$)$
(magenta), $(10^{12},$1\text{\,}\mathrm{y}\mathrm{r}$)$ (green), and
$(10^{20},$1\text{\,}\mathrm{y}\mathrm{r}$)$ (cyan) are assumed. The black
dash-dotted line represents the combination of the current best constraints on
$|g_{app}|$ and $|g_{ann}|$, including constraints on $|g_{ann}|$ from neutron
star cooling [67], $\ce{K}$–$\ce{{}^{3}He}$ comagnetometer [68], and ChangE
[69] and ChangE NMR [70] experiments, and a constraint on $|g_{app}|$ from
SN1987A [71]. The black dotted line represents the prospect of constraints,
including constraints on $|g_{ann}|$ from future comagnetometers [72], the
electrostatic storage ring [73], the CASPEr-gradient experiment [74], and the
homogeneous precession domain of the superfluid $\ce{{}^{3}He}$ [75], and
constraints on $|g_{app}|$ from the proton storage ring [76], the CASPEr-
gradient experiment [74], and the nuclear magnon in $\ce{MnCO3}$ [77]. The
limit data is adopted from [57].

Figure 5:  The top (bottom) panel shows the limits on $|g_{ann}|$
($|g_{app}|$) from the Ramsey setup under an artificial assumption that
$|g_{app}|\ll|g_{ann}|$ ($|g_{ann}|\ll|g_{app}|$). We assume the dephasing
time scale of $T_{2}^{*}=$7.25\text{\,}\mathrm{m}\mathrm{s}$$. The color
conventions and the meaning of the black lines for the existing constraints
and prospects are the same as in Fig. 4.

The $95\,\%$ exclusion limits on $\tilde{f}_{a}$ defined in Eq. 31 from the
Ramsey sequence are shown in Fig. 4 with the assumed relaxation time and the
free precession time $T_{2N}^{*}=$7.25\text{\,}\mathrm{m}\mathrm{s}$$,
choosing $\tau=T_{2N}^{*}/2$.777 Precisely speaking, there are periodic $O(1)$
fluctuations of the sensitivity due to the discrete binning of the frequency
with the bin width $2\pi/{t_{\mathrm{obs}}}$. In Fig. 4, we smooth out these
fluctuations to focus on the larger-scale frequency dependence of the
sensitivity. Furthermore, the small step of the magenta line at
$\sim$1\text{\,}\mathrm{H}\mathrm{z}$$ is due to the difference between Eqs.
71 and 72.  Three colored lines correspond to the most conservative setup with
an already-achieved number of NV centers $N=10^{12}$ [63] and
${t_{\mathrm{obs}}}=$1\text{\,}\mathrm{s}$$ (magenta), the same $N=10^{12}$
but with ${t_{\mathrm{obs}}}=$1\text{\,}\mathrm{y}\mathrm{r}$$ (green), and a
rather optimistic choice of $N=10^{20}$ with
${t_{\mathrm{obs}}}=$1\text{\,}\mathrm{y}\mathrm{r}$$ (cyan). Given the
available NV concentration of
$1.6\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$ [49],
$N=10^{20}$ is achievable with a total volume of
$\sim${10}^{3}\text{\,}\mathrm{c}\mathrm{m}^{3}$$, which can be realized by a
synthesis technique for a large diamond [78] or by preparing multiple
measurement setups. Also shown by the black lines are the combined constraints
on $\tilde{f}_{a}$ from the existing experimental results (dash-dotted) and
the prospects (dotted). See [57] for details.

Since the Ramsey sequence delivers full performance when the frequency $f$
satisfies $f\lesssim 1/\tau$, the choice of
$T_{2N}^{*}=2\tau=$7.25\text{\,}\mathrm{m}\mathrm{s}$$ leads to a frequency
coverage $m_{a}/2\pi\lesssim$200\text{\,}\mathrm{H}\mathrm{z}$$
($m_{a}\lesssim 8\times 10^{-13}\,\mathrm{eV}$), outside of which the
sensitivity is rapidly lost. Another remarkable feature of our sensitivities
is the kinks of the green and cyan lines at
$m_{a}\sim$2\text{\times}{10}^{-17}\text{\,}\mathrm{e}\mathrm{V}$$ with the
corresponding axion coherence time
$\tau_{a}\sim$1\text{\,}\mathrm{y}\mathrm{r}$$. Both lines below this point
(and also the magenta line) correspond to ${t_{\mathrm{obs}}}>\tau_{a}$, so
the axion signal maintains coherence during the observation. Thus, the Ramsey
sequence has frequency-independent sensitivities for this mass range. For
higher masses, on the other hand, we need to account for a slower sensitivity
improvement $\propto(\tau_{a}{t_{\mathrm{obs}}})^{1/4}$ shown in Eq. 46, so
the detection limit plots have slopes. Finally, compared with solid lines that
show detection limit prospects, the dash-dotted part of the magenta lines,
which corresponds to the mass range $2\pi/m_{a}\ll$1\text{\,}\mathrm{s}$$,
needs special care. In this mass range, the signal strengths for all repeated
measurements are proportional to $\cos{\phi}$ with a randomly chosen phase
factor ${\phi}$; thus, it is always possible that no signal is observed
irrespective of the size of $\tilde{f}_{a}$. The magenta dash-dotted lines
should then be interpreted as a $1\sigma$ lower bound on $\tilde{f}_{a}$ when
no signal is observed.

Comparison between our results and the existing constraints or prospects shows
that our approach is promising for a broad mass range with
$m_{a}/2\pi\lesssim$200\text{\,}\mathrm{H}\mathrm{z}$$. It should be noted,
however, that the exclusion limits in Fig. 4 need to be carefully interpreted
since experimental results that lead to the orange lines have a dominant
constraint on either $|g_{ann}|$ or $|g_{app}|$, contrary to our approach
where $\tilde{f}_{a}$, a linear combination shown in Eq. 31, is directly
constrained. Due to the expression Eq. 31, both constraints on $|g_{ann}|$ and
$|g_{app}|$ in principle affect the orange lines. However, practically, either
the $|g_{ann}|$ or $|g_{app}|$ that is less constrained at a chosen $m_{a}$
determines how strongly $\tilde{f}_{a}$ is constrained at that value of
$m_{a}$. To disentangle the mixed effect of $|g_{ann}|$ and $|g_{app}|$
constraints, we also demonstrate the detection limits of our setup with
$T_{2N}^{*}=$7.25\text{\,}\mathrm{m}\mathrm{s}$$ on an individual $|g_{ann}|$
($|g_{app}|$) coupling in the top (bottom) panel of Fig. 5 under an artificial
assumption that the corresponding coupling is much larger than the other one.
Fig. 5 is useful for comparison; in particular, the optimistic setup (cyan)
shows remarkable sensitivities to $|g_{ann}|$ for
$m_{a}\lesssim$2\text{\times}{10}^{-19}\text{\,}\mathrm{e}\mathrm{V}$$
comparable to the future comagnetometer prospect [72] and to $|g_{app}|$ for
$${10}^{-15}\text{\,}\mathrm{e}\mathrm{V}$\lesssim
m_{a}\lesssim${10}^{-12}\text{\,}\mathrm{e}\mathrm{V}$$ corresponding to a gap
between the proton storage ring [76] and the CASPEr-gradient prospects [74].
However, it should be remembered that Fig. 4 is a more fundamental result of
our approach obtained without any artificial assumptions on physics
parameters.

Figure 6:  Same as Fig. 4 but with the Hahn-echo sequence. The decoherence
times of $T_{2}=$10\text{\,}\mathrm{m}\mathrm{s}$$ (top) and
$1\text{\,}\mathrm{s}$ (bottom) are assumed.

In Fig. 6, we show the $95\,\%$ exclusion limits from the Hahn-echo sequence
at room temperature with $T_{2N}=$10\text{\,}\mathrm{m}\mathrm{s}$$ (top) and
a cryogenic environment with $T_{2N}=$1\text{\,}\mathrm{s}$$ (bottom). As is
clearly shown in the plots, this approach is a narrow-band search targeted at
the frequencies
$1/\tau\sim$1\text{\,}\mathrm{H}\mathrm{z}$$–$100\text{\,}\mathrm{H}\mathrm{z}$
depending on the choice. Although the frequency coverage is limited in this
approach, the sensitivity around the target frequency is much better than the
Ramsey setup under the cryogenic environment when $T_{2N}\gg T_{2N}^{*}$ Note
that these exclusion limits can also be reinterpreted as limits on $|g_{ann}|$
and $|g_{app}|$ under certain assumptions similar to Fig. 5.

Figure 7:  Same as Fig. 4 but with the DD sequence. The decoherence time of
$T_{2}=$1\text{\,}\mathrm{s}$$ and the number of $\pi$-pulses $N_{\pi}=63$ are
assumed.

In Fig. 7, we show the $95\,\%$ exclusion limits from the DD sequence with
$T_{2}=$1\text{\,}\mathrm{s}$$ and $N_{\pi}=63$. Despite the improved
sensitivity at the peak due to the prolonged $T_{2}$, the peak width becomes
narrower for a larger $N_{\pi}$, which makes this sequence generally not
suitable for dark matter searches with unknown signal frequency.

## VI. Discussion and conclusion

We proposed a novel method to use the ${\ce{{}^{1}4N}}$ spin of NV centers in
diamond for axion dark matter searches. Our nuclear spin magnetometry
metrology approach is based on new types of protocols from Figs. 2 and 3 aimed
at dark matter searches, and provides constraints on the axion-nucleus
couplings $g_{ann}$ and $g_{app}$, which are completely independent of those
on $g_{aee}$ obtained with conventional magnetometry protocols in [29]. This
opens up a new direction for quantum sensing techniques based on NV centers
and motivates further investigation into the properties of the
${\ce{{}^{1}4N}}$ spin, including the relaxation time scales $T_{2N}^{*}$ and
$T_{2N}$, for an ensemble of NV centers under various conditions.

One of the benefits of our approach compared with other proposed ideas to
constrain $g_{ann}$ and $g_{app}$ is its broad frequency coverage arising from
the wide dynamic range of NV center magnetometry. We found several frequency
windows in which our approach has a relatively high potential within the
overall target frequency range $\lesssim$100\text{\,}\mathrm{H}\mathrm{z}$$.
Another remarkable feature is the natural sensitivity to a roughly $1:1$
linear combination of two coupling constants $g_{ann}$ and $g_{app}$ shown in
Eq. 31. Accordingly, our approach is sensitive not only to an individual
$g_{ann}$ or $g_{app}$ coupling under existence of a large hierarchy between
them but also to a relative phase between them when they have comparable
sizes, which enables us to explore the axion model after its discovery.

When preparing the diamond sample, it is possible to have the majority of the
NV centers contain the nitrogen isotope $\ce{{}^{1}5N}$ by creation via
implantation [79] or by doping during chemical vapor deposition synthesis of
diamond. In this case, we primarily obtain constraints on $g_{app}$ because
the nuclear spin of $\ce{{}^{1}5N}$ is predominantly influenced by proton
contributions. Indeed, the nuclear shell model indicates that the expression
of the axion-induced magnetic field Eq. 30 is replaced by

$\displaystyle\gamma_{N}B_{N}$
$\displaystyle\simeq-\frac{1}{3}\frac{g_{app}}{m_{p}}\sqrt{2\rho_{a}}v_{a},$
(47)

for $\ce{{}^{1}5N}$. Therefore, NV center metrology based on $\ce{{}^{1}5N}$
spins provides yet another independent piece of information about axion-
nucleon couplings.

Finally, we reflect shortly on the current state of experimental NV center
sensors with respect to our proposal. As mentioned before, the nuclear spin is
generally not used for magnetic field sensing, as it is much less sensitive
compared to the electron spin. Therefore, there is limited information about
their properties, which gives a lot of opportunities to investigate them.
Noise-wise, single-shot readout via the electron spin, which overhead is
relatively short given the long nuclear spin coherence times, makes that the
noise in the system can be rather close to the spin-projection noise [37, 80],
as investigated in this work.

###### Acknowledgements.

SC thanks Joseph Bramante for useful discussions. This work was supported by
World Premier International Research Center Initiative (WPI), MEXT, Japan.
This work was also supported by the DOE, Office of Science under contract DE-
AC02-05CH11231, partially through Quantum Information Science Enabled
Discovery (QuantISED) for High Energy Physics (KA2401032). This work was also
supported by JST Moonshot (Grant Number JPMJMS226C), JSPS KAKENHI (Grant
Number 20H05661, 23H04390, 20H05661, 24K07010), and CREST (JPMJCR23I5), JST.

## References

  * [1] F. Zwicky, “Die Rotverschiebung von extragalaktischen Nebeln,” Helvetica Physica Acta 6 (Jan., 1933) 110–127.
  * [2] S. Smith, “The Mass of the Virgo Cluster,” Astrophys. J.  83 (Jan., 1936) 23.
  * [3] V. C. Rubin and J. Ford, W. Kent, “Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions,” Astrophys. J.  159 (Feb., 1970) 379.
  * [4] K. C. Freeman, “On the Disks of Spiral and S0 Galaxies,” Astrophys. J.  160 (June, 1970) 811.
  * [5] J. P. Ostriker and P. J. E. Peebles, “A Numerical Study of the Stability of Flattened Galaxies: or, can Cold Galaxies Survive?,” Astrophys. J.  186 (Dec., 1973) 467–480.
  * [6] H. Hoekstra, H. K. C. Yee, and M. D. Gladders, “Properties of galaxy dark matter halos from weak lensing,” in Dark Matter in Galaxies, S. Ryder, D. Pisano, M. Walker, and K. Freeman, eds., vol. 220, p. 439.  July, 2004.  arXiv:astro-ph/0310756 [astro-ph].
  * [7] D. Clowe, M. Bradač, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones, and D. Zaritsky, “A direct empirical proof of the existence of dark matter*,” The Astrophysical Journal 648 no. 2, (Aug, 2006) L109. https://dx.doi.org/10.1086/508162.
  * [8] C. L. Bennett, D. Larson, et al., “Nine-year wilkinson microwave anisotropy probe (wmap) observations: Final maps and results,” The Astrophysical Journal Supplement Series 208 no. 2, (Sep, 2013) 20. https://dx.doi.org/10.1088/0067-0049/208/2/20.
  * [9] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [10] Particle Data Group Collaboration, R. L. Workman and Others, “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [11] J. Billard et al., “Direct detection of dark matter—APPEC committee report*,” Rept. Prog. Phys. 85 no. 5, (2022) 056201, arXiv:2104.07634 [hep-ex].
  * [12] M. Pato, L. Baudis, G. Bertone, R. Ruiz de Austri, L. E. Strigari, and R. Trotta, “Complementarity of Dark Matter Direct Detection Targets,” Phys. Rev. D 83 (2011) 083505, arXiv:1012.3458 [astro-ph.CO].
  * [13] J. E. Kim, “Weak-interaction singlet and strong $\mathrm{CP}$ invariance,” Phys. Rev. Lett. 43 (Jul, 1979) 103–107. https://link.aps.org/doi/10.1103/PhysRevLett.43.103.
  * [14] M. Shifman, A. Vainshtein, and V. Zakharov, “Can confinement ensure natural cp invariance of strong interactions?” Nuclear Physics B 166 no. 3, (1980) 493–506. https://www.sciencedirect.com/science/article/pii/0550321380902096.
  * [15] M. Dine, W. Fischler, and M. Srednicki, “A simple solution to the strong cp problem with a harmless axion,” Physics Letters B 104 no. 3, (1981) 199–202. https://www.sciencedirect.com/science/article/pii/0370269381905906.
  * [16] A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31 (1980) 260.
  * [17] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [18] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * [19] F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [20] E. Witten, “Some properties of o(32) superstrings,” Physics Letters B 149 no. 4, (1984) 351–356. https://www.sciencedirect.com/science/article/pii/0370269384904222.
  * [21] P. Svrcek and E. Witten, “Axions In String Theory,” JHEP 06 (2006) 051, arXiv:hep-th/0605206.
  * [22] J. P. Conlon, “The QCD axion and moduli stabilisation,” JHEP 05 (2006) 078, arXiv:hep-th/0602233.
  * [23] K.-S. Choi, H. P. Nilles, S. Ramos-Sanchez, and P. K. S. Vaudrevange, “Accions,” Phys. Lett. B 675 (2009) 381–386, arXiv:0902.3070 [hep-th].
  * [24] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, “String Axiverse,” Phys. Rev. D 81 (2010) 123530, arXiv:0905.4720 [hep-th].
  * [25] B. S. Acharya, K. Bobkov, and P. Kumar, “An M Theory Solution to the Strong CP Problem and Constraints on the Axiverse,” JHEP 11 (2010) 105, arXiv:1004.5138 [hep-th].
  * [26] M. Cicoli, M. Goodsell, and A. Ringwald, “The type IIB string axiverse and its low-energy phenomenology,” JHEP 10 (2012) 146, arXiv:1206.0819 [hep-th].
  * [27] J. Halverson, C. Long, and P. Nath, “Ultralight axion in supersymmetry and strings and cosmology at small scales,” Phys. Rev. D 96 no. 5, (2017) 056025, arXiv:1703.07779 [hep-ph].
  * [28] M. Demirtas, C. Long, L. McAllister, and M. Stillman, “The Kreuzer-Skarke Axiverse,” JHEP 04 (2020) 138, arXiv:1808.01282 [hep-th].
  * [29] S. Chigusa, M. Hazumi, E. D. Herbschleb, N. Mizuochi, and K. Nakayama, “Light Dark Matter Search with Nitrogen-Vacancy Centers in Diamonds,” arXiv:2302.12756 [hep-ph].
  * [30] T. Fukui, Y. Doi, et al., “Perfect selective alignment of nitrogen-vacancy centers in diamond,” Applied Physics Express 7 no. 5, (2014) 055201.
  * [31] T. Miyazaki, Y. Miyamoto, T. Makino, H. Kato, S. Yamasaki, T. Fukui, Y. Doi, N. Tokuda, M. Hatano, and N. Mizuochi, “Atomistic mechanism of perfect alignment of nitrogen-vacancy centers in diamond,” Applied Physics Letters 105 no. 26, (2014) 261601.
  * [32] J. Michl, T. Teraji, et al., “Perfect alignment and preferential orientation of nitrogen-vacancy centers during chemical vapor deposition diamond growth on (111) surfaces,” Applied Physics Letters 104 no. 10, (2014) 102407.
  * [33] M. Lesik, J.-P. Tetienne, A. Tallaire, J. Achard, V. Mille, A. Gicquel, J.-F. Roch, and V. Jacques, “Perfect preferential orientation of nitrogen-vacancy defects in a synthetic diamond sample,” Applied Physics Letters 104 no. 11, (2014) 113107.
  * [34] J. F. Barry, J. M. Schloss, E. Bauch, M. J. Turner, C. A. Hart, L. M. Pham, and R. L. Walsworth, “Sensitivity optimization for nv-diamond magnetometry,” Rev. Mod. Phys. 92 (Mar, 2020) 015004. https://link.aps.org/doi/10.1103/RevModPhys.92.015004.
  * [35] “Appendix ii - table of constants and conversion factors,” in Handbook of MRI Pulse Sequences, M. A. BERNSTEIN, K. F. KING, and X. J. ZHOU, eds., pp. 960–962.  Academic Press, Burlington, 2004.  https://www.sciencedirect.com/science/article/pii/B9780120928613500261.
  * [36] G. Vidal and C. M. Dawson, “Universal quantum circuit for two-qubit transformations with three controlled-not gates,” Phys. Rev. A 69 (Jan, 2004) 010301. https://link.aps.org/doi/10.1103/PhysRevA.69.010301.
  * [37] P. Neumann, J. Beck, M. Steiner, F. Rempp, H. Fedder, P. R. Hemmer, J. Wrachtrup, and F. Jelezko, “Single-shot readout of a single nuclear spin,” Science 329 no. 5991, (2010) 542–544, https://www.science.org/doi/pdf/10.1126/science.1189075. https://www.science.org/doi/abs/10.1126/science.1189075.
  * [38] N. F. Ramsey, “A molecular beam resonance method with separated oscillating fields,” Phys. Rev. 78 (Jun, 1950) 695–699. https://link.aps.org/doi/10.1103/PhysRev.78.695.
  * [39] A. Jarmola, V. M. Acosta, K. Jensen, S. Chemerisov, and D. Budker, “Temperature- and magnetic-field-dependent longitudinal spin relaxation in nitrogen-vacancy ensembles in diamond,” Phys. Rev. Lett. 108 (May, 2012) 197601. https://link.aps.org/doi/10.1103/PhysRevLett.108.197601.
  * [40] T. Rosskopf, A. Dussaux, K. Ohashi, M. Loretz, R. Schirhagl, H. Watanabe, S. Shikata, K. M. Itoh, and C. L. Degen, “Investigation of surface magnetic noise by shallow spins in diamond,” Phys. Rev. Lett. 112 (Apr, 2014) 147602. https://link.aps.org/doi/10.1103/PhysRevLett.112.147602.
  * [41] N. Aslam, M. Pfender, et al., “Nanoscale nuclear magnetic resonance with chemical resolution,” Science 357 no. 6346, (2017) 67–71, https://www.science.org/doi/pdf/10.1126/science.aam8697. https://www.science.org/doi/abs/10.1126/science.aam8697.
  * [42] N. Mizuochi, P. Neumann, et al., “Coherence of single spins coupled to a nuclear spin bath of varying density,” Physical Review B 80 no. 49, (Jul, 2009) 041201 (R). https://doi.org/10.1103/PhysRevB.80.041201.
  * [43] G. Waldherr, J. Beck, P. Neumann, R. S. Said, M. Nitsche, M. L. Markham, D. J. Twitchen, J. Twamley, F. Jelezko, and J. Wrachtrup, “High-dynamic-range magnetometry with a single nuclear spin in diamond,” Nature Nanotechnology 7 (2012) 105–108.
  * [44] A. Jarmola, A. Berzins, J. Smits, K. Smits, J. Prikulis, F. Gahbauer, R. Ferber, D. Erts, M. Auzinsh, and D. Budker, “Longitudinal spin-relaxation in nitrogen-vacancy centers in electron irradiated diamond,” Applied Physics Letters 107 no. 24, (12, 2015) 242403.
  * [45] E. Herbschleb, H. Kato, et al., “Ultra-long coherence times amongst room-temperature solid-state spins,” Nature communications 10 no. 1, (2019) 3766.
  * [46] E. L. Hahn, “Spin echoes,” Phys. Rev. 80 (Nov, 1950) 580–594. https://link.aps.org/doi/10.1103/PhysRev.80.580.
  * [47] H. Y. Carr and E. M. Purcell, “Effects of diffusion on free precession in nuclear magnetic resonance experiments,” Phys. Rev. 94 (May, 1954) 630–638. https://link.aps.org/doi/10.1103/PhysRev.94.630.
  * [48] S. Meiboom and D. Gill, “Modified spin-echo method for measuring nuclear relaxation times,” Review of scientific instruments 29 no. 8, (1958) 688–691.
  * [49] T. Wolf, P. Neumann, K. Nakamura, H. Sumiya, T. Ohshima, J. Isoya, and J. Wrachtrup, “Subpicotesla diamond magnetometry,” Phys. Rev. X 5 (Oct, 2015) 041001. https://link.aps.org/doi/10.1103/PhysRevX.5.041001.
  * [50] J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the invisible axion,” Physics Letters B 120 no. 1, (1983) 127–132. https://www.sciencedirect.com/science/article/pii/0370269383906378.
  * [51] L. Abbott and P. Sikivie, “A cosmological bound on the invisible axion,” Physics Letters B 120 no. 1, (1983) 133–136. https://www.sciencedirect.com/science/article/pii/037026938390638X.
  * [52] M. Dine and W. Fischler, “The not-so-harmless axion,” Physics Letters B 120 no. 1, (1983) 137–141. https://www.sciencedirect.com/science/article/pii/0370269383906391.
  * [53] M. Kawasaki and K. Nakayama, “Axions: Theory and Cosmological Role,” Ann. Rev. Nucl. Part. Sci. 63 (2013) 69–95, arXiv:1301.1123 [hep-ph].
  * [54] D. J. E. Marsh, “Axion Cosmology,” Phys. Rept. 643 (2016) 1–79, arXiv:1510.07633 [astro-ph.CO].
  * [55] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, “The landscape of QCD axion models,” Phys. Rept. 870 (2020) 1–117, arXiv:2003.01100 [hep-ph].
  * [56] W. Hu, R. Barkana, and A. Gruzinov, “Cold and fuzzy dark matter,” Phys. Rev. Lett. 85 (2000) 1158–1161, arXiv:astro-ph/0003365.
  * [57] C. O’HARE, “cajohare/axionlimits: Axionlimits.” July, 2020.  https://doi.org/10.5281/zenodo.3932430.
  * [58] Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [59] Y. Ema, K. Hamaguchi, T. Moroi, and K. Nakayama, “Flaxion: a minimal extension to solve puzzles in the standard model,” JHEP 01 (2017) 096, arXiv:1612.05492 [hep-ph].
  * [60] L. Calibbi, F. Goertz, D. Redigolo, R. Ziegler, and J. Zupan, “Minimal axion model from flavor,” Phys. Rev. D 95 no. 9, (2017) 095009, arXiv:1612.08040 [hep-ph].
  * [61] Y. Ema, M. Pospelov, and A. Ray, “Probing Earth-Bound Dark Matter with Nuclear Reactors,” arXiv:2402.03431 [hep-ph].
  * [62] W. W. True, “Nitrogen-14 and the shell model,” Phys. Rev. 130 (May, 1963) 1530–1537. https://link.aps.org/doi/10.1103/PhysRev.130.1530.
  * [63] J. F. Barry, M. H. Steinecker, S. T. Alsid, J. Majumder, L. M. Pham, M. F. O’Keefe, and D. A. Braje, “Sensitive AC and DC Magnetometry with Nitrogen-Vacancy Center Ensembles in Diamond,” arXiv:2305.06269 [quant-ph].
  * [64] J. A. Dror, S. Gori, J. M. Leedom, and N. L. Rodd, “Sensitivity of Spin-Precession Axion Experiments,” Phys. Rev. Lett. 130 no. 18, (2023) 181801, arXiv:2210.06481 [hep-ph].
  * [65] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, “Asymptotic formulae for likelihood-based tests of new physics,” Eur. Phys. J. C 71 (2011) 1554, arXiv:1007.1727 [physics.data-an]. [Erratum: Eur.Phys.J.C 73, 2501 (2013)].
  * [66] E. D. Herbschleb, S. Chigusa, R. Kawase, H. Kawashima, M. Hazumi, K. Nakayama, and N. Mizuochi, “Robust quantum sensing via the standard deviation,” to appear .
  * [67] M. Buschmann, C. Dessert, J. W. Foster, A. J. Long, and B. R. Safdi, “Upper Limit on the QCD Axion Mass from Isolated Neutron Star Cooling,” Phys. Rev. Lett. 128 no. 9, (2022) 091102, arXiv:2111.09892 [hep-ph].
  * [68] G. Vasilakis, J. Brown, T. Kornack, and M. Romalis, “Limits on new long range nuclear spin-dependent forces set with a k- he 3 comagnetometer,” Physical review letters 103 no. 26, (2009) 261801.
  * [69] K. Wei et al., “Dark matter search with a strongly-coupled hybrid spin system,” arXiv:2306.08039 [hep-ph].
  * [70] Z. Xu et al., “Constraining Ultralight Dark Matter through an Accelerated Resonant Search,” arXiv:2309.16600 [hep-ph].
  * [71] A. Lella, P. Carenza, G. Co’, G. Lucente, M. Giannotti, A. Mirizzi, and T. Rauscher, “Getting the most on supernova axions,” Phys. Rev. D 109 no. 2, (2024) 023001, arXiv:2306.01048 [hep-ph].
  * [72] I. M. Bloch, Y. Hochberg, E. Kuflik, and T. Volansky, “Axion-like Relics: New Constraints from Old Comagnetometer Data,” JHEP 01 (2020) 167, arXiv:1907.03767 [hep-ph].
  * [73] C. Brandenstein, S. Stelzl, E. Gutsmiedl, W. Schott, A. Weiler, and P. Fierlinger, “Towards an electrostatic storage ring for fundamental physics measurements,” EPJ Web Conf. 282 (2023) 01017, arXiv:2211.08439 [hep-ex].
  * [74] D. F. Jackson Kimball et al., “Overview of the Cosmic Axion Spin Precession Experiment (CASPEr),” Springer Proc. Phys. 245 (2020) 105–121, arXiv:1711.08999 [physics.ins-det].
  * [75] C. Gao, W. Halperin, Y. Kahn, M. Nguyen, J. Schütte-Engel, and J. W. Scott, “Axion Wind Detection with the Homogeneous Precession Domain of Superfluid Helium-3,” Phys. Rev. Lett. 129 no. 21, (2022) 211801, arXiv:2208.14454 [hep-ph].
  * [76] P. W. Graham, S. Hacıömeroğlu, D. E. Kaplan, Z. Omarov, S. Rajendran, and Y. K. Semertzidis, “Storage ring probes of dark matter and dark energy,” Phys. Rev. D 103 no. 5, (2021) 055010, arXiv:2005.11867 [hep-ph].
  * [77] S. Chigusa, T. Moroi, K. Nakayama, and T. Sichanugrist, “Dark matter detection using nuclear magnetization in magnet with hyperfine interaction,” Phys. Rev. D 108 no. 9, (2023) 095007, arXiv:2307.08577 [hep-ph].
  * [78] M. Schreck, S. Gsell, R. Brescia, and M. Fischer, “Ion bombardment induced buried lateral growth: the key mechanism for the synthesis of single crystal diamond wafers,” Scientific reports 7 no. 1, (2017) 44462.
  * [79] J. R. Rabeau, P. Reichart, G. Tamanyan, D. N. Jamieson, S. Prawer, F. Jelezko, T. Gaebel, I. Popa, M. Domhan, and J. Wrachtrup, “Implantation of labelled single nitrogen vacancy centers in diamond using N15,” Applied Physics Letters 88 no. 2, (01, 2006) 023113.
  * [80] A. Dréau, P. Spinicelli, J. R. Maze, J.-F. Roch, and V. Jacques, “Single-shot readout of multiple nuclear spin qubits in diamond under ambient conditions,” Phys. Rev. Lett. 110 (Feb, 2013) 060502.

##  Appendix A Spin synthesis

In this appendix, we focus on the synthesis of the spin $S=1/2$ and the
orbital angular momentum $L=\ell$ and derive how spin operators act on the
eigenstates of the total spin $J$. Let $\ket{\uparrow}$ and $\ket{\downarrow}$
be the spin-up and spin-down states, respectively, and $\ket{m}$
($m=-\ell,\dots,\ell$) be the eigenstates of the orbital angular momentum with
$L_{z}=m$. The synthesized states of these two quantum numbers decompose into
two groups with total angular momenta $J=\ell\pm\frac{1}{2}$. By
parameterizing these states as $\ket{J,M}$ with $M=-J,\dots,J$ representing
the $z$-component of the total angular momentum, the matrix elements of the
total spin operators are characterized by the Clebsch-Gordan coefficients

$\displaystyle\Braket{J,M\pm 1}{J^{\pm}}{J,M}=\sqrt{(J\mp M)(J\pm M+1)},$ (48)

with $J^{\pm}\equiv J^{x}\pm iJ^{y}$. Using this expression, we can relate the
eigenstates of various spins as follows:

$\displaystyle\Ket{J=\ell+\frac{1}{2},M}$
$\displaystyle=\frac{1}{\sqrt{2J}}\left(\sqrt{J+M}\ket{\uparrow}\ket{M-1/2}+\sqrt{J-M}\ket{\downarrow}\ket{M+1/2}\right),$
(49) $\displaystyle\Ket{J=\ell-\frac{1}{2},M}$
$\displaystyle=\frac{1}{\sqrt{2J+2}}\left(\sqrt{J-M+1}\ket{\uparrow}\ket{M-1/2}-\sqrt{J+M+1}\ket{\downarrow}\ket{M+1/2}\right),$
(50)

or equivalently,

$\displaystyle\ket{\uparrow}\ket{m}$
$\displaystyle=\frac{1}{\sqrt{2\ell+1}}\left(\sqrt{\ell+m+1}\Ket{\ell+\frac{1}{2},m+\frac{1}{2}}+\sqrt{\ell-m}\Ket{\ell-\frac{1}{2},m+\frac{1}{2}}\right),$
(51) $\displaystyle\ket{\downarrow}\ket{m}$
$\displaystyle=\frac{1}{\sqrt{2\ell+1}}\left(\sqrt{\ell-m+1}\Ket{\ell+\frac{1}{2},m-\frac{1}{2}}-\sqrt{\ell+m}\Ket{\ell-\frac{1}{2},m-\frac{1}{2}}\right).$
(52)

Using Eqs. 51 and 52, we can calculate all the non-zero matrix elements of the
spin operators as follows:

$\displaystyle\Braket{J=\ell+\frac{1}{2},M\pm 1}{S^{\pm}}{J,M}$
$\displaystyle=\frac{\sqrt{(J\mp M)(J\pm M+1)}}{2J},$ (53)
$\displaystyle\Braket{J=\ell+\frac{1}{2},M}{S^{z}}{J,M}$
$\displaystyle=\frac{M}{2J},$ (54)
$\displaystyle\Braket{J=\ell-\frac{1}{2},M\pm 1}{S^{\pm}}{J,M}$
$\displaystyle=-\frac{\sqrt{(J\mp M)(J\pm M+1)}}{2J+2},$ (55)
$\displaystyle\Braket{J=\ell-\frac{1}{2},M}{S^{z}}{J,M}$
$\displaystyle=-\frac{M}{2J+2},$ (56)

with $S^{\pm}\equiv S^{x}\pm iS^{y}$. From the above equations, we see that,
for a fixed value of $J$, the spin operators effectively act as $SU(2)$
generators in the spin-$J$ representation with a non-trivial factor,

$\displaystyle\Braket{J,M^{\prime}}{\vec{S}}{J,M}=\pm\frac{1}{2\ell+1}\Braket{J,M^{\prime}}{\vec{J}}{J,M},$
(57)

for $J=\ell\pm\frac{1}{2}$.

Now the calculation so far can be applied to the case of the ${\ce{{}^{1}4N}}$
spin, which has $I=1$ composed of a neutron and a proton in individual
$(1p)_{1/2}$ orbitals. First, each neutron and proton resides in the $1p$
orbital with $\ell=1$, resulting in the total angular momentum
$J_{\chi}=\ell-1/2=1/2$ (${\chi}=n,p$). According to Eq. 57, we obtain the
effective relationships among operators

$\displaystyle\vec{S}_{\chi}\sim-\frac{1}{3}\vec{J}_{\chi}.~{}~{}(f=n,p)$ (58)

Since either one of the total angular momentum operators, say $\vec{J}_{n}$,
has spin $J_{n}=1/2$, we can repeat the same estimation, combining it with
$\vec{J}_{p}$ to obtain $I=1$ states. Again according to Eq. 57, the angular
momentum operators of nucleons are related to the ${\ce{{}^{1}4N}}$ spin
operator $\vec{I}$ as

$\displaystyle\vec{J}_{n}\sim\vec{J}_{p}\sim\frac{1}{2}\vec{I}.$ (59)

Therefore, the coefficient $-1/6$ in Eq. 27 is successfully reconstructed.

##  Appendix B Calculation of the power spectral density and the quantum noise

In this appendix, we provide a detailed calculation of the PSD defined in Eq.
37 starting from Eq. 35. Firstly, the density matrix of the nuclear spin state
in an ensemble of $N$ NV centers before the $j$-th measurement, $\rho_{j}$, is
expressed as

$\displaystyle\rho_{j}$
$\displaystyle\simeq\bigotimes_{\ell=1}^{N}\rho_{j\ell},$ (60)
$\displaystyle\rho_{j\ell}$
$\displaystyle\equiv\frac{1}{2}\left\\{(1-2{F}_{j})\outerproduct{+}{+}+\sqrt{1-4{F}_{j}^{2}}\outerproduct{+}{0}+\sqrt{1-4{F}_{j}^{2}}\outerproduct{0}{+}+(1+2{F}_{j})\outerproduct{0}{0}\right\\},$
(61)

for fixed values of ${\phi}$ and $\vec{v}_{a}$, where ${F}_{j}$ is given by
Eqs. 14, 16 and 17 for the Ramsey, the Hahn echo, and the dynamic decoupling
sequences, respectively, with the replacement ${\phi}\to m_{a}t_{j}+{\phi}$.
For notational simplicity, we omit the indices $j$ and $\ell$ for each bra and
ket, but they are assumed implicitly. Thus, $\rho_{j}$ is a
$2^{N}$-dimensional density matrix.

Our next task is to evaluate the trace factor
${\mathrm{Tr}\,\left[{\rho_{jj^{\prime}}M_{jj^{\prime}}^{z}}\right]}$ for
various choices of $j$ and $j^{\prime}$, where
$\rho_{jj^{\prime}}\equiv\rho_{j}\otimes\rho_{j^{\prime}}$ and
$M_{jj^{\prime}}^{z}\equiv M_{j}^{z}\otimes M_{j^{\prime}}^{z}$ for $j\neq
j^{\prime}$ and $\rho_{jj}\equiv\rho_{j}$ and $M_{jj}^{z}\equiv
M_{j}^{z}M_{j}^{z}$. Firstly, when $j=j^{\prime}$, we obtain

$\displaystyle{\mathrm{Tr}\,\left[{\rho_{j}M_{j}^{z}M_{j}^{z}}\right]}$
$\displaystyle=\frac{1}{4N^{2}}\sum_{\ell}1+\frac{1}{4N^{2}}\sum_{\ell\neq\ell^{\prime}}{\mathrm{Tr}\,\left[{\rho_{j\ell}\otimes\rho_{j\ell^{\prime}}\sigma_{j\ell}^{z}\sigma_{j\ell^{\prime}}^{z}}\right]}$
$\displaystyle=\frac{1}{4N}+\frac{N(N-1)}{N^{2}}F_{j}^{2}.$ (62)

For a different time two-point function with $j\neq j^{\prime}$, we can
instead decompose the trace in two parts for the time $t_{j}$ and
$t_{j^{\prime}}$ and obtain

$\displaystyle{\mathrm{Tr}\,\left[{(\rho_{j}\otimes\rho_{j^{\prime}})(M_{j}^{z}\otimes
M_{j^{\prime}}^{z})}\right]}$
$\displaystyle={\mathrm{Tr}\,\left[{\rho_{j}M_{j}^{z}}\right]}\,{\mathrm{Tr}\,\left[{\rho_{j^{\prime}}M_{j^{\prime}}^{z}}\right]}$
$\displaystyle={F}_{j}{F}_{j^{\prime}}.$ (63)

We can combine these expressions in a compact form

$\displaystyle{\mathrm{Tr}\,\left[{\rho_{jj^{\prime}}M_{j}^{z}M_{j^{\prime}}^{z}}\right]}\simeq\frac{1}{4N}\delta_{jj^{\prime}}+{F}_{j}{F}_{j^{\prime}},$
(64)

where subleading terms of $N\gg 1$ are neglected. Recalling that the ensemble
average of a single data is zero, $\Braket{M_{j}}=0$, the two-point function
defined in the main text,

$\displaystyle C_{jj^{\prime}}=$ $\displaystyle\left.\frac{1}{2\pi}\int
d{\phi}\,\frac{1}{4\pi}\int d\hat{v}_{a}\,\frac{1}{2\pi}\int
d{\phi}^{\prime}\,\frac{1}{4\pi}\int
d\hat{v}_{a}^{\prime}\,{\mathrm{Tr}\,\left[{\rho_{jj^{\prime}}M_{jj^{\prime}}^{z}}\right]}\right|_{v_{a}=10^{-3}}$
$\displaystyle\times\left[\Theta(|t_{j}-t_{j^{\prime}}|-\tau_{a})+8\pi^{2}\delta({\phi}-{\phi}^{\prime})\delta(\hat{v}_{a}-\hat{v}_{a}^{\prime})\Theta(\tau_{a}-|t_{j}-t_{j^{\prime}}|)\right],$
(35)

is thus calculated as

$\displaystyle C_{jj^{\prime}}$
$\displaystyle\simeq\frac{1}{4N}\delta_{jj^{\prime}}+\frac{1}{2\pi}\int
d{\phi}\,\frac{1}{4\pi}\int
d\hat{v}_{a}\,{F}_{j}{F}_{j^{\prime}}\Big{|}_{{\phi}^{\prime}={\phi},\hat{v}_{a}^{\prime}=\hat{v}_{a},v_{a}=10^{-3}}\Theta(\tau_{a}-|t_{j}-t_{j^{\prime}}|)$
$\displaystyle=\frac{1}{4N}\delta_{jj^{\prime}}+\mathcal{A}\cos\left[m_{a}(t_{j}-t_{j^{\prime}})\right]\Theta(\tau_{a}-|t_{j}-t_{j^{\prime}}|),$
(65)

with the protocol-dependent coefficient $\mathcal{A}$ defined as

$\displaystyle\mathcal{A}\equiv\begin{cases}\dfrac{\rho_{a}v_{0}^{2}}{27\tilde{f}_{a}^{2}m_{a}^{2}}\sin^{2}\dfrac{m_{a}\tau}{2},&(\text{Ramsey})\\\\[10.0pt]
\dfrac{4\rho_{a}v_{0}^{2}}{27\tilde{f}_{a}^{2}m_{a}^{2}}\sin^{4}\dfrac{m_{a}\tau}{4},&(\text{Hahn
echo})\\\\[10.0pt]
\dfrac{\rho_{a}v_{0}^{2}}{27\tilde{f}_{a}^{2}m_{a}^{2}}\sin^{2}\dfrac{m_{a}\tau}{2}\tan^{2}\dfrac{m_{a}\tau}{2(N_{\pi}+1)},&(\text{DD})\\\\[10.0pt]
\end{cases}$ (66)

where $v_{0}=10^{-3}$ denotes the typical axion velocity.

Next, we calculate the PSD ${\mathcal{P}_{k}}\equiv\Braket{O_{k}}$ using the
operator $\mathcal{O}_{k}$ defined as follows:

$\displaystyle\mathcal{O}_{k}\equiv\frac{\tau^{2}}{{t_{\mathrm{obs}}}}\sum_{j,j^{\prime}}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}M_{j}^{z}M_{j^{\prime}}^{z}.$ (36)

Using a modified expression

$\displaystyle{\mathcal{P}_{k}}=\frac{\tau^{2}}{{t_{\mathrm{obs}}}}\sum_{j,j^{\prime}}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}C_{jj^{\prime}},$ (37)

an easy way to accomplish this task is to consider the continuum limit as
follows:

$\displaystyle{\mathcal{P}_{k}}$
$\displaystyle\simeq\frac{1}{{t_{\mathrm{obs}}}}\int_{0}^{{t_{\mathrm{obs}}}}dt\,\int_{0}^{{t_{\mathrm{obs}}}}dt^{\prime}\,e^{i\omega_{k}(t-t^{\prime})}C(t,t^{\prime}),$
(67)

where $\omega_{k}\equiv 2\pi k/{t_{\mathrm{obs}}}$ and the function
$C(t,t^{\prime})$ is defined as a natural extension of $C_{jj^{\prime}}$ to
the continuous choice of time. By substituting Eq. 65 into the above
expression, we obtain

$\displaystyle{\mathcal{P}_{k}}\simeq\frac{\tau}{4N}+\begin{cases}\dfrac{2\mathcal{A}}{{t_{\mathrm{obs}}}\Delta\omega_{k}^{2}}\sin^{2}\dfrac{{t_{\mathrm{obs}}}\Delta\omega_{k}}{2},&({t_{\mathrm{obs}}}<\tau_{a})\\\\[10.0pt]
\dfrac{2\mathcal{A}}{{t_{\mathrm{obs}}}\Delta\omega_{k}^{2}}\sin^{2}\dfrac{\tau_{a}\Delta\omega_{k}}{2}+\dfrac{{t_{\mathrm{obs}}}-\tau_{a}}{{t_{\mathrm{obs}}}\Delta\omega_{k}}\mathcal{A}\sin\left[\tau_{a}\Delta\omega_{k}\right],&({t_{\mathrm{obs}}}>\tau_{a})\end{cases}$
(68)

with $\Delta\omega_{k}\equiv\omega_{k}-m_{a}$, where we have neglected the
fast oscillation terms.

Finally, we evaluate the quantum noise on the PSD without the axion effect,

$\displaystyle{\mathcal{B}_{k}}\equiv\left.\sqrt{\Braket{\mathcal{O}_{k}^{2}}-\Braket{\mathcal{O}_{k}}^{2}}\right|_{B_{N}=0}.$
(40)

The first term in the square root can be deformed as

$\displaystyle\Braket{\mathcal{O}_{k}^{2}}$
$\displaystyle=\Braket{\left(\frac{\tau^{2}}{{t_{\mathrm{obs}}}}\sum_{j,j^{\prime}}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}M_{j}^{z}M_{j^{\prime}}^{z}\right)^{2}}$
$\displaystyle=\Braket{\left(\frac{\tau^{2}}{{t_{\mathrm{obs}}}}\sum_{j,j^{\prime}}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}\frac{1}{4N^{2}}\sum_{\ell,\ell^{\prime}}\sigma_{j\ell}^{z}\sigma_{j^{\prime}\ell^{\prime}}^{z}\right)^{2}}$
$\displaystyle=\frac{\tau^{4}}{16{t_{\mathrm{obs}}^{2}}N^{4}}\Braket{\left(\sum_{j,\ell}\bm{1}+\sum_{(j,\ell)\neq(j^{\prime},\ell^{\prime})}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}\sigma_{j\ell}^{z}\sigma_{j^{\prime}\ell^{\prime}}^{z}\right)^{2}},$
(69)

where $\bm{1}$ is the identity operator. To go further, we note that the odd
number of Pauli matrices for a certain combination of $(j,\ell)$ leads to the
vanishing contribution
${\mathrm{Tr}\,\left[{\rho_{j\ell}\sigma_{j\ell}^{z}}\right]}\big{|}_{v_{a}=0}=0$.
Thus, the only remaining contribution comes from the terms proportional to the
identity matrix. In the parenthesis of the third line of the previous
equation, the first term trivially leads to such a contribution with size
${N_{\mathrm{obs}}^{2}}N^{2}$, while the second term also contributes as

$\displaystyle\left(\sum_{(j,\ell)\neq(j^{\prime},\ell^{\prime})}e^{2\pi
ik(j-j^{\prime})/{N_{\mathrm{obs}}}}\sigma_{j\ell}^{z}\sigma_{j^{\prime}\ell^{\prime}}^{z}\right)^{2}$
$\displaystyle=\sum_{(j_{1},\ell_{1})\neq(j_{2},\ell_{2})}\sum_{(j_{3},\ell_{3})\neq(j_{4},\ell_{4})}e^{2\pi
ik(j_{1}-j_{2}+j_{3}-j_{4})/{N_{\mathrm{obs}}}}\sigma_{j_{1}\ell_{1}}^{z}\sigma_{j_{2}\ell_{2}}^{z}\sigma_{j_{3}\ell_{3}}^{z}\sigma_{j_{4}\ell_{4}}^{z}$
$\displaystyle=\sum_{(j_{1},\ell_{1})\neq(j_{2},\ell_{2})}\left(1+e^{4\pi
ik(j_{1}-j_{2})/{N_{\mathrm{obs}}}}\right)+\cdots$
$\displaystyle={N_{\mathrm{obs}}}N({N_{\mathrm{obs}}}N-1)+\sum_{j_{1},j_{2}}\sum_{\ell_{1},\ell_{2}}e^{4\pi
ik(j_{1}-j_{2})/{N_{\mathrm{obs}}}}-\sum_{j,\ell}1+\cdots$
$\displaystyle={N_{\mathrm{obs}}}N({N_{\mathrm{obs}}}N-2)+{N_{\mathrm{obs}}^{2}}N^{2}\delta_{k,0}+\cdots,$
(70)

where the identity operator $\bm{1}$ is implicit, while dots represent terms
with remnant Pauli matrices. Substituting this result in the original
definition, we obtain

$\displaystyle{\mathcal{B}_{0}}$ $\displaystyle\simeq\frac{\tau}{2\sqrt{2}N},$
(71) $\displaystyle{\mathcal{B}_{k\neq 0}}$
$\displaystyle\simeq\frac{\tau}{4N},$ (72)

where we neglect the subleading terms of ${N_{\mathrm{obs}}}$ and $N$.