# Supplementary Information for:
“Search for axion-like dark matter with spin-based amplifiers”

Min Jiang  These authors contributed equally to this work Hefei National
Laboratory for Physical Sciences at the Microscale and Department of Modern
Physics, University of Science and Technology of China, Hefei 230026, China
CAS Key Laboratory of Microscale Magnetic Resonance, University of Science and
Technology of China, Hefei 230026, China  Synergetic Innovation Center of
Quantum Information and Quantum Physics, University of Science and Technology
of China, Hefei 230026, China  Haowen Su  These authors contributed equally to
this work Hefei National Laboratory for Physical Sciences at the Microscale
and Department of Modern Physics, University of Science and Technology of
China, Hefei 230026, China  CAS Key Laboratory of Microscale Magnetic
Resonance, University of Science and Technology of China, Hefei 230026, China
Synergetic Innovation Center of Quantum Information and Quantum Physics,
University of Science and Technology of China, Hefei 230026, China  Antoine
Garcon  Helmholtz-Institut, GSI Helmholtzzentrum für Schwerionenforschung,
Mainz 55128, Germany  Johannes Gutenberg University, Mainz 55128, Germany
Xinhua Peng<EMAIL_ADDRESS>Hefei National Laboratory for Physical
Sciences at the Microscale and Department of Modern Physics, University of
Science and Technology of China, Hefei 230026, China  CAS Key Laboratory of
Microscale Magnetic Resonance, University of Science and Technology of China,
Hefei 230026, China  Synergetic Innovation Center of Quantum Information and
Quantum Physics, University of Science and Technology of China, Hefei 230026,
China  Dmitry Budker  Helmholtz-Institut, GSI Helmholtzzentrum für
Schwerionenforschung, Mainz 55128, Germany  Johannes Gutenberg University,
Mainz 55128, Germany  Department of Physics, University of California,
Berkeley, CA 94720-7300, USA

###### Contents

  1. I Experimental apparatus
  2. II Analysis of spin-based amplifier
    1. II.1 Linear and nonlinear response
    2. II.2 Amplification factor
    3. II.3 Bandwidth
  3. III Derivation of the effective magnetic field induced by axion-like dark matter
    1. III.1 Axion dark-matter halo
    2. III.2 Coordinate transformation
  4. IV Data analysis for the axion-like dark matter search
    1. IV.1 Calibration experiments
    2. IV.2 Axion-like dark matter search experiment
      1. IV.2.1 The timing of axion-like dark matter searches
      2. IV.2.2 Detection threshold determination
  5. V Constraint plots for other possible bosonic dark matter fields
  6. VI Projected sensitivity reach
  7. VII Testing the data analysis procedure with simulated axion signals

Figure S1: Experimental setup for axion-like dark matter searches. The details
are present in the text. DAQ, data Acquisition; PID,
proportional–integral–derivative controller; AWG, arbitrary wavefunction
generator; $\lambda/4$, quarter-wave plate.

##  I Experimental apparatus

This section presents the details of our experimental setup for axion-like
dark matter searches. Figure S1 shows the schematic of the experimental
apparatus. 0.5 cm3 volume cubic vapor cell made from pyrex glass contains 5
torr isotopically enriched 129Xe, 250 torr N2 as buffer gas, and a droplet
(several miligrams) of isotopically enriched 87Rb metal. The vapor cell is
heated to $140$ ∘C with twisted wire using AC current. The cell is
magnetically shielded with a five-layer cylindrical $\mu$-metal shield
(shielding factor of 106). Two sets of three pairs orthogonal coils are placed
around the vapor cell and can provide bias and oscillating magnetic fields in
an arbitrary direction. A bias field $B_{0}^{z}$ (above 200 $\rm{nT}$) is
applied along $z$ with current driven with a precision current source (Krohn-
Hite Model 523) and an oscillating field $B_{\rm{ac}}^{y}$ is applied along
$y$ with an arbitrary waveform generator (Keysight Model 33210A). The 87Rb
atoms in the vapor cell are pumped with a circularly polarized laser light
along ${z}$. The 129Xe atoms are polarized by spin-exchange collisions with
optically pumped 87Rb atoms. The pump laser light is tuned to the D1
transition at 795 nm and its optical power is $\sim$26 mW. The magnetic field
is measured via optical rotation of the linearly polarized probe laser beam
propagating along the $x$ direction. The probe beam is blue-detuned by 110 GHz
from the D2 transition at 780 nm and its power is $\sim$0.8 mW. The optical
rotation of the probe beam after passing the vapor cell is (see, for example,
refs. [1, 2, 3])

$\theta=\frac{1}{4}lr_{e}cfnP_{x}^{e}D(v),$ (S1)

where $P_{x}^{e}$ is the electron spin polarization of 87Rb atoms along the
$x$ axis; its explicit form (which depends on the unknown magnetic fields) is
described below, $l\approx 8$ mm is the optical path length, $r_{e}=2.8\times
10^{-13}$ cm is the classical radius of the electron, $c$ is the speed of
light, $f$ is the oscillator strength (about 1/3 for D1 light and 2/3 for D2
light), $D(v)=(v-v_{0})/[(v-v_{0})^{2}+(\Delta v/2)^{2}]$, $v$ is the
frequency of the probe laser light, and $\Delta v$ is the full-width at half-
maximum (FWHM) of the optical D2 transition of frequency $v_{0}$. To suppress
the influence of low-frequency noise of the probe beam, the polarization of
the probe beam is modulated with a photoelastic modulator (PEM) at 50
$\rm{kHz}$. The signal is demodulated with a lock-in amplifier (SRS Model 830)
and is acquired with a 24-bit acquisition card (NI 9239).

As discussed below, when an oscillating external magnetic field is resonant
with 129Xe spins, the induced 129Xe transverse magnetization becomes
considerable. The transverse nuclear magnetization can generate an effective
magnetic field on 87Rb atoms and such magnetic field can be detected with a
87Rb magnetometer (the details are present in Sec. II). We show below that the
amplitude of the effective field could be significantly larger than that of
the external oscillating field. Thus, 129Xe spins act as a spin-based
amplifier, see Sec. II, to boost the signal from the external oscillating
magnetic field. Using the 129Xe spin-based amplifier, the magnetic sensitivity
of the 87Rb magnetometer can be enhanced by two orders of magnitude and can
reach 18.5 fT/$\sqrt{\rm{Hz}}$. This is in contrast to other state-of-the-art
87Rb magnetometers based on 87Rb-129Xe vapor cells, where their magnetic
sensitivities are typically on the order of pT/$\sqrt{\rm{Hz}}$. Benefiting
from the high sensitivity of spin-based amplifier, we use this device to
search for axion-like dark matter, which can be seen as a coherently
oscillating magnetic field (see Sec. IV.2).

##  II Analysis of spin-based amplifier

This section presents a detailed analysis of our spin-based amplifier. The
deveice makes use of 129Xe noble gas, which spatially overlaps with 87Rb gas
in the same vapor cell. The spin dynamics of the overlapping 129Xe and 87Rb
can be described by the coupled Bloch equations [4],

$\displaystyle\frac{\partial\textbf{P}^{e}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{e}}{Q}(\textbf{B}^{0}+\lambda{M}^{n}\textbf{P}^{n})\times\textbf{P}^{e}+\frac{P^{e}_{0}\bm{z}-\textbf{P}^{e}}{\\{T_{2e},T_{2e},T_{1e}\\}Q},$
(S2) $\displaystyle\frac{\partial\textbf{P}^{n}}{\partial t}$ $\displaystyle=$
$\displaystyle\gamma_{n}(\textbf{B}^{0}+\lambda{M}^{e}\textbf{P}^{e})\times\textbf{P}^{n}+\frac{P^{n}_{0}\bm{z}-\textbf{P}^{n}}{\\{T_{2n},T_{2n},T_{1n}\\}},$
(S3)

where $\gamma_{e}=2\pi\times 28~{}\textrm{Hz/nT}$ and $\gamma_{n}=2\pi\times
0.0118~{}\textrm{Hz/nT}$ are, respectively, the gyromagnetic ratio of a bare
electron and 129Xe nuclear spin, $\textbf{P}^{e}$ and $\textbf{P}^{n}$ are,
respectively, the polarization of 87Rb electrons and 129Xe atoms, $P_{0}^{e}$
and $P_{0}^{n}$ are the equilibrium polarization of 87Rb and 129Xe, $M^{e}$
and $M^{n}$ are, respectively, the magnetization of 87Rb atoms and 129Xe atoms
with unity polarization, $T_{1n},T_{2n}\approx 20$ s are the longitudinal and
transverse relaxation times of 129Xe spins, $\textbf{B}^{0}$ is the applied
magnetic field, and $\lambda M\textbf{P}$ is the effective field induced by P
polarization. $T_{1e},T_{2e}$ are the longitudinal and transverse relaxation
times of 87Rb spins. The relaxation times of 87Rb spins satisfy
$1/T_{2e}=1/T_{1e}+1/T_{2}^{\rm{SE}}$, where $1/T_{2}^{\rm{SE}}$ is the
relaxation caused by spin-exchange collisions. In our experiment, $1/T_{1e}$
is the dominant relaxation term. Accordingly, $T_{1e}\approx T_{2e}$ can be
regarded as the common relaxation time without distinction here. The factor
$Q$ is the slowing-down factor of 87Rb atoms, which depends on the 87Rb
polarization. Because $P_{z}^{e}$ is much larger than transverse $P_{x}^{e}$
and $P_{y}^{e}$, $Q$ primarily depends on $P_{z}^{e}$. Moreover, as
$P_{z}^{e}$ is approximated as a constant, $Q$ can be approximated as a
constant.

For spatially overlapping 87Rb and 129Xe spins, their Fermi-contact
interaction introduces an effective field $\lambda M^{n}\textbf{P}^{n}$ on
87Rb spins generated by 129Xe nuclear magnetization [5]. Such an effective
field can be described by

$\textbf{B}_{\rm{eff}}=\lambda
M^{n}\textbf{P}^{n}=\frac{8\pi}{3}\kappa_{0}M^{n}\textbf{P}^{n},$ (S4)

where the Fermi-contact enhancement factor $\kappa_{0}$ is about 540 and
$\mu_{0}$ is the permeability of vacuum. Similarly, 87Rb magnetization also
generates an effective field $\lambda M^{e}\textbf{P}^{e}$ on 129Xe spins. In
our experiment, $\lambda M^{e}\textbf{P}^{e}$ is on the order of nT and
$\lambda M^{n}\textbf{P}^{n}$ is at least one order greater than $\lambda
M^{e}\textbf{P}^{e}$. In Sec. VI, we present some numerical calculations of
$\textbf{B}_{\rm{eff}}$.

We now consider how to solve the coupled Bloch equations (Eqs. S2 and S3). At
first, we simplify them. Unlike the near-zero field case where noble gas spins
and electron spins are strongly coupled together and operate as a
comagnetometer [4], we consider a different situation where 87Rb and 129Xe
spins are under a large bias field $B_{0}^{z}$ along $z$ (above 200 nT). In
this situation, 87Rb and 129Xe spins are weakly coupled together, greatly
simplifying Eqs. S2 and S3. The details are explained as follows. (1)
$P^{n}_{z}$ can be approximated as a constant, which is much larger than
$P^{n}_{x}$ and $P^{n}_{y}$. Based on Eq. S3, the 129Xe spins evolve under the
static field of the effective field $\lambda M^{e}{P}_{z}^{e}$ along $z$
generated by 87Rb spins and the bias field $B^{z}_{0}$. Because the applied
$B_{0}^{z}$ is much larger than the effective field $\lambda
M^{e}{P}_{z}^{e}$, we can first neglect the term $\lambda M^{e}{P}_{z}^{e}$ in
Eq. S3 for simplicity (but the 87Rb effective field along $z$ should be
considered when the Larmor frequency of 129Xe is calibrated, see Sec. IV.1).
As a result, Eq. S3 can be independently solved at first and $\textbf{P}^{n}$
can be evaluated. (2) Then, we can take the solution of $\textbf{P}^{n}$ into
Eq. S2 and solve for $\textbf{P}^{e}$.

We write the total magnetic field experienced by the 87Rb spins as
$\textbf{B}=\textbf{B}^{0}+\lambda{M}^{n}\textbf{P}^{n}$ and obtain the
simplified Bloch equations,

$\displaystyle\frac{\partial\textbf{P}^{e}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{e}}{Q}\textbf{B}\times\textbf{P}^{e}+\frac{P^{e}_{0}\bm{z}-\textbf{P}^{e}}{T_{e}Q},$
(S5) $\displaystyle\frac{\partial\textbf{P}^{n}}{\partial t}$ $\displaystyle=$
$\displaystyle\gamma_{n}\textbf{B}^{0}\times\textbf{P}^{n}+\frac{P^{n}_{0}\bm{z}-\textbf{P}^{n}}{\\{T_{2n},T_{2n},T_{1n}\\}}.$
(S6)

We first solve the evolution of 129Xe spins under a bias field $B^{0}_{z}$
along $z$ and an oscillating magnetic field $B_{\textrm{ac}}^{y}\cos(2\pi\nu
t)$. The Larmor frequency of 129Xe spins can be written as
$\nu_{0}=\gamma_{n}B_{0}^{z}/(2\pi)$. With rotating-wave approximation, we can
rewrite the Bloch equation of 129Xe spins in the rotating frame,

$\displaystyle\dfrac{\partial}{\partial
t}\tilde{\textbf{P}}^{n}=\gamma_{n}\tilde{\textbf{B}}^{0}\times\tilde{\textbf{P}}^{n}-\dfrac{\tilde{P}_{x}^{n}\bm{x}+\tilde{P}_{y}^{n}\bm{y}}{T_{2n}}-\dfrac{(P_{z}^{n}-P_{0}^{n})\bm{z}}{T_{1n}},$
(S7)

where
$\tilde{\textbf{B}}^{0}=\dfrac{2\pi(\nu_{0}-\nu)}{\gamma_{n}}\bm{z}+\dfrac{B_{\textrm{ac}}^{y}}{2}\bm{y}$
is the equivalent magnetic field in the rotating frame. We derive the three
components of $\tilde{\textbf{P}}^{n}$ by solving Eq. S7 in rotating frame and
transforming back to the laboratory frame. Finally, the transverse
polarization components are

$\displaystyle
P^{n}_{x}=\dfrac{1}{2}P^{n}_{0}\gamma_{n}B_{\textrm{ac}}^{y}\dfrac{T_{2n}\cos(2\pi\nu
t)+2\pi(\nu-\nu_{0})T_{2n}^{2}\sin(2\pi\nu
t)}{1+({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}+[2\pi(\nu-\nu_{0})]^{2}T_{2n}^{2}}+Ce^{-t/T_{2n}}\cos(2\pi\nu_{0}t),$
(S8) $\displaystyle
P^{n}_{y}=\dfrac{1}{2}P^{n}_{0}\gamma_{n}B_{\textrm{ac}}^{y}\dfrac{T_{2n}\sin(2\pi\nu
t)-2\pi(\nu-\nu_{0})T_{2n}^{2}\cos(2\pi\nu
t)}{1+({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}+[2\pi(\nu-\nu_{0})]^{2}T_{2n}^{2}}+Ce^{-t/T_{2n}}\sin(2\pi\nu_{0}t),$

where $C$ is a constant determined by initial conditions. When the oscillating
field $B^{y}_{\rm{ac}}$ is suddenly applied to the vapor cell, the transient
response of 129Xe spins can be created, as described in Eq. S8. According to
$\textbf{B}_{\rm{eff}}=\lambda M^{n}\textbf{P}^{n}$, we can derive the
effective field experienced by 87Rb atoms (or 87Rb magnetometer) as

$\displaystyle\textbf{B}_{\rm{eff}}=\textbf{B}_{\rm{eff}}^{\rm{CW}}+\textbf{B}_{\rm{eff}}^{\rm{Tran}},$
(S9)
$\displaystyle\textbf{B}_{\rm{\rm{eff}}}^{\rm{CW}}=\overbrace{\dfrac{1}{2}\lambda
M^{n}P^{n}_{0}\gamma_{n}B_{\textrm{ac}}^{y}\dfrac{T_{2n}\cos(2\pi\nu
t)+2\pi(\nu-\nu_{0})T_{2n}^{2}\sin(2\pi\nu
t)}{1+({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}+[2\pi(\nu-\nu_{0})]^{2}T_{2n}^{2}}\bm{x}}^{\textrm{
effective field generated by ${}^{129}$Xe $x$ magnetization}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\overbrace{\dfrac{1}{2}\lambda
M^{n}P^{n}_{0}\gamma_{n}B_{\textrm{ac}}^{y}\dfrac{T_{2n}\sin(2\pi\nu
t)-2\pi(\nu-\nu_{0})T_{2n}^{2}\cos(2\pi\nu
t)}{1+({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}+[2\pi(\nu-\nu_{0})]^{2}T_{2n}^{2}}\bm{y}}^{\textrm{effective
field generated by ${}^{129}$Xe $y$ magnetization}},$
$\displaystyle\textbf{B}_{\rm{eff}}^{\rm{Tran}}=\overbrace{\lambda
M^{n}Ce^{-t/T_{2n}}\cos(2\pi\nu_{0}t)\bm{x}+\lambda
M^{n}Ce^{-t/T_{2n}}\sin(2\pi\nu_{0}t)\bm{y}}^{\textrm{transient response}}.$

To measure the steady-state response signal (i.e.,
$\textbf{B}_{\rm{\rm{eff}}}^{\rm{CW}}$) from $\textbf{B}_{\rm{eff}}$, we set
50 s waiting time to ensure that $\textbf{B}_{\rm{eff}}^{\textrm{Tran}}$
decays to zero and thus the transient term can be neglected before data
sampling. Hence, in the following, $\textbf{B}_{\rm{eff}}^{\rm{Tran}}$ is
neglected and $\textbf{B}_{\rm{eff}}\approx\textbf{B}_{\rm{eff}}^{\rm{CW}}$.
The total oscillating magnetic field $\textbf{B}_{\rm{tot}}$ on 87Rb detected
can be written as

$\textbf{B}_{\rm{tot}}=\textbf{B}_{\rm{eff}}^{\textrm{CW}}+B^{y}_{\rm{ac}}\cos(2\pi\nu
t)\bm{y}.$ (S10)

Figure S2: Response to an oscillating field ($B_{y}$) of the 87Rb magnetometer
as a function of the bias magnetic field along with a fit to
$\dfrac{-B_{y}{\Delta}B}{(B_{z}^{0})^{2}+(\Delta B)^{2}}$, where the parameter
is $\Delta B\approx 2700\pm 13$ nT. In experiment, a 7.11 nT field with 320 Hz
frequency is applied along ${y}$. The bias field $B^{0}_{z}$ is scanned from 0
to 8000 $\rm{nT}$, and the amplitude of the oscillating signal is recorded.

We now solve for the evolution of 87Rb spins (Eq. S5). In our experiment, the
${x}$ component of 87Rb polarization $P^{e}_{x}$ is detected with a probe beam
(see Eq. S1). Thus, we only need obtain the explicit expression of
$P^{e}_{x}$. Under the quasi-static field condition, we can obtain the steady-
state solution

$\displaystyle
P_{x}^{e}\propto\dfrac{B_{x}B_{z}-B_{y}{\Delta}B}{|\textbf{B}|^{2}+(\Delta
B)^{2}}\approx\dfrac{B_{x}B_{z}-B_{y}{\Delta}B}{(B_{z}^{0})^{2}+(\Delta
B)^{2}},$ (S11)

where ${\bigtriangleup}B=1/(\gamma_{e}T_{e})$. To calibrate the parameter
${\bigtriangleup}B$, we apply a 7.11 nT field with 320 Hz frequency along
${y}$ and scan the bias field from 0 to 8000 $\rm{nT}$. Figure S2 provides the
experimental data, which is fitted with Eq. S11 and yields $\Delta B\approx
2700\pm 13$ nT. In our experiment, a static bias field $B^{z}_{0}$ above 200
nT is applied along the ${z}$ axis. In this situation, the 87Rb magnetometer
becomes simultaneously sensitive to the magnetic fields along ${x}$ and ${y}$.
This is in contrast to near-zero-field magnetometers, where 87Rb magnetometer
is only sensitive to the magnetic field along $y$. Therefore, 87Rb
magnetometer can measure the effective magnetic fields along ${x}$ and ${y}$
produced by 129Xe nuclear $x,y$-magnetization.

###  II.1 Linear and nonlinear response

The effective field $\textbf{B}_{\textrm{eff}}$ induced by 129Xe magnetization
has a close relation with the measured $B_{\textrm{ac}}^{y}$ (for example,
axion-like dark matter field). Thus, it is important to analyse the response
of detectable $\textbf{B}_{\textrm{eff}}$ to unknown $B_{\textrm{ac}}^{y}$. We
consider three situations.

  * (1)

When the applied oscillating field $B^{y}_{\textrm{ac}}$ is weak and satisfies
$({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}\ll 1$, based on Eq. S9,
the effective field is proportional to the strength of $B^{y}_{\rm{ac}}$. In
this situation, the 129Xe spin-based amplifier has linear response to the
measured $B_{\textrm{ac}}^{y}$. According to our experimental parameters, the
$B_{\textrm{ac}}^{y}$ should satisfy

$B^{y}_{\textrm{ac}}\ll\frac{2}{\gamma_{n}\sqrt{T_{1n}T_{2n}}}\approx
8.5~{}\textrm{nT}.$ (S12)

  * (2)

When the applied oscillating field $B^{y}_{\rm{ac}}$ is comparable to that of
Eq. S12, the $({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}$ term
becomes significant. In this situation, the response is nonlinear.

  * (3)

When the applied oscillating field $B^{y}_{\textrm{ac}}$ is much larger, the
term $({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}$ cannot be neglected
and causes nuclear magnetization-saturation effects. The effective field
$\textbf{B}_{\rm{eff}}$ is negligible. In this situation, the linear response
signal of the 87Rb magnetometer is dominant.

We verify the above discussions (1)-(3) with experiments. Figure S3a provides
experimental linear and nonlinear response signals, scanned as a function of
the oscillating field amplitude. In contrast to the far off-resonant case (320
Hz), the signals for resonant ($8.96$ Hz) and near-resonant cases ($8.80$ Hz,
$9.05$ Hz) first linearly increase, then decrease, and lastly linearly
increase again, originating from nuclear magnetization saturation. For
example, the resonant signal (red circles) has a large response slope when the
oscillating field amplitude is weak. Due to the large slope, the signal from
an external oscillating field can be enhanced. With higher oscillating field
amplitude, the response signal decreases. When the oscillating field amplitude
is increased further, the response signal from the 87Rb magnetometer becomes
linear again. In the large oscillating field amplitude, the slopes for
resonant ($8.96$ Hz) and near-resonant ($8.80$ Hz, $9.05$ Hz) cases are the
same as that of the far off-resonant case (320 Hz). This is because the
response signals are from the rubidium magnetometer. The solid curves in Fig.
S3b represent our theoretical calculations based on Eq. S9, which agree well
with the experiment.

Figure S3: Response of the spin-based amplifier as a function of the
oscillating field amplitude. (a) Experimental response as a function of
oscillating magnetic field amplitude. Each experiment (with the oscillation
frequency of 7.00, 8.80, 8.96, 9.05, 320 Hz, respectively) is performed in the
same bias field $B_{z}^{0}\approx 759$ nT (corresponding to that 129Xe Larmor
frequency is 8.96 Hz). (b) Theoretical profiles corresponding to the
experiment shown in (a). The experimental data agrees well with the theoretcal
profiles.

A particularly sensitive window for measuring external oscillating fields
corresponds to the amplitude of oscillating fields below several nT (see Eq.
S12). In practice, we are interested in a sensitive measurement of small
magnetic fields (for example, axion-like dark matter field), whose amplitudes
are within the sensitive window demonstrated in this work. To maintain the
operation in the sensitive window, we should suppress ambient electromagnetic
interference, for example, by using a five-layer magnetic shield and an
ultralow-noise current source.

Apart from the oscillating field amplitude, the effective magnetic field
$\textbf{B}_{\textrm{eff}}$ also depends on the oscillating field frequency.
We discuss the frequency dependence, as shown in Fig. S4: (1) On-resonance
case: when $\nu\approx\nu_{0}$, the effective field $\textbf{B}_{\rm{eff}}$
reaches a maximum and could be much larger than the oscillating field
amplitude $B^{y}_{\rm{ac}}$. We present how to calculate and experimentally
measure the amplification factor, as discussed in Sec. II.2; (2) Near-
resonance case: the effective field increases when the frequency of the
oscillating field frequency $\nu$ becomes close to the Larmor frequency
$\nu_{0}$. Thus, there is a frequency bandwidth for the spin-based amplifier,
as discussed in Sec. II.3. (3) Far-off-resonance case: when $\nu\gg\nu_{0}$ or
$\nu\ll\nu_{0}$, the term $[2\pi(\nu-\nu_{0})]^{2}T_{2n}^{2}$ is dominant in
Eq. S9. In this situation, the effective field $\textbf{B}_{\rm{eff}}$
generated by 129Xe spins is negligible and the applied oscillating field
$B^{y}_{\rm{ac}}\cos(2\pi\nu t)\bm{y}$ is dominant.

Figure S4: Amplification performance of the spin-based amplifier. (a) On- and
near-resonance cases. When the frequency of the applied oscillating field
$\nu$ is close to the 129Xe Larmor frequency $\nu_{0}$, the effective field
$\textbf{B}_{\rm{eff}}$ is considerable. Such an effective field can be
measured with 87Rb magnetometer. (b) Far-off-resonance cases. When the
frequency of applied oscillating field $\nu$ is away from resonance frequency
$\nu_{0}$, the effective field $\textbf{B}_{\rm{eff}}$ generated from 129Xe
spins is negligible.

###  II.2 Amplification factor

The resonant oscillating field can induce an oscillating 129Xe nuclear
magnetization, which can generate a considerable effective magnetic field
$\textbf{B}_{\rm{eff}}$ on 87Rb spins. As discussed below, the strength of
$\textbf{B}_{\rm{eff}}$ could be much larger than the applied oscillating
field $B^{y}_{\rm{ac}}$. Accordingly, we define an amplification factor:

$\eta=|\textbf{B}_{\rm{eff}}/B^{y}_{\rm{ac}}|.$ (S13)

We first derive the amplification factor $\eta$ on resonance. Here we assume
that the oscillating field strength is small and thus the term
$({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}$ can be neglected in Eq.
S9. In this situation, the spin-based amplifier works in the sensitive linear-
response regime and $\textbf{B}_{\rm{eff}}$ can be written as

$\textbf{B}_{\textrm{eff}}(\nu=\nu_{0})=\dfrac{1}{2}\lambda
M^{n}P^{n}_{0}\gamma_{n}T_{2n}[\cos(2\pi\nu t)\bm{x}+\sin(2\pi\nu
t)\bm{y}]B_{\textrm{ac}}^{y}.$ (S14)

Based on this result, the effective field $\textbf{B}_{\textrm{eff}}$ is a
circularly polarized field (see Fig. S5) and its amplitude is equal to
$\dfrac{1}{2}\lambda M^{n}P^{n}_{0}\gamma_{n}T_{2n}\cdot B_{\textrm{ac}}^{y}$.
As a result, the amplification factor is

$\eta=\dfrac{1}{2}\lambda M^{n}P^{n}_{0}\gamma_{n}T_{2n}.$ (S15)

Based on Eq. S15, many ways can be used to increase the amplification factor,
such as prolonging relaxation time $T_{2n}$ and improving equilibrium
polarization $P^{n}_{0}$. The amplification factor can be as large as 104 in a
3He-K magnetometer (see Sec. VI for details). The sensitivity of 3He-K
magnetometer can be improved by four orders of magnitude and can potentially
reach a few aT/$\sqrt{\rm{Hz}}$.

Figure S5: Effective magnetic field $\textbf{B}_{\textrm{eff}}$. (a) An
oscillating field $B_{\textrm{ac}}^{y}\bm{y}$ is applied along the $y$ axis.
(b) The magnetic field generated from 129Xe magnetization traces a circle in
the $xy$ plane.

Although the amplification factor can be theoretically calculated based on Eq.
S15, this requires knowledge of experimental parameters
$\\{M^{n},P_{0}^{n},T_{2n}\\}$. Alternatively, we provide a simple
experimental method to measure the amplification factor $\eta$. The method
consists of the following steps:

  * (1)

A resonant oscillating field is applied along $y$; the output signal of 87Rb
magnetometer is an oscillating voltage signal and we record its amplitude as
$A(B_{\textbf{eff}})$. By taking $\textbf{B}_{\textrm{eff}}$ (see Eq. S14)
into $P_{x}^{e}$ (Eq. S11), we achieve the output signal of the 87Rb
magnetometer

$\displaystyle A(B_{\textbf{eff}})\propto
P^{e}_{x}(\textbf{B}_{\textrm{eff}})\propto\dfrac{1}{2}\lambda
M^{n}P_{0}^{n}\gamma_{n}T_{2n}\dfrac{1}{\sqrt{[{(B^{0}_{z})^{2}}+({\Delta
B})^{2}]}}B_{\rm{ac}}^{y}.$ (S16)

In a practical measurement, we usually do not know the exact Larmor frequency
without prior calibration. In this situation, an alternative way is to scan
the oscillating field frequency over a small frequency range, corresponding to
the narrow bandwidth of the spin-based amplifier (see Sec. II.3). The
oscillating field frequency is automatically changed through an AWG (keysight
33210A) with a 0.04 Hz step applied. We fit the experimental data and find the
maximum as the resonant signal amplitude.

  * (2)

A far-off-resonant oscillating field is applied along $y$; similarly, we
record the amplitude $A(B_{\textbf{ac}}^{y}\textbf{y})$ of the 87Rb
magnetometer output signal. In this situation, we only need to consider
$B^{y}_{\textrm{ac}}\bm{y}$ without considering $\textbf{B}_{\textrm{eff}}$.
The direct response of the 87Rb magnetometer to $B^{y}_{\textrm{ac}}\bm{y}$
can be expressed as

$\displaystyle A(B_{\textbf{ac}}^{y}\bm{y})\propto
P^{e}_{x}(B^{y}_{\textrm{ac}}\bm{y})\propto\dfrac{\Delta
B}{[{(B^{0}_{z})^{2}}+({\Delta B})^{2}]}B_{\rm{ac}}^{y}.$ (S17)

Here, the frequency of the far-off-resonant field should be within the
bandwidth of the 87Rb magnetometer, otherwise a frequency-response factor
should be taken into account.

  * (3)

By comparing the above two amplitudes, we have

$\displaystyle\Phi(\nu_{0})$
$\displaystyle=\dfrac{A(B_{\textbf{eff}})}{A(B_{\textbf{ac}}^{y}\bm{y})}=\dfrac{1}{2}\lambda
M^{n}P_{0}^{n}\gamma_{n}T_{2n}\cdot\sqrt{1+\left(\dfrac{\nu_{0}}{\gamma_{n}\Delta
B}\right)^{2}},$ (S18)

and further, the amplification factor can be evaluated by the following
relation

$\displaystyle\eta$
$\displaystyle=\Phi(\nu_{0})/\sqrt{1+\left(\dfrac{\nu_{0}}{\gamma_{n}\Delta
B}\right)^{2}}.$ (S19)

Figure S6a shows the experimental $\Phi(\nu_{0})$ as a function of the
resonant frequency, which agrees well with our theoretical model (shown with
green line) described in Eq. S18. The amplitude $A(B_{\textbf{eff}})$ scales
as $1/B_{z}^{0}$, while the amplitude of Eq. S17 scales as
$1/(B_{z}^{0})^{2}$. This is eventually leads to the dependence on Fig. S6a.
The amplification factors calculated according to Eq. S19 are shown in Fig.
S6b. The experimental result shows that the achieved amplification factors are
nearly independent of $\nu_{0}$, which agrees with our theoretical model. The
average of such amplification factors is $128\pm 0.3$.

Figure S6: Measurements of amplification factor. (a) $\Phi$ as a function of
resonant frequency (corresponding to different applied bias field). The solid
line is a fit based on Eq. S18 with $\eta\approx 128.1$. (b) The measured
amplification factor $\eta$ at different resonant frequencies. The averaged
$\eta$ is measured to be $\eta=128\pm 0.3$.

In the axion-like dark matter searches (Secs. III and IV), we focus on the
search via the axion-nucleon interactions known as ‘the axion gradient effect’
[6, 7, 8, 9, 10, 11, 12] (see Sec. III.1). In this model, axion-like dark
matter acts as a time-oscillating magnetic field which couples to nuclear
spins. Through spin-based amplifier, the effect of axion-like field can be
enhanced and the measurement sensitivity can be greatly improved.

###  II.3 Bandwidth

In Sec. II.2, we only consider the amplification factor in the resonant case.
We now analyse the amplification performance in the near-resonant case. Based
on the analysis of the frequency dependence above, the amplification factor
reaches a maximum on resonance, and rapidly decreases when the frequency of
applied oscillating magnetic field $\nu$ is far off-resonance. Therefore,
there is a bandwidth for spin-based amplifier. Based on Eq. S9, we can rewrite
the equation in the amplitude spectrum,

$\displaystyle|\textbf{B}_{\textrm{eff}}^{\textrm{CW}}|\propto\dfrac{\Lambda/2}{\sqrt{(\nu-\nu_{0})^{2}+(\Lambda/2)^{2}}}.$
(S20)

We assume that the oscillating field amplitude is weak and thus the term
$({\gamma_{n}B_{\textrm{ac}}^{y}}/2)^{2}T_{1n}T_{2n}$ can be neglected in Eq.
S9. As a demonstration, we set the bias field as $B_{0}^{z}=759$ nT and the
oscillating field strength $B_{\textrm{ac}}^{y}=30$ pT. The frequency response
is measured from 8 to 10 Hz, as shown in Fig. S7. The experimental data are
well described as a single-pole, band-pass filter (solid red curve). The full-
width at half-maximum (FWHM) is $\sqrt{3}\Lambda\approx 0.052$ Hz; the spin-
based amplifier can enhance the signal in a correspondingly narrow frequency
range.

Figure S7: Bandwidth measurement. With a fit to
$A/\sqrt{(\nu-\nu_{0})^{2}+(\Lambda/2)^{2}}$, where the bandwidth
$\sqrt{3}\Lambda=$ 0.052 Hz.

##  III Derivation of the effective magnetic field induced by axion-like dark
matter

This section derives the effective magnetic field induced by axion-like dark
matter.

###  III.1 Axion dark-matter halo

Axion-like particles could be hypothetically produced in the early Universe by
non-thermal mechanisms, such as “vacuum misalignment” [13]. Then, they form a
coherently oscillating classical field $a(t)=a_{0}\cos(2\pi\nu_{a}t)$ [6]. The
axion-like dark matter field $a(t)$ oscillates at its Compton frequency

$\displaystyle\nu_{a}=\dfrac{m_{a}c^{2}}{h},$ (S21)

where $m_{a}$ is the mass of the axion-like particle, $c$ is the speed of
light, and $h$ is the Planck constant. The coherence time determined by the
random fluctuations of the field corresponds to the de Broglie wavelength
$\lambda_{\rm{dB}}$,

$\displaystyle\tau_{a}\approx\dfrac{\lambda_{\rm{dB}}}{v}\approx\dfrac{h}{m_{a}v^{2}},$
(S22)

where $|\bm{v}|\sim 10^{-3}c$ is the galactic virial velocity of the axion
dark-matter field. Assuming that axion field comprises the totality of the
local dark matter energy density $\rho_{\rm{DM}}\approx 0.4~{}\rm GeV/cm^{3}$,
the amplitude $a_{0}$ of the field obeys the equation,

$\displaystyle\rho_{\rm{DM}}\approx\dfrac{c^{2}}{2\hbar^{2}}m_{a}^{2}a_{0}^{2}.$
(S23)

When axion-like dark matter field interacts with nuclear spins, the axion-like
dark matter can be considered as a time-oscillating pseudo-magnetic field
directed along the $\bm{v}$. Thus, the dark-matter couplings to nuclear spins
can be seen as a Zeeman interaction between nuclear spins and this pseudo-
magnetic field generated by axion-like particles. The resultant Hamiltonian is
[6, 7, 8, 9, 10, 11, 12]

$\displaystyle\mathcal{H}_{\rm{spin}}\approx\textrm{g}_{\rm{aNN}}\sqrt{2\hbar^{3}c\rho_{\rm{DM}}}\sin(2\pi\nu_{a}t)\bm{v}\cdot\textbf{I}_{\rm{N}},$
(S24)

where $\textrm{g}_{\rm{aNN}}$ is the strength of the axion-nucleon coupling
and $\textbf{I}_{\rm{N}}$ is the spin-matrix vector of the nuclear spin. For a
nucleus with a particular gyromagnetic ratio
$\gamma_{n}=g_{n}\mu_{\rm{N}}/\hbar$, the $g_{n}$ is the nuclear
Land$\acute{e}$ factor and $\mu_{\rm{N}}$ is the nuclear magneton. As a
result, we can derive the effective magnetic field induced by axion-like dark
matter,

$\displaystyle\textbf{B}_{a}\approx
10^{-3}\dfrac{\textrm{g}_{\rm{aNN}}}{\hbar\gamma_{n}}\sqrt{2\hbar^{3}c^{3}\rho_{\rm{DM}}}\sin(2\pi\nu_{a}t)\times\dfrac{\bm{v}}{|\bm{v}|},$
(S25) $\displaystyle\textbf{B}_{a}[\textrm{T}]\approx
10^{-7}\dfrac{\rm{g}_{\rm{aNN}}[\rm{GeV}^{-1}]}{\rm{g}_{n}}\sin(2\pi\nu_{a}t)\times\dfrac{\bm{v}}{|\bm{v}|}.$

According to Eq. S25, we can build the relation between the detection of
magnetic field and the coupling strength $\textrm{g}_{\rm{aNN}}$. Here,
$\textrm{g}_{\rm{aNN}}$ is the important parameter that we should measure or
constrain in experiments. In Sec. IV, we present the details how to use our
experimental data to constrain the coupling strength $\textrm{g}_{\rm{aNN}}$
over an axion mass range, from 10 feV to 1 peV or 2 Hz to 200 Hz.

###  III.2 Coordinate transformation

Based on Eq. S25, the effective magnetic field is along the galactic virial
velocity $\bm{v}$. For terrestrial experiments, we need to calculate the
projection component of $\textbf{B}_{a}$ along the sensitive axis of our setup
shown in Fig. S8b. To this end, we should calculate the projection of the
virial velocity along sensitive axis. In lab coordinates, 129Xe spins are
polarized along the ${z}$ axis and the transverse magnetization is detected.
As discussed in Sec. II, the sensitive direction of the spin-amplifier-based
magnetometer is perpendicular to the bias field (along ${z}$), see Fig. S8c.
Thus, we need to calculate the projection of pseudo-magnetic field in the
plane perpendicular to the ${z}$ axis in laboratory frame.

There are three coordinate systems, which are galactic coordinates, the
celestial coordinates and laboratory coordinates. The virial velocity
direction is defined as a direction of 270∘ longitude and 0∘ latitude in
galactic coordinates. We use the coordinate convert tool from the official
NASA website (NASA LAMBDA–Tools), to convert the direction of the axion
velocity from galactic coordinates to celestial coordinates. The axion
velocity in celestial coordinates can be written as
$\bm{v}=|\bm{v}|(\cos\delta\cos\eta\bm{X}+\cos\delta\sin\eta\bm{Y}+\sin\delta\bm{Z})$,
where declination is $\delta=-48^{\circ}$ and right ascension is
$\eta=138^{\circ}$ [8, 9].

Next, we transform the axion velocity in celestial coordinates into laboratory
coordinates. Based on ref. [14], we can use the following matrix to convert
coordinates,

$\begin{bmatrix}x\\\ y\\\ z\end{bmatrix}=\begin{bmatrix}\cos\chi\cos\Omega
t&\cos\chi\sin\Omega t&-\sin\chi\\\ -\sin\Omega t&\cos\Omega t&0\\\
\sin\chi\cos\Omega t&\sin\chi\sin\Omega
t&\cos\chi\end{bmatrix}\begin{bmatrix}X\\\ Y\\\ Z\end{bmatrix},$ (S26)

where $\begin{bmatrix}x&y&z\end{bmatrix}^{\rm{T}}$ and
$\begin{bmatrix}X&Y&Z\end{bmatrix}^{\rm{T}}$ are lab and celestial
coordinates, $\Omega=2\pi\times 1.16\times 10^{-5}$ s-1 is the sidereal
frequency or the Earth’s rotation frequency, $\chi$ is the angle between the
${z}$ axis in lab and in celestial coordinates. We can derive this angle
through the equation $\cos\chi=\cos\alpha\cos\beta$, where $\alpha$ is the
latitude of our lab and $\beta$ is the angle between the ${z}$ axis of lab and
the Geographic North (the axis ${z}$ of celestial coordinates projected in the
lab plane). Our apparatus for axion-like dark matter searches is located in
University of Science and Technology of China (Hefei, China), where
$\alpha=30^{\circ}$ and $\beta=90^{\circ}$. From the matrix Eq. S26, we can
obtain the relation between coordinates transformation $z=\sin\chi\cos\Omega
tX+\sin\chi\sin\Omega tY+\cos\chi Z$. Finally, we derive the axion velocity
$v_{z}$ in lab coordinates,

$v_{z}=\bm{v}\cdot\bm{z}=|\bm{v}|[\cos\chi\sin\delta+\sin\chi\cos\delta\cos(\Omega
t-\eta)].$ (S27)

Because $\cos\chi$ is close to 0 based on the relation
$\cos\chi=\cos\alpha\cos\beta$, we can ignore the $\cos\chi\sin\delta$ term.
Based on this simplification, we derive the transverse velocity:

$\sqrt{v_{x}^{2}+v_{y}^{2}}=\sqrt{\bm{v}^{2}-v_{z}^{2}}=|\bm{v}|\sqrt{1-[\sin\chi\cos\delta\cos^{2}(\Omega
t-\eta)]}=|\bm{v}|\sqrt{1-3/4\cos^{2}(\Omega t-\eta)}.$ (S28)

Substituting the velocity term back into Eq. S24, we derive the equation
utilized to search for axion-like dark matter in our experiment:

${B}_{a}[\rm{T}]\approx
10^{-7}\dfrac{\textrm{g}_{\rm{aNN}}[\rm{GeV}^{-1}]}{g_{n}}\sin(2\pi\nu_{a}t)\times\sqrt{1-3/4\cos^{2}(\Omega
t-\eta)}.$ (S29)

According to Eq. S29, axion-like dark matter induced pseudo-magnetic field is
modulated by the Earth’s rotation. Because the sampling time for each search
measurement is much shorter than the Earth’s rotation period,
$\sqrt{1-3/4\cos^{2}(\Omega t-\eta)}$ during the sampling time can be
appromixated as the value at the initial time of each search experiment.

Figure S8: Coordinates transformation. (a) Relative relation between the
celestial coordinates and the lab plane in celestial coordinates. (b) Relation
between the celestial North and the lab coordinates. (c) The magnetic
sensitive plane of the spin-based amplifier.

##  IV Data analysis for the axion-like dark matter search

This section presents the detailed experimental procedure and data analysis of
axion-like dark matter searches. Our search experiment includes two parts,
i.e., calibration and axion search.

###  IV.1 Calibration experiments

As demonstrated in Sec. II.3, the 87Rb magnetometer assisted with 129Xe spin-
based amplifier can be highly sensitive to the oscillating magnetic field
whose frequency is within the bandwidth of the spin-based amplifier. To use
such a magnetometer in the searches of axion-like dark matter, it is the
prerequisite to calibrate the frequency bandwidth of magnetometer and the
magnetic sensitivity over the bandwidth. To this end, we should perform three
calibration experiments:

  * (1)

Calibration of the relation between 129Xe Larmor frequency and applied coil
current. The bias magnetic field is generated with a commercial current source
(Krohn-Hite Model 523). This current source drives a set of coils (along $z$).
Due to the existence of 87Rb effective field (see Sec. II) and magnetic shield
residual field, we need to carefully calibrate the relation between 129Xe
Larmor frequency and the applied current. At first, to measure the 129Xe
Larmor frequency, we apply a $\pi/2$ magnetic-field pulse along the ${x}$ axis
to rotate the 129Xe spins from ${z}$ axis to ${y}$ axis, and then the free-
decay of 129Xe signals are measured for 100 s. The free-decay signal can be
fitted with a single-exponential-decay trigonometric function, which gives
corresponding Larmor frequency. Figure S9a shows experimental Larmor
frequencies (black circles) as a function of the applied current. The linear
fit for the experimental data yields

$\nu_{0}(I)=\frac{\gamma_{n}}{2\pi}(761\cdot I+0.66),$ (S30)

where $I$ is the applied current (mA) and the conversion coefficient between
magnetic field and applied current is about $\approx 761$ nT/mA. The intercept
of 0.66 nT represents the sum of 87Rb effective field and magnetic shield
residual field. Based on our calibration, we determine the dark-matter search
frequency range when the applied current is set.

Figure S9: Results of experimental calibrations. (a) Calibration of relation
between 129Xe Larmor frequency and applied coil current. (b) Calibration of
87Rb magnetometer magnetic sensitivity assisted the 129Xe spin-based
amplifier. Note that $18~{}\textrm{fT}/\textrm{Hz}^{1/2}$ is achieved at 129Xe
Larmor frequency ($\nu_{0}\approx 8.96$ Hz), which is beyond the sensitivity
($\approx$2 pT/Hz1/2) of the rubidium magnetometer itself. (c) Calibration of
the enhanced magnetic sensitivity over different bias fields. There are some
spurious noise peaks (e.g., around 80 Hz and 180 Hz) that originate from the
probe-beam intensity noise, which can be removed by properly tuning and
stabilizing the probe laser.

  * (2)

Calibration of 87Rb magnetometer magnetic sensitivity over different bias
fields. Here the amplification effect is not taken into account and we only
calibrate the 87Rb magnetometer sensitivity itself. A certain bias field is
first set and then we wait for 50 s before 100 s sampling. The reason of 50-s
waiting is to make sure that the set current becomes stable. During the
procedure, a 320 Hz calibration field with 7.11 nT amplitude is applied along
$y$. Based on the information of the calibration field, we determine the
conversion coefficient between voltage (the 87Rb magnetometer output signal is
voltage) and applied oscillating field. Then, we calculate the power spectrum
of 100 s data, which can be converted to magnetic field sensitivity. In
particular, for each bias field $B^{0}_{z}$, we record the 87Rb magnetometer
sensitivity at the corresponding 129Xe Larmor frequency
$\nu_{0}=\gamma_{n}B^{0}_{z}$, which is the most sensitive frequency because
of spin-based amplifier (see Sec. II.3). Further, we scan bias fields to
repeat the above step. In a small bias field, for example at $B^{0}_{z}=1125$
nT, the magnetic sensitivity can reach
$\Upsilon_{\rm{Rb}}(\gamma_{n}B^{0}_{z}\approx 13.3)\approx
1.5~{}\rm{pT}/\sqrt{\rm{Hz}}$. With increasing bias field $B^{0}_{z}$, the
sensitivity gradually decreases due to the loss of $\frac{\partial
P_{x}^{e}}{\partial B_{y}}\propto[(B^{0}_{z})^{2}+(\Delta B)^{2}]^{-1}$ (see
Eq. S11). In a large bias field, for example at $B^{0}_{z}=15000$ nT, the
sensitivity is $\Upsilon_{\rm{Rb}}(\gamma_{n}B^{0}_{z}\approx 177.0)\approx
58~{}\rm{pT}/\sqrt{\rm{Hz}}$.

  * (3)

Calibration of amplification factor and enhanced magnetic sensitivity. The
detailed procedure to measure amplification factor is described in Sec. II.2.
The amplification factor is measured to be $\eta\approx 128$. Further, we can
calculate the enhanced magnetic sensitivity at scanned Larmor frequency

$\Upsilon_{\rm{amp}}(\nu_{0},\nu_{0})=\frac{\Upsilon_{\rm{Rb}}(\nu_{0})}{\Phi(\nu_{0})},$
(S31)

where $\Phi(\nu_{0})$ represents the ratio between the resonant response
signal amplitude and far off-resonant response signal amplitude, as seen in
Eq. S19:

$\Phi(\nu_{0})=\eta\sqrt{1+\left(\frac{\nu_{0}}{\gamma_{n}\Delta
B}\right)^{2}}.$

Figure S6a shows the experimental result for $\Phi(\nu_{0})$. Figure S9b
signals the significant improvement of magnetic sensitivity of the spin-based
amplifier, comparing with the 87Rb magnetometer itself. Figure S9c shows the
87Rb magnetometer magnetic sensitivity at 129Xe Larmor frequency over
different bias fields. However, the scanned Larmor frequency is discrete. For
other frequencies $\nu$ nearby each resonant frequency $\nu_{0}$, as
demonstrated in Sec. II.3, the spin-based amplifier can operate within the
bandwidth of 0.052 Hz centered at 129Xe Larmor frequency. Therefore, it is
necessary to derive the explicit form of the magnetic sensitivity of the spin-
based amplifier within the bandwidth

$\Upsilon_{\rm{amp}}(\nu,\nu_{0})=\Upsilon_{\rm{amp}}(\nu_{0},\nu_{0})\times\dfrac{\sqrt{(\nu-\nu_{0})^{2}+(\Lambda/2)^{2}}}{\Lambda/2},$
(S32)

The magnetic sensitivity of 87Rb magnetometer assisted with the 129Xe spin-
based amplifier can be achieved over different bias fields.

###  IV.2 Axion-like dark matter search experiment

####  IV.2.1 The timing of axion-like dark matter searches

The spin-based amplifier can operate with high amplification factor within the
bandwidth of 0.052 Hz centered at 129Xe Larmor frequency, as discussed in Sec.
II.3. Due to the narrow bandwidth, the spin-based amplifier is only highly
sensitive to a narrow range of axion-like particle masses. For example, the
current 0.052 Hz bandwidth corresponds to 0.2 feV bandwidth of axion-like
particle mass (see Eq. S21). To search for a broad range of axion mass, we can
move the high-sensitivity bandwidth window by changing bias field, which can
accordingly change the 129Xe Larmor frequency of the spin-based amplifier. The
detailed procedure of our search experiment is shown in Fig. S10. In our
search experiment, the scanning step of 129Xe Larmor frequency is 0.0358 Hz,
which is slightly smaller than the spin-based amplifier bandwidth (i.e., 0.052
Hz). In principle, we could choose a smaller scanning step, but it would take
more time to complete the entire search. The entire 129Xe Larmor frequency
range is from 2 Hz to 180 Hz, corresponding to axion-like particle masses from
8.3 feV to 0.74 peV. In each search, we take 100 s of data and calculate the
detection threshold of axion-like dark matter (see below). Last, we combine
the all independent search results to achieve the entire detection threshold
from 2 Hz to 180 Hz. The detailed data processing is described in the
following. As an initial search, we also report experimental search for 100
probed axion-like particle mass windows in the range of 10 feV to 1 peV, where
each result is from 5 h of data.

Figure S10: The timing of axion-like dark matter search. A single search
consists of 50 s waiting time and 100 s sampling time. Before each search, the
129Xe Larmor frequency is increased by 0.0358 Hz per step by changing the bias
field. Figure S11: Data analysis of axion-like dark matter search. (a) A small
window of the experimental power spectrum. In this situation, the bias field
is set to 1681.3 nT and the corresponding 129Xe Larmor frequency is
$\nu_{0}\approx 19.839$ Hz. The searched axion mass window is centered at
82.05 feV with the width of 0.15 feV. (b) White-noise power histogram, which
can be modelled as a chi-squared distribution with 2 degrees of freedom. (c)
Cumulative distribution of chi-square distribution and fit to the exponential
distribution $1-e^{-\chi/2\sigma^{2}}$. The fit yields a white-noise variance
$\sigma\approx 4.241\times 10^{-7}$ V. (d) Limits ($95\%$ confidence level) on
the coupling strength of nucleons with the axion-like dark matter within a
mass range centered at $m_{a}=82.05$ feV and with the width of 0.15 feV. The
limits are ahieved by five hour measurement at $\nu_{0}\approx 19.839$ Hz of
the spin-based amplifier.

####  IV.2.2 Detection threshold determination

To define the detection threshold, we use the standard $p$-value hypothesis
test. Detection threshold for a hypothesis test at a frequency $\nu$ is
defined by the probability to find data which are less compatible than the
observed $z_{0}$. In the following, we use one of 5-h search data
($\nu_{0}\approx 19.839$ Hz, corresponding to $m_{a}=82.05$ feV) as an example
to explain our data procedure. We calculate the power spectrum of the 5-h
data. Figure S11a shows the amplitude spectrum (the square root of the power
spectral density). We analyze the power spectral density in $\pm 2$ Hz region
centered at the corresponding Larmor frequency. The histogram of the power
spectral density in this frequency is shown in Fig. S11b with the bin width of
$1.62\times 10^{-13}$ V2. The histogram can be modeled as the chi-squared
distribution with two degrees of freedom. The corresponding cumulative
distribution function (CDF) follows an exponential distribution [15], as shown
in Fig. S11c. Specifically, the probability that the measured power $P(\nu)$
is smaller than $z_{0}$ is

$\textrm{Pr}\\{\textrm{P}(\nu)<z_{0}\\}=1-e^{-z_{0}/2\sigma^{2}},$ (S33)

where $\sigma$ is the standard deviation. Hence, we can fit the experimental
data to obtain the standard deviation $\sigma$. The red curve is the fit and
agrees well with the experiment, yielding $\sigma\approx 4.2\times 10^{-7}$ V.
For 100 s sampling, we adopt the same data processing method.

We now define a signal power threshold $z_{0}$ such that, if at the probed
frequency $\nu$, the power of a spectrum containing both the signal and noise,
is greater than $z_{0}$, then peak has a $p_{0}$ probability of being induced
by noise fluctuation,

$p_{0}:=\textrm{Pr}\\{\textrm{P}(\nu)>z_{0}\\}=e^{-z_{0}/2\sigma^{2}}.$ (S34)

Based on the above equation, we can achieve

$z_{0}=-2\sigma^{2}\textrm{ln}(p_{0}).$ (S35)

Thus, when $p_{0}=0.05$ (corresponding to 95$\%$ confidence level), we find
the signal detection threshold is

$z_{0}=5.99\sigma^{2}.$ (S36)

Because our analysis is done on the amplitude spectrum, we define the 95$\%$
confidence detection threshold from the amplitude spectrum as

$A_{\textrm{th}}=2.44\sigma.$ (S37)

Here, we should note that the unit of $A_{\textrm{th}}$ is voltage and can be
converted into magnetic field according to the calibration (2) in Sec. IV.1.
Further, to obtain the constraint on the coupling parameter
$\textrm{g}_{\rm{aNN}}$ (the coupling of axion-like particles to nucleons), we
can use Eq. S29 and $A_{\rm{th}}$ to achieve the corresponding
$\textrm{g}_{\rm{aNN}}$. For example, our experiment constrains the parameter
space describing the coupling to axion-like dark matter,
$|\textrm{g}_{\rm{aNN}}|<3.5\times 10^{-9}$ GeV-1 at 82.05 feV axion mass.

In each search experiment (corresponding to a bias field), not only
$\textrm{g}_{\rm{aNN}}$ can be constrained at the corresponding Larmor
frequency, but also $\textrm{g}_{\rm{aNN}}$ can be constrained at near-
resonance frequencies, because the spin-based amplifier still has a high
sensitivity to axion-like signals whose oscillating frequencies are within the
amplifier bandwidth. We have

$\textrm{g}_{\textrm{aNN}}(\nu,\nu_{0})=\textrm{g}_{\textrm{aNN}}(\nu_{0},\nu_{0})\times\dfrac{\sqrt{(\nu-\nu_{0})^{2}+(\Lambda/2)^{2}}}{\Lambda/2},$
(S38)

where $\textrm{g}_{\textrm{aNN}}(\nu_{0},\nu_{0})$ is the axion-nucleon
coupling limit at the Larmor frequency $\nu_{0}$ and
$\textrm{g}_{\textrm{aNN}}(\nu,\nu_{0})$ is the axion-nucleon limit at the
near-resonance frequency. Figure S11d provides the constraint limit on
$\textrm{g}_{\rm{aNN}}$ over the spin-based amplifier bandwidth. After
obtaining the threshold for the spin-based amplifier bandwidth range, we scan
the bias field to obtain the threshold of axion-like dark matter from 2 Hz and
180 Hz. The limit for $\textrm{g}_{\textrm{aNN}}$ is shown in Fig. 4b of the
main text.

##  V Constraint plots for other possible bosonic dark matter fields

This section presents three other possible couplings between bosonic dark
matter fields and nuclear spins [16, 17, 11, 8, 10]: the quadratic wind
coupling with the axions, the couplings to dark electric field and magnetic
field mediated by spin-1 bosons such as dark photons [8, 10]. The Hamiltonian
of the three couplings can be written as

$\displaystyle\mathcal{H}_{\textrm{quad}}\approx
2\textrm{g}_{\textrm{quad}}^{2}\hbar^{2}c^{2}\dfrac{\rho_{\textrm{DM}}}{2\pi\nu_{a}}\sin(2\pi\nu_{a}t)\bm{v}\cdot\textbf{I}_{\rm{N}},$
(S39)
$\displaystyle\mathcal{H}_{\rm{dMDM}}\approx\textrm{g}_{\rm{dMDM}}\sqrt{2\hbar^{3}c\rho_{\rm{DM}}}\sin(2\pi\nu_{a}t)\bm{v}\cdot\textbf{I}_{\rm{N}},$
(S40)
$\displaystyle\mathcal{H}_{\rm{dEDM}}\approx\textrm{g}_{\rm{dEDM}}\sqrt{2\hbar^{3}c\rho_{\rm{DM}}}\sin(2\pi\nu_{a}t)\bm{v}\cdot\textbf{I}_{\rm{N}}/|\bm{v}|,$
(S41)

where $\textrm{g}_{\textrm{quad}}$ is the quadratic coupling strength between
axions and nuclear spins, $\textrm{g}_{\rm{dMDM}}$ and
$\textrm{g}_{\rm{dEDM}}$ respectively parameterizes the coupling strengths
between nuclear spins with the dark magnetic fields and electric fields.
$\bm{v}$ is taken as the same $\bm{v}$ in Sec. III. Based on the same analysis
and calculation in Sec. III, we can derive the effective magnetic field
induced by the bosonic dark matter fields

$\displaystyle\textbf{B}_{\textrm{quad}}[\textrm{T}]$ $\displaystyle\approx
10^{-3}\dfrac{\textrm{g}_{\rm{quad}}^{2}}{\gamma_{n}}2\hbar
c^{3}\dfrac{\rho_{\rm{DM}}}{2\pi\nu_{a}}\sin(2\pi\nu_{a}t)\times\dfrac{\bm{v}}{|\bm{v}|},$
$\displaystyle\approx 2.9\cdot
10^{-4}\dfrac{\rm{g}_{\rm{quad}}^{2}[(\rm{GeV}^{-1})^{2}]}{\rm{g}_{n}}\dfrac{\sin(2\pi\nu_{a}t)}{2\pi\nu_{a}}\times\sqrt{1-3/4\cos^{2}(\Omega
t-\eta)},$ (S42) $\displaystyle\textbf{B}_{\textrm{dMDM}}[\textrm{T}]$
$\displaystyle\approx
10^{-3}\dfrac{\textrm{g}_{\rm{dMDM}}}{\hbar\gamma_{n}}\sqrt{2\hbar^{3}c^{3}\rho_{\rm{DM}}}\sin(2\pi\nu_{a}t)\times\dfrac{\bm{v}}{|\bm{v}|},$
$\displaystyle\approx
10^{-7}\dfrac{\rm{g}_{\rm{dMDM}}[\rm{GeV}^{-1}]}{\rm{g}_{n}}\sin(2\pi\nu_{a}t)\times\sqrt{1-3/4\cos^{2}(\Omega
t-\eta)},$ (S43) $\displaystyle\textbf{B}_{\textrm{dEDM}}[\textrm{T}]$
$\displaystyle\approx\dfrac{\textrm{g}_{\rm{dEDM}}}{\hbar\gamma_{n}}\sqrt{2\hbar^{3}c\rho_{\rm{DM}}}\sin(2\pi\nu_{a}t)\times\dfrac{\bm{v}}{|\bm{v}|},$
$\displaystyle\approx
10^{-4}\dfrac{\rm{g}_{\rm{dEDM}}[\rm{GeV}^{-1}]}{\rm{g}_{n}}\sin(2\pi\nu_{a}t)\times\sqrt{1-3/4\cos^{2}(\Omega
t-\eta)}.$ (S44)

Thus, based on the same procedures to derive the constraints level of the
axion-nucleon coupling (see Sec. IV), the constraint plots for the couplings
between the nuclear spins and the bosonic dark matter fields are shown in Fig.
S12 ($\textrm{g}_{\textrm{dMDM}}$), Fig. 4c ($\textrm{g}_{\textrm{dEDM}}$) and
Fig. 4d ($\textrm{g}_{\textrm{quad}}$) of the main text.

Figure S12: Results of axion-like dark matter search. Limits on dark
photon–nucleon dMDM coupling in the mass range from 8.3 feV to 744 feV. The
blue-shade is excluded by our measurements (100 s for each run) at the $95\%$
confidence level. The red lines show our advanced sensitivity (five hours for
each run) for dark photon-nucleon interactions at 100 probed dark-photon
masses with a window width of 0.15 feV. The grey line shows the limit given by
the CASPEr-ZULF experiment [10]. The horizontal black line shows the
astrophysical limit from supernova SN1987A cooling [18, 19].

##  VI Projected sensitivity reach

Our experimental results demonstrate the feasibility of using spin-based
amplifier to search for axion-like dark matter. We anticipate that the
experimental sensitivity to axion-like dark matter can be further improved by
several orders of magnitude. The main factors limiting the sensitivity of our
search for axion-like dark matter are: (1) the amplification factor of the
spin-based amplifier, and (2) the magnetic sensitivity of the alkali-metal
magnetometer itself. Considering this, we plan to employ 3He-K systems as a
spin-based amplifier, because 3He spins have much longer coherence time
($T_{2}\sim 1000$ s) and larger gyromagnetic ratio than those of 129Xe
($T_{2}\sim 20$ s in this work). In the following, we estimate the
corresponding amplification factor and axion-like dark matter search
sensitivity.

We theoretically consider the amplification effect on 3He-K system. Compared
to 87Rb magnetometer, K magnetometer has demonstrated higher sensitivity [20,
4] and a 3He-39K spin-based amplifier can achieve larger amplification factor
as stated below. Based on the Eq. S18, we predict the amplification factor and
higher magnetic sensitivity of such a spin-based amplifier. To calculate the
amplification factor, we need to calculate the term $\lambda M^{n}P^{n}_{0}$.
Using the SI units, $\lambda M^{n}P^{n}_{0}$ can be expressed as [5, 4]

$\displaystyle\dfrac{2\mu_{0}}{3}\kappa_{0}M^{n}P_{0}^{n}=\dfrac{2\mu_{0}}{3}\kappa_{0}g\mu
P_{0}^{n}[\rm{N}],$ (S45)

where $\kappa_{0}$ is the enhancement factor for the given alkali metal and
noble gas pairs (for example, for 3He and K, $\kappa_{0}\approx 5.9$), $g$ is
the Land$\acute{e}$ factor for noble gas nucleus, $P^{n}_{0}$ is the
polarization of noble gas nucleus, $[\rm{N}]$ is the number density of the
noble gas, and $\mu=1.07\times 10^{-26}~{}\textrm{J/T}$ is the magnetic
moments of 3He nucleus. We use Eq. S45 to estimate the 3He effective magnetic
field and numerically calculate the amplification factor. Based on the
demonstrated experimental parameters in ref. [4], $P_{0}^{n}\approx 0.02$ is
the typical polarization of the 3He gas, $[\rm{He}]\approx 1.88\times
10^{26}~{}\textrm{m}^{-3}$ corresponds to the number density with 7 amg 3He.
We obtain $\lambda M^{n}P_{0}^{n}\approx 2.0~{}\textrm{mG}$, which is in good
agreement with 2.2 mG demonstrated in ref. [4]. According to Eq. S15, the
amplification factor of 3He-K spin-based amplifier is estimated to be
$\eta\approx 10^{4}$. Compared to 129Xe-Rb systems, the amplification factor
of 3He-K systems yields about two orders of magnitude improvement.

Based on the above discussion, one can evaluate the experimental sensitivity
to axion-like dark matter. We still need to know the sensitivity of the K
magnetometer itself. The 3He-K system has five orders smaller spin-destruction
cross section than that of 129Xe-Rb sytems, and thus the K magnetometer can
still achieve a femtotesla sensitivity
$\sim$$1~{}\textrm{fT}/\textrm{Hz}^{1/2}$, as demonstrated in ref. [4]; thus,
the magnetic sensitivity based on 3He spin amplifier could probably reach
$1~{}\textrm{aT}/\textrm{Hz}^{1/2}$ within the amplifier bandwidth. We note
that the experimental sensitivity to $\rm{g}_{\rm{aNN}}$ scales as $t^{-1/2}$
as a function of measurement time $t$. When the measurement time is five
hours, the search sensitivity of $|\textrm{g}_{\textrm{aNN}}|$ can reach
$\sim$$10^{-13}~{}\textrm{GeV}^{-1}$, which improves over the astrophysical
limit (from supernova SN1987A cooling [18, 19]) by at least four orders of
magnitude. For dark photon-nucleon interactions, the experimental sensitivity
could reach $\rm{g}_{\rm{dEDM}}\sim 10^{-17}$ GeV-1 and
$\rm{g}_{\rm{dMDM}}\sim 10^{-13}$ GeV-1, which are far beyond the
astrophysical limits. Using such a spin-based amplifier, it is possible to
explore new physics beyond the standard model.

If an alarm signal is observed, experiment would involve re-scans of the
candidate frequencies to verify if they are statistical deviations or real
axion-like signals. The re-scans would be particularly time-demanding if the
experiment is based on long-term search data. Alternatively, we plan to build
a network of synchronized spin-based amplifiers to differentiate false alarms
from true detection events. Because the apparatus used in our experiment is
small-scale and inexpensive, we could quickly build such a network comprising
of, for example, ten spin-based amplifiers. Each sensor in the network
independently performs the search for axion-like signals in the same frequency
range from 1 Hz to 200 Hz. Using the network, claims of detecting axion-like
dark matter could be made for candidates that are all above the thresholds for
every sensor measurement at the searched frequency. This leads to a
significant decrease in false-alarm probability. For example, with a network
of ten spin-based amplifiers, the false-alarm probability becomes
$p_{0}^{10}\approx 10^{-13}$ and the false-alarm rate is $\sim 10^{-7}$ alarms
per experiment. In addition to reducing false-alarm rate, such a sensor
network is promising to compose an exotic field telescope array for multi-
messenger astronomy, as recently proposed [21, 16], and address the stochastic
fluctuations of bosonic dark matter [22].

##  VII Testing the data analysis procedure with simulated axion signals

Our data analysis procedure is tested by adding simulated axion-like dark
matter signals into our experimental data. As a demonstration, an on-resonance
axion-like signal and a signal with a small frequency deviation from resonance
frequency (corresponding to near-resonant case but still within the bandwidth)
are injected into experimental data. For example, we add two simulated axion-
like signals at the Compton frequency of 19.84 Hz (on resonance) and 19.858 Hz
(near-resonance) and with the coupling strength
$\textrm{g}_{\rm{aNN}}$=25.18$\times 10^{-9}$ GeV-1. These simulated signals
can generate an effective magnetic field of 1.63$\times 10^{-3}$ pT. The
corresponding power spectrum is shown in Fig. S13a, where the axion-like
signals are clearly visible at 19.84 Hz and 19.858 Hz. As shown in Fig. S13b,
there are two notable points (blue and red circles) that are significantly
above the 95$\%$ threshold in the the histogram. Based on the response of the
spin-based amplifier, we recover the effective magnetic field as 1.63$\times
10^{-3}$ pT and the coupling strength $\rm{g}_{\rm{aNN}}$=25.18$\times
10^{-9}$ GeV-1, which agrees well with the injected simulated axion-like
signals. It indicates that our experiment can distinguish the simulated axion-
like dark matter signal within the spin-based amplifier bandwidth.
Furthermore, we inject simulated signals from 1$\sigma$ to 20$\sigma$
significance corresponding coupling strength from 1.26$\times 10^{-9}$ GeV-1
to 25.18$\times 10^{-9}$ GeV-1. The recovery axion-nucleon coupling strength
are shown with blue squares in Fig. S13c, which are fitted with a linear fit
(red line). The slope of the linear fit is 0.996, yielding that the couplings
match the injected signals.

Figure S13: Injecting simulated axion-like dark matter signals into the
experimental data. (a) A small window of the experimental power spectrum, with
two numerically-generated synthetic axion-like signals added. The Compton
frequencies of such axion-like signals are set to 19.84 Hz and 19.858 Hz and
the axion-nucleon coupling strength is set to
$\textrm{g}_{\rm{aNN}}$=25.18$\times 10^{-9}$ GeV-1. The resonant frequency of
the spin-based amplifier is 19.84 Hz. There are two visible axion-like signals
at their corresponding Compton frequencies. The peak amplitude at resonant
19.84 Hz is larger than that at near-resonant 19.858 Hz. (b) The histogram of
the power spectral data within the bin shown in (a). The vertical dashed line
marks the 5.99$\sigma^{2}$ threshold ($95\%$ confidence level) for an event
candidate. The injected synthetic axion-like signal was detected with
20$\sigma$ significance. The inset shows that the magnetic sensitivity is
about 9.95 fT/Hz1/2 at the resonant frequency and 15.7 fT/Hz1/2 at the near-
resonant frequency. (c) Verification of detection and coupling strength
recovery of synthetic axion-like signals, with coupling strengths between
$1.26\times 10^{-9}$GeV-1 and $25.18\times 10^{-8}$ GeV-1. The coupling
strengths recovered from detected signals are shown as blue squares, and the
coupling matching the injected signal is shown as the red line.

## References

  * [1] Wu, Z., Kitano, M., Happer, W., Hou, M. & Daniels, J.  Optical determination of alkali metal vapor number density using faraday rotation.  _Appl. Opt._ 25, 4483–4492 (1986).
  * [2] Opechowski, W.  Magneto-optical effects and paramagnetic resonance.  _Rev. Mod. Phys._ 25, 264 (1953).
  * [3] Jiang, M. _et al._ Interference in atomic magnetometry.  _Adv. Quantum Technol._ 3, 2000078 (2020).
  * [4] Kornack, T., Ghosh, R. & Romalis, M.  Nuclear spin gyroscope based on an atomic comagnetometer.  _Phys. Rev. Lett._ 95, 230801 (2005).
  * [5] Walker, T. G. & Happer, W.  Spin-exchange optical pumping of noble-gas nuclei.  _Rev. Mod. Phys._ 69, 629 (1997).
  * [6] Kimball, D. F. J. _et al._ Overview of the cosmic axion spin precession experiment (CASPEr).  In _Microwave Cavities and Detectors for Axion Research_ , 105–121 (Springer, 2020).
  * [7] Abel, C. _et al._ Search for axionlike dark matter through nuclear spin precession in electric and magnetic fields.  _Phys. Rev. X_ 7, 041034 (2017).
  * [8] Wu, T. _et al._ Search for axionlike dark matter with a liquid-state nuclear spin comagnetometer.  _Phys. Rev. Lett._ 122, 191302 (2019).
  * [9] Smorra, C. _et al._ Direct limits on the interaction of antiprotons with axion-like dark matter.  _Nature_ 575, 310–314 (2019).
  * [10] Garcon, A. _et al._ Constraints on bosonic dark matter from ultralow-field nuclear magnetic resonance.  _Sci. Adv._ 5, eaax4539 (2019).
  * [11] Graham, P. W. _et al._ Spin precession experiments for light axionic dark matter.  _Phys. Rev. D_ 97, 055006 (2018).
  * [12] Graham, P. W. & Rajendran, S.  New observables for direct detection of axion dark matter.  _Phys. Rev. D_ 88, 035023 (2013).
  * [13] Marsh, D. J.  Axion cosmology.  _Phys. Rep._ 643, 1–79 (2016).
  * [14] Kosteleckỳ, V. A. & Lane, C. D.  Constraints on lorentz violation from clock-comparison experiments.  _Phys. Rev. D_ 60, 116010 (1999).
  * [15] Walck, C.  Hand-book on statistical distributions for experimentalists. Internal Report SUF–PFY/96–01 (Univ. Stockholm, 2007).
  * [16] Pospelov, M. _et al._ Detecting domain walls of axionlike models using terrestrial experiments.  _Phys. Rev. Lett._ 110, 021803 (2013).
  * [17] Graham, P. W., Irastorza, I. G., Lamoreaux, S. K., Lindner, A. & van Bibber, K. A.  Experimental searches for the axion and axion-like particles.  _Annu. Rev. Nucl. Part. Sci._ 65, 485–514 (2015).
  * [18] Vysotsskii, M., Zel’Dovich, Y. B., Khlopov, M. Y. & Chechetkin, V.  Some astrophysical limitations on the axion mass.  _JETP Lett._ 27, 533 (1978).
  * [19] Raffelt, G. G.  Astrophysical axion bounds.  In _Axions_ , 51–71 (Springer, 2008).
  * [20] Kominis, I., Kornack, T., Allred, J. & Romalis, M. V.  A subfemtotesla multichannel atomic magnetometer.  _Nature_ 422, 596–599 (2003).
  * [21] Dailey, C. _et al._ Quantum sensor networks as exotic field telescopes for multi-messenger astronomy.  _Nat. Astron.,_ 1–9 (2020).
  * [22] Centers, G. P. _et al._ Stochastic fluctuations of bosonic dark matter. Preprint at https://arxiv.org/abs/1905.13650 (2019).