Search for dark-photon dark matter in the SuperMAG geomagnetic field dataset

In our recent companion paper [arXiv:2106.00022], we pointed out a novel signature of ultralight kinetically mixed dark-photon dark matter. This signature is a quasi-monochromatic, time-oscillating terrestrial magnetic field that takes a particular pattern over the surface of the Earth. In this work, we present a search for this signal in existing, unshielded magnetometer data recorded by geographically dispersed, geomagnetic stations. The dataset comes from the SuperMAG Collaboration and consists of measurements taken with one-minute cadence since 1970, with $\mathcal{O}(500)$ stations contributing in all. We aggregate the magnetic field measurements from all stations by projecting them onto a small set of global vector spherical harmonics (VSH) that capture the expected vectorial pattern of the signal at each station. Within each dark-photon coherence time, we use a data-driven technique to estimate the broadband background noise in the data, and search for excess narrowband power in this set of VSH components; we stack the searches in distinct coherence times incoherently. Following a Bayesian analysis approach that allows us to account for the stochastic nature of the dark-photon dark-matter field, we set exclusion bounds on the kinetic mixing parameter in the dark-photon dark-matter mass range $2\times10^{-18}\,\text{eV} \lesssim m_{A'} \lesssim 7\times10^{-17}\,\text{eV}$ (corresponding to frequencies $6\times 10^{-4}\,\text{Hz}\lesssim f_{A'} \lesssim 2\times 10^{-2}\,\text{Hz}$). These limits are complementary to various existing astrophysical constraints. Although our main analysis also identifies a number of candidate signals in the SuperMAG dataset, these appear to either fail or be in tension with various additional robustness checks we apply to those candidates. We report no robust and significant evidence for a dark-photon dark-matter signal in the SuperMAG dataset.


I. INTRODUCTION
Over an enormous range of scales from the dwarfgalactic to the cosmological, there is overwhelming evidence for the existence of dark matter (DM) via its gravitational effects. However, despite a broad and decades-long experimental program to detect the effects of any non-gravitational interactions which the dark matter may possess, either in the laboratory or via astrophysical probes, the identity of the dark matter remains elusive. The difficulty of the search for the nature of dark matter stems in part from its extremely broad range of allowed masses, spanning some ∼ 80 orders of magnitude, from ultralight fuzzy dark matter around 10 −21 eV [2][3][4][5][6], up to macroscopic primordial black hole dark matter around ∼ 10 56 eV [7]. Moreover, the various possible DM candidates that populate this allowed mass range give rise to a diverse array of potential phenomenological effects that cannot all be searched for using a single experimental approach. In the past decade or so, there has in particular been a rapid growth of interest in novel experimental techniques aiming to detect bosonic DM candidates that admit a classical wave description, 1 of which * mfedderke@jhu.edu † pwgraham@stanford.edu ‡ derek.jacksonkimball@csueastbay.edu § saarik@stanford.edu 1 To be precise, we mean here that given the local galactic abundance of the DM, ρdm ∼ 0.3 GeV/cm 3 , excitations of the bosonic one well-motivated example is the kinetically mixed [8] dark photon [9], sometimes also referred to as the 'hidden photon'; see, e.g., Refs. [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Such dark-photon dark matter (DPDM) can be produced in the early Universe in a variety of model-dependent and -independent ways; e.g., Refs. [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. In a recent companion paper [1], we pointed out the existence of a new signature of ultralight kinetically mixed dark-photon dark matter: a spatially and temporally coherent, oscillating, terrestrial magnetic field signal that is narrowband in frequency, and that takes a particular vectorial field pattern over the whole surface of the Earth. This signal arises because of the same photon-darkphoton mixing effects responsible for the generation of the signal in, e.g., DM Radio [14]. In Ref. [1], we provided a high-level summary and the results of an experimental search for this novel signal that we undertook using a publicly available geomagnetic field dataset maintained by the SuperMAG Collaboration [43][44][45]. This dataset, which exists primarily for geophysical metrology and solar activity research purposes, consists of time-series magnetic field measurements obtained with unshielded three-axis magnetometers located at O(500) ground stations that are widely dispersed over the surface of the Earth and that, collectively, have been recording data dark-matter quantum field have expected local occupancy numbers (i.e., number of particles per cubic de Broglie wavelength) in the vicinity of the Earth that are greater than 1. This occurs for DM masses lighter than ∼ 10 eV assuming vdm ∼ 10 −3 . arXiv:2108.08852v2 [hep-ph] 30 Nov 2021 continuously since the early 1970s with a sampling rate of (at least) once per minute [43][44][45]. We reported no significant evidence for the existence of a robust dark-photon dark-matter signal in these data in the dark-photon mass range 2 × 10 −18 eV m A 7 × 10 −17 eV corresponding to frequencies 6 × 10 −4 Hz f A 2 × 10 −2 Hz. We thus placed direct observational bounds on the kinetic mixing parameter ε that are complementary to various existing astrophysical constraints [46][47][48][49][50]. In this paper, we supplement Ref. [1] by providing a detailed technical description of this experimental search.
The remainder of this paper is structured as follows: in Sec. II, we briefly summarize the main features of the signal we described in detail in Ref. [1]. We then give a description of the SuperMAG dataset in Sec. III, before giving a high-level description of our analysis strategy for this dataset in Sec. IV. With this high-level overview as a guidepost, we give a detailed technical description of the analysis in Sec. V. The results of this analysis in the form of exclusion bounds on dark-photon dark-matter parameter space are shown at Fig. 4 and discussed in Sec. V G. Our analysis in Sec. V also identifies a number of naïve signal candidates in the SuperMAG data, in addition to placing bounds on parameter space; we test these candidates for robustness in Sec. VI. On the basis of those tests and other indicia, we find no naïve signal candidate for which there is robust evidence of a real signal, although a handful of these candidates would be of potential interest to examine in follow-up work. We discuss our results and conclude in Sec. VII. There are a number of appendices that provide additional information, conventions, or details. Appendices A and B give our conventions for the Fourier transform and vector spherical harmonics, respectively. Appendix C gives some additional technical details of the signal as it appears in the Super-MAG dataset in our analysis construction. Appendix D contains some derivations of important statistical results used in our analysis construction in Sec. V. Finally, Appendix E contains a series of detailed validation checks on the data-driven noise estimation procedures applied in our analysis.

II. SIGNAL
In our recent companion paper [1], we showed that kinetically mixed dark-photon dark matter generates a coherent magnetic field signal across the surface of the Earth of the form as measured in the rotating Earth-fixed frame, where Ω = (θ, φ) is the location on the surface of the Earth (in the geographic co-ordinate system referenced to True Geographic North), ε is the kinetic mixing parameter (as defined in Ref. [1]); m A ≡ 2πf A is the DPDM mass (with f A being the corresponding cycles-per-second frequency); 2 R is the radius of the Earth; A m are the (complex) amplitudes describing the (amplitude and phase of the) three different polarization modes of the dark photon in the vicinity of the Earth (as measured in a nonrotating fixed inertial frame); 3 Φ m are vector spherical harmonics (see Appendix B for conventions); and the additional frequency f d = (sidereal day) −1 appears in the m = ±1 modes owing to the rotation of the Earth [1]. As discussed in detail in Ref. [1], Eq. (1) is a good description of the signal within a single DM coherence time T coh ∼ 2π/(m A v 2 dm ) ∼ 10 6 f −1 A , where we have taken v dm ∼ 10 −3 to be a representative value for the galactic DM velocity dispersion. Within that coherence time, the complex amplitudes A 0,± characterizing the local DM field remain approximately constant, but in general they evolve significantly from one coherence time to the next. 4 Indeed, since the local DPDM field can be thought of as being comprised of the sum of a large number of independent plane waves with frequencies f ∼ f A 1 + O(v 2 dm ) , each with its own random phase, within a single coherence time each of the real and imaginary parts of the A m can by virtue of the central limit theorem be described as a random draw from a zero-mean normal distribution with standard deviation √ ρ dm /( √ 3m A ), such that together they satisfy where |A | 2 = (3) 2 We work in natural units where = c = 1. The mass-frequency conversion is thus f A ≈ 24 mHz × (m A /10 −16 eV). 3 As discussed in Ref. [1], the A m technically describe the amplitudes of the polarization modes of the sterile component (in the interaction basis) of the DPDM, as measured in the vicinity of the Earth but well outside the atmosphere, and in the inertial frame. Our convention for the A m is such that A ± = ∓ 1 √ 2 A x ∓ iA y so that the Cartesian components (in the inertial frame) of the dark vector potential are given by We employ the shorthand A ± ≡ A ±1 . 4 Generically, this implies that the 'dark electric field' E ∼ m A A of the dark-photon dark matter is not simply a vector of fixed direction in 3D space with an amplitude oscillating at frequency f A . Instead, each of the components of E executes oscillations at frequency f A , which implies that E has a periodic (with period T A = f −1 A ) variation of both its instantaneous amplitude and its direction in 3D space. In the generic case, |E | does not vanish instantaneously at any moment in time, and the tip of the unit vectorÊ traces out a closed periodic curve with period f −1 A ; that curve evolves secularly on characteristic timescales of order the coherence time.
the angle-brackets · · · describe an average over times τ much longer than the coherence time, τ T coh ; and ρ dm is the average local (to the Earth) dark-matter mass-density, which we will take to be fixed at ρ dm = 0.3 GeV/cm 3 throughout this paper. We pause to note that there is some discussion in the literature regarding the appropriate treatment of the DPDM polarization state; see, e.g., discussion in Refs. [29,51]. We have assumed above, and will continue to do so throughout this paper, that the DPDM field is a sum of plane waves with random individual phases and randomly oriented individual polarization states; this guarantees that the overall polarization state will necessarily randomize over a coherence time. However, certain production mechanisms (e.g., misalignment) may possibly give rise to a polarization state that does not evolve (significantly) in time today, because every individual mode is produced in the early Universe with the same (or similar) polarization; provided that structure formation and interactions with matter do not spoil this, the DPDM field today would then consist of a sum of plane waves with individual random phases but all the same (or similar) polarization states. Such a DPDM field would still exhibit phase-decoherence over a coherence time and thus amplitude fluctuations from one such time to the next, but the polarization state would of course not randomize significantly from one coherence time to the next.
A crucial feature of the signal is the factor of (m A R) in Eq. (1), which encodes that the signal suffers a suppression in the ratio of the radius of the Earth to the (Compton) wavelength of the DM. Naïvely, however, one might be tempted to think that the depth of the atmosphere L atmos R would instead be the relevant length scale governing this suppression (see, e.g., brief comments in Ref. [46]), which would have implied that this factor would instead be replaced by a factor of m A L atmos ∼ 10 −2 m A R, dramatically weakening the signal prediction. We discuss at length in Ref. [1] why this is not in fact the case.
Numerically, the signal Eq. (1) is expected to have an amplitude on the order of where we assumed for the purposes of this rough estimate that |A | ∼ √ 2ρ dm /m A with ρ dm = 0.3 GeV/cm 3 , and evaluated the maximum value of the field on the Earth's surface. While this signal is very small in amplitude (many orders of magnitude smaller than the static geomagnetic field, which is of order B ⊕ ∼ 0.5 G [52]), it is at nonzero frequency, effectively monochromatic with a long coherence time and has a very particular global field pattern over the surface of the Earth; it can thus be meaningfully distinguished from many noise sources via techniques that are tailored to search for the specific spatial and frequency structure of the signal.
The magnetic field amplitude of the signal Eq. (4) is also potentially much smaller than the individual point-in-time, single-station, single-field-component digital measurement resolution of the magnetic field stations that contribute to SuperMAG, which are typically in the 10-100 pT ≈ 100-1000 nG range [53][54][55][56][57][58][59][60][61]. However, given that the instantaneous fluctuating random noise in the detectors greatly exceeds 5 this measurement resolution [43][44][45] (see generally Sec. V C and Appendix E), the fact that the signal amplitude is sub-readout-resolution does not degrade the sensitivity for our signal search. This can be understood from the following intuitive argument: suppose the station digital measurement readout resolution is ρ, and the signal has amplitude A in some given field component. If A < ρ, then on average only a fraction ξ ∼ A/ρ of single-station, single-fieldcomponent measurements will be impacted by the signal being present. This is because only that fraction of pointin-time noise realizations lie close enough to the breakpoint in the digital readout rounding for the presence of the signal to alter the reading of the magnetometer. However, for those fraction ξ of measurements, the readout is changed by a full resolution unit ρ > A, which is ρ/A ∼ 1/ξ larger than the signal amplitude A. These two effects thus effectively cancel out. And indeed, we have verified numerically as well in simple cognate examples that, in the presence of the readout digitization noise, a narrowband signal of sub-resolution amplitude added to super-resolution instantaneous random noise remains narrowband and visible in the finite-resolution data, provided of course that the signal amplitude is larger than the averaged-down noise level.

III. SUPERMAG DATA
We now turn to a general description [Sec. III A] of the magnetic field dataset we have analyzed in this paper to search for the signal shown at Eq. (1), before turning to a more extensive discussion of various salient details in Secs. III B-III D.

A. Overview
The SuperMAG Collaboration [44,45] maintains and makes available for research purposes a large archival dataset of three-axis geomagnetic field measurements sourced from around O(500) individual measurement stations 6 which are geographically dispersed across the surface of the Earth; see in a common format, in a well-defined co-ordinate system [45], with common temporal measurement resolution and synchronization, and (where relevant) are preprocessed in a common manner.
The data product with which we will be primarily interested in this paper is their 'low fidelity' dataset, which encompasses measurements made with one-minute resolution, beginning in 1970 [43,45]. The SuperMAG Collaboration has also recently released a 'high fidelity' dataset of measurements taken by O(100) stations with one-second temporal resolution in the time-frame 2012-2020 [43]. We defer analysis of the one-second resolution data to future work.
While various individual stations contributing to this dataset have come online and/or gone offline since 1970, and even otherwise operational stations do not have 100% uptime (so that data from individual stations are unavailable during certain periods of time), the cumulative number of stations in this dataset for which at least some amount of data are available and analyzed in this work is 494, and the number of stations operational in recent years has fluctuated between around 150 and 250 at any given time [43,45]; see Fig. 2. Note however that for technical reasons, 7 we restrict our attention to the 48 years of data taken starting at the beginning of 1972 and concluding at the end of 2019.

B. SuperMAG co-ordinate system
The three-axis magnetic field measurements supplied by SuperMAG are reported for every station in locally well-defined co-ordinate systems whose definitions vary from station to station [45]. As described in this subsection, these local co-ordinate systems must be rotated to 7 Insufficiently many stations are operative in 1970 and 1971 for our analysis to be applied; see also footnote 39 for a similar point. reporting as a function of the date. We discuss the manifest annual trends in Sec. III D.
obtain the magnetic field components in the global geographic co-ordinate system that we will require for our analysis.
As discussed in detail in Ref. [45], SuperMAG assumes that each station has correctly reported the orientation of the vertical component (i.e., the component radially directed at the center of the Earth) B z of their field, into the ground being positive. However, owing to a proliferation of possible conventions in use by individual station operators to report the two orthogonal components of the field in the plane perpendicular to the vertical ('the horizontal plane'), SuperMAG undertakes a data-driven procedure to rotate the horizontal field components reported by every station into a co-ordinate system whose orthogonal basis vectors are oriented along Local Magnetic North (LMN) and Local Magnetic East (LME); the corresponding magnetic field components along these directions are B lmn and B lme , respectively.
The instantaneous rotation angle about the vertical axis that would be required to perform this transformation is obtained in an unambiguous way from the field data themselves by demanding that the 'typical value' (as defined in Ref. [45]) of B lme in a sliding 17-day window period centered on the observation time is zero by definition [45]; the rotation angles that are actually applied to the data are smoothed versions of these instantaneous rotation angles, where the smoothing is performed over the same 17-day period [45,63].
For the purposes of the discussion here, it is only relevant to note that the 17-day (∼ 1.5 × 10 6 s) time windows used in these procedures are much longer than the intrinsic timescales associated with the oscillation of the DPDM in our mass range of interest: 6×10 −4 Hz f A 2 × 10 −2 Hz, corresponding roughly to 2 × 10 −18 eV m A 7 × 10 −17 eV. Since our analysis procedure will construct an observable linear in the magnetic field, the principle of superposition implies that the rotation pro-cedure cannot induce or mask signals that are in-band. 8 While the local co-ordinate systems described above have the advantage of being unambiguous, their definitions are inherently local and time-dependent. In order to obtain the field components in the rigid, timeindependent, global geographic co-ordinate system required for our analysis, the field components in these local co-ordinate systems must be rotated to the global frame using the per-station time-dependent true magnetic declination angles 9 δ(Ω, t). The SuperMAG data products include the time-dependent declination angles δ(Ω, t) for each station, assuming the International Geomagnetic Reference Field model of the appropriate epoch [45,63]; these declination angles again vary only on timescales that are out-of-band for our DPDM mass range of interest. A simple 2D rotation about the vertical axis can thus be applied to yield the magnetic field components B tgn and B tge reported along the directions of True Geographic North (TGN) and True Geographic East (TGE), respectively: Finally, note that the field components in the global geographic (i.e., Earth-fixed) co-ordinate system (r,θ,φ) are related to those discussed above by B θ = −B tgn (recall,θ points to the geographic South Pole), B φ = B tge , and B r = −B z (r points locally out of the ground, 10 whereas the SuperMAG convention is to measure B z positive when pointing down).

C. Post-processing by SuperMAG
In addition to the 17-day windowing procedure outlined in Sec. III B that is used to obtain field compo- 8 One might however also naïvely be concerned that, even if this procedure might not induce or mask in-band signals, it might somehow impact the coherence of the DPDM signal over timescales longer that 17 days (provided that T coh > 17 days). This is however not the case: we remind the reader the statement that the coherence time T coh ∼ 10 6 T A [where T A = 1/f A ] is simply the statement that the width of the DPDM signal peak in Fourier space is ∼ 10 −6 of the carrier frequency f A set by the DPDM mass. For our range of interest, 6 × 10 −4 Hz f A 2 × 10 −2 Hz, and the above considerations imply that all the information regarding the signal coherence is similarly restricted to (approximately) the same frequency range. On the other hand, the 17-day smoothing and rotation effects discussed here will only impact frequencies at or below (17 days) −1 ∼ 7 × 10 −7 Hz. As such, these modifications do not impact the coherence of the signal in the data. 9 Conventionally δ is defined as the angle between True Geographic North (TGN) and LMN, with the sign convention chosen such that δ is positive when the direction of LMN lies eastward (i.e., clockwise on a compass vane) of TGN, and is negative when LMN lies westward (i.e., counterclockwise on a compass vane) of TGN. 10 We assume an exactly spherical surface for the Earth through this paper; this approximation is accurate at the 0.3% level [64].
nents in the SuperMAG co-ordinate system, the default SuperMAG data product magnetic field measurements have also been post-processed to remove a 'baseline' field component consisting of a sum of time-varying diurnal and slower annual components (as well as a constant offset irrelevant for the purposes of searching for a timedependent signal, as here).
The procedure used to perform this subtraction is detailed in Ref. [45]; for the present purposes we simply note that the procedure utilized to remove the diurnal component involves examining magnetic field data in discrete coarse-grained intervals of 30 minutes in length. While this diurnal baseline subtraction can result in the removal of even quite monochromatic signals with periods longer than 30 minutes (see, e.g., Fig. 6 of Ref. [45], where a strong six-hour signal is removed from the B lmn data from a single station), the coarse-graining of the data into 30 minute intervals implies that this procedure should not significantly impact any frequencies somewhat higher than (30 min) −1 ∼ 5 × 10 −4 Hz, although it could impact frequencies around or below this. Since the lower end of our frequency range of interest is f A ∼ 6 × 10 −4 Hz, the effects of this diurnal baseline subtraction are mostly out-of-band for the DPDM mass range of greatest interest to us. We have explicitly re-run our analysis pipeline on the non-baseline-subtracted dataset that is also available from SuperMAG [43] and verified that the diurnal subtraction is not observed to remove any DM-likeline features in the results near our frequency range of interest.
The annual baseline subtraction procedure makes use of data which are aggregated using the same 17-day sliding window that was employed to determine the coordinate system rotations, and is thus also well out-ofband. Therefore, it would be consistent to use the data either with or without these two time-dependent baseline subtractions for the purposes of our analysis. However, in order to determine the appropriate data-driven weighting to give the measurements from each station in our analysis [see Eqs. (13) and (14) below], it is more appropriate to utilise the data with the time-dependent baseline subtracted, as this disregards some noise below our frequency range of interest.
Later in our analysis treatment we also assume that a station that is not reporting data reports exactly zero field (instead of the mean field). For consistency with this treatment, we must work with data whose mean DC value is zero in order to avoid introducing artificial discontinuities of order the size of the mean field. Therefore, we will work with the fully baseline subtracted data throughout.

D. Temporal features in SuperMAG data
It is typical for the number of stations that are active to change significantly at the beginning of a calendar year; see Fig. 2. This is likely due to a tendency for stations to report data associated with specific calendar years. Importantly for us, this means that the amount of available data fluctuates relatively little within a calendar year, but may fluctuate significantly between calendar years. In our noise analysis, we will therefore compute separate noise spectra for each calendar year. This, of course, assumes that the noise level remains relatively constant within a calendar year. We assess the validity of this assumption in Appendix E 1.

IV. ANALYSIS STRATEGY
We now turn to a high-level description of the analysis procedures we have utilized in order to search for the signal (described in Sec. II) in the SuperMAG data (described in Sec. III). Details of the implementation of this analysis follow in Sec. V.
Station i, located at geographic co-ordinates Ω i = (θ i , φ i ), reports a time series of three-axis magnetic field measurements; we denote the field measured at time t j as B i (t j ). We denote the set of sampling times at which station i reports valid measurements as T i . As already described, the T i vary among stations, making a straightforward analysis of the individual stations fairly complicated.
Two possible analysis approaches suggest themselves: (A) from the predicted signal Eq. (1), one could construct the expected per-station signals B signal i (t j ) at every time t j ∈ T i . One could then perform a simultaneous joint search in the observed B observed i (t j ) over all stations for the correlated signal predictions in all the stations, and thereby extract a signal or limits on the value of ε as a function of m A ; or (B) one could exploit the observation that our signal Eq. (1) is predicted to be in only a small number of global vector spherical harmonics (VSH). Therefore, one could instead first extract from all the individual station measurements B observed i (t j ) a small number of time series that give the projections of the entire collection of station measurements onto the field components of the independent global vector spherical harmonic modes of interest. One could then perform a search on these distilled VSH component time series for the signal Eq. (1) and thereby extract a signal or limits on the value of ε as a function of m A . For technical reasons, we find it simpler to utilize approach (B).
Our analysis strategy will thus be to first identify the appropriate global VSH components of interest to find a signal of the form Eq. (1); we find that there are five such components of interest, which we denote X (1) , . . . , X (5) . We will then combine all station measurements B i (t j ) available at a given time t j to extract the values of X (1) (t j ), . . . , X (5) (t j ), which are the distilled time series of the VSH components previously mentioned. As noted, the signal Eq. (1) is very narrow in frequency, so it is most appropriate to search for the signal in the frequency domain; however, since the total duration of the data-taking for the available Super-MAG data is in many cases significantly longer than the signal coherence time, a straightforward Fourier transform of the full time series would in general result in a non-monochromatic peak in frequency space if a DM signal were present; extracting a rigorous limit or signal amplitude estimate would then require knowledge of the exact shape of this peak, which is a fairly nontrivial problem that relies on detailed knowledge of the velocity dispersion of the DM (see, e.g., Refs. [65][66][67][68][69][70][71] on this and related points applicable to dark-photon and axion-like dark matter). In order to avoid this issue, whenever the coherence time of the signal is shorter than the available data-taking duration, we instead break up our full time series X (1) (t j ), . . . , X (5) (t j ) into a number of shorter subseries, each of which has a duration of (approximately; see next section) a single coherence time for the frequency of interest. We then Fourier transform each of these subseries to the frequency domain, and because the signal is then coherent in each of the subseries by construction, we can conduct independent searches for a monochromatic (i.e., single frequency bin) signal in the frequency domain subseries. As a second step, we then incoherently stack the results obtained from the searches on subseries into a single result for the frequency of interest, which we do taking into account that the signal phase and DM amplitude vary stochastically from one coherence time to the next. The requisite estimates for the noise in the VSH time series are obtained in a data-driven fashion within each coherence time, as detailed in the next section.
We utilize a Bayesian analysis framework: assuming a reparametrization-invariant (Jeffreys) prior on ε at each mass m A , we use the SuperMAG data to extract the fully marginalized Bayesian posterior distribution on ε, from which we extract upper limits on ε at each mass m A ; see Refs. [69,72] for a similar approach.
In the cases where our search indicates the possible existence of a signal at some threshold significance in the full dataset, we first verify that it has the appropriate frequency-domain width for a DM signal. If the candidate peak passes this test, we then perform subsampling checks to test whether or not the full-dataset signal is consistent with being the expected global DM signal, or whether the degree either of geographical variation in the signal between different random selections of stations, or of temporal variation in the signal between different disjoint subsets of the data broken up over time, are too great to be consistent with the DM interpretation, perhaps indicating that the signal is being driven by some large noise fluctuation in a small number of stations, or for some finite duration of time.

V. ANALYSIS DETAILS
In the previous section, we provided a high-level overview of our analysis strategy; in this section we will provide a detailed description of the analysis. We begin in Sec. V A, with a discussion of the selection of the Expected value of X k in the presence of the signal Eq. (1) x (m) (tj) Sec. V C Hypothetical realization of the data X (m) (tj) over a duration τ contained entirely within a calendar year a five time series X (1) (t j ), . . . , X (5) (t j ) and the procedure by which we combine the different station measurements into these five time series. Additionally, we discuss the breaking of these time-series data into single-coherencetime subseries. In Sec. V B, we discuss how the same procedure affects a hypothetical signal Eq. (1), as this will determine the expectation values of the X (n) (t j ) that enter in the likelihood function we utilize to construct the posterior distribution on ε. A noise estimate is also required to construct the likelihood function; in Sec. V C, we outline our data-driven noise estimation procedure (with validation checks discussed in Appendix E). In Sec. V D, we construct the likelihood function, and then derive the posterior distribution on ε given the observed SuperMAG data. Finally, in Sec. V E, we address a technical point related to the choice of frequencies in our analysis and the way in which they relate to an approximation we make for the coherence time of the signal of interest. For the convenience and reference of the reader, we collect in Tab. I a variety of the analysis variables we will define in this section, a cross-reference to where they are defined, and a brief description of each.
A. Time series

Selection of VSH components
The first point to address is the selection of the appropriate time-series VSH components on which to per-form our analysis. We would like to extract a set of variables defined on the available magnetic field measurements which keep only the information in the measured fields which could have a spatial overlap 11 with the signal; that is, the dot-product B(Ω i , t j ) · B i (t j ) of the expected signal B(Ω i , t j ) and the observed fields B i (t j ) should form the basis of the information we wish to extract. The signal Eq. (1) is proportional to the expression for some complex a m (t j ) that are related to A m . Now, noting that the B i (t j ) are real, it is easy to see that using the explicit expressions for the Φ 1m (Ω i ) VSH that are displayed at Eqs. (B9)-(B11), it can be shown that this sum can be written as a linear combination of the variables These X i (t j ) hold all the information about the observed fields at station i that could possibly spatially overlap with the expected signal we wish to constrain or observe, and we can thus structure our analysis around these variables.

Combination of stations
To combine the results from all the stations into a small number of time series, we simply take the weighted averages over all stations of these projections on the VSH components: 11 The temporal overlap is considered by going to the frequency domain later.
where the notation '{i|t j ∈ T i }' indicates that the sum is over the set of stations i such that there is a valid field measurement from station i at time t j . 12 The weights w (n) i (t j ) we choose for the station at location Ω i will be taken to be constant within the time span over which we will assume stationarity of the noise distributions [one calendar year; see Secs. III D and V C], so that the noise distributions we estimate from the data are informed entirely from the magnetic field noise and from the fluctuations in which stations are reporting, and we do not inject additional temporal variation via explicitly timedependent weights.
The choice of weights for each period of assumed noise stationarity (i.e., one calendar year) could in principle be arbitrary; however, a reasonable assumption is to take the weights to be informed by the per-station white noise levels, assuming the noise between stations is uncorrelated (see discussion in Sec. V C below). That is, we will take w (n) i for n = 1, 2 to be the inverse of the station-i white noise level in B φ i over a given year, while w (n) i for n = 3, 4, 5 will be taken to be the inverse of the station-i white noise level in B θ i over a given year. More specifically, for all t within year a, we take where T a i is the subset of T i [see Sec. IV and below Eq. (12)] contained entirely within year a, and N a i is the corresponding number of samples in T a i . The normalizing total weight is then simply defined by note that even though in our analysis all the w i (t j ) are themselves constant within a year, W (n) (t j ) may still change on more rapid timescales because the number of stations reporting generically changes over time.

Coherent-signal data subsets
The X (n) (t j ) time series contain all the relevant information we need to proceed with our data analysis in Sec. V. 12 We could equivalently formulate this as saying that the sum is over all stations i, but that the weights are zero-ed out at all times when a station is not reporting valid data: w (n) However, as we have already explained in Sec. IV, the coherence time of the signal can be shorter than the full duration of available SuperMAG data, and we wish to avoid having to search for signals that have a resolvable width in frequency space, as this complicates the analysis significantly (and depends in part on the exact DPDM lineshape). Instead, we will perform our search over total durations that are longer than the intrinsic signal coherence time by first analyzing the data coherently within each separate, disjoint coherence time, and then incoherently combining the results from these single-coherencetime searches. By way of concrete example, we mean that if we are confronted with searching for a signal with a sixyear coherence time (f A ∼ 5.3 mHz for v dm ∼ 10 −3 ), we would separate the total 48 years of available SuperMAG data into eight separate consecutive data subsets. We then perform eight independent fully coherent searches for the signal, one in each of these subsets; finally, we stack these search results incoherently to obtain a final search result.
To this end, let us denote by X (n) k the subseries of the time series X (n) that contains only the data from the k-th disjoint interval [k = 1, . . . , K] of duration T within the full SuperMAG data taking duration. As we discuss below in some detail in Sec. V E, we will take T to be approximately equal to the coherence time for the signal: 13 We will analyze each of the subseries X (n) k independently in the frequency domain; we denote the Fourier transform (FT) of the subseries X  k (f ). We can see from Eq. (1) that a signal in the data would contribute power not only at the cycles-per-second frequency corresponding to the DPDM mass, f A = m A /(2π); to the extent that the data contain a signal spatially oriented such that there is some m = ±1 contribution, there will also be power at f = f A ± f d where, as before, f d = (sidereal day) −1 . Therefore it will be relevant to consider the Fourier transformsX 14 Actually, since we obtain the FT of a time-domain signal of total duration T and discrete sampling cadence ∆t (for a total of N = T /∆t samples) via the Discrete Fourier Transform (DFT) [or, more precisely, by the Fast Fourier Transform (FFT) implementation of the DFT], we only obtain independent frequency information at a discrete set of predetermined frequencies f k = k∆f 13 Unless more precision is required to avoid confusion, in order to avoid the repetitive incantation of 'approximate coherence time' and other such caveats, we will hereinafter simply refer to the duration of time T as 'the coherence time', and to any such time period of duration T as a 'coherence time', with this approximation implicitly understood. 14 Note however that it is possible that some of theX where ∆f = 1/T and k = 0, . . . , N − 1. We are thus not generically able to obtain (at least not within the context of the FFT) the FT valueX (n) k (f ) at exactly all the frequencies f A and f A ± f d . Instead, we will consider the FT at f A , which we will always by construction take to be an exact DFT frequency (f A = m∆f for some integer m; see Sec. V E); and at f A ±f d , wheref d is the closest multiple of ∆f to f d (i.e.,f d ≡ n∆f for the integer value of n such that |n∆f − f d | is minimized). We note that a refinement of this approach, in particular one that considers the full lineshape in the Fourier domain, would likely be possible at additional computational expense.
With this in mind, we define a new 15-dimensional 'analysis vector' X k that contains the values of the FT X These analysis vectors X k are the central objects we use in Sec. V D to construct the likelihood function on which our analysis is based; we will consider them to be multivariate Gaussian variables, an assumption we validate in Appendix E 3. As such we will need to know their expectation values given an injected signal [Sec. V B] and their covariance matrices [Sec. V C].

B. Signal
Having described the general procedures we utilize to search for signals appearing in the relevant VSH compo-nents of the SuperMAG data in the previous subsection, in this subsection we will derive the expected values of the X k variables which arise under the dark-photon darkmatter signal model, Eq. (1); we denote the expected X k under the signal hypothesis with parameter ε by X k .
In principle, this derivation amounts to simply substituting Eq. (1) into the definitions of the time series in Eqs. (7)-(11); however, it is useful to examine intermediate results here, so we will develop this section pedagogically.
As a first step, it is useful to more explicitly understand the signal Eq. (1) in the case where the DPDM is polarized along any of the three inertial Cartesian axes [i.e., the set of rigid, mutually orthogonal axes fixed in space with respect to the (average) positions of a field of distant stars, not the body-fixed axes rigidly attached to the rotating Earth]. To this end, we define variables c i which define the orientation of the DPDM field along the inertial i-axis for i = x, y, z: or, in terms of the A ±,0 , we have These are normalized such that for a (hypothetical) linearly polarized DPDM signal that is oriented along the unit vectorn (in the inertial frame) and that satisfies 1 2 m 2 A |A | 2 = ρ dm , we would have c i =n i where i = x, y, z andn i denotes the i-th Cartesian component of n. Expanding Eq. (1) in terms of these variables, we can write where we have defined 15 Here, and throughout, we use x to denote a vector x with 15 components, and y to indicate a vector y with three (usually spatial) components.
Here Ω ≡ (θ, φ) denotes the spherical co-ordinates of the observation point on the surface of the Earth, andθ and φ are the associated unit vectors (recall: θ, φ,θ, andφ are all defined in the body-fixed frame that an observer co-rotating with the surface of the Earth would naturally use). Substituting these expressions for B into Eqs. (7)-(11) and Fourier transforming yields the contribution of each polarization to X k .
We define several auxiliary time-dependent functions H (n) (t j ) [n = 1, . . . , 7] which will allow us to more compactly express our results provided that the per-station weights are taken to be 16 w for n = 3, 4, 5: where W (θ) and W (φ) are defined as in Eq. (15). Armed with these functions, we are in a position to write down the contributions to X k . For instance, in the case that the signal is a dark photon oriented entirely in the z-direction (and has a phase such that c z ∈ R), we have 16 From Eqs. (7)-(11), we see that the X (n) i for n = 1, 2 depend only on B θ i , while the X (n) i for n = 3, 4, 5 depend only on B φ i . Therefore, since we use per-station data-driven estimates for the station weights (see Sec. V A 2), it stands to reason that the weights applied for X (n) i for n = 1, 2 should be common, while those for X (n) i for n = 3, 4, 5 should also be common, but with the latter common value distinct from the former.
Here, H (m) k represents the subseries of H (m) consisting of the same sampling times as the subseries X is its Fourier transform. The series 1 k is a series of 1's at these same sampling times, and1 k its Fourier transform. 17 Note that we have made the approximation in moving from Eq. (35) to Eq. (36) thatH (m) k decays rapidly with increasing frequency, so that we may discard high frequency contributions (i.e., those at f = 2f A , and f = 2f A ±f d ). We verified that, with our choices of weightings, this is a valid approximation: for instance, H is typically of order a few 17 Therefore typically1 k (0) = T . However, we do adjust this to account for the situation in which no stations report valid measurements for some subset of times within the k-th coherence time, or the analysis duration for any one k happens to be shorter than T (e.g., for the last interval).

percent, andH
(m) This approximation has the computational advantage that we need only know the Fourier transforms of the H (m) k (t j ) at three frequencies, and so we can avoid performing an FFT.
Under this same assumption, it is not difficult to see that if the signal is oriented along the z-direction but with ic z ∈ R, we have We may likewise define µ xk and µ yk as the contributions to X k coming from the x-and y-polarizations, and an exact analog of Eq. (37) holds for these too; the full expressions for µ xk and µ yk are shown in Appendix C. It follows immediately that the full expression for the expectation value of the X k under the signal hypothesis Eq. (1) can be written in terms of the µ ik (i = x, y, z) as where the c ik (i = x, y, z) encode the inertial-frame polarization state of DPDM during the k-th coherence time, and * denotes complex conjugation.

Signal in other VSH modes
The signal Eq. (1) was derived in Ref. [1] under the assumption of an exactly spherical geometry; i.e., assuming that the spherical ionosphere acts as the outer boundary for the lower atmospheric cavity in which the dark-photon signal is sourced. In this case, the dark photon sources only a Φ 1m component of the magnetic field. However, as we discussed at length in Sec. II.B of Ref. [1], details of the ionosphere call this assumption into question, and suggest that the aspherical magnetopause may instead act as the outer boundary of the geometry. We showed in Sec. III.C of Ref. [1] that when the spherical ionospheric outer boundary assumption is relaxed, the signal Eq. (1) in general receives additional contributions from other VSH (e.g., Ψ m and Y m contributions; see Appendix B for definitions); however the Φ 1m component shown at Eq. (1) remains correct to leading order in an m A R( 1) expansion.
In principle, these additional field contributions are distinguishable from the Φ 1m component at Eq. (1) due to the global orthogonality of the VSHs; see Eq. (B8). If the station locations Ω i were uniformly distributed over the Earth's surface and the weights were taken to be w (n) i (t j ) = 1 at all stations i and times t j , then the definition of X (n) in Eq. (12) would approximate a uniform integral over the sphere in the limit of many stations. This would project out any Ψ m or Y m contributions to the observed magnetic field B i (t j ), leaving only the contributions from Eq. (1). Following the analysis through, this would imply that Eq. (38) would give the exact signal expectation for the X k .
However, due to the nonuniformity of the station distribution, differing noise levels among stations, and variations in the number of stations reporting at a given time, Eq. (12) for X (n) deviates from approximating a uniform integral over the sphere. This implies that field contributions from other VSH modes arising from the magnetospheric asphericity could give unsuppressed contributions to the time series X (n) ; we estimate that this 'leakage' of other VSH components into X (n) could be at the level of tens of percent. However, while such contributions in principle enter the X (n) in such a way that the expected X k in the presence of the full signal that includes these other VSH modes would deviate from Eq. (38) at the level of an O(1) factor, it would require a highly unlikely environmental fine-tuning for these modifications to completely cancel the signal contribution Eq. (38) that we search for. For instance, the asphericity in the magnetopause is variable with Solar activity as its shape is strongly sculpted by the radial outflow of the variable Solar wind, and other Solar activity (Coronal Mass Ejection events, etc.); the Earth also rotates inside of it. It would be exceedingly surprising for some conspiracy between the stochastically varying DPDM field and the evolving magnetopause shape in which the Earth rotates to somehow engineer cancellation of all three components of the vectorial signal Eq. (1) as it enters the X (n) at Eq. (12), and for that cancellation to be maintained precisely for O(50) years when considered over all O(500) stations that switch on and off over time and have varying noise levels completely uncorrelated with the DPDM signal.
While a more refined future analysis may hope to deal with these signal additional contributions more precisely, we are satisfied that these considerations imply that our search is still accurate at the level of (at worst) O(1) factors even when they are present and not explicitly accounted for.

C. Noise spectra
The statistical analysis of the SuperMAG magnetic field dataset-as expressed in terms of the variables X k [see Sec. V A]-in order to search for a signal of the form X k [see Sec. V B] requires a quantitative estimate of the noise; we utilize a data-driven noise estimation procedure, which we detail in this subsection.
Our analysis is constructed around the assumptions that the noise in the data time series X k is (1) Gaussian, and (2) statistically stationary within each calendar year. We quantify the extent to which (1) and (2) are acceptable assumptions in detail in Appendix E.
Let x (m) (t j ) [m = 1, . . . , 5] represent a single hypothetical realization of the data time series which we have denoted as X (m) (t j ), taken over some time span of duration τ contained entirely within a single calendar year a. Assume the data are taken with a measurement cadence ∆t, such that τ ≡ N ∆t with N ∈ Z and, here, t j = j∆t for j = 0, . . . , N − 1; additionally, we consider x (m) (t j ) to be obtained under the assumption that no DPDM signal is present in the data. 18 Then, x (m) (t j ) is simply a single hypothetical duration-τ realization of the noise in the 18 Note that even if any true dark-photon signal were present in the data, it would have to be very large to invalidate this approach. The DPDM signal line has a width of order σ f ∼ 10 −6 f A . However, the spacing of the DFT frequencies in Eq. (39) is approximately (∆f ) ∼ 1 × 10 −6 Hz if τ = 16834 min. Because our frequency range of interest is 6 × 10 −4 Hz f A 2 × 10 −2 Hz, this means that (∆f ) lies in the range 1700 (∆f ) /σ f 50. A true DPDM signal would thus have to be huge, at least 50 times larger than the noise level in neighboring bins, to make even an O(1) impact on the noise estimate. For signals smaller than this, the estimate we have outlined here is acceptably accurate. For a large signal, the noise estimate outlined here would be formally incorrect; however, we would still see an obvious signal candidate in this case, but further analysis would be required to extract an accurate noise estimate; see, for instance, our signal injection analysis in Sec. VI C and Fig. 6. data time series X (m) (t) in year a. We define the twosided cross-power spectral density of the noise for year a by where · · · ε=0 denotes the expectation taken over all possible noise realizations [i.e., with no signal, ε = 0], evaluated at one of the set of DFT frequencies f p,q (see below for discussion and definition of f p ), and δ pq is the Kronecker delta. Our data-driven noise estimate of year a is constructed from S a mn (f ), which we wish to estimate from our single realization of the actual data time series, X (m) (t). One of our fundamental analysis assumptions is that the noise properties of the data are statistically stationary within each calendar year period; see Appendix E 1 for validation of this assumption. Therefore, we divide each calendar year of data X (m) (t) (with t entirely within year a) into many temporal 'chunks', each of duration τ , and treat each chunk as an independent noise realization [that is, we convert the ensemble average in Eq. (39) over hypothetical noise realizations to a straight average over chunks of the actual data, under the assumption of noise stationarity].
Since the length of calendar years varies between leap and non-leap years and we wish to use as much of our data as possible, we do not fix the length of τ universally, but instead choose a universal minimum value τ min , and divide each individual calendar year evenly into chunks whose durations exceed τ min . We choose the shortest such duration that allows us to evenly divide the entire year. Namely for a year of length T a , we use N chunks chunks of length τ , where and where the second expression assumes a unit of time measurement of minutes (i.e., the 'floor function' notation in the second expression is abused to mean 'round this result to the nearest minute'). Generically, computing the DFT of a time series of duration τ can be computationally difficult if the number of sample points in the duration τ is not a power of 2 (since the measurement cadence of SuperMAG data is ∆t = 1 min, this means that τ itself should be a power of 2 when measured in minutes). We therefore pad our time series x (m) with zeros to extend the number of data points in the chunk to the next power of 2 (i.e., we add additional values of x (m) (t j ) = 0 at assumed sample times t j = j∆t with j = N, . . . , 2 p − 1 for some p ∈ Z). 19 We therefore find it convenient to choose τ min to be a power of 2, and thus take the extended, padded chunk duration to be 2τ min . 20 The frequencies f p at which the DFTx (m) is computed will thus be multiples of (∆f ) = 1/(2τ min ). We find τ min = 16384 min = 2 14 min to be an adequate choice. [This implies τ = 16425 min for non-leap years and τ = 16470 min for leap years. Additionally, the DFT frequencies f p will be multiples of (∆f ) = 1/(32768 min) ∼ 5 × 10 −7 Hz.] We justify this choice in Appendix E 2, and show that our results do not depend strongly on the specific choice we have made.
For the i-th chunk of actual data X (m) in year a, we compute the quantity 21 and average over all M chunks within year a in order to estimate S a mn (f ): This process allows us to estimate S a mn (f p ) at the discrete frequencies f p = p(∆f ) for p ∈ Z. However, in the course of analyzing the data over durations longer than 2τ min , we will have access to a finer frequency spacing than (∆f ) , and so we really need access to S a mn (f ) sampled over this finer frequency range; since it is not possible to directly estimate S a mn (f ) on that finer grid with only our single data realization, our analysis interpolates the S a mn (f p ) estimated as at Eq. (42) to intermediate frequencies. Although this is approximate, there is no obvious superior approach.
Armed with the estimate Eq. (42) for the noise crosspower spectra S a mn (f p ), which yields the covariances be-tweenX (m) within a given year, we then compute the the padded data and the ensemble average of the PSD from the unpadded data agree statistically with their respective standard deviations of the mean. This step is purely for computational advantage. 20 Naive application of the definition of the PSD at Eq. (A6) taking the padded duration and padded number of data points yields the incorrect normalization for the desired PSD in this case because of the dead time associated with the padding. However, since we pad in such a way as to maintain the same ∆t in both the padded and un-padded data, the normalization of the FFT given at Eq. (A5) is correct, and the only modification we must make is to re-scale the PSD computed per Eq. (A6) by a factor of (2τ min )/τ ; cf. Eq. (41) and the comments in footnote 21. 21 The value of τ appearing in the denominator of Eq. (41) is actually taken to be τ ≡ N i data ∆t, where N i data is the number of data sampling points within chunk i for which at least one station has a valid measurement to allow the construction of x (m) (which necessarily is none of the points that have been padded with zeros). Generically, there is at least one station reporting at every time throughout the i-th chunk, and this procedure has no effect, yielding a value for τ that matches the value discussed in the main text; however, for the small number of cases where no stations happen to report data for some duration of the i-th chunk, τ as appearing in Eq. (41) is proportionally re-scaled to a smaller value.
covariances of the analysis variables X k as defined at Eq. (16). Suppose that N a k is the number of data points in the subseries X (m) k which were obtained in year a, so that a N a k ≡ ℵ where ℵ is the number of data points in the subseries X (m) k , then we have where T a k = N a k ∆t is the duration of time corresponding to the number of data samples in the subseries X (m) k in year a, assuming a measurement cadence of ∆t, such that in turn we have a T a k = T , the total duration of the k-th coherence time (except for the situations already noted in footnote 17, which are also handled appropriately here); see also Sec. V A 3 and the more detailed discussion in Sec. V E.
We may then write the covariance matrix for the X k schematically as for the appropriate values of m and n in the relevant locations; this matrix takes a block diagonal form because we assume the DFT results at distinct frequencies are uncorrelated variables, and X k is constructed in such a way that the successive blocks of entries all refer to the same frequency.

D. Bayesian statistical analysis
In the previous two subsections, we computed the expected X k under the signal hypothesis Eq. (1), and dis-cussed our data-driven noise estimation procedure. We can now synthesize these developments to construct a likelihood function for our model in terms of the expected signal vectors µ ik and the estimated covariance matrix Σ k . We can then use that likelihood function to construct the marginalized Bayesian posterior for ε given the data.

Likelihood function
Up to normalization, the likelihood function for the kth coherence time given the signal hypothesis Eq. (1) with a kinetic-mixing parameter ε is (we set the normalization factor for the likelihood to 1 arbitrarily) 22 where c k is the 3-vector with entries c ik for i = x, y, z [i.e., the variables defined in Eq. (38) which specify the arbitrary phase and spatial orientation of the DPDM polarization vector in the inertial frame for the k-th coherence time]. Assuming that all coherence times are treated as independent 'experiments', the full likelihood function L over all the available data will be taken to be the product 22 The normalization of the RHS of this equation (that is the normalization of ln L k , not L k ) cannot be chosen arbitrarily. It is set by demanding that X k X † k = Σ, or equivalently that Y k as defined by Eq.
Before proceeding to utilize this likelihood to construct a Bayesian posterior on ε, it will be advantageous and simplifying to make some changes of notation. Since Σ k is by construction a Hermitian, positive-definite matrix, it is possible to decompose it as Σ k = A k A † k , for some invertible A k . If we then define it can be shown that Eq. (45) can be expressed as Now, let N k be the 15 × 3 matrix whose first, second, and third columns take entries equal to the corresponding components of ν ik for i = x, y, z, respectively. Eq. (49) can then be rewritten as The singular value decomposition of N k can be written as 23 where U k is a 15×3 matrix with orthonormal columns [so that, specifically, U † k U k = 1 3 ], S k is a real 3 × 3 diagonal matrix, and V k is a 3 × 3 unitary matrix. We also define the 3-vector variables We can then re-write Eq. (50) as where to obtain the second expression we have expanded out, used U † k U k = 1 3 , added and subtracted |Z k | 2 , and simplified.
Our immediate goal now is to use Eq. (54) to define a likelihood function in terms of the variables Z k , which will be central to our analysis going forward.
To this end, consider the following preparatory argument. The matrix P k = U k U † k is an orthogonal projection operator: Armed with that knowledge, consider now the term in ( · · · )-brackets in Eq. (54). This term is (a) independent of the parameters ε and c k , and (b) equal to 23 Our convention is that of Ref. [73]; an alternative convention would take U k to be a square unitary matrix (here, 15 × 15), and S k to be rectangular diagonal (here, 15 × 3).
it thus does not depend on Z k . These observations imply, respectively, that (a ) the ( · · · )-term can simply be dropped from Eq. (54) in constructing a likelihood for ε and c k in terms of the Z k : where we dropped an additional irrelevant constant offset; and (b ) the resulting likelihood at Eq. (55) is still also interpretable in the usual way (again up to a constant offset) as the probability density for the Z k given the parameters, which we will see is necessary for our arguments in Sec. VI. 24 Physically, what has happened here is that the full 15dimensional analysis vectors X k that we constructed at Eq. (16) hold much more information about the measured magnetic fields than just the pieces necessary to find the signal Eq. (1), as is clear from the fact that the signal expectations X k are expressible as a sum over only three linearly independent vectors in the 15-dimensional space; see Eq. (38). What the foregoing mathematical manipulations have succeeded in identifying is the relevant part of the data X k to keep in the likelihood, Z k ; and the part that is superfluous to the signal search, B k .
The cognate full likelihood combining all the coherence times (assuming they are independent 'experiments') is given by

Marginalized likelihood function
Our goal is to construct the posterior distribution for ε in a Bayesian analysis framework. In constructing this 24 Indeed, for the purposes of Sec. V only, observation (a) would have sufficed. This is because Eq. (54) gives a likelihood, which is a function of parameters for fixed data, and we only use this in Sec. V to construct a marginalized posterior on ε in Eq. (63) below. Any term in Eq. (54) that is a function of the data only and independent of the parameters gives no useful information about those parameters, and constitutes a piece of the parameterindependent normalization constant for that marginalized posterior; but the structure of that parameter-independent normalization constant is irrelevant, since the posterior gets re-normalized to a give a unit integral. This leads to conclusion (a ). The reason that this argument is insufficient is that in Sec. VI we again interpret the likelihood Eq. (55) [or, really, the marginalized likelihood Eq. (60)] expressed in terms of the Z k as the probability density for Z k to be observed given theb parameters. Naturally, this is usually exactly what a likelihood like Eq. (45) is, by definition: L k (ε, c k X k ) ≡ α · p( X k |ε, c k ) with α a numerical constant. However, had we dropped a parameter-independent but Z k -dependent term in Eq. (54), we could no longer make the cognate identification for Eqs. (55) or (60). As such, it is important for the arguments in Sec. VI that observation (b) is true.
posterior however, we must account for the fact that our model for the DPDM field is such that we may not treat the d k simply as arbitrary model parameters which can be specified by us: instead, the statistical behavior of the DPDM field that emerges from the field being the sum of a large number of interfering plane waves (see discussion in Sec. II and Ref. [1]) dictates that the individual d k should themselves be treated as random variables that must be drawn from the appropriate distribution; see, e.g., Refs. [66,67,69,70,72]. Within the Bayesian framework, the appropriate procedure to fold that information into the ε posterior is to marginalize the combined likelihood Eq. (56) over the d k .
In this subsection we discuss the appropriate likelihood that describes the distribution of the d k , and then construct the marginalized combined likelihood.
Given the discussion in Sec. II [in particular Eq. (2)], and the definitions of the c k in Eqs. (18)-(20), both the real and imaginary parts of c k are independent normally distributed variables with mean zero which satisfy |c k | 2 = 1. Since V k is a unitary matrix, the same is true for the derived d k : |d k | 2 = 1. Therefore, the appropriate auxiliary likelihoods for the d k should be taken to be up to an irrelevant overall normalization. 25 The combined auxiliary likelihood for the d k is thus again assuming that the polarization vectors in distinct coherence times are independent random draws. The marginalized combined likelihood, defined as is thus given by (see Appendix D 1 for a detailed derivation) where z ik is the i-th component of Z k , and s ik is the diagonal (i, i)-element of the matrix S k (i.e., the i-th singular value of N k ); in both cases, i = 1, 2, 3. 25 The numerical factor of 3 in the exponent arises from assuming that the probability density function for each of the Re d i k and Im d i k for i = 1, 2, 3 takes the (common) form of a zero-mean normal distribution with unknown width, f (x) ∝ exp[−αx 2 ] for x = Re d 1 k , Im d 1 k , . . . , Im d 3 k , and then finding the value of α such that the normalization condition |d k | 2 = 1 is satisfied. See also footnote 22.

Priors and posteriors
Bayes' theorem constructs the marginalized posterior for ε, denoted by p(ε|{Z k }), from the marginalized likelihood given by Eq. (60), and the prior on ε, denoted by p(ε): We must thus specify a choice of prior on ε. Following Ref. [69], we will take the (reparametrization-invariant) objective Jeffreys prior [74] for ε; in a similar context, this choice of prior has the additional feature that it yields limits from a Bayesian analysis which are broadly in agreement with an alternative, frequentist approach [69].
The Jeffreys prior is defined formally in terms of the Fisher information matrix [74]; applying the formal definition, we show in Appendix D 2 that, for our analysis, this prior takes the form The posterior for ε is thus where N is a normalization factor. We can without loss of generality 26 restrict ε ≥ 0, and demand that N is set such that ∞ 0 dε p(ε|{Z k }) = 1. 27 With the appropriately normalized posterior, we can then set upper bounds on ε; for instance, the 95% credible upper limit (local significance)ε will be given by solving 26 Since the kinetic mixing term is the only term in the Lagrangian (see Ref. [1]) that is odd in A (in the interaction basis), a trivial field definition A → −A maps ε → −ε. Moreover, both the prior and posterior are even in ε. 27 Note that in the kinetically mixed basis in which L ⊃ − 1 4 F 2 − 1 4 (F ) 2 − 1 2 F F , there is a bound on | | < 1 for the physical region of parameter space that is smoothly connected to = 0; at = ±1, one or other of the two linear combinations F ±F becomes a non-propagating degree of freedom (i.e., the kinetic term vanishes). However, in the interaction basis we use in this work, we have ε = / √ 1 − 2 , so ε is unbounded above in the physical region of parameter space. Note however that our computation of the signal Eq. (1) is only valid for ε 1 [i.e., we have neglected terms at O(ε 2 )] [1]. However, our posterior distributions have little support for ε 1; our results are thus self-consistent.

E. Coherence time approximation and choice of frequencies
Our analysis to this point has been constructed to obtain a bound at a single frequency f A , but we have not yet specified how this frequency was chosen. One would ideally simply scan this frequency in some range. However, computationally we require the use of an FFT which can only evaluate bounds at discrete frequencies, and the specific set of frequencies depends on the duration of data we choose to analyze coherently. We discuss these issues further in this subsection.
In this subsection, we will be more precise about our usage of the term 'coherence time'; cf. footnote 13. Let T refer to the length of the data subseries analyzed in a coherent fashion under the analysis procedures thus far outlined in Sec. V, and denote by the shorter of the actual DPDM signal coherence time, and the total duration T tot = 48 yrs of SuperMAG data available for analysis.
For the following reasons, it has been implicit in our analysis construction to this point that T ≈ T coh (f A ): (1) beginning in Sec. V A 3 we split the SuperMAG data into subseries of length T , and assumed for the purposes of constructing the likelihood in Sec. V D that the polarization of the signal was constant for the duration of each [i.e., that the d k for each k were a single random draw from the expected distribution, Eq. (57)]. For that to be a consistent assumption, each signal subseries must not extend beyond a single actual DPDM coherence time, because the polarization wanders randomly on the latter timescale: T T coh (f A ); and (2) in Sec. V D, we explicitly constructed the joint likelihood over all the duration-T intervals by treating the polarization in each interval as having a distinct random orientation uncorrelated with that in the neighboring intervals, if any (i.e., each of the d k is a distinct random draw from the distribution defined by Eq. (57), uncorrelated with previous or future draws). Because we additionally analyse our data in contiguous blocks of duration T , this is only a good assumption if the duration of each block is long enough that, at the start of the subsequent block, the polarization has effectively been randomised by phase drifts. That latter time period is however again simply the definition of the DPDM coherence time, so we have: T T coh (f A ). 28 28 We note that if our analysis were over non-contiguous blocks, the criterion is simply that the start times of consecutive subseries are spaced by at least T coh (f A ), and not that the subseries' durations themselves must be at least T coh (f A ) long. However, we have mandated that there is no gap between consecutive subseries in order to maximize data usage, so the criterion for our analysis as constructed is indeed as stated in the text.
Since both T T coh (f A ) and T T coh (f A ) are needed, we have to take T ≈ T coh (f A ) for consistency. 29 Temporarily setting aside that T coh (f A ) itself is only known approximately (because the DM velocity profile is not known exactly, and the entire concept of the DPDM coherence time arises precisely because of the velocity dispersion in the interfering constituent plane waves), there is a computational problem in assuming T = T coh (f A ) exactly: it makes the duration of the signal to be analyzed an explicit function of the frequency at which the analysis is being performed [at least for all frequencies such that 1/(f A v 2 dm ) < T tot ]. That would preclude the application, necessary here owing to the multi-gigabyte volume of the full SuperMAG dataset, of the FFT algorithm to process the analysis of many frequencies simultaneously, because the FFT relies on having a fixedduration signal to transform. Having to either perform the slow DFT for each frequency, or indeed having to reperform the FFT for every frequency of interest, would be computationally prohibitive given available resources.
At a high-level, our solution to this computational issue seeks a trade-off between implementing the condition T ≈ T coh , 30 and still being able to exploit the computational speedup of the FFT algorithm to process many frequencies simultaneously. We will break up our full range of frequencies of interest into some small number of narrower frequency ranges, indexed by n, and perform our analysis for the frequencies in each such range n using a fixed duration of the data subseries, T = T n . We chose T n to be independent of frequency within each individual frequency range n (but varying for different n) such that T n ≈ T coh is satisfied, up to some fixed tolerance (which we take to be 3%), for all of the frequencies that lie within range n. Since T n is fixed within each frequency range n, we can then utilize the FFT algorithm to obtain results simultaneously for the whole set of FFT frequencies that lie within range n, {f ni } [where the f ni are multiples of 1/T n ; see below]. As we will not actually require too many different n (indeed we need only 56 such ranges) to cover our whole frequency range of interest in this way, this strategy allows us to construct results under the assumption that T ≈ T coh (f ) up to 29 Technically, without additional assumptions, only point (1) holds for the case where 1/(f A v 2 dm ) Ttot, such that T coh (f A ) = Ttot per Eq. (65). In that case, for point (2) to hold, the necessary assumption is simply that we wish to analyze all the available data to maximize the statistical power of the search; we implement this assumption by analyzing the whole dataset coherently with T = Ttot in this case. Since we define T coh (f ) such that max [T coh (f )] = Ttot, this case is automatically handled correctly by the discussion in the main text. 30 Note that, in some sense, this solution exploits the existing inherent uncertainty in the exact length of the DPDM coherence time to our advantage: we are not honor-bound to implement an inefficient analysis strategy to obtain T = T coh (f A ) exactly when the latter is only approximately known. We have some freedom to instead design an efficient analysis strategy that obtains T ≈ T coh (f A ). some known controllable tolerance, while also exploiting the FFT computational speedup, at only the modest cost of having to run the FFT algorithm 56 times. More precisely, we choose T n to be 31 where q = 0.03 fixes the aforementioned 3% tolerance, and the consecutive set of integers n = 0, . . . , 55 is chosen such that T n ranges from T tot down to approximately 10 6 minutes [i.e., the coherence time, assuming v dm = 10 −3 , corresponding to the sampling rate of the SuperMAG data, which is 1/(1 min)]. The set of frequencies {f ni } that we will consider to fall within range n will be f ni = i/T n for i = i min n , . . . , i max n . For n = 0, we take i min n = 10 6 /(1 + q) , while for the special case n = 0 (i.e., when the entire dataset is treated coherently), we have i min n = 0. The value of i max n is defined iteratively starting with the highest-frequency set, and for each n is taken to be the largest integer such that max [{f ni }] < min [{f n+1,i }]; this means that, approximately, i max n ≈ 10 6 (1+q) . 32 Because i max n must be iteratively constructed beginning with the set of frequencies containing the highest frequency, we must specify the highest frequency in the construction: this is taken to be one DFT frequency bin below the SuperMAG sampling frequency (i.e., twice the Nyquist frequency), such that i max 55 = T 55 /(1 min) − 1.
Defined this way, the individual sets of frequencies {f ni } cover non-overlapping ranges of frequencies. Moreover, for frequencies f ni such that T coh (f ni ) < T tot , we have Meanwhile it is easy to show that frequencies with T coh (f ni ) = T tot necessarily have n = 0 and so T coh (f ni ) = T n trivially by Eq. (66). Therefore, we indeed approximate the coherence time (or total data duration) to within a fixed percentage for every frequency f ni within every range n. 31 Let M be the number of data points corresponding to the time interval Tn. For computational purposes in the FFT, it is preferable for all the prime factors of M to be small. Therefore, we actually choose Tn such that M is the integer within 10 of the estimate implied by Eq. (66) that has the minimal largest prime factor. 32 While the FFT algorithm run on each duration-Tn dataset will also generally yield results for frequencies f ni with i outside the range shown in the text, for those frequencies the coherence time approximation tolerance will not be satisfied. We thus discard those results and utilize a different n for the construction of the results at the corresponding frequencies.  We show a graphical representation of this approximation scheme in Fig. 3.

F. Correction for finite signal width
Our analysis construction to this point has operated on the assumption that the dark-photon signal is exactly monochromatic within a coherence time, so that the entirety of the signal power appears in a single DFT frequency bin; in order words, we assumed exact coherence of the signal for a full coherence time T coh ∼ (f A v 2 dm ) −1 . Indeed, in the preceding subsection we matched the DFT frequency bin width to the coherence time to within 3% over the entire frequency range we consider in order to preserve this property. 33 However, this is a slight oversimplification of the situation: the DPDM signal is actually σ f ∼ 1/T coh wide in frequency space, so while we do expect the majority of the signal power to appear in the 33 For f A 6.4 × 10 −4 Hz, preservation of this property begins to fail because the signal coherence time begins to exceed (3% more than) the available data duration; see left edge of Fig. 3. As the frequency is decreased further, the coherence time further exceeds the data duration, and the signal therefore begins to become much narrower than a single DFT bin. This concentration of signal power in a single bin more closely matches our analysis construction, which implies that the degradation factor should be smoothly tapered to 1 (i.e., no degradation) for f A 6.4 × 10 −4 Hz. However, the lowest frequency that we explicitly present limits for in this work (see Fig. 4) is f A = 6 × 10 −4 Hz; at this frequency, the coherence time is still within 10% of the available data duration, and so we find it unnecessary to implement any such tapering of the degradation factor in presenting our results.
DFT bin corresponding to f A , some power will appear in the neighboring (few) bins as well. Given the way our analysis is constructed, if we did not account for this, we would set limits that are too aggressive.
While a more sophisticated approach to this analysis would have considered this spreading of the signal power from the beginning of the analysis construction, we leave such an improvement to future work. Instead, precisely because we have matched the DFT bin-width to the signal width to high accuracy over the whole frequency range, we can apply a simple frequency-independent rescaling factor to approximately correct for this in a post hoc fashion. That is, we proceed by simply degrading the limit on the kinetic mixing parameter from Eq. (64): with ζ > 1. In all of our results to follow, we present the degraded limitsε unless otherwise explicitly noted. It remains to estimate ζ. For the purposes of this estimate, we ignore the vectorial nature of the DPDM field, and focus only on the frequency-space spreading (this is equivalent to considering each vectorial component of the DPDM field independently); see Sec. VI C for the cognate signal injection that accounts for the vectorial nature of the signal and that validates this approach.
Assume that the DPDM field (component) is a sum of a large number of plane waves (see, e.g., Sec. II.A of Ref. [1]): where v n are samples drawn from an assumed galacticframe Maxwellian velocity distribution with an rms speed v 0 ∼ 10 −3 , and φ n is a random phase. Now analyse this field in the Fourier domain via the DFT, and let the largest single-bin value of the resulting PSD be S * . 34 Let the sum over the whole PSD of this field (i.e., the total signal power) be Σ S . Because our analysis is very roughly constructed so as to compare single-bin signal power to single-bin noise power, and because the PSD of the resulting magnetic field signal Eq. (1) arising from 34 Note that the average frequency of the DPDM field constructed in this fashion is 2πf A = m A √ 1 + v 2 ≈ m A (1+v 2 0 /2), which differs from the standard relationship we have employed to this point, 2πf A = m A , by a frequency shift of order (half) the DFT bin spacing. The correct way to interpret this shift is to identify the physical mean frequency of the DPDM field with the frequency at which we set limits, and consider this shift to be a modification to the relationship between f A and m A ; the correction is however negligible everywhere except for the frequency-mass identification. This procedure guarantees that the highest-power DFT bin is (except for fluctuations) the bin centered on f A . See the discussion in Sec. VI C. the DPDM is proportional to ε 2 , the appropriate degradation factor would then be ζ ∼ Σ S /S * . Averaging over 100 distinct random realizations of DPDM fields of this type constructed from sums of 2N + 1 = 5001 plane waves, we estimate numerically that the degradation factor would be ζ ∼ 1.24 (1). We therefore set ζ = 1.25 as the degradation factor. This 25% degradation factor is comparable to the uncertainties on many of our noise properties (see Appendix E), and so proceeding in this way is consistent with the overall accuracy of our full analysis.
We discuss this degradation factor further in Sec. VI C, where we verify that an injected signal would be correctly reconstructed.

G. Results
We now have in place all the relevant tools to set upper bounds on the kinetic mixing parameter ε; the results of our analysis are shown as the blue band in Fig. 4 as 95% credible upper limits (local significance) on ε as a function of the dark-photon mass m A . These constraints are complementary to the existing astrophysical limits also presented in Fig. 4, which arise from dark-photon heating of gas in the interstellar medium in the Milky Way (dotted orange) [46], the intergalactic medium around the time of helium reionization (short-dashed red) [48], and in the Leo T dwarf galaxy (dash-dotted purple) [49]; 35,36 or from dark-photon-photon conversion depleting dark matter (long-dashed green) [48]. Future bounds based on 21 cm observations are expected to become strong in this mass range [50], but we do not show current limits or projections here in light of the EDGES global 21 cm anomaly [78].
We note the existence of some clearly visible sharp peaks in the exclusion bounds shown in Fig. 4. These peaks are among some 30 naïve signal candidates that we identify in the data. We discuss these naïve signal candidates in detail in the next section, where we conclude that none of them clearly survive robustness checks on their consistency with the expected signal properties. Because we dismiss all these signal candidates, we can reasonably also plot in Fig. 4 as a guide to the eye a smoothed version of our limits (light blue solid line) that is obtained by averaging our limits over the ±25000 neighboring frequency bins. 35 Per Ref. [75], the limits in Ref. [49] are mildly weaker than those in the arXiv v1 and v2 preprints of that paper, on account of inter alia updated gas metallicity measurements of Leo T that were incorporated in the published version of Ref. [49]. 36 Limits similar to those in Ref. [49] appear also in Ref. [47]. The latter reference also gives a stronger preliminary bound based on a gas cloud of anomalously low (and disputed) temperature, which we do not show here; see the discussion in Refs. [47,76,77] and our comments in Ref. [1]. Exclusion bounds on the kinetic mixing parameter ε of the dark-photon dark matter as a function of the dark-matter mass m A (frequency f A ). The darker blue line (appearing as a band owing to frequency-to-frequency limit fluctuations) shows our 95% credible upper limit (local significance) [cf. Eqs. (64) and (70)] on the kinetic mixing parameter ε as a function of the dark-photon dark-matter mass (corresponding Compton frequency noted on upper axis), assuming that the dark photon constitutes all of the local dark-matter density, ρdm = 0.3 GeV/cm 3 , but taking into account the stochastic variations expected for classical-field dark matter (see, e.g., Refs. [66,67,69,70,72]). These limits include the effect of the 25% degradation factor discussed in Sec. V F. To guide the eye and give a sense of the relative density of stronger vs. weaker limits in narrow frequency bands, we also show as the lighter blue solid line the sliding average of the limit taken over the neighboring ±25000 frequencies.
Various sharply rising narrow spikes in our limits provide a variety of potential candidate signals; we examine these in detail in Sec. VI, where we conclude that none constitute robust evidence for a real signal. The various other lines show a variety of existing astrophysical limits arising from dark-photon dark-matter heating of gas in a number of astrophysical environments: the ionized interstellar medium in the Milky Way (dotted orange) [46]; the intergalactic medium around helium reionization (short-dashed red, labeled 'He ++ ') [48]; and gas in the Leo T dwarf galaxy (dot-dashed purple) [49]. A DM-depletion limit from nonresonant dark-photon-photon conversion [48] is also shown (long-dashed green, labeled '∆ρcdm'). Our limits are complementary to these existing bounds as they arise from terrestrial experimental data (analogous to 'direct detection'), and are thus subject to completely different sources of systematic uncertainty as compared to the other bounds shown (which are analogous to 'indirect detection').
Finally, we note that for m A 3 × 10 −17 eV our limits scale with increasing mass m A faster than m −1 A [cf. Eq. (1)], which is a manifestation of the decreasing noise at higher frequency in the SuperMAG magnetic field data. This trend is only terminated at the upper end of the plotted mass range owing to decreased sensitivity around and above the Nyquist frequency (f Nyq corresponds to a mass m A ∼ 3 × 10 −17 eV). This observation is highly encouraging because, assuming that this noise trend were to be maintained in the highercadence (i.e., one-second) SuperMAG data currently being released, it is plausible that this search method would allow access to kinetic-mixing parameter space at higher frequency that is currently unconstrained by astrophysical observations; we have not however undertaken any analysis of the higher-cadence data to check whether this is the case-this is deferred to future work. In any event, we note that even our existing limits are subject to completely independent systematics as compared to the existing astrophysical constraints in the mass range where we have presented limits in Fig. 4, and are already therefore complementary.

VI. CANDIDATES, VALIDATION, AND REJECTION
Our results in Sec. V G are phrased as exclusions (upper bounds) on the value of the parameter ε as a function of the DPDM mass. However, our analysis would be incomplete without also considering whether there are any indicia in the data of nonzero DPDM signals; indeed, this is the logical prior step. Even casual observation of Fig. 4 indicates the existence of multiple 'peaks' in the limits: frequencies at which the bounds are considerably weaker than those at neighboring frequencies. While this behavior may be the product of statistical fluctuations or other real non-DM-signal features in the data, it would also be expected behavior for the upper limit on ε to fluctuate upward for any specific frequency or frequencies at which a real DPDM signal(s) were present in the data (with a 'true' value of ε somewhat smaller than the value of the limit we have placed on ε for the respective frequencies).
In this section, we therefore complete our analysis by first developing in Sec. VI A formal criteria for identifying what we call 'naïve signal candidates': we find 30 such candidates in the data. Then, in Sec. VI B, we develop and apply tests to analyze whether or not the identified candidates are fully consistent with the expected properties of a DPDM signal, Eq. (1): on the basis of the discussion there, we conclude that none of the 30 naïve signal candidates can be considered robust evidence for a real DPDM signal in the SuperMAG data. Finally, we show in Sec. VI C that a mock signal of the form Eq. (1) injected into the (partially processed) SuperMAG data would be identified by our analysis, and not rejected by the robustness tests we develop, which validates our analysis approach. We offer discussion in Sec. VI D.

A. Naïve signal candidates
We begin by developing the formal criterion for declaring a feature in the data to be a naïve signal candidate.
As a first step, we must determine the statistical significance of any such feature under the zero-signal, null hypothesis: ε = 0. We work with the quantities z ik , defined in Sec. V D 2, whose (marginalized) likelihood for a given ε is given by Eq. (60).
From Eq. (60), we can see that under the null hypothesis (zero-signal; ε = 0), the real and imaginary parts of the quantities z ik are described by a zero-mean multivariate normal distribution. Therefore to determine which frequencies in our original analysis are inconsistent with the absence of a signal, we can simply compute the χ 2 statistic 37 and its 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 duration-T subseries into which we partitioned the full time series (i.e., the number of values through which the index k ranges; see Sec. V A 3). The relevant number of degrees of freedom is 6K 0 since each of the three Cartesian components of Z k has independent real and imaginary parts; note that we have made no correction for parameter estimation to the naïve number of degrees of freedom. Fig. 5 shows the value of p 0 for each frequency in the range of interest, as well as a histogram of all p 0 values taken over the whole range of frequencies we analyse. We consider a data feature at a certain frequency to be a naïve signal candidate (with 95% confidence) if p 0 is below the threshold p crit defined by where N f is the number of frequencies we consider in the range of interest; i.e., this threshold takes into account a trials factor, so the 95% confidence is global. For the frequency range of interest, 6 × 10 −4 Hz < f A < (1 min) −1 − 6 × 10 −4 Hz, 38 we have N f ∼ 3.3 × 10 6 , and the corresponding threshold is p crit = 1.6 × 10 −8 ; this threshold is shown as the horizontal (respectively, vertical) orange line in the left (right) panel of Fig. 5.
Using the criterion Eq. (74), we identify 30 naïve DPDM candidates in our frequency range of interest; see Tab. II. All of these naïve candidates are sufficiently narrow (i.e., they are only one-to-two frequency bins wide, consistent with ∆f /f A ∼ v 2 dm ∼ 10 −6 ) to be a potential DPDM signal. However, we cannot yet declare any of these naïve candidates to be a DPDM signal, as we must first verify that they pass further checks on their spatial and/or temporal characteristics.

B. Tests of candidates
Having identified 30 naïve signal candidates (see Tab. II) on the basis of the criterion specified in Sec. VI A, it is important to develop tests for the robustness of those naïve candidates. Candidates which fail these robustness checks can be rejected as being inconsistent with a real DPDM signal.
In particular, although the naïve signal candidates are indeed real magnetic field signal patterns in the data that have strong overlap with the VSH field pattern expected from a DPDM signal, it must also be the case that a DPDM signal should be present in all stations, and at all times. Therefore, re-analysis of any subdivision of the SuperMAG dataset should, for a real signal, yield parameter determinations consistent with the analysis of the full dataset. If, on the other hand, analysis of subdivisions of the full dataset yield inconsistent DPDM parameter determinations, that is strong evidence that the naïve signal candidate is not a real DM candidate; instead, this could be evidence for strong in-band local (in time or space) fluctuations driving the identification of the naïve candidate.
In this subsection, we develop these resampling tests and apply them to the 30 signal candidates we have identified.
Our first task is to identify the appropriate subdivisions of the full dataset to analyze independently for this resampling analysis. We perform two types of divisions of the data: by geographical location of the station, and by epoch of data acquisition. For the geographical division, we randomly partition the stations into four disjoint subsets. 39 For the temporal division, we divide the full dataset into four consecutive, non-overlapping temporal intervals, each with a duration of 12 years.
For each subset of data, we re-performed the analysis described in Sec. V for the frequencies corresponding to 39 Because we require each disjoint subset to have at least three active stations at all times in order to construct five linearly independent time series X (n) [j = 1, . . . , 5] for each subdivision, we are forced for this part of the analysis only to ignore the first six years of available data: insufficiently many stations are continuously operative during this time. Thus, for the geographical subdivisions in this resampling analysis, we analyse only the last 42 years of data available.
each of the 30 naïve signal candidates only. In performing this analysis, our choice of the time interval T described at length in Sec. V E is however inherited from the analysis of the full dataset instead of being re-adjusted on the basis of the subdivided data. Note however that this does mean that the number of duration-T intervals, K j , in each subdivision j of the data may be different from the number of such intervals in the analysis of the full dataset, K 0 . For the analysis of each division of the data j = 1, . . . , 8 (with j = 1, . . . , 4 being the temporal splits and j = 5, . . . , 8 being the geographical splits), this reanalysis procedure yields the quantities z ik,j and s ik,j ; these have the same definitions as the z ik and s ik in Sec. V, with the exception that they are evaluated only on the j-th subdivision of the data. For a kinetic mixing parameter of size ε, Eq. (60) indicates that the relevant test statistic to compute on the data in each subdivision j is the following χ 2 statistic: Given a fixed value of ε, it is easy to compute the p-value for this test statistic: it is simply found from the CDF of the χ 2 distribution as 40 where F χ 2 (6Kj ) is again the CDF for a χ 2 distribution with 6K j degrees of freedom [see below Eq. (73)]; we have once again made no correction for parameter estimation to the naïve number of degrees of freedom.
As we wish to test for consistency of the subdivisions of the data with analysis of the full dataset, we then weight these p-values by the posterior on ε computed in the analysis of the full dataset, p(ε|{z ik }) as defined at Eq. (63). That is, we assign to subdivision j the p-value p j = dε p(ε|{z ik }) · p j (ε).
We utilize Fisher's method [79][80][81][82] to combine the p j into a single p-value for the resampling checks. That is, we construct a joint test statistic over all the data subsamples by summing of the logarithms of the p j , and then compare it to a χ 2 -distribution. Our tests here must however be two-tailed as both large and small p j indicate disagreement with the original analysis. Therefore, the appropriate quantity whose logarithm must be summed is the minimum of p j and 1 − p j , rather than just p j . In other words, we will combine these p j into the single test statistic which has the corresponding joint p-value for n = 8 tests: with the relevant number of degrees of freedom in Fisher's method being 2n [79][80][81][82]. We additionally examine the p-values that arise from considering the temporal-only and geographical-only splits, p time and p geo , respectively; these are defined in the same fashion as p full , but with the appropriate restriction on the sum over the data subsets j in Eq. (78) where N f ≈ 3.3 × 10 6 ; see discussion below Eq. (74). Before proceeding to the interpretation of these pvalues, we note an important caveat. The procedure for the construction of p full (but not for p time or p geo ) given above is approximate in the following sense. Each data subset j = 1, . . . , 4 from the temporal split of the data is fully independent of each of the other distinct data subsets from the temporal split, and the same is true of the subsets j = 5, . . . , 8 from the geographic split. However, the temporal-split data subsets j = 1, . . . , 4 are not fully independent of the geographical-split data subsets j = 5, . . . , 8. For instance, the j = 1 data subset consists of data at all available stations for a fixed temporal duration, while the j = 5 data subset consists of data at all available times over a fixed subset of stations; this obviously implies that data from some station that is considered in the set of stations in j = 5 and that were taken during a time that is included in the temporal duration considered for j = 1, will appear in both the j = 1, 5 data subsets. If the data were exactly evenly distributed among data subsets, this would mean that there is an approximately 1/16-th overlap between each temporalspatial pair of data subsets; of course, the temporal split yields unevenly distributed data subsets (see Fig. 2), so this is only a rough estimate. Naturally, the overlap of data in the different data subsets correlates their p j values in some complicated way, although given the relatively modest data overlap that we estimate, we do not expect this effect to be large. A detailed accounting for this effect would, given the complexity of the SuperMAG dataset and the nontrivial analysis manipulations we apply to it, however likely require significant Monte Carlo modeling, which we consider to be beyond the scope of this work.
The preceding caveat notwithstanding, we proceed as follows. In the absence of the correlation in the data subsets, we would reject any naïve signal candidate with p full < 0.05 (95% confidence); because we expect that the correlation effect is mild, we will continue to automatically reject any naïve signal candidate with p full < 0.01 0.05. This automatically rejects 23 of the 30 candidates in Tab. II. For the remaining seven candidates that are not automatically rejected, we find that four have 0.01 < p full < 0.05 (candidates numbered 3, 8, 13, and 24 in Tab. II). We consider these to be in strong tension with the robustness checks; they are likely ruled    (77). Moreover, pj for j = 1, . . . , 8 that lie either close to 0 or close to 1 indicate disagreement with the original analysis. ptime (respectively, pgeo) is the combined significance of the temporal (geographical) split of the data in the resampling analysis. p full is the overall combined significance in the resampling analyses. We discuss these results at length in the text.
out, but we cannot give a definitive statement absent a quantitative accounting for the correlation caveat noted above. The remaining three candidates have p full > 0.05 (candidates numbered 6, 17, and 19 in Tab. II); these are not formally excluded, but they too have issues that put them in tension with an interpretation as a robust DPDM signal.
We discuss all seven of the candidates that are not automatically excluded in detail: Candidate 3 -This candidate is in strong tension with the combined temporal and geographical robustness test (p full = 0.026), but cannot be ruled out definitively on these grounds. It is however severely excluded by the temporal robustness test alone: p time = 2.3 × 10 −3 . We reject this candidate.
Candidates 8, 13, and 24 -All three of these candi-dates are also in strong tension with the combined temporal and geographical robustness test, but again cannot be definitively ruled out: 0.01 < p full < 0.05. Examining the temporal and geographical robustness tests separately, none of these candidates is clearly rejected by either alone: 0.05 < p time , p geo < 0.11. We note however that these candidates only exhibit moderate global significances: σ(p 0 ) = 2.6, 3.8, and 2.1, respectively. We consider these to be weak candidates which are heavily disfavored by our resampling analysis.
Candidate 6 -This candidate is perhaps the most interesting of the seven. It has a strong global significance, with σ(p 0 ) = 6.7, and is in good agreement both with the combined robustness test (p full = 0.27), and the geographical robustness test (p geo = 0.97). It is however in tension (but not definitively) with the temporal ro-bustness test: p time = 0.034. Because of this tension, we do not believe that there is clear or robust evidence that this candidate constitutes a signal; however it may warrant follow-up, for instance, in an analysis of the highercadence SuperMAG data.
Candidates 17 and 19 -While these candidates are in good agreement with the combined robustness test, as well as the separate temporal and geographical tests, we note that they lie exactly one DFT frequency bin above and below the Nyquist frequency for a one-minute data cadence. Although not dispositive, this constitutes good reason to believe that these peaks are systematic artefacts of the analysis. Future analysis of the highercadence data would settle this point definitively, as the Nyquist frequency for the one-second data would differ.
We conclude that none of the 30 naïve signal candidates should be considered to be a robust dark-photon dark-matter candidate signal, but that definitive exclusion of candidate 6 would require follow-up work.
Having already reached this conclusion, we have terminated our checks at this point; however, additional checks would have been possible had any candidate survived without demonstrating tension or inconsistency with the existing checks. For instance, we could also check whether the signal shows evidence of enough variation in the DPDM field from one coherence time to the next. The idea here would be to first perform parameter estimation to fix ε and m A for the signal using the full dataset. Then, we would examine subsets of the data with duration equal to the signal coherence time, and for each such data subset, we would estimate: (1) the DPDM polarization state, and (2) the DPDM field amplitude. We would then test whether the polarization state indeed randomizes on coherence-time timescales (i.e., whether the polarization state is too highly correlated between coherence times), and whether the DPDM field amplitude shows the appropriate distribution for a true DPDM field consisting of a sum of plane-waves with random phases; see, e.g., Eq. (81) below. Of course, because of the Earth's rotation, a re-orientation of the DPDM field polarization state in the galactic frame not only causes a change to the vectorial orientation of the magnetic field to be expected at each station, but it also shifts the relative power of the signal between the three frequencies f = f A , f A ± f d ; both effects are of course captured in Eq. (1), but one could additionally specifically test to ensure that the variation in the extracted DPDM polarization state is occasioned by the expected shift in the relative amount of the signal power at each of these frequencies, holding m A fixed.

C. Validation of analysis pipeline
In order to confirm that our analysis would identify and not reject a true DPDM signal in the data, we injected a mock physical DPDM signal with the expected spatial and temporal coherence properties into the Super-MAG dataset, and tested whether our analysis pipeline would correctly reconstruct it. We confirm that (a) the extracted limit curve on the kinetic mixing parameter ε is unaffected at frequencies other than in the expected vicinity of the injected signal frequency, (b) the limit on ε at the injection frequency is close to the value used to construct the mock signal (up to expected deviations owing to idealized assumptions in our pipeline), and (c) our resampling analysis that tests for signal robustness on the basis of spatial and temporal coherence of the signal correctly does not reject the injected mock signal candidates that appear at the injection frequency.
As we discussed in Ref. [1], the dark photon vector potential A can be considered to be a sum over a collection of plane waves with random vectorial orientations (although see also Ref. [51] for discussion on this point) and phase offsets, and frequencies centered on where v 0 ∼ 10 −3 is the DM velocity dispersion in the Milky Way (MW). Because the exact lineshape of this quasimonochromatic DPDM signal in the Fourier domain is actually not known as it depends on the detailed properties of the collection of plane waves being summed over (e.g., the exact velocity distribution in the MW), we choose to construct a straw-man artificial injected DPDM signal as follows.
Let c j (t) [j = x, y, z] denote the time-dependent orientation of the DPDM signal as defined in Eqs. (18)- (20) but with A m = A m (t) now time-dependent. Write c(t) as the 3-vector those j-th component is c j (t). We construct c(t) as follows: where N 1 is a large number of random plane waves (see Tab. III);â n andb n are randomly selected real unit 3-vectors specifying the polarization state of each such plane wave; φ n is a random phase; f A is the physical average frequency of the signal; v n is a 3-velocity selected from an isotropic Gaussian velocity distribution with a single-velocity-component standard deviation 3 such that the root-mean-square DM speed is v dm = 10 −3 (note: we neglect the difference between the Earth/Solar System and galactic frames for this purposes of this construction); and ξ ≈ 1 − v 2 dm /2 + O(v 4 dm ) is an otherwise-negligible correction factor to the massfrequency relationship chosen such that ξ √ 1 + v 2 ≡ 1, with the average taken over the velocity distribution (i.e., we inject a signal with a slightly corrected mass m A = 2πξf A such that f A is the average signal frequency; this is done for the technical reason that we wish to inject the signal on average at an exact DFT frequency, and not shifted upward by order of a DFT bin half-width-see also footnote 34). These c(t) have the appropriate normalization (i.e., |c(t)| 2 = 1 with the TABLE III. Parameters used to generate the mock signal injected into the SuperMAG data for the purposes of analysis pipeline validation; parameters are defined in detail in the text (see Sec. VI C). Also shown for comparison is the value of ∆f , the DFT frequency spacing corresponding to the value of the approximate coherence time used in the vicinity of f * (see discussion in Sec. V E).

Parameter
Symbol Value Central frequency f * 7.5 × 10 −3 Hz Signal 'width' (≡ f * v 2 dm ) σ f 7.5 × 10 −9 Hz Kinetic mixing parameter ε * 10 −3 Number of plane waves summed 2N + 1 1001 DFT frequency spacing at f = f * ∆f 7.45 × 10 −9 Hz average being an ensemble average over random phases, orientations, and velocities; or, equivalently, a temporal average over times much longer than the coherence time) and coherence properties for a realistic DPDM signal (technically, in the galactic rest-frame), assuming a non-truncated Standard Halo Model (see, e.g., Ref. [83]) for the DM velocity distribution. Note that the injected signal Eq. (81) has a frequency-space 'width' of order σ f ≡ f A v 2 dm . As in Sec. V B, we can express the effect of a DPDM signal on the data time series X (n) (t j ) [as defined in Eqs. (7)-(11)] using various combinations of the time series H (n) (t j ) [as defined in Eqs. (28)- (34)], weighted by the time-dependent orientations c i (t j ). Specifically, we inject our mock DPDM signal into the partiallyprocessed SuperMAG time-series dataset by making the following substitution: with t denoting transpose; see Eq. (36) and Eqs. (C1)-(C2) for similar expressions for the Fourier transform of the signal. We refer to the dataset consisting of the SuperMAG data plus this injected mock signal as the 'mock-signal dataset'; the injected signal parameters are shown in Tab. III along with other relevant data. We rerun the full analysis detailed in Secs. V, VI A, and VI B on the mock-signal dataset. The exclusion bounds resulting from this analysis are shown in red in Fig. 6, with our exclusion bound from Fig. 4 derived from the SuperMAG dataset superimposed in blue. As expected, the exclusion bounds derived from the mock-signal dataset and the unadulterated SuperMAG dataset agree well everywhere except in the vicinity of the injected frequency, 41 and at its reflection across the Nyquist frequency (1 min) −1 − f * ; there are strong narrowband peaks in the exclusion bounds at those two frequencies, which are diagnostic of a signal.
We note that the exclusion bound which is set at the injected frequency actually appears slightly weaker than the value of the kinetic mixing parameter used to construct the signal; cf. Tab. III and Fig. 6. This is of course the expected behavior for an upper limit in the presence of a signal and additive noise. Note however that the degradation factor of ζ = 1.25 that was discussed in Sec. V F proves crucial in obtaining this result; indeed, one can clearly see in the inset plot of Fig. 6 that injected signal power has appeared in more than one DFT bin, in line with our arguments in Sec. V F. The limit set at the Nyquist reflection is slightly stronger than the injected signal value of ε * , but this is traceable to the fact that the Nyquist reflection of the signal appears not at an exact DFT frequency bin, and so more signal power leaks to neighboring bins than is the case for the signal at the 41 The attentive reader will note from the inset axes in Fig. 6 that even at some significant distance (as compared to σ f ) from the strong signal peaks in the mock-signal dataset (but still in its vicinity), the exclusion bounds from the two datasets do not agree exactly. This is because the noise level that is used to set the exclusion bounds at any one frequency is estimated in a datadriven way using all the data, including any injected signal, in some nearby frequency range; see Sec. V C and Appendix E. As a result, the bounds away from the region f = f * ±(few)·σ f that are derived using the mock-signal dataset are weakened slightly as compared to the unadulterated dataset, owing to the higher noise estimate in the former. A refinement of this approach would be possible, but we note from Fig. 6 that even an injected signal with a very high signal-to-noise ratio (SNR) weakens the exclusion bounds in the vicinity of that signal by a factor of only O(1). See also footnote 18.
injected frequency; we have not degraded the limits so significantly as to account for this effect, as our limits are to be formally interpreted as correct only for injected signals that lie at exact DFT frequencies.
As a separate check, we also verified that when injecting into our analysis pipeline an idealized exactly monochromatic signal (at an exact DFT frequency, so that spectral leakage effects can be ignored; see, e.g., Refs. [84,85]) with A aligned to the Earth's rotational axis, the limit on ε correctly appears slightly weaker than the injected signal size, even without the ζ = 1.25 degradation factor applied.
Taken together, these checks confirm that the analysis operates correctly for the exactly monochromatic signal it is formally constructed to search for without any post hoc correction, and that the amplitude of the post hoc correction applied in Sec. V F for a real signal with the appropriate frequency-space width is of an appropriate magnitude.
We also point out that our validation here has been phrased entirely in terms of exclusion bounds; formally, we should perform parameter estimation on the mock signal to estimate the recovered frequency and kinetic mixing parameter. However, since the injected signal is assumed to have a large SNR in our tests, that level of detail in this validation analysis would be unwarranted.
Additionally, we reran the resampling analysis detailed in Secs. VI A and VI B on the mock-signal dataset to ensure that we both identified candidate peaks at the injected signal frequency, and did not reject them on the basis of the robustness tests. All 30 of the alreadydiscussed candidate peaks listed in Tab. II appear again when we analyse the mock-signal dataset, with the same p 0 -values as in the original analysis. We also identify 15 additional naïve signal candidate peaks: two of these are clustered at f * + f d , five at f * , and one at f * − f d ; the remaining seven peaks appear clustered around the reflections of these peaks across the Nyquist frequency with two at (1 min) −1 − (f * + f d ) and five at (1 min) −1 − f * . Our resampling analysis correctly rejects (p full < 0.01) all five of the naïve signal candidates at f * ± f d and (1 min) −1 − (f * ± f d ), along with one at f * . All the other naïve signal candidates that appear around f * and its Nyquist reflection at (1 min) −1 − f * are not rejected, indicating that the resampling analysis correctly does not rule out potential signals with the correct spatial and temporal properties. We stress that this non-rejection of multiple naïve signal candidates around f * (and around its reflection through the Nyquist frequency) should not be read to indicate that our pipeline identified some unexpected multi-peaked structure; rather, this is an expected result given that σ f ∼ ∆f and the SNR of the injected signal is large. Specifically, the non-rejected candidates all appear in contiguous ranges of DFT frequencies; see the inset of Fig. 6. It is to be expected that, evaluated on a bin-by-bin basis, some number of consecutive DFT bins (with widths similar to the signal width parameter) in the immediate vicinity of a large injected signal would each contain sufficient signal power in their own right to qualify as candidates and not be rejected, since they have the appropriate signal properties. The results in this subsection validate that our analysis pipeline functions as expected.

D. Discussion
In this section, we examined in detail the 30 naïve signal candidates that we identified in the data that exceed a 95% confidence global significance threshold. Applying a resampling analysis to temporal and spatial subsets of the data to test for the robustness of these candidates and their consistency with the expected persistent global nature of the signal, we found that we could automatically reject 23 of the 30 candidates on the grounds of their failing by a large margin the combined temporal and spatial checks. Of the remaining seven candidates, one fails the temporal check severely; three weak candidates are in strong tension with the combined robustness test (but cannot be definitely ruled out owing to the unaccountedfor correlation issue discussed in Sec. VI B); two pass all formal robustness tests and are reasonably globally significant, but they appear respectively one DFT frequency bin above or below the Nyquist sampling frequency, and must therefore be viewed with appropriate skepticism with regard to analysis systematics; and the final candidate is statistically significant globally, and passes the combined test and the spatial test, but is in strong tension with the temporal test. As such, we do not consider any of these naïve signal candidates to be robust candidates for a real DPDM signal in the data on the basis of the analysis of the one-minute SuperMAG dataset presented in this work. Nevertheless, it would be worthwhile for future work to perform an analysis of the higher-cadence SuperMAG data, as this would likely definitively settle some questions regarding those candidates that were not automatically rejected on the basis of the formal statistical criteria we applied here.
In this section we also presented a validation of our analysis pipeline by showing that a fake signal injected at the level of the X (n) variables defined at Eq. (12) is (a) recovered by the analysis at (b) the appropriate frequency with (c) an appropriate value of the kinetic mixing parameter (up to expected deviations; see discussion in Sec. VI C), and that (d) this injected signal survived the spatial and temporal robustness checks we apply to naïve signal candidates. This verifies that our analysis pipeline performs as expected, and would correctly identify and not reject a real signal in the data.
In conclusion, we find no robust statistical evidence for the existence of the DPDM signal in the SuperMAG data, and we are confident that our analysis pipeline would have identified such a signal had it been present.

VII. CONCLUSION
In this work, we presented the details of our analysis of the one-minute-cadence SuperMAG geomagnetic field dataset [43][44][45] for the quasi-monochromatic (fractional linewidth σ f /f ∼ 10 −6 ) global magnetic field signal of dark-photon dark matter that we recently proposed in a companion paper [1]. Because the size of the magnetic field signal is B ∝ εm A R √ ρ dm , suffering only a geometrical suppression by the radius of the Earth R, we were able to place competitive limits on the parameter space for kinetically mixed dark-photon dark matter. These limits are shown in Fig. 4, and cover the mass range 2 × 10 −18 eV m A 7 × 10 −17 eV corresponding to frequencies 6×10 −4 Hz f A 2×10 −2 Hz, with the upper end of this mass reach being limited by the one-minute sampling cadence of the SuperMAG data. Our analysis made use of a Bayesian framework (using a Jeffreys prior on our kinetic mixing parameter ε) [69,72] in order to incorporate the effects of the statistically varying local dark-photon dark-matter field amplitude, which assumes values of A such that ρ dm ∼ 0.3 GeV/cm 3 only on average; this is a conservative value for the DM density (recent ADMX limits [86] assume a value 50% larger; see also Ref. [87]).
In the course of our search, we initially identified 30 naïve signal candidates for dark-photon dark matter that exceeded a global 95% confidence threshold on the basis of our main analysis. In order to ensure that we did not miss a signal, we performed robustness cross-checks on these candidates, testing them for both the spatial and temporal coherence characteristics we expected from our signal. On the basis of these cross-checks and other indicia, we concluded that none of these signal candidates provide robust evidence for the existence of a DPDM signal in the data.
We also verified our analysis pipeline by injecting both (a) an exactly monochromatic signal aligned along the Earth's rotational axis and appearing exactly at one of the DFT frequency values where we set limits, and (b) a signal with the correct statistical properties for a dark-photon dark-matter signal whose phase-and polarization-coherence are fixed by DM velocity dispersion (i.e., the σ f /f * ∼ 10 −6 linewidth of the signal). In the former case, we correctly recovered a limit at the expected value of the kinetic mixing parameter without any post hoc degradation to our results. In the latter case, without any correction factor applied, we would have recovered limits slightly stronger than the injected signal owing to a mild violation of the assumptions used in the construction of the analysis by a signal with the full coherence properties of the dark-photon dark-matter signal (e.g., our analysis assumed an exactly monochromatic signal with exact phase coherence for a full coherence time, which in reality is only approximately true as the phase and polarization of the signal actually evolves fractionally by O(1) over a full coherence time). We applied a 25% degradation to our limits to compensate for this in a post hoc fashion. The resulting degraded limits correctly fail to exclude the injected signal, as expected. We also verified in both cases that the temporal-spatial robustness checks we used to dismiss naïve signal candidates in the real data, did not dismiss these injected signals.
The dark-photon dark-matter exclusion bounds we set in Fig. 4 are complementary to existing astrophysical bounds in the same mass range that arise from gas heating in various astrophysical settings [46,48,49], and a DM-depletion bound arising from nonresonant darkphoton-photon conversion [48]. Importantly, as we noted in Ref. [1], the scaling of our limits with mass is steeper than m −1 A owing to falling noise in the SuperMAG dataset as a function of increasing frequency; this raises the prospect, assuming that this falling noise trend continues to hold, that higher-cadence magnetic data would allow the search for this signal to access currently un-constrained dark-photon dark-matter parameter space. SuperMAG is currently in the process of releasing onesecond-cadence data, and we defer analysis of those data to future work (such an analysis would also provide a separate check on our dismissal of the naïve signal candidates that we identified near the Nyquist frequency in the one-minute-cadence data).
Additionally, we set our limits under the assumption of a Standard Halo Model (see, e.g., Ref. [83]) DM velocity abundance and dispersion. If stream-like structures actually dominate the local DM abundance (see, e.g., Refs. [88][89][90][91][92][93]), then the signal would be narrower and the DM abundance increased; these effects would make a signal more easily discernible in the data.
Finally, we note that while the analysis approach we presented here exploits the power of this large dataset well, it is not necessarily optimal. We leave to future work refinement and optimization of the analysis.
In this appendix, we give our conventions for the continuous Fourier Transform (FT) and Discrete Fourier Transform (DFT).
The continuous FTF (ω) [alternatively,F (f ) with ω = 2πf ] of a continuous signal F (t) is defined by The that are equally spaced and taken with a cadence ∆t for a total duration T ≡ N ∆t is defined by where t n ≡ n∆t = nT /N , and f k ≡ k∆f ≡ k/T . The (discrete) two-sided power spectral density (PSD) is defined in terms of the DFT: Note that a monochromatic signal with frequency has a DFT given bŷ where k = 0, . . . , N −1; the corresponding two-sided PSD isŜ (A9) we have where at Eq. (D5) we completed the square in the exponent to obtain Gaussian integrals, and simplified. Up to an arbitrary normalization, Eq. (60) follows.

Jeffreys prior
In this subsection, we derive the Jeffreys prior, Eq. (62).
The Fisher information matrix, I is defined as [74] where Θ is the parameter vector, E · · · | Θ is the expectation value over data realizations drawn assuming the values of the parameters Θ, and L = L Θ {x} is the likelihood considered as a function of the model parameters given the data realization {x}. The Jeffreys prior is the unique reparametrizationinvariant prior, and is defined in terms of I [74]: For the case of a one-dimensional parameter vector, as is our case after marginalizing over the d k , this simplifies: for notational simplicity, we leave the '|ε' implicit in what follows. Therefore, Interpreting the likelihood L given by Eq. (60) via its definition as the probability density function of the data given the model parameter ε (see discussion at footnote 24), we see that (when viewed as random variables rather than as the specific data realizations we have) the real and imaginary parts of the z ik are all independent, zero-mean normally distributed variables satisfying |z ik | 2 = 1 + ε 2 s 2 ik /3. Therefore the expectation values of all the cross terms (i.e., those with differing i and k) in Eq. (D11) vanish since they factor into two quantities, each with expectation value zero. The only remaining terms are thus those where i and k are the same for both factors. Therefore,

Noise variation within a calendar year
Our analysis in Sec. V C computed noise spectra S a mn (f p ) associated with particular calendar years, under the assumption that the noise remained statistically stationary within a calendar year. We evaluate this assumption by computing the same quantity S a mn (f p ) on timescales shorter than a full year and comparing the results to the full-year estimate. In particular, we divide a given year evenly into four quarters and recompute S a mn (f p ) as in Sec. V C using the data from each quarter (again taking τ min = 16384 min). Due to the finite number of samples used in Eq. (42), the estimate S a mn (f p ) that we compute has a large variance from one frequency to the next. For the purposes of this comparison it is useful to examine and compare the moving average of S a mn (f p ) taken over a range of frequency bins in a sliding window centered on each frequency. That is, we computē where we take the window half-width to be w = 512, corresponding to a top-hat sliding window with width 5.2 × 10 −4 Hz. Likewise it is useful to compute the standard deviation of the values of S a mn (f p ) within the sliding window, as a statistic to quantify the spread: In Figs. 8 and 9, we compare the quarterly estimates of S a mn (f p ) to the full-year moving averageS a mn (f p ), for diagonal (m = n) and off-diagonal (m = n) elements, respectively. Specifically, in the top panels of Fig. 8 appear in Eq. (E1) for a few representative frequencies.
In Fig. 9, on the other hand, we show scatter plots of the quarterly estimates S a mn (f p ) that appear in Eq. (E1), along with their corresponding 68% coverage ellipses; the full-year averageS a mn (f p ) is marked by a black cross. The variety of years, components m, n, and frequencies f p displayed in Fig. 9 are broadly representative. It is clear that for a wide range of frequencies and component choices (both diagonal and off-diagonal), the full-year average is consistent with the distribution of quarterly estimates, indicating (within the precision of the shorter-timescale estimates) that the assumption of statistical stationarity is satisfied. We do however note that there is variation in the degree to which the individual quarterly results are consistent with each other within a year [e.g., the 1989 results for the (m, n) = (5, 5) component show some mild tension between the first and third quarters; whereas, e.g., the 2007 results for the (m, n) = (4, 2) component are in better agreement]. It is possible that a more sophisticated analysis than that presented here could account for this.

Choice of τ min
Our analysis in Sec. V C has one arbitrary parameter: τ min , the lower bound for the duration of the chunks of data that are used in our noise analysis. There is a trade-off in the selection of this parameter: it should be as small as possible to allow for more independent chunks, and thus a better overall statistical estimate of the noise. On the other hand, taking the chunks length too short would increase the relative impact of possible correlations between data points near the edges of consecutive chunks, which could systematically bias our noise estimate as the chunks would no longer be sufficiently statistically independent. It is therefore important to understand the dependence of S a mn (f p ) on τ min , and to select τ min large enough that the dependence of S a mn (f p ) on that parameter becomes subdominant to the uncertainty on the estimate of S a mn (f p ) itself. To quantify this, we re-compute S a mn (f p ) as in Sec. V C, but using several different choices of τ min . For each choice of τ min , we compute the frequency-space moving average spectrum,S a mn (f p ), as in Eq. (E1). The results are shown in Fig. 10: the panels on the left of Fig. 10 show the full spectraS a mn (f p ) obtained using values of τ min = 2 j min for various integers j = 11, . . . , 15, again for various representative choices of the year a and components m, n. The panels on the right in Fig. 10 show the dependence ofS a mn (f p ) on τ min for the same years and components at three particular representative frequencies chosen in the low-, mid-, and high-frequency regions of our range of interest, along with their associated standard deviations (shaded bands), σ a mn (f p ) [defined in Eq. (E2)].
The figure demonstrates that for well-behaved components (typically the larger, diagonal m = n components in most years), the dependence ofS a mn (f p ) on τ min is smaller than the spread σ a mn (f p ) in the values of S a mn (f p ), for τ min = 16384 min or greater. The less well-behaved components that vary somewhat more with changing τ min are typically off-diagonal (i.e., m = n) components, which tend to be smaller by a factor of O(3-10) than the on-diagonal components; these thus have less impact on our results, regardless of the choice of τ min . Generically, the picture is that, for most years, a choice of τ min = 16384 min is reasonable, and gives results that do not strongly depend on τ min ; this choice also allows sufficiently many independent chunks to obtain a noise estimate within a O(10%), which is sufficient for the level of precision in our analysis.
Nevertheless, we do note that for certain specific years, there are exceptions to this generic picture. For instance, Fig. 10 shows an example of a less well-behaved offdiagonal component [(m, n) = (2, 1)] from an early year (a = 1976), when much fewer stations were active as compared to later years; cf. Fig. 2. This result clearly shows a sharp change in behavior between τ min = 8192 min and τ min = 16384 min that might be worrisome. However, anomalous behavior in years with very few stations running would not be entirely unexpected, as a single station turning on/off would have a larger impact on the overall results than in years with many stations running. While we still include results from such 'anomalous' years in our analysis, we note that-precisely because they have much fewer stations running-they will be naturally down-weighted in our analysis owing to the larger noise obtained with fewer stations running, and will thus have little influence on our final results.
FIG. 9. Noise stationarity validation; off-diagonal elements of S a mn (f ). Because the off-diagonal components of the Hermitian matrices S a mn (f ) are complex, we cannot show the quarter-by-quarter agreement as a function of frequency as simply for the off-diagonal elements as we did in Fig. 8 for the real, on-diagonal components. In this figure, at three selected representative frequencies (the same ones indicated by the vertical grey dashed lines in Fig. 8), we select some representative off-diagonal components (m, n) for some representative years a, and show scatter plots in the complex plane of the values of S a mn (fq) for the fq that lie within the corresponding averaging windows used to determine the quarterly values ofS a mn (f ) [see discussion around Eq. (E1)] (shaded colored circular points; see legend). Note that the density of points, and not the depth of shading, indicates the clustering of values within each quarter (with some loss of resolution in the denser regions); we have fixed the depth of shading for each quarter to be independent of the density of points in order to make the differences in the clustering of points from quarter to quarter clearer. Also shown are the 68% coverage ellipses for two-dimensional Gaussian fits to the scattered points for each quarter (like-colored solid ellipses; see legend), along with the full-year average value ofS a mn (f ) (black cross). The left-hand column shows a case where there is some mild tension (within factors of ∼ 2-3) between the intra-year statistical stationary of the noise assumed in our analysis, and the realised noise (i.e., the fitted 68% coverage ellipses for different quarters vary somewhat), while the right-hand column shows a case where intra-year noise stationarity is realised well.  (m, n). The colors as annotated in the legend on the lower left distinguish the various cases. Note that the sliding window used to compute the frequency-space average is varied (w = 2 j−5 for the integers j defined above), so that the sliding window width used to compute the frequency-space averages is maintained at 5.2 × 10 −4 Hz in all cases. Note also that, in this plot, where relevant and in contrast to Fig. 9, we show the absolute value of the complex off-diagonalS a mn (f ). Right column: Dependence of |S a mn (f )| on τmin for three representative fixed frequencies: f = 1.04 × 10 −3 Hz, f = 2.95 × 10 −3 Hz, and f = 8.33 × 10 −3 Hz (the locations of these frequencies are also indicated as dashed vertical grey lines in the left panels). The grey bands represent the spread in the values in the sliding windows, σ a mn (f ); see Eq. (E2). It is clear that, in most cases, for τmin ≥ 16384 min, the dependence of |S a mn (f )| on τmin is much smaller than the spread σ a mn (f ). Note that the top row shows an anomalous off-diagonal component from an early year (a = 1976), whose behavior is likely strongly influenced by the small number of active stations. Additionally, the panels in the third row show a result for an off-diagonal component [(m, n) = (3, 1)], which is somewhat smaller in normalization than the other on-diagonal components from the same year; residual variation of these results with τmin will thus have little effect on our results. . This is realised (note that the vertical axis covers only a small range of values) for |z ik | 2 to within ∼ 8% and for |z ik | 4 /(2 |z ik | 2 2 ) to within ∼ 2%, over the entire frequency range of interest.

Gaussianity of variables
In Sec. V D, we constructed a likelihood function for our Bayesian analysis under the assumption of Gaussianity on the real and imaginary components of the z ik ; see, e.g., Eq. (60). We evaluate the validity of this assumption by computing the four-point function of z ik and comparing it to the expected result assuming Gaussianity: |z ik | 4 Gaussian = 2 |z ik | 2 2 . We do this at each frequency by treating the z ik that we compute for each value of i and k as independent samples of Eq. (60) in the absence of a signal (ε = 0), also assuming that values of z ik at different frequencies are independent. 43 We again take moving averages in frequency-space of these independent samples and study the result as a function of frequency. That is, we compute for m = 2, 4, where f p runs over the full set {f ni } described in Sec. V E, z ik (f p ) are our analysis variables calculated for f A = f p , and K(f p ) is the number of subseries [cf. the definition of K in Sec. V A 3] used in the analysis for f A = f p (which will vary for f ni of different n). Here, we use w = 25000.
In Fig. 11 we show in orange the ratio of |z ik | 4 /2 |z ik | 2 2 , which is expected to be 1 in the case of exact Gaussianity, as a function of frequency. For comparison, we also show |z ik | 2 in blue, which per Eq. (60) should be 1 under our analysis assumptions. Deviations of |z ik | 2 from 1 would potentially stem from mis-estimation of the noise S a mn due to the uncertainties σ a mn referenced in the above subsections. It is clear from Fig. 11 that the deviation from Gaussianity is smaller than the mis-estimation of the noise. There is a small deviation of |z ik | 2 = 1, within the 10% level, across the entire frequency range; on the other hand, the assumption of Gaussianity as tested in this fashion is good to within 2%. These levels are acceptably accurate for the purposes of our analysis.