Earth as a transducer for axion dark-matter detection

We demonstrate that ultralight axion dark matter with a coupling to photons induces an oscillating global terrestrial magnetic field signal in the presence of the background geomagnetic field of the Earth. This signal is similar in structure to that of dark-photon dark matter that was recently pointed out and searched for in [arXiv:2106.00022] and [arXiv:2108.08852]. It has a global vectorial pattern fixed by the Earth's geomagnetic field, is temporally coherent on long time scales, and has a frequency set by the axion mass $m_a$. In this work, we both compute the detailed signal pattern, and undertake a search for this signal in magnetometer network data maintained by the SuperMAG Collaboration. Our analysis identifies no strong evidence for an axion dark-matter signal in the axion mass range $2\times10^{-18}\text{eV} \lesssim m_a \lesssim 7\times10^{-17}\text{eV}$. Assuming the axion is all of the dark matter, we place constraints on the axion-photon coupling $g_{a\gamma}$ in the same mass range; at their strongest, for masses $3\times 10^{-17}\text{eV} \lesssim m_a \lesssim 4\times 10^{-17}\text{eV}$, these constraints are comparable to those obtained by the CAST helioscope.

The identity of the dark matter (DM) [1,2] remains one of the most prominent unsolved puzzles about the Universe. Based on its gravitational effects, we know that the dark matter composes about 26% of the total energy density of the Universe, the rest consisting of 4% ordinary baryonic matter and 70% dark energy [3][4][5][6]. To date, no non-gravitational interactions of the dark matter with Standard Model (SM) particles have been observed; if these interactions exist, they must thus be feeble. While the ongoing and multi-faceted weakly interacting massive particle (WIMP) detection program continues its long-term [7] search for ever-weaker couplings of particle dark matter to the SM, there has recently been increased interest in other dark-matter candidates. In particular, much attention has been devoted to the study of ultralight classical-field bosonic dark matter, which exhibits distinct and highly varied phenomenology. The most popular candidates of this type are the QCD axions [8][9][10], axionlike particles (ALPs) [11,12], and dark photons [13][14][15]. In this paper, we focus our attention on axions.
The QCD axion was originally proposed in order to solve the strong CP problem [16][17][18]. Its mass and couplings with SM particles are defined in terms of a single parameter F a , the axion decay constant. On the other hand, many extensions to the SM, as well as generic string compactifications, imply the existence of light pseudoscalars with properties similar to the QCD ax-ion [19,20], but with the crucial difference that their mass and their SM couplings are independent parameters. These are usually called axionlike particles (ALPs). Both QCD axions and ALPs are produced in the early Universe by non-thermal mechanisms (e.g., misalignment, decays of topological defects, and others [8][9][10][21][22][23][24][25][26][27][28][29][30][31][32]) in sufficient abundance that they can potentially constitute all of the dark-matter energy density (for QCD axions this is true for F a all the way up to ∼ the Planck scale [23,31]). As these non-thermal production mechanisms create nonrelativistic particles, both QCD axions and ALPs are excellent dark-matter candidates [8][9][10]21]. In this work, we restrict our attention to ALPs, which we will from now on simply refer to as 'axions'.
The magnetic fields involved in many of these searches are produced in the laboratory by ferromagnets or solenoids carrying a strong electric current. Another possibility is to use instead the natural geomagnetic field of the Earth. 1 An important aspect of this field is that its spatial extent is much larger than the length scales that can be achieved in laboratory experiments. Axionphoton conversion depends both on the magnitude of the external magnetic field and on the spatial extent over which the interaction takes place. Having axion dark matter interacting with the Earth's magnetic field over a length scale of order the Earth's radius thus has the potential to boost axion-photon conversion to a level competitive with experiments that use stronger magnetic fields in a smaller, laboratory-scale volume.
In this paper, we exploit this observation and point out a novel signal of axion dark matter. The signal is an oscillating magnetic field generated near the surface of the Earth, and applies for axion masses m a 10 −14 eV. It is very similar to the one recently pointed out in Refs. [61,62] for kinetically mixed dark-photon dark matter, which in turn has the same conceptual origin as the signal in the DM Radio experiment [63].
Consider the axion-photon coupling term in the La-grangian, in the presence of a static background magnetic field B 0 : L ⊃ +g aγ aE · B ∼ g aγ (∂ t a)B 0 · A ∼ −im a g aγ aB 0 · A, where we used E ∼ −∂ t A (ignoring the EM scalar potential) integrated by parts dropping the boundary term, and used ∂ t a ∼ −im a a as relevant for a non-relativistic axion. This term has the same mathematical structure as the mass-mixing term that appears in a Lagrangian describing a dark photon that is kinetically coupled to the SM photon, when expressed in the interaction basis [61]: L ⊃ εm 2 A A · A ∼ −εm 2 A A · A (again ignoring the scalar potential), where m A is the dark-photon mass, and ε is the kinetic mixing parameter. The axion coupling must thus source similar physical effects to the dark-photon coupling. In particular, the well-known effect of the mass-mixing term in the darkphoton case is to create a small mixing of the relevant sterile dark-matter field with the SM electric field in such a way as to drive free-charge motion. If one considers a shielding conducting box, this results in charges being driven in the conducting walls of the shielding box, and this charge motion will in turn generate real observable electromagnetic fields: in particular, the shield forces the dominant electromagnetic field generated inside the box to be a magnetic field (assuming that the box is of a spatial extent smaller than the dark-photon Compton wavelength). In exactly the same way, if that box is permeated by a static background magnetic field and the dark matter is instead an oscillating axion field, an oscillating magnetic field is again generated inside the box. An approximate translation from the dark-photon case to the axion case is simple: εm 2 A A → ig aγ m a aB 0 , or, if we assume in each case that the relevant DM candidate is all of the DM, εm A Â → ig aγ B 0B0 . In particular, the role of the dark-photon polarization state is replaced by the background magnetic-field direction for the axion case, and the axion-induced signal amplitude is obtained from the dark-photon-induced signal amplitude by the replacement εm A → g aγ B 0 .
As for the dark-photon case considered in Ref. [61], our axion signal however does not arise from considering a human-engineered shielding box (as employed, e.g., in other axion searches, such as resonant cavity experiments [38,[64][65][66][67][68] and LC circuits [39,63,69]); instead, in our mass range of interest, Nature provides us with a ready-made shield. The near-Earth environment itself can be modeled as a conducting spherical cavity with a vacuum gap [61]: the lower layer of the atmosphere is a poor conductor sandwiched on one side by the conductive innermost layers of the Earth, and on the other side by the conductive ionosphere and/or interplanetary medium. The conductive inner-Earth, ionosphere, and interplanetary medium are thick enough to damp the electromagnetic active mode, at least in our mass range of interest. In the case of dark-photon dark matter, this is sufficient to give rise to a dark-matter induced magnetic field at the surface of the Earth [61]. Provided that m A 1/R, where R ∼ (3 × 10 −14 eV) −1 is the radius of the Earth, the electric field was suppressed by a factor (m A R) 2 because of the large natural shield. In this paper, we show that because the lower atmospheric vacuum gap is also permeated by the geomagnetic field of the Earth B 0 , a similar axion-induced magnetic field is generated if the dark matter is instead composed of axions. As in the dark-photon case, the accompanying electric field is suppressed by (m a R) 2 1, and the magnetic field oscillates at an angular frequency equal to the axion mass m a ; it also has a predictable vectorial pattern over the whole surface of the Earth (albeit one that differs from the cognate pattern for the case of darkphoton dark matter), and it inherits the temporal phasecoherence properties of the axion field.
In addition to demonstrating the existence of this novel signal of axion dark matter, we propose to search for it in a way that is conceptually identical to the approach described in Refs. [61,62], which advanced and applied this technique for dark-photon dark-matter searches: expose a geographically dispersed network of sensitive magnetometers to the ambient magnetic environment at the surface of the Earth, and record the magnetic field as a function of time over long time periods. Using distributed networks of sensors 2 in this fashion carries many advantages (see, e.g., Refs. [75,76] for some recent discussions). This search can make use of the same public database maintained by the SuperMAG Collaboration [77,78] that was previously analyzed in Refs. [61,62]. SuperMAG collates data from hundreds of unshielded three-axis magnetometers that are widely dispersed over the surface of the Earth and that have been measuring geomagnetic activity since the early 1970s with a time resolution (for the relevant dataset) of one minute.
Because this dataset has previously been analyzed for the dark-photon signal, we can easily motivate why the cognate axion search is interesting. Ignoring for the purposes of this argument that the Earth's magnetic field takes a non-trivial spatial pattern, we can employ the rough parametric signal-amplitude mapping εm A → g aγ B 0 on the existing dark-photon limits set in Ref. [61,62]. For instance, for a dark-photon mass of m A ∼ 4 × 10 −17 eV, the corresponding smoothed 95%credible upper limit on ε was found in Ref. [61,62] to be ε ∼ 1.4 × 10 −5 . The Earth's geomagnetic field ranges from B 0 ∼ 25-65 µT across its surface, 3 so the rough parametric mapping indicates that a bound in the range g aγ ∼ (4.4-11) × 10 −11 GeV −1 is potentially achievable for m a ∼ 4 × 10 −17 eV. At this mass, this esti- 2 We note that while there is a geographically distributed array of atomic magnetometers (the Global Network of Optical Magnetometers for Exotic physics searches, GNOME [70][71][72]) specifically designed to search for evidence of beyond-the-Standard-Model physics (such as couplings of axion dark-matter fields to nuclear spins [73]), the GNOME magnetometers are enclosed in meter-scale, multi-layer magnetic shields that effectively cancel the signatures [74] searched for in this work and that discussed in Refs. [61,62]. 3 Recall that 1 T ≈ 195.4 eV 2 . mate brackets the existing low-mass CAST helioscope 95%-confidence bound on the axion-photon coupling: g aγ 6.6 × 10 −11 GeV −1 [42]. This indicates that it is worthwhile to undertake this analysis carefully, accounting fully for the differing spatial patterns of the darkphoton and axion signals.
Proceeding with a careful analysis of the SuperMAG data, we find no robust evidence for a statistically significant axion-induced oscillating magnetic-field signals. Our search in the axion mass range 2 × 10 −18 eV m a 7×10 −17 eV initially identifies some 27 naïve signal candidates that appear globally significant at the 95%confidence level on the basis of our analysis pipeline. However, further robustness checks performed on these candidates cleanly eliminate the majority of them. The only candidates that are not eliminated are either in some tension with some subset of the robustness checks, or have relatively low global significance that would be insufficient to robustly claim anything more than some tension with the background-only model. Because we find no robust evidence for an axion signal, we proceed to set limits: following a Bayesian analysis procedure that accounts for stochastic fluctuations of the amplitude of the axion dark matter [79,80], we derive a posterior on the axion-photon coupling and place 95%-credible upper limits on g aγ in the same mass range as for the signal search. Assuming that the axion is all of the dark matter, our limits indeed reach the current CAST bound at their most sensitive: we set the constraint g aγ 6.5 × 10 −11 GeV −1 for 3 × 10 −17 eV m a 4 × 10 −17 eV (the limits weaken outside this range). Nevertheless, our limits have distinct systematics as compared to CAST, and future improvements using other existing archival datasets as well as via dedicated searches with new experiments are possible.
The rest of this paper is structured as follows: in Sec. II, we derive the axion dark-matter induced magnetic-field signal, beginning with a derivation of the axion effective current in Sec. II A, then discussing the IGRF-13 geomagnetic-field model in Sec. II B, giving a quick signal derivation argument in Sec. II C with many details deferred to the appendices, and then comparing it to the cognate dark-photon signal [61] in Sec. II D. In Sec. III, we summarize our signal search at a high level, again deferring details to the appendices, and present a set of axion-photon coupling exclusion bounds in Fig. 1. We conclude in Sec. IV. There are a number of appendices that expand on the main text with further detail: Appendix A gives our conventions for the vector spherical harmonics and a number of relevant identities that we utilize in this work. In Appendix B, we present a thorough and detailed derivation of the signal under two different sets of assumptions regarding the modeling of the near-Earth conductivity environment; this reinforces the quicker derivation given in the main text. The details of our analysis are presented in Appendix C; we mirror the discussion in Ref. [62] and give a minute accounting of all relevant differences. Included in Appendix C 4 is a detailed investigation of some anomalies ('naïve signal candidates') that we identified in the data, but which we do not consider to be strong, robust signals of axion dark matter for reasons also discussed in that appendix.

II. SIGNAL
In this section, we describe the observable magneticfield signal sourced by axion dark matter at the surface of the Earth. The signal described in this work is analogous to the one described in Ref. [61], but with axion dark matter replacing the role of dark-photon dark matter.
It can be shown that the Earth itself acts as an effective conducting shield for electromagnetic waves of frequencies 10 −21 eV ω 3 × 10 −14 eV (see Sec. II B of Ref. [61]). At these frequencies, the lower atmosphere just above the surface of the Earth, however, has negligible damping effects on electromagnetic waves [61]. Even further from the surface of the Earth, either the ionosphere or interplanetary medium act as an effective shield in this frequency range [61] owing to their large conductivity or plasma frequency, respectively. 4 The near-Earth environment can thus be treated as an inner conducting sphere (the Earth) and a surrounding conducting shield (the ionosphere/interplanetary medium), separated by a vacuum region (the lower atmosphere).
In Ref. [61], the effect of kinetically mixed dark-photon dark matter in this environment was parametrized by an effective background current, oriented in the direction of the local dark-photon field. It was shown that the effect of this current would be to generate an oscillating magnetic-field signal at the surface of the Earth which exhibits a particular global spatial pattern, given by Φ m vector spherical harmonics (VSH) [see Appendix A for VSH conventions]. Importantly, the leading Φ m contributions did not depend on the details of the conducting boundaries.
The derivation in this work of the signal of axion dark matter in this near-Earth conductivity environment proceeds similarly to the derivation for dark-photon dark matter in Ref. [61], but with some important differences. Specifically, we can also employ an effective current approach for electromagnetically coupled axion dark matter. One crucial difference between the axion and dark-photon cases however is that the axion requires the Earth's static magnetic field in order to convert into an observable electromagnetic signal. The axion-induced effective current will thus not only depend on the local axion field value but also on the local geomagnetic field.
In particular, the effective current will inherit its direction from the Earth's quasi-static magnetic field. Once translated into the language of an effective current, the axion calculation proceeds similarly to the dark-photon calculation. We thus defer a detailed calculation of the axion dark-matter signal to Appendix B, and instead in this section present a simpler argument (which could also apply to the dark-photon case).
We begin this section with a description of the effective current for axion dark matter in the presence of a magnetic field. Next, we describe the International Geomagnetic Reference Field (IGRF) model for the Earth's static magnetic field. Then we present a simple derivation for the Φ m component of the axion dark-matter signal. Finally, we conclude by comparing the properties of this axion signal with the previously described darkphoton signal.

A. Effective current
In this work, we consider an axion a coupled to electromagnetism with strength g aγ , described by the Lagrangian In the presence of such an axion, Maxwell's equations are modified as [33,81] ∇ · E = g aγ ∇a · B, If the axion a comprises the dark matter, then it is nonrelativistic and so ∂ t a ∼ −im a a and |∇a| ∼ m a v dm a, where v dm ∼ 10 −3 . Therefore, all ∇a terms will be parametrically suppressed compared to the ∂ t a terms, and the leading effect of the axion dark matter will thus come from the −g aγ (∂ t a)B term in Eq. (5). We can see from the form of Eq. (5) that this term behaves similarly to a current in the usual Ampère-Maxwell law, In the presence of a static background magnetic field B 0 , this effective current is given by Note that in the case of dark-photon dark matter, the direction of J eff is set by the direction of the dark photon polarization [61] and is thus spatially uniform. By contrast, in the axion case, the direction of J eff is not determined by any property of the axion, but rather by the direction of the static magnetic field B 0 . For the derivation of our signal, B 0 will be the Earth's DC magnetic field. Thus, the effective current in the axion case will not be uniform in space, but will have an (approximately) dipolar angular dependence and decay with radial distance from the Earth's center.

B. IGRF model
In this work, we employ the IGRF-13 model [82] of the Earth's magnetic field, which provides coefficients for a multipole expansion of the field. As the geomagnetic field drifts slowly over time, the IGRF model provides coefficients for the field at five-year intervals and specifies an interpolation procedure on these coefficients to obtain the field at intermediate times. The most recent generation, IGRF-13, provides values dating from 1900 up to 2020. The IGRF model parametrizes the geomagnetic field in terms of a scalar potential 5 B 0 = −∇V 0 , which is then expanded as where R = 6371.2 km is the exact reference value for the Earth's radius that is used in the specification of the IGRF model [82], θ and φ are geographic co-latitude and longitude, 6 P m are the Schmidt-normalized associated Legendre polynomials, 7 and Y m are the scalar spherical 5 Recall that in the magnetoquasistatic limit and in the absence of free currents, the Ampère-Maxwell law becomes ∇ × B = 0. We can therefore define a scalar potential V so that B = −∇V . See Sec. 5.9.B of Ref. [83] for a more detailed discussion. 6 Note that the coordinates θ and φ co-rotate with the Earth; that is, they describe a rotating frame with coordinates fixed to the Earth, not the inertial frame with coordinates fixed to the average positions of distant stars. Henceforth, all VSH and spherical harmonics will implicitly use coordinates in this corotating frame as well. As the rotational speed at the surface of the Earth is non-relativistic, there are no relativistic fieldmixing effects which need to be taken into account when switching between these frames. Thus, it remains consistent to apply Maxwell's equations as in Eqs.
(2)-(5) even in the co-rotating frame. Specifically, because both the Earth's magnetic field and the observation points for the magnetic observatories co-rotate with the Earth, this is the natural co-ordinate frame to use. 7 These are given for m ≥ 0 by [84] where P are the Legendre polynomials. Note that the Schmidtnormalized P m used in the IGRF specification differ from those defined at Eq. (3.49) in Ref. [83] in several ways; however, our scalar Y m are normalized to agree with those of Ref. [83].
harmonics. The IGRF model provides the 'Gauss coefficients' g m and h m for ≥ m ≥ 0 at five-year intervals (see Tab. 2 of Ref. [82]), and their values at intermediate times are to be calculated by linear interpolation.
Here we adopt the conventions g ,−m = (−1) m g m and h ,−m = (−1) m+1 h m to extend the coefficients to negative m. From B 0 = −∇V 0 , we can then write a multipole expansion for the geomagnetic field in terms of the VSH as where C m are related to the Gauss coefficients by Note that our phase conventions for g m and h m ensure that C ,−m = (−1) m C * m , in analogy to the VSH phase conventions Eqs. (A2)-(A4). As the Earth's magnetic field is approximately dipolar, with the dipole axis oriented relatively close to the Earth's rotational axis, the largest of these coefficients will be C 10 ; however, subsequent terms can provide corrections of O(10%). The IGRF-13 model provides values up to = 13 for the most recent coefficients. In our analysis (summarized in Sec. III and detailed in Appendix C), we find it sufficient to utilize the Gauss coefficients up to = 4. We have verified explicitly that the addition of higher-modes has no significant impact on our results.

C. Signal derivation
Here we give a simple derivation for the leading order Φ m contribution to the magnetic-field signal of axion dark matter at the Earth's surface (see Appendix B for a more detailed calculation). This derivation relies only on the effective current approach, and so a similar derivation can also be applied to the dark-photon case computed in Ref. [61]. As in Ref. [61], we model the near-Earth environment as a perfectly conducting sphere of radius R (the Earth), surrounded by some vacuum region (the lower atmosphere), which is further surrounded by some perfectly conducting boundary. Here we do not assume any particular shape for the outer boundary, only that it has a longest length scale L m −1 a . For this reason, our model allows for the outer boundary to be either the ionosphere, which is approximately spherical and located only ∼ 100 km from the Earth's surface, or the magnetopause (beyond which lies the interplanetary medium), which is highly aspherical and extends to ∼ 200R from the Earth's surface in the 'downwind' direction of the solar wind. 8 As this longest length scale L is smaller than the de Broglie wavelength 9 of the axion dark matter λ dB ∼ (m a v DM ) −1 , we can take the axion field value to be constant over the entire geometry, and write it as Due to the stochasticity of the axion field [79,85], a 0 is not uniquely determined by the DM density, although it is generically of order |a 0 | ∼ √ 2ρ dm /m a (see discussion at the end of this subsection).
From Eqs. (7) and (11), the effective current which this axion dark matter sources is then given by Now we argue that the ∂ t E term in Eq. (6) can be neglected. 10 This is because E vanishes both deep within the Earth and within a skin depth of the outer boundary (as they are both good enough conducting shields to effectively damp all electromagnetic waves). Moreover, as our geometry has longest length scale L m −1 a , these two surfaces on which E vanishes are separated by a subwavelength distance. We therefore only expect that E can grow quadratically in m a L between them. In particular, we expect parametrically E ∼ (g aγ a 0 )(m a L) 2 C m . Comparing to Eq. (14), we see that ∂ t E is parametrically smaller than J eff . Therefore, up to corrections at order O((m a L) 2 ), it suffices to only consider the first and last terms of Eq. (6).
Given the form of J eff in Eq. (14), we can apply the VSH curl properties Eqs. (A13)-(A15) to solve Eq. (6). Namely, we find that B must be of the form damping effects of the ionosphere are uncertain, yet the furthest point of the magnetopause may not be a sub-wavelength distance from the Earth's surface; i.e., 200R ∼ m −1 a . The validity of our assumptions are therefore questionable around ma ∼ (few) × 10 −16 eV; see also the discussion in Sec. II B (and in particular Sec. II B 6) of Ref. [61]. This mass lies outside of the range explicitly constrained in our analysis (see Sec. III). 9 In principle, there are two length scales of the dark matter which could be relevant here: the de Broglie wavelength λ dB ∼ (mav rel ) −1 and the coherence length λ coh ∼ (ma∆v) −1 , where v rel is the mean relative velocity between the dark-matter rest frame and the Earth, and ∆v is the magnitude of the local DM velocity dispersion. In order for Eq. (13) to remain valid, L must be smaller than both of these length scales. We note however that for virialized dark matter following a Maxwell-Boltzmann velocity distribution, v rel ∼ ∆v ∼ 10 −3 in the Earth's rest frame, so that λ dB ∼ λ coh . 10 See the end of Sec. III C in Ref. [61] for a similar discussion.
where V is some scalar function (so that ∇V is curl-free). From the spherical-harmonic gradient relation Eq. (A9), we can also note that ∇V consists entirely of Y m and Ψ m modes. Therefore, the leading order Φ m contribution to the magnetic field is precisely given by the first line of Eq. (15). In particular, at the surface of the Earth (r = R), to leading order in m a L, the Φ m contribution to the magnetic field signal of axion dark matter is 11 Note that due to the VSH orthogonality properties Eqs. (A6)-(A8), any vectorial function on the sphere can be decomposed into VSH (much like any scalar function on the sphere can be decomposed into scalar spherical harmonics). Thus, when searching for our signal in global magnetic-field data across the Earth, we can project onto the particular combination of Φ m modes appearing in Eq. (16). This allows us to neglect the Y m and Ψ m contributions coming from ∇V , which may generically depend on the shape of the outer boundary, and instead focus on the Φ m contributions which we know to be present regardless of details of the outer boundary. Finally, we comment on the temporal coherence of our signal in Eq. (16). The monochromatic description of the axion given in Eq. (13) remains valid only on timescales less than the coherence time T coh ∼ 2π/(m a v 2 dm ) of the axion. For the mass range relevant to our analysis (summarized in Sec. III and detailed in Appendix C), we have T coh ∼ 2-50 yr. On timescales longer than the coherence time, a 0 will vary stochastically in both amplitude and phase. Therefore, Eq. (16) only remains valid for times t T coh : the magnetic-field signal's phase offset and amplitude randomize on longer timescales, with the phase offset within each coherence time being uniformly distributed on [0, 2π), and the amplitude being set by a 0 drawn from a distribution [79,85] and satisfying |a 0 | 2 τ = 2ρ dm /m 2 a , on average over timescales τ T coh . See also the more detailed discussion of this point in the context of dark-photon dark matter in Ref. [61].

D. Comparison with dark-photon signal
The signal described by Eq. (16) takes a very similar form to the signal described in Ref. [61]. In particular, if the Earth's magnetic field is assumed to be exactly dipolar (C m = 0 for > 1), then Eq. (16) takes precisely the same form as the signal from a dark photon expressed in inertial coordinates, with the role of the dark-photon polarization in setting the signal orientation replaced by the Earth's magnetic dipole (cf. Eq. (38) of Ref. [61]). Here we highlight three important differences that allow an axion dark-matter signal to be distinguished from a dark-photon signal.
The first is due to the fact that the geomagnetic field is not exactly dipolar, and so > 1 modes will contribute to Eq. (16), giving it a slightly different angular dependence than a dark-photon signal. As mentioned before, this correction in angular dependence will be at the level of O(10%).
Secondly, the dark-photon signal in co-rotating coordinates receives a shift in frequency by f d = (sidereal day) −1 due to the rotation of the Earth (see Eq. (42) of Ref. [61]), whereas the axion dark-matter signal does not. This is because the effective current (and thus the angular dependence of the magnetic-field signal) inherits its direction from the geomagnetic field in the axion case but from the dark-photon field itself in the dark-photon case. Because the geomagnetic field corotates with the Earth, the angular dependence of the axion signal is constant in geographic coordinates (on timescales short enough that the geomagnetic field does not drift significantly). On the other hand, since the dark-photon direction is fixed in inertial coordinates (on timescales shorter than the coherence time of the darkphoton field), the dark-photon signal precesses in geographic coordinates. Thus, in frequency space, the axion signal (as measured by magnetometers fixed on the Earth's surface) only appears at f = f a ≡ m a /2π (assuming v dm = 0; see footnote 34 in Ref. [62]), while the dark-photon signal also exhibits sidebands at f = f a ±f d .
Finally, the stochastic properties of the axion darkmatter signal could differ from those of the dark-photon signal. Classical-field dark-matter candidates (both axions and dark photons) are comprised of a sum (really, an integral) over constituent Fourier modes, each of which has a random phase and, for the dark-photon case, a vectorial orientation (which is in general complex). As a result, classical-field dark-matter exhibits amplitude and overall phase-offset fluctuations from one coherence time to the next; see, e.g., Refs. [79,85]. In addition, for the dark-photon case, there can be a fluctuation of the polarization state of the field, but this depends on the assumed underlying structure of the individual Fourier modes' vectorial orientations. Depending on the formation model and subsequent cosmological evolution of the dark matter, it is an open question (see, e.g., Ref. [86]) whether these individual Fourier modes' vectorial orientations are all the same, or whether they are effectively random. In the former case, the dark-photon polarization state does not randomize from one coherence time to the next; in the latter case, it does. 12 Since the direction of the effective current that gives rise to the dark-photon signal is set by the dark-photon polarization state, the polarization-state fluctuation that arises in this latter case results in a fluctuation of the relative phases appearing between the different components of the effective current from one coherence time to the next. Because the magnetic-field signal is determined by the effective current, this results in an O(1) change to the global angular dependence of the darkphoton-induced magnetic-field signal from one coherence time to the next in this case.
By contrast, for axion dark matter, the direction of the effective current is determined by the geomagnetic field, so no fluctuation in the relative phases of the components of the effective current appears from one coherence time to the next. The direction of the effective current and the global angular dependence of the magnetic-field signal instead drift only as the geomagnetic field drifts on very long timescales (which is independent of the coherence properties of the axion).
Over the 50-year duration of the SuperMAG dataset utilized in our analysis (summarized in Sec. III and detailed in Appendix C), the direction of the geomagnetic field drifts O(5%). Thus, the change in global angular dependence of the axion dark-matter signal is significantly smaller over the duration of our analysis than the O(1) change that the dark-photon signal experiences over timescales of order the coherence time, in the case assumed above. Moreover, the Earth's background magnetic-field drift is modeled, and we include this effect in our analysis of the axion signal, whereas the possible coherence-time to coherence-time drift in the global angular dependence of the signal in the dark-photon case is inherently stochastic, although accounted for in Ref. [62].

III. SEARCH FOR SIGNAL IN SUPERMAG DATASET
The axion dark-matter signal described by Eq. (16) is an oscillating magnetic field at the surface of the Earth of magnitude assuming the axion constitutes all of the dark matter (we take ρ dm = 0.3 GeV/cm 3 ). It is temporally coherent over a long time period; it is also spatially coherent, taking a known pattern across the entire globe. As such, ter case. However, because the network of magnetometers contributing to the SuperMAG dataset has reasonably isotropic directional sensitivity (being a network of O(500) three-axis magnetometers that are widely distributed on the rotating and orbiting Earth), we expect that even if the dark-photon dark matter behaved according to the former case, our limits in Refs. [61,62] would be changed by only an O(1) factor. it could be detected by any global array of unshielded magnetometers taking data over several decades. The SuperMAG Collaboration [77,78,87] maintains a public database of measurements from precisely such an array of magnetometers. In particular, they report three-axis magnetic field measurements from O(500) stations with a one-minute time resolution, with the measurements from some stations dating back to 1970.
In this work, we perform a search of the SuperMAG dataset 13 for the signal described by Eq. (16), similar to the dark-photon dark-matter search undertaken in Ref. [62]. As the analysis proceeds in a similar way to the one in Ref. [62], we reserve the details of this work's analysis (and, in particular, how they differ from those of the analysis in Ref. [62]) to Appendix C. Instead, in this section, we summarize the results of this work's axion dark-matter search. We first enumerate some naïve candidate signals which we identified but through further robustness checks dismissed. Having dismissed all such naïve signal candidates, we then present an exclusion limit on axion dark-matter parameter space. Finally, we discuss a higher-resolution dataset also maintained by the SuperMAG collaboration, a future analysis of which could extend the results of this work to new parameter space.
As in Ref. [62], we performed a search for our signal at O(10 6 ) discrete frequencies in the frequency range 6 × 10 −4 Hz f a 2 × 10 −2 Hz, corresponding to the mass range 2 × 10 −18 eV m a 7 × 10 −17 eV. For each frequency, we constructed analysis variables (see Appendix C 1) and wrote down a likelihood function for the photon-axion coupling g aγ , given the observed values of these variables (see Appendix C 3). Using this likelihood, we determined whether the analysis variables at each frequency were consistent with the lack of a signal g aγ = 0. We declared any frequency at which the data were inconsistent with g aγ = 0 at 95% confidence global significance to be a 'naïve signal candidate' for axion dark matter.
Based on this initial analysis, we identified 27 such candidates (some of which can be seen as narrow peaks above the dark blue exclusion band in Fig. 1). We then reevaluated each naïve signal candidate for robustness, to test if it exhibits features of a physical axion dark-matter signal. In particular, a physical axion dark-matter signal should be present for the entire duration of time over which SuperMAG has collected data, and should appear in the data from all stations across the globe. 14 Because of this, we partitioned the full SuperMAG dataset into 13 In this work, we use a slightly updated dataset as compared to the one used in Ref. [62], which includes a few additional stations.
The temporal duration of the dataset remains the beginning of 1970 through the end of 2019. As this updated dataset contains minimal additional data, we do not expect to gain significant sensitivity. 14 The manner in which our main analysis proceeds does attempt to project onto this global mode, but the exact fashion in which four temporal subsets, consisting of data only from certain time periods, and four geographical subsets, consisting of data only from certain stations. We re-performed our analysis on each of these subsets and searched for each naïve signal candidate to see if it re-appeared in the subset. In particular, we checked if the analysis variables constructed from each data subset were consistent with the signal size implied by the original analysis of the full SuperMAG dataset, as characterized by the Bayesian posterior on g aγ constructed from the latter (see Appendix C 3). We combined the results of all eight of these resampling checks and immediately rejected any candidates with a combined p-value of p full < 0.01 (see Appendix C 4). We note three candidates of remaining potential interest that we were unable to immediately reject based on this criterion: (1) one candidate at f a ≈ 4.2 mHz with a high (6.7σ) global significance that is in strong but not definitive tension (0.01 < p full < 0.05) with the combined spatio-temporal robustness tests and also the geographical tests alone (0.01 < p geo < 0.05); (2) one candidate at f a ≈ 5.5 mHz with 3.6σ global significance that is however in strong tension with the temporal resampling checks (0.01 < p time < 0.05); and (3) one candidate (actually a pair symmetrically arranged around the Nyquist frequency) at f a ≈ 8.3 mHz with a low (2.3σ) global significance that is nevertheless consistent with all resampling checks. While we do not consider these candidates to constitute strong and robust evidence for dark matter, they would require further work to definitively exclude.
With no robust axion dark-matter candidates identified, we set 95%-credible exclusion limits (local significance) on the axion-photon coupling g aγ based on the Bayesian posterior on g aγ derived in our analysis. Fig. 1 shows our exclusion limit on g aγ as a function of m a , assuming that the axion is all of the dark matter and that ρ dm = 0.3 GeV/cm 3 . Also shown in Fig. 1 are limits on g aγ set by the CAST solar axion search [42], and limits based on non-observation of gamma-rays in coincidence with supernova SN1987A by the Gamma-Ray Spectrometer instrument on the Solar Maximum Mission satellite [88]. The latter limit arises because axions produced via the Primakoff process in SN1987A would convert to gamma rays in the Milky Way's magnetic field. 15 We do note that neither the CAST nor the SN1987A limits must assume that the axion is all of the dark matter. Nevertheless, from Fig. 1, it can be seen that our limits are competitive with CAST bounds in some mass ranges; this is done leaves the main analysis vulnerable to a false-positive signal identification if a small subset of the stations exhibit a very large signal oriented in the appropriate direction. The point of the re-evaluation/robustness tests is to exclude this vulnerability. 15 It was recently pointed out [89] that accounting for the turbulent component of the Milky Way magnetic field complicates the computation of the conversion of astrophysical axions into gamma rays. Accounting for this turbulent component can change the axion-photon conversion probability by up to a factor of two.  16), as summarized in Sec. III and detailed in Appendix C. Our exclusion limit is shown as a function of the axion mass ma in solid dark blue and appears as a wide band due to the density of masses at which limits are plotted. The solid light blue line shows the sliding average of our exclusion limit over nearby frequencies. Our exclusion limit exhibits several narrow spikes (each at most a few frequency bins wide), which correspond to potential signal candidates. We investigate these candidates further in Appendix C 4 and show that none constitute robust evidence for dark matter. Also shown are existing limits from the CAST helioscope search for axions produced in the Sun [42] (dashed orange), and a constraint due to non-observation of a gamma-ray signal from axions in coincidence with SN1987A [88] (dotted green).
they also rely on independent systematics.
A comment is in order on the mass dependence of our limits as shown in Fig. 1: because our magneticfield signal has B ∝ m a a 0 and because the root-meansquare (rms) value of a 0 ∝ √ ρ dm /m a , the rms amplitude of the magnetic field signal is independent of the axion mass once the rms axion field amplitude is normalized to the dark-matter density. The m a -dependence of our limits is thus driven by the underlying noise behavior of the SuperMAG dataset as a function of frequency [f a = m a /(2π)]. This contrasts with the dark-photon case [61,62], for which the signal amplitude itself also still depended linearly on the dark-photon mass, even after normalizing the dark-photon field amplitude to the dark-matter abundance (a well-known decoupling effect of the massless dark-photon limit).
In addition to the one-minute resolution dataset analyzed in this work, SuperMAG also maintains a onesecond resolution dataset from a smaller number of stations. A similar axion dark-matter search in this higher time-resolution dataset would allow for sensitivity to higher axion masses. Given that our limit shown in Fig. 1 improves with increasing mass (owing to lower noise at higher frequencies in the one-minute SuperMAG dataset), we anticipate that the constraint from a search in the one-second resolution SuperMAG dataset could potentially outperform existing constraints from CAST in some mass range (and perhaps the SN1987A constraint, although this is less clear), assuming the noise in that dataset continues to behave similarly. Addition-ally, a search in the higher time-resolution dataset would present an opportunity to re-evaluate the three naïve signal candidates of interest that we discussed above to see if they appear in that dataset. We intend to undertake such a search in future work.

IV. CONCLUSION
In this paper, we described a novel signature of ultralight axion dark matter with a coupling to photons g aγ . This signal is similar to that of dark-photon dark matter that was recently discussed in Ref. [61] and searched for as described in Ref. [62]. Namely, we pointed out that such an axion field converts off the quasi-static geomagnetic field of the Earth B 0 , to produce at ground-level all across the Earth's surface an observable magnetic-field signal B a (t). This signal oscillates at a frequency f a set by the axion mass m a ≈ 2πf a , a fundamental physics parameter; moreover, it is narrowband in the sense that the bandwidth ∆f ∼ v 2 dm f a ∼ 10 −6 f a , implying a long phase-coherence time for these oscillations. The signal amplitude is |B a | ∼ g aγ |B 0 |R √ ρ dm where R is the radius of the Earth, implying that it is detectably large (note the appearance of R here, and not some other length scale, such as the height of the atmosphere). Finally, the signal has a global vectorial pattern that is set by the Earth's quasi-static geomagnetic field. As such, this signal is an ideal candidate to be searched for using a network of globally distributed, terrestrial magnetic-field metrology stations. There is an existing publicly available dataset of measurements of this type maintained by the SuperMAG Collaboration [77,78], which consists of 50 years' worth of one-minuteresolution, three-axis magnetometer readings taken at 508 geographically dispersed stations in total (although not all stations report data at all times). We made use of this dataset to search for our axion-induced magneticfield signal in the axion mass range 2 × 10 −18 eV m a 7 × 10 −17 eV, constructing our analysis around projections of this large dataset onto a small number of vector spherical harmonic coefficients in which our signal is expected to appear. Our search initially identified 27 naïve signal candidates in the data that at a global 95% confidence level were inconsistent with a background-only hypothesis. However, applying further robustness checks to test these candidates for spatial consistency and temporal uniformity, we definitively eliminated all but three of them. The three candidates that were not definitively eliminated however still exhibited strong tension with (at least some subset of) our robustness tests, or had weak global significance. As such, we do not consider any of them to be strong and robust signals of axion dark matter on the basis of this analysis.
Having dismissed all the anomalies in the data as either relatively weak and/or exhibiting of some defect and thus not robust, we turned to setting limits on the axionphoton coupling g aγ . We made use of a Bayesian analysis procedure that folded in the effects of the stochastic fluctuations of the axion dark-matter field amplitude from one coherence time to the next. Assuming that the axion is all of the DM, we set 95%-credible upper limits (local significance) on the axion-photon coupling g aγ as a function of the axion mass m a in the same mass-range as our signal search; see Fig. 1. These limits are strongest for 3 × 10 −18 eV m a 4 × 10 −18 eV: smoothed over frequency-to-frequency fluctuations, the mean limit in this mass-range reaches g aγ 6.5 × 10 −11 GeV −1 , which is comparable to limits on axions set by the CAST helioscope [42]; they are however about an order of magnitude weaker than astrophysical limits set by observations of SN1987A [88] (but see also Ref. [89]).
The already impressive reach for this search, which has dramatically different systematics as compared to the other constraints in this axion mass range and is thus both competitive and complementary, could be further improved by future analysis of a higher resolution (onesecond) dataset also maintained by the SuperMAG Collaboration. Although this dataset has data from fewer stations and over a shorter total temporal duration as compared to the data analyzed in this work, if the decrease in the noise moving to higher frequencies that is evident in Fig. 1 persists also in that other dataset, it is possible that such an analysis could further probe for signals below existing CAST bounds at frequencies up to a factor of 60 higher than those searched in the present work. Moreover, this would afford the opportunity to revisit some of the weak or non-robust anomalies observed in this work for further analysis. In future planned work that will be undertaken in collaboration with members of the SuperMAG Collaboration, we will apply the analysis techniques we have developed in this work and in Refs. [61,62] to this one-second resolution SuperMAG dataset. Some of the computing for this project was performed on the Sherlock cluster. We thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results.
The work of M.A.F. was performed in part at the Aspen Center for Physics, which is supported by NSF Grant No. PHY-1607611.
We gratefully acknowledge the SuperMAG Collaboration for maintaining and providing the database of ground magnetometer data that were analyzed in this work, and we thank Jesper W. Gjerloev for helpful correspondence regarding technical aspects of the Super-MAG data. SuperMAG receives funding from NSF Grant Nos. ATM-0646323 and AGS-1003580, and NASA Grant No. NNX08AM32G S03.
We acknowledge those who contributed data to the Su- This appendix, which defines the VSH conventions used in this work, is reproduced from Ref. [61] with minor modifications for the convenience of the reader.
The VSH are defined in terms of the scalar spherical harmonics Y m by the relations, wherer is the radial unit vector. Thus Y m points radially, while Ψ m and Φ m point tangentially. Some of their relevant properties (and our phase conventions) are For any radially dependent function f (r), the gradient of the scalar spherical harmonics can be related to the VSH by Additionally, the divergences and curls of the VSH are given by with the Laplacians then being In this appendix, we derive in more detail the axion dark-matter induced magnetic-field signal, Eq. (16), that would be measured in the lower-atmospheric air gap just above the surface of the Earth. The calculation in this appendix closely follows the calculations in Secs. III B and III C of Ref. [61], but using the effective current given by Eq. (14). As in Ref. [61], we will derive the axion dark-matter signal using two different models of the near-Earth conductivity environment. In the first, we will take the inner and outer boundaries of the loweratmospheric air gap to be spherical perfectly conducting shells at r = R (corresponding to the Earth) and r = R+L (corresponding to the ionosphere), respectively, with the lower-atmospheric region separating them assumed to be vacuum. In this case, the magnetic-field signal at r = R, to leading order in m a R, will be exactly Eq. (16). In the second model, we will still take the inner boundary of the lower-atmospheric air-gap to be a spherical Earth (this is accurate to 0.3% [91]), but we will allow the outer boundary to have an arbitrary shape; this corresponds to the scenario where the outer boundary is the interplanetary medium. In this case, Eq. (16) will give only the leading Φ m contribution to the magnetic-field signal at r = R; in general, other Y m and Ψ m contributions may also be present, but these contributions may depend on the details of the outer boundary.
In either model, we can decompose the electric field in the vacuum region into two contributions (cf. Eqs. (14)-(16) of Ref. [61]), where E inh is chosen to satisfy and E hom is chosen so that E fulfills the boundary conditions on the full solution, while satisfying Both contributions must also satisfy Because of this last criterion, both contributions must be composed of 'transverse electric' (TE) and 'transverse magnetic' (TM) modes, whose electric fields are of the form for some scalar functions f m and g m to be determined below. We denote the TE and TM contributions to E inh by E inh,TE and E inh,TM , with relevant scalar functions f inh, m and g inh, m (and likewise for E hom ).
Comparing Eq. (14) with the forms of the RHS of Eqs. (B5) and (B6), it can be shown straightforwardly that J eff has the same form as the RHS of Eq. (B6), with where x ≡ m a r and x 0 ≡ m a R. Substituting Eq. (B6) into Eq. (B2) and making use of the VSH Laplacian properties Eqs. (A16)-(A18), it follows that These equations, one for each ( , m) pair, are solved by 16 Likewise it can be readily seen that Eq. (B3) implies so that g hom, m are linear combinations of spherical Bessel functions 16 This solution can be read off by inspection after noting that where we have extracted a factor of x 2 +1 0 from B m so that A m and B m have the same power-counting in an expansion in x 0 .

Spherical boundary conditions
Now to solve for A m and B m and obtain the full solution, we must fix boundary conditions. Let us first consider the near-Earth conductivity model with spherical boundaries, where the inner boundary lies at r = R and the outer boundary lies at r = R + L. This means that the Ψ m components of the full electric field must vanish at these radii (as they are parallel to the boundaries). This implies the boundary conditions, (B13) Let us assume that the axion Compton wavelength is much larger than the largest physical length-scales in the problem: m a R, m a (R + L) 1. We can then employ the small-x expansion of the spherical Bessel functions where n!! = n · (n − 2) · (n − 4) · · · k n , where k n is the smallest positive integer with the same parity as n. The leading terms in Eqs. (B12)-(B13) in this expansion are (B17) This is solved by A m = 0 and with corrections appearing at higher orders in m a R and m a (R+L); note that this expansion has not assumed anything as to the relative sizes of R and L. Interestingly, this is in a sense the opposite outcome to that of the dark-photon calculation in Ref. [61], where the B m were the coefficients that vanished to leading order. As in the dark-photon case, all contributions to the full electric field vanish everywhere in the cavity to order O(x 2 0 ), where at the last equality, we used the expansion Eq. (B15) in the ( · · · ) bracket, which causes the leading (x 0 /x) l+2 term to vanish leaving a leading correction at O(x 2 0 (x 0 /x) l ); since x ≥ x 0 , the ( · · · ) bracket is parametrically O(x 2 0 ) too. The magnetic-field contributions, on the other hand, can be derived by Faraday's law, Eq. (4). Because the inhomogeneous contributions have no associated magnetic field the full magnetic field comes solely from the homogeneous contributions: which agrees with Eq. (16).

Aspherical boundary conditions
Now consider the case where the boundaries have arbitrary shape (but we are still interested in the magnetic field at r = R). As shown in Eq. (B21), the electric field in the above solution vanishes everywhere in the cavity in all directions up to corrections at O[(m a R) 2 ]. This means that one could perturb the boundary shape in any arbitrary fashion, so long as the largest length-scale associated with boundary remains smaller than the Compton wavelength of the axion, and the correct electric-field boundary condition would still be satisfied on that perturbed boundary up to corrections at O[(m a R) 2 ]. Therefore, regardless of the boundary conditions, the above solution will remain the correct solution for the homogeneous electric field up to corrections at order O[(m a R) 2 ]. It then only remains to determine what this implies for the magnetic-field solution.
There is a power-counting argument for the magnetic field solution that can be followed and which proceeds similarly to the argument in Sec. III C of Ref. [61]. Namely, if we Taylor expand E hom,TE in x 0 = m a R as for some scalar functions f m , where ∇ ξ denotes that the radial co-ordinate in the derivative is taken to be ξ. Importantly, B where Thus the full homogeneous electric and magnetic fields can be expanded as where As noted above, the homogeneous electric-field solution in the non-spherical case will be precisely the same as in the spherical case up to additional corrections at O(x 2 0 ). Therefore, E hom is as in the spherical case, and E and the leading order magnetic field will be hom,TE + B hom,TM , While it is clear from Eq. (B38) that B hom,TE can give a leading-order contribution to the magnetic field, it will contribute to different VSH than B (1) hom,TM . Namely TE contributions to the magnetic field are comprised of Y m and Ψ m modes, while TM contributions are comprised of Φ m modes. Therefore, Eq. (16) [or Eq. (B24)] indeed gives the correct leading order Φ m contributions to the axion dark-matter magnetic-field signal, regardless of the boundary shape (so long as the largest physical scale in the problem is shorter than the axion Compton wavelength).

Appendix C: Analysis details
In this appendix, we explain the details of the signal search whose results are summarized in Sec. III. The analysis in this work proceeds similarly to the dark-photon analysis described in Ref. [62]. Here we therefore focus on the aspects of the axion analysis described in this work that differ from the dark-photon analysis, and we refer the reader to Ref. [62] for the parts of the analysis that are identical.
The analysis described in this appendix searches global magnetometer data maintained by the SuperMAG Collaboration [77,78,87] for the signal Eq. (16). The Super-MAG dataset that we analyze consists of time series of three-axis magnetic-field measurements from each of 508 stations which together cover a 50-year-long time period from the beginning of 1970 to the end of 2019. Since many of the stations began reporting data later than the beginning of 1970, shut down prior to the end of 2019, or underwent periods of inactivity, the time series that each station reports is not continuous over the entire 50-year duration 17 of the SuperMAG dataset. We denote 17 For technical reasons related to the number of independent stations' measurements required to perform the analysis for a vectorial dark-photon dark-matter signal, the analysis in Ref. [62] restricted its attention to the data taken from the beginning of 1972 to the end of 2019; see footnotes 7 and 39 of Ref. [62]. These reasons do not apply to the analysis in this work, as only one active station is required to produce the two linearly independent time series X θ and X φ required to perform the analysis for the scalar axion dark-matter search. Therefore, the present analysis also utilizes the data from the years 1970 and 1971.
the geographic coordinates of station i by Ω i = (θ i , φ i ) and the three-axis magnetic-field measurement it reports 18 We also denote the set of sampling times at which station i reports valid measurements by T i . Importantly, T i differs between stations, making a straightforward station-bystation analysis difficult.
Our search for an axion signal at frequency f a = m a /2π (assuming v dm = 0; see footnote 34 in Ref. [62]) proceeds roughly as follows. We combine the B θ i (t j ) measurements from all 508 stations into one time series X θ (t j ), based on theθ-component of the signal in Eq. (16), and we likewise combine the B φ i (t j ) measurements into another time series X φ (t j ) based on theφcomponent of the signal. We then partition each time series into segments X θ k (t j ) and X φ k (t j ) of duration roughly equal to the coherence time T coh of the axion dark-matter field (which depends on f a ). We Fourier transform each segment to findX θ k (f a ) andX φ k (f a ), and combine them into a two-dimensional vector X k . These X k are the primary variables that we use below to construct the likelihood function for our analysis. As such, we must compute both their expectation values under the signal hypothesis and variances under the background-only hypothesis. Their expectation values under the signal hypothesis are computed by performing the same time-series combination on the signal Eq. (16) as we performed on the SuperMAG data. Their variances under the background-only hypothesis are computed by a data-driven noise estimation procedure, identical to the one outlined in Ref. [62].
With the statistics of X k computed, we construct a likelihood function under the assumption of a signal with coupling g aγ . Assuming an objective Jeffreys prior on g aγ , we use this likelihood in a Bayesian analysis framework to compute the posterior for g aγ . Finally, this posterior is used to set 95%-credible upper limits on g aγ .
Of course, setting exclusion limits on a parameter should not be the main goal of a signal search; logically, the prior and more interesting question is whether the data support the inference of a non-zero signal above background. Therefore, in addition to setting exclusion limits, we search for naïve signal candidates. We perform 18 SuperMAG reports the magnetic-field measurements from each station in locally defined coordinates, oriented along Local Magnetic North and Local Magnetic East. Using the IGRF model, these measurements can be rotated on a station-by-station basis to globally defined coordinates, oriented along True Geographic North and True Geographic East. In what follows, we will work solely with the measurements in these geographic coordinates. In particular, we denote the components of B i (t j ) by B θ i (t j ) (oriented towards geographic South) and B φ i (t j ) (oriented towards geographic East). The vertical component of B i (t j ) will not be relevant for our analysis. See Sec. III B of Ref. [62] for more details on the SuperMAG coordinate systems and the rotations required to achieve field measurements aligned to geographic coordinates.
this part of the analysis in a frequentist fashion, searching for isolated frequencies at which the data variables X k are inconsistent with the absence of a signal as determined using the likelihood function under the assumption of g aγ = 0. We identify 27 such candidate signals in our analysis that were inconsistent with the null hypothesis at 95% global significance. We re-evaluate each such identified naïve signal candidate to test it for robustness: i.e., temporal consistency and spatial uniformity, which are required properties of a physical axion dark-matter signal. These tests consist of splitting the full SuperMAG dataset into a number of smaller data subsets, either by restricting the temporal duration of the data to create temporally disjoint subsets, or by restricting the stations whose data are utilized to create sets of data recorded by disjoint sets of stations. We check if the X k computed from the data subsets are consistent with the posterior on g aγ derived from the Bayesian analysis of the full dataset.
The subsections of this appendix will follow the above structure. Namely, Appendix C 1 will discuss the construction of the X θ and X φ time series and the X k variables; Appendix C 2 will compute the expectation and variances of the X k variables under the appropriate hypotheses; Appendix C 3 will derive the likelihood function for the X k variables and use it to set exclusion limits on g aγ ; and finally, Appendix C 4 will identify the naïve signal candidates in our analysis and test them for robustness.

Time series construction
In Ref. [62], five time series X (n) were constructed based on the five distinctθ-andφ-components of the Φ 1m modes. Here we construct our time series in a similar manner, but instead only require two time series X θ and X φ based on the components of the signal Eq. (16). In particular, we define 19 where Φ θ m (Ω i ) and Φ φ m (Ω i ) are theθ-andφcomponents of Φ m as evaluated at the location Ω i of 19 We note that since C m has units of nT, then X θ and X φ in this work have units of (nT) 2 . This is in contrast to Ref. [62], where the X (n) time series have units of nT. station i. The notation '{i|t j ∈ T i }' indicates that the outer sum is over all stations i which recorded a valid measurement at time t j . The inner sum is taken over 20 ≤ 4 and − ≤ m ≤ . The coefficients C m are computed from the Gauss coefficients g m and h m of the IGRF model (see Sec. II B), per Eq. (12). The IGRF model provides the values for g m and h m at five-year intervals from 1900 to 2020, and the values at all intermediate times are computed via linear interpolation of these coefficients. Thus C m (t j ) exhibits a gradual time dependence in Eqs. (C1) and (C2).
Motivated by the stationarity of the noise in our time series over any given calendar year (see Appendix E 1 in Ref. [62]), we take the weights w θ i (t j ) and w φ i (t j ) in Eqs. (C1) and (C2) to be constant over each calendar year. In principle, these weights could be arbitrary; however, as in Ref. [62] we set the weights within a calendar year in a data-driven fashion using the variances of the measured magnetic fields within each given year. In particular, for α = θ, φ, where T a i is the subset of T i contained entirely within year a, and N a i is the corresponding number of data points in T a i . Note that we make use of the baseline-subtracted (zero-mean) SuperMAG field data, so that Var [B α i ] = [B α i ] 2 ; see the detailed discussion of the appropriateness of the latter choice in Sec. III C of Ref. [62]. The normalizing total weights are then defined by The rest of our analysis works solely with the time series X θ (t j ) and X φ (t j ), rather than the station-bystation data. As mentioned at the end of Sec. II C, the axion dark matter has a finite coherence time T coh (m a ), which may be shorter than the 50-year duration of the SuperMAG dataset. It is therefore convenient to partition our full time series X θ and X φ into shorter segments X θ k and X φ k , each roughly the length of the coherence time T coh ∼ 2π/(m a v 2 dm ) ∼ 10 6 f −1 a . Proceeding in this fashion allows us to analyze each individual segment k coherently [i.e., Eq. (16) can be assumed to be accurate over the whole duration of the segment] and then combine 20 Because the Gauss coefficients g m and h m are largest for low , higher terms in the sums in Eqs. (C1) and (C2) become increasingly negligible. We choose to truncate the sums at = 4 and have explicitly verified that this choice has negligible effect on our analysis. In particular, the results we get by truncating at = 3 or = 5 differ negligibly from those we present with = 4. the individual segments' results incoherently. As we will be interested in setting a bound at a particular axion frequency f a , we take the Fourier transform of each X θ k and X φ k at f a and combine them into several two-dimensional 'analysis vectors' 21 We note that while in Ref. [62] it was necessary to include the Fourier transforms at f a ± f d [where f d = (sidereal day) −1 ] in the analysis vector, this is not necessary in the present analysis, as the axion signal has no f d -dependence (see the discussion in Sec. II D). Therefore, the dimensionality of the analysis vector X k in this work is significantly reduced from its 15 dimensions in the dark-photon dark-matter case considered in Ref. [62] to only two dimensions in the axion dark-matter case considered here. Finally, we note that in order to efficiently make use of the Fast Fourier Transform (FFT), we must approximate T coh in such a way that many frequencies f a at which we set bounds have the same approximate coherence time.
The framework by which we choose our frequencies f a and approximate T coh in a computationally efficient manner is identical to the framework used in Ref. [62]. We refer the interested reader to Sec. V E of Ref. [62] for more details on our frequency choice and coherence-time approximation.

Statistics of X k
Now that we have constructed the primary variables X k for our analysis, we compute their statistics in this subsection. 22 Let us begin with the expectation X k under the hypothesis of a signal with axion-photon coupling g aγ . We parametrize the axion amplitude as so that |c| 2 = 1. In the case of a true axion dark-matter signal with amplitude a 0 and coupling g aγ , the physical field which would be measured by the SuperMAG magnetometers would be the real part of Eq. (16). Because the VSH sum in Eq. (16) is manifestly real, the measured 21 Here and throughout we use x to denote a vector x with 2 components. 22 As in Ref. [62], we assume the variables X k to be Gaussian, so that their statistics are entirely described by their expectation value and variance. See Appendix E 3 of Ref. [62] for validation of this assumption.
field at r = R can simply be written as 23 Let us also define the time series for α = θ, φ, and let H α k represent the subseries of these with the same sampling times as X α k , and letH α k be their Fourier transforms. Substituting Eq. (C7) into Eq. (C5), we find that the expectation value of our analysis vector, under the signal hypothesis, is where c k denotes the normalized axion amplitude [Eq. (C6)] in the k-th coherence time, and where at the ≈ sign we have assumed thatH θ k andH φ k decay rapidly with frequency, so that we may ignore higherfrequency contributions (see discussion after Eq. (36) in Ref. [62]). 24 Now we turn to the variance of the X k variables under the zero-signal hypothesis. As in Ref. [62], we estimate the noise in X k under the assumption that it is stationary within any given calendar year. (This assumption was validated in Appendix E 1 of Ref. [62] by showing that 23 As noted in Sec. II C, this expression for the magnetic field only accounts for the Φ m contributions to the axion dark-matter signal, but Y m and Ψ m may in principle also be present. As in Ref. [62], we ignore such additional contributions in our analysis. If the station locations Ω i were uniformly distributed over the Earth's surface, and the weights were taken to be w θ i (t j ) = w φ i (t j ) = 1 for all stations i and times t j , then Eqs. (C1) and (C2) would approximate a uniform integral over the sphere and so project out any Y m and Ψ m contributions due to the VSH orthogonality relations Eqs. (A6)-(A8). Since this is not exactly the case, it is possible for other VSH modes to 'leak' into our time series X θ and X φ . However, this leakage would at worst affect our search at the level of O(1) factors, and so we neglect such contributions. See Sec. V B 1 of Ref. [62] for further discussion on this point. 24 Note that the coefficients C m that appear in Eq. (C8) exhibit a time dependence. These coefficients thus introduce an additional time dependence which was not present in the dark-photon darkmatter case considered in Ref. [62]. As these coefficients drift by O(10%) on the century timescale, their time dependence affects significantly lower frequencies than are relevant for our analysis.
the noise estimates from different quarters of a calendar year agreed with the full-year estimate.) In particular, let x α (t j ) denote a hypothetical realization of the noise in X α (i.e., under the assumption of no signal g aγ = 0) over a duration τ entirely contained within calendar year a. Then we define the two-sided cross-power spectral density S a αβ for the year a by wherex α (f p ) is the Discrete Fourier Transform (DFT) of x α (t j ) evaluated at one of the DFT frequencies f p,q , and δ pq is the Kronecker delta. We compute the power spectral density S a αβ in a data-driven manner identical to that used in Sec. V C of Ref. [62]. Namely, we partition the calendar year into several chunks, each of which we treat as an independent noise realization, in order to evaluate Eq. (C11). Once the power spectral density has been calculated for each calendar year, the covariance matrix for X k can be computed by combining the noise spectra from different years falling within the same coherence time: where T a k is the duration of the subseries X (m) k contained within calendar year a, so that a T a k is the full duration of the subseries X (m) k , approximately given by the coherence time T coh . The quantities µ k and Σ k are then sufficient to describe the statistics of X k .

Bayesian statistical analysis
Having computed the statistics of the X k variables, we can now analyze them to set bounds on g aγ . We will do so in a Bayesian framework: we will write down a likelihood function for the X k ; then beginning with a prior on g aγ , we will utilize this likelihood function to derive a posterior on g aγ . Apart from a few O(1) numbers, the formulae in this subsection will be almost identical to those derived in Sec. V D of Ref. [62], but we reproduce them here for completeness. Let us begin by writing down the likelihood function for an axionphoton coupling g aγ and normalized axion amplitudes c k , in terms of the observed analysis vectors X k Here we treat each coherence time as independent, so that we may sum over the individual log-likelihoods from each coherence time. Since Σ k is an Hermitian, positivedefinite matrix, we may write Σ k = A k A † k for some invertible A k , and define In terms of these new variables, the likelihood Eq. (C13) can then be rewritten as We can further simplify Eq. (C16) by projecting 25 each term onto the direction of ν k . Namely if we make a fur- 25 It can be shown that the component orthogonal to this projection does not depend on gaγ , c k , or z k . It can therefore be neglected when restricting our attention to the likelihood in terms of z k . See the detailed discussion of this point in the same context in Sec. V D 1 of Ref. [62].
ther change of variables we can write the likelihood as We are interested in setting constraints on g aγ , given data in the form of z k . We are therefore not concerned with the c k amplitudes and so marginalize over them in the likelihood Eq. (C18). The real and imaginary parts of c k are each Gaussian with variance 1/2 [see discussion after Eq. (C6)]. Therefore, their likelihoods are given by Using these likelihoods, we can marginalize over c k in Eq. (C18) to find the marginalized likelihood Now that we have the likelihood for g aγ in terms of z k , we can derive a posterior for g aγ given the observed values of z k . First we must begin with a prior for g aγ . As in Refs. [62] and [79], we take the (reparametrizationinvariant) objective Jeffreys prior [92]. Formally this is defined in terms of the Fisher information matrix [92]. In our context, it takes the form, Then after observing z k , the posterior for g aγ becomes The normalization for Eq. (C23) can be calculated by requiring ∞ 0 dg aγ p(g aγ |{z k }) = 1. Finally with the appropriate normalization computed, we can set a 95%credible upper limit (local significance)ĝ aγ by solving ĝaγ 0 dg aγ p(g aγ |{z k }) = 0.95. (C24) We apply a 25% degradation factor to our upper limit, g aγ →ĝ aγ ≡ 1.25 ·ĝ aγ (C25) to correct for the finite width of a physical axion darkmatter signal. The reason that this degradation factor is necessary is that, thus far, we have assumed that our signal would be exactly monochromatic within a coherence time. However, even within a coherence time, a physical signal would have finite width, comparable to the DFT frequency resolution. This would lead to a suppression of power at the central frequency. In other words, the limitĝ aγ is too strong because it assumes a signal size slightly larger than what a physical signal of finite width would exhibit at its central frequency. In Ref. [62], we estimated that a 25% degradation factor is appropriate to correct for this assumption. Because that estimate was based on a single component of the vectorial dark-photon dark-matter signal, the same numerical estimate applies in the axion case. Additionally, following a prescription similar to that described in Sec. VI C of Ref. [62] (with the appropriate simplifications made for the axion case), we have confirmed explicitly by injecting a physical axion dark-matter signal that a 25% degradation factor is sufficient for our analysis to return an upper limit on g aγ that is consistent with the injected signal parameters. Our results in Fig. 1 show the corrected limitĝ aγ , as defined in Eq. (C25).

Naïve signal candidate re-evaluation
In this subsection, we identify any naïve axion darkmatter signal candidates in our analysis and test them for robustness. Our candidate identification and reevaluation procedure is essentially the same as in Ref. [62]. Here we briefly review the details of the procedure, as well as discuss the results of the re-evaluation.
We define a naïve signal candidate in a frequentist fashion as a frequency for which the observed z k are inconsistent with the absence of a signal. According to the likelihood Eq. (C20) with g aγ = 0, each z k should have Gaussian real and imaginary parts, each with variance 1/2. Therefore, the statistic, should follow a χ 2 -distribution. We can compute an associated p-value where F χ 2 (ν) is the cumulative distribution function (CDF) of the χ 2 -distribution with ν degrees of freedom, and K 0 is the number of distinct segments into which we break our time series (or in other words, the number of values through which the index k ranges). We declare a frequency to be a 'naïve signal candidate' (with 95% global confidence) if p 0 is below the threshold, where N f ≈ 3.3 × 10 6 is the number of distinct frequencies in our range of interest 6 × 10 −4 Hz < f A < (1 min) −1 − 6 × 10 −4 Hz. Based on this criterion, we identify 27 naïve signal candidates; see Table I. The global significance of each candidate can also be characterized by its equivalent one-sided Gaussian-standarddeviation We however do not immediately consider these 27 naïve signal candidates to necessarily be promising axion dark-matter signals: we first re-evaluate each naïve signal candidate by performing several robustness checks to test their temporal consistency and spatial uniformity. Any axion dark-matter signal should be present over the entire 50-year duration of the SuperMAG dataset and should be present in all 508 stations. We therefore reperform our analysis on eight subsets of the full Super-MAG dataset: four temporal subsets and four geographical subsets. The four temporal subsets consist only of the data from the ranges of years 1970-1982, 1983-1994, 1995-2007, and 2008-2019, respectively. Meanwhile, the geographical subsets consist of the data from four ran- σ(p0) p1 p2 p3 p4 ptime p5 p6 p7 p8 pgeo p full 1 2.777777 6.6 × 10 −17 6.2 TABLE I. Naïve signal candidates and associated p-values. p0 indicates the local significance of the candidate in the original analysis, under the hypothesis of no signal. This is also translated into a global one-sided Gaussian-standard-deviation significance σ(p0). The values p1, . . . , p4 indicate how well the signal sizes implied by four temporal subsets agree with the signal size implied by the full dataset. The values p5, . . . , p8 indicate a similar agreement, but for four geographical subsets. Note that pj near either 0 or 1 indicates poor agreement with the full dataset. ptime (respectively, pgeo) indicates the combined significance of all the temporal (geographical) checks. p full represents the combined significance of all eight tests. dom disjoint subsets of stations. 26 For each frequency and each subset, the robustness test consists of checking whether the signal appears in a way that is consistent with the posterior on g aγ derived from the original analysis.
More specifically, let s k,j and z k,j denote the quantities defined in Eq. (C17), but computed using subset j of the full SuperMAG dataset (where j = 1, . . . , 4 refers to the temporal subsets and j = 5, . . . , 8 refers to the geographical subsets). From the likelihood Eq. (C20), we can see that the statistic, should follow a χ 2 -distribution under the assumption of a signal with coupling g aγ . We can again compute an associated p-value where K j is the number of time-series segments in the j th resampling analysis (which may differ from K 0 due to the reduced total duration of the subsets). To test for consistency with the original analysis, we then weight the p-values for each resampling test according to the posterior p(g aγ |{z k }) derived in the original analysis p j = dg aγ p(g aγ |{z k }) · p j (g aγ ).
(C32) Note that p j near either 0 or 1 indicates disagreement with the original analysis, as they imply that the inferred signal present in the subset is much smaller or larger, respectively, than the signal in the full dataset. Finally, we can combine the p-values from all n = 8 tests into a single statistic using a two-tailed version of Fisher's method [93][94][95][96]. Specifically we define the combined χ 2statistic and its associated p-value In addition to p full , we also compute combined p-values p time and p geo for the temporal-only and geographicalonly tests, respectively, by restricting the sum in Eq. (C33) to the appropriate tests and setting n = 4 in Eq. (C34). As noted and discussed at greater length in Sec. VI B of Ref. [62], the temporal and geographical tests here are not entirely independent, as the temporal and geographical subsets generally contain overlapping data. Thus our naïve procedure of combining the results of all the tests as if they were independent in Eqs. (C33) and (C34) is not exactly accurate. Therefore, rather than immediately rejecting (at 95% confidence) all candidates with p full < 0.05, we instead only reject candidates with p full < 0.01, and deem candidates with 0.01 < p full < 0.05 to be in 'strong tension' with our resampling checks.
Tab. I shows the results of our resampling analysis for all 27 naïve signal candidates. Of the 27 candidates, 23 have p full < 0.01, and so we reject them as axion darkmatter signals on the basis of our combined resampling checks. Although none of the other four candidates constitute strong evidence for axion dark matter, they require further discussion. Here we review each one in detail.
Candidate 5 -This candidate has σ(p 0 ) = 6.7, but 0.01 < p full < 0.05. While this signal is thus large, it exhibits strong (but not necessarily definitive) tension with the combined spatio-temporal robustness check. In addition, this candidate is also in strong tension with the geographical tests alone (p geo = 0.033 < 0.05). We also note in passing that this candidate appears at the DFT frequency closest to half of the Nyquist frequency. Based on the tensions with the robustness tests, we do not consider this a robust axion dark-matter candidate.
Candidate 7 -This candidate has σ(p 0 ) = 3.6 and exhibits good agreement with both the geographical tests (p geo = 0.74) and the combined tests (p full = 0.091). It however has 0.01 < p time < 0.05, so that it is in strong tension with the temporal tests. We therefore consider this candidate to be in strong tension with some of our robustness tests, but we cannot definitively rule it out.
Candidates 15 and 17 -These two candidates are reflections of each other across the Nyquist frequency (and thus should be thought of as a single candidate). They exhibit good agreement with all the robustness tests (p time = 0.14, p geo = 0.11, p full = 0.068). However, they exhibit a fairly weak global significance of σ(p 0 ) = 2.3, and so we do not consider these to be strong axion darkmatter candidates. These candidates would be interesting to revisit in an analysis of the one-second Super-MAG dataset to see if they increase in significance. They also motivate the development of further checks in future work to verify the physicality of our naïve signal candidates.
In summary, we find that none of the 27 naïve signal candidates constitute both strong and robust evidence for an axion dark-matter signal. However, four of the candidates warrant further investigation in follow-up work, such as an analysis of the one-second resolution Super-MAG dataset.