Search for transient ultralight dark matter signatures with networks of precision measurement devices using a Bayesian statistics method
B. M. Roberts G. Blewitt C. Dailey A. Derevianko
 









	Department of Physics, University of Nevada, Reno, 89557, USA

	Department of Physics, University of Nevada, Reno, 89557, USA 

	Department of Physics, University of Nevada, Reno, 89557, USA



  



We analyze the prospects of employing a distributed global network of precision measurement devices as a dark matter and exotic physics observatory. 
In particular, we consider the atomic clocks of the Global Positioning System (GPS), consisting of a constellation of 32 medium-Earth orbit satellites equipped with either Cs or Rb microwave clocks and a number of Earth-based receiver stations, some of which employ highly-stable H-maser atomic clocks.
High-accuracy timing data is available for almost two decades.
By analyzing the satellite and terrestrial atomic clock data, it is possible to search for transient signatures of exotic physics, such as “clumpy” dark matter and dark energy, effectively transforming the GPS constellation into a 50,000 km aperture sensor array.
Here we characterize the noise of the GPS satellite atomic clocks, describe the search method based on Bayesian statistics, and test the method using simulated clock data.
We present the projected discovery reach using our method, and demonstrate that it can surpass the existing constrains by several order of magnitude for certain models.
Our method is not limited in scope to GPS or atomic clock networks, and can also be applied to other networks of precision measurement devices.




Search for transient ultralight dark matter signatures with networks of precision measurement devices using a Bayesian statistics method
    A. Derevianko
     December 30, 2023
========================================================================================================================================








§ INTRODUCTION


Astrophysical observations suggest that ordinary luminous and baryonic matter contributes only about 5% to the total energy density of the Universe, with the rest due to dark matter (DM)  at ∼ 25%, and dark energy (DE) at ∼ 70%.
Despite the overwhelming cosmological evidence for the existence of DM, and the considerable effort of the scientific community over several decades, there is as of yet no definitive evidence for DM in terrestrial experiments.

Currently all the evidence for DM comes from observations carried out over distances greater than or comparable to  galactic  scales <cit.>.  
In general, in order to perform  a direct DM detection experiment, these vast 10 kpc (∼10^13m)  distances must be extrapolated down to scales that are accessible in laboratory settings (∼ 1 m). 
This extrapolation leads to a  variety  of plausible theoretical possibilities for DM models, ranging from elementary particles to black holes. 
Considering the broad variety of models and the associated assortment of non-gravitational interactions of DM with ordinary matter, it is important to constrain DM models by creatively reanalyzing archival data <cit.>.
Compared to investments into dedicated experiments, this is a relatively low-cost strategy with potential for important discovery.
Here we develop a method based on Bayesian statistics for a time-domain DM search using data accumulated by networks of precision measurements devices.

The field of low-energy precision measurements has proven to be an important area for probing fundamental laws and searching for new physics, that is often complementary to collider experiments <cit.>.
The idea of using a distributed network of precision measurement devices to search for DM and other exotic physics signatures is one promising approach <cit.>.
The particular network considered here is the Global Positioning System (GPS), a satellite constellation of nominally 32 satellites in medium-Earth orbit (altitude ∼20,000 km) housing atomic clocks, as well as a large number of atomic clocks on ground-based receiver stations. 
Following the proposal of Ref. <cit.>, we use the GPS constellation as a ∼50,000km aperture sensor array, analyzing the satellite and terrestrial atomic clock data for transient signatures of exotic physics, such as DM and DE.
High-quality timing data from the GPS network exists for the past 18 years, and is made freely available by, e.g., the Jet Propulsion Laboratory (JPL), NASA <cit.>. This dataset is routinely augmented with more recent data.



The global scale of the GPS network offers a unique opportunity to search for spatially-extended DM objects (or “clumps”), such as topological defects (TDs) <cit.>, which are otherwise not detectable by most  ongoing and planned DM searches.
The large number of clocks and the very large aperture of the network increase both the chance of an interaction and the sensitivity of the search, since we seek the correlated propagation of new physics signals throughout the entire network.
The large network diameter also increases the overall interaction time.
Therefore, by analyzing the  GPS timing data, one can perform a sensitive search for transient signals of exotic physics, and if no sought signals are found, stringent limits on the relevant interaction strengths can be placed.


Recently, our GPS.DM collaboration carried out an initial analysis <cit.> of the archival GPS data, looking for signatures of a particular type of TDs (domain walls, quasi-2D cosmic structures). 
While no such signatures were found, we placed limits on certain DM couplings to atoms that are many orders of magnitude more stringent than the previous constraints.
Here, we present a search method based on Bayesian statistics. We demonstrate that compared to our initial search, the Bayesian approach greatly increases the search sensitivity. This approach also broadens the discovery reach to more general DM models and to lower DM field masses.

Our approach is not limited in scope to the GPS network, but applies equally to other networks of precision measurement devices.
In principle, timing data from any other atomic clocks as well as data from other precision measurement devices can be included in the analysis.
In particular, there are similarities to another experiment, the Global Network of Optical Magnetometers for Exotic physics (GNOME), which employs a geographically-distributed Earth-based network of magnetometers to search for transient signatures of exotic physics, including topological defect DM <cit.>.
Techniques described in this work may prove useful for such experiments.





Beyond “clumpy” DM models, one can use  networks to search for other types of DM, such as non-self-interacting virialized ultralight fields (VULFs), that lead to signals that oscillate at the DM Compton frequency.
Such a search would rely on a multi-node spatio-temporal correlation function <cit.>.
One may also search for both transient and oscillating effects due to ultralight DM (and other exotic physics) with laser interferometers and gravitational wave detectors
<cit.>, 
by directly exploiting the scalar–photon coupling
<cit.>,
 atomic spectroscopy <cit.> and
noise statistics <cit.>,
 electric dipole moment searches
<cit.>,
and even pulsar timing <cit.>.








The structure of this paper is as follows.
Section <ref> reviews the background theory for topological defect DM, the DM-induced transient variations of fundamental constants, and how atomic clocks can be used to search for DM signatures.
Section <ref> discusses aspects of the GPS network relevant to our search.
In Section <ref> we describe the Bayesian statistics method for the data analysis and the search, and in Section <ref> we use this method with simulated data to demonstrate its efficacy.
Finally, in Section <ref> we present the projected sensitivity and the discovery reach of the search.

This paper has three appendices, which include the derivation of the  velocity distributions for macroscopic DM objects, a brief characterization of the noise properties of the clock data relevant to our search, and the expected signals for a few specific DM models.
The supplementary information <cit.>
 presents a detailed analysis of noise characteristics such as Allan variance, power spectrum, and autocorrelation  for individual GPS satellite clocks. 
Since the intended audience includes both atomic and particle physics communities, we restore ħ and c in the formulas. We use the rationalized Heaviside-Lorentz units for electromagnetism.









§ THEORY





 §.§ Ultralight dark matter and topological defects



Despite the extensive searches, both laboratory direct detection and high-energy collider experiments have so far failed to yield convincing evidence for the existence of weakly interacting massive particles (WIMPs) with masses ∼ 10 – 10^4GeV, see, e.g., 
Refs. <cit.>.
While WIMPs are theoretically well-motivated, they are by no means the only DM candidate.
The null WIMP searches have partially motivated searches for ultralight bosonic DM, such as axions <cit.>.
While direct DM searches with particle detectors rely on measuring energy deposition by individual DM particles, precision measurement techniques are well suited for detecting candidates that act as coherent entities on the scale of individual devices or their networks. 
In other words, precision measurement devices can be used for detecting ultralight DM and this approach probes the mass region that is
complementary to particle detectors.


Ultralight fields may form coherent (on certain time-scales) oscillating fields, or they may form stable macroscopic-scale objects 
<cit.>. The formation of macroscopic “clumpy” DM objects requires self-interactions in the dark sector.
An example of macroscopic DM are topological defects, which may have various dimensionalities: monopoles (0D), strings (1D), and domain walls (2D).  Depending on their cosmological fluid equation of state, these objects can contribute to both DM and DE. 


The interactions of light scalar fields with standard model (SM) fields can be phenomenologically parameterized as a sum of effective interaction Lagrangians (portals) <cit.>

    Ł_ int = Ł^ PS + Ł^ S^1 + Ł^ S^2  +  …  ,

where Ł^ PS represents the pseudoscalar (axionic) portal, and Ł^ S^1 and Ł^ S^2 are the linear and quadratic scalar portals, respectively.
The linear and quadratic scalar portals can lead to changes in the effective values of certain fundamental constants and thus cause shifts in atomic transition frequencies.  Atomic clocks in particular are sensitive probes of varying fundamental constants.
The axionic portal leads to interactions that mimic spin-dependent shifts due to fictitious magnetic fields, and  thus are well suited for magnetometry searches <cit.>.
We also note that there are stringent limits on the interaction strength for the linear scalar interaction coming from astrophysics and gravitational experiments (see, e.g., <cit.>).
However, the constraints on the quadratic portal are substantially weaker <cit.>.
For concreteness, here we will focus on the quadratic scalar portal.


While we  refer to specific models, namely topological defect DM with quadratic scalar couplings, it is important to note that the search technique is not limited in scope to this possibility.
Any large (on laboratory scales), “clumpy” object that interacts with standard model particles is detectable using this scheme.
Examples of such other models include Q-balls <cit.>, solitons <cit.>, axion stars <cit.>, and other stable objects formed due to self-interactions in the DM sector.








 §.§ Searching for dark matter with atomic clocks



Since the microscopic nature of DM is unknown, we take a phenomenological approach for the non-gravitational interactions with ordinary matter (see, e.g., <cit.>).
Explicitly for the quadratic scalar portal, we have 

    -Ł^ S^2=  ϕ^2(
    Γ_fm_f c^2 ψ̅_f ψ_f    
     +Γ_αF_μν^2 /4
    +  …),

where ϕ is the scalar DM field (measured in units of energy), m_f are the fermion masses,
ψ_f and F_μν
are the SM fermion fields and the electromagnetic Faraday tensor, respectively, 
and  Γ  are coupling constants that quantify the strength of the DM–SM interaction. 
There is an implicit sum over the SM fermions f in the above equation.
The above Lagrangian leads to the effective redefinition of  fundamental masses and coupling constants,

    α^ eff(ř,t)    = [1+ Γ_α ϕ^2(ř,t)]  α,
    
    m_f^ eff(ř,t)    = [1+ Γ_f ϕ^2(ř,t)] m_f,

where α≈1/137 is the electromagnetic fine-structure constant 
and m_f are the fermion (electron m_e and light quark m_q≡[m_u+m_d]/2) masses.
The coupling constants Γ have units of  [ Energy]^-2 and to aid the comparison with previous literature we also define the effective energy scales  Λ_X ≡ 1/√(Γ_X)  with X=α, m_e, m_q.

Considering TDs, the DM field ϕ^2→0 outside the defect, hence the effective couplings are only realized inside the defect.[Strictly speaking, this condition requires an auxiliary DM field, see Ref. <cit.> for the mechanism. ]
The  field amplitude inside the defect, A, can be linked to the average energy density inside the defect as ρ_ inside= A^2/( c   d^2), where d is the spatial extent or width of the defect. 
In TD models, the width d is set naturally by the field Compton wavelength, d=ħ/(m_ϕ c), where m_ϕ is the mass of the DM field particles; in general, we treat d as a free observational parameter.
Further, in the assumption that these objects saturate the local DM energy density, one can link A and d to the local DM energy density ρ_ DM,

    A^2 = (ħ c)  ρ_ DM v_g   d  ,

where 

is the average time between close encounters of the DM objects with a point-like device, and the galactic velocity v_g=⟨v|∼⟩300kms^-1∼ 10^-3 c is the average relative velocity of DM objects that cross paths with the Earth.


From Eqs. (<ref>) – (<ref>), we may relate the observable DM-induced atomic frequency shift to the transient variation of fundamental constants (and thus to the DM field parameters).
The fractional shift in the  frequency ω_0 of a particular clock transition can be expressed as

    δω(ř, t)/ω_0 
    = ∑_X _̨XΓ_Xϕ^2(ř, t)≡ ϕ^2(ř, t),

where X runs over relevant fundamental constants, and _̨X are dimensionless sensitivity coefficients.
For convenience, we introduced the effective constant, Γ_ eff≡∑_X κ_X Γ_X, which depends on the specific clock.

The dimensionless sensitivity coefficients _̨X are known from atomic and nuclear structure calculations.
For example, considering only the variation in the fine-structure constant α and ignoring relativistic effects, the optical and microwave transitions frequencies scale as
ω^ opt_c ∝α^2, and
ω^ mw_c ∝α^4, respectively.
Relativistic atomic-structure effects add small corrections to these scalings <cit.>.
For the microwave  Rb, Cs, and H clocks of the GPS network, the effective coupling constants read (using computations <cit.>)

    ( ^87Rb)   =4.34 Γ_-0.069 Γ_m_q+2 Γ_m_e,
    
    ( ^133Cs)   =4.83 Γ_-0.048 Γ_m_q+2 Γ_m_e,
    
    ( ^1H)   =4 Γ_α-0.150 Γ_m_q+2 Γ_m_e.

The values of _̨m_q come from a combination of shifts in the nuclear magnetic moment and  in the nuclear size <cit.>, and from the variation in the proton mass with δ m_p/m_p = 0.05 δ m_q/m_q <cit.>.

Although each clock type is sensitive to a combination of three coupling constants,
by combining results for three types of clocks within the network one can unfold individual coupling constants Γ_X or, equivalently, individual energy scales Λ_X. 
Until recently, the existing constraints came from observations of supernova emission <cit.>:
Λ_m_e,≳3TeV, and 
Λ_m_p≳15TeV.
More stringent constraints for certain regions of the (d,) parameter space have recently been placed on Λ_ using a laboratory optical Sr clock by the Toruń group <cit.>. 
Using 16 years of archival GPS data, our GPS.DM collaboration constrained Λ_, Λ_m _e, and Λ_m_q <cit.>; 
that initial search focused on domain walls. 
These newly-established constraints reach the ∼ 10^7TeV level depending on the size of the objects and the frequency of their encounters with the Earth.

With the model-specific theoretical background established, now we proceed to developing a method for a sensitive search for  macroscopic DM objects. 
We will demonstrate that  compared to our initial search <cit.>, the developed method improves the sensitivity by several orders of magnitude, and also substantially increases the range of probed DM field masses. 
It is also sufficiently general to enable mining for signatures of all the prototypical topological defects: monopoles, strings, and walls.
The method is Bayesian in nature and we start with describing known DM halo properties, velocity distribution and directionality, that serve as priors to the search.









 §.§ Priors on velocity distribution and event rate









We form our search priors based on the standard halo model (SHM), see, e.g., Ref. <cit.>.
Within the SHM framework, the velocity distribution of DM objects in the galactic rest frame is isotropic and quasi-Maxwellian; further details are given in Appendix <ref>.
The Milky Way rotates through the DM halo, with the Sun moving at v≃ 220kms^-1 in the direction towards the Cygnus constellation.
This defines the most probable incident direction for a collision with a DM object; in fact, more than 90% of events are expected to come from the forward-facing hemisphere, as shown below. 
We define the unit-vector, ň̂̌_g, that points from the Earth center along this direction.
Further, as shown in Fig. <ref>, we define the incident direction of the DM object, ň̂̌, to be pointing away from the Earth center, so that ň̂̌ = -v̌/v̌, where v̌ is the velocity vector of the DM object.
The  angle of incidence ψ is defined as ň̂̌·ň̂̌_g=-cosψ. According to this definition, the forward-facing angle ψ=π points in the direction of the galactic motion of the Solar system, towards the Cygnus constellation.



We consider three topological defect “templates”: domain walls, strings, and monopoles.
We assume that over the length scales of the GPS constellation a string/wall can be modeled to be straight/flat.
A domain wall that crosses the GPS constellation incident with a velocity v̌ that is at an angle η to the vector normal to the wall, would be indistinguishable from a wall (of the same width) incident with a (slower) velocity 

    v_⊥=vcosη

that is normal to the wall, see Fig. <ref>. 
We will refer to the v̌_⊥ component of the relative velocity v̌ as the “perpendicular” velocity, and define ň̂̌_⊥≡ -v̌_⊥/v̌_⊥.
The same argument applies to strings (for strings, v̌_⊥ is defined to lie in the plane containing v̌ and the string symmetry axis).
Therefore, in these cases, the more relevant quantity is the distribution of the perpendicular velocities, f_v_⊥; this distribution is derived in Appendix <ref>.


Since we are focusing on macroscopic DM objects, it is also instructive to consider the distribution of transit durations. 
The transit duration, =d/v, is defined as the time it takes a DM object of width d to sweep through a point in space (or a single device). 
Similarly we can consider _ GPS=D_ GPS/v, the time for the center of the DM object to pass the entire GPS constellation.
Our derived speed, incident angle, and transit duration distributions are shown for monopole- and wall-like objects in Fig. <ref>.


We treat the expected event rate, 1/ [see Eq. (<ref>)], as a free parameter.
This parameter can be linked to the number density of DM objects in the galaxy. 
For monopole-like objects (including non-topological solitons, Q-balls, bubbles etc.), the relevant quantity is the volume number density, while for strings and domain walls it is the areal and linear number densities, respectively.
Thereby,  can be related to the energy density inside the DM object as

    = ρ_ inside/ρ_ TDMd/v_g,

where ρ_ TDM is the galactic energy density of the considered DM objects.
In the assumption that these objects saturate the local DM density,  we have ρ_ TDM = ρ_ DM. 
Direct measurements of the local DM density give 0.3±0.1GeVcm^-3 <cit.>; we take ρ_ DM≈ 0.4GeVcm^-3 for definitiveness and to be consistent with recent literature.
Note that Eq. (<ref>) is model independent and applies to any DM object of characteristic size d.






§ GPS ARCHITECTURE AND CLOCK SOLUTIONS




A  detailed description of modern GPS data acquisition and processing techniques and their application in precision geodesy can be found in Ref. <cit.>. Details relevant to DM searches with GPS constellation
are given in Ref. <cit.>. Here, we briefly review the main aspects of GPS and introduce  relevant concepts and terminology.


GPS works by broadcasting microwave signals from nominally 32 satellites in medium-Earth orbit (altitude ∼ 20,000 km).
The transmissions are driven by an atomic clock (either based on Rb or Cs atoms) on board each satellite. 
It is namely the carrier phase of these microwave signals that is measured by the specialized GPS receivers and is used in deriving the GPS clock solutions.
Typically, each satellite houses four atomic clocks, only one of which is broadcasting at any given time. 
Clock swaps are marked in databases supporting the archival GPS dataset.
There are also a large number of ground-based receiver stations, several of which employ highly-stable H-maser clocks.
The more recent satellites predominantly employ Rb clocks as it has become clear that unpredictable variations in clock phases for the Cs-clock satellites are significantly worse than for Rb. 
As of early September 2017, there were 30 Rb satellites and only two Cs satellites in orbit.

The GPS satellites are grouped into several generations, called blocks: II, IIA, IIR, and IIF <cit.>, and each satellite is assigned a unique identifier known as the Space Vehicle Number (SVN).
Each subsequent block was built with significant improvements, and the effect of these improvements can be seen in the noise characteristics of the satellite clocks, as discussed below.
Block III satellites are currently under development, and are  to be launched from mid-2018. 
Table <ref> presents a summary, including the number of days worth of data that is available for each satellite block in the archival data set.
Further, the network can be extended to incorporate the network of Earth-based receiver clocks, as well as clocks from other Global Navigation Satellite Systems, such as the European Galileo, Russian GLONASS, and Chinese BeiDou, and networks of laboratory clocks <cit.>.





Here we analyze data generated by the Jet Propulsion Laboratory (JPL) <cit.>, in which clock time-series are given at τ_0=30  s intervals (epochs).
The data are  clock biases, that is the difference in the time readings (clock phases) between the given clock and a reference clock.
The same reference clock is used for the entire GPS network for each day.
The biases are generated using data from a global network of ∼100 GPS station receivers <cit.> by a mature analysis system that is used routinely for purposes of centimeter-level satellite orbit determination, and millimeter-level positioning for scientific purposes, such as plate tectonics, Earth rotation, and geodynamics. 

We also note that while the currently available clock time-series are sampled every 30 seconds, the raw GPS data is sampled every second for some stations. 
It is therefore possible to re-process the GPS data to generate higher-rate 1s clock solutions. 
This work is currently underway in our group. 
Notice that a fiducial DM object sweep through the entire constellation takes about 170 seconds, thereby it lasts for just 6 epochs for the currently available 30s sampling intervals.  Clearly, the resolution would improve for the 1s data.  

In the initial GPS data processing (performed by JPL <cit.>), there is effectively no restriction on the allowed behavior of the clocks from one epoch to the next.
Crucially, if a clock were to have a real transient that far exceeded engineering expectations, the data over that time window would not have been removed as outliers. 





The clock biases from JPL <cit.> also come with a “formal error”, σ_F.
The formal error, typically on the order of σ_F∼0.02-0.03  ns, quantifies uncertainty in the determination of the clock bias, and does not directly incorporate the intrinsic clock noise or slowly varying biases due to correlated orbit errors and higher-order general relativistic effects (∼0.1ns).
Only the most recent satellite clocks (Rb clocks on board the block IIF satellites) have observed temporal variations from one epoch to the next that are at a similar level as the formal error, indicating that temporal variations in older clocks are actually due to clock behavior rather than estimation error.
In fact, the observed variances in the data from the most modern Rb IIF satellite clocks are significantly better than suggested by the formal error, see Appendix <ref>.


Due to frequency drifts and other long time-scale (∼ hours) effects, it is typical for a second-order polynomial (y_2) to be subtracted from the raw GPS time series data before the analysis <cit.>.
One may form y_2(j) for each clock, for each individual day, using a weighted least-squares approach, taking the weights as the inverse of the formal error.
Then the polynomial-reduced data (residuals) are defined
d_j^(0) = x_j-y_2(j),
where {x_j} are the raw time-series data, and j denotes the same-time (epoch).
This procedure is useful for visualizing the data, however,  it is not necessary for our analysis. 
Unless noted otherwise, we take  d_j^(0)≡ x_j in this paper.




The relative phase of an atomic clock (bias) d_j^(0) is a non-stationary time series, dominated by random walk noise.
To perform the analysis, we must first “whiten” the data.
To this end, we employ (depending on the clock type, as discussed below) either a first- or second-order differencing, and define

    d_j^(1)   ≡ d_j^(0)-d_j-1^(0), 
    
    d_j^(2)   ≡ d_j^(0)-2d_j-1^(0)+d_j-2^(0).

In general, first-order differencing is sufficient for Rb clocks, while second-order differencing is required for Cs clocks. 
Since d_j^(1) is proportional  to the discreet derivative of the clock biases, we refer to it as a pseudo-frequency.
Further discussion of the clock noise characteristics is presented in Appendix <ref>; see also the Supplementary Information, where we quantify the noise characteristics of each satellite clock individually. 








§ BAYESIAN SEARCH FOR DM EVENTS




 §.§ Likelihoods and odds ratio



A DM-induced perturbation in the  device data would be indistinguishable from a perturbation caused by other external non-DM factors or random statistical processes.
The key then, is to rely on the correlated propagation of clock “glitches” across a network caused by the sweep of a DM object through the network. 
Based on the standard halo model (see Sec. <ref>), DM objects are expected to travel relative to the Earth with galactic-scale speeds, v_g∼300kms^-1, incident from a certain direction.
Thus the speed and the directionality of the sweeps serve as DM signatures.  
There are similarities between the method we describe and those employed in gravitational wave detection by the LIGO collaboration, see, e.g., Refs. <cit.>.


Consider a candidate model, denoted M, that predicts a DM signal across the network (for example, the passing of a  domain wall). 
In order to determine whether such a model is supported by the data D, we employ a Bayesian technique, see, e.g., Ref. <cit.>.
In Bayesian statistics, model selection is based on forming the odds ratio of two probabilities (likelihoods)

    O_M,M̅(j_0)=p(D_j_0|M,I)/p(D_j_0|M̅,I).

Here, M̅ denotes the proposition that no signal is present in the data, i.e. the data is purely random,
 and I encodes the knowledge of the SHM priors discussed in Sec. <ref>. 
The data stream D is sampled at intervals of _0 (_0=30s for our current GPS data set).
Since we search for transient events of finite duration, the odds ratio is tested in a time window of length J_W points, centered at epoch j_0. 
The value of J_W is determined by the maximum duration of the transient signals to be tested.
In the analysis, we scan over j_0 for a fixed value of J_W, so the odds ratio is an explicit function of j_0. 

The likelihoods p(D_j_0|M,I) entering Eq. (<ref>) are described by the Gaussian multi-variate distributions marginalized over model parameters,

	
    p(D_j_0|M,I)   =
    		K∫^̣3v   p(v̌|M,I)
    		∫ḥ  p(h|M) ∫x̣  p(x|M) ∫_ (j_0-1)τ_0^ j_0τ_01/_0 ṭ_0 exp(
    			-1/2χ^2(s)
    		) .


Here,
v̌ is the velocity of the incident DM object in the Earth-centered inertial (ECI) frame,
h quantifies the amplitude of the DM signal,
x stands for the remaining model-specific parameters of the DM object,
and K is a normalization factor.  
Further, t_0 is the moment of time at which the DM object passes by the center of the Earth. 
It is assumed to occur in the time interval ( (j_0-1)_0, j_0 _0 ], and we marginalize over t_0 in the last integral (<ref>).
Note that compared to our initial search <cit.>, the single-device sweep time (d/v) may last longer than _0.



The data and the model-prescribed DM signal are combined in the argument of the exponential, 

    χ^2(s)
    	=∑_ab^∑_jl^J_W[d^a_j-s^a_j](E^-1)^ab_jl[d^b_l-s^b_l],

where E is the covariance matrix discussed in the following section. 
Here and below we use the “upstairs” indices to label devices, and the “downstairs” indices to denote epochs (sampling times).
The indices a and b run over all  devices in the network, and the indices j and l run over the J_W data points in the  time window.   The device data d notation is generic and it can stand for the singly– or doubly–differenced clock bias data, Eq. (<ref>).
Finally, s^a_j=s^a_j(M,t_0,v̌,h,x) is the model-prescribed DM signal in device a at epoch j, discussed in Sec. <ref>.


Continuing with the discussion of factors entering the likelihood, p(v̌|M,I) is the (normalized) probability density for the velocity distribution of DM objects in the ECI frame.
In the case of monopoles, for example, it is reasonable to take this to be given by the SHM.
Likewise, the function p(h|M) is the normalized probability density for the DM signal amplitude in the time series, and is described by a flat prior.[We note that the normalization for the h prior is arbitrary. For our purposes it is not important, since we do not rely on the actual value of the likelihoods function but rather define some threshold, above which false-positives are sufficiently rare, as discussed in the following sections.]
To calculate the likelihoods, we perform the integral over h analytically (possible because s is linear in h, see below), and use a randomized Monte-Carlo integration for the other parameters.


Finally, the likelihood that no signal is present in the data
is given simply by

    p(D_j_0|M̅,I)=K
    		exp(
    			-1/2χ^2(0)
    		),

where the DM signal is set to zero. 
The window size J_W is kept the same as in the  p(D_j_0|M,I) computations.


The likelihood functions (<ref>) and (<ref>) are calculated for every available epoch j_0, and the odds ratios (<ref>) are formed. Large spikes in the odds ratio as  a function of j_0 can indicate potential DM events.






 §.§ Correlations and covariance




The covariance matrix entering Eq. (<ref>) is defined as 

    E^ab_jl≡⟨d^a_jd^b_l|,⟩

where ⟨⋯|$⟩ denotes averaging.
To compute its elements, one requires a stationary time series, for which (depending on the clock type) we use either the first- or second-order differenced data (<ref>).
Note that for pure uncorrelated white noise, the covariance matrix is completely diagonal, with elements given by the variances.
In this case, the matrix inversion required for computingχ^2(s)(<ref>)  is trivial.
Realistic device noise is, however, correlated. 
Specific to the GPS clocks, additional short-range anti-correlation for individual clocks is introduced by the propagation of the formal error (which is roughly white noise ind^(0)) through the differencing procedure (<ref>).
Moreover, the underlying clock biasesd^(0)are the differences between the phases of the given clock and a reference clock. 
Since the reference clock is common to all clocks, biases and the differenced datad^(1)andd^(2)are correlated between different clocks. 

It is convenient to split the covariance matrix into two contributions,E=A+B, where 

    A^ab_jl   ≡ E^ab_jlδ^ab,
    
    B^ab_jl   ≡ E^ab_jl(1-δ^ab).

The first term,A, represents the correlation between data points for a single clock, i.e.,  auto-correlation.
TheBcontribution describes the correlations between different clocks, and is referred to herein as the cross-correlation.
The autocorrelation part of the covariance matrix is block diagonal, built fromindependent symmetricJ_W ×J_Wmatrices.
The elements ofAdepend only on the distance from the diagonal, and can be related to the autocorrelation functionA^a(Δt)asA^aa_jl = (^a)^2 A^a(Δt_jl),whereΔt_jl = j-lτ_0is the lag and^ais the standard deviation (see Appendix <ref>).


Each clock in the network is referenced against a common reference clock.
This adds a common noise component to all the data streams, and is main source of cross-correlations.
Therefore, the time series for each clock can be decomposed as

    d^a_j = e^a_j + c_j,

wherec_jis the component due to the shared reference clock, ande^a_jis the component unique to clocka(⟨e^a_je^b_l|=⟩0fora≠b).
Then, it is clear thatBdepends only on the reference clock, and is independent ofa,b:

    B^a,b≠ a_jl = ⟨d^a_jd^b_l|=⟩⟨c_jc_l|≡⟩b_jl.






To calculate the likelihoods, we need to invert the covariance matrix.
First, we note that the Earth-based H-maser clocks used as reference in the JPL data processing are typically much quieter than the  satellite clocks.
 Therefore, the cross-correlation contributionBis typically smaller thanA, soBcan be treated perturbatively. 
Further, we may neglect the even smaller termsB^ab_jlwithj≠l, definingb_0≡B^ab_jj.

Thus, we express the inverse of the covariance matrix asE^-1 = H + W,where

    H^aa_jl   = (A^-1)^aa_jl ,  
    
    W^ab_jj   ≈-b_0/(^a^b)^2(1-δ^ab).

Then, Eq. (<ref>) can be expressed (withη≡d-s) as

    χ^2(s)=
    ∑_a^∑_jl^J_Wη^a_j H_jl^aa η^a_l
    -
     ∑_a ≠ b^∑_j^J_Wb_0 η^a_j η^b_j /(^a^b)^2.


The described approximation holds when the clock noises far exceed that of the reference clock. 
This approximation breaks down if the network includes clocks with noise levels similar to that of  the reference clock.
For example, when including multiple station, Rb-IIF, or laboratory clocks,  Eq. (<ref>) is no longer valid.
In this case, we define the weighted mean of all (other) clocks

    d̅^a̅_j = ∑_b≠ a d^b_j   (^b)^-2/∑_b≠ a(^b)^-2≈ c_j ±σ/√() ,

which is subtracted from each time series [Eq. (<ref>)]:

    d^a_j-d̅^a̅_j ≈ e^a_j ±σ/√().

Here,is the typical standard deviation of the clock data.
Each new data stream still contains a common component,∼σ/√(), however this is small enough so that the above approximation (<ref>) holds true.
In these cases, the same procedure must be applied also to the expected signalss^a_j →s^a_j-s̅^a̅_j .











 §.§ Transient dark matter signals


The likelihood function in Eq. (<ref>) requires a model-prescribed DM signal,s^a_j, for the data streams to be compared against.
The DM signal depends on the assumed coupling strength to the device, and on the kinematics and spatial structure of the DM object (monopole, domain wall, etc.). To quantify the transient signal we need to specify the collision geometry. 
We work in the ECI (Earth-centered inertial) J2000  frame, which has its origin (denoted ECI0) at the center of mass of the Earth, andz-axis aligned with Earth's spin axis. 
Thex-axis is aligned with the mean equinox  at 12:00 Terrestrial Time on 1 January 2000.
The important aspect is that the ECI frame orientation remains fixed in the galactic rest frame, i.e., it does not rotate with the Earth.

Here, we consider three generic and geometrically unique templates: walls, strings, and monopoles. 
Albeit more complex geometries are plausible, such as walls closing on themselves forming cosmic bubbles, the presentation below is sufficient for extending the formalism to such more complex object geometries. 
While the field profile inside the DM object can be arbitrary, we focus on Gaussian profiles.  
Beyond qualitative arguments, the reasons for Gaussian-profiled objects can be also supported by Bayesian logic. 
Indeed, application of the maximum entropy principle to a distribution with the  mean and variance (determined by the defect sizedin our case) as the only given information yields the Gaussian distribution <cit.>.
In any case, the presented formalism can be applied to arbitrarily-shaped DM object profiles.

We  assume that the linear trajectory and velocity of the DM object are not affected by the gravitational pull of the Earth or the portal couplings to the Earth constituents, and that the shape of the DM object is preserved through the encounter. 
Another assumption is that the DM encounters are well separated, i.e. DM objects do not overlap and at most one of them interacts with the entire network at any given time. 
Finally, we assume that objects lacking spherical symmetry do not rotate.







Consider an event in which the center of a DM object moving with velocity v̌ crosses the plane that is perpendicular tov̌and contains ECI0 at timet_0, as shown in Fig. <ref>.
The accumulated time bias between a clock thats frequency is perturbed byδωand an unaffected clock (ω_0) is given by∫_-∞^tδω(t')/ω_0ṭ'.
Therefore, at timet, the DM-induced clock phase bias in clockareads

    s^a^(0)(t) =  ∫_-∞^t 
    		[ 
    			h^a φ_M^2(t^a,ρ^a,t') - h^R φ_M^2(t^R,ρ^R,t')
    		] ṭ' 	,

whereφ_M^2is the normalized profile[The DM “profile”φ differs from the field ϕ [Eq. (<ref>)] only by normalization, and is defined for convenience; see Appendix <ref>.] of the DM object (for specific modelM),ρ^ais the impact parameter,h ∝A^2 Γ_effis a clock-specific constant that determines the magnitude of the signal in the data (see Appendix <ref>), andt^a(t^R) is the time of encounter for clocka(reference clock).
The time of encounter is defined as the moment the DM object passes by clocka.
More precisely, it is the time at which the center of the DM object (central plane for walls, or central axis for strings) crosses the plane that is perpendicular to v̌ and contains the given clock:

    t^a = t_0-ř^a·ň̂̌/v,

whereň̂̌is the unit vector that points from ECI0 parallel to the incident direction of the DM object
(v̌=-vň̂̌, see Fig. <ref>), andř^ais coordinate of clockain the ECI frame. 
The satellite and ground station positionsř^aare a part of the JPL GPS dataset, and are known with∼cmand∼mmaccuracies, respectively.
Note that the impact parameters are zero for domain walls, but are, in general, non-zero for strings and monopoles; see Appendix <ref>.








For domain walls and strings, we usev̌_⊥andň̂̌_⊥, see Fig. <ref> and the discussion around Eq. (<ref>). 
The discreet matrixs^a_jis generated by integrating to the specific values oftthat correspond to the GPS epochs (data sampling times).
Then we form either the first- or second-order differenced DM signals as in Eq. (<ref>), withd→s.


The particular form ofφ^2_Mdepends  on the spatial structure and the kinematics  of the DM object.
In Appendix <ref>, we present explicit signals for domain walls, monopoles, and strings, and linkφ_Mandh^ato the field parameters for these templates.

In our discussion of DM signals, we neglected the Earth orbital motion about the Sun at∼30 km/s, orbital velocities of satellites about the Earth (∼4 km/s), and the ground station rotational velocities (∼0.5 km/s). 
While these velocities are much smaller than the galactic velocities, the motional effects can become important for large-scale or slowly-moving objects. 
For example, the motional effects become relevant if the overall duration of an encounter is comparable to the 12-hour satellite orbit.
The modification of the DM signal templates to account for clock motion is straightforward, as the satellite and ground station positions are known. 
We leave this generalization for future work.






 §.§ Mixed networks



There are several different clock types (Cs, Rb, H-maser) in the GPS network.
As our search is expanded to include other laboratory clocks (and other high-precisions sensors) the diversity will increase further.
Each clock species may respond differently to the interaction with the DM field, see Eq. (<ref>).
Therefore, we cannot assumehto be uniform  across the network.

There are several approaches for inhomogeneous networks.
One approach, as per Ref. <cit.>, is to consider separately the homogeneous sub-networks (e.g., consider only the Rb clocks).
The major drawbacks of this approach is that we lose the benefit of the highly-stable H-maser reference clocks (none of the GPS satellite clocks have H-masers), and that we also limit the total number of clocks that are considered at any given time.

The simplest approach is to assume that one of the couplings in Eq. (<ref>) dominates, and carry out the analysis separately for each case.
For example, we may assume that|Γ_α|≫|Γ_m_e|, |Γ_m_q|in Eq. (<ref>).
The drawback of this approach is that it does not account for the possibility that several couplings may produce effects that are of a similar magnitude.


Furthermore, a Bayesian-like approach is to introduce additional marginalization parameters for each extra free parameter in place ofh.
The number of such free parameters is equal to the smaller of either the number of distinct clock species in the network, or the number of distinct couplings we consider.
For example, considering a network of Rb, Cs, and H clocks (as per GPS), we can substitute∫ḥ→∫ḥ_Rb∫ḥ_Cs∫ḥ_Hin Eq. (<ref>).









 §.§ Directional signatures




A possible scenario is that  a large number of small events are flagged by the Bayesian search (by “small” we mean the magnitude of the signal in the data compared to the clock noise, or the small magnitude of the spikes in the odds ratio).
Of course, such events may be  simply due to random statistical fluctuations, or other conventional-physics non-DM perturbations.
Here we consider signatures unique to DM (or other galactic sources) allowing us to exclude non-DM signals. 
While these signatures are  included in the Bayesian approach through the priors (e.g., the likelihood are suppressed for velocities outside the SHM range through the prescribed velocity distribution prior), we could also examine inferred values of collision geometry parameters through the Bayesian parameter estimation, as discussed in Sec. <ref>. 
Being able to resolve the event velocity magnitude and directionality is a powerful feature of geographically distributed networks.


First of all, the distributed nature of the network offers the direct sensitivity to the magnitude of DM object velocities. 
If the observed incident velocity falls too far outside of the bounds allowed by the standard halo model, then a DM origin can be excluded. 
There is also sensitivity to the directionality of the DM object velocity. 
The most probable incident direction is from the average forward direction of the Sun's motion through the galaxy (roughly from the direction of the Cygnus constellation), see Fig. <ref>.
We are only aware of one external systematic effect that has propagation speeds comparable tov_g, which is the solar wind <cit.>.
This effect, however, can be vetoed out on the basis of distinct directionality from the Sun, and by the fact that the solar wind does not affect the satellites in the Earth's shadow.



In addition to individual event signatures, one can also focus on the overall event statistics, provided the event rates are sufficiently high on the yearly basis <cit.>. 
For example, due to the∼10%annual variation in the relative velocities of the Earth and Sun in the galaxy <cit.>, one would expect to observe the annual modulation in the event rate.
This approach parallels the method employed in WIMP searches, e.g., Refs. <cit.>. 
Unlike WIMP searches, where the event rate may depend strongly on the DM velocity <cit.> (due to energy dependence of the cross section), here the rate is linear inv.
Also unlike (most) WIMP searches, the distributed network approach is additionally sensitive to the annual modulation in the average incident velocity  direction, which varies by∼20, as shown in Fig. <ref>.
(A WIMP-detection scheme that does have directional sensitivity is presented in  Ref. <cit.>.)








§ BENCHMARKING THE METHOD



 §.§ Simulating realistic clock time series





We generate simulated time series data that have the same noise characteristics as the real clock data for each individual GPS satellite.
This is achieved by “coloring” pseudo-random white noise with the known power spectral densities for each clock <cit.>.
We calculate the power spectral densities for each specific SVN using the clock data provided by JPL, as the clock performance may degrade over time, and the clocks can perform differently when in orbit than when tested in a laboratory environment. We also simulate cross-correlations (correlations between different clocks).
This is achieved by simulating a reference clock, which adds a common noise stream to all the clocks in the network. 


In Fig. <ref>, we plot several arbitrarily selected real JPL clock solutions,d^(0), alongside the simulated clock solutions for the corresponding SVNs (denotedz^(0)) to demonstrate the quality of simulated data.
The standard deviations of the simulated data (after first- or second-order differencing) match exactly those of the real data for the given SVNs.  
Further, the longer-scale noise characteristics also match  – in Fig. <ref>, we plot the autocorrelation functions and Allan variances for both the simulated and real data for the same clocks. 
These figures demonstrate that the simulated clock data do indeed have the same noise characteristics as the real data.

Having generated simulated time series, we can test our Bayesian search code in a number of distinct ways:

  * 
To gauge the prevalence of statistical false-positives, we run the code for the event-free simulated data.

  * 
We inject DM event signals into the simulated data streams to gauge the efficacy of our technique to pick out true-positive events.

  * 
We inject “bad” events (i.e., signals that are not properly correlated) into the simulated data as a test of the robustness of the search technique.

  * 
We use parameter estimation to extract the observed parameters of the injected DM event, and compare the results to those used to generate the injected DM signal as a test of the method accuracy and efficacy.



 §.§ Prevalence of statistical false-positives




We wish to define a threshold for the odds ratio,O_M,M̅(j_0), Eq.(<ref>). 
If the spike in the odds ratio is larger than this threshold, such an event can be investigated as a potential DM event.
In order to do this, we need to calibrate the rate of statistical false positives. 
To this end, we ran multiple simulations of various GPS clock network configurations, and computed the odds ratio (<ref>) for each epoch.
For each combination of clocks, we considered 2048 realizations of 2048 30s-epochs, amounting to approximately 4 years of simulated data for each simulation.
A plot of the rate of false positives  as a function of the threshold is presented in Fig. <ref>.
This plot is for the specific DM model of “thin” (d ≪10^4km, see Appendix <ref>) domain walls. 


Note that for this exercise, a false positive is counted whenever an epoch has a value for the odds ratio above the given threshold.
This is a conservative definition, since the “width” of the odds-ratio spike (due to the imperfect resolution) may lead to the same false-positive event appearing in more than one neighboring epoch. 
By our definition, this will be counted several times.


We can also drastically reduce the number of false positives that occur by introducing a minimum value (magnitude),h_min, for the integral over signal magnitudes (<ref>).
Of course, this also means we can only detect positive events with|h|>h_min.
We can then perform the analysis in several sweeps, systematically reducingh_mineach time until signals of a given magnitude can no longer be excluded.












 §.§ Detecting injected DM events



To determine the sensitivity of the method, we must know the probability of positively detecting DM events of a given magnitude.
To this end, we generate clock data per Sec. <ref>, inject randomized DM signals into the data streams and compute the odds ratios. In Fig. <ref>, we present one such simulation as an example. 
Here we show the first 8 (of 30) simulated time series' for a 1.5 hr window.
In this example, for simplicity, the clock noise was taken to be white (ind^(1)).
Then we injected a single “thick” domain-wall event for a wall of sized=10^4  km; the velocity and incident direction were chosen randomly.
The odds ratio was calculated for each epoch.
The spike in the odds ratio at the event is apparent, while the event is not discernible by eye in the data streams.
Note that the search routine is isolated from the simulation – it is not made aware of the event time, speed, direction, magnitude, wall width (or if there was an event at all).





Figure <ref> shows the fraction of injected thin domain wall events that are correctly identified, as a function of the signal magnitude.
The velocity and incident direction for each wall was chosen randomly (according to the SHM distributions, Fig. <ref>), and we assumed all clocks were affected by the DM in the same way (i.e., all clocks have the sameΓ_eff).
For this analysis, the odds threshold was set to allow fewer than 10 false positive events per year (O_thresh ∼10^3, see Fig. <ref>).
Note, for 30 s data, there are over10^6epochs in a year.
We count an event as “found” if there was a spike in the odds ratio above the determined threshold that appears within±1epoch from the injected incident timet_0.
(Of course, occurrences where a single event leads to an odds-ratio spike for more than one epoch are not double-counted, only one event is injected per trial, and it is either found or not.)
Also shown in Fig. <ref> is the average of the log odds ratio for each of these simulated networks as a function of the magnitude of the injected domain wall signal.
The large “gap” in the sensitivity that occurs around 2010 is due to the introduction of the Rb-IIF satellite clocks, which are substantially more stable than the older generation satellite clocks; see Appendix <ref>.


In Fig. <ref>, we show the same true- and false-positive test results, but for networks of a varying number of identical pure white frequency noise devices (d^(1)equivalent).
This is to demonstrate the general efficacy of the method, without specific reference to the properties of the GPS data.











Note that the ground receiver clocks contribute only minimally, even though they are significantly more precise than the GPS satellite clocks.
That is because our current data is sampled only every 30 s, which is about the time it would take for a DM object to cross the Earth, meaning many of the Earth-bound clocks will be affected by the DM during the same data acquisition interval.
As discussed in Sec. <ref>, it is possible to re-process the existing raw GPS data to produce 1 s sampled time series.
In addition to the statistical improvement from the larger data set, this would also further allow us to take full advantage of the highly-stable Earth-based receiver and laboratory clocks.
Of course, this advantage comes at the cost of significantly increased computation time, which scales (roughly) quadratically with the number of data pointsJ_Wdue to the correlations, see Eq. (<ref>).












We also check the “robustness” of the method, to ensure incorrectly correlated events (that may exist in the data due to Earth-sourced or other non-galactic perturbations) are not flagged as potential events.
To do this, we inject a single perturbation of a specified magnitude into each satellite data stream at a random epoch, all within the same 5 minute time window.
This simulates a domain wall crossing, except in the important fact that the network perturbations are not correctly correlated between different satellites.
Injecting a large2σperturbation of this kind (σis the typical standard deviation of the clock noise) into the simulated data streams for a 30-clock network, fewer than 1% present odds ratios anywhere within the 5 minute window that are above the threshold.







 §.§ Parameter estimation



When a spike in the odds ratio is  above the pre-determined threshold value, we can investigate this region of data as a potential event.
For example, by finding the set of “best-fit” parameters that maximize the un-marginalized likelihood, 
we can estimate the properties of the possible event (e.g., the time of arrival, size of the object, coupling strength etc.). 
In Fig. <ref>, we show histograms of the parameter estimation for a number of simulated trials where event signals were randomly injected into simulated data.
Shown in the plots is the difference between the injected value and the extracted best-fit value for the crossing timet_0, speedv, and incident polar angleθ, for simulated domain wall crossings.
These parameters are representative of the spatial and temporal resolution of the method.
Note that Fig. <ref> was generated for 
 30 s sampled data – using the re-processed 1 s data (as discussed above) will lead to a substantially improved resolution in the arrival time and velocity.

We also perform the parameter estimation for the false positive trials, where the analysis is performed on simulated event-free data.
The resultant histograms are presented in Fig. <ref>.
In this case, when neglecting the priors, the histograms are flat, with a  slight bias of more false-positives towards higher velocities.
When including the priors, the distribution of extracted parameters from the false-positives match the priors, as expected.

This means that there is a potential to search for events even below the “false positive floor”.
Reducing the odds ratio threshold will allow us to detect much smaller DM events, but will also lead to a larger number of false positives.
The true positive results, however, are expected to follow the distribution of velocities and incident directions predicted by the standard halo model.
This is relevant for the part of the parameter space with≪1year.
There would also be expected annual modulations in the event rate, average event speed, and most-common incident direction, see Sec. <ref>.
In this case, the analysis would need to be performed without the priors (i.e., assuming flat priors) to avoid biasing the false positives.
















§ SEARCH SENSITIVITY AND DISCOVERY REACH



Combining Eqs. (<ref>), (<ref>), and (<ref>), we find the maximum signal amplitude observable in a given clock (a, with reference clockR) for a  domain wall crossing to be

    s^(1)_ max≃  (ħ c)  ρ_ DM√(π) d τ̃v_g 𝒯[   Γ_ eff^a  -  Γ_ eff^Rexp(-L^2/d^2)  ],

where 
the interaction duration is givenτ̃= d/vford/v<τ_0andτ̃= τ_0otherwise (τ_0=30 s is the time period between data sample points for GPS),
andL∼10^4 kmis the  distance between the clock  and the reference clock.




The subtraction of two terms in square brackets in Eq. (<ref>) is due to the fact that when the maximum of DM field affects the clock, the reference clock is affected by its exponentially-suppressed tail.
When employing a network of identical clocks, this term leads to a fast decline in sensitivity for larged. 
This is because the clock and reference clock are affected in the same way, so no bias is built up between them.
In contrast, when employing clocks with significantly different effective couplingsΓ_eff(particularly, combining microwave and optical clocks) this suppression is not realized.




Statistically, the minimum detectable signal is proportional to

    s^(1)_ min∝_y(_0)_0/√(N_ clk N_ pts),

whereN_clkis the number of affected clocks, andN_pts ∼d/v_g_0is the number of data samples per clock for which the DM-induced signal is appreciable, and_y(_0)is the Allan deviation.
The proportionality constant depends on the efficacy of the search technique, and onO_thresh, the odds ratio threshold required to eliminate false positives.
Therefore, we should have sensitivity to the region

    Γ_X ≳ε σ_y(τ_0)τ_0   (_̨X^a-_̨X^Re^-L^2/d^2)^-1/ c ρ_ DM √(N_ clk N_ pts) v_g   d  ,

where the factorε∼O(1)is the efficiency factor determined from the simulations, and depends onO_thresh.
From the results presented in Fig. <ref>, for a 90% detection confidence level, and when requiring fewer than ten false-positives per year, we haveε≈5for the existing GPS data.
Future improvements in the search method should allow us to further decreaseε.



Theεfactor depends only fairly weakly on the search parameters.
For example, as shown in Figs. <ref> and Fig. <ref>, increasing the odds threshold by a factor of 10 decreases the number of false positives by a factor of 10, while only increasings_min, the smallest detectable signal magnitude, by∼10%.
Therefore, we may estimate that for a 90% detection confidence level, and when requiring fewer than one false-positive every 10 years,ε≈6.




The average time between consecutive encounters with a DM object,, is considered a free parameter in our model (set by the number density of the DM objects).
The dependence of Eq. (<ref>) oncomes via the DM field amplitude (<ref>), and the requirement to not oversaturate the galactic DM density; the higher the number density of objects, the lower the field amplitude must be per object to compensate.
In order to determine the maximumthat one can have sensitivity to, we assume the sequence of DM events can be modeled as a Poissonian process.
For example, if we expect one DM object to cross the Earth every period ofon average, then in order to be∼90%confident that an event would have occurred in the observation timeT_obs, we must requireT_obs ≳2.3 .




We present the projected sensitivity of our search in Fig. <ref>, along with the existing constraints.
To be consistent with existing literature, we present the sensitivity in terms of the effective energy scales,Λ_X ≡1/√(Γ_X).
Specifically, we show the projections forΛ_α; the projections forΛ_m_eandΛ_m_qare essentially the same, the only difference arising from the different sensitivity coefficientsκ_X, see Eq. (<ref>).


The reduction in sensitivity aboved≃10^4kmfor the homogeneous clock networks is due to the fact that large DM objects will interact with the clock and reference clock at the same time, see Eq. (<ref>).
That is, above this value, we are only sensitive to the gradient in the DM field when using a homogeneous network.
The limits from our previous work <cit.> have a sharp cut-off above this value, since in that work, we required that the DM signal would be present for just a single data point (see Appendix <ref>).
The Bayesian method presented in this work does not suffer this constraint.


Performing the simulations for strings and monopoles is substantially more computationally demanding, due to the number of extra free parameters that must be marginalized over (see Appendix <ref>).
However, the sensitivity can be approximated by analogy with the domain wall case.
Ford≫R_GPS, the monopole and domain wall cases are essentially the same.
Ford < R_GPS, the sensitivity of the search can be estimated by noting the typical number of clocks that would be affected in a monopole crossing,N_eff ≃⌈N_clk   d^2/R_GPS^2⌉.
A similar equation exists for strings,N_eff ∝N_clk   d/R_GPS.
For strings and monopoles, we required that at least 3 satellite clocks are affected during the DM sweep,N_eff≥3, which leads to a sharp drop in sensitivity for smalld, as shown in Fig. <ref>.














§ CONCLUSION

We have described a method to use data from a distributed global network of precision measurement devices to search for transient signals that may be associated with sweeps by macroscopic-scale dark matter.
In particular, we considered the network of microwave atomic clocks on board the GPS satellites and ground stations, for which nearly two decades of archival data is available.
The method was demonstrated using simulated atomic clock data, and the prospects and discovery reach for topological defect dark matter was presented.
This approach can be extended in a straightforward fashion to other networks of high-precision measurement devices.






This work was supported by the U.S. National Science Foundation grant PHY-1506424.
We thank Chris Pankow, Derek Jackson Kimball, and Tigran Kalaydzhyan for discussions.
BMR is grateful to the CIERA institute and Northwestern University for hospitality during the CIERA Data Analysis Workshop,
and to the Perimeter Institute for Theoretical physics for support to attend the New Directions in Dark Matter and Neutrino Physics workshop and acknowledges the many helpful discussions that took place there.








§ APPENDIX



§ VELOCITY DISTRIBUTION AND EVENT RATE



Assuming the standard halo model, the velocity distribution of DM objects in the galactic rest frame is isotropic and quasi-Maxwellian, with dispersion of290kms^-1and a threshold above the galactic escape velocity ofv_ge≃544kms^-1.
The vector velocity distribution for DM objects that cross paths with the Earth can be expressed in the Earth-centered inertial (ECI) frame as 

    f_v(v̌)= C v exp[- (v̌+v̌_g)^2/v_c^2]Θ(v_ e(ψ)-v),

whereΘis the Heaviside step function with

    v_ e(ψ)=√(v_ ge^2-v_g^2sin^2ψ)-v_gcosψ

being the effective escape velocity (the maximum allowable relative DM velocity as a function ofψ),ψis the angle betweenň̂̌_g=v̌_g/v_g(the direction of Earth's motion through the galaxy) andň̂̌=-v̌/v(the vector of the incident DM object,cosψ=-ň̂̌·ň̂̌_g),v_cis the speed of the Sun in the galactic rest frame,v_gis the galactic speed of the local reference frame (ECI),
andCis a normalization constant.
For the purposes of this work, we can neglect the smaller relative velocity of the Earth in its orbit around the Sun, and takev_g≈v_c≈220 km/s.
For a more detailed overview, see, e.g., Ref. <cit.>.



The direction of motion of the solar system through the galaxy points towards the Cygnus constellation; in the ECI frameň̂̌_g≈(0.46,-0.49,0.74)^T.
The angular distribution function for events is obtained by integrating over  velocities

    f_ψ(ψ) = 2π∫_0^∞ f_v(v,ψ) v^2sinψ ṿ .

In fact, more than 90% of events would be expected to come from the forward-facing hemisphere, see Fig. <ref>.
Similarly, the scalar velocity distribution can be found by integrating over the angles

    f_v(v)=2π∫_0^π f_v(v,ψ)  v^2 sinψ ψ̣
    ≈ C v^2[ exp(-(v-v_c)^2/v_c^2) - exp(-(v+v_c)^2/v_c^2)].





Since we are focusing on macroscopic DM objects, it is also instructive to constructive to consider the distribution of crossing durations.
Define=d/vto be the time it takes a DM object of widthdto pass by a point in space (similarly we can consider_GPS=D_GPS/v, time to sweep the entire GPS constellation by the DM object central point).
It is convenient to define the inverse velocityu≡1/v, and its corresponding distribution function

    f_u(u)≡P̣_u(u)/ụ=ṿ/ụP̣_v(v)/ṿ=v^2f_v(v),

where P̣_u(u)is the infinitesimal probability for the DM object to have inverse velocityu.
Then, the distribution for the crossing times is

    f_()≡P̣_()/=ụ/P̣_u(u)/ụ=(d/τ)^2/df_v(d/τ).



In the case of domain walls, we are actually interested in the distribution of perpendicular velocitiesv̌_⊥, see Fig. <ref>.
Note that the Earth is more likely to cross paths with walls that have velocities close to the normal (such objects sweep out a greater volume per unit time).
Ifηis the angle betweenv̌andv̌_⊥, see Eq. (<ref>), then 
the probability of encountering a wall with this angle is proportional tocosη, and can be expressed as

    P̣_η(cosη)=2cosη (̣cosη) = 2 v_⊥/v^2ṿ_⊥.

Therefore, we havef_v_⊥(v̌_⊥)  ṿ_⊥≡P̣_v̌_⊥(v̌_⊥), with 

    P̣_v̌_⊥(v̌_⊥)=∫_v_⊥^∞∫_-π/2^π/2 P̣_η(cosη) P̣_v(v̌) ṿ η̣,

which implies

    f_v_⊥(v̌_⊥)=2 ∫_v_⊥^∞∫_-π/2^π/2 f_v(v̌)v_⊥/v^2 ṿ η̣.

We can further find the angular, scalar, and crossing-time distributions as above, which are also presented in Fig. <ref>.






















§ CLOCK NOISE PROFILES





Here we present a brief overview of the noise characteristics of the GPS satellite clocks.
For more detail, including the analysis for each individual SVN, see the Supplementary Information.







In Table <ref>, we present the average standard deviations of each clock and satellite combination for both first- and second-order differenced data, averaged over all available SVNs and reference clocks.
We also form the autocorrelation function (ACF),

    A^a(j _0)=∑_l=0^J-j-1d^a_l  d^a_l+j/(J-j) (σ^a)^2,

for each clock, where the time-series data{d^a_j}is assumed to be centered around 0, 
andJis the total number of data points for each clock per day. 
For the 30s sampled data,_0=30sandJ=2880.
(Here,σ^ais the standard deviation of the clock data, not the formal error.)
We calculate ACFs for the first- and second-order differenced data (d^(1)andd^(2)).
In Fig. <ref> we show the ACF averaged over all clocks of a specific type between July 2004 and June 2016 for first- and second-order differenced data.
For pure white dataA(0)=1andA()→0for≠0; other noise profiles have distinct ACF forms (see, e.g., Ref. <cit.>).
First-order differencing is sufficient to ensure all the Rb and the block IIF Cs clock time-series are sufficiently stationary, while the block II and IIA Cs clocks require second-order differencing. 











We can also compute the Allan variance for each clock 

    σ_y^2(Δ t)=∑_j=0^J-2l-1( d_j-2d_j+l+d_j+2l)^2/2 l^2_0^2(J-2l),

which is a widely utilized time-domain measure of the frequency stability <cit.>.
Note that the Allan variance is a function of the averaging time, which in our case can be writtenΔt=l _0.
The Allan variance is shown is Fig. <ref>, where we use the non-differenced data,d^(0).








For a given clock,a, we can form the power spectral density (PSD),

    S^a(k)=_0/Jd^a(k)^2,

whered^a(k)=∑_j d^a_j exp(-i2π j k/J)is the discreet Fourier transform  of the time-series data for the clocka. 
The PSD units ares^2/Hz.
Plots of the PSD for each of the clock/satellite combinations are shown in Fig. <ref>, for which we use the singly-differenced data,d^(1).












The periodic spikes that appear in the power spectrum and autocorrelation function (particularly visible for the Rb-IIF satellite clocks, see Figs. <ref> and <ref>) correspond to a 5-minute period, and are technical artefacts traceable to the partitioning method used in the initial JPL data processing.
This has been addressed in recent updates to their software.







§ SPECIFIC DARK-MATTER SIGNALS









Here, we present the specific DM profiles and resultant signals for domain walls, monopoles, and strings.
We then link the generalhparameters back to the specific field parameters for those models.




  
*Thin walls—

The simplest case to consider is thin domain walls.
By “thin”, we mean that the width of the wall is sufficiently small such that it will pass through any clock within the 30 s sampling period,

    d≪∼ 300kms^-1× 30s≈ 10^4 km.

In this case, the profile can be considered to be a delta-functionφ^2(t^a,t)=δ(t^a-t),so that

    s^a_j^(0)=
    	
    		0 		    t ≤ t^a, t^R  
    
    		h^a 		    t^a ≤ t < t^R 
    
    		-h^R 	    t^a > t ≥ t^R 
    
    		h^a-h^R	    t≥ t^a, t^R

(t=j_0).
From the normalization defined in Eq. (<ref>), in the thin wall case, the parameterhcan be linked back to the field parameters as

    h = A^2∑_X_̨X Γ_X.



  
*Gaussian profile walls—

For walls of finite thickness, we assume a Gaussian density profile, with root-mean-square widthd, such that

    φ^2(t^a,t')=v_⊥/d √(π)exp( -v_⊥^2/d^2(t^a-t')^2 ).

The normalization coefficient, which includesv_⊥andd, is chosen purely for convenience so that the maximum accumulated clock bias will beh, in order to be consistent with the thin wall case and because it ishthat is the directly observable parameter.
Then, the integral in Eq. (<ref>) can be expressed in terms of error functions, 

    s^a_j^(0)=
    1/2{ h^a - h^R  +   h^a [ v_⊥(j-t^a)/d]
    
    
    -h^R [ v_⊥(j-t^R)/d] }.

In the Gaussian-profile wall case, the parameterhrelates to the field parameters as

    h = A^2 d √(π)/v_⊥∑_X_̨X Γ_X.



  
*Monopoles—




For monopoles, we assume Gaussian profile spherical objects, and also have to consider the impact parameter,ρ, the distance between the clock and the center of the DM object in the plane perpendicular to ň̂̌, the incident direction of the object.
In this plane, the distance of a clock from ECI0 is given by

    r^a_⊥^2=ř^a^2 - (ň̂̌·ř^a)^2.

If, in this plane, the DM object enters with a perpendicular distance ofRfrom ECI0, and at an angleαmeasured fromz', the projection ofzdown to the plane perpendicular ton̂[ẑ̌' = ẑ̌ - (ẑ̌·ň̂̌)ň̂̌],[4] 

 then the impact parameter is

    ρ^a = √(r^a_⊥^2+R^2-2r^a_⊥Rcosγ^a),

whereγ^a=α-β^ais the angle betweenRandr^a_⊥, andβ^ais the angle thatr^a_⊥makes in the plane perpendicular toň̂̌also measured from thez'-axis, and is given by

    tanβ^a=(ř^a×ẑ̌)·ň̂̌/ř^a·ẑ̌-(ř^a·ň̂̌)(ň̂̌·ẑ̌),

 as shown in Fig. <ref>.
The profile can be expressed as

    φ^2(t^a,t')=v/d √(π)exp( -v^2/d^2(t^a-t')^2 -ρ^a^2/d^2),

and the parameterhis linked to the field parameters as

    h = A^2 d √(π)/v∑_X_̨X Γ_X.



  
*Strings—


The string case is similar to the monopole case, except here the impact parameter is set by the perpendicular distance from each clock to the string.
We assume that on the scale of the GPS network, the string can be modelled as a straight line segment.
For a string that enters from incident directionň̂̌_⊥=-v̌_̌⊥̌/v_⊥(we are interested in the velocity perpendicular to the string), with a perpendicular distance ofRfrom ECI0, at an angleα(measured from thez'-axis toRas above), the impact parameter for each satellite is

    ρ^a =R-r^a_⊥cosγ^a,

where, as above,γ^a=α-β^ais the angle betweenRandr^a_⊥, andβ^ais the angle thatr^a_⊥makes in the plane perpendicular toň̂̌_⊥also measured from thez'-axis. 


Assuming a 2D Gaussian profile with radial widthd, the string profile can be expressed

    φ^2(t^a,t')=v_⊥/d √(π)exp( -v_⊥^2/d^2(t^a-t')^2 -ρ^a^2/d^2).

The parameterhis linked back to the field parameters in the same way as for the Gaussian-profile wall case (<ref>).





[4]When ẑ̌·ň̂̌≈ 1, we instead define ẑ̌' = x̌̂̌ - (x̌̂̌·ň̂̌)ň̂̌.