On the Sensitivity of Spin-Precession Axion Experiments

A leading direction in the hunt for axion dark matter is to search for its influence on nuclear spins. The detection scheme involves polarizing a sample of nuclei within a strong static magnetic field and then looking for a spin precession induced by the oscillating axion field. We study the axion signal and background contributions that arise in such experiments (a prominent example being CASPEr), finding key differences with the existing literature. Most importantly, in the limit where the transverse spin-relaxation time of the material is the largest timescale of the problem, we show that the induced signal continues to grow even beyond the coherence time of the axion field. As a result, 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 at an axion mass of 100 neV using a Xenon-129 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.

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 [1][2][3][4], axions have since been appreciated for both their ubiquity in string theory [5][6][7] and the generic expectation that they contribute to cold dark matter through, for example, misalignment production [8][9][10]. 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 1 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, In the presence of an axion-wave background, these interactions source an oscillating magnetic dipole (MD), weighted by g N ∝ 1/f a , and an oscillating electric dipole (ED), weighted by g d ∝ 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 µνG µν , and is fixed to g d = (3.7 ± 1.5) × 10 −18 GeV −2 (10 15 GeV/f a ) [11]. For a nonrelativistic nucleus, these interactions lead to the Hamiltonian, H int = −2(g N ∇a + g d aE * ) · S. Here, S is the nucleon spin operator, and E * 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 axionnucleus interaction can be characterized as an effective axion magnetic field, with γ being the gyromagnetic ratio of the nucleon. This axion-induced magnetic field can be detected using spinprecession techniques such as nuclear magnetic resonance (NMR), as originally proposed in the seminal CASPEr papers [12,13], and it was further developed experimentally and theoretically in Refs. [14][15][16][17][18][19][20][21][22][23]. 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. ). In this Letter, we revisit the sensitivity of spinprecession 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, τ a ∼ (4 neV/m a ) 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 τ 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, B 0 . The field induces a bulk magnetization M 0 that is parallel to B 0 . A perpendicular magnetic fieldsuch 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 B 0 at the Larmor frequency, ω 0 ≡ γB 0 , generating an oscillating transverse magnetization, which can then be detected with a sensitive magnetometer. For m a ∼ ω 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 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 M 0 ∝ B 0 ∝ẑ, the Bloch equations read The longitudinal relaxation time commonly satisfies T 1 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 B = B 0ẑ + B a (t). As |B a (t)| B 0 , we study the magnetization perturbatively, taking M M 0 + M a , and we will work only to linear order in the axion-induced fields throughout [45]. Working to this order, M z (t) = M 0 , leaving the dynamics to the transverse magnetizations, These equations can be decoupled. If we measure thex component of the magnetization, the equation to solve is Here and throughout, we neglect terms of O(1/ω 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 ω 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 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 = √ 2ρ DM /m a , fixed by the local dark-matter density, and oscillates with frequency ω a m a (1 + v 2 /2). The statistics of the field are then encoded by requiring the field, which carries velocity v ∼ 10 −3 , obtain a new random phase uniformly sampled on [0, 2π), every coher-ence time, τ a = 2π/m a v 2 [47]. Each time the phase is updated, we further update the direction of the axion field's momentum, k (where |k| = m a v), though, parametrically, our results are insensitive to the stochastic nature of the momentum vector. In summary, where we focus on the dominant term in the driving force, which carries a phase ϕ(t), uniformly sampled on [0, 2π), but shifted from the axion phase. Furthermore, we assumed E * ∝x for the ED operator. From Eq. (5), a resonant response is induced in the system . We will assume a scan strategy such that this condition is always satisfied, and thereby assume ω a ω 0 .
Whenever t τ a , the axion behaves as if it were a perfectly coherent driving force: ϕ and k are constant, so that Eq. (6) yields the following oscillating solution For t T 2 , the amplitude grows linearly in time. Beyond T 2 , however, the growth saturates.
We next extend these results to finite τ a . Due to the stochastic variation of ϕ and k for t > τ a , it is useful to compute the autocorrelation function of the induced magnetization, The expectation value vanishes unless the random phases are identical, which requires |t −t | < τ a , so that in terms of the step function Θ. (Here, we have also used that ϕ and k are uncorrelated.) The remaining expectation value can be determined by averaging over the incident direction of the axion, yielding At short times (t, t τ a ) where the axion is continuous, we have (cf. Eq. (8)), At longer times, the stochastic fluctuations in the axion associated with τ a must be accounted for. Combining Eqs. (9) and (10), for t, t τ a we can evaluate the integral by rotating coordinates tot ±t , yielding ×e −(t+t )/T2 e 2min(t,t )/T2 − 1 .
Eqs. (12) and (13) allow us to infer the growth of the magnetization in the presence of a finite τ a , as C(t, t) = M 2 x (t) . Firstly, for t 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 τ a . For t T 2 , however, the behaviors differ, The first result, that the amplitude of the magnetization grows linearly with time for t τ a , T 2 , is consistent with Eq. (8). However, for τ a < t < T 2 , we see that the amplitude continues to grow, leading to an ever-increasing size, albeit as √ t. Intuitively, this transition in behavior can be understood as the magnetization executing a random walk. For t > τ 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 τ a (given the growth for t < τ a ), and with a number of steps t/τ a , so that we expect |M x (t)| ∝ √ tτ a , exactly as found (cf. Ref. [48]). In Fig. 1 (Left), we show the growth of the magnetization for various axion masses, computed directly from Eq. (9). The three regimes (t < τ a , τ 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∆t)} of measurements of the magnetization collected at a frequency 1/∆t, for an integration time T = N ∆t. We can then compute the power spectral density (PSD) as [49], We can compute the PSD exactly for arbitrary T , T 2 , and τ a [50]. The result is particularly simple for T T 2 , taking the form with ∆ω k ≡ 2πk/T − ω 0 . Due to the resonant response of the sample, the signal peaks for ∆ω k 0. For T 2 T τ a , the signal falls dominantly in a single k-bin, or exactly if ω 0 T /2π ∈ N. Once T > τ 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.
Noise Sources. The axion signal must be detected on top of three relevant Gaussian background contributions: thermal noise, SQUID noise, and spin-projection noise [51]. 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. [23]). The transverse magnetic field, related to the transverse magnetization by an 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 10 Hz [52, 53] and can be modeled as, where by default, we take P SQ ΦΦ (µΦ 0 ) 2 /Hz, with Φ 0 the magnetic flux quantum, and A eff is the effective area of the sample sensed by a pickup loop, which we take to be 0.3 cm 2 , following Ref. [13].
The final background source we consider is spinprojection noise, which originates directly from the quantum nature of the nuclear spins in the sample [54]. We can determine its magnitude from the autocorrelation function. Consider two successive measurements of the spin operator along thex direction, S x , taken at t and then t . The operators are related through the time evolution operator, 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, 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 T 2 is given by, This result determines the spin-projection noise for an arbitrary J, assuming a hyperpolarized sample.
Experimental Sensitivity. We now combine the above At short times the magnetization grows ∝ t, saturating at T2. If τa < T2, the amplitude grows ∝ √ t for τa < t < T2. The magnitudes are compared to those from a SQUID and spin-projection noise, evaluated at k * ≡ ω0T /2π, and assuming an integration time T = T2. Right: Projected sensitivity to gN (solid) in comparison to previous projections of Ref. [12] (dashed), for an almost identical scan strategy. Strengthened sensitivity for ma 2π/v 2 T2 arises primarily due to the signal growth we account for when τa < t < T2, 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.
results to forecast the expected sensitivity to an axioninduced signal. We use the signal PSD in the τ a T and τ a 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. [48] and inserting the Asimov dataset [55]). Above the transition region of τ 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 τ a , T 2 , and T . Specifically, 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 T τ a , where the signal is dominantly in a single bin, k * ≡ ω 0 T /2π, we see that P a k * ∝ T from Eq. (16), whereas P SQ k * and P SP k * are independent of T . Estimating sensitivity by matching the signal to the background, we find the limit scales as To provide quantitative projections, we match the parameters specified in the CASPEr papers [13,24], focusing on the magnetic dipole operator. The accessible Larmor frequencies set the mass range we consider, and we take 10 −14 eV < ω 0 < γB max , with B max the maximum magnetic field. We fix T 2 = 100 sec and for each mass, adopt a variable integration time, T = max[τ 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π max[τ −1 a , T −1 2 ]. 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 µ N and a nuclear spin density of 1.3 × 10 22 cm −3 , and we assume B max = 10 T. The second is a more optimistic projection using Helium-3 -µ = 2.12 µ N and n = 2.8 × 10 22 cm −3 -and B max = 20 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. [12]. The discrepancy has at least three sources: 1. the additional growth the signal we have accounted for when τ a < T < T 2 , see Eq. (14); 2. a different treatment of the spin-projection noise; and 3. a different calculation of the Larmor fre-quency given B max .
[56] 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 ω 0 , as in Eq. (A2) in Ref. [13]. 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 spinprojection 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 [57].
Discussion. In this work we have derived the axioninduced 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 [58]. 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 [59] -are more detectable with spin-precession instruments than may otherwise have been concluded [60].  (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 gN prevent this from occurring. We can estimate when they must enter by determining when Mx(t) ∼ M0. Setting gN to the largest allowed value (set by the SN 1987A bound), we find t ∼ 30 years, dramatically larger than T1 for any sample we consider, which represents the longest time we can interrogate the sample for. Including a finite T2 or τa only strengthens this conclusion.

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. [12,13] 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 E * . To minimize the shielding, one needs to carefully choose the sample and Refs. [12,13] 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 → √ 3g d E * /ω 0 v, see Eq. (11). The experiment proposed in Refs. [12,13] has two phases, with both utilizing a lead sample of diameter 10 cm and E * = 3 × 10 10 V/m. In the first phase (P1), the polarization p, relaxation time T 2 , and maximum magnetic field B max are taken to be 10 −3 , 1 msec, and 10 T, respectively. In the second phase (P2), these same parameters are given values p = 1, T 2 = 1 sec, and B max = 20 T. For each phase, we envision scanning from 10 −14 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 ω 0 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. [13]. 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 SP k 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 msec, T 2 < τ a across the entire mass range. Therefore, at no point can the magnetization enter the additional period of growth we predict for τ a < t < T 2 . For P2, τ 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.

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 γ, in a uniform magnetic field B 0 = B 0ẑ . The operator for the total spin of the system is the sum of the individual nuclear spins, For each nuclei, we use a basis spanned by the total spin and the projection along the z-axis, so that Neglecting interactions between nuclei, the Hamiltonian of the system is given by, H = −µ·B 0 = −ω 0 S z , leading to a time evolution operator, U(t) = exp[iω 0 S z t]. We now follow Ref. [12] and assume that the system is hyperpolarized such that all sites are in the highest S First, we compute the autocorrelation for a single spin in the sample, finding, The result for the full system follows, We can now derive Eqs. (18) and (19). Firstly, the spin-projection autocorrelation function is, 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 T 2 is, as shown in the main text. For J = 1/2 this agrees parametrically with existing literature [12,13,23,54].

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 ∆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, where the angular frequency for a given k mode is ω k ≡ 2πk/T . To evaluate the integrals, we switch variables to x = (t − t )/2 and y = (t + t )/2. Additionally, except for the SQUID noise, the other expressions for C(t, t ) that we studied are of the form f (x, y) cos[2ω 0 x], allowing us to write, 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 τ a , T, T 2 ω −1 0 , which holds for almost our entire parameter space, the expressions simplify, and we have In the limit where T 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 SQ k , 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.
[48], 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∆t), n = 0, 1, . . . , N − 1, taken over the integration time T = N ∆t. From the experimentally measured values, we can compute the PSD for the data in each bin, P (d) k . If we then model the PSD in each bin by P k , the likelihood can be written as We will decompose our PSD model into a signal and background component as P k ≡ S k + B k . The result for the signal is P a k , as given in Eq. (16), whereas for the background, we combine the SQUID and spin-projection noises as B k = P SQ k + P SP k , 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 [55] rather than generating Monte Carlo simulations. In particular, we compute 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 τ a , S k will be non-zero in all bins except k * ω 0 T /2π. 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 Here W 0 is the principal branch of the Lambert W function, and from this we can compute κ 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 B k , we will need multiple bins, and therefore T T 2 , τ a , a case we return to below.
From Eq. (S-13), we see that the 95% limit is obtained for S k * κ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), As A 2 ∝ 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 τ a , T 2 , so that the signal becomes resolved into many bins. From the above, we see that P a k 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, where S(ω) = S k=ωT /2π , and similarly for B(ω). We note that, in actuality, the frequency integration will extend to 2π(N − 1)/T rather than ∞. 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.
[48], 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(ω) B(ω), then the limit is set when In particular, as S 2 (ω) ∝ g 4 , we see that in this limit, the sensitivity is scaling as g ∝ T −1/4 , as claimed in Eq. (20).
As a final comment, we note that sensitivity to S 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 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 τ a 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 [48,59,62], For each state in the sum, we evaluate a random velocity v i and phase φ i ∈ [0, 2π), where the former is drawn from the underlying dark-matter velocity distribution, and used to form the frequency ω i = m a (1 + |v i | 2 /2) and wavenumber k i = m a v i . In this formalism, we can straightforwardly include any velocity distribution, f v (v), that we choose. (Note we add the subscript v to clearly indicate this is the velocity and not speed distribution.) For simplicity we will here take an isotropic distribution, with v ∼ 10 −3 . For the magnetic dipole operator, the axion driving force arises from the spatial gradient, which we can take to be evaluated at x = 0 and takes the form, The PSD of G(t) has been studied in Refs. [41,43,44,61], 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 thex component of M, if we assume that ∆t → 0 while T = N ∆t remains finite, we obtain where again ω k = 2πk/T , andG s k is the Fourier transform ofŝ · 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 ω 0 and approximate width T −1 2 ; the larger T 2 , the narrower the response of the sample. Even when τ a T , if T 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 τ a , theG i k have an essentially Dirac δ-function support in frequency, and the signal again lies in a single bin.
Let us now focus on the high mass parameter space where τ a T = T 2 and show that the fluctuating phase and plane wave models agree for the signal PSD up to O(1) factors. Assuming an isotropic dark-matter distribution, near ω k ω 0 , we have -22) and the third term in Eq. (S-21) averages to zero. Then we find This matches the second expression in Eq. (16) up to a factor of 1.7.
The plane wave model further provides an alternative explanation for the √ t growth exhibited by the magnetization when τ a t 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 → ∞, as it plays no role in the physics of interest. Further, we adopt a simplified version of the plane wave model, where k i = m a vx for all i, so that If we take ω 0 = m a , and write ω i = m a + δω i , with δω i m a , the axion-induced magnetization is given by where c i = cos φ i and s i = sin φ i . We see two frequencies that enter the sums, the oscillations at the resonant frequency m a and the oscillations at the beat frequency δω i . Accordingly, This result exhibits the relevant behavior. When t τ a , we cannot resolve the various frequencies of the axion field, and therefore δω i t 1. We recover the early M x (t) ∝ t growth in this case. Once t τ a , there will be frequencies for which δω 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 τ a /t for which δω i t 1, and which remain in the linear growth regime. Accordingly, we expect M x (t) ∝ √ t τ 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 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 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: wheren is the unit vector perpendicular to the surface, and we evaluate the magnetization at the retarded time t ret = t − |r − r |. 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 |r − r | 1 for all points r such that, 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, (S- 29) In particular, if the magnetic field is measured close to the sample, and aligned such that we measure the field at r = Rx, 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 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. [13] 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 = max[τ a , T 2 ]. The second ingredient in the search strategy is how to spread the ω 0 values scanned. At high masses, when τ 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, ∆ω = 2π/τ a . At lower masses, the response of the instrument is broader than the axion, and therefore we should shift by this larger value ∆ω = 2π/T 2 . In summary, our integration time and ∆ω both change on either side of τ a = T 2 .
We can now determine the total integration time required in the high and low mass regimes. Consider high masses first, where τ a < T 2 . Here we integrate by a constant time, but we have a variable ∆ω. To describe this, it is convenient to rewrite the coherence time in terms of the equivalent quality factor for the axion, τ a = 2π/m a v 2 = 2πQ a /m a . For dark matter, Q a ∼ v −2 ∼ 10 6 , however, this allows us to generalize the discussion to more general driving sources (see, e.g., Ref. [59]). If we scan a frequency range [ω min , ω max ], then the total number of bins required is .

(S-31)
At low masses, where τ a > T 2 , we have a constant ∆ω, but a variable T . Accordingly, the number of bins to cover a range [ω min , ω max ] can be immediately written as, (ω max − ω min ). (S-32) The total integration time can then be written in terms of the digamma function, ψ(x), as 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 sec, implying that in both cases, the transition between the low and high mass regime occurs at m 41.4 peV. At the high end, the maximum frequency varies given the different assumed B max in addition to the larger Helium-3 magnetic moment. Combining these details, we confirm the integration times quoted in the main body, These times are significant, although they follow from the strategy suggested in Ref. [12]. 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 neV (10 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. [13]. The first phase has T 2 = 1 msec, whereas the second has T 2 = 1 sec. The corresponding masses where τ a = T 2 are 4.14 µeV and 4.14 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, The integration times here are significantly shorted than those in Eq. (S-34). We note, however, that this is assuming integration times of T = max[τ a , T 2 ]. In Ref. [13], they claim to use T = 100 sec × max[1, τ 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.