# On the Sensitivity of Spin-Precession Axion Experiments

Jeff A. Dror  Department of Physics, University of California Santa Cruz, 1156
High St., Santa Cruz, CA 95064, USA
and Santa Cruz Institute for Particle Physics, 1156 High St., Santa Cruz, CA
95064, USA  Stefania Gori  Department of Physics, University of California
Santa Cruz, 1156 High St., Santa Cruz, CA 95064, USA
and Santa Cruz Institute for Particle Physics, 1156 High St., Santa Cruz, CA
95064, USA  Jacob M. Leedom  Deutsches Elektronen-Synchrotron DESY, Notkestr.
85, 22607 Hamburg, Germany  Nicholas L. Rodd  Theoretical Physics Department,
CERN, 1 Esplanade des Particules, CH-1211 Geneva 23, Switzerland

###### Abstract

We study the signal and background that arise in nuclear magnetic resonance
searches for axion dark matter, finding key differences with the existing
literature. We find that spin-precession instruments are much more sensitive
than what has been previously estimated in a sizable range of axion masses,
with sensitivity improvement of up to a factor of 100 using a ${}^{129}{\rm
Xe}$ sample. This improves the detection prospects for the QCD axion, and we
estimate the experimental requirements to reach this motivated target. Our
results apply to both the axion electric and magnetic dipole moment operators.

††preprint: CERN-TH-2022-163††preprint: DESY 22-155

The axion, a light scalar whose leading interactions exhibit a shift symmetry,
is one of the most compelling extensions to the Standard Model of particle
physics. Originally proposed as an elegant solution to the Strong CP problem
Peccei:1977hh ; Peccei:1977ur ; Weinberg:1977ma ; Wilczek:1977pj , axions have
since been appreciated for both their ubiquity in string theory Svrcek:2006yi
; Arvanitaki:2009fg ; Halverson:2019cmy  and the generic expectation that they
contribute to cold dark matter through, for example, misalignment production
Abbott:1982af ; Preskill:1982cy ; Dine:1982ah . The continuous shift symmetry
of the axion leads to a natural expectation that the axion should be an
extraordinarily light state, with mass $m_{a}\ll 1~{}{\rm eV}$ suppressed by
its decay constant $f_{a}$, which arises from breaking the continuous symmetry
to a discrete one via instanton effects. Such a small mass implies that axion
dark matter on Earth should be well described by a classical wave.

At low energies, the possible interactions of an axion, $a$, with a nucleon
$N$ take the following form,

${\cal L}\supset
g_{N}(\partial_{\mu}a)\bar{N}\gamma^{\mu}\gamma_{5}N-\frac{i}{2}g_{d}a\bar{N}\sigma_{\mu\nu}\gamma_{5}NF^{\mu\nu}.$
(1)

In the presence of an axion-wave background, these interactions source an
oscillating magnetic dipole (MD), weighted by $g_{N}\propto 1/f_{a}$, and an
oscillating electric dipole (ED), weighted by $g_{d}\propto 1/m_{N}f_{a}$,
with $m_{N}$ the nucleon mass. For the QCD axion, $g_{N}$ depends on the
charge assignments and field content of the UV theory, whereas $g_{d}$ depends
on the term that resolves the strong CP problem,
$(a/f_{a})G^{\mu\nu}\tilde{G}_{\mu\nu}$, and is fixed to $g_{d}=(3.7\pm
1.5)\times 10^{-18}~{}{\rm GeV}^{-2}(10^{15}~{}{\rm GeV}/f_{a})$
Pospelov:1999mv . For a non-relativistic nucleus, these interactions lead to
the Hamiltonian, $H_{\rm int}=-2(g_{N}\nabla a+g_{d}a{\bf E}^{\ast})\cdot{\bf
S}$. Here, ${\bf S}$ is the nucleon spin operator, and ${\bf E}^{\ast}$ is the
effective electric field felt by a nucleus (and differs from the applied field
by at least two orders of magnitude due to shielding by the atomic electrons).
By drawing an analogy between the interaction of spins and a magnetic field,
the axion-nucleus interaction can be characterized as an effective axion
magnetic field,

${\bf B}_{a}(t)=-\frac{2}{\gamma}\left\\{\begin{array}[]{cr}g_{N}\nabla
a(t)&\text{(MD)},\\\\[2.84544pt] g_{d}{\bf
E}^{\ast}a(t)&\text{(ED)},\end{array}\right.$ (2)

with $\gamma$ being the gyromagnetic ratio of the nucleon. This axion-induced
magnetic field can be detected using spin-precession techniques such as
nuclear magnetic resonance (NMR), as originally proposed in the seminal CASPEr
papers Graham:2013gfa ; Budker:2013hfa , and it was further developed
experimentally and theoretically in Refs. Stadnik:2013raa ; Abel:2017rtm ;
Wang:2017ixp ; Garcon:2019inh ; Smorra:2019qfx ; Roussy:2020ily ;
Jiang:2021dby ; Aybas:2021nvn ; JacksonKimball:2017elr ; Aybas:2021cdk . The
CASPEr approach is opening up a frontier for axion dark-matter direct
detection beyond the widely exploited axion-photon coupling (for related
proposals, see Refs. Graham:2011qk ; Sikivie:2014lha ; Graham:2017ivz ;
Chang:2017ruk ; Chang:2018dvm ; Wu:2019exd ; Terrano:2019clh ; Bloch:2019lcy ;
Yang:2019xdz ; Abusaif:2019gry ; Flambaum:2019rvu ; Stephenson:2020jzx ;
Gramolin:2020ict ; Graham:2020kai ; Gaul:2020bdq ; Arvanitaki:2021wjk ;
Kim:2021pld ; Gao:2022nuq ; Berlin:2022mia ; Lisanti:2021vij ; Lee:2022vvb ).

In this Letter, we revisit the sensitivity of spin-precession experiments,
clarifying fundamental aspects of the behavior of the expected axion signal
and noise sources. The system depends on three fundamental timescales: the
axion coherence time, $\tau_{a}\sim(4~{}{\rm neV}/m_{a})~{}{\rm sec}$, the
transverse spin-relaxation time, $T_{2}$, and the experimental integration
time, $T$. We demonstrate that there are two previously overlooked effects
that enhance the growth of the signal when $\tau_{a}<T_{2}$, a realization
that improves the prospects for QCD axion detection with spin-precession
instruments, and we comment on the requirements to achieve this goal. We
further reconsider a dominant noise source – spin-projection noise – and
demonstrate that it is larger at high axion masses, thereby reducing the
utility of using materials with large magnetic moments.

The NMR Axion Signal. We begin with a brief review of spin-precession axion
experiments. Fundamentally, they involve a macroscopic sample of atoms with
non-zero nuclear spin placed within a static magnetic field, ${\bf B}_{0}$.
The field induces a bulk magnetization ${\bf M}_{0}$ that is parallel to ${\bf
B}_{0}$. A perpendicular magnetic field – such as that induced by the axion –
will rotate the nuclear spins by a small amount. However, as soon as they do,
the spins will precess around ${\bf B}_{0}$ at the Larmor frequency,
$\omega_{0}\equiv\gamma B_{0}$, generating an oscillating transverse
magnetization, which can then be detected with a sensitive magnetometer. For
$m_{a}\sim\omega_{0}$, the effect is resonantly enhanced, with the width of
the system’s response controlled by the transverse relaxation time, $T_{2}$,
which is a macroscopic property of the sample and quantifies how long the
precession of the transverse spins can be maintained coherently. By varying
${\bf B}_{0}$, the instrument can scan a range of axion masses.

To study the dynamics of the magnetization, we turn to the Bloch equations.
Preparing the sample with ${\bf M}_{0}\propto{\bf B}_{0}\propto\hat{\bf z}$,
the Bloch equations read

$\frac{d{\bf M}}{dt}={\bf M}\times\gamma{\bf B}-\frac{M_{x}\hat{\bf
x}+M_{y}\hat{\bf y}}{T_{2}}-\frac{(M_{z}-M_{0})\hat{\bf z}}{T_{1}}.$ (3)

The longitudinal relaxation time commonly satisfies $T_{1}\gg T_{2}$, and we
will work in the limit where $T_{1}$ is much longer than any other timescale
of interest, so that it will play no further role in our discussion. We
decompose the magnetic field as ${\bf B}=B_{0}\hat{\bf z}+{\bf B}_{a}(t)$. As
$|{\bf B}_{a}(t)|\ll B_{0}$, we study the magnetization perturbatively, taking
${\bf M}\simeq{\bf M}_{0}+{\bf M}_{a}$, and we will work only to linear order
in the axion-induced fields throughout 111This is an excellent approximation.
If we take $T_{2},\tau_{a}\to\infty$, then $M_{x}(t)\sim
g_{N}v\sqrt{2\rho_{\scriptscriptstyle\textrm{DM}}}M_{0}\sin(\omega_{0}t)t$ for
the MD interaction, see Eq. (8). This suggests that with enough time an
infinite magnetization can be generated, which is inconsistent with the finite
number of spins in the system. Terms higher order in $g_{N}$ prevent this from
occurring. We can estimate when they must enter by determining when
$M_{x}(t)\sim M_{0}$. Setting $g_{N}$ to the largest allowed value (set by the
SN 1987A bound), we find $t\sim 30~{}{\rm years}$, dramatically larger than
$T_{1}$ for any sample we consider, which represents the longest time we can
interrogate the sample for. Including a finite $T_{2}$ or $\tau_{a}$ only
strengthens this conclusion. . Working to this order, $M_{z}(t)=M_{0}$,
leaving the dynamics to the transverse magnetizations,

$\displaystyle\dot{M}_{x}$
$\displaystyle=\omega_{0}M_{y}-T_{2}^{-1}M_{x}-M_{0}\gamma B_{ay},$ (4)
$\displaystyle\dot{M}_{y}$
$\displaystyle=-\omega_{0}M_{x}-T_{2}^{-1}M_{y}+M_{0}\gamma B_{ax}.$

These equations can be decoupled. If we measure the $\hat{\bf x}$ component of
the magnetization, the equation to solve is

$\displaystyle\ddot{M}_{x}+2T_{2}^{-1}\dot{M}_{x}+\omega_{0}^{2}M_{x}\simeq
F(t),$ (5) where, $\displaystyle F(t)=\gamma
M_{0}[\omega_{0}B_{ax}-\dot{B}_{ay}].$

Here and throughout, we neglect terms of ${\cal O}(1/\omega_{0}T_{2})$, as
they are significantly suppressed. Eq. (5) has reduced the system to the form
of a simple harmonic oscillator with a resonant frequency $\omega_{0}$ and
bandwidth $1/T_{2}$, that is being driven by the axion wave. Assuming
$M_{x}(0)=M_{y}(0)=0$, the solution is given by 222The solution in Eq. (6) is
for the Bloch equation before terms of $\mathcal{O}(1/\omega_{0}T_{2})$ are
dropped; cf. Eq. (5).

$M_{x}(t)=\frac{1}{\omega_{0}}\int_{0}^{t}dt^{\prime}e^{(t^{\prime}-t)/T_{2}}\sin\left[\omega_{0}(t-t^{\prime})\right]F(t^{\prime}).$
(6)

To complete our solution for the magnetization, we require a model for the
axion field. Here we treat the axion as a field with a fluctuating phase—we
show that our results are reproduced when the axion is modeled as a sum over
plane waves in the Supplementary Material (SM). The axion is taken to have
constant amplitude $a_{0}=\sqrt{2\rho_{\scriptscriptstyle\textrm{DM}}}/m_{a}$,
fixed by the local dark-matter density, and oscillates with frequency
$\omega_{a}\simeq m_{a}(1+v^{2}/2)$. The statistics of the field are then
encoded by requiring the field, which carries velocity $v\sim 10^{-3}$, obtain
a new random phase uniformly sampled on $[0,2\pi)$, every coherence time,
$\tau_{a}=2\pi/m_{a}v^{2}$ 333The coherence time is a measure of how long the
axion field can be approximated as a perfectly coherent driving force, or in
the frequency domain how long of a measurement is required to resolve the
intrinsic width of the axion. It is an approximate notion, as the transition
from perfectly coherent to incoherent happens gradually near $t\sim
2\pi/m_{a}v^{2}$, not exactly at $\tau_{a}$. For further discussion of the
coherence time, see for instance Ref. Foster:2020fln .. Each time the phase is
updated, we further update the direction of the axion field’s momentum, ${\bf
k}$ (where $|{\bf k}|=m_{a}v$), though, parametrically, our results are
insensitive to the stochastic nature of the momentum vector. In summary,

$\displaystyle F(t)$ $\displaystyle=A\cos[\omega_{a}t+\varphi(t)],$ (7)
$\displaystyle A$ $\displaystyle\equiv
2M_{0}a_{0}\left\\{\begin{array}[]{l}g_{N}\sqrt{\omega_{0}^{2}({\bf
k}\cdot\hat{\bf x})^{2}+\omega_{a}^{2}({\bf k}\cdot\hat{\bf
y})^{2}},\\\\[2.84544pt] g_{d}\omega_{0}E^{\ast},\end{array}\right.$

where we focus on the dominant term in the driving force, which carries a
phase $\varphi(t)$, uniformly sampled on $\left[0,2\pi\right)$, but shifted
from the axion phase. Furthermore, we assumed ${\bf E}^{\ast}\propto\hat{\bf
x}$ for the ED operator. From Eq. (5), a resonant response is induced in the
system when $\left|\omega_{0}-\omega_{a}\right|\lesssim\pi\,{\rm
max}[\tau_{a}^{-1},T_{2}^{-1}]$. We will assume a scan strategy such that this
condition is always satisfied, and thereby assume
$\omega_{a}\simeq\omega_{0}$.

Whenever $t\ll\tau_{a}$, the axion behaves as if it were a perfectly coherent
driving force: $\varphi$ and ${\bf k}$ are constant, so that Eq. (6) yields
the following oscillating solution

$M_{x}(t)\simeq\left(1-e^{-t/T_{2}}\right)\frac{AT_{2}}{2\omega_{0}}\sin[\omega_{0}t+\varphi].$
(8)

For $t\ll T_{2}$, the amplitude grows linearly in time. Beyond $T_{2}$,
however, the growth saturates.

We next extend these results to finite $\tau_{a}$. Due to the stochastic
variation of $\varphi$ and ${\bf k}$ for $t>\tau_{a}$, it is useful to compute
the autocorrelation function of the induced magnetization,
$C(t,t^{\prime})\equiv\langle M_{x}(t)M_{x}(t^{\prime})\rangle$,

$\displaystyle C(t,t^{\prime})$
$\displaystyle=\frac{1}{\omega_{0}^{2}}\int_{0}^{t}d\bar{t}\int_{0}^{t^{\prime}}d\bar{t}^{\prime}e^{-(t-\bar{t})/T_{2}}e^{-(t^{\prime}-\bar{t}^{\prime})/T_{2}}$
(9) $\displaystyle\times$
$\displaystyle\sin\left[\omega_{0}(t-\bar{t})\right]\sin\left[\omega_{0}(t^{\prime}-\bar{t}^{\prime})\right]\langle
F(\bar{t})F(\bar{t}^{\prime})\rangle.$

The expectation value vanishes unless the random phases are identical, which
requires $|\bar{t}-\bar{t}^{\prime}|<\tau_{a}$, so that

$\langle F(\bar{t})F(\bar{t}^{\prime})\rangle=\frac{1}{2}\langle
A^{2}\rangle\cos\left[\omega_{0}(\bar{t}-\bar{t}^{\prime})\right]\Theta(\tau_{a}-|\bar{t}-\bar{t}^{\prime}|),$
(10)

in terms of the step function $\Theta$. (Here, we have also used that
$\varphi$ and ${\bf k}$ are uncorrelated.) The remaining expectation value can
be determined by averaging over the incident direction of the axion, yielding

$\displaystyle\langle
A^{2}\rangle\simeq\left(2M_{0}a_{0}\omega_{0}\right)^{2}\left\\{\begin{array}[]{l}(g_{N}\omega_{0}v)^{2}/3,\\\\[2.84544pt]
(g_{d}E^{\ast})^{2}.\end{array}\right.$ (11)

At short times ($t,t^{\prime}\ll\tau_{a}$) where the axion is continuous, we
have (cf. Eq. (8)),

$\displaystyle C(t,t^{\prime})$ $\displaystyle=\frac{\langle A^{2}\rangle
T_{2}^{2}}{8\omega_{0}^{2}}\cos\left[\omega_{0}(t-t^{\prime})\right]$ (12)
$\displaystyle\times$
$\displaystyle\left(1-e^{-t/T_{2}}\right)\left(1-e^{-t^{\prime}/T_{2}}\right)\\!.$

At longer times, the stochastic fluctuations in the axion associated with
$\tau_{a}$ must be accounted for. Combining Eqs. (9) and (10), for
$t,t^{\prime}\gg\tau_{a}$ we can evaluate the integral by rotating coordinates
to $\bar{t}\pm\bar{t}^{\prime}$, yielding

$\displaystyle C(t,t^{\prime})$ $\displaystyle=\frac{\langle A^{2}\rangle
T_{2}\tau_{a}}{16\omega_{0}^{2}}\cos\left[\omega_{0}(t-t^{\prime})\right]$
(13) $\displaystyle\times$ $\displaystyle
e^{-(t+t^{\prime})/T_{2}}\left(e^{2{\rm
min}(t,t^{\prime})/T_{2}}-1\right)\\!.$

Eqs. (12) and (13) allow us to infer the growth of the magnetization in the
presence of a finite $\tau_{a}$, as $C(t,t)=\langle M_{x}^{2}(t)\rangle$.
Firstly, for $t\gg T_{2}$, we see that the growth saturates even with a finite
coherence time in the driving force, implying saturation occurs in this limit
regardless of the size of $\tau_{a}$. For $t\ll T_{2}$, however, the behaviors
differ,

$\lim_{t\ll T_{2}}\langle M_{x}^{2}(t)\rangle=\frac{\langle
A^{2}\rangle}{8\omega_{0}^{2}}\left\\{\begin{array}[]{lc}t^{2}&t\ll\tau_{a},\\\
\tau_{a}t&t\gg\tau_{a}.\end{array}\right.$ (14)

The first result, that the amplitude of the magnetization grows linearly with
time for $t\ll\tau_{a},\,T_{2}$, is consistent with Eq. (8). However, for
$\tau_{a}<t<T_{2}$, we see that the amplitude continues to grow, leading to an
ever-increasing size, albeit as $\sqrt{t}$. Intuitively, this transition in
behavior can be understood as the magnetization executing a random walk. For
$t>\tau_{a}$, $M_{x}(t)$ in Eq. (6) is now a sum of contributions from the
axion field at different coherence times, all of which are out of phase. The
sum is analogous to a 2D random walk with steps of length $\tau_{a}$ (given
the growth for $t<\tau_{a}$), and with a number of steps $t/\tau_{a}$, so that
we expect $|M_{x}(t)|\propto\sqrt{t\tau_{a}}$, exactly as found (cf. Ref.
Foster:2017hbq ). In Fig. 1 (Left), we show the growth of the magnetization
for various axion masses, computed directly from Eq. (9). The three regimes
($t<\tau_{a}$, $\tau_{a}<t<T_{2}$, $T_{2}<t$) can be clearly observed.

To compute the experimental sensitivity to a highly coherent axion signal, it
is convenient to move to the frequency domain. We imagine a dataset
$\\{M_{n}=M_{x}(n\Delta t)\\}$ of measurements of the magnetization collected
at a frequency $1/\Delta t$, for an integration time $T=N\Delta t$. We can
then compute the power spectral density (PSD) as 444In Eq. (15) we specify the
expectation value of the PSD, rather than the PSD itself. Given a single
experimental dataset, the average PSD cannot be computed, and it will vary
between realizations as the axion field which gave rise to it is stochastic.
Nevertheless, our goal is to determine what value we expect to generate, and
so we will compute the expected PSD throughout. A more careful distinction is
provided in the SM,

$P_{k}=\frac{\Delta t^{2}}{T}\langle|\tilde{M}_{k}|^{2}\rangle,\hskip
14.22636pt\tilde{M}_{k}\equiv\sum_{n=0}^{N-1}e^{-i2\pi kn/N}M_{n}.$ (15)

We can compute the PSD exactly for arbitrary $T$, $T_{2}$, and $\tau_{a}$
555This is discussed in detail in the SM. Also see Lisanti:2021vij ;
Gramolin:2021mqv ; Lee:2022vvb ; Gao:2022nuq  for related considerations.. The
result is particularly simple for $T\gg T_{2}$, taking the form

$P_{k}^{a}\simeq\frac{\langle A^{2}\rangle
T_{2}^{2}}{16\omega_{0}^{2}}\left\\{\begin{array}[]{lc}\displaystyle\frac{4}{\Delta\omega^{2}_{k}T}\sin^{2}\left[\tfrac{1}{2}\Delta\omega_{k}T\right]&T\ll\tau_{a},\\\\[11.38092pt]
\displaystyle\frac{\tau_{a}}{1+\Delta\omega^{2}_{k}T_{2}^{2}}&T\gg\tau_{a},\end{array}\right.$
(16)

with $\Delta\omega_{k}\equiv 2\pi k/T-\omega_{0}$. Due to the resonant
response of the sample, the signal peaks for $\Delta\omega_{k}\simeq 0$. For
$T_{2}\ll T\ll\tau_{a}$, the signal falls dominantly in a single $k$-bin, or
exactly if $\omega_{0}T/2\pi\in\mathbb{N}$. Once $T>\tau_{a},\,T_{2}$, the
signal becomes resolved into multiple bins. This will impact the signal-to-
noise scaling, as we will discuss after considering the relevant background
contributions.

Figure 1: Projections for the sensitivity of a prototypical spin-precession
experiment to the magnetic dipole moment operator,
$g_{N}(\partial_{\mu}a)\bar{N}\gamma^{\mu}\gamma_{5}N$. Left: The axion-
induced growth of the magnetization determined from the Bloch equations. At
short times the magnetization grows $\propto t$, saturating at $T_{2}$. If
$\tau_{a}<T_{2}$, the amplitude grows $\propto\sqrt{t}$ for
$\tau_{a}<t<T_{2}$. The magnitudes are compared to those from a SQUID and
spin-projection noise, evaluated at $k^{*}\equiv\omega_{0}T/2\pi$, and
assuming an integration time $T=T_{2}$. Right: Projected sensitivity to
$g_{N}$ (solid) in comparison to previous projections of Ref. Graham:2013gfa
(dashed), for an almost identical scan strategy. Strengthened sensitivity for
$m_{a}\gg 2\pi/v^{2}T_{2}$ arises primarily due to the signal growth we
account for when $\tau_{a}<t<T_{2}$, whereas the suppression at lower masses
arises partially from a refined estimate of the spin-projection noise. Blue
and red lines correspond to Xenon and Helium targets, and for the latter we
label a “QCD axion targeted” for the results that could be obtained with
negligible spin-projection noise and five years of integration time.

Noise Sources. The axion signal must be detected on top of three relevant
Gaussian background contributions: thermal noise, SQUID noise, and spin-
projection noise 666In both the fluctuating phase and plane wave models, the
axion field itself forms a Gaussian random field, and hence one can show so
too does any linear functional of the axion (see, e.g., Ref. Weinberg:2008zzc
, Appendix E).. Thermal noise arises in the readout circuit and can be
suppressed by cooling the apparatus, and we assume this can be done
sufficiently for this noise source to be neglected (for details, see Ref.
Aybas:2021cdk ). The transverse magnetic field, related to the transverse
magnetization by an ${\cal O}(1)$ factor which we take to be unity (see the
SM), is read out using a SQUID. This magnetometer noise is frequency-
independent for $f\gtrsim 10~{}{\rm Hz}$ PhysRevLett.110.147002 ; Kahn:2016aff
and can be modeled as,

$C_{\rm SQ}(t,t^{\prime})=\delta(t-t^{\prime})\frac{1}{A_{\rm
eff}^{2}}P_{\Phi\Phi}^{\rm SQ}\hskip 2.84544pt\Rightarrow\hskip
2.84544ptP_{k}^{\rm SQ}=\frac{1}{A_{\rm eff}^{2}}P_{\Phi\Phi}^{\rm SQ},$ (17)

where by default, we take $P_{\Phi\Phi}^{\rm SQ}\simeq(\mu\Phi_{0})^{2}/{\rm
Hz}$, with $\Phi_{0}$ the magnetic flux quantum, and $A_{\rm eff}$ is the
effective area of the sample sensed by a pickup loop, which we take to be
$\simeq 0.3~{}{\rm cm}^{2}$, following Ref. Budker:2013hfa .

The final background source we consider is spin-projection noise, which
originates directly from the quantum nature of the nuclear spins in the sample
PhysRevB.75.085310 . We can determine its magnitude from the autocorrelation
function. Consider two successive measurements of the spin operator along the
$\hat{\bf x}$ direction, $S_{x}$, taken at $t$ and then $t^{\prime}$. The
operators are related through the time evolution operator,
$\mathcal{U}(t)=\exp(i\omega_{0}S_{z}t)$,
$S_{x}(t^{\prime})=\mathcal{U}^{\dagger}(t^{\prime}-t)S_{x}(t)\,\mathcal{U}(t^{\prime}-t).$
The magnetization can be determined from the sum over all nuclear spins, and
assuming the sample is hyperpolarized (i.e. unit polarization fraction), we
obtain,

$\displaystyle C_{\rm SP}(t,t^{\prime})$
$\displaystyle=\frac{\gamma^{2}}{2V^{2}}e^{-|t-t^{\prime}|/T_{2}}\sum_{i}\big{\langle}S_{x}^{(i)}(t)S_{x}^{(i)}(t^{\prime})\big{\rangle}+{\rm
h.c.}$
$\displaystyle=\frac{\gamma^{2}nJ}{2V}e^{-|t-t^{\prime}|/T_{2}}\cos\left[\omega_{0}(t-t^{\prime})\right]\\!,$
(18)

where $J$ is the nuclear spin, $n$ is the number density of spins, and $V$ is
the volume of the sample. The exponential factor is included to account for
transverse-spin relaxation, and for further details, see the SM. Eq. (18)
exhibits a $V^{-1}$ scaling, which suggests that for a large enough sample,
this noise source can be suppressed. The corresponding PSD for
$t,t^{\prime}\gg T_{2}$ is given by,

$P_{k}^{\rm
SP}=\frac{\gamma^{2}nJ}{2V}\frac{T_{2}}{1+\Delta\omega^{2}_{k}T_{2}^{2}}\,.$
(19)

This result determines the spin-projection noise for an arbitrary $J$,
assuming a hyperpolarized sample.

Experimental Sensitivity. We now combine the above results to forecast the
expected sensitivity to an axion-induced signal. We use the signal PSD in the
$\tau_{a}\gg T$ and $\tau_{a}\ll T$ limits and compare it to the background
PSD, using all $k$-bins and the likelihood framework presented in the SM
(employing the formalism of Ref. Foster:2017hbq  and inserting the Asimov
dataset Cowan:2010js ). Above the transition region of $\tau_{a}=T$, we
interpolate between the two regimes with a horizontal line. In principle, one
could extend our analysis to handle the intermediate regime. Detailed
projections require us to specify explicit experimental parameters. Even
before this, we can determine the sensitivity scaling with integration time,
which varies depending on the hierarchy between $\tau_{a}$, $T_{2}$, and $T$.
Specifically,

$\displaystyle T\ll\tau_{a},\,T_{2}$ $\displaystyle\hskip
5.69046pt\Rightarrow\hskip 5.69046ptg\propto T^{-3/2}$ (20)
$\displaystyle\tau_{a}\ll T\ll T_{2}$ $\displaystyle\hskip
5.69046pt\Rightarrow\hskip 5.69046ptg\propto T^{-1}$ $\displaystyle T_{2}\ll
T\ll\tau_{a}$ $\displaystyle\hskip 5.69046pt\Rightarrow\hskip
5.69046ptg\propto T^{-1/2}$ $\displaystyle\tau_{a},\,T_{2}\ll T$
$\displaystyle\hskip 5.69046pt\Rightarrow\hskip 5.69046ptg\propto T^{-1/4},$

where $g=g_{N},g_{d}$. The growth becomes slowest once $T$ is the largest
timescale, and the signal is resolved into multiple bins. These scalings are
derived in the SM, however, they arise from comparing the growth of the signal
and background. For instance, for $T_{2}\ll T\ll\tau_{a}$, where the signal is
dominantly in a single bin, $k^{*}\equiv\omega^{0}T/2\pi$, we see that
$P_{k^{*}}^{a}\propto T$ from Eq. (16), whereas $P_{k^{*}}^{\rm SQ}$ and
$P_{k^{*}}^{\rm SP}$ are independent of $T$. Estimating sensitivity by
matching the signal to the background, we find the limit scales as $g\propto
T^{-1/2}$.

To provide quantitative projections, we match the parameters specified in the
CASPEr papers Graham:2011qk ; Budker:2013hfa , focusing on the magnetic dipole
operator. The accessible Larmor frequencies set the mass range we consider,
and we take $10^{-14}~{}{\rm eV}<\omega_{0}<\gamma B_{\rm max}$, with $B_{\rm
max}$ the maximum magnetic field. We fix $T_{2}=100~{}{\rm sec}$ and for each
mass, adopt a variable integration time, $T={\rm max}[\tau_{a},T_{2}]$, to
ensure we always run until the $T^{-1/4}$ growth sets in from Eq. (20). This
implies that the signal, and therefore the likelihood, remains dominated by a
single frequency bin. To ensure each mass is covered only once, resonant
frequencies are adjusted by $2\pi\,{\rm max}[\tau_{a}^{-1},T_{2}^{-1}]$.

We consider two different spin-1/2 samples, and both assumed to be
hyperpolarized. The first consists of pure Xenon-129, which has a nuclear
magnetic moment of $0.78~{}\mu_{N}$ and a nuclear spin density of $1.3\times
10^{22}~{}{\rm cm}^{-3}$, and we assume $B_{\rm max}=10~{}{\rm T}$. The second
is a more optimistic projection using Helium-3 – $\mu=2.12~{}\mu_{N}$ and
$n=2.8\times 10^{22}~{}{\rm cm}^{-3}$ – and $B_{\rm max}=20~{}{\rm T}$, as
well as assuming the SQUID noise can be decreased by two orders of magnitude
below that discussed around Eq. (17). As shown in the SM, to cover the full
mass range, these two experiments would require a total integration time of
56.1 and 61.5 years, respectively.

Our results are shown in Fig. 1 (Right), and contrasted with the projections
of Ref. Graham:2013gfa . The discrepancy has at least three sources: 1. the
additional growth the signal we have accounted for when $\tau_{a}<T<T_{2}$,
see Eq. (14); 2. a different treatment of the spin-projection noise; and 3. a
different calculation of the Larmor frequency given $B_{\rm max}$. 777For
Helium our result is exactly a factor of two larger, whereas for Xenon the
difference is slightly larger than a factor of 4. For Xenon, this difference
is attributed to the use of a combination of isotopes, whereas we considered a
pure sample of Xenon-129. For the spin-projection noise, we compute the value
for each $k$ with Eq. (19) – recall the signal is dominantly peaked in a
single bin – rather than integrating a result over a range of frequencies near
$\omega_{0}$, as in Eq. (A2) in Ref. Budker:2013hfa . Further discussion is
provided in the SM, where we also contrast our results for the electric dipole
operator. As a final benchmark, in Fig. 1, we also show the Helium-3
sensitivity assuming that spin-projection noise could be evaded, and the mass
range is set by assuming five years of integration time. This benchmark cuts
into the QCD axion parameter space, illustrated by the yellow band 888The band
is defined as the region bounded from below by the KSVZ axion model and from
above by the DFSZ model with the down-type Higgs doublet vacuum expectation
value set to zero. (Note that the full range of DFSZ Higgs’ vacuum expectation
values allows the QCD axion parameter space to fill the entire lower right
corner of Fig. 1.).

Discussion. In this work we have derived the axion-induced signal in spin-
precession experiments. These instruments remain one of the most promising
paths to measuring axion-induced magnetic and electric dipoles, and our
calculations demonstrate that their sensitivity is significantly different to
what has previously been estimated. Arguably our most important finding is the
enhanced detection prospects for the QCD axion at high masses, a result which
arises from the continued growth of the axion-induced signal when integrating
beyond the axion coherence time 999In light of our results, one may wonder if
its possible to build an experiment using an element that carries a larger
gyromagnetic ratio than Helium-3 to increase the maximum testable axion mass
and improve the coverage of the QCD axion line. Surveying all measured
gyromagnetic moments of nuclear isotopes, the only element with a
substantially larger value is Thallium-200, although since it carries a half-
life of $\sim$1 day it does not appear viable. Tritium, atomic Hydrogen, and
Fluorine-19 have larger gyromagnetic ratios than Helium-3, but by less than
40%.. Our findings have broader implications. To name one, they demonstrate
that signals that are less coherent than dark matter – for instance, a Cosmic
axion Background Dror:2021nyr  – are more detectable with spin-precession
instruments than may otherwise have been concluded CaB-Nuclear:upcoming .

Acknowledgements. Our work benefited from conversations with Hendrik Bekker,
Yonatan Kahn, Alexander Sushkov, and Arne Wickenbrock. Further, we thank the
anonymous referee for useful feedback. J.M.L. is supported by the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy - EXC 2121 “Quantum
Universe” - 390833306. The research of JD and SG is supported in part by NSF
CAREER grant PHY-1915852 and in part by the U.S. Department of Energy grant
number DE-SC0023093. Part of this work was performed at the Aspen Center for
Physics, which is supported by National Science Foundation grant PHY-1607611.

## References

  * (1) R. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38 (1977) 1440–1443.
  * (2) R. Peccei and H. R. Quinn, Constraints Imposed by CP Conservation in the Presence of Instantons, Phys. Rev. D 16 (1977) 1791–1797.
  * (3) S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40 (1978) 223–226.
  * (4) F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279–282.
  * (5) P. Svrcek and E. Witten, Axions In String Theory, JHEP 06 (2006) 051, [hep-th/0605206].
  * (6) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String Axiverse, Phys. Rev. D 81 (2010) 123530, [arXiv:0905.4720].
  * (7) J. Halverson, C. Long, B. Nelson, and G. Salinas, Towards string theory expectations for photon couplings to axionlike particles, Phys. Rev. D 100 (2019), no. 10 106010, [arXiv:1909.05257].
  * (8) L. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120 (1983) 133–136.
  * (9) J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120 (1983) 127–132.
  * (10) M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120 (1983) 137–141.
  * (11) M. Pospelov and A. Ritz, Theta vacua, QCD sum rules, and the neutron electric dipole moment, Nucl. Phys. B 573 (2000) 177–200, [hep-ph/9908508].
  * (12) P. W. Graham and S. Rajendran, New Observables for Direct Detection of Axion Dark Matter, Phys. Rev. D 88 (2013) 035023, [arXiv:1306.6088].
  * (13) D. Budker, P. W. Graham, M. Ledbetter, S. Rajendran, and A. Sushkov, Proposal for a Cosmic Axion Spin Precession Experiment (CASPEr), Phys. Rev. X 4 (2014), no. 2 021030, [arXiv:1306.6089].
  * (14) Y. V. Stadnik and V. V. Flambaum, Axion-induced effects in atoms, molecules, and nuclei: Parity nonconservation, anapole moments, electric dipole moments, and spin-gravity and spin-axion momentum couplings, Phys. Rev. D 89 (2014), no. 4 043522, [arXiv:1312.6667].
  * (15) C. Abel et al., Search for Axionlike Dark Matter through Nuclear Spin Precession in Electric and Magnetic Fields, Phys. Rev. X 7 (2017), no. 4 041034, [arXiv:1708.06367].
  * (16) T. Wang, D. F. Jackson Kimball, A. O. Sushkov, D. Aybas, J. W. Blanchard, G. Centers, S. R. O’ Kelley, A. Wickenbrock, J. Fang, and D. Budker, Application of Spin-Exchange Relaxation-Free Magnetometry to the Cosmic Axion Spin Precession Experiment, Phys. Dark Univ. 19 (2018) 27–35, [arXiv:1701.08082].
  * (17) A. Garcon et al., Constraints on bosonic dark matter from ultralow-field nuclear magnetic resonance, Sci. Adv. 5 (2019), no. 10 eaax4539, [arXiv:1902.04644].
  * (18) C. Smorra et al., Direct limits on the interaction of antiprotons with axion-like dark matter, Nature 575 (2019), no. 7782 310–314, [arXiv:2006.00255].
  * (19) T. S. Roussy et al., Experimental Constraint on Axionlike Particles over Seven Orders of Magnitude in Mass, Phys. Rev. Lett. 126 (2021), no. 17 171301, [arXiv:2006.15787].
  * (20) M. Jiang, H. Su, A. Garcon, X. Peng, and D. Budker, Search for axion-like dark matter with spin-based amplifiers, Nature Phys. 17 (2021), no. 12 1402–1407, [arXiv:2102.01448].
  * (21) D. Aybas et al., Search for Axionlike Dark Matter Using Solid-State Nuclear Magnetic Resonance, Phys. Rev. Lett. 126 (2021), no. 14 141802, [arXiv:2101.01241].
  * (22) 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].
  * (23) D. Aybas, H. Bekker, J. W. Blanchard, D. Budker, G. P. Centers, N. L. Figueroa, A. V. Gramolin, D. F. J. Kimball, A. Wickenbrock, and A. O. Sushkov, Quantum sensitivity limits of nuclear magnetic resonance experiments searching for new fundamental physics, Quantum Sci. Technol. 6 (2021), no. 3 034007, [arXiv:2103.06284].
  * (24) P. W. Graham and S. Rajendran, Axion Dark Matter Detection with Cold Molecules, Phys. Rev. D 84 (2011) 055013, [arXiv:1101.2691].
  * (25) P. Sikivie, Axion Dark Matter Detection using Atomic Transitions, Phys. Rev. Lett. 113 (2014), no. 20 201301, [arXiv:1409.2806]. [Erratum: Phys.Rev.Lett. 125, 029901 (2020)].
  * (26) P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, W. A. Terrano, L. Trahms, and T. Wilkason, Spin Precession Experiments for Light Axionic Dark Matter, Phys. Rev. D 97 (2018), no. 5 055006, [arXiv:1709.07852].
  * (27) S. P. Chang, S. Haciomeroglu, O. Kim, S. Lee, S. Park, and Y. K. Semertzidis, Axionlike dark matter search using the storage ring EDM method, Phys. Rev. D 99 (2019), no. 8 083002, [arXiv:1710.05271].
  * (28) S. Chang, O. Kim, Y. Semertzidis, S. Haciomeroglu, S. Lee, and S. Park, Axion Searches with the Storage Ring EDM Method, PoS ICHEP2018 (2019) 842.
  * (29) T. Wu et al., Search for Axionlike Dark Matter with a Liquid-State Nuclear Spin Comagnetometer, Phys. Rev. Lett. 122 (2019), no. 19 191302, [arXiv:1901.10843].
  * (30) W. A. Terrano, E. G. Adelberger, C. A. Hagedorn, and B. R. Heckel, Constraints on axionlike dark matter with masses down to $10^{-23}$ eV/c2, Phys. Rev. Lett. 122 (2019), no. 23 231301, [arXiv:1902.04246].
  * (31) 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].
  * (32) Q. Yang, Probe Dark Matter Axions using the Hyperfine Structure Splitting of Hydrogen Atoms, arXiv:1912.11472.
  * (33) F. Abusaif et al., Storage Ring to Search for Electric Dipole Moments of Charged Particles - Feasibility Study, arXiv:1912.07881.
  * (34) V. V. Flambaum, D. Budker, and A. Wickenbrock, Oscillating nuclear electric dipole moments inside atoms, arXiv:1909.04970.
  * (35) JEDI Collaboration, E. Stephenson, A Search for Axion-like Particles with a Horizontally Polarized Beam In a Storage Ring, PoS PSTP2019 (2020) 018.
  * (36) A. V. Gramolin, D. Aybas, D. Johnson, J. Adam, and A. O. Sushkov, Search for axion-like dark matter with ferromagnets, Nature Phys. 17 (2021), no. 1 79–84, [arXiv:2003.03348].
  * (37) 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 (2021), no. 5 055010, [arXiv:2005.11867].
  * (38) K. Gaul, M. G. Kozlov, T. A. Isaev, and R. Berger, Chiral Molecules as Sensitive Probes for Direct Detection of P -Odd Cosmic Fields, Phys. Rev. Lett. 125 (2020), no. 12 123004, [arXiv:2005.02429].
  * (39) A. Arvanitaki, A. Madden, and K. Van Tilburg, The Piezoaxionic Effect, arXiv:2112.11466.
  * (40) O. Kim and Y. K. Semertzidis, New method of probing an oscillating EDM induced by axionlike dark matter using an rf Wien filter in storage rings, Phys. Rev. D 104 (2021), no. 9 096006, [arXiv:2105.06655].
  * (41) 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, arXiv:2208.14454.
  * (42) A. Berlin and K. Zhou, Discovering QCD-Coupled Axion Dark Matter with Polarization Haloscopes, arXiv:2209.12901.
  * (43) M. Lisanti, M. Moschella, and W. Terrano, Stochastic properties of ultralight scalar field gradients, Phys. Rev. D 104 (2021), no. 5 055037, [arXiv:2107.10260].
  * (44) J. Lee, M. Lisanti, W. A. Terrano, and M. Romalis, Laboratory Constraints on the Neutron-Spin Coupling of feV-scale Axions, arXiv:2209.03289.
  * (45) This is an excellent approximation. If we take $T_{2},\tau_{a}\to\infty$, then $M_{x}(t)\sim g_{N}v\sqrt{2\rho_{\scriptscriptstyle\textrm{DM}}}M_{0}\mathop{sin}\nolimits(\omega_{0}t)t$ for the MD interaction, see Eq. (8). This suggests that with enough time an infinite magnetization can be generated, which is inconsistent with the finite number of spins in the system. Terms higher order in $g_{N}$ prevent this from occurring. We can estimate when they must enter by determining when $M_{x}(t)\sim M_{0}$. Setting $g_{N}$ to the largest allowed value (set by the SN 1987A bound), we find $t\sim 30~{}{\rm years}$, dramatically larger than $T_{1}$ for any sample we consider, which represents the longest time we can interrogate the sample for. Including a finite $T_{2}$ or $\tau_{a}$ only strengthens this conclusion.
  * (46) The solution in Eq. (6) is for the Bloch equation before terms of $\mathcal{O}(1/\omega_{0}T_{2})$ are dropped; cf. Eq. (5).
  * (47) The coherence time is a measure of how long the axion field can be approximated as a perfectly coherent driving force, or in the frequency domain how long of a measurement is required to resolve the intrinsic width of the axion. It is an approximate notion, as the transition from perfectly coherent to incoherent happens gradually near $t\sim 2\pi/m_{a}v^{2}$, not exactly at $\tau_{a}$. For further discussion of the coherence time, see for instance Ref. Foster:2020fln .
  * (48) J. W. Foster, N. L. Rodd, and B. R. Safdi, Revealing the Dark Matter Halo with Axion Direct Detection, Phys. Rev. D 97 (2018), no. 12 123006, [arXiv:1711.10489].
  * (49) In Eq. (15) we specify the expectation value of the PSD, rather than the PSD itself. Given a single experimental dataset, the average PSD cannot be computed, and it will vary between realizations as the axion field which gave rise to it is stochastic. Nevertheless, our goal is to determine what value we expect to generate, and so we will compute the expected PSD throughout. A more careful distinction is provided in the SM.
  * (50) This is discussed in detail in the SM. Also see Lisanti:2021vij ; Gramolin:2021mqv ; Lee:2022vvb ; Gao:2022nuq  for related considerations.
  * (51) In both the fluctuating phase and plane wave models, the axion field itself forms a Gaussian random field, and hence one can show so too does any linear functional of the axion (see, e.g., Ref. Weinberg:2008zzc , Appendix E).
  * (52) S. M. Anton, J. S. Birenbaum, S. R. O’Kelley, V. Bolkhovsky, D. A. Braje, G. Fitch, M. Neeley, G. C. Hilton, H.-M. Cho, K. D. Irwin, F. C. Wellstood, W. D. Oliver, A. Shnirman, and J. Clarke, Magnetic flux noise in dc squids: Temperature and geometry dependence, Phys. Rev. Lett. 110 (Apr, 2013) 147002.
  * (53) Y. Kahn, B. R. Safdi, and J. Thaler, Broadband and Resonant Approaches to Axion Dark Matter Detection, Phys. Rev. Lett. 117 (2016), no. 14 141801, [arXiv:1602.01086].
  * (54) M. Braun and J. König, Faraday-rotation fluctuation spectroscopy with static and oscillating magnetic fields, Phys. Rev. B 75 (Feb, 2007) 085310.
  * (55) 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]. [Erratum: Eur.Phys.J.C 73, 2501 (2013)].
  * (56) For Helium our result is exactly a factor of two larger, whereas for Xenon the difference is slightly larger than a factor of 4. For Xenon, this difference is attributed to the use of a combination of isotopes, whereas we considered a pure sample of Xenon-129.
  * (57) The band is defined as the region bounded from below by the KSVZ axion model and from above by the DFSZ model with the down-type Higgs doublet vacuum expectation value set to zero. (Note that the full range of DFSZ Higgs’ vacuum expectation values allows the QCD axion parameter space to fill the entire lower right corner of Fig. 1.).
  * (58) In light of our results, one may wonder if its possible to build an experiment using an element that carries a larger gyromagnetic ratio than Helium-3 to increase the maximum testable axion mass and improve the coverage of the QCD axion line. Surveying all measured gyromagnetic moments of nuclear isotopes, the only element with a substantially larger value is Thallium-200, although since it carries a half-life of $\sim$1 day it does not appear viable. Tritium, atomic Hydrogen, and Fluorine-19 have larger gyromagnetic ratios than Helium-3, but by less than 40%.
  * (59) J. A. Dror, H. Murayama, and N. L. Rodd, Cosmic axion background, Phys. Rev. D 103 (2021), no. 11 115004, [arXiv:2101.09287].
  * (60) J. A. Dror, S. Gori, J. M. Leedom, and N. L. Rodd, to appear.
  * (61) A. V. Gramolin, A. Wickenbrock, D. Aybas, H. Bekker, D. Budker, G. P. Centers, N. L. Figueroa, D. F. J. Kimball, and A. O. Sushkov, Spectral signatures of axionlike dark matter, Phys. Rev. D 105 (2022), no. 3 035029, [arXiv:2107.11948].
  * (62) J. W. Foster, Y. Kahn, R. Nguyen, N. L. Rodd, and B. R. Safdi, Dark Matter Interferometry, Phys. Rev. D 103 (2021), no. 7 076018, [arXiv:2009.14201].
  * (63) J. A. Dror, S. Gori, J. M. Leedom, and N. L. Rodd, See Supplemental Material at [URL will be inserted by publisher] for further details., .
  * (64) S. Weinberg, Cosmology.  2008\.

SUPPLEMENTARY MATERIAL
“On the Sensitivity of Spin-Precession Axion Experiments”
Jeff A. Dror, Stefania Gori, Jacob M. Leedom, and Nicholas L. Rodd

##  S-I Electric Dipole Operator

In the main text, we focused on the impact of our results on the magnetic
dipole operator in Eq. (1). To measure the electric dipole operator, Refs.
Graham:2013gfa ; Budker:2013hfa  proposed using a hyperpolarized sampled with
an additional applied electric field. Note that as the electrons adjust
themselves to shield the field, the nucleus experiences a smaller effective
field ${\bf E}^{\ast}$. To minimize the shielding, one needs to carefully
choose the sample and Refs. Graham:2013gfa ; Budker:2013hfa  proposed
Lead-207. While there are significant differences in the instrument required
to measure the electric dipole operator, in terms of the calculation we
presented in the main text, the primary difference is the $A$ coefficient
controlling the magnitude of the driving force (defined in Eq. (7)).
Accordingly, we can readily extend our results to this operator, simply by
making the replacement $g_{N}\to\sqrt{3}g_{d}E^{\ast}/\omega_{0}v$, see Eq.
(11).

The experiment proposed in Refs. Graham:2013gfa ; Budker:2013hfa  has two
phases, with both utilizing a lead sample of diameter $10~{}{\rm cm}$ and
$E^{\ast}=3\times 10^{10}~{}{\rm V/m}$. In the first phase (P1), the
polarization $p$, relaxation time $T_{2}$, and maximum magnetic field $B_{\rm
max}$ are taken to be $10^{-3}$, $1~{}{\rm msec}$, and $10~{}{\rm T}$,
respectively. In the second phase (P2), these same parameters are given values
$p=1$, $T_{2}=1~{}{\rm sec}$, and $B_{\rm max}=20~{}{\rm T}$. For each phase,
we envision scanning from $10^{-14}~{}{\rm eV}$ up to the maximum accessible
Larmor frequency, and as discussed in App. S-VII, the total scan time is 5.41
hours (213 days) for P1 (P2).

Using the expressions provided in the main text, we compute the expected
sensitivity shown as the solid lines in Fig. S-1. The different scaling
between the mass and coupling in Fig. 1 (Right) and Fig. S-1 arises from the
appearance of $\omega_{0}\simeq m_{a}$ in the amplitude of the driving force
for the magnetic dipole operator, see Eq. (7). More importantly, we again see
clear differences to the projections from the earlier results, here taken from
Ref. Budker:2013hfa . Firstly, we find that the maximum mass that can be
obtained is a factor of two larger than the earlier calculation. Next, note
that the spin-projection noise for Lead-207 is subdominant to the background
from the SQUID, and so our different treatment of $P_{k}^{\rm SP}$ has no
impact in this figure (cf. Fig. 1). For P1, aside from the maximum mass, we
find very good agreement. This arises as for $T_{2}=1~{}{\rm msec}$,
$T_{2}<\tau_{a}$ across the entire mass range. Therefore, at no point can the
magnetization enter the additional period of growth we predict for
$\tau_{a}<t<T_{2}$. For P2, $\tau_{a}<T_{2}$ at higher axion masses, implying
a range where we predict additional growth, and indeed our sensitivity
projection begins to parametrically separate from the earlier result.

Figure S-1: Projected sensitivity for the electric dipole operator. The
experimental proposal has two phases corresponding to the different parameters
discussed in the text. For Phase 1, the experiment never reaches the
transition for the stochastic regime while for Phase 2, the transition region
is shown by the vertical dashed line. For this operator the QCD axion line is
well predicted, and forms a strongly motivated target. The bound from SN1987A
is taken from Ref. Graham:2013gfa .

##  S-II Spin-Projection Noise

In this appendix, we derive the expression for the spin-projection noise given
in Eqs. (18) and (19) in the main text. The aim is to determine the
autocorrelation of spins measured in the $x$-direction, given a sample with
$N_{s}$ nuclei of spin $J$ and gyromagnetic ratio $\gamma$, in a uniform
magnetic field ${\bf B}_{0}=B_{0}\hat{\bf z}$. The operator for the total spin
of the system is the sum of the individual nuclear spins,

$S_{z}=\sum_{i=1}^{N_{s}}S_{z}^{(i)}.$ (S-1)

For each nuclei, we use a basis spanned by the total spin and the projection
along the $z$-axis, so that

$S_{z}^{(i)}\ket{i;J,m}=m\ket{i;J,m},\hskip 14.22636ptm=-J,-J+1,\ldots,J.$
(S-2)

Neglecting interactions between nuclei, the Hamiltonian of the system is given
by, $H=-\bm{\mu}\cdot{\bf B}_{0}=-\omega_{0}S_{z}$, leading to a time
evolution operator, $\mathcal{U}(t)=\exp[i\omega_{0}S_{z}t]$. We now follow
Ref. Graham:2013gfa  and assume that the system is hyperpolarized such that
all sites are in the highest $S_{z}^{(i)}$ eigenstate with $m=J$. The state of
the sample is

$\ket{\uparrow}\equiv\bigotimes_{i=1}^{N_{s}}\ket{i;J,J}\\!.$ (S-3)

First, we compute the autocorrelation for a single spin in the sample,
finding,

$\bra{i;J,J}S^{(i)}_{x}(t)S^{(i)}_{x}(t^{\prime})\ket{i;J,J}=\frac{1}{2}Je^{-i\omega_{0}(t-t^{\prime})}.$
(S-4)

The result for the full system follows,

$\bra{\uparrow}S_{x}(t)S_{x}(t^{\prime})\ket{\uparrow}=\sum_{i=1}^{N_{s}}\bra{i;J,J}S^{(i)}_{x}(t)S^{(i)}_{x}(t^{\prime})\ket{i;J,J}=\frac{1}{2}N_{s}Je^{-i\omega_{0}(t-t^{\prime})}.$
(S-5)

We can now derive Eqs. (18) and (19). Firstly, the spin-projection
autocorrelation function is,

$C_{\rm
SP}(t,t^{\prime})=\frac{\gamma^{2}}{2V^{2}}e^{-|t-t^{\prime}|/T_{2}}\bra{\uparrow}S_{x}(t)S_{x}(t^{\prime})+S_{x}(t^{\prime})S_{x}(t)\ket{\uparrow}=\frac{\gamma^{2}nJ}{2V}e^{-|t-t^{\prime}|/T_{2}}\cos\left[\omega_{0}(t-t^{\prime})\right]\,.$
(S-6)

As discussed in the main text, the exponential factor involving $T_{2}$ is
phenomenological and added to account for the transverse spin relaxation. From
here, we can directly compute the PSD, which in the limit of $t,t^{\prime}\gg
T_{2}$ is,

$P_{k}^{\rm
SP}=\frac{\gamma^{2}nJ}{2V}\frac{T_{2}}{1+\Delta\omega^{2}_{k}T_{2}^{2}}\,,$
(S-7)

as shown in the main text. For $J=1/2$ this agrees parametrically with
existing literature PhysRevB.75.085310 ; Graham:2013gfa ; Budker:2013hfa ;
Aybas:2021cdk .

##  S-III Power Spectral Densities

In the main text, we presented expressions for the autocorrelation functions
for the magnetization induced by the axion, as well as for several noise
sources. Nevertheless, our eventual sensitivities were determined in the
frequency domain, and in particular, we made use of the PSD defined in Eq.
(15), as discussed in App. S-IV. Here we briefly expand on the procedure used
to compute $P_{k}$.

To begin with, we assume that $\Delta t$ is sufficiently finely spaced
compared to the scale over which $M_{x}(t)$ varies that we can approximate the
sums in the discrete Fourier transform with the following integrals,

$P_{k}\simeq\frac{1}{T}\int_{0}^{T}dt\int_{0}^{T}dt^{\prime}e^{i\omega_{k}(t-t^{\prime})}C(t,t^{\prime}),$
(S-8)

where the angular frequency for a given $k$ mode is $\omega_{k}\equiv 2\pi
k/T$. To evaluate the integrals, we switch variables to $x=(t-t^{\prime})/2$
and $y=(t+t^{\prime})/2$. Additionally, except for the SQUID noise, the other
expressions for $C(t,t^{\prime})$ that we studied are of the form
$f(x,y)\cos[2\omega_{0}x]$, allowing us to write,

$P_{k}=\frac{2}{T}\int_{0}^{T/2}dx\cos\left[2(\omega_{0}-\omega_{k})x\right]\int_{x}^{T-x}dyf(x,y),$
(S-9)

where we have dropped a rapidly oscillating term and exploited
$f(x,y)=f(-x,y)$.

Following the procedure above, we can exactly evaluate $P_{k}$ for each
expression we consider, and we use these expressions when computing the limits
we show in Figs. 1 and S-1. If we assume
$\tau_{a},\,T,\,T_{2}\gg\omega_{0}^{-1}$, which holds for almost our entire
parameter space, the expressions simplify, and we have

$\displaystyle\lim_{T\ll\tau_{a}}P_{k}^{a}$ $\displaystyle=\frac{\langle
A^{2}\rangle
T_{2}^{2}}{16\omega_{0}^{2}}\frac{1}{\Delta\omega_{k}^{2}T}\frac{1}{1+\Delta\omega_{k}^{2}T_{2}^{2}}\left[4\sin^{2}\left[\tfrac{1}{2}\Delta\omega_{k}T\right]+\left(\sin[\Delta\omega_{k}T]-\left[1-e^{-T/T_{2}}\right]\Delta\omega_{k}T_{2}\right)^{2}-\sin^{2}[\Delta\omega_{k}T]\right]\\!,$
$\displaystyle\lim_{T\gg\tau_{a}}P_{k}^{a}$ $\displaystyle=\frac{\langle
A^{2}\rangle
T_{2}^{2}}{16\omega_{0}^{2}}\frac{\tau_{a}}{1+\Delta\omega^{2}_{k}T_{2}^{2}}\left[1-\frac{T_{2}}{T}\left(\frac{1}{2}\left[1+e^{-2T/T_{2}}\right]+\frac{1-\Delta\omega_{k}^{2}T_{2}^{2}}{1+\Delta\omega_{k}^{2}T_{2}^{2}}\right.\right.$
(S-10) $\displaystyle\hskip
113.81102pt\left.\left.-2e^{-T/T_{2}}\frac{\cos[\Delta\omega_{k}T]-\Delta\omega_{k}T_{2}\sin[\Delta\omega_{k}T]}{1+\Delta\omega^{2}_{k}T^{2}_{2}}\right)\right]\\!,$
$\displaystyle P_{k}^{\rm SP}$
$\displaystyle=\frac{\gamma^{2}nJ}{2V}\frac{T_{2}}{1+\Delta\omega_{k}^{2}T_{2}^{2}}\left[1-\frac{T_{2}}{T}\left(\frac{1-\Delta\omega_{k}^{2}T_{2}^{2}}{1+\Delta\omega_{k}^{2}T_{2}^{2}}-e^{-T/T_{2}}\frac{(1-\Delta\omega_{k}^{2}T_{2}^{2})\cos[\Delta\omega_{k}T]-2\Delta\omega_{k}T_{2}\sin[\Delta\omega_{k}T]}{1+\Delta\omega_{k}^{2}T_{2}^{2}}\right)\right]\\!.$

In the limit where $T\gg T_{2}$, these expressions reduce to those in Eqs.
(16) and (19), and in App. S-IV, we will use these to confirm the scalings in
Eq. (20). We have not restated $P_{k}^{\rm SQ}$, as the exact expression was
already given in Eq. (17).

##  S-IV Likelihood Framework

Here we expand on the likelihood framework used to compute the projected
limits throughout this work. The formalism we employ is similar to that of
Ref. Foster:2017hbq , and we refer there for additional details. The PSD for
each contribution we consider is exponentially distributed, and therefore the
exponential distribution provides the appropriate likelihood. Our dataset,
$d$, is given by the time series of magnetization measurements, $M(n\Delta
t)$, $n=0,1,\ldots,N-1$, taken over the integration time $T=N\Delta t$. From
the experimentally measured values, we can compute the PSD for the data in
each bin, $P_{k}^{(d)}$. If we then model the PSD in each bin by $P_{k}$, the
likelihood can be written as

${\cal L}(d\,|P_{k})=\sum_{k=1}^{N-1}\frac{1}{P_{k}}e^{-P_{k}^{(d)}/P_{k}}.$
(S-11)

We will decompose our PSD model into a signal and background component as
$P_{k}\equiv S_{k}+B_{k}$. The result for the signal is $P_{k}^{a}$, as given
in Eq. (16), whereas for the background, we combine the SQUID and spin-
projection noises as $B_{k}=P_{k}^{\rm SQ}+P_{k}^{\rm SP}$, using Eqs. (17)
and (19). In order to determine the projected 95% exclusion limits, we will
use a log-likelihood ratio test statistic, and further make use of the Asimov
dataset Cowan:2010js  rather than generating Monte Carlo simulations. In
particular, we compute

$q=2\sum_{k=1}^{N-1}\left[\left(1-\frac{B_{k}}{S_{k}+B_{k}}\right)-\ln\left(1+\frac{S_{k}}{B_{k}}\right)\right]\\!,$
(S-12)

and then determine the amplitude of the signal where $q=-2.71$, which will
correspond to the limiting value.

When computing estimated limits, we make use of the full expression in Eq.
(S-12). However, we can gain intuition from approximate versions of this
result. To begin with, as noted in the main body, until $T$ is larger than
$T_{2}$ and $\tau_{a}$, $S_{k}$ will be non-zero in all bins except
$k^{*}\simeq\omega_{0}T/2\pi$. This is the case for our scan strategy when we
never allow $T$ to become parametrically larger. Accordingly, only a single
bin in (S-12) contributes, and the 95% limit occurs when

$\frac{S_{k^{*}}}{B_{k^{*}}}=-\frac{1}{W_{0}(-e^{-2.71/2-1})}-1\equiv\kappa.$
(S-13)

Here $W_{0}$ is the principal branch of the Lambert $W$ function, and from
this we can compute $\kappa\simeq 8.48$. Note this implies that in the single
bin limit, the properties of the exponential distribution imply that we can
only exclude a value of the signal noticeably larger than the background. To
probe $S_{k}\ll B_{k}$, we will need multiple bins, and therefore $T\gg
T_{2},\tau_{a}$, a case we return to below.

From Eq. (S-13), we see that the 95% limit is obtained for
$S_{k^{*}}\simeq\kappa B_{k^{*}}$ whenever $T$ is not the largest timescale.
We can then use this result to confirm the first three scalings in Eq. (20).
In particular, from Eq. (S-10),

$P_{k^{\ast}}^{a}\simeq\frac{\langle
A^{2}\rangle}{16\omega_{0}^{2}}\left\\{\begin{array}[]{lc}\displaystyle\tfrac{1}{4}T^{3}&T\ll\tau_{a},\,T_{2},\\\\[8.5359pt]
\displaystyle\tfrac{1}{3}T^{2}\,\tau_{a}&\tau_{a}\ll T\ll T_{2},\\\\[8.5359pt]
\displaystyle T\,T_{2}^{2}&T_{2}\ll T\ll\tau_{a},\\\\[8.5359pt]
\displaystyle\tau_{a}\,T_{2}^{2}&\tau_{a},\,T_{2}\ll T,\end{array}\right.$
(S-14)

As $\langle A^{2}\rangle\propto g^{2}$, the claimed scalings follow
immediately from the first three results above.

Although we did not consider this case in our scan strategy, let us take
$T\gg\tau_{a},\,T_{2}$, so that the signal becomes resolved into many bins.
From the above, we see that $P_{k}^{a}$ is no longer growing with $T$ in this
limit. As larger values of $T$ correspond to a finer frequency spacing of the
$k$ values in the PSD, for large enough $T$, we can replace the sum in Eq.
(S-12) with an integral,

$q\simeq\frac{T}{\pi}\int_{0}^{\infty}d\omega\left[\left(1-\frac{B(\omega)}{S(\omega)+B(\omega)}\right)-\ln\left(1+\frac{S(\omega)}{B(\omega)}\right)\right]\\!,$
(S-15)

where $S(\omega)=S_{k=\omega T/2\pi}$, and similarly for $B(\omega)$. We note
that, in actuality, the frequency integration will extend to $2\pi(N-1)/T$
rather than $\infty$. As the signal we consider will be strongly peaked near
the resonant frequency, this distinction is not important. Further, we note
that this same approach was used by default in Ref. Foster:2017hbq , as
instruments searching for the axion-photon coupling commonly work in the
regime where $T$ is the largest timescale. If we assume that we can now probe
$S(\omega)\ll B(\omega)$, then the limit is set when

$q\simeq-\frac{T}{2\pi}\int_{0}^{\infty}d\omega\,\frac{S^{2}(\omega)}{B^{2}(\omega)}.$
(S-16)

In particular, as $S^{2}(\omega)\propto g^{4}$, we see that in this limit, the
sensitivity is scaling as $g\propto T^{-1/4}$, as claimed in Eq. (20).

As a final comment, we note that sensitivity to $S\ll B$ – and therefore the
validity of Eq. (S-15) – requires a large number of bins and also $T$ to be
parametrically larger than any other timescale. We can demonstrate this with a
simple example. Consider a signal present in $n_{S}$ bins, with a PSD value of
$S$ in each, and similarly assume a flat background PSD, $B$. Using Eq.
(S-12), we find the contribution to $q$ is the same from every bin, so in
direct analogy to Eq. (S-13) we obtain the following sensitivity

$\frac{S}{B}=-\frac{1}{W_{0}(-e^{-2.71/2n_{S}-1})}-1.$ (S-17)

For $n_{S}=1$, we recover the single bin result of $S/B=8.48$. To obtain
$S/B=1$ requires $n_{S}\simeq 7$, whereas $S/B=0.1$ occurs for $n_{S}\simeq
300$. For Eq. (S-15) to be percent level accurate requires over $1,000$ bins.

##  S-V Plane Wave Model

In this section, we frame the discussion from the main text in the language of
the plane wave model for axion dark matter, thereby providing an alternative
to the fluctuating phase model used in the main body. All results can be
restated in this alternative approach, however, we focus on the magnetic
dipole operator for simplicity. In particular, we obtain the same result for
the magnetization PSD induced by the axion in the limit $\tau_{a}\ll T$,
demonstrating both models can account for the finite axion coherence time. In
this plane wave model, motivated by the macroscopic occupation number of
axions locally, we treat the dark-matter wave as a sum of $N_{a}$ plane waves
Foster:2017hbq ; Foster:2020fln ; Dror:2021nyr ,

$a(t,{\bf
x})=\frac{a_{0}}{\sqrt{N_{a}}}\sum_{i=1}^{N_{a}}\cos\left[\omega_{i}t-{\bf
k}_{i}\cdot{\bf x}+\phi_{i}\right]\\!.$ (S-18)

For each state in the sum, we evaluate a random velocity ${\bf v}_{i}$ and
phase $\phi_{i}\in[0,2\pi)$, where the former is drawn from the underlying
dark-matter velocity distribution, and used to form the frequency
$\omega_{i}=m_{a}(1+|{\bf v}_{i}|^{2}/2)$ and wavenumber ${\bf
k}_{i}=m_{a}{\bf v}_{i}$. In this formalism, we can straightforwardly include
any velocity distribution, $f_{\bf v}({\bf v}),$ that we choose. (Note we add
the subscript ${\bf v}$ to clearly indicate this is the velocity and not speed
distribution.) For simplicity we will here take an isotropic distribution,

$f_{\bf v}(\mathbf{v})=\frac{1}{\pi^{3/2}v^{3}}e^{-|\mathbf{v}|^{2}/v^{2}}\,,$
(S-19)

with $v\sim 10^{-3}$. For the magnetic dipole operator, the axion driving
force arises from the spatial gradient, which we can take to be evaluated at
${\bf x}=0$ and takes the form,

${\bf G}(t)\equiv\nabla a(t)=\frac{a_{0}}{\sqrt{N_{a}}}\sum_{i=1}^{N_{a}}{\bf
k}_{i}\cos\left[\omega_{i}t+\phi_{i}\right]\\!.$ (S-20)

The PSD of ${\bf G}(t)$ has been studied in Refs. Lisanti:2021vij ;
Gramolin:2021mqv ; Lee:2022vvb ; Gao:2022nuq , and we refer there for further
discussion, although we will not require any details from those works here. We
can now use this expression to model the axion driving force in the Bloch
equations, and determine the resulting magnetization, as we did for the
fluctuating phase model in the main body. Focusing on the $\hat{\bf x}$
component of ${\bf M}$, if we assume that $\Delta t\to 0$ while $T=N\Delta t$
remains finite, we obtain

$|\tilde{M}_{k}|^{2}=\frac{4g_{N}^{2}M_{0}^{2}}{(\omega_{0}^{2}-\omega_{k}^{2})^{2}+4\omega_{k}^{2}/T_{2}^{2}}\bigg{(}\omega_{0}^{2}|\tilde{G}^{x}_{k}|^{2}+\omega_{k}^{2}|\tilde{G}^{y}_{k}|^{2}-2\omega_{0}\omega_{k}\text{Im}[\tilde{G}^{x}_{k}(\tilde{G}^{y}_{k})^{*}]\bigg{)}\,,$
(S-21)

where again $\omega_{k}=2\pi k/T$, and $\tilde{G}_{k}^{s}$ is the Fourier
transform of $\hat{\bf s}\cdot{\bf G}(t)$. This result encapsulates much of
the physics described in the main text. The transfer function factor
multiplying the gradient PSDs has a Lorentzian-like shape with resonant
frequency $\omega_{0}$ and approximate width $T_{2}^{-1}$; the larger $T_{2}$,
the narrower the response of the sample. Even when $\tau_{a}\ll T$, if $T\ll
T_{2}$, then the axion is resolved in the frequency domain, but the
magnetization will still lie within a single bin given the narrow sample
response. In the limit where $T\ll\tau_{a}$, the $\tilde{G}^{i}_{k}$ have an
essentially Dirac $\delta$-function support in frequency, and the signal again
lies in a single bin.

Let us now focus on the high mass parameter space where $\tau_{a}\ll T=T_{2}$
and show that the fluctuating phase and plane wave models agree for the signal
PSD up to ${\cal O}(1)$ factors. Assuming an isotropic dark-matter
distribution, near $\omega_{k}\simeq\omega_{0}$, we have

$\langle\lvert\tilde{G}_{k}^{x}\rvert^{2}\rangle=\langle\lvert\tilde{G}_{k}^{y}\rvert^{2}\rangle\simeq\frac{2\pi^{2}}{3}m_{a}a_{0}^{2}v^{3}f_{\bf
v}(v),$ (S-22)

and the third term in Eq. (S-21) averages to zero. Then we find

$P_{k}^{a}=\frac{4\sqrt{\pi}}{3e}\frac{g_{N}^{2}M_{0}^{2}m_{a}a_{0}^{2}T_{2}^{2}}{1+\Delta\omega_{k}^{2}T_{2}^{2}}.$
(S-23)

This matches the second expression in Eq. (16) up to a factor of $\simeq 1.7$.

The plane wave model further provides an alternative explanation for the
$\sqrt{t}$ growth exhibited by the magnetization when $\tau_{a}\ll t\ll
T_{2}$, seen in Eq. (14), and underlying our enhanced sensitivity at higher
masses seen in Fig. 1. To make the result maximally transparent, we take
$T_{2}\to\infty$, as it plays no role in the physics of interest. Further, we
adopt a simplified version of the plane wave model, where ${\bf
k}_{i}=m_{a}v\hat{\bf x}$ for all $i$, so that

$\nabla a=\frac{a_{0}m_{a}v\hat{\bf
x}}{\sqrt{N}}\sum_{i=1}^{N_{a}}\cos\left[\omega_{i}t+\phi_{i}\right]\\!.$
(S-24)

If we take $\omega_{0}=m_{a}$, and write $\omega_{i}=m_{a}+\delta\omega_{i}$,
with $\delta\omega_{i}\ll m_{a}$, the axion-induced magnetization is given by

$M_{x}(t)\simeq
g_{N}a_{0}m_{a}vM_{0}\frac{1}{\sqrt{N_{a}}}\sum_{i=1}^{N_{a}}\,[c_{i}\sin(m_{a}t)+s_{i}\cos(m_{a}t)]\frac{\sin\left[\tfrac{1}{2}\delta\omega_{i}t\right]}{\tfrac{1}{2}\delta\omega_{i}},$
(S-25)

where $c_{i}=\cos\phi_{i}$ and $s_{i}=\sin\phi_{i}$. We see two frequencies
that enter the sums, the oscillations at the resonant frequency $m_{a}$ and
the oscillations at the beat frequency $\delta\omega_{i}$. Accordingly,

$\langle
M_{x}^{2}(t)\rangle\simeq(g_{N}a_{0}m_{a}vM_{0})^{2}\frac{1}{2N_{a}}\sum_{i=1}^{N_{a}}\,\left\langle\left(\frac{\sin\left[\tfrac{1}{2}\delta\omega_{i}t\right]}{\tfrac{1}{2}\delta\omega_{i}}\right)^{2}\right\rangle\\!.$
(S-26)

This result exhibits the relevant behavior. When $t\ll\tau_{a}$, we cannot
resolve the various frequencies of the axion field, and therefore
$\delta\omega_{i}t\ll 1$. We recover the early $M_{x}(t)\propto t$ growth in
this case. Once $t\gg\tau_{a}$, there will be frequencies for which
$\delta\omega_{i}t$ is no longer small, we have passed the timescale of their
associated beat. The magnetization contributed by these terms will turn over,
and quickly they will no longer contribute to the sum. However, if we expect
the frequencies are uniformly distributed, there will still be a fraction of
frequencies of size $\tau_{a}/t$ for which $\delta\omega_{i}t\ll 1$, and which
remain in the linear growth regime. Accordingly, we expect
$M_{x}(t)\propto\sqrt{t\,\tau_{a}}$ growth to begin, exactly as in Eq. (14).
Although this behavior was derived using a simplified model, the full plane
wave model obeys the same scalings, as can be confirmed numerically.

##  S-VI Magnetic Flux from Axion Dark Matter

In our results, we assumed a simple relation between the magnetic field ${\bf
B}$ outside the sample that is measured at the pickup loop and the
magnetization within the sample, in particular, $B_{x}=M_{x}$. In general, an
${\cal O}(1)$ factor will correct this relation between the magnetization and
the magnetic field. In this appendix, we demonstrate this by explicitly
computing that factor assuming a spherical sample.

The vector potential from such a setup is the sum of the contributions from a
bulk and surface magnetization:

${\bf A}(t,{\bf r})=\frac{1}{4\pi}\int d^{3}{\bf
r}^{\prime}\frac{[\nabla^{\prime}\times{\bf M}({\bf r}^{\prime},t)]_{\rm
ret}}{|{\bf r}-{\bf
r}^{\prime}|}+\frac{1}{4\pi}\oint_{S}dS^{\prime}\frac{[{\bf M}({\bf
r}^{\prime},t)]_{\rm ret}\times\hat{\bf n}^{\prime}}{|{\bf r}-{\bf
r}^{\prime}|},$ (S-27)

where $\hat{\bf n}^{\prime}$ is the unit vector perpendicular to the surface,
and we evaluate the magnetization at the retarded time $t_{\rm
ret}=t-\lvert{\bf r}-{\bf r}^{\prime}\rvert$. Since the coherence length of
axion dark matter for the range of masses we consider is much larger than the
characteristic length scale of the experiment, the magnetization the axion
induces in the sample will be spatially uniform. In this case, the first term
in Eq. (S-27) is parametrically suppressed, and the magnetic field induced by
the sample will be sourced primarily by the surface term.

If the pickup loop is placed near the surface of the sample then $m_{a}|{\bf
r}-{\bf r}^{\prime}|\ll 1$ for all points ${\bf r}^{\prime}$ such that,

${\bf A}(t,{\bf r})\simeq\frac{1}{4\pi}\oint_{S}dS^{\prime}\frac{{\bf
M}(t)\times\hat{\bf n}^{\prime}}{|{\bf r}-{\bf r}^{\prime}|}.$ (S-28)

For a sample composed of a sphere of radius $R$, the vector potential and
magnetic field outside the sample are those of a magnetic dipole,

$\displaystyle{\bf A}(t,{\bf r})$ $\displaystyle=\frac{R^{3}}{3}\frac{{\bf
M}(t)\times\hat{\bf r}}{|{\bf r}|^{3}},$ (S-29) $\displaystyle{\bf B}(t,{\bf
r})$ $\displaystyle=\frac{R^{3}}{3|{\bf r}|^{3}}\left[3({\bf
M}(t)\cdot\hat{\bf r})\hat{\bf r}-{\bf M}(t)\right]\\!.$

In particular, if the magnetic field is measured close to the sample, and
aligned such that we measure the field at ${\bf r}=R\,\hat{\bf x}$, then the
measured magnetic field will be proportional to the magnetization:
$B_{x}(t)=(2/3)M_{x}(t)$. The relationship between the two is an ${\cal O}(1)$
value as claimed.

##  S-VII Integration Time and Scan Strategy

As shown in the main text, a resonant response is one of the characteristic
features of the interaction of the axion with a sample of nuclear spins. Given
the narrow bandwidth of the response, in order to cover a range of possible
axion masses, one must adjust the Larmor frequency by varying the external
magnetic field. One then needs to specify a scanning strategy for how long to
integrate at each frequency and how far apart those frequencies should be
placed. Here we will provide further details of the scan strategy adopted in
the main text and, in particular, detail the calculation of the total
integration time. In this discussion, we will continue to ignore the impact of
the longitudinal relaxation time, $T_{1}$. However, after $T_{1}$ the spins
will lose their coherence in the direction proportional to the external
magnetic field, and the sample will need to be repolarized. We defer to Ref.
Budker:2013hfa  for further discussion.

By default, the scan strategy we adopt is to run until the sensitivity scaling
slows to its lowest parametric growth. From Eq. (20), this occurs when $T={\rm
max}[\tau_{a},T_{2}]$. The second ingredient in the search strategy is how to
spread the $\omega_{0}$ values scanned. At high masses, when $\tau_{a}>T_{2}$,
the width of the axion is larger than the resonant response of the sample.
Therefore, to ensure no axion is missed, we need to shift by a mass-dependent
quantity, $\Delta\omega=2\pi/\tau_{a}$. At lower masses, the response of the
instrument is broader than the axion, and therefore we should shift by this
larger value $\Delta\omega=2\pi/T_{2}$. In summary, our integration time and
$\Delta\omega$ both change on either side of $\tau_{a}=T_{2}$.

We can now determine the total integration time required in the high and low
mass regimes. Consider high masses first, where $\tau_{a}<T_{2}$. Here we
integrate by a constant time, but we have a variable $\Delta\omega$. To
describe this, it is convenient to rewrite the coherence time in terms of the
equivalent quality factor for the axion, $\tau_{a}=2\pi/m_{a}v^{2}=2\pi
Q_{a}/m_{a}$. For dark matter, $Q_{a}\sim v^{-2}\sim 10^{6}$, however, this
allows us to generalize the discussion to more general driving sources (see,
e.g., Ref. Dror:2021nyr ). If we scan a frequency range $[\omega_{\rm
min},\omega_{\rm max}]$, then the total number of bins required is

$N_{\rm high}=\frac{\ln(\omega_{\rm max}/\omega_{\rm
min})}{\ln(1+Q_{a}^{-1})},$ (S-30)

implying the high-mass integration time is,

$T_{\rm high}=T_{2}\frac{\ln(\omega_{\rm max}/\omega_{\rm
min})}{\ln(1+Q_{a}^{-1})}.$ (S-31)

At low masses, where $\tau_{a}>T_{2}$, we have a constant $\Delta\omega$, but
a variable $T$. Accordingly, the number of bins to cover a range $[\omega_{\rm
min},\omega_{\rm max}]$ can be immediately written as,

$N_{\rm low}=\frac{T_{2}}{2\pi}(\omega_{\rm max}-\omega_{\rm min}).$ (S-32)

The total integration time can then be written in terms of the digamma
function, $\psi(x)$, as

$T_{\rm low}=\frac{\omega_{\rm
max}T_{2}^{2}}{2\pi}\left[\psi\left(\frac{T_{2}\omega_{\rm
max}}{2\pi}\right)-\psi\left(\frac{T_{2}\omega_{\rm
min}}{2\pi}\right)\right]\\!.$ (S-33)

These expressions can now be used to evaluate the total integration time. For
the magnetic dipole operator, we considered two results in the main body, one
for Xenon and the other for Helium. In both cases, we took $T_{2}=100~{}{\rm
sec}$, implying that in both cases, the transition between the low and high
mass regime occurs at $m\simeq 41.4~{}{\rm peV}$. At the high end, the maximum
frequency varies given the different assumed $B_{\rm max}$ in addition to the
larger Helium-3 magnetic moment. Combining these details, we confirm the
integration times quoted in the main body,

${\rm magnetic~{}dipole\\!:}\hskip 14.22636ptT_{\rm low}+T_{\rm
high}\simeq\left\\{\begin{array}[]{lc}26.4+29.7~{}{\rm years}\simeq
56.1~{}{\rm years}&({\rm Xenon}),\\\\[2.84544pt] 26.4+35.1~{}{\rm years}\simeq
61.5~{}{\rm years}&({\rm Helium}).\end{array}\right.$ (S-34)

These times are significant, although they follow from the strategy suggested
in Ref. Graham:2013gfa . Given that we find the sensitivity is enhanced at
high masses and reduced at lower values, if we instead only run from the
maximum mass down to $1~{}{\rm neV}$ ($10~{}{\rm neV}$), then the total
integration time would reduce to 19.6 years (12.3 years) for Xenon, and 25.0
years (17.7 years) for Helium. We note that the “QCD axion targeted” scan we
show in Fig. 1 only used a decade of integration time.

We can further determine the total integration time required for the electric
dipole operator, considered in App. S-I. As discussed there, we consider the
two phases suggested in Ref. Budker:2013hfa . The first phase has
$T_{2}=1~{}{\rm msec}$, whereas the second has $T_{2}=1~{}{\rm sec}$. The
corresponding masses where $\tau_{a}=T_{2}$ are $4.14~{}\mu{\rm eV}$ and
$4.14~{}{\rm neV}$, respectively. This transition regime for phase one is at
such a large mass that, in fact, the entire integration occurs in the low-mass
regime. The integration times are now,

${\rm electric~{}dipole\\!:}\hskip 14.22636ptT_{\rm low}+T_{\rm
high}\simeq\left\\{\begin{array}[]{lc}5.41+0~{}{\rm hours}\simeq 5.41~{}{\rm
hours}&({\rm Phase}~{}1),\\\\[2.84544pt] 152+60.3~{}{\rm days}\simeq
213~{}{\rm days}&({\rm Phase}~{}2).\end{array}\right.$ (S-35)

The integration times here are significantly shorted than those in Eq. (S-34).
We note, however, that this is assuming integration times of $T={\rm
max}[\tau_{a},T_{2}]$. In Ref. Budker:2013hfa , they claim to use
$T=100~{}{\rm sec}\times{\rm max}[1,\tau_{a}/2T_{2}]$, and if this was
adopted, we find that Phases 1 and 2 would run for roughly 120 and 40 years,
respectively. Nevertheless, we suspect this is not what was adopted in their
results, as with this scan strategy our results in Fig. S-1 would look
significantly different and we would no longer find close agreement for Phase
1.