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

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 ≪ 1 eV/c2 ) such as hidden photons, axions, or axion-like particles (ALPs) [1][2][3].If ultralight bosons are the dark matter, under reasonable assumptions 1 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 * ibrahim.sulai@bucknell.edu† kalias@umn.edu‡ derek.jacksonkimball@csueastbay.edu 1 Here 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].
In this paper we describe initial results of the "Search for Non-Interacting Particles Experimental Hunt" (SNIPE Hunt [31]): a campaign to search for axion 2 and hidden-photon dark matter using magnetome-ters located in the "wilderness" (away from the high levels of magnetic noise associated with urban environments [32,33]).This work extends to higher axion/hiddenphoton 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 ′ 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 wavelike 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]: where only terms up to first order in the kinetic mixing parameter ε ≪ 1 are retained.In Eq. ( 1), F µν is the field-strength tensor for the "interacting" mode of the electromagnetic field that couples to charges, (F ′ ) µν is the field-strength tensor for the "sterile" mode that does not interact with charges, A µ is the four-potential for the interacting mode, (A ′ ) µ is the four-potential for the sterile mode, and J µ EM is the electromagnetic four-current density.In our case of interest, the hidden-photon darkmatter field in the vicinity of Earth is a coherently oscillating vector field with random polarization: 4 where A ′ is the sterile vector potential, ρ DM ≈ 0.3 GeV/cm 3 is the local dark-matter density [48], ni are a set of orthonormal unit vectors, ξ i (r, t) are slowly varying O(1) amplitudes, and ϕ i (r, t) are slowly varying random phases.Both the amplitudes ξ i (r, t) and phases ϕ i (r, t) of the hidden-photon dark-matter field change stochastically on length scales given by the dark-matter coherence length, and time scales given by the coherence time of the field, where v DM ∼ 10 −3 is the characteristic dispersion (virial) velocity of the dark matter in the vicinity of Earth [7,49].Note that the timelike component of the four-potential (A ′ ) µ is suppressed relative to the spacelike component 3 Throughout, we use natural units where ℏ = c = 1. 4 In 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.
(the vector potential A ′ ) by ∼ v DM ∼ 10 −3 .From inspection of Eq. ( 1), it can be seen that the physical effects due to the hidden-photon dark-matter field (A ′ ) µ are to leading order the same as those generated by an effective current density Inside a good conductor, the interacting mode vanishes, F µν = 0 and A µ = 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 nonzero 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 (≪ 1 eV) and couplings to Standard Model particles and fields [2].Like hidden photons, axions are ubiquitous features of beyond-the-Standard-Model theories [53,[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: where a is the axion field, m a is the axion mass, g aγ parameterizes the axion-photon coupling, and F µν is the dual field-strength tensor.The last term appearing in Eq. ( 6) describes the interaction between the axion and electromagnetic fields: where E and 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: It follows that the physical effects of the axion-photon coupling in the presence of a magnetic field B, as in the case of hidden photons [Eq.( 5)], manifest as an effective current: where a(r, t) = a 0 (r, t)e −imat (10) is the axion field with a stochastically (slowly) varying amplitude |a 0 | ∼ √ 2ρ DM /m a , with coherence length ℓ coh and coherence time τ coh analogous to those for hidden photons described by Eqs. ( 3) and ( 4), with the replacement m A ′ → m a .The interaction of an axion darkmatter 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 ≈ 0.5 Hz and there is ongoing analysis of SuperMAG data covering frequencies up to ≈ 1 Hz that is expected to surpass the sensitivity of this experiment.For the upper bound of 5 Hz, we are limited by the well-studied Schumann resonances of the Earth-ionosphere cavity [61,62].We cannot make a robust prediction for frequencies corresponding to the Schumann resonances because of finite conductivity effects and 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 ≤ 5 Hz a region where the dark-matter-induced magnetic-field pattern can be reliably derived (see Sec. IV C 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: λ A ′ ≈ 1/m A ′ ≫ R and λ a ≈ 1/m a ≫ R.This sets an upper limit on the hidden-photon mass m A ′ and axion mass m a of ∼ 3 × 10 −14 eV and, correspondingly, for their Compton frequencies: f A ′ and f a must be ≪ 7 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 ∼ 1/R.We write the magnetic and electric fields in terms of vector spherical harmonics (VSH; see Appendix D of [26]) Y ℓm , Ψ ℓm , Φ ℓm as where ω is the oscillation angular frequency of the darkmatter effective current.For the dark-matter effective current J which stands for both hidden photons and axion-like particles, we use the fact that it satisfies ∇ × J = 0 to write Inserting the above ansatz into Maxwell's equations, we get and the other components are determined by ℓm + J (1) ℓm ( 16) B ℓm .
This system is solved with boundary conditions such that E (1) ℓm and E (2) ℓ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 ωh ≪ 1, the boundary condition for E (2) ℓm implies immediately that it is zero everywhere; it follows that B (r) ℓm and B (1) ℓm also vanish identically.
Writing B (2) ℓm = u ℓm /r, in the limit in which h ≪ R we find where We write the solution for Notice that the magnetic field signal at Earth's surface (r = R) is simply given by From the boundary condition u ℓm at r = R and r = R + h, we find at zeroth order in h/R A. Hidden-Photon Signal In terms of vector spherical harmonics, the hiddenphoton effective current, given in Eq. ( 5), is written as Here , where f d is the frequency associated to the sidereal day, 5 and the hidden-photon amplitudes A ′ m (for polarizations m = 0, ±1) appearing in Eq. ( 22) are normalized via where ρ DM = 0.3 GeV/cm 3 is the local dark-matter density.Extracting J 1m from Eq. ( 22), we find

B. 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 [63], which parameterizes Earth's magnetic field B ⊕ in terms of a scalar potential V 0 , such that B ⊕ = −∇V 0 , where V 0 is expanded as where P m ℓ are the Schmidt-normalized associated Legendre polynomials.The Gauss coefficients g ℓm and h ℓm are specified by the IGRF model at five-year intervals (see Tab. 2 of Ref. [63]).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 ℓ = 4) forward to July 23, 2022 using these time derivatives, and adopt the conventions 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γ generates can be written as [28] where a 0 is the (complex) axion amplitude, normalized by 1  2 m 2 a ⟨|a 0 | 2 ⟩ = ρ DM , and Now, by identifying J (1) (r) in Eq. ( 26), the magneticfield signal from axion dark matter is found to be

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 pT/ √ Hz over a frequency range of 0.1-100 Hz.Prior to deploying the sensors in the field, we verified the calibration of the magnetometers with well-known external oscillating fields applied to the sensors within a magnetically shielded environment.An accurate determination of the oscillating magnetic fields used for calibration was independently attained by observing and measuring magneto-optical resonances in alkali-metal vapor magnetometers [64][65][66].
In addition to the magnetic field sensors, 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 A • 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 ≲ 100 ms.
The location of the three stations is shown in Table I.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 ≲ 1 • .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.A threeaxis 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.

A. 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 ampli- tude spectral density in most of the band of interest is flat and corresponds to approximately 300 pT/ √ Hz, the noise floor of the sensors.
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 ∼ 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 C, we assign an uncertainty on the quoted HPDM and axion limits due to temperature drifts.
Between hours ∼ 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 ∼ 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 [67,68].The storm led to the modulation of Earth's magnetic field which we detected.Including data from this window in 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.

A. 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: and B ϕ (Ω 3 , t j ). 6We model these time series as being given by (the real part of) the signal in Eq. ( 24) plus Gaussian white noise.Our goal is then to extract a bound on ε.As the exact amplitudes A ′ m are unknown, we utilize a Bayesian framework and treat these as nuisance parameters.We also take a Gaussian distribution for them, 7 normalized by Eq. (23).
The signal in Eq. ( 24) indicates that all relevant information is contained at the frequencies f A ′ and f A ′ ± f d .Thus we Fourier transform the six time series B α (Ω i ), and construct an 18-dimensional data vector 8 ⃗ X which contains all information which may be relevant to setting a bound at f A ′ .Namely, ⃗ X consists of the six values Bα Ω i , f A ′ − fd , followed by the six values Bα (Ω i , f A ′ ), followed by the six values Bα Ω i , f A ′ + fd .In our analysis, we compute bounds only at discrete Fourier transform (DFT) frequencies f A ′ = 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 fd , which we define as the nearest DFT frequency to f d .With these choices, ⃗ X can be computed via a fast Fourier transform (FFT).(This allows us to compute ⃗ 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 ⃗ X, namely its expectation and variance.
First, let us compute the expectation of ⃗ X.As mentioned above, we model our measurements as being Gaussian noise on top of the signal in Eq. (24).Since the expectation of the noise vanishes, the expectation of ⃗ X simply comes from Fourier transforming Eq. ( 24) and assembling its relevant components into a vector.To remove the normalization from the amplitudes A ′ m , let us define These now have m ⟨|c m | 2 ⟩ = 1.In the case c ± = 0, (the real part of) Eq. ( 24) takes the simple form (30) and the only nonzero components of ⟨ ⃗ X⟩ are 7 A ′ can be written as a sum of several independent plane wave solutions of different velocities vn ∼ O(v DM ).These have corresponding frequencies fn ∼ f A ′ 1 + O(v 2 DM ) .On timescales longer than τ coh ∼ 1/(f A ′ v 2 DM ), the value of A ′ is thus a sum of many contributions with random phases.By the central limit theorem, it is thus distributed as a Gaussian variable. 8We use ⃗ x to denote a vector x with 18 components (or six components in Sec.IV B), and y to indicate a vector y with three components.
On the other hand, if c 0 = c − = 0, then the signal becomes and so the expectation of ⃗ X is where and ∆t = (1/160) s is the time resolution.Note that, in principle, Eq. ( 35) should have an additional term proportional to c + , which contains factors of Q f , these will all be significantly smaller than the Q factors appearing in Eq. ( 35).Thus we are safe to neglect this additional term.Similarly, ⟨ ⃗ X⟩ − ≡ c * − ε⃗ µ − can be computed (for the case when c 0 = c + = 0).Then generically, the full expectation of ⃗ X is Now that we have computed the expectation of ⃗ X, let us consider its variance.In this analysis, we consider the frequency range 0.5 Hz ≤ f A ′ ≤ 5 Hz, over which the noise is roughly frequency independent [see Fig. (4)].Therefore, we may consider each instance of ⃗ 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 ⃗ X as where f k indexes the DFT frequencies between 0.5 Hz and 5 Hz (for k = 1, . . ., N ∼ 10 5 ).9 Now that we understand the statistics of ⃗ X, we can write down its likelihood From this likelihood, the computation of the bound on ε proceeds as in Sec.V D of Ref. [27], but we reproduce it here for completeness.Let us write Σ = LL † and then define If we let N be the 18 × 3 matrix whose columns are ⃗ ν m , then Eq. ( 39) becomes Now if we perform a singular value decomposition N = U SV † (where U is a 18 × 3 matrix with orthonormal columns, S is a 3 × 3 diagonal matrix, and V is a 3 × 3 unitary matrix) and further define then the likelihood in Eq. ( 42) can be reduced to 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 Marginalizing over d, the likelihood Eq. ( 45) reduces to where z m are the components of Z and s m are the diagonal entries of S [see Appendix D 1 of Ref. [27] for a derivation of Eq. ( 47)].In order to turn this into a posterior on ε, we must assume some prior.We take a Jeffreys prior over, the three diagonal blocks should be identical, since they correspond to the same averages in Eq. ( 38) (only with the frequency f k shifted by fd ).Thus it suffices to only compute Σ ij for 7 ≤ i, j ≤ 12.
again see Appendix D 1 of Ref. [27].The posterior for ε is thus where N must be calculated to normalize the integral of p(ε|Z) to 1.We then set a 95% credible upper limit ε by solving 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.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. ( 45), we see that under the null hypothesis of no dark matter signal (ε = 0), the vector Z should be distributed as a multivariate Gaussian of mean zero.Specifically, the statistic should follow a χ 2 -distribution with six degrees of freedom.We may therefore compute the corresponding local p-value where F χ 2 (ν) denotes the cumulative distribution function for a χ 2 -distribution with ν degrees of freedom.Fig. (7) shows the local p-values at each frequency f A ′ 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 crit defined by 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 frequency bins 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 2 DM , where v DM is the dark matter velocity dispersion.Since the frequency bin size for our analysis is roughly 10 −5 Hz and each cluster spans multiple bins, these clusters represent signal candidates with widths of roughly 10 −5 f a , corresponding to large velocity dispersions of v DM ∼ 1000 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 Hz ≤ f A ′ ≤ 5 Hz range.
We have verified our entire analysis by injecting artificial HPDM signals into our data set and ensuring that the analysis correctly identified them.For example, when we added a monochromatic signal of the form in Eq. ( 24) with ε = 10 −5 and m A ′ = 10 −14 eV to the time series data from each station, and re-ran our analysis, we found the resulting limit only changed in the vicinity of m A ′ = 10 −14 eV, where it became ε ∼ 1.4 × 10 −5 .(Note that the limit is slightly weaker than the injected signal, as expected.)Moreover, the candidate analysis correctly identified DM candidates near the injected masses with high significance.We applied a similar verification process to the axion analysis described in the next section.

B. 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 ⃗ X consisting of Fourier transforms of the measured magnetic field at each location.Since the axion signal in Eq. ( 28) contains no f d dependence, however, the only relevant information is contained at frequency f a .Therefore in this analysis, we only take ⃗ X to be a six-dimensional vector, consisting of the measurements: where Φ θ ℓm and Φ ϕ ℓm denote the θ-component and φcomponents of the VSH Φ ℓm , and The covariance matrix Σ of ⃗ X can again be determined by averaging over independent frequencies, as in Eq. (37) [except that Σ will now be a 6 × 6 matrix].If we define ⃗ Y and ⃗ ν as in Eqs. ( 40) and (41) [without the m index], and further define we can write the likelihood function for the axion signal as Again marginalizing over c (which we take to have a Gaussian distribution with ⟨|c| 2 ⟩ = 1), and utilizing a Jeffreys prior for g aγ , we arrive at the posterior distribution ) Note that Eq. ( 59) is properly normalized, which is possible because its integral over g aγ can be taken analytically.The 95% credible limit ĝaγ can then be defined, as in Eq. (50).In this case, we can solve for it analytically to find ĝaγ = 1 s − |z| 2 log 0.95 + 0.05e −|z| 2 − 1.
(60) 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 ĝaγ → 4.36/s in the limit z → 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γ asymptotes to a  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.
finite floor. 10s 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 pvalue at a particular frequency under the null hypothesis (g aγ = 0) as (The χ 2 -distribution only has two degrees of freedom now, since the likelihood in Eq. ( 58) 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 crit , as defined in Eq. ( 53).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 Hz ≤ f a ≤ 5 Hz range.

C. Error Budget
The results of this science run and analysis are summarized in Figs. 6 and 9.They show upper limits on ε, the HPDM kinetic mixing parameter, and on g aγ , the axion-photon coupling constant, respectively.Below, we discuss the impact of uncertainties in the signal model and experimental conditions on the quoted limits.

Signal model uncertainty
The signals in Eqs. ( 24) and ( 28) 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 [69] and the ionosphere [70,71] achieve conductivities of at least 10 −4 S/m at certain depths/heights, which translate to skin depths of ∼ 50 km for frequencies f ∼ 1 Hz.Given that the only relevant length scale appearing in Eqs. ( 24) and ( 28) is the radius of Earth R ∼ 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. ( 24) and ( 28), however, reveals that our model predicts resonances in the signal at mR = ℓ(ℓ + 1) (for ℓ = 1 in the HPDM case, and ℓ ≥ 1 in the axion case).These are the Schumann resonances of the Earth-ionosphere cavity [61,62].Our simplified spherical model predicts the first of these resonances to occur at ∼ 10 Hz, but the central frequency of this resonance has been measured to be ∼ 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 ≤ 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 ∼ 4 Hz [62].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 A were at f ∼ 0.5, 0.75 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.

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 ≲ 1 • .By adjusting the orientation of the sensor in the analysis, we estimate that the impact of such an orientation error is to change the ε and g aγ upper limits by ≲ 1%.

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 0.5-5.0Hz band, we estimate the impact of a possibly drifting calibration on the upper limits of ε and g aγ by  running analyses where we independently scaled the sensor readings by up to 10 percent for Hayward, and up to 3 percent for the other two stations.We then determined the resulting limits, concluding that a drifting calibration of the magnitude we observed would change the limits on ε and g aγ by ≲ 3%.

Timing synchronization
As discussed in Sec.III, the magnetic-field measurements were digitized at 160 samples per second.An onsensor 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 ∼ 100 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 ≈ 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 ≳ 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 ≲ 1 km from one another, it is possible to compute the numerical derivatives of B, and particularly components of ∇ × B. In the envisioned measurement scheme [72], 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 where J eff encapsulates the effect of the dark matter [see Eqs. ( 5) and ( 9)].Since E is negligible in directions tangent to the ground, a measurement of ∇ × B in a tangent direction gives a direct measurement of the dark matter, which is insensitive to the atmospheric bound-ary conditions.Moreover, we expect this scheme to reduce sensitivity to geomagnetic noise, as physical geomagnetic fields in the lower atmosphere should have (∇ × B) ∥ = J ∥ = 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 ε for hidden-photon dark matter and on the axion-photon coupling constant g aγ .
Figure 11 displays constraints on ε as a function of hidden-photon mass m A ′ obtained in our experiment as well as those from other experiments [27,74], derived from planetary science [75,76], and based on astrophysical observations [46,[77][78][79][80][81][82].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γ as a function of axion mass m a .
We are actively pursuing further measurements based on this concept, but instead using induction-coil magnetometers [89][90][91].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.

FIG. 1 .
FIG.1.Block diagram of SNIPE station setup.A threeaxis 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.

FIG. 3 .FIG. 4 .
FIG. 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.

FIG. 6 .
FIG. 6. 95% credible upper limit on ε, 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.

FIG. 7 .
FIG. 7. The local p0-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. 53).The left panels show p0 as a function of frequency with candidates having p-values below the threshold.The right panels show histograms of p0 for the two different scans and candidates as outliers to the right of the threshold.

FIG. 8 .
FIG. 8.The local p0-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.

FIG. 10 .
FIG. 10.The local p0-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).pcrit, 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 p0-values for each scan.Signal candidates would appear as outliers to the right of the threshold.

TABLE I .
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.
FIG. 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 heavyduty plastic tent screws.