# A Hunt for Magnetic Signatures of Hidden-Photon and Axion Dark Matter in the
Wilderness

Ibrahim A. Sulai<EMAIL_ADDRESS>Department of Physics &
Astronomy, Bucknell University, Lewisburg, Pennsylvania 17837, USA  Saarik
Kalia<EMAIL_ADDRESS>School of Physics & Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  Ariel Arza  Tsung-Dao Lee Institute (TDLI),
Shanghai Jiao Tong University, Shanghai 200240, China  Itay M. Bloch  Berkeley
Center for Theoretical Physics, University of California, Berkeley, CA 94720,
USA  Theory Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA  Eduardo Castro Muñoz  Department of Physics & Astronomy, Oberlin College,
Oberlin, Ohio 44074, USA  Christopher Fabian  Department of Physics &
Astronomy, Bucknell University, Lewisburg, Pennsylvania 17837, USA  Michael A.
Fedderke  The William H. Miller III Department of Physics and Astronomy, The
Johns Hopkins University, Baltimore, MD 21218, USA  Madison Forseth
Department of Physics, California State University – East Bay, Hayward,
California 94542-3084, USA  Brian Garthwaite  Department of Physics &
Astronomy, Bucknell University, Lewisburg, Pennsylvania 17837, USA  Peter W.
Graham  Stanford Institute for Theoretical Physics, Department of Physics,
Stanford University, Stanford, CA 94305, USA  Kavli Institute for Particle
Astrophysics & Cosmology, Stanford University, Stanford, CA 94305, USA  Will
Griffith  Department of Computing, Math and Physics, Messiah University,
Mechanicsburg, PA 17055, USA  Erik Helgren  Department of Physics, California
State University – East Bay, Hayward, California 94542-3084, USA  Katie
Hermanson  Department of Physics, California State University – East Bay,
Hayward, California 94542-3084, USA  Andres Interiano-Alvarado  Department of
Physics, California State University – East Bay, Hayward, California
94542-3084, USA  Brittany Karki  Department of Physics, California State
University – East Bay, Hayward, California 94542-3084, USA  Abaz Kryemadhi
Department of Computing, Math and Physics, Messiah University, Mechanicsburg,
PA 17055, USA  Andre Li  Department of Physics, California State University –
East Bay, Hayward, California 94542-3084, USA  Ehsanullah Nikfar  Department
of Physics & Astronomy, Oberlin College, Oberlin, Ohio 44074, USA  Jason E.
Stalnaker  Department of Physics & Astronomy, Oberlin College, Oberlin, Ohio
44074, USA  Yicheng Wang  Department of Physics & Astronomy, Bucknell
University, Lewisburg, Pennsylvania 17837, USA  Derek F. Jackson Kimball
<EMAIL_ADDRESS>Department of Physics, California State
University – East Bay, Hayward, California 94542-3084, USA

###### Abstract

Earth can act as a transducer to convert ultralight bosonic dark matter
(axions and hidden photons) into an oscillating magnetic field with a
characteristic pattern across its surface. Here we describe the first results
of a dedicated experiment, the Search for Non-Interacting Particles
Experimental Hunt (SNIPE Hunt), that aims to detect such dark-matter-induced
magnetic-field patterns by performing correlated measurements with a network
of magnetometers in relatively quiet magnetic environments (in the wilderness
far from human-generated magnetic noise). Our experiment constrains parameter
space describing hidden-photon and axion dark matter with Compton frequencies
in the 0.5-5.0 Hz range. Limits on the kinetic-mixing parameter for hidden-
photon dark matter represent the best experimental bounds to date in this
frequency range.

##  I Introduction

Understanding the nature of dark matter is of paramount importance to
astrophysics, cosmology, and particle physics. A well-motivated hypothesis is
that the dark matter consists of ultralight bosons (masses $\ll$ 1 eV$/c^{2}$)
such as hidden photons, axions, or axion-like particles (ALPs) [1, 2, 3]. If
ultralight bosons are the dark matter, under reasonable assumptions111Here we
assume models where the self-interactions among the bosons are sufficiently
feeble that they do not collapse into large composite structures (such as
boson stars [4]). Therefore, the bosons can be treated as an ensemble of
independent particles described by the standard halo model (SHM) of dark
matter [5, 6, 7]. the ensemble of virialized bosons constituting the dark
matter halo has extremely large mode-occupation numbers and can be well
described as a stochastic classical field [8, 9, 10, 11, 12].

Ultralight bosonic fields can couple to Standard Model particles through
various “portals” [13, 14], one of which is the interaction between the
ultralight bosonic dark matter (UBDM) and the electromagnetic field. Several
ongoing laboratory experiments employ sensitive magnetometers located within
controlled magnetic environments to search for electromagnetic signatures of
UBDM; see, for example, Refs. [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. As
noted in Refs. [26, 27, 28], the conceptual framework for UBDM-to-photon
conversion upon which these aforementioned laboratory searches are based also
applies to Earth as a whole. For hidden-photon dark matter, the non-conducting
atmosphere sandwiched between the conductive Earth interior and the ionosphere
acts as a transducer to convert the hidden photon field into a real magnetic
field, just as laboratory-scale shields act as transducers in lumped-element
or resonant-cavity experiments [23, 24, 25]. For axion dark matter, Earth’s
geomagnetic field causes axion-to-photon conversion via the inverse Primakoff
effect [29, 30], playing the role of the applied magnetic field in laboratory-
scale axion haloscope experiments [15, 16, 17, 18, 19, 20, 21, 22]. Thus,
unshielded magnetometers can be used to search for ambient oscillating
magnetic fields generated by UBDM.

In this paper we describe initial results of the “Search for Non-Interacting
Particles Experimental Hunt” (SNIPE Hunt [31]): a campaign to search for
axion222We use the term “axion” as a generic descriptor of both QCD axions
(that solve the strong-CP problem) and axion-like particles (ALPs). and
hidden-photon dark matter using magnetometers located in the “wilderness”
(away from the high levels of magnetic noise associated with urban
environments [32, 33]). This work extends to higher axion/hidden-photon
Compton frequencies (covering the range from 0.5-5 Hz) than earlier analyses
of archival data from the SuperMAG network of magnetometers [34, 35, 36]
published in Refs. [27, 28]. In this frequency range, the dominant magnetic
field noise sources are anthropogenic [37], so we anticipate that the
sensitivity to UBDM can be drastically enhanced by measuring in a remote
location.

The rest of this paper is structured as follows. Section II reviews the model
developed in Refs. [26, 28] to predict the global magnetic field patterns
induced by hidden-photon and axion dark matter and used to interpret our data.
In Sec. III, we discuss the experimental setup for the magnetometers that
measured the magnetic fields at three different locations in July 2022 as well
as the time and frequency characteristics of the acquired data. In Sec. IV,
the data analysis procedure is described, which is closely based on that
presented in Refs. [27, 28]. Section IV is subdivided into one subsection on
the hidden-photon dark-matter analysis and another on the axion dark-matter
analysis; in both cases no evidence of a dark-matter-induced magnetic signal
was discovered, so each subsection concludes by summarizing the constraints
obtained on relevant parameters. In Sec. V, we summarize the next steps for
the SNIPE Hunt research program, namely developing and carrying out an
experiment for higher Compton frequencies with more sensitive magnetometers.
Finally, in our conclusion we summarize results and compare them to other
experiments and observational limits.

##  II Dark-Matter Signal

First, we review relevant features of the theory motivating our hidden-photon
dark-matter search. The hidden photon is associated with an additional $U(1)$
symmetry, beyond that corresponding to electromagnetism, which is a common
feature of beyond-the-Standard-Model theories, such as string theory [38]. In
our case, we are interested in hidden photons that kinetically mix with
ordinary photons [39]. This allows hidden and ordinary photons to interconvert
via a phenomenon akin to neutrino mixing [40]; i.e., the mass (propagation)
and interaction eigenstates are misaligned. Hidden photons possess a non-zero
mass $m_{A^{\prime}}$ and can be generated in the early universe (see, for
example, Refs. [41, 42, 43, 44]), which means that they have the right
characteristics to be wave-like dark matter [45]. A useful way to understand
the impact of the existence of hidden-photon dark matter on electrodynamics is
to write the Lagrangian describing real and hidden photons in the
“interaction” basis [24, 26]:333Throughout, we use natural units where
$\hbar=c=1$.

$\displaystyle\mathcal{L}\supset-\frac{1}{4}{\left[F_{\mu\nu}F^{\mu\nu}+{\left(F^{\prime}\right)}_{\mu\nu}{\left(F^{\prime}\right)}^{\mu\nu}\right]}+\frac{1}{2}m_{A^{\prime}}^{2}{\left(A^{\prime}\right)}_{\mu}{\left(A^{\prime}\right)}^{\mu}+\varepsilon
m_{A^{\prime}}^{2}{\left(A^{\prime}\right)}^{\mu}A_{\mu}-J_{\mathrm{EM}}^{\mu}A_{\mu}\>,$
(1)

where only terms up to first order in the kinetic mixing parameter
$\varepsilon\ll 1$ are retained. In Eq. (1), $F_{\mu\nu}$ is the field-
strength tensor for the “interacting” mode of the electromagnetic field that
couples to charges, ${\left(F^{\prime}\right)}_{\mu\nu}$ is the field-strength
tensor for the “sterile” mode that does not interact with charges, $A_{\mu}$
is the four-potential for the interacting mode,
${\left(A^{\prime}\right)}_{\mu}$ is the four-potential for the sterile mode,
and $J_{\mathrm{EM}}^{\mu}$ is the electromagnetic four-current density. In
our case of interest, the hidden-photon dark-matter field in the vicinity of
Earth is a coherently oscillating vector field with random polarization:444In
this work, we assume that both the hidden-photon phase and its polarization
state randomize on the coherence timescale. It is also possible, depending on
the production mechanism and subsequent structure-formation processing, that
the hidden-photon polarization state could be fixed in inertial space; see,
e.g., the discussions in Refs. [46, 47]. We do not explicitly consider this
case in this work; a closely related, but different, analysis would need to be
undertaken. However, absent accidental geometrical cancellations that are made
unlikely by virtue of the length of the data-taking period compared to Earth’s
sidereal rotational period and the widely separated geographical locations of
the magnetic-field stations on which we report, limits in that case are
expected to be of the same order of magnitude as those we obtain.

$\displaystyle\bm{A}^{\prime}(\bm{r},t)\approx\frac{\sqrt{2\rho_{\mathrm{DM}}}}{m_{A^{\prime}}}e^{-im_{A^{\prime}}t}\sum_{i=1}^{3}\xi_{i}{\left(\bm{r},t\right)}\bm{\hat{n}}_{i}e^{i\phi_{i}(\bm{r},t)}~{},$
(2)

where $\bm{A}^{\prime}$ is the sterile vector potential,
$\rho_{\mathrm{DM}}\approx 0.3~{}{\rm{GeV/cm^{3}}}$ is the local dark-matter
density [48], $\bm{\hat{n}}_{i}$ are a set of orthonormal unit vectors,
$\xi_{i}(\bm{r},t)$ are slowly varying $\mathcal{O}(1)$ amplitudes, and
$\phi_{i}(\bm{r},t)$ are slowly varying random phases. Both the amplitudes
$\xi_{i}{\left(\bm{r},t\right)}$ and phases $\phi_{i}(\bm{r},t)$ of the
hidden-photon dark-matter field change stochastically on length scales given
by the dark-matter coherence length,

$\displaystyle\ell{{}_{\mbox{\scriptsize
coh}}}\approx\frac{2\pi}{m_{A^{\prime}}v{{}_{\mbox{\scriptsize DM}}}}~{},$ (3)

and time scales given by the coherence time of the field,

$\displaystyle\tau{{}_{\mbox{\scriptsize
coh}}}\approx\frac{\ell{{}_{\mbox{\scriptsize coh}}}}{v{{}_{\mbox{\scriptsize
DM}}}}\approx\frac{2\pi}{m_{A^{\prime}}v_{\mathrm{DM}}^{2}}~{},$ (4)

where $v{{}_{\mbox{\scriptsize DM}}}\sim 10^{-3}$ is the characteristic
dispersion (virial) velocity of the dark matter in the vicinity of Earth [49,
7]. Note that the timelike component of the four-potential
${\left(A^{\prime}\right)}^{\mu}$ is suppressed relative to the spacelike
component (the vector potential $\bm{A}^{\prime}$) by $\sim
v{{}_{\mbox{\scriptsize DM}}}\sim 10^{-3}$. From inspection of Eq. (1), it can
be seen that the physical effects due to the hidden-photon dark-matter field
${\left(A^{\prime}\right)}^{\mu}$ are to leading order the same as those
generated by an effective current density

$\displaystyle\bm{J}_{A^{\prime}}=-\varepsilon
m_{A^{\prime}}^{2}\bm{A}^{\prime}~{}.$ (5)

Inside a good conductor, the interacting mode vanishes, $F_{\mu\nu}=0$ and
$A_{\mu}=0$, whereas the sterile mode can propagate into a conducting region
with essentially no perturbation. Outside a conducting region, the effective
current density due to the sterile mode acts to generate a non-zero
interacting mode. These effects, where Earth’s conducting interior and the
conducting ionosphere provide relevant boundary conditions, give rise to the
oscillating magnetic-field pattern we seek to measure in our experiment, as
described in detail in Ref. [26].

The second theoretical scenario we consider is the hypothesis that the dark
matter consists primarily of axions [50, 51, 52, 53, 54, 55]. Axions are
pseudoscalar particles arising from spontaneous symmetry breaking at a high
energy scale associated, for example, with grand unified theories (GUTs) or
even the Planck scale [56]. Combined with explicit symmetry breaking at lower
energy scales, such pseudoscalar particles acquire small masses ($\ll
1~{}{\rm{eV}}$) and couplings to Standard Model particles and fields [2]. Like
hidden photons, axions are ubiquitous features of beyond-the-Standard-Model
theories [57, 58, 59, 60], and have all the requisite characteristics to be
the dark matter [1, 2, 3]. The focus of our experiment is the axion-to-photon
coupling which is described by the Lagrangian:

$\displaystyle\mathcal{L}\supset-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\frac{1}{2}{\left(\partial_{\mu}a\right)}^{2}-\frac{1}{2}m_{a}^{2}a^{2}+\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}~{},$
(6)

where $a$ is the axion field, $m_{a}$ is the axion mass, $g_{a\gamma}$
parameterizes the axion–photon coupling, and $\tilde{F}^{\mu\nu}$ is the dual
field-strength tensor. The last term appearing in Eq. (6) describes the
interaction between the axion and electromagnetic fields:

$\displaystyle\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}=-g_{a\gamma}a\bm{E}\cdot\bm{B}~{},$
(7)

where $\bm{E}$ and $\bm{B}$ are the electric and magnetic fields. In the non-
relativistic limit, the leading-order correction to Maxwell’s equations
arising from the existence of the axion–photon coupling described by Eq. (7)
appears in the Ampère–Maxwell Law:

$\displaystyle\bm{\nabla}\times\bm{B}-\partial_{t}\bm{E}=\bm{J}-g_{a\gamma}{\left(\partial_{t}a\right)}\bm{B}~{}.$
(8)

It follows that the physical effects of the axion–photon coupling in the
presence of a magnetic field $\bm{B}$, as in the case of hidden photons [Eq.
(5)], manifest as an effective current:

$\displaystyle\bm{J}_{a}=-g_{a\gamma}{\left(\partial_{t}a\right)}\bm{B}=ig_{a\gamma}m_{a}a{\left(\bm{r},t\right)}\bm{B}~{},$
(9)

where

$\displaystyle
a{\left(\bm{r},t\right)}=a_{0}{\left(\bm{r},t\right)}e^{-im_{a}t}$ (10)

is the axion field with a stochastically (slowly) varying amplitude
$\left|a_{0}\right|\sim\sqrt{2\rho{{}_{\mbox{\scriptsize DM}}}}/m_{a}$, with
coherence length $\ell{{}_{\mbox{\scriptsize coh}}}$ and coherence time
$\tau{{}_{\mbox{\scriptsize coh}}}$ analogous to those for hidden photons
described by Eqs. (3) and (4), with the replacement $m_{A^{\prime}}\rightarrow
m_{a}$. The interaction of an axion dark-matter field with the geomagnetic
field of Earth thus generates an oscillating magnetic-field pattern, which is
discussed in detail in Ref. [28].

In this work, we aim to analyze the first dedicated measurements of the SNIPE
Hunt experiment in the frequency range 0.5–5 Hz. The lower frequency bound of
0.5 Hz for our analysis was chosen for practical reasons: $1/f$ noise begins
to reduce our sensitivity below $\approx 0.5~{}{\rm Hz}$ and there is ongoing
analysis of SuperMAG data covering frequencies up to $\approx 1~{}{\rm Hz}$
that is expected to surpass the sensitivity of this experiment. For the upper
bound of 5 Hz, we have considered the fact that we do not have a robust
prediction for Schumann resonances because of finite conductivity effects and
also inhomogeneities in the ionosphere refractive index [61]. Indeed, the
first Schumann resonance occurs at a frequency around 7.8 Hz with time-
dependent fluctuations of the order of 0.5 Hz. Most importantly, its width is
about 2 Hz, which makes $f\leq 5~{}\text{Hz}$ a region where the dark-matter-
induced magnetic-field pattern can be reliably derived (see Sec. IV.3.1 for
further discussion). The analyses carried out in Refs. [26, 28] considered a
quasi-static limit valid only when the UBDM Compton wavelengths are much
larger than Earth’s radius $R$: $\lambda_{A^{\prime}}\approx
1/m_{A^{\prime}}\gg R$ and $\lambda_{a}\approx 1/m_{a}\gg R$. This sets an
upper limit on the hidden-photon mass $m_{A^{\prime}}$ and axion mass $m_{a}$
of $\sim 3\times 10^{-14}~{}{\rm eV}$ and, correspondingly, for their Compton
frequencies: $f_{A^{\prime}}$ and $f_{a}$ must be $\ll 7~{}{\rm Hz}$. As we
are working at frequencies up to 5 Hz, the formulas used in Refs. [26, 28] are
only marginally correct, and therefore more robust formulas are needed here.

In the following we calculate a more general signal for dark-matter masses
close to $\sim 1/R$. We write the magnetic and electric fields in terms of
vector spherical harmonics (VSH; see Appendix D of [26]) $\bm{Y}_{\ell m}$,
$\bm{\Psi}_{\ell m}$, $\bm{\Phi}_{\ell m}$ as

$\displaystyle\bm{B}(\bm{x},t)$ $\displaystyle=e^{-i\omega
t}\sum_{\ell,m}\left(B^{(r)}_{\ell m}(r)\bm{Y}_{\ell m}+B^{(1)}_{\ell
m}(r)\bm{\Psi}_{\ell m}+B^{(2)}_{\ell m}(r)\bm{\Phi}_{\ell m}\right)$ (11)
$\displaystyle\bm{E}(\bm{x},t)$ $\displaystyle=e^{-i\omega
t}\sum_{\ell,m}\left(E^{(r)}_{\ell m}(r)\bm{Y}_{\ell m}+E^{(1)}_{\ell
m}(r)\bm{\Psi}_{\ell m}+E^{(2)}_{\ell m}(r)\bm{\Phi}_{\ell m}\right),$ (12)

where $\omega$ is the oscillation angular frequency of the dark-matter
effective current. For the dark-matter effective current $\bm{J}$ which stands
for both hidden photons and axion-like particles, we use the fact that it
satisfies $\bm{\nabla}\times\bm{J}=0$ to write

$\displaystyle\bm{J}(\bm{x},t)$ $\displaystyle=e^{-i\omega
t}\sum_{\ell,m}\left(J^{(r)}_{\ell m}(r)\bm{Y}_{\ell m}+J^{(1)}_{\ell
m}(r)\bm{\Psi}_{\ell m}\right).$ (13)

Inserting the above ansatz into Maxwell’s equations, we get

$\displaystyle\left({1\over r^{2}}{d\over dr}\left(r^{2}{d\over
dr}\right)+\omega^{2}-{\ell(\ell+1)\over
r^{2}}\right)\left(\begin{array}[]{cc}B^{(2)}_{\ell m}\\\ E^{(2)}_{\ell
m}\end{array}\right)=0~{},$ (16)

and the other components are determined by

$\displaystyle E^{(r)}_{\ell m}$ $\displaystyle={1\over
i\omega}\left({\ell(\ell+1)\over r}B^{(2)}_{\ell m}+J^{(r)}_{\ell m}\right)$
(17) $\displaystyle E^{(1)}_{\ell m}$ $\displaystyle={1\over
i\omega}\left({1\over r}{d\over dr}\left(rB^{(2)}_{\ell
m}\right)+J^{(1)}_{\ell m}\right)$ (18) $\displaystyle B^{(r)}_{\ell m}$
$\displaystyle=-{1\over i\omega}{\ell(\ell+1)\over r}E^{(2)}_{\ell m}$ (19)
$\displaystyle B^{(1)}_{\ell m}$ $\displaystyle=-{1\over i\omega}{1\over
r}{d\over dr}\left(rE^{(2)}_{\ell m}\right).$ (20)

This system is solved with boundary conditions such that $E^{(1)}_{\ell m}$
and $E^{(2)}_{\ell m}$ vanish at both Earth’s surface $r=R$ and ionosphere
$r=R+h$, where $h$ is the ionosphere height. Because we work in the regime
$\omega h\ll 1$, the boundary condition for $E^{(2)}_{\ell m}$ implies
immediately that it is zero everywhere; it follows that $B^{(r)}_{\ell m}$ and
$B^{(1)}_{\ell m}$ also vanish identically.

Writing $B^{(2)}_{\ell m}=u_{\ell m}/r$, in the limit in which $h\ll R$ we
find

$\displaystyle u_{\ell m}^{\prime\prime}-\lambda_{\ell}^{2}u_{\ell m}=0,$ (21)

where $\lambda_{\ell}^{2}=\ell(\ell+1)/R^{2}-\omega^{2}$. We write the
solution for $u_{\ell m}$ as $u_{\ell m}=\alpha_{\ell
m}\cosh(\lambda_{\ell}(r-R))+\beta_{\ell m}\sinh(\lambda_{\ell}(r-R))$. Notice
that the magnetic field signal at Earth’s surface ($r=R$) is simply given by

$\displaystyle\bm{B}=\sum_{\ell,m}{\alpha_{\ell m}\over R}\bm{\Phi}_{\ell m}.$
(22)

From the boundary condition $u_{\ell m}^{\prime}=-rJ^{(1)}_{\ell m}$ at $r=R$
and $r=R+h$, we find at zeroth order in $h/R$

$\displaystyle\alpha_{\ell m}=-{J^{(1)}_{\ell m}(R)+RJ^{(1)}_{\ell
m}{}^{\prime}(R)\over\lambda_{\ell}^{2}}.$ (23)

###  II.1 Hidden-Photon Signal

In terms of vector spherical harmonics, the hidden-photon effective current,
given in Eq. (5), is written as

$\displaystyle\bm{J}_{A^{\prime}}=-\sqrt{4\pi\over 3}\varepsilon
m_{A^{\prime}}^{2}\sum_{m=-1}^{1}A_{m}(\bm{Y}_{1m}+\bm{\Psi}_{1m})e^{-i\omega_{m}t}\enspace.$
(24)

Here $\omega_{m}=m_{A^{\prime}}-2\pi f_{d}m$, where $f_{d}$ is the frequency
associated to the sidereal day,555The appearance of $f_{d}$ here is due to the
rotation of Earth. While the direction of the hidden photon is fixed in the
inertial celestial frame, our measurements are performed by magnetometers
which are fixed to the rotating Earth. Transforming the hidden-photon
amplitude from the inertial to co-rotating frame, introduces an additional
time dependence related to Earth’s rotational frequency. and the hidden-photon
amplitudes $A^{\prime}_{m}$ (for polarizations $m=0,\pm 1$) appearing in Eq.
(24) are normalized via

$\frac{1}{2}m_{A^{\prime}}^{2}\langle|A^{{}^{\prime}}|^{2}\rangle=\rho_{\mathrm{DM}},$
(25)

where $\rho_{\mathrm{DM}}=0.3\,\mathrm{GeV/cm}^{3}$ is the local dark-matter
density. Extracting $J^{(1)}_{1m}$ from Eq. (24), we find

$\displaystyle\bm{B}_{A^{\prime}}=\sqrt{4\pi\over
3}{\varepsilon\,m_{A^{\prime}}^{2}R\over
2-m_{A^{\prime}}^{2}R^{2}}\sum_{m=-1}^{1}A_{m}\bm{\Phi}_{1m}e^{-i\omega_{m}t}.$
(26)

###  II.2 Axion Signal

For axion dark matter, the orientation of the effective current is determined
by Earth’s dc magnetic field [see Eq. (9)]. As in Ref. [28], we utilize the
IGRF-13 model [62], which parameterizes Earth’s magnetic field
$\bm{B}_{\oplus}$ in terms of a scalar potential $V_{0}$, such that
$\bm{B}_{\oplus}=-\nabla V_{0}$, where $V_{0}$ is expanded as

$V_{0}=\sum_{\ell=1}^{\infty}\sum_{m=0}^{\ell}\frac{R^{\ell+2}}{r^{\ell+1}}(g_{\ell
m}\cos(m\phi)+h_{\ell m}\sin(m\phi))P_{\ell}^{m}(\cos\theta),$ (27)

where $P_{\ell}^{m}$ are the Schmidt-normalized associated Legendre
polynomials. The Gauss coefficients $g_{\ell m}$ and $h_{\ell m}$ are
specified by the IGRF model at five-year intervals (see Tab. 2 of Ref. [62]).
The last of these coefficients correspond to the year 2020, with time
derivatives provided for their subsequent evolution. In this work, we
extrapolate the 2020 values (up to $\ell=4$) forward to July 23, 2022 using
these time derivatives, and adopt the conventions $g_{\ell,-m}=(-1)^{m}g_{\ell
m}$ to $h_{\ell,-m}=(-1)^{m+1}h_{\ell m}$ to extend to negative $m$.

Once Earth’s dc field has been parametrized in this way, the effective current
that axion dark matter of mass $m_{a}$ and axion–photon coupling $g_{a\gamma}$
generates can be written as [28]

$\bm{J}_{a}=ig_{a\gamma}a_{0}m_{a}\sum_{\ell,m}C_{\ell m}\left(R\over
r\right)^{\ell+2}((\ell+1)\bm{Y}_{\ell m}-\bm{\Psi}_{\ell m})e^{-im_{a}t},$
(28)

where $a_{0}$ is the (complex) axion amplitude, normalized by
$\frac{1}{2}m_{a}^{2}\langle|a_{0}|^{2}\rangle=\rho_{\mathrm{DM}}$, and

$C_{\ell m}=(-1)^{m}\sqrt{\frac{4\pi(2-\delta_{m0})}{2\ell+1}}\frac{g_{\ell
m}-ih_{\ell m}}{2}.$ (29)

Now, by identifying $J^{(1)}(r)$ in Eq. (28), the magnetic-field signal from
axion dark matter is found to be

$\displaystyle\bm{B}_{a}=-iga_{0}m_{a}R\sum_{\ell,m}{(\ell+1)C_{\ell
m}\over\ell(\ell+1)-m_{a}^{2}R^{2}}\bm{\Phi}_{\ell m}e^{-im_{a}t}.$ (30)

##  III Experimental Details

From 21 July 2022 to 24 July 2022, we conducted the first coordinated SNIPE
Hunt science run. Measurements were made with battery-operated magnetometers
located at three sites which were chosen to have minimal magnetic-field
interference from power lines, traffic, and other anthropic sources. A block
diagram of the experimental setup at an individual station is shown in Fig. 1.
The magnetometers were Vector Magnetoresistive (VMR) sensors manufactured by
Twinleaf LLC. The VMRs use three mutually perpendicular giant magnetoresitive
(GMR) field sensors to measure all three components of the magnetic field. The
sensitivity of the GMR sensors is specified to be 300
$\mathrm{pT}/\sqrt{\mathrm{Hz}}$ over a frequency range of 0.1–100 Hz. In
addition to the magnetic field, the VMR also has a three-axis gyroscope, a
three-axis accelerometer, a barometer, and a thermometer. The measurements
from all of these sensors were recorded during the course of the science run
on a laptop computer which also provided power to the VMR via a USB
connection. The sample rate for the data acquisition was set to 160 samples/s.
In order to limit the influence of magnetic noise from the laptop on the VMR,
the laptop was located in a camping tent 9–12 m from the sensor, depending on
the station. The laptops were powered by 50 ${\rm{A\cdot hr}}$ powerbanks,
which were swapped with fully charged powerbanks every 6–10 hours and
recharged using a solar generator. Fig. 3 shows the operation times for the
three stations.

The data were time stamped using the computer clocks, which were steered to
GPS time using a receiver antenna and synchronization software. To account for
the software lag present in the timing calibration, the timing offset
correction was set prior to the science run using a time server from the
National Institute for Standards and Technology. The accuracy of the timing
was tested in the laboratory by applying magnetic-field signals that were
triggered by an external GPS receiver before and after the science run. Based
on these tests, we estimate the accuracy of the timing to be $\lesssim 100$
ms.

The location of the three stations is shown in Table 1. The magnetometers were
aligned so that the $y$ axis of the magnetometers was vertical, relative to
local gravity, and the $z$ axis of the detectors was pointing to true north as
determined by smart-phone compasses. We estimate the pointing accuracy of the
detectors to be $\lesssim 1^{\circ}$. An example of one of the mounts used for
the alignment of the magnetometers is shown in Fig. 2. The sensors and mounts
were covered with a plastic container that was secured to the ground to guard
against rain.

Station | Location | Latitude | Longitude | Elevation
---|---|---|---|---
|  | (deg) | (deg) | (meters)
Hayward | Auburn State Recreation Area | 39.1017 | -120.924 | 355.0
Lewisburg | Penn Roosevelt State Park | 40.7404 | -77.7113 | 692.2
Oberlin | Findley State Park | 41.1303 | -82.2069 | 277.4
Table 1: Locations of sensors used in the 2022 SNIPE hunt. The stations are
referred to by the location of the home institution for the groups in charge
of each station. Figure 1: Block diagram of SNIPE station setup. A three-axis
GMR magnetometer was connected via USB to a laptop located 9–12 m from the
sensor. The data were recorded with a laptop and time stamped using the laptop
computer time, which was steered to GPS time using a GPS timing receiver. The
laptop was powered with battery power banks that were swapped out every 6–10
hours. Figure 2: Mount for the detector. The pitch, roll, and yaw can be
adjusted. A smart phone fits onto the table that holds the sensor for
alignment. The phone is removed during data collection. The mount was attached
to the ground using heavy-duty plastic tent screws. Figure 3: Activity for the
2022 SNIPE science run. The horizontal bars indicate when the Hayward,
Lewisburg, and Oberlin stations were operational. Two subsets of the data were
analyzed independently: Scan-1 covering the interval shown as the light blue
shaded region on the left, and Scan-2, the grey shaded region on the right.
Figure 4: Amplitude spectral densities of the North-South and East-West
components of the magnetic field measurements from the three measurement
sites. The shaded band 0.5–5.0 Hz shows the range of frequencies probed in
this work. In this band, the noise floor is limited by the instrumental
sensitivity of $\sim 300\,\mathrm{pT}/\sqrt{\mathrm{Hz}}$.

###  III.1 Noise Characteristics

For the three sites, we show in Fig. 4 the amplitude spectral density for the
East-West and North-South components of the magnetic field – the components
relevant for this search. A couple of features are evident. The Hayward
station had noticeably smaller power-line noise at 60 Hz than the Lewisburg
and Oberlin stations. The Lewisburg station had a significant $1/f$ pedestal
in the 0.1 to 0.5 Hz band that was absent in the other two stations. Also, the
Oberlin station had narrow peaks at 0.25, 0.5, and 0.75 Hz suggesting a common
origin as harmonics of some fundamental frequency. As the local magnetic
environments are distinct, this difference in noise profile between the
stations is expected even though we have not identified the origins of the
particular features noted above. However, for the three stations, the
amplitude spectral density in most of the band of interest is flat and
corresponds to approximately $300\,\mathrm{pT/\sqrt{Hz}}$, the noise floor of
the sensors.

Figure 5: Time series of magnetic fields made at the Hayward, Lewisburg, and
Oberlin measurement stations. The North-South, East-West, and Vertical (normal
to Earth’s surface) directions are shown. Scan-1 begins at time $t=0$, and
covers the first 24 hours of the data shown. The red dashed boxes correspond
to the occurrence of a geomagnetic storm. During that time, we noticed
correlated low-frequency oscillations in all three stations. Data from this
period were not included in Scan-1, as discussed in the main text. The blue
dashed line shows the sensor temperature measured at the different locations.

In Fig. 5, we plot time series of the sensor temperature (shown as the blue
dashed lines on the right), and of the temperature-corrected measurements of
the magnetic field covering the first $\sim 30$ hours of the observing run.
The rows correspond to the different sites, and the columns to the North-
South, East-West, and Vertical components of the field. We apply the
temperature correction purely for plotting purposes, as we noticed a
temperature-dependent drift in the sensor calibration at dc of up to 10
percent in the case of the Hayward station and about 2 percent for the other
two stations. However, in the analysis band – 0.5 to 5.0 Hz – we do not make
any temperature correction. Instead, as we discuss in Sec. IV.3, we assign an
uncertainty on the quoted HPDM and axion limits due to temperature drifts.

Between hours $\sim 13$ and 20 of the time series, we observe increased
fluctuations in the North and East components of the Lewisburg data –
fluctuations which were not present in the other stations. This interval
coincides with an overnight thunderstorm, during which mechanical agitation of
the sensor or lightning occurring nearby may have led to the fluctuations.
However, in the temporal window between hours $\sim 25$ and $32$ (shown
enclosed in the red dashed boxes of Fig. 5), we notice features which are
clearly correlated across all three stations, and which we believe are due to
a geomagnetic storm associated with the eruption of sunspot AR3060. This
produced a C5-class solar flare and a coronal mass ejection directed toward
Earth [63, 64]. The storm led to the modulation of Earth’s magnetic field
which we detected. Including data from this window in the the analysis
presented below led to noticeable non-gaussianities in the test statistic used
for setting limits on the HPDM and axion parameters. For this reason, we
excluded the time interval containing the geomagnetic storm in the analysis
and instead separate the data into two independently analyzed measurement
periods: Scan-1 and Scan-2. These time periods are shown as shaded regions in
Fig. 3.

##  IV Data Analysis

In this section, we outline how the SNIPE Hunt data is analyzed to search for
both a hidden-photon dark-matter (HPDM) and axion dark-matter signal.

###  IV.1 Hidden-Photon Analysis

We begin with the HPDM signal. Our analysis follows a similar (but simplified)
methodology to that described in Ref. [27]. In this search, our data consist
of six time series, corresponding to the south-directed and east-directed
magnetic field components measured at each of the three SNIPE Hunt measurement
locations: $B_{\theta}(\Omega_{1},t_{j})$, $B_{\phi}(\Omega_{1},t_{j})$,
$B_{\theta}(\Omega_{2},t_{j})$, $B_{\phi}(\Omega_{2},t_{j})$,
$B_{\theta}(\Omega_{3},t_{j})$, and $B_{\phi}(\Omega_{3},t_{j})$.666Here
$\Omega_{i}=(\theta_{i},\phi_{i})$ denotes the geographic coordinates of each
station. Note that while $\phi_{i}$ is exactly the longitude of each station,
the latitude of each station is given by $\frac{\pi}{2}-\theta_{i}$. Likewise,
$\bm{\hat{\phi}}$ points east, while $\bm{\hat{\theta}}$ points south. We
model these time series as being given by (the real part of) the signal in Eq.
(26) plus Gaussian white noise. Our goal is then to extract a bound on
$\varepsilon$. As the exact amplitudes $A^{\prime}_{m}$ are unknown, we
utilize a Bayesian framework and treat these as nuisance parameters. We also
take a Gaussian distribution for them,777$A^{\prime}$ can be written as a sum
of several independent plane wave solutions of different velocities
$v_{n}\sim\mathcal{O}(v_{\mathrm{DM}})$. These have corresponding frequencies
$f_{n}\sim f_{A^{\prime}}\left(1+\mathcal{O}(v_{\mathrm{DM}}^{2})\right)$. On
timescales longer than $\tau_{\mathrm{coh}}\sim
2\pi/(f_{A^{\prime}}v_{\mathrm{DM}}^{2})$, the value of $A^{\prime}$ is thus a
sum of many contributions with random phases. By the central limit theorem, it
is thus distributed as a Gaussian variable. normalized by Eq. (25).

The signal in Eq. (26) indicates that all relevant information is contained at
the frequencies $f_{A^{\prime}}$ and $f_{A^{\prime}}\pm f_{d}$. Thus we
Fourier transform the six time series $B_{\alpha}(\Omega_{i})$, and construct
an 18-dimensional data vector888We use $\vec{x}$ to denote a vector $x$ with
18 components (or six components in Sec. IV.2), and $\bm{y}$ to indicate a
vector $y$ with three components. $\vec{X}$ which contains all information
which may be relevant to setting a bound at $f_{A^{\prime}}$. Namely,
$\vec{X}$ consists of the six values
$\tilde{B}_{\alpha}\left(\Omega_{i},f_{A^{\prime}}-\hat{f}_{d}\right)$,
followed by the six values
$\tilde{B}_{\alpha}\left(\Omega_{i},f_{A^{\prime}}\right)$, followed by the
six values
$\tilde{B}_{\alpha}\left(\Omega_{i},f_{A^{\prime}}+\hat{f}_{d}\right)$. In our
analysis, we compute bounds only at discrete Fourier transform (DFT)
frequencies $f_{A^{\prime}}=n/T$ (where $T$ is the total duration of the time
window in consideration). Note that $f_{d}$ may not generically be a DFT
frequency, and so we have instead used $\hat{f}_{d}$, which we define as the
nearest DFT frequency to $f_{d}$. With these choices, $\vec{X}$ can be
computed via a fast Fourier transform (FFT). (This allows us to compute
$\vec{X}$ at all frequencies simultaneously, and perform the subsequent
analysis for all frequencies in parallel.) The first step of our analysis is
to characterize the statistics of $\vec{X}$, namely its expectation and
variance.

First, let us compute the expectation of $\vec{X}$. As mentioned above, we
model our measurements as being Gaussian noise on top of the signal in Eq.
(26). Since the expectation of the noise vanishes, the expectation of
$\vec{X}$ simply comes from Fourier transforming Eq. (26) and assembling its
relevant components into a vector. To remove the normalization from the
amplitudes $A^{\prime}_{m}$, let us define

$c_{m}=\frac{\sqrt{2}\pi
f_{A^{\prime}}A^{\prime}_{m}}{\sqrt{\rho_{\mathrm{DM}}}}.$ (31)

These now have $\sum_{m}\langle|c_{m}|^{2}\rangle=1$. In the case $c_{\pm}=0$,
(the real part of) Eq. (26) takes the simple form

$\bm{B}_{0}(\Omega,t)=-\frac{2\pi f_{A^{\prime}}R}{2-(2\pi
f_{A^{\prime}}R)^{2}}\varepsilon\sqrt{2\rho_{\mathrm{DM}}}\sin\theta\cdot\text{Re}\left[c_{0}e^{-2\pi
if_{A^{\prime}}t}\right]\bm{\hat{\phi}},$ (32)

and the only nonzero components of $\langle\vec{X}\rangle$ are

$\displaystyle\langle X_{8}\rangle_{0}$
$\displaystyle=\tilde{B}_{0,\phi}(\Omega_{1},f_{A^{\prime}})=-\frac{2\pi
f_{A^{\prime}}R}{2-(2\pi f_{A^{\prime}}R)^{2}}c_{0}^{*}\varepsilon
T\sqrt{\frac{\rho_{\mathrm{DM}}}{2}}\sin\theta_{1}\equiv
c_{0}^{*}\varepsilon\mu_{0,8}$ (33) $\displaystyle\langle X_{10}\rangle_{0}$
$\displaystyle=\tilde{B}_{0,\phi}(\Omega_{2},f_{A^{\prime}})=-\frac{2\pi
f_{A^{\prime}}R}{2-(2\pi f_{A^{\prime}}R)^{2}}c_{0}^{*}\varepsilon
T\sqrt{\frac{\rho_{\mathrm{DM}}}{2}}\sin\theta_{2}\equiv
c_{0}^{*}\varepsilon\mu_{0,10}$ (34) $\displaystyle\langle X_{12}\rangle_{0}$
$\displaystyle=\tilde{B}_{0,\phi}(\Omega_{3},f_{A^{\prime}})=-\frac{2\pi
f_{A^{\prime}}R}{2-(2\pi f_{A^{\prime}}R)^{2}}c_{0}^{*}\varepsilon
T\sqrt{\frac{\rho_{\mathrm{DM}}}{2}}\sin\theta_{3}\equiv
c_{0}^{*}\varepsilon\mu_{0,12},$ (35)

On the other hand if $c_{0}=c_{-}=0$, then the signal becomes

$\bm{B}_{+}(\Omega,t)=\frac{2\pi f_{A^{\prime}}R}{2-(2\pi
f_{A^{\prime}}R)^{2}}\varepsilon\sqrt{\rho_{\mathrm{DM}}}\cdot\text{Re}\left[c_{+}\left(i\bm{\hat{\theta}}-\cos\theta\bm{\hat{\phi}}\right)e^{-2\pi
i(f_{A^{\prime}}-f_{d})t+i\phi}\right],$ (36)

and so the expectation of $\vec{X}$ is

$\langle\vec{X}\rangle_{+}\approx-\frac{\pi f_{A^{\prime}}R}{2-(2\pi
f_{A^{\prime}}R)^{2}}c_{+}^{*}\varepsilon\Delta
t\sqrt{\rho_{\mathrm{DM}}}\begin{pmatrix}ie^{-i\phi_{1}}Q(f_{d}-\hat{f}_{d})\\\
\cos\theta_{1}e^{-i\phi_{1}}Q(f_{d}-\hat{f}_{d})\\\
ie^{-i\phi_{2}}Q(f_{d}-\hat{f}_{d})\\\
\cos\theta_{2}e^{-i\phi_{2}}Q(f_{d}-\hat{f}_{d})\\\
ie^{-i\phi_{3}}Q(f_{d}-\hat{f}_{d})\\\
\cos\theta_{3}e^{-i\phi_{3}}Q(f_{d}-\hat{f}_{d})\\\ ie^{-i\phi_{1}}Q(f_{d})\\\
\cos\theta_{1}e^{-i\phi_{1}}Q(f_{d})\\\ ie^{-i\phi_{2}}Q(f_{d})\\\
\cos\theta_{2}e^{-i\phi_{2}}Q(f_{d})\\\ ie^{-i\phi_{3}}Q(f_{d})\\\
\cos\theta_{3}e^{-i\phi_{3}}Q(f_{d})\\\ ie^{-i\phi_{1}}Q(f_{d}+\hat{f}_{d})\\\
\cos\theta_{1}e^{-i\phi_{1}}Q(f_{d}+\hat{f}_{d})\\\
ie^{-i\phi_{2}}Q(f_{d}+\hat{f}_{d})\\\
\cos\theta_{2}e^{-i\phi_{2}}Q(f_{d}+\hat{f}_{d})\\\
ie^{-i\phi_{3}}Q(f_{d}+\hat{f}_{d})\\\
\cos\theta_{3}e^{-i\phi_{3}}Q(f_{d}+\hat{f}_{d})\end{pmatrix}\equiv
c_{+}^{*}\varepsilon\vec{\mu}_{+},$ (37)

where

$Q(f)=\frac{1-e^{-2\pi ifT}}{1-e^{-2\pi if\Delta t}},$ (38)

and $\Delta t=(1/160)\,\mathrm{s}$ is the time resolution. Note that, in
principle, Eq. (37) should have an additional term proportional to $c_{+}$,
which contains factors of $Q(2f_{A^{\prime}}-f_{d}-\hat{f}_{d})$,
$Q(2f_{A^{\prime}}-f_{d})$, and $Q(2f_{A^{\prime}}-f_{d}+\hat{f}_{d})$. Since
$f_{d}\ll f_{A^{\prime}}$ and $Q(f)\sim 1/f$, these will all be significantly
smaller than the $Q$ factors appearing in Eq. (37). Thus we are safe to
neglect this additional term. Similarly, $\langle\vec{X}\rangle_{-}\equiv
c_{-}^{*}\varepsilon\vec{\mu}_{-}$ can be computed (for the case when
$c_{0}=c_{+}=0$). Then generically, the full expectation of $\vec{X}$ is

$\langle\vec{X}\rangle=\varepsilon(c_{+}^{*}\vec{\mu}_{+}+c_{0}^{*}\vec{\mu}_{0}+c_{-}^{*}\vec{\mu}_{-}).$
(39)

Now that we have computed the expectation of $\vec{X}$, let us consider its
variance. In this analysis, we consider the frequency range
$0.5\,\text{Hz}\leq f_{A^{\prime}}\leq 5\,\text{Hz}$, over which the noise is
roughly frequency independent [see Fig. (4)]. Therefore, we may consider each
instance of $\vec{X}$ for different frequencies as independent realizations of
the noise, and use these to estimate the noise. In particular, we can compute
the covariance matrix for $\vec{X}$ as

$\Sigma_{ij}\equiv\langle
X_{i}X_{j}^{*}\rangle=\frac{1}{N}\sum_{k=1}^{N}X_{i}(f_{k})X_{j}(f_{k})^{*},$
(40)

where $f_{k}$ indexes the DFT frequencies between $0.5\,\text{Hz}$ and
$5\,\text{Hz}$ (for $k=1,\ldots,N\sim 10^{5}$).999Note that since the first
six elements, the middle six elements, and the final six elements of $\vec{X}$
correspond to different frequencies, then covariances between elements from
these different groups should vanish, i.e. $\Sigma$ should be block diagonal.
Moreover, the three diagonal blocks should be identical, since they correspond
to the same averages in Eq. (40) (only with the frequency $f_{k}$ shifted by
$\hat{f}_{d}$). Thus it suffices to only compute $\Sigma_{ij}$ for $7\leq
i,j\leq 12$.

Now that we understand the statistics of $\vec{X}$, we can write down its
likelihood

$-\ln\mathcal{L}\left(\varepsilon,\bm{c}|\vec{X}\right)=\left(\vec{X}-\varepsilon\sum_{m}c_{m}^{*}\vec{\mu}_{m}\right)^{\dagger}\Sigma^{-1}\left(\vec{X}-\varepsilon\sum_{m}c_{m}^{*}\vec{\mu}_{m}\right).$
(41)

From this likelihood, the computation of the bound on $\varepsilon$ proceeds
as in Sec. V D of Ref. [27], but we reproduce it here for completeness. Let us
write $\Sigma=LL^{\dagger}$ and then define

$\displaystyle\vec{Y}$ $\displaystyle=L^{-1}\vec{X},$ (42)
$\displaystyle\vec{\nu}_{m}$ $\displaystyle=L^{-1}\vec{\mu}_{m}.$ (43)

If we let $N$ be the $18\times 3$ matrix whose columns are $\vec{\nu}_{m}$,
then Eq. (41) becomes

$-\ln\mathcal{L}\left(\varepsilon,\bm{c}|\vec{Y}\right)=\left|\vec{Y}-\varepsilon
N\bm{c}^{*}\right|^{2}.$ (44)

Now if we perform a singular value decomposition $N=USV^{\dagger}$ (where $U$
is a $18\times 3$ matrix with orthonormal columns, $S$ is a $3\times 3$
diagonal matrix, and $V$ is a $3\times 3$ unitary matrix) and further define

$\displaystyle\bm{d}$ $\displaystyle=V^{\dagger}\bm{c}^{*},$ (45)
$\displaystyle\bm{Z}$ $\displaystyle=U^{\dagger}\vec{Y},$ (46)

then the likelihood in Eq. (44) can be reduced to

$-\ln\mathcal{L}\left(\varepsilon,\bm{d}|\bm{Z}\right)=\left|\bm{Z}-\varepsilon
S\bm{d}\right|^{2}.$ (47)

As mentioned earlier, the polarization amplitudes $c_{m}$, and thus also the
parameters $d_{m}$, are nuisance parameters over which we need to marginalize.
We take them to have a Gaussian likelihood

$\mathcal{L}(\bm{d})=\exp(-3|\bm{d}|^{2}).$ (48)

Marginalizing over $\bm{d}$, the likelihood Eq. (47) reduces to

$\mathcal{L}\left(\varepsilon|\bm{Z}\right)\propto\prod_{m}\frac{1}{3+\varepsilon^{2}s_{m}^{2}}\exp\left(-\frac{3|z_{m}|^{2}}{3+\varepsilon^{2}s_{m}^{2}}\right),$
(49)

where $z_{m}$ are the components of $\bm{Z}$ and $s_{m}$ are the diagonal
entries of $S$ [see Appendix D 1 of Ref. [27] for a derivation of Eq. (49)].
In order to turn this into a posterior on $\varepsilon$, we must assume some
prior. We take a Jeffreys prior

$p(\varepsilon)\propto\sqrt{\sum_{m}\frac{4\varepsilon^{2}s_{m}^{4}}{(3+\varepsilon^{2}s_{m}^{2})^{2}}};$
(50)

again see Appendix D 1 of Ref. [27]. The posterior for $\varepsilon$ is thus

$p(\varepsilon|\bm{Z})=\mathcal{N}\sqrt{\sum_{m}\frac{4\varepsilon^{2}s_{m}^{4}}{(3+\varepsilon^{2}s_{m}^{2})^{2}}}\prod_{m}\frac{1}{3+\varepsilon^{2}s_{m}^{2}}\exp\left(-\frac{3|z_{m}|^{2}}{3+\varepsilon^{2}s_{m}^{2}}\right),$
(51)

where $\mathcal{N}$ must be calculated to normalize the integral of
$p(\varepsilon|\bm{Z})$ to 1. We then set a 95% credible upper limit
$\hat{\varepsilon}$ by solving

$\int_{0}^{\hat{\varepsilon}}d\varepsilon\,p(\varepsilon|\bm{Z})=0.95.$ (52)

By performing this analysis at all DFT frequencies between 0.5 Hz and 5 Hz, we
arrive at a bound over a range of HPDM masses. Fig. (6) shows the results of
our analysis for both Scan-1 and Scan-2.

Figure 6: 95% credible upper limit on $\varepsilon$, the HPDM kinetic-mixing
parameter. The top figure shows the results for Scan-1, and the bottom figure
shows the results for Scan-2. The orange traces on both plots are smoothed
versions of the limits obtained by averaging over 100 adjacent frequency bins.

Following the methodology in Sec. VI of Ref. [27], we evaluate our data at
each frequency for evidence of a significant dark-matter candidate. From Eq.
(47), we see that under the null hypothesis of no dark matter signal
($\varepsilon=0$), the vector $\bm{Z}$ should be distributed as a multivariate
Gaussian of mean zero. Specifically, the statistic

$Q=2\sum_{m}|z_{m}|^{2}$ (53)

should follow a $\chi^{2}$-squared distribution with six degrees of freedom.
We may therefore compute the corresponding local $p$-value

$p_{0}=1-F_{\chi^{2}(6)}(Q),$ (54)

where $F_{\chi^{2}(\nu)}$ denotes the cumulative distribution function for a
$\chi^{2}$-distribution with $\nu$ degrees of freedom. Fig. (7) shows the
local $p$-values at each frequency $f_{A^{\prime}}$ for both Scan-1 and
Scan-2. We consider there to be evidence for a DM candidate at a given
frequency (with 95% global significance) if its local $p$-value is below the
threshold $p_{\text{crit}}$ defined by

$(1-p_{\text{crit}})^{N}=0.95.$ (55)

This threshold is shown as a dotted line in Fig. (7). Scan-1 exhibits seven
frequency bins which cross the threshold. Four of these are clustered around
0.5 Hz, while the other three are clustered around 0.75 Hz. Scan-2, likewise,
exhibits three candidate frequencies clustered around 0.5 Hz, and one at 0.75
Hz. We expect these candidates are associated with the narrow peaks observed
in the Oberlin station data. We have re-performed our analysis using only the
Hayward and Lewisburg data, and find that these peaks do not cross the
threshold for significance in either scan when restricting to these two
stations [see Fig. (8)]. Since dark matter should be present in all locations
at all times, this strongly suggests that these signal candidates do not
correspond to dark matter. Moreover, we note that the width of a dark-matter
signal is given by $f_{a}v_{\mathrm{DM}}^{2}$, where $v_{\mathrm{DM}}$ is the
dark matter velocity dispersion. Since the frequency bin size for our analysis
is roughly $10^{-5}\,\mathrm{Hz}$, these signal candidates have widths of
roughly $10^{-5}f_{a}$, corresponding to a large velocity dispersion of
$v_{\mathrm{DM}}\sim 1000\,\mathrm{km/s}$ (which is far above the escape
velocity of the Milky Way). We therefore rule out these dark-matter candidates
and conclude that our analysis finds no evidence for HPDM in the
$0.5\,\text{Hz}\leq f_{A^{\prime}}\leq 5\,\text{Hz}$ range.

Figure 7: The local $p_{0}$-values for each of the $N=414572$ frequency bins
analyzed in the Scan-1, shown in the top (blue) figure, and each of the
$N=340291$ bins searched in Scan-2, shown in the lower (grey) figure. The
threshold value for declaring a dark-matter candidate at 95% global confidence
is shown by the dotted line (after accounting for the trials factor given by
the multiplicity of frequencies searched; see Eq. 55). The left panels show
$p_{0}$ as a function of frequency with candidates having $p$-values below the
threshold. The right panels show histograms of $p_{0}$ for the two different
scans and candidates as outliers to the right of the threshold.  Figure 8: The
local $p_{0}$-values for each frequency bin when only data from the Hayward
and Lewisburg stations are considered. No beyond-threshold candidates appear
in common in _both_ Scan-1 and Scan-2. Also, the peaks at 0.50 and 0.75 Hz
evident in Fig. (7) are not present in this subset of stations. This indicates
that those candidates were due to artefacts in the Oberlin data.

###  IV.2 Axion Analysis

Now we move to the analysis for an axion dark-matter signal. This analysis
proceeds similarly to the HPDM analysis, but is slightly simpler. As in the
HPDM analysis, we construct a data vector $\vec{X}$ consisting of Fourier
transforms of the measured magnetic field at each location. Since the axion
signal in Eq. (30) contains no $f_{d}$ dependence, however, the only relevant
information is contained at frequency $f_{a}$. Therefore in this analysis, we
only take $\vec{X}$ to be a six-dimensional vector, consisting of the
measurements: $\tilde{B}_{\theta}(\Omega_{1},f_{a})$,
$\tilde{B}_{\phi}(\Omega_{1},f_{a})$, $\tilde{B}_{\theta}(\Omega_{2},f_{a})$,
$\tilde{B}_{\phi}(\Omega_{2},f_{a})$, $\tilde{B}_{\theta}(\Omega_{3},f_{a})$,
and $\tilde{B}_{\phi}(\Omega_{3},f_{a})$. The expectation of $\vec{X}$ is now
given by

$\langle\vec{X}\rangle=ic^{*}g_{a\gamma}RT\sqrt{\frac{\rho_{\mathrm{DM}}}{2}}\sum_{\ell
m}\frac{(\ell+1)C_{\ell m}}{\ell(\ell+1)-(2\pi
f_{a}R)^{2}}\begin{pmatrix}\Phi^{\theta}_{\ell m}(\Omega_{1})\\\
\Phi^{\phi}_{\ell m}(\Omega_{1})\\\ \Phi^{\theta}_{\ell m}(\Omega_{2})\\\
\Phi^{\phi}_{\ell m}(\Omega_{2})\\\ \Phi^{\theta}_{\ell m}(\Omega_{3})\\\
\Phi^{\phi}_{\ell m}(\Omega_{3})\end{pmatrix}\equiv
c^{*}g_{a\gamma}\vec{\mu},$ (56)

where $\Phi^{\theta}_{\ell m}$ and $\Phi^{\phi}_{\ell m}$ denote the
$\bm{\hat{\theta}}$-component and $\bm{\hat{\phi}}$-components of the VSH
$\bm{\Phi}_{\ell m}$, and

$c=\frac{\sqrt{2}\pi f_{a}a_{0}}{\sqrt{\rho_{\mathrm{DM}}}}.$ (57)

The covariance matrix $\Sigma$ of $\vec{X}$ can again be determined by
averaging over independent frequencies, as in Eq. (39) [except that $\Sigma$
will now be a $6\times 6$ matrix]. If we define $\vec{Y}$ and $\vec{\nu}$ as
in Eqs. (42) and (43) [without the $m$ index], and further define

$\displaystyle s$ $\displaystyle=|\vec{\nu}|,$ (58) $\displaystyle z$
$\displaystyle=\frac{\vec{\nu}^{\dagger}\vec{Y}}{s},$ (59)

we can write the likelihood function for the axion signal as

$-\ln\mathcal{L}(g_{a\gamma},c|z)=\left|z-g_{a\gamma}c^{*}s\right|^{2}.$ (60)

Again marginalizing over $c$ (which we take to have a Gaussian distribution
with $\langle|c|^{2}\rangle=1$), and utilizing a Jeffreys prior for
$g_{a\gamma}$, we arrive at the posterior distribution

$p(g_{a\gamma}|z)=\frac{|z|^{2}}{1-e^{-|z|^{2}}}\cdot\frac{2g_{a\gamma}s^{2}}{(1+g_{a\gamma}^{2}s^{2})^{2}}\exp\left(-\frac{|z|^{2}}{1+g_{a\gamma}^{2}s^{2}}\right).$
(61)

Note that Eq. (61) is properly normalized, which is possible because its
integral over $g_{a\gamma\gamma}$ can be taken analytically. The 95% credible
limit $\hat{g}_{a\gamma}$ can then be defined, as in Eq. (52). In this case,
we can solve for it analytically to find

$\hat{g}_{a\gamma}=\frac{1}{s}\sqrt{-\frac{|z|^{2}}{\log\left(0.95+0.05e^{-|z|^{2}}\right)}-1}.$
(62)

Fig. (9) shows the resulting limit as a function of frequency, for both Scan-1
and Scan-2. Note that the lower edge of the limit appears as a smooth curve.
This is due to the fact that $\hat{g}_{a\gamma}\rightarrow 4.36/s$ in the
limit $z\rightarrow 0$. Therefore, even when the measured data at a particular
frequency becomes arbitrarily small (compared to the estimated noise level),
the limit on $g_{a\gamma}$ asymptotes to a finite floor.101010This floor
exhibits a slight frequency dependence because of the $f_{a}$-dependence in
Eq. (56).

Figure 9: 95% CL upper limit on $g_{a\gamma\gamma}$ for Scan-1 and Scan-2. The
orange traces on both plots show smoothed versions of the limits obtained by
averaging over 100 adjacent frequency bins.

As in the HPDM case, we evaluate our data at each frequency in order to
determine whether there is evidence for a significant DM signal. We may
compute the local $p$-value at a particular frequency under the null
hypothesis ($g_{a\gamma\gamma=0}$) as

$p_{0}=1-F_{\chi^{2}(2)}(2|z|^{2}).$ (63)

(The $\chi^{2}$-distribution only has two degrees of freedom now, since the
likelihood in Eq. (60) only has one $z$ variable.) Fig. (10) shows these
$p$-values as a function of frequency for both Scan-1 and Scan-2, along with
the threshold value $p_{\text{crit}}$, as defined in Eq. (55). Neither scan
shows any significant signal candidates, and so we again conclude that our
data contains no evidence for axion dark matter in the $0.5\,\text{Hz}\leq
f_{a}\leq 5\,\text{Hz}$ range.

Figure 10: The local $p_{0}$-values for each of the $N=414572$ frequency bins
analyzed in Scan-1 (top), and each of the $N=340291$ frequency bins searched
in Scan-2 (bottom). $p_{\text{crit}}$, the threshold value for declaring a
candidate signal at $95\%$ confidence is shown as the dotted line on each of
the plots. The right panel shows a histogram of all the $p_{0}$-values for
each scan. Signal candidates would appear as outliers to the right of the
threshold.

###  IV.3 Error Budget

The results of this science run and analysis are summarized in Figs. 6 and 9.
They show upper limits on $\varepsilon$, the HPDM kinetic mixing parameter,
and on $g_{a\gamma}$, the axion–photon coupling constant, respectively. Below,
we discuss the impact of uncertainties in the signal model and experimental
conditions on the quoted limits.

####  IV.3.1 Signal model uncertainty

The signals in Eqs. (26) and (30) assume a simplified model of Earth and the
ionosphere, where both are treated as spherical perfect conductors. In Ref.
[26], it is argued that this model holds to a high degree of accuracy in the
frequency range relevant to this work. In particular, both Earth’s crust [65]
and the ionosphere [66, 67] achieve conductivities of at least $10^{-4}$ S/m
at certain depths/heights, which translate to skin depths of $\sim 50\,$km for
frequencies $f\sim 1\,$Hz. Given that the only relevant length scale appearing
in Eqs. (26) and (30) is the radius of Earth $R\sim 6000\,$km, finite-
conductivity effects only modify the geometry of the system at the percent
level. In the absence of resonances, we conclude that the signal should also
only be affected at the percent level.

Close examination of Eqs. (26) and (30), however, reveals that our model
predicts resonances in the signal at $mR=\sqrt{\ell(\ell+1)}$ (for $\ell=1$ in
the HPDM case, and $\ell\geq 1$ in the axion case). These are the well-studied
Schumann resonances of the Earth-ionosphere cavity [61, 68]. Our simplified
spherical model predicts the first of these resonances to occur at $\sim 10$
Hz, but the central frequency of this resonance has been measured to be $\sim
8$ Hz [61], indicating that our spherical model does not accurately account
for environmental effects on the Schumann resonances. Moreover, since the
signal nominally diverges at the Schumann resonances, small deviations in
their central frequency can have a large impact on the predicted signal. For
this reason, we limit our analysis to $f\leq 5$ Hz, in order to remain below
the measured Schumann resonances.

We note that the measured width of the Schumann resonances can, however, be
quite large at certain times. In the summer, during the day, the first
Schumann resonance can reach widths as large as $\sim 4$ Hz [68]. The upper
end of our frequency range may therefore be mildly affected by the first
Schumann resonance for certain portions of the runtime. Such an effect would
result in a slight _enhancement_ of the signal, beyond what our model
predicted. Therefore our exclusion limits are still conservative. In
principle, the effect of the Schumann resonances may, however, invalidate our
signal-candidate rejection procedure. This is because environmental effects
could influence each station differently, meaning we cannot accurately
characterize the spatial dependence of a true signal. To this point, we simply
note that our only signal candidates presented at the end of Sec. IV.1 were at
$f\sim 0.5,0.75\,\mathrm{Hz}$, and so are too low frequency to be affected by
the Schumann resonances. We therefore conclude that both our exclusion
analysis and our candidate rejection are robust to signal-model uncertainties.

####  IV.3.2 Sensor orientation

As discussed in section III, we orient the magnetometers at each site such
that the N-S, and E-W axes of each sensor lie in a horizontal plane with North
indicating True (i.e., geographic) north, and the Normal (Up-Down) axis lies
in the direction of the local force of gravity. We are able to achieve this
orientation with repeatability $\lesssim 1^{\circ}$. By adjusting the
orientation of the sensor in the analysis, we estimate that the impact of such
an orientation error is to change the $\varepsilon$ and $g_{a\gamma}$ upper
limits by $\lesssim 1\%$.

####  IV.3.3 Calibration drift

A temperature-dependent sensor calibration will lead to systematic errors in
magnetic-field measurements. As shown in Fig. 5, we observed that the
temperature swing over the course of a day at the Hayward station was
significantly greater than that in the Oberlin and Lewisburg stations. In that
period, we recorded changes in the dc magnetic-field readings that tracked the
sensor temperature of up to $10\%$ for the Hayward station, and less than
$3\%$ for the Oberlin and Lewisburg stations. In the the 0.5–5.0 Hz band, we
estimate the impact of a drifting calibration on the upper limits of
$\varepsilon$ and $g_{a\gamma}$ by running analyses where we scale the sensor
readings by up to 10 percent of their values. We then determined the resulting
limits, concluding that a drifting calibration of the magnitude we observed
would change the limits on $\varepsilon$ and $g_{a\gamma}$ by $\lesssim 3\%$.

####  IV.3.4 Timing synchronization

As discussed in Sec. III, the magnetic-field measurements were digitized at
160 samples per second. An on-sensor real-time clock ensured sample-to-sample
timing to better than 1 ppm and a GPS-referenced computer clock provided the
absolute time reference for the time stamps. The absolute timing accuracy
between sensors was limited to $\sim 100\,\mathrm{ms}$ due to latencies in the
steering of the DAQ clock to GPS. This can be significantly improved. However,
such an accuracy was adequate for an analysis covering the 0.5 to 5 Hz window.
We estimate the systematic on the derived limits due to this error to be
neglible.

##  V Future Directions

The current experiment is limited by the sensitivity of the magnetometers,
rather than by the geomagnetic noise, and our model only accurately describes
signals at frequencies below $\approx 5$ Hz. In the next generation of the
experiment, we plan to use more sensitive magnetometers to reach the limit
imposed by geomagnetic noise. In addition, we propose to employ a novel
experimental geometry to avoid model uncertainties in interpretation of our
data.

At frequencies $\gtrsim 5$ Hz, the DM-induced magnetic field signal becomes
sensitive to the details of Earth’s atmosphere, which would require more
careful modelling than that needed for the lower-frequency analysis presented
in this paper. In order to be sensitive to higher-mass ALPs and hidden
photons, we are investigating the prospect of measuring spatial derivatives of
the magnetic field. By measuring components of the magnetic field across
multiple stations which are positioned $\lesssim 1~{}{\rm km}$ from one
another, it is possible to compute the numerical derivatives of ${\bm{B}}$,
and particularly components of $\nabla\times{\bm{B}}$. In the envisioned
measurement scheme, we do not expect to have significant local electric
currents, so the modified Ampère–Maxwell law describing the sought-after
effect of DM fields is

$\nabla\times{\bm{B}}-\partial_{t}{\bm{E}}={\bm{J}}_{\rm eff},$ (64)

where ${\bm{J}}_{\rm eff}$ encapsulates the effect of the dark matter [see
Eqs. (5) and (9)]. Since ${\bm{E}}$ is negligible in directions tangent to the
ground, a measurement of $\nabla\times\bm{B}$ in a tangent direction gives a
direct measurement of the dark matter, which is insensitive to the atmospheric
boundary conditions. Moreover, we expect this scheme to reduce sensitivity to
geomagnetic noise, as physical geomagnetic fields in the lower atmosphere
should have $(\nabla\times\bm{B})_{\parallel}=\bm{J}_{\parallel}=0$. However,
it is important to note that, unlike the low-frequency measurements whose
signal is enhanced by the full radius of Earth, the effective enhancement here
would only be the separation between stations. SNIPE Hunt is currently
carrying out an investigation of the expected background and signal, while
simultaneously taking steps to perform a search based on this new methodology.

##  VI Conclusions

In this work, we reported on a search for axion and hidden-photon dark matter
using a network of unshielded vector magnetoresistive (VMR) magnetometers
located in relatively quiet magnetic environments, in wilderness areas far
from anthropogenic magnetic noise. The magnetic signal pattern targeted by our
search could, in principle, be generated by the interaction of axion or hidden
photon dark matter with Earth, which can act as a transducer to convert the
dark matter into oscillating magnetic fields as described in Refs. [26, 27,
28]. Analysis of the data acquired over the course of approximately three days
in July 2022 revealed no evidence of a persistent oscillating magnetic field
matching the expected characteristics of a dark-matter-induced signal.
Consequently, we set upper limits on the kinetic-mixing parameter
$\varepsilon$ for hidden-photon dark matter and on the axion–photon coupling
constant $g_{a\gamma}$.

Figure 11: Constraints on the hidden-photon kinetic-mixing parameter
$\varepsilon$ as a function of hidden-photon mass $m_{A^{\prime}}$. The plot
was created based on Refs. [69] and [47] and includes the SuperMAG limit [26,
27] and the recent measurement using a network of magnetometers in meter-scale
shielded rooms [70], which we denote the “Synchronized Quantum Sensor Network”
(SQSN). The results reported in Refs. [26, 27, 70] are the only other
laboratory measurements in this mass range. In additional to the laboratory
constraints, the plot also shows various astrophysical bounds, including the
geomagnetic limit obtained from satellite measurements of the Earth’s magnetic
field [71], the hidden-photon limits from magnetic-field measurements in
Jupiter’s magnetosphere [72], limits from cold gas clouds at the Milky Way
center [73], heating of the ionized interstellar medium in the galaxy from
hidden photons [74], and the limit on heating/cooling due to DM in the Leo T
dwarf galaxy [75]. Cosmological bounds on hidden photons from COBE/FIRAS data
estimated from potential hidden-photon interactions with plasmas in the
universe are from Refs. [76, 77, 46]. Finally, the figure also displays
cosmological/astrophysical bounds on hidden photons from He II reionization
[78].

Figure 11 displays constraints on $\varepsilon$ as a function of hidden-photon
mass $m_{A^{\prime}}$ obtained in our experiment as well as those from other
experiments [27, 70], derived from planetary science [71, 72], and based on
astrophysical observations [73, 74, 75, 76, 77, 46, 78]. We note that, in the
studied frequency range, the results of the SNIPE Hunt experiment are the most
stringent experimental bounds, and can be regarded as complementary to the
more severe observational constraints. Fig. 12 shows bounds on the
axion–photon coupling constant parameter $g_{a\gamma}$ as a function of axion
mass $m_{a}$.

Figure 12: Constraints on the axion–photon coupling constant parameter
$g_{a\gamma}$ as a function of axion mass $m_{a}$. The plot was created based
on Ref. [69], and includes the relevant experimental bounds based on the
SuperMAG analysis [28] in maroon and the CAST result [79] in grey.
Additionally, the plot displays astrophysical limits on the axion–photon
interaction, represented in various shades of green, including (Diffuse SNe)
[80], (Hydra A) [81], (Super star clusters) [82], M87 [83], and (H1821+643)
[84].

We are actively pursuing further measurements based on this concept, but
instead using induction-coil magnetometers [85, 86, 87]. We anticipate an
improvement in sensitivity to dark-matter-induced magnetic signals of several
orders of magnitude. Furthermore, as discussed in Sec. V, we will use local
multi-sensor arrays to measure the curl of the local magnetic field at the
various sites and thereby extend the frequency range probed up to about a kHz.

###### Acknowledgements.

The Oberlin group thank Michael Miller for his work on the construction of the
sensor mount for the Oberlin station. This work was supported by the U.S.
National Science Foundation under grants PHY-2110370, PHY-2110385, and
PHY-2110388. S.K. and M.A.F. are supported by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers,
Superconducting Quantum Materials and Systems Center (SQMS) under contract
number DE-AC02-07CH11359. M.A.F. is also supported by the Simons Investigator
Award No. 827042. S.K. and M.A.F. thank the Aspen Center for Physics for
hospitality during the final stages of this work, supported by NSF Grant No.
PHY-2210452.

## References

  * Jackson Kimball and van Bibber [2022] D. F. Jackson Kimball and K. van Bibber,  _The Search for Ultralight Bosonic Dark Matter_ (Springer, 2022).
  * Graham  _et al._ [2015a] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Experimental searches for the axion and axion-like particles, Annu. Rev. Nucl. Part. Sci. 65, 485 (2015a).
  * Arias  _et al._ [2012a] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, WISPy Cold Dark Matter, JCAP 06, 013, arXiv:1201.5902 [hep-ph] .
  * Braaten and Zhang [2019] E. Braaten and H. Zhang, Colloquium: The physics of axion stars, Rev. Mod. Phys. 91, 041002 (2019).
  * Freese  _et al._ [2013] K. Freese, M. Lisanti, and C. Savage, Colloquium: Annual modulation of dark matter, Rev. Mod. Phys. 85, 1561 (2013).
  * Pillepich  _et al._ [2014] A. Pillepich, M. Kuhlen, J. Guedes, and P. Madau, The distribution of dark matter in the Milky Way’s disk, Astrophys. J. 784, 161 (2014).
  * Evans  _et al._ [2019] N. W. Evans, C. A. O’Hare, and C. McCabe, Refinement of the standard halo model for dark matter searches in light of the Gaia Sausage, Phys. Rev. D 99, 023012 (2019).
  * Hui  _et al._ [2017] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Ultralight scalars as cosmological dark matter, Phys. Rev. D 95, 043541 (2017).
  * Foster  _et al._ [2018] J. W. Foster, N. L. Rodd, and B. R. Safdi, Revealing the dark matter halo with axion direct detection, Phys. Rev. D 97, 123006 (2018).
  * Lin  _et al._ [2018] S.-C. Lin, H.-Y. Schive, S.-K. Wong, and T. Chiueh, Self-consistent construction of virialized wave dark matter halos, Phys. Rev. D 97, 103523 (2018).
  * Centers  _et al._ [2021] G. P. Centers, J. W. Blanchard, J. Conrad, N. L. Figueroa, A. Garcon, A. V. Gramolin, D. F. Jackson Kimball, M. Lawson, B. Pelssers, J. A. Smiga, _et al._ , Stochastic fluctuations of bosonic dark matter, Nature Comm. 12, 7321 (2021).
  * Lisanti  _et al._ [2021] M. Lisanti, M. Moschella, and W. Terrano, Stochastic properties of ultralight scalar field gradients, Phys. Rev. D 104, 055037 (2021).
  * Graham  _et al._ [2016a] P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, and W. A. Terrano, Dark matter direct detection with accelerometers, Phys. Rev. D 93, 075029 (2016a).
  * Safronova  _et al._ [2018] M. Safronova, D. Budker, D. DeMille, D. F. Jackson Kimball, A. Derevianko, and C. W. Clark, Search for new physics with atoms and molecules, Rev. Mod. Phys. 90, 025008 (2018).
  * Sikivie [1983] P. Sikivie, Experimental tests of the “invisible” axion, Phys. Rev. Lett. 51, 1415 (1983).
  * Asztalos  _et al._ [2010] S. J. Asztalos, G. Carosi, C. Hagmann, D. Kinion, K. Van Bibber, M. Hotz, L. Rosenberg, G. Rybka, J. Hoskins, J. Hwang, _et al._ , SQUID-based microwave cavity search for dark-matter axions, Phys. Rev. Lett. 104, 041301 (2010).
  * Braine  _et al._ [2020] T. Braine, R. Cervantes, N. Crisosto, N. Du, S. Kimes, L. Rosenberg, G. Rybka, J. Yang, D. Bowring, A. Chou, _et al._ , Extended search for the invisible axion with the axion dark matter experiment, Phys. Rev. Lett. 124, 101303 (2020).
  * Zhong  _et al._ [2018] L. Zhong, S. Al Kenany, K. Backes, B. Brubaker, S. Cahn, G. Carosi, Y. Gurevich, W. Kindel, S. Lamoreaux, K. Lehnert, _et al._ , Results from phase 1 of the HAYSTAC microwave cavity axion experiment, Phys. Rev. D 97, 092001 (2018).
  * Backes  _et al._ [2021] K. M. Backes, D. A. Palken, S. A. Kenany, B. M. Brubaker, S. Cahn, A. Droster, G. C. Hilton, S. Ghosh, H. Jackson, S. K. Lamoreaux, _et al._ , A quantum enhanced search for dark matter axions, Nature 590, 238 (2021).
  * Salemi  _et al._ [2021] C. P. Salemi, J. W. Foster, J. L. Ouellet, A. Gavin, K. M. Pappas, S. Cheng, K. A. Richardson, R. Henning, Y. Kahn, R. Nguyen, _et al._ , Search for low-mass axion dark matter with ABRACADABRA-10 cm, Phys. Rev. Lett. 127, 081801 (2021).
  * Gramolin  _et al._ [2021] A. V. Gramolin, D. Aybas, D. Johnson, J. Adam, and A. O. Sushkov, Search for axion-like dark matter with ferromagnets, Nature Physics 17, 79 (2021).
  * Andrew  _et al._ [2023] K. Y. Andrew, S. Ahn, Ç. Kutlu, J. Kim, B. R. Ko, B. I. Ivanov, H. Byun, A. F. van Loo, S. Park, J. Jeong, _et al._ , Axion dark matter search around 4.55 $\mu$ ev with dine-fischler-srednicki-zhitnitskii sensitivity, Phys. Rev. Lett. 130, 071002 (2023).
  * Wagner  _et al._ [2010] A. Wagner, G. Rybka, M. Hotz, L. Rosenberg, S. Asztalos, G. Carosi, C. Hagmann, D. Kinion, K. Van Bibber, J. Hoskins, _et al._ , Search for hidden sector photons with the ADMX detector, Phys. Rev. Lett. 105, 171801 (2010).
  * Chaudhuri  _et al._ [2015] S. Chaudhuri, P. W. Graham, K. Irwin, J. Mardon, S. Rajendran, and Y. Zhao, Radio for hidden-photon dark matter detection, Phys. Rev. D 92, 075012 (2015).
  * Phipps  _et al._ [2020] A. Phipps, S. E. Kuenstner, S. Chaudhuri, C. S. Dawson, B. A. Young, C. T. FitzGerald, H. Froland, K. Wells, D. Li, H. Cho, _et al._ , Exclusion limits on hidden-photon dark matter near 2 neV from a fixed-frequency superconducting lumped-element resonator, in  _Microwave Cavities and Detectors for Axion Research: Proceedings of the 3rd International Workshop_ (Springer, 2020) p. 139.
  * Fedderke  _et al._ [2021a] M. A. Fedderke, P. W. Graham, D. F. Jackson Kimball, and S. Kalia, Earth as a transducer for dark-photon dark-matter detection, Phys. Rev. D 104, 075023 (2021a).
  * Fedderke  _et al._ [2021b] M. A. Fedderke, P. W. Graham, D. F. Jackson Kimball, and S. Kalia, Search for dark-photon dark matter in the supermag geomagnetic field dataset, Phys. Rev. D 104, 095032 (2021b).
  * Arza  _et al._ [2022] A. Arza, M. A. Fedderke, P. W. Graham, D. F. Jackson Kimball, and S. Kalia, Earth as a transducer for axion dark-matter detection, Phys. Rev. D 105, 095007 (2022).
  * Primakoff [1951] H. Primakoff, Photo-production of neutral mesons in nuclear electric fields and the mean life of the neutral meson, Phys. Rev. 81, 899 (1951).
  * Raffelt and Seckel [1988] G. Raffelt and D. Seckel, Bounds on exotic-particle interactions from SN1987A, Phys. Rev. Lett. 60, 1793 (1988).
  * [31] https://en.wikipedia.org/wiki/Snipe_hunt.
  * Bowen  _et al._ [2019] T. A. Bowen, E. Zhivun, A. Wickenbrock, V. Dumont, S. D. Bale, C. Pankow, G. Dobler, J. S. Wurtele, and D. Budker, A network of magnetometers for multi-scale urban science and informatics, Geosci. Instrum. Methods Data Syst. 8, 129 (2019).
  * Dumont  _et al._ [2022] V. Dumont, T. A. Bowen, R. Roglans, G. Dobler, M. S. Sharma, A. Karpf, S. D. Bale, A. Wickenbrock, E. Zhivun, T. Kornack, J. S. Wurtele, and D. Budker, Do cities have a unique magnetic pulse?, J. of Appl. Phys. 131, 204902 (2022).
  * [34] http://supermag.jhuapl.edu.
  * Gjerloev [2009] J. W. Gjerloev, A Global Ground-Based Magnetometer Initiative, Eos 90, 230 (2009).
  * Gjerloev [2012] J. W. Gjerloev, The SuperMAG data processing technique, J. Geophys. Res. Space Phys. 117, A09213 (2012).
  * Constable and Constable [2004] C. G. Constable and S. C. Constable, Satellite magnetic field measurements: Applications in studying the deep earth, in  _The State of the Planet: Frontiers and Challenges in Geophysics_ (American Geophysical Union, AGU, 2004) pp. 147–159.
  * Cvetic and Langacker [1996] M. Cvetic and P. Langacker, Implications of Abelian extended gauge structures from string models, Phys. Rev. D 54, 3570 (1996).
  * Holdom [1986] B. Holdom, Two U(1)’s and $\epsilon$ charge shifts, Phys. Lett. B 166, 196 (1986).
  * Graham  _et al._ [2014] P. W. Graham, J. Mardon, S. Rajendran, and Y. Zhao, Parametrically enhanced hidden photon search, Phys. Rev. D 90, 075017 (2014).
  * Graham  _et al._ [2016b] P. W. Graham, J. Mardon, and S. Rajendran, Vector Dark Matter from Inflationary Fluctuations, Phys. Rev. D 93, 103520 (2016b).
  * Ahmed  _et al._ [2020] A. Ahmed, B. Grzadkowski, and A. Socha, Gravitational production of vector dark matter, J. High Energy Phys. 2020 (8), 59.
  * Kolb and Long [2021] E. W. Kolb and A. J. Long, Completely dark photons from gravitational particle production during the inflationary era, J. High Energy Phys. 2021 (3), 283.
  * Adshead  _et al._ [2023] P. Adshead, K. D. Lozanov, and Z. J. Weiner, Dark photon dark matter from an oscillating dilaton, (2023), arXiv:2301.07718 [hep-ph] .
  * Nelson and Scholtz [2011] A. E. Nelson and J. Scholtz, Dark light, dark matter, and the misalignment mechanism, Phys. Rev. D 84, 103501 (2011).
  * Arias  _et al._ [2012b] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, Wispy cold dark matter, J. Cosmol. Astropart. Phys. 06 (06), 013.
  * Caputo  _et al._ [2021] A. Caputo, C. A. J. O’Hare, A. J. Millar, and E. Vitagliano, Dark photon limits: a cookbook, (2021), arXiv:2105.04565 [hep-ph] .
  * Read [2014] J. I. Read, The local dark matter density, J. Phys. G: Nucl. Part. Phys. 41, 063101 (2014).
  * Bland-Hawthorn and Gerhard [2016] J. Bland-Hawthorn and O. Gerhard, The galaxy in context: structural, kinematic, and integrated properties, Annu. Rev. Astron. Astrophys. 54, 529 (2016).
  * Peccei and Quinn [1977] R. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg [1978] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Preskill  _et al._ [1983a] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120, 127 (1983a).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120, 137 (1983).
  * Graham and Scherlis [2018] P. W. Graham and A. Scherlis, Stochastic axion scenario, Phys. Rev. D 98, 035017 (2018).
  * Preskill  _et al._ [1983b] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the invisible axion, Phys. Lett. B 120, 127 (1983b).
  * Svrcek and Witten [2006] P. Svrcek and E. Witten, Axions in string theory, J. High Energy Phys. 2006 (06), 051.
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String axiverse, Phys. Rev. D 81, 123530 (2010).
  * Graham  _et al._ [2015b] P. W. Graham, D. E. Kaplan, and S. Rajendran, Cosmological relaxation of the electroweak scale, Phys. Rev. Lett. 115, 221801 (2015b).
  * Sentman [2017] D. D. Sentman, Schumann resonances, in  _Handbook of atmospheric electrodynamics_ (CRC Press, 2017) pp. 267–295.
  * Alken  _et al._ [2021] P. Alken, E. Thébault, C. D. Beggan, H. Amit, J. Aubert, J. Baerenzung, T. N. Bondar, W. J. Brown, S. Califf, A. Chambodut, A. Chulliat, G. A. Cox, C. C. Finlay, A. Fournier, N. Gillet, A. Grayver, M. D. Hammer, M. Holschneider, L. Huder, G. Hulot, T. Jager, C. Kloss, M. Korte, W. Kuang, A. Kuvshinov, B. Langlais, J.-M. Léger, V. Lesur, P. W. Livermore, F. J. Lowes, S. Macmillan, W. Magnes, M. Mandea, S. Marsal, J. Matzka, M. C. Metman, T. Minami, A. Morschhauser, J. E. Mound, M. Nair, S. Nakano, N. Olsen, F. J. Pavón-Carrasco, V. G. Petrov, G. Ropp, M. Rother, T. J. Sabaka, S. Sanchez, D. Saturnino, N. R. Schnepf, X. Shen, C. Stolle, A. Tangborn, L. Tøffner-Clausen, H. Toh, J. M. Torta, J. Varner, F. Vervelidou, P. Vigneron, I. Wardinski, J. Wicht, A. Woods, Y. Yang, Z. Zeren, and B. Zhou, International Geomagnetic Reference Field: the thirteenth generation, Earth Planets Space 73, 49 (2021).
  * [63] U. G. Survey, Usgs geomagnetic survey, https://geomag.usgs.gov/ws/docs, accessed: 2023-03-10.
  * [64] Solar storm forecast, https://www.space.com/moderate-solar-storm-forecast-july-23, accessed: 2023-03-10.
  * atl [2015] World Atlas of Ground Conductivities, Rec. ITU-R P.832-4, International Telecommunication Union (2015).
  * Takeda and Araki [1985] M. Takeda and T. Araki, Electric conductivity of the ionosphere and nocturnal currents, J. Atmos. Terr. Phys. 47, 601 (1985).
  * Richmond and Thayer [2000] A. Richmond and J. Thayer, Ionospheric Electrodynamics: A Tutorial, in  _Magnetospheric Current Systems (Geophysical Monograph 118)_ , edited by S. Ohtani, R. Fujii, M. Hesse, and R. L. Lysak (American Geophysical Union, Washington, DC, 2000) pp. 131–146.
  * Rodríguez-Camacho  _et al._ [2022] J. Rodríguez-Camacho, A. Salinas, M. C. Carrión, J. Portí, J. Fornieles-Callejón, and S. Toledo-Redondo, Four year study of the schumann resonance regular variations using the sierra nevada station ground-based magnetometers, Journal of Geophysical Research: Atmospheres 127, e2021JD036051 (2022).
  * [69] https://cajohare.github.io/AxionLimits/.
  * Jiang  _et al._ [2023] M. Jiang, T. Hong, D. Hu, Y. Chen, F. Yang, T. Hu, X. Yang, J. Shu, Y. Zhao, and X. Peng, Search for dark photons with synchronized quantum sensor network, arXiv:2305.00890 10.48550/arXiv.2305.00890 (2023).
  * Fischbach  _et al._ [1994] E. Fischbach, H. Kloor, R. A. Langel, A. T. Y. Lui, and M. Peredo, New geomagnetic limits on the photon mass and on long-range forces coexisting with electromagnetism, Phys. Rev. Lett. 73, 514 (1994).
  * Marocco [2021] G. Marocco, Dark photon limits from magnetic fields and astrophysical plasmas, (2021), arXiv:2110.02875 [hep-ph] .
  * Bhoonah  _et al._ [2019] A. Bhoonah, J. Bramante, F. Elahi, and S. Schon, Galactic center gas clouds and novel bounds on ultralight dark photon, vector portal, strongly interacting, composite, and super-heavy dark matter, Phys. Rev. D 100, 023001 (2019).
  * Dubovsky and Hernández-Chifflet [2015] S. Dubovsky and G. Hernández-Chifflet, Heating up the galaxy with hidden photons, J. Cosmol. Astropart. Phys. 2015 (12), 054.
  * Wadekar and Farrar [2021] D. Wadekar and G. R. Farrar, Gas-rich dwarf galaxies as a new probe of dark matter interactions with ordinary matter, Phys. Rev. D 103, 123028 (2021).
  * Caputo  _et al._ [2020a] A. Caputo, H. Liu, S. Mishra-Sharma, and J. T. Ruderman, Dark photon oscillations in our inhomogeneous universe, Phys. Rev. Lett. 125, 221303 (2020a).
  * McDermott and Witte [2020] S. D. McDermott and S. J. Witte, Cosmological evolution of light dark photon dark matter, Phys. Rev. D 101, 063030 (2020).
  * Caputo  _et al._ [2020b] A. Caputo, H. Liu, S. Mishra-Sharma, and J. T. Ruderman, Dark photon oscillations in our inhomogeneous universe, Phys. Rev. Lett. 125, 221303 (2020b).
  * CAST Collaboration  _et al._ [2017] CAST Collaboration, V. Anastassopoulos, S. Aune, K. Barth, A. Belov, H. Bräuninger, G. Cantatore, J. M. Carmona, J. F. Castel, S. A. Cetin, F. Christensen, J. I. Collar, T. Dafni, M. Davenport, T. A. Decker, A. Dermenev, K. Desch, C. Eleftheriadis, G. Fanourakis, E. Ferrer-Ribas, H. Fischer, J. A. García, A. Gardikiotis, J. G. Garza, E. N. Gazis, T. Geralis, I. Giomataris, S. Gninenko, C. J. Hailey, M. D. Hasinoff, D. H. H. Hoffmann, F. J. Iguaz, I. G. Irastorza, A. Jakobsen, J. Jacoby, K. Jakovčić, J. Kaminski, M. Karuza, N. Kralj, M. Krčmar, S. Kostoglou, C. Krieger, B. Lakić, J. M. Laurent, A. Liolios, A. Ljubičić, G. Luzón, M. Maroudas, L. Miceli, S. Neff, I. Ortega, T. Papaevangelou, K. Paraschou, M. J. Pivovaroff, G. Raffelt, M. Rosu, J. Ruz, E. R. Chóliz, I. Savvidis, S. Schmidt, Y. K. Semertzidis, S. K. Solanki, L. Stewart, T. Vafeiadis, J. K. Vogel, S. C. Yildiz, and K. Zioutas, New CAST limit on the axion–photon interaction, Nature Physics 13, 584 (2017).
  * Caputo  _et al._ [2022] A. Caputo, H.-T. Janka, G. Raffelt, and E. Vitagliano, Low-energy supernovae severely constrain radiative particle decays, Phys. Rev. Lett. 128, 221103 (2022).
  * Wouters and Brun [2013] D. Wouters and P. Brun, Constraints on axion-like particles from x-ray observations of the hydra galaxy cluster, Astrophys. J. 772, 44 (2013).
  * Dessert  _et al._ [2020] C. Dessert, J. W. Foster, and B. R. Safdi, X-ray searches for axions from super star clusters, Phys. Rev. Lett. 125, 261102 (2020).
  * Marsh  _et al._ [2017] M. D. Marsh, H. R. Russell, A. C. Fabian, B. R. McNamara, P. Nulsen, and C. S. Reynolds, A new bound on axion-like particles, J. Cosmol. Astropart. Phys. 2017 (12), 036.
  * Sisk Reynés  _et al._ [2021] J. Sisk Reynés, J. H. Matthews, C. S. Reynolds, H. R. Russell, R. N. Smith, and M. C. D. Marsh, New constraints on light axion-like particles using _Chandra_ transmission grating spectroscopy of the powerful cluster-hosted quasar H1821+643, Mon. Not. R. Astron. Soc. 510, 1264 (2021).
  * Votis  _et al._ [2018] C. I. Votis, G. Tatsis, V. Christofilakis, S. K. Chronopoulos, P. Kostarakis, V. Tritakis, and C. Repapis, A new portable ELF Schumann resonance receiver: Design and detailed analysis of the antenna and the analog front-end, EURASIP Journal on Wireless Communications and Networking 2018, 155 (2018).
  * Poliakov  _et al._ [2017] S. Poliakov, B. Reznikov, A. Shchennikov, E. Kopytenko, and B. Samsonov, The range of induction-coil magnetic field sensors for geophysical explorations, Seismic instruments 53, 1 (2017).
  * Hospodarsky [2016] G. B. Hospodarsky, Spaced-based search coil magnetometers, Journal of Geophysical Research: Space Physics 121, 12068 (2016).