Search for domain wall dark matter with atomic clocks on board global positioning system satellites

Cosmological observations indicate that dark matter makes up 85% of all matter in the universe yet its microscopic composition remains a mystery. Dark matter could arise from ultralight quantum fields that form macroscopic objects. Here we use the global positioning system as a ~ 50,000 km aperture dark matter detector to search for such objects in the form of domain walls. Global positioning system navigation relies on precision timing signals furnished by atomic clocks. As the Earth moves through the galactic dark matter halo, interactions with domain walls could cause a sequence of atomic clock perturbations that propagate through the satellite constellation at galactic velocities ~ 300 km s−1. Mining 16 years of archival data, we find no evidence for domain walls at our current sensitivity level. This improves the limits on certain quadratic scalar couplings of domain wall dark matter to standard model particles by several orders of magnitude.

D espite the overwhelming cosmological evidence for the existence of dark matter (DM), there is as of yet no definitive evidence for DM in terrestrial experiments. Multiple cosmological observations suggest that ordinary matter makes up only about 15% of the total matter in the universe, with the remaining portion composed of DM 1 . All the evidence for DM (e.g., galactic rotation curves, gravitational lensing, cosmic microwave background) comes from galactic or larger scale observations through the gravitational pull of DM on ordinary matter 1 . Extrapolation from the galactic to laboratory scales presents a challenge because of the unknown nature of DM constituents. Various theories postulate additional nongravitational interactions between standard model (SM) particles and DM. Ambitious programs in particle physics have mostly focused on (so far unsuccessful) searches for weakly interacting massive particle (WIMP) DM candidates with 10-10 3 GeV c −2 masses (c is the speed of light) through their energy deposition in particle detectors 2 . The null results of the WIMP searches have partially motivated an increased interest in alternative DM candidates, such as ultralight fields. These fields, in contrast to particle candidates, act as coherent entities on the scale of an individual detector.
Here we focus on ultralight fields that may cause apparent variations in the fundamental constants of nature. Such variations in turn lead to shifts in atomic energy levels, which may be measurable by monitoring atomic frequencies [3][4][5] . Such monitoring is performed naturally in atomic clocks, which tell time by locking the frequency of externally generated electromagnetic radiation to atomic frequencies. Here, we analyse time as measured by atomic clocks on board global positioning system (GPS) satellites to search for DM-induced transient variations of fundamental constants 6 . In effect we use the GPS constellation as a5 0,000 km-aperture DM detector. Our DM search is one example of using GPS for fundamental physics research. Another recent example includes placing limits on gravitational waves 7 .
GPS works by broadcasting microwave signals from nominally 32 satellites in medium-Earth orbit. The signals are driven by an atomic clock (either based on Rb or Cs atoms) on board each satellite. By measuring the carrier phases of these signals with a global network of specialised GPS receivers, the geodetic community can position stations at the 1 mm level for purposes of investigating plate tectonics and geodynamics 8 . As part of this data processing, the time differences between satellite and station clocks are determined with <0.1 ns accuracy 9 . Such high-quality timing data for at least the past decade are publicly available and are routinely updated. Here we analyse data from the Jet Propulsion Laboratory 10 . A more detailed overview of the GPS architecture and data processing relevant to our search is given in Supplementary  Topological defects may be formed during the cooling of the early universe through a spontaneous symmetry breaking phase transition 11,12 . Technically, this requires the existence of hypothesised self-interacting DM fields, φ. While the exact nature of TDs is model-dependent, the spatial scale of the DM object, d, is generically given by the Compton wavelength of the particles that make up the DM field d ¼ h=ðm φ cÞ, where m φ is the field particle mass, and h is the reduced Plank constant. The fields that are of interest here are ultralight: for an Earth-sized object the mass scale is m φ $ 10 À14 eV c À2 , hence the probed parameter space is complementary to that of WIMP searches 2 , as well as searches for other DM candidates [20][21][22] . Searches for TDs have been performed via their gravitational effects, including gravitational lensing [23][24][25] . Limits on TDs have been placed by the Planck 26 and Background Imaging of Cosmic Extragalactic Polarization 2 (BICEP2) 27 collaborations from fluctuations in the cosmic microwave background. So far the existence of TDs is neither confirmed nor ruled out. The past few years have brought several proposals for TD searches via their non-gravitational signatures 6,[28][29][30][31][32] .
Here, we report the results of the search for domain walls, quasi-2D cosmic structures. As a result, we improve the limits on certain quadratic scalar couplings of domain wall DM to standard model particles by several orders of magnitude.

Results
Domain wall theory and expected signal. We focus on the search for domain walls, since they would leave the simplest DM signature in the data. General signature matching for the vast set of GPS data has proven to be computationally expensive and is in progress. While we interpret our results in terms of domain wall DM, we remark that our search applies equally to the situation where walls are closed on themselves, forming a bubble that has transverse size significantly exceeding the terrestrial scale. The galactic structure formation in that case may occur as per conventional cold dark matter theory 33 , since from the large distance perspective the bubbles of domain walls behave as point-like objects.
We employ the known properties of the DM halo to model the statistics of encounters of the Earth with TDs. Direct measurements 34 of the local dark matter density give 0.3 ± 0.1 GeV cm −3 , and we adopt the value of ρ DM ≈ 0.4 GeV cm −3 for definitiveness. According to the standard halo model, in the galactic rest frame the velocity distribution of DM objects is isotropic and quasi-Maxwellian, with dispersion 35 v ' 290 km s −1 and a cut-off above the galactic escape velocity of v esc ' 550 km s −1 . The Milky Way rotates through the DM halo with the Sun moving at~220 km s −1 towards the Cygnus constellation. For the goals of this work we can neglect the much smaller orbital velocities of the Earth around the Sun (~30 km s −1 ) and GPS satellites around the Earth (~4 km s −1 ). Thereby one may think of a TD wind impinging upon the Earth, with typical relative velocities v g $ 300 km s −1 .
Assuming the standard halo model, the vast majority of events (~95%) would come from the forward-facing hemisphere centred about the direction of the Earth's motion through the galaxy, with typical transit times through the GPS constellation of about 3 min. An example of a domain wall crossing is shown in Fig. 1. Note that we make an additional assumption that the distribution of wall velocities is similar to the standard halo model, which is expected if the gravitational force is the main force governing wall dynamics within the galaxy. However, even if this distribution is somewhat different, the qualitative feature of a TD wind is not expected to change.
A positive DM signal can be visualised as a coordinated propagation of clock glitches at galactic velocities through the GPS constellation, see Fig. 2. The powerful advantage of working with the network is that non-DM clock perturbations do not mimic this signature. The only systematic effect that has propagation velocities comparable to v g is the solar wind 36 , an effect that is simple to exclude based on the distinct directionality from the Sun and the fact that the solar wind does not affect the satellites in the Earth's shadow.
As the nature of non-gravitational interactions of DM with ordinary matter is unknown, we take a phenomenological approach that respects the Lorentz and local gauge invariances. We consider quadratic scalar interactions between the DM objects and clock atoms that can be parameterised in terms of shifts in the effective values of fundamental constants 6 . The relevant combinations of fundamental constants include α ¼ e 2 = hc % 1=137, the dimensionless electromagnetic finestructure constant (e is the elementary charge), the ratio m q /Λ QCD of the light quark mass to the quantum chromodynamics (QCD) energy-scale, and m e and m p , the electron and proton masses. With the quadratic scalar coupling, the relative change in the local value for each such fundamental constant is proportional to the square of the DM field where Γ X is the coupling constant between dark and ordinary matter, with X = α,m e ,m p ,m q /Λ QCD (see Supplementary Note 2 for further details).
As the DM field vanishes outside the TD, the apparent variations in the fundamental constants occur only when the TD overlaps with the clock. This temporary shift in the fundamental constants leads in-turn to a transient shift in the atomic energy levels referenced by the clocks, which may be measurable by ARTICLE monitoring atomic frequencies [3][4][5] . The frequency shift can be expressed as where ω c is the unperturbed clock frequency and K X are known coefficients of sensitivity to effective changes in the constant X for a particular clock transition 37 . It is worth noting that the values of the sensitivity coefficients K X depend on experimental realisation.
Here we compare spatially separated clocks (to be contrasted with the conventional frequency ratio comparisons [3][4][5], and thus our used values of K X somewhat differ from other places in the literature 37 ; full details are presented in Supplementary Note 2. For example, for the microwave frequency 87 Rb clocks on board the GPS satellites, the sensitivity coefficients are where we have introduced the short-hand notation Γ q Γ mq=ΛQCD and Γ e=p 2Γ me À Γ mp , and the effective coupling constant Γ eff ≡ ∑ X K X Γ X . From Eqs. (1) and (2), the extreme TD-induced frequency excursion, δω ext , is related to the field amplitude φ max inside the defect as δω ext ¼ Γ eff ω c φ 2 max . Further, assuming that a particular TD type saturates the DM energy density, we have 6 Here, T is the average time between consecutive encounters of the clock with DM objects, which, for a given ρ DM , depends on the energy density inside the defect 6 which is valid for TDs of any type (monopoles, walls and strings). The frequency excursion is positive for Γ eff > 0, and negative for Γ eff < 0.
The key qualifier for the preceding Eq. (4) is that one must be able to distinguish between the clock noise and DM-induced frequency excursions. Discriminating between the two sources relies on measuring time delays between DM events at network nodes. Indeed, if we consider a pair of spatially separated clocks (Fig. 2), the DM-induced frequency shift Eq. (2) translates into a distinct pattern. The velocity of the sweep is encoded in the time delay between two DM-induced spikes and it must lie within the boundaries predicted by the standard halo model. Generalisation to the multi-node network is apparent (see GPS-specific discussion below). The distributed response of the network encodes the spatial structure and kinematics of the DM object, and its coupling to atomic clocks.
Analysis and search. Working with GPS data introduces several peculiarities into the above discussion (see Supplementary Note 1 for details). The most relevant is that the available GPS clock data are clock biases (i.e., time differences between the satellite and reference clocks) S (0) (t k ) sampled at times (epochs) t k every 30 s. Thus we cannot access the continuously sampled clock frequencies as in Fig. 2. Instead, we formed discretised pseudofrequencies S ð1Þ ðt k Þ S ð0Þ ðt k Þ À S ð0Þ ðt kÀ1 Þ. Then the signal is especially simple if the DM object transit time through a given clock, d/v g , is smaller than the 30-s epoch interval (i.e., thin DM objects with d≲10 4 km, roughly the size of the Earth), since in this case S (1) collapses into a solitary spike at t k if the DM object was encountered during the (t k−1 , t k ) interval. The exact time of interaction within this interval is treated as a free parameter.
One of the expected S (1) signatures for a thin domain wall propagating through the GPS constellation is shown in Fig. 3a. This signature was generated for a domain wall incident with v = 300 km s −1 from the most probable direction. The derivation of the specific expected domain wall signal is presented in Supplementary Note 3, which includes Supplementary Fig. 3. As the DM response of Rb and Cs satellite clocks can differ due to their distinct effective coupling constants Γ eff , we treated the Cs and Rb satellites as two sub-networks, and performed the analysis separately. Within each sub-network we chose the clock on board the most recently launched satellite as the reference because, as a rule, such clocks are the least noisy among all the clocks in orbit.
To search for domain wall signals, we analysed the S (1) GPS data streams in two stages. At the first stage, we scanned all the data from May 2000 to October 2016 searching for the most general patterns associated with a domain wall crossing, without taking into account the order in which the satellites were swept.
cut , and the blue depict S ð1Þ < À S cut , with S ð1Þ cut ¼ 0:18 ns and 0.13 ns, respectively. At the 0.13 ns level, c, this data window would be flagged as a potential event, but not at the 0.18 ns level shown in b. In this case, the potential event c is excluded because the reference clock experiences a much larger perturbation than the rest of the clock network We required at least 60% of the clocks to experience a frequency excursion at the same epoch, which would correspond to when the wall crossed the reference clock (vertical blue line in Fig. 3a). This 60% requirement is a conservative choice based on the GPS constellation geometry, and ensures sensitivity to walls with relative speeds of up to v≲700 km s À1 . Then, we checked if these clocks also exhibit a frequency excursion of similar magnitude (accounting for clock noise) and opposite sign anywhere else within a given time window (red tiles in Fig. 3a). Any epoch for which these criteria were met was counted as a potential event. We considered time windows corresponding to sweep durations through the GPS constellation of up to 15,000 s, which is sufficiently long to ensure sensitivity to walls moving at relative velocities v ≲ 4 km s À1 (given that <0.1% of DM objects are expected to move with velocities outside of this range). Further details of the employed search technique are presented in the Methods section and Supplementary Note 4.
The tiled representation of the GPS data stream depends on the chosen signal cut-off S ð1Þ cut (see Fig. 3). We systematically decreased the cut-off values and repeated the above procedure. Above a certain threshold, S ð1Þ thresh , no potential events were seen. This process is demonstrated for a single arbitrarily chosen data window in Fig. 3b, c. The thresholds for the Rb and Cs subnetworks above which no potential events were seen are S The second stage of the search involved analysing the potential events in more detail, so that we may elevate their status to candidate events if warranted by the evidence. We examined a few hundred potential events that had S (1) magnitudes just below S ð1Þ thresh , by matching the data streams against the expected patterns; one such example is shown in Fig. 3a. At this second stage, we accounted for the ordering and time at which each satellite clock was affected. The velocity vector and wall orientation were treated as free parameters within the bounds of the standard halo model. As a result of this pattern matching, we found that none of these events were consistent with domain wall DM, thus we have found no candidate events at our current sensitivity. Analysing numerous potential events well below S ð1Þ thresh has proven to be substantially more computationally demanding, and is beyond the scope of the current work.

Discussion
As we did not find evidence for encounters with domain walls at our current sensitivity, there are two possibilities: either DM of this nature does not exist, or the DM signals are below our sensitivity. In the latter case we may constrain the possible range of the coupling strengths Γ eff . For the discrete pseudofrequencies, and considering the case of thin domain walls, Eq. (4) becomes Our technique is not equally sensitive to all values for the wall widths, d, or average times between collisions, T . This is directly This work Optical Sr Fig. 5 Constraints on the coupling of dark matter to electromagnetism. Limits (90% confidence level) on the energy scale Λ α as a function of the wall width d and average time between encounters T . The shaded yellow region shows the Global Positioning System limits from this work (assuming Γ α ) Γ q;e=p ), the shaded green region shows the limits derived from an optical Sr clock 38 , and the shaded blue region shows the astrophysical bounds 39 . The solid red line shows the potential discovery reach using the global network of Global Positioning System microwave atomic clocks. For T ≲7 yr, the Global Positioning System reach is limited by the modern Rb block IIF satellite clocks 46 (σ y ð30 sÞ $ 10 À11 ), and for T ≲7 yr, the reach is limited by the older Rb (block IIR, IIA and II) clocks (σ y ð30 sÞ $ 10 À10 ). Compared to more accurate optical clocks, microwave clocks provide additional sensitivity to Λ q and Λ e/p (optical clocks only have sensitivity to Λ α ) NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-01440-4 ARTICLE NATURE COMMUNICATIONS | 8: 1195 | DOI: 10.1038/s41467-017-01440-4 | www.nature.com/naturecommunications taken into account by introducing a sensitivity function, s (d)∈[0,1], that is included in Eq. (5) to determine the final limits at the 90% confidence level. For example, the smallest width is determined by the servo-loop time of the GPS clocks, i.e., by how quickly the clock responds to the changes in atomic frequencies.
In addition, we are sensitive to events that occur less frequently than once every~150 s (so the expected patterns do not overlap), which places the lower bound on T . Further, we incorporate the expected event statistics into Eq. (5). Details are presented in the Methods section.
Our results are presented in Fig. 4. To be consistent with previous literature 6,38 , the limits are presented for the effective energy scale Λ eff 1= ffiffiffiffiffiffiffiffiffiffi Γ eff j j p . Further, on the assumption that the coupling strength Γ α dominates over the other couplings in the linear combination in Eq. (3), we place limits on Λ α . The resulting limits are shown in Fig. 5, together with existing constraints 38,39 . For certain parameters, our limits exceed the 10 7 TeV level; astrophysical limits 39 on Λ α , which come from stellar and supernova energy-loss observations 40,41 , have not exceeded 10 TeV. The derived constraints on Λ α can be translated into a limit on the transient variation of the fine-structure constant, which for d = 10 4 km corresponds to δα=α≲10 À12 . Because of the scaling of the constraints on Λ X , this result is independent of T , and scales inversely with d (within the region of applicability). It is worth contrasting this constraint with results from the searches for slow linear drifts of fundamental constants. For example, the search 5 resulting in the most stringent limits on long-term drifts of α was carried out over a year and led to δα α ≲3 10 À17 . Such long-term limits apply only for very thick walls of thickness d ) v g 1 yr $ 10 10 km, which are outside our present discovery reach.
Further, by combining our results from the Rb and Cs GPS sub-networks with the recent limits on Λ α from an optical Sr clock 38 , we also place independent limits on Λ e/p , and Λ q ; for details, see Supplementary Note 5. These limits are presented in Fig. 6 as a function of the average time between events. For certain values of the d and T parameters, we improve current bounds on Λ e/p by a factor of~10 5 and for the first time establish limits on Λ q .
While we have improved the current constraints on DMinduced transient variation of fundamental constants by several orders of magnitude, it is possible that DM events remain undiscovered in the data noise. Our current threshold S ð1Þ thresh is larger than the GPS data noise by a factor of~5-20, depending on which clocks/time periods are examined. By applying a more sophisticated statistical approach with greater computing power, we expect to improve our sensitivity by up to two orders of magnitude. Indeed, the sensitivity of the search is statistically determined by the number of clocks in the network, N clocks , and the Allan deviation 42 , σ y (τ 0 ), evaluated at the data sampling interval τ 0 = 30 s reads, or, combining with Eq. (5), Note that this estimate differs from a previous estimate 6 , since while arriving at Eq. (7), we assumed a more realistic white frequency noise (instead of white phase noise). The projected discovery reach of GPS data analysis is presented in Fig. 5. Prospects for the future include incorporating substantially improved clocks on next-generation satellites, increasing the network density with other Global Navigation Satellite Systems, such as European Galileo, Russian Global Navigation Satellite System (GLONASS), and Chinese BeiDou, and including networks of laboratory clocks 38,43 . Such an expansion can bring the total number of clocks to~100. Moreover, the GPS search can be extended to other TD types (monopoles and strings), as well as different DM models, such as virialized DM fields 32,44 .
In summary, by using the GPS as a dark matter detector, we have substantially improved the current limits on DM domain wall induced transient variation of fundamental constants. Our approach relies on mining an extensive set of archival data, using existing infrastructure. As the direct DM searches are widening to include alternative DM candidates, it is anticipated that the mining of time-stamped archival data, especially from laboratory precision measurements, will play an important role in verifying or excluding predictions of various DM models 45 . In the future, our approach can be used for a DM search with nascent networks  38 to place assumption-free limits on Λ e/p and Λ q . The blue region shows the astrophysical bounds for Λ e/p ; note that Λ q was previously unconstrained 39 of laboratory atomic clocks that are orders of magnitude more accurate than the GPS clocks 43 .

Methods
Crossing duration distribution. Before placing limits on Γ eff , we must account for the fact that we do not have equal sensitivity to each domain wall width, d, or equivalently crossing durations, where v ⊥ is the component of the velocity perpendicular to the face of the DM wall. This is due in part to aspects of the clock hardware and operation, the timeresolution (data sampling frequency), and the employed search method. Therefore, for a given d, we must determine the proportion of events that have crossing durations within the range τ min to τ max , where τ min(max) is the minimum (maximum) crossing duration for which this method is sensitive.
There are two factors that determine τ min . The first is the servo-loop time-the fastest perturbation that can be recorded by the clock. This servo-loop time is manually adjusted by military operators and is not available to us at a given epoch, however, it is known 46,47 to be within 0.01 and 0.1 s. As such, we consider the best and worse case scenarios: The second condition that affects τ min is the clock degeneracy: the employed GPS data set has only 30 s resolution, so any clocks which are affected within 30 s of the reference clock will not exhibit any DM-induced frequency excursion in their data; see satellites 15 and 16 in Fig. 3a. For <60% of the clocks to experience the jump, the (thin-wall) DM object would have to be travelling at over 700 km s −1 , which is close to the galactic escape velocity (for head-on collisions), so the degeneracy does not affect the derived limits in a substantial way. (In fact, assuming the standard halo model, <0.1% of events are expected to have v ⊥ > 700 km s −1 .) This velocity corresponds to a crossing duration for the entire network of~70 s. Transforming this to crossing duration for a single clock, τ, amounts to multiplying by the ratio d/(D GPS ): For walls thicker than d~100 km, τ min is determined by the d/(D GPS )70 s term.
As to the maximum crossing duration, there are also two factors that affect τ max . First, the wall must pass each clock in less than the sampling interval of 30 s-this is the condition for the wall to be considered thin: If a wall takes longer than 30 s to pass by a clock, the simple single-data-point signals shown in Fig. 3 would become more complicated, and would require a more-detailed pattern-matching technique. Second, we only consider time windows, J w , of a certain size in our analysis (see Supplementary Note 4). If a wall moves so slowly that it does not sweep all the clocks within this window, the event would be missed: Therefore, the overall expression for τ max is: Making J w large, however, also tends to increase S ð1Þ thresh (since there is a higher chance that a large window will satisfy the condition for a potential event). By performing the analysis for multiple values for J w , we can probe the largest portion of the parameter space; for further details see Supplementary Note 4 and Supplementary Tables 2 and 3. In this work, we consider windows of J w up to 500 epochs (15,000 s), which corresponds to a minimum velocity of~4 km s −1 , which is roughly the orbit speed of the satellites. This has a negligible effect on our sensitivity, since <0.1% of walls are expected to have v ⊥ < 4 km s −1 .
Domain wall width sensitivity. Assuming the standard halo model, the relative scalar velocity distribution of DM objects that cross the GPS network is quasi-