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

TL;DR

This study uses 16 years of GPS atomic clock data to search for domain wall dark matter, setting new constraints on their interactions with standard model particles.

Contribution

It introduces a novel method of using GPS satellite atomic clocks as a large-scale detector for topological defect dark matter, improving existing limits.

Findings

No evidence of domain wall dark matter was found.

Constraints on quadratic scalar couplings were significantly improved.

GPS data can serve as a sensitive detector for macroscopic dark matter objects.

Abstract

Cosmological observations indicate that 85% of all matter in the Universe is dark matter (DM), yet its microscopic composition remains a mystery. One hypothesis is that DM arises from ultralight quantum fields that form macroscopic objects such as topological defects. Here we use GPS as a ~ 50,000 km aperture DM detector to search for such defects in the form of domain walls. GPS navigation relies on precision timing signals furnished by atomic clocks hosted on board GPS satellites. As the Earth moves through the galactic DM halo, interactions with topological defects could cause atomic clock glitches that propagate through the GPS satellite constellation at galactic velocities ~ 300 km/s. Mining 16 years of archival GPS data, we find no evidence for DM in the form of domain walls at our current sensitivity level. This allows us to improve the limits on certain quadratic scalar couplings of domain wall DM to standard model particles by several orders of magnitude.

Tables

Table 1. Table 1: Typical standard deviations (in ns) of the S ( 1 ) superscript 𝑆 1 S^{(1)} clock data over the 16 years analysed in this work. The standard deviations vary significantly between the satellite generations  1 (II, IIA, IIR, and IIF), and also depend on whether the ground based H-maser clock was used as reference (in this case USN3  2 ), or one of the other GPS satellite clocks.

GPS Block	Reference
	USN3	RbIIR
RbII	0.047	0.070
RbIIA	0.038	0.074
RbIIR	0.073	0.097
RbIIF	0.013	0.067
GPS Block	Reference
	USN3	CsIIA
CsII	0.081	0.112
CsIIA	0.088	0.124
CsIIF	0.098	0.128

Table 2. Table 2: An illistrative subset of the results of the stage 1 analysis, showing the number of potential events for the given cut-off values S cut ( 1 ) subscript superscript 𝑆 1 cut S^{(1)}_{\rm cut} (using a 90% C.L. range) for the Rb satellite GPS sub-network for various time window lengths J w subscript 𝐽 𝑤 J_{w} . S thresh ( 1 ) subscript superscript 𝑆 1 thresh S^{(1)}_{\rm thresh} is defined as the smallest S cut ( 1 ) subscript superscript 𝑆 1 cut S^{(1)}_{\rm cut} value with zero potential events. The analysis covers a total effective period of 16.5 years, and includes a total of 131 353 clock-days of data.

	$J_w$ (s)
$S_{\mathrm{cut}}^{\left(1\right)}$ (ns)	90	270	450	630	810	990	1170
0.1	4295864	4387849	4387849	4387849	4387849	4387849	4387849
0.15	682150	973116	973133	973133	973133	973133	973133
0.2	69619	210761	211120	211127	211127	211127	211127
0.25	4969	46234	47214	47310	47322	47328	47330
0.3	363	10233	11727	11906	11955	11974	11994
0.35	23	1331	2290	2594	2680	2724	2737
0.4	1	97	235	345	385	413	437
0.45	0	5	20	29	39	45	51
0.5	0	1	3	4	4	5	5
0.55	0	0	1	1	1	1	1
0.6	0	0	0	1	1	1	1
0.65	0	0	0	0	0	0	0

Table 3. Table 3: As in Table  2 , but for the Cs GPS sub-network. Here, the analysis covers a total effective period of 10.5 years, and includes a total of 37 099 clock-days of data.

	$J_w$ (s)
$S_{\mathrm{cut}}^{\left(1\right)}$ (ns)	90	270	450	630	810	990	1170
0.1	6807411	6885078	6885078	6885078	6885078	6885078	6885078
0.15	2283645	3123230	3123326	3123326	3123326	3123326	3123326
0.2	476928	1490303	1492203	1492225	1492225	1492225	1492225
0.25	60719	552654	566282	566538	566555	566559	566559
0.3	6525	148529	171379	173184	173369	173395	173402
0.35	512	28240	40620	43181	43822	43993	44050
0.4	39	4264	7660	8982	9513	9742	9845
0.45	3	495	1087	1442	1655	1777	1883
0.5	0	52	125	182	222	257	285
0.55	0	8	16	24	28	33	41
0.6	0	1	2	3	3	4	4
0.65	0	1	1	1	1	1	1
0.7	0	0	0	0	0	0	0

Equations

$$\frac{\delta X(\mathbf{r},t)}{X}={\Gamma }_X\phi (\mathbf{r},t{)}^2,$$

$$\frac{\delta \omega (\mathbf{r},t)}{{\omega }_c}=\sum _X{K}_X\frac{\delta X(\mathbf{r},t)}{X},$$

$$\frac{\delta \omega }{{\omega }_c}(\mathrm{Rb})=(4.34{\Gamma }_{\alpha }-0.019{\Gamma }_q+{\Gamma }_{e/p}){\phi }^2\equiv {\Gamma }_{\mathrm{eff}}^{(\mathrm{Rb})}{\phi }^2,$$

$$\frac{\delta {\omega }_{\mathrm{ext}}}{{\omega }_c}={\Gamma }_{\mathrm{eff}}\hslash c{\rho }_{\mathrm{DM}}{v}_g\mathcal{T}d,$$

$$\left|{\Gamma }_{\mathrm{eff}}\right|<\frac{S_{\mathrm{thresh}}^{(1)}}{\hslash c\sqrt{\pi }{\rho }_{\mathrm{DM}}\mathcal{T}{d}^2}.$$

$$\frac{\delta \alpha }{\alpha }=\hslash c{\rho }_{\mathrm{DM}}{v}_g\frac{\mathcal{T}d}{{\Lambda }_{\alpha }^2{K}_{\alpha }},$$

$$S^{(1)}\gtrsim \frac{{\sigma }_y({\tau }_0){\tau }_0}{\sqrt{N_{\mathrm{clocks}}}},$$

$${\Lambda }_X\lesssim d\sqrt{\frac{\hslash c{\rho }_{\mathrm{DM}}\mathcal{T}{K}_X\sqrt{N_{\mathrm{clocks}}}}{{\sigma }_y({\tau }_0){\tau }_0}}.$$

$$S^{(1)}(t_k)=S^{(0)}(t_k)-S^{(0)}(t_{k-1}),$$

$${\mathcal{L}}_{\mathrm{int}}={\mathcal{L}}^{\mathrm{PS}}+{\mathcal{L}}^{{\mathrm{S}}^1}+{\mathcal{L}}^{{\mathrm{S}}^2}+\dots ,$$

$$-{\mathcal{L}}^{\mathrm{DM}-\mathrm{SM}}=$$

$${\phi }^2\left(\mathbf{r},t\right)\left({\Gamma }_f{m}_f{c}^2{\bar{\psi }}_f{\psi }_f+\frac{{\Gamma }_{\alpha }}{4}{F}_{\mu \nu }{F}^{\mu \nu }+\dots \right),$$

$$-{\mathcal{L}}^{\mathrm{SM}}=m_f{c}^2{\bar{\psi }}_f{\psi }_f+\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }+\dots ,$$

$${\alpha }^{\mathrm{eff}}=\alpha \left(1+{\Gamma }_{\alpha }{\phi }^2\right),$$

$$m_{e,p}^{\mathrm{eff}}=m_{e,p}\left(1+{\Gamma }_{m_{e,p}}{\phi }^2\right),$$

$${\left(\frac{m_q}{{\Lambda }_{\mathrm{QCD}}}\right)}^{\mathrm{eff}}=\frac{m_q}{{\Lambda }_{\mathrm{QCD}}}\left(1+{\Gamma }_{m_q/{\Lambda }_{\mathrm{QCD}}}{\phi }^2\right),$$

$${\rho }_{\mathrm{inside}}={\phi }_{\max }^2/(\hslash c{d}^2),$$

$$d=\hslash /(m_{\phi }c).$$

$$\frac{\delta \omega (t)}{{\omega }_c}=\sum _X{K}_X\frac{\delta X(t)}{X}\equiv {\Gamma }_{\mathrm{eff}}{\phi }^2,$$

$$\frac{\delta X(t)}{X}={\Gamma }_X{\phi }^2,$$

$${\Gamma }_{\mathrm{eff}}\equiv \sum _X{K}_X{\Gamma }_X,$$

$$\frac{\delta \omega }{{\omega }_c}(\mathrm{Rb})={\phi }^2(4.34{\Gamma }_{\alpha }-0.019{\Gamma }_q+{\Gamma }_{e/p}),$$

$$\frac{\delta \omega }{{\omega }_c}(\mathrm{Cs})={\phi }^2(4.83{\Gamma }_{\alpha }+0.002{\Gamma }_q+{\Gamma }_{e/p}),$$

$$\frac{\delta \omega }{{\omega }_c}(\mathrm{S}{\mathrm{r}}_{\mathrm{optical}})={\phi }^{2}1.06{\Gamma }_{\alpha },$$

$$S^{(0)}_{i}(t)=\varphi_{\rm max}^{2}\int_{-\infty}^{t}\Bigg{[}{\Gamma_{{\rm eff,}i}}\,\exp\left(\frac{-(t^{\prime}-t^{\times}_{i})^{2}}{\tau^{2}}\right)\qquad\qquad$$

$$\qquad\qquad-{\Gamma_{{\rm eff,}R}}\,\exp\left(\frac{-(t^{\prime}-t^{\times}_{R})^{2}}{\tau^{2}}\right)\Bigg{]}dt^{\prime},$$

$$\tau \equiv \frac{d}{v_{\perp }},$$

$$S^{(1)}={\phi }_{\max }^2{\Gamma }_{\mathrm{eff}}\sqrt{\pi }\tau ,$$

$${\phi }_{\max }^2=\hslash c{\rho }_{\mathrm{DM}}{d}^2\frac{\mathcal{T}}{{\tau }_{\mathrm{avg}}},$$

$$S_{\mathrm{avg}}^{(1)}=\hslash c{\Gamma }_{\mathrm{eff}}\sqrt{\pi }{\rho }_{\mathrm{DM}}{d}^2\mathcal{T}.$$

Full text

Supplementary Information

Here we present details of the GPS architecture, data acquisition, and data analysis. Section S.1 presents an overview of the relevant aspects of GPS. Section S.2 outlines the relevant theoretical background and presents the details of the atomic clock response to passing topological defect dark matter. Section S.3 discusses the sought domain wall dark matter signals, and section S.4 presents the details of our search method for their transient signatures in the GPS atomic clock data. Section S.5 describes the sensitivity of our approach to different regions of the parameter space, and section S.6 presents the resulting limits.

Despite 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 non-gravitational interactions between Standard Model (SM) particles and DM. Ambitious programs in particle physics have mostly focused on (so far unsuccessful) searches for WIMP (Weakly Interacting Massive Particle) DM candidates with $10-10^{3}\,\mathrm{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 of 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 analyze 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 a $\sim$ 50,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\,\mathrm{ns}$ accuracy 9. Such high-quality timing data for at least the past decade are publicly available and are routinely updated. Here we analyze 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 section S.1 of the Supplementary Material.

The large aperture of the GPS network is well suited to search for macroscopic DM objects or “clumps”. Examples of clumpy DM are numerous: topological defects (TDs) 11, 12, $Q$-balls 13, 14, 15, solitons 16, 17, axion stars 18, 19, and other stable objects formed due to dissipative interactions in the DM sector. For concreteness, we interpret our results in terms of TDs. Each TD type (monopoles, strings, or domain walls) would exhibit a transient in GPS data with a distinct signature. General signature matching for the vast set of GPS data has proven to be computationally expensive and is in progress. Here, we report the results of the search for domain walls, quasi-2D cosmic structures, since domain walls would leave the simplest DM signature in the data. An example of a domain wall crossing is shown in Fig. 1. 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 theory20, since from the large distance perspective the bubbles of domain walls behave as point-like objects.

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, $\varphi$. 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=\hbar/(m_{\varphi}c)$, where $m_{\varphi}$ is the field particle mass, and $\hbar$ is the reduced Plank constant. The fields that are of interest here are ultralight: for an Earth-sized object the mass scale is $m_{\varphi}\sim 10^{-14}\,\mathrm{eV}\,c^{-2}$, hence the probed parameter space is complementary to that of WIMP searches 2, as well as searches for other DM candidates 21, 22, 23. Searches for TDs have been performed via their gravitational effects, including gravitational lensing 24, 25, 26. Limits on TDs have been placed by Planck 27 and BICEP2 28 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, 29, 30, 31, 32, 33.

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\pm 0.1\,\mathrm{GeV}\,\mathrm{cm}^{-3}$, and we adopt the value of $\rho_{\mathrm{DM}}\approx 0.4\,\mathrm{GeV}\,\mathrm{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\simeq 290\,$km s*-1* and a cut-off above the galactic escape velocity of $v_{\rm esc}\simeq 550\,$km s*-1*. The Milky Way rotates through the DM halo with the Sun moving at $\sim$ $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 ($\sim$ $30\,{\rm km\,s}^{-1}$) and GPS satellites around the Earth ($\sim$ $4\,{\rm km\,s}^{-1}$). Thereby one may think of a TD “wind” impinging upon the Earth, with typical relative velocities $v_{g}\sim 300\,$km s*-1*. Assuming the standard halo model, the vast majority of events ($\sim$ 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 three minutes. A positive DM signal can be visualised as a coordinated propagation of clock “glitches” at galactic velocities through the GPS constellation, see Fig. 2. 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. 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 $\alpha=e^{2}/\hbar c\approx 1/137$, the dimensionless electromagnetic fine-structure constant ($e$ is the elementary charge), the ratio $m_{q}/\Lambda_{\rm 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

[TABLE]

where $\Gamma_{X}$ is the coupling constant between dark and ordinary matter, with $X=\alpha,\,m_{e},\,m_{p},\,m_{q}/\Lambda_{\rm QCD}$ (see section S.2 of the Supplementary Information for 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 frequency referenced by the clocks, which may be measurable by monitoring atomic frequencies 3, 4, 5. The frequency shift can be expressed as

[TABLE]

where $\omega_{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 Ref. 37; full details are presented in section S.2 of the Supplementary Information. For example, for the microwave frequency 87Rb clocks on board the GPS satellites, the sensitivity coefficients are

[TABLE]

where we have introduced the short-hand notation $\Gamma_{q}\equiv\Gamma_{m_{q}/\Lambda_{\rm QCD}}$ and $\Gamma_{e/p}\equiv 2\Gamma_{m_{e}}-\Gamma_{m_{p}}$, and the effective coupling constant $\Gamma_{\mathrm{eff}}\equiv\sum_{X}K_{X}\Gamma_{X}$.

From Eqs. (1) and (2), the extreme TD-induced frequency excursion, $\delta{\omega}_{\mathrm{ext}}$, is related to the field amplitude $\varphi_{\rm max}$ inside the defect as $\delta{\omega}_{\mathrm{ext}}=\Gamma_{\mathrm{eff}}{\omega_{c}}{\varphi_{\rm max}^{2}}$. Further, assuming that a particular TD type saturates the DM energy density, we have 6 $\varphi_{\rm max}^{2}=\hbar c\rho_{\mathrm{DM}}\mathcal{T}v_{g}d$. Here, $\mathcal{T}$ is the average time between consecutive encounters of the clock with DM objects, which, for a given $\rho_{\mathrm{DM}}$, depends on the energy density inside the defect 6 $\rho_{\rm inside}=\rho_{\rm DM}{\mathcal{T}v_{g}}/{d}.$ Thus the expected DM-induced fractional frequency excursion reads

[TABLE]

which is valid for TDs of any type (monopoles, walls, and strings). The frequency excursion is positive for $\Gamma_{\mathrm{eff}}>0$, and negative for $\Gamma_{\mathrm{eff}}<0$.

The key qualifier for the preceding equation (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 (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. Generalization 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.

Working with GPS data introduces several peculiarities into the above discussion (see section S.1 of the Supplementary Information 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 seconds. Thus we cannot access the continuously sampled clock frequencies as in Fig. 2. Instead, we formed discretised pseudo-frequencies $S^{(1)}(t_{k})\equiv 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-second epoch interval (i.e., “thin” DM objects with $d\lesssim 10^{4}\,\mathrm{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. 3(a). 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 section S.3 of the Supplementary Information. Since the DM response of Rb and Cs satellite clocks can differ due to their distinct effective coupling constants $\Gamma_{\mathrm{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 analyzed 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. 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. 3(a)). 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\lesssim 700\,{\rm 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. 3(a)). 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\gtrsim 4\,{\rm km\,s}^{-1}$ (given that $<$ 0.1% of DM objects are expected to move with velocities outside of this range). The full details of the employed search technique are presented in section S.4 of the Supplementary Information.

The tiled representation of the GPS data stream depends on the chosen signal cut-off $S^{(1)}_{\mathrm{cut}}$ (see Fig. 3). We systematically decreased the cut-off values and repeated the above procedure. Above a certain threshold, $S^{(1)}_{\mathrm{thresh}}$, no potential events were seen. This process is demonstrated for a single arbitrarily chosen data window in Figs. 3(b) and 3(c). The thresholds for the Rb and Cs subnetworks above which no potential events were seen are $S^{(1)}_{\mathrm{thresh}}(\mathrm{Rb})=0.48\,{\rm ns}$ and $S^{(1)}_{\mathrm{thresh}}(\mathrm{Cs})=0.56\,{\rm ns}$ for $v\approx 300\,$km s*-1* sweeps.

The second stage of the search involved analyzing 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)}_{\mathrm{thresh}}$, by matching the data streams against the expected patterns; one such example is shown in Fig. 3(a). 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)}_{\rm thresh}$ has proven to be substantially more computationally demanding, and is beyond the scope of the current work.

Since we did not find evidence for encounters with domain walls at our current sensitivity, there are two possibilities: (i) DM of this nature does not exist, or (ii) the DM signals are below our sensitivity. In the latter case we may constrain the possible range of the coupling strengths $\Gamma_{\mathrm{eff}}$. For the discrete pseudo-frequencies, and considering the case of thin domain walls, Eq. (4) becomes

[TABLE]

Our technique is not equally sensitive to all values for the wall widths, $d$, or average times between collisions, $\mathcal{T}$. This is directly taken into account by introducing a sensitivity function that is combined with 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 $\sim$ 150 s (so the expected patterns do not overlap), which places the lower bound on $\mathcal{T}$. Further, we incorporate the expected event statistics into Eq. (5). Details can be found in section S.5 of the Supplementary Information.

Our results are presented in Fig. S.6. To be consistent with previous literature 6, 38, the limits are presented for the effective energy scale $\Lambda_{\mathrm{eff}}\equiv 1/\sqrt{\left|\Gamma_{\mathrm{eff}}\right|}$. Further, on the assumption that the coupling strength $\Gamma_{\alpha}$ dominates over the other couplings in the linear combination (S.12), we place limits on $\Lambda_{\alpha}$. The resulting limits are shown in Fig. 5, together with existing constraints 38, 39. For certain parameters, our limits exceed the $10^{7}\,{\rm TeV}$ level; astrophysical limits 39 on $\Lambda_{\alpha}$, which come from stellar and supernova energy-loss observations 41, 42, have not exceeded $\sim$ $10\,{\rm TeV}$.

The derived constraints on $\Lambda_{\alpha}$ can be translated into a limit on the transient variation of the fine-structure constant,

[TABLE]

which for $d=10^{4}\,{\rm km}$ corresponds to ${\delta\alpha}/{\alpha}\lesssim 10^{-12}$. Due to the scaling of the constraints on $\Lambda_{X}$, this result is independent of $\mathcal{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 $\alpha$ was carried out over a year and led to $\frac{\delta\alpha}{\alpha}\lesssim 3\times 10^{-17}$. Such limits may apply only for very thick walls $d\gg v_{g}\times 1\,{\rm yr}\sim 10^{10}\,{\rm 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 $\Lambda_{\alpha}$ from an optical Sr clock 38, we also place independent limits on $\Lambda_{e/p}$, and $\Lambda_{q}$; for details, see section S.6 of the Supplementary Information. These limits are presented in Fig. 6 as a function of the average time between events. For certain values of the $d$ and $\mathcal{T}$ parameters, we improve current bounds on $\Lambda_{e/p}$ by a factor of $\sim 10^{5}$ and for the first time establishes limits on $\Lambda_{q}$.

While we have improved the current constraints on DM-induced 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)}_{\rm thresh}$ is larger than the GPS data noise by a factor of $\sim 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_{\rm clocks}$ and the Allan deviation 43, $\sigma_{y}(\tau_{0})$, evaluated at the data sampling interval $\tau_{0}=30\,{\rm s}$ reads,

[TABLE]

or, combining with Eq. (5),

[TABLE]

Note that this estimate differs from that in Ref. 6, since while ariving 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 GLONASS, and Chinese BeiDou, and including networks of laboratory clocks 38, 44. Such an expansion can bring the total number of clocks to $\sim$ 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 33, 45.

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 46. In the future, our approach can be used for a DM search with nascent networks of laboratory atomic clocks that are orders of magnitude more accurate than the current GPS clocks 44.

• Bertone et al. 2005 G. Bertone, D. Hooper,  and J. Silk, Phys. Rep. 405, 279 (2005).
• Liu et al. 2017 J. Liu, X. Chen, and X. Ji, Nat. Phys. 13, 212 (2017).
• Rosenband et al. 2008 T. Rosenband et al., Science 319, 1808 (2008).
• Huntemann et al. 2014 N. Huntemann, B. Lipphardt, C. Tamm, V. Gerginov, S. Weyers, and E. Peik, Phys. Rev. Lett. 113, 210802 (2014).
• Godun et al. 2014 R. M. Godun et al., Phys. Rev. Lett. 113, 210801 (2014).
• Derevianko and Pospelov 2014 A. Derevianko and M. Pospelov, Nat. Phys. 10, 933 (2014).
• Aoyama et al. 2014 S. Aoyama, R. Tazai, and K. Ichiki, Phys. Rev. D 89, 067101 (2014).
• Blewitt 2015 G. Blewitt, in Treatise on Geophysics (Elsevier, Oxford, 2015) pp. 307–338.
• Ray and Senior 2005 J. Ray and K. Senior, Metrologia 42, 215 (2005).
• D. Murphy et al. 2015 D. Murphy et al., JPL Analysis Center Technical Report (in IGS Technical Report) 2015, 77 (2015).
• Kibble 1980 T. Kibble, Phys. Rep. 67, 183 (1980).
• Vilenkin 1985 A. Vilenkin, Phys. Rep. 121, 263 (1985).
• Coleman 1985 S. Coleman, Nucl. Phys. B 262, 263 (1985).
• Kusenko and Steinhardt 2001 A. Kusenko and P. J. Steinhardt, Phys. Rev. Lett. 87, 141301 (2001).
• Lee et al. 1989 K. Lee, J. A. Stein-Schabes, R. Watkins, and L. M. Widrow, Phys. Rev. D 39, 1665 (1989).
• Marsh and Pop 2015 D. J. E. Marsh and A.-R. Pop, Mon. Not. R. Astron. Soc. 451, 2479 (2015).
• Schive et al. 2014 H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nat. Phys. 10, 496 (2014).
• Hogan and Rees 1988 C. Hogan and M. Rees, Phys. Lett. B 205, 228 (1988).
• Kolb and Tkachev 1993 E. W. Kolb and I. I. Tkachev, Phys. Rev. Lett. 71, 3051 (1993).
• Blumenthal et al. 1984 G. R. Blumenthal, S. M. Faber, J. R. Primack, and M. J. Rees, Nature 311, 517 (1984).
• Bozek et al. 2015 B. Bozek, D. J. E. Marsh, J. Silk, and R. F. G. Wyse, Mon. Not. R. Astron. Soc. 450, 209 (2015).
• Arvanitaki and Dubovsky 2011 A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83, 044026 (2011).
• Arvanitaki and Geraci 2014 A. Arvanitaki and A. A. Geraci, Phys. Rev. Lett. 113, 161801 (2014).
• Schneider et al. 1992 P. Schneider, J. Ehlers, and E. E. Falco, Gravitational Lenses (Springer-Verlag, Berlin Heidelberg, 1992).
• Vilenkin and Shellard 1994 A. Vilenkin and E. Shellard, Cosmic Strings and Other Topological Defects (Cambridge University Press, Cambridge, 1994).
• Cline 2001 D. B. Cline, ed., Sources and Detection of Dark Matter and Dark Energy in the Universe. Fourth International Symposium Held at Marina del Rey, CA, USA February 23–25, 2000 (Springer, Berlin Heidelberg, 2001).
• The Planck Collaboration 2014 The Planck Collaboration, Astron. Astrophys. 571, A25 (2014).
• The BICEP2 Collaboration 2014 The BICEP2 Collaboration, Phys. Rev. Lett. 112, 241101 (2014).
• Pospelov et al. 2013 M. Pospelov, S. Pustelny, M. P. Ledbetter, D. F. J. Kimball, W. Gawlik, and D. Budker, Phys. Rev. Lett. 110, 021803 (2013).
• Pustelny et al. 2013 S. Pustelny, D. F. J. Kimball, C. Pankow, M. P. Ledbetter, P. Wlodarczyk, P. Wcisło, M. Pospelov, J. R. Smith, J. Read, W. Gawlik, and D. Budker, Ann. Phys. 525, 659 (2013).
• Hall et al. 2016 E. D. Hall, T. Callister, V. V. Frolov, H. Müller, M. Pospelov, and R. X. Adhikari,  (2016), arXiv:1605.01103.
• Stadnik and Flambaum 2014 Y. V. Stadnik and V. V. Flambaum, Phys. Rev. Lett. 113, 151301 (2014).
• Kalaydzhyan and Yu 2017 T. Kalaydzhyan and N. Yu,  (2017), arXiv:1705.05833.
• Bovy and Tremaine 2012 J. Bovy and S. Tremaine, Astrophys. J. 756, 89 (2012).
• Freese et al. 2013 K. Freese, M. Lisanti, and C. Savage, Rev. Mod. Phys. 85, 1561 (2013).
• Meyer-Vernet 2007 N. Meyer-Vernet, Basics of the Solar Wind (Cambridge University Press, Cambridge, 2007).
• Flambaum and Dzuba 2009 V. V. Flambaum and V. A. Dzuba, Can. J. Phys. 87, 25 (2009).
• Wcisło et al. 2016 P. Wcisło, P. Morzyński, M. Bober, A. Cygan, D. Lisak, R. Ciuryło, and M. Zawada, Nat. Astron. 1, 0009 (2016).
• Olive and Pospelov 2008 K. A. Olive and M. Pospelov, Phys. Rev. D 77, 043524 (2008).
• Griggs et al. 2015 E. Griggs, E. R. Kursinski, and D. Akos, Radio Sci. 50, 813 (2015).
• Raffelt 1999 G. G. Raffelt, Annu. Rev. Nucl. Part. Sci. 49, 163 (1999).
• K. S. Hirata et al. 1988 K. S. Hirata et al., Phys. Rev. D 38, 448 (1988).
• Barnes et al. 1971 J. A. Barnes et al., IEEE Trans. Instrum. Meas. IM-20, 105 (1971).
• Riehle 2017 F. Riehle, Nat. Photonics 11, 25 (2017).
• Derevianko 2016 A. Derevianko,  (2016), arXiv:1605.09717.
• Budker and Derevianko 2015 D. Budker and A. Derevianko, Phys. Today 68, 10 (2015).
• 47 See also the Supplementary Information, which contains figures S.1 – S.10, tables 1 – 3, and references 1 – 35.

**Acknowledgements We thank A. Sushkov and Y. Stadnik for discussions, and J. Weinstein for his comments on the manuscript. We acknowledge the International GNSS Service for the GPS data acquisition. We used JPL’s GIPSY software for the orbit reference frame conversion. This work was supported in part by the US National Science Foundation grant PHY-1506424. **

Data Availability— We used publicly available GPS timing and orbit data for the past 16 years from the Jet Propulsion Laboratory 10.

∗andrei@unr.edu.

S.1 GPS data processing and clock estimation

The Global Positioning System (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. While each satellite may host multiple Rb and Cs backup clocks, at any given time only one of these clocks drives the GPS signals, which are transmitted on carrier waves at both L1 ($1.57542\,{\rm GHz}$) and L2 ($1.2276\,{\rm GHz}$) bands. Superimposed on the carrier waves are streams of pseudo-random bits generated by flipping the sign of the electric field. A geodetic GPS receiver can sample the dual-frequency signals simultaneously from all GPS satellites in view (typically 8 to 10) at user-specified intervals (typically 1 to 30 seconds). At every such interval, for each satellite in view, timing measurements are made (according to the receiver clock) of the peak cross-correlation of the incoming signal with the receiver’s replica model of the signal. At both L1 and L2 frequencies, two types of data are generated including a pseudorange (group delay, using the bits) and a carrier phase (phase delay, using the carrier wave). The term pseudorange is used because it is a measure of delay that is biased by the receiver clock. This bias cancels when differencing data between pairs of satellites. The ionospheric delay of $\sim$ $10\,{\rm ns}$ is calibrated with $\sim$ $0.02\,{\rm ns}$ precision by forming a specific “ionosphere-free” linear combination of the data at L1 and L2 frequencies, hence the purpose of the dual-frequency system.

Measurement precision in such cross-correlation systems tends to scale with the relevant wavelength of the signal. In the case of pseudorange, precision scales with the time interval between bit transitions. As a consequence, the pseudorange precision is typically $\sim 2\,{\rm ns}$, whereas in contrast, carrier phase is measured with $\sim 0.02\,{\rm ns}$ precision, but suffers from a constant integer cycle bias that is initially unknown. Resolving this bias is known as integer ambiguity resolution. Combinations of the 4 observables (pseudorange and carrier phase on both frequencies) allow for robust detection of data outliers and cycle slips in the integer bias, and enable robust integer ambiguity resolution 3, 4. With integer ambiguities resolved, carrier phase data can then be modelled as pseudorange data, but with two orders of magnitude more precision.

Here we analyse data from the Jet Propulsion Laboratory (JPL) 5, in which the clock biases are given at $30\,{\rm s}$ intervals 10. These clock biases are generated using data from a global network of $\sim 100$ GPS geodetic receivers by a mature analysis system that is used routinely for purposes of centimetre-level satellite orbit determination, and millimetre-level positioning for scientific purposes, such as plate tectonics, Earth rotation, and geodynamics. The analysis standards that are applied are consistent with models and conventions specified by the International Earth Rotation and Reference Frames Service (IERS), in line with resolutions of the International Astronomical Union (IAU). These standards ensure best practices and consistency between various geodetic techniques, including very long baseline interferometry (VLBI), and satellite laser ranging (SLR), both of which have been instrumental in developing the IERS conventions. We note that similar types of analyses abiding by the IERS conventions are performed by several analysis centres around the world, and are routinely compared by the International GNSS Service. Results of these comparisons, along with the variety of scientific applications that depend on such data, provide an abundance of evidence that the clock biases are determined relative to each other with an accuracy $<0.1\,{\rm ns}$.

Analysis of GPS receiver data includes the effects of special and general relativity. All ideal clocks stationary on one of the Earth’s equipotential surfaces such as the geoid (“sea level”) have no relative phase drift. For purposes of discussion, let us define coordinate time as proper time on the geoid. (Actually a different coordinate time is chosen in geodesy, but that does not change the argument here.) For each satellite, the combination of velocity and gravitational potential causes satellite proper time to vary with respect to coordinate time. This is dominated by a positive drift at the level of $0.45$ parts per billion, which would be $\sim 13\,{\rm ns}$ per $30\,{\rm s}$ epoch. The hardware in the GPS satellite is designed to set the effective frequency of the atomic clocks such that the drift is zero on average, using the known semi-major axis of each satellite’s orbit. A residual effect results from the satellite moving in an ellipse, and thus with a time varying velocity and gravitational potential. The eccentricity of the GPS satellites is typically small $\lesssim 0.02$; nevertheless, the effect is a periodic variation of satellite proper time with an amplitude of $\sim 30\,{\rm ns}$ at $\sim 12\,{\rm hr}$ orbit period. This is modelled to first order by knowing the satellite’s position and velocity, and by considering the Earth as a spherically symmetric mass. All higher order effects effectively go into the definition of the satellite clock bias provided by JPL, consistent with international conventions. These residual effects are less than $1\,{\rm ns}$ over the $\sim 12\,{\rm hr}$ orbit period, and so are completely negligible over the time periods investigated here.

In JPL’s global GPS data analysis, the results of which are employed as input data to our dark matter (DM) search, the clock biases are estimated using GPS ionosphere-free carrier phase and pseudorange data combinations, as part of a multi-parameter least-squares estimation process using a square-root information filter. Other parameters that are estimated along with the clock biases include GPS orbits, station positions, Earth rotation, and atmospheric delay. Correlated errors between the clocks and other parameters are formally computed by the square-root information filter to be at the level $<0.1\,{\rm ns}$, consistent with the inferred level of accuracy.

In this process, there is effectively no restriction on the allowed behaviour of the clocks from one $30\,{\rm s}$ 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. Since only relative clock bias can be estimated, one clock is always held fixed in the estimation procedure. The choice of reference clock is irrelevant, in that our search algorithms look at differences in estimated clock biases.

In order to analyse the estimated clock biases, we further apply a first-order differencing procedure to the data, and define the pseudo-frequencies

[TABLE]

where $S^{(0)}(t_{k})$ is the original clock bias for the $t_{k}$ epoch (data point), and $t_{k}-t_{k-1}=30\,{\rm s}$ is the sampling interval. This is equivalent to taking a discrete derivative (up to a multiplicative factor), and acts to whiten the data (since the clock bias noise is dominated by random walk processes). This differencing procedure also removes any constant bias offsets, and transforms the linear frequency drifts into constant offsets in $S^{(1)}$. In practice, these residual offsets are small, and are removed by subtracting the mean of $S^{(1)}$ over the given day. The effect of this procedure is demonstrated for a single arbitrarily-chosen clock in Fig. S.1.

The clock biases from JPL 5 also come with a “formal error”. The formal error quantifies uncertainty in the determination of the clock bias, and does not directly incorporate the intrinsic clock noise. It varies slowly over time, and is typically on the order of $\sim 0.02-0.03\,{\rm ns}$; the formal error is dominated by the uncertainty in the satellite orbit determination. Only the most recent satellite clocks have observed temporal variations in $S^{(1)}$ that are at a similar level as the formal error, indicating that temporal variations in older clocks are actually due to clock behaviour rather than estimation error. A snapshot of the standard deviations for the $S^{(1)}$ clock data used in our analysis is presented in Table 1. A plot showing a few hours of $S^{(1)}$ data for an arbitrarily selected few satellite clocks is shown in Fig. S.2. This shows quite clearly how the more modern block IIF satellites 1 are substantially less noisy than the older-generations of satellites (by orders of magnitude).

Since the clocks are noisy, in order to discern the DM-induced signal from the intrinsic clock noise we rely on signals correlated across the entire GPS network, as discussed in the following sections. Moreover, there already exists more than 15 years of high-accuracy timing data that can be exploited in the search. This data stream is being routinely updated, and, in principal, timing data from any other atomic clocks that are synchronised with GPS can be included in the analysis (for a discussion on this point, see Ref. 46). Therefore, by analysing the new and existing GPS timing data, we can perform a sensitive search for transient DM signals, and if no signals are found we can then place limits on the DM–ordinary matter interaction strengths 6. The GPS network is particularly well suited for this type of search, for a number of reasons. The large number of clocks, and the very large ($\sim$ $50\,000\,{\rm km}$) diameter of the network increases both the chance of an interaction, and the sensitivity of the search, since we rely on correlated signal throughout the entire network. Similar arguments underpin motivations for searches using global network of atomic magnetometers 29, 30. Such magnetometry searches are sensitive to different types of interactions (as discussed below), and are therefore complementary to atomic clock searches.

S.2 Using GPS to search for topological defects

The null results from recent WIMP (weakly-interacting massive particle) direct-detection experiments have partly contributed to increased attention to ultralight field DM, such as axions 6, 7, 8, 9, 10. Ultralight fields may form stable topological defects (TDs), such as monopoles, strings, or domain walls, which can be a dominant or subdominant contribution to both DM and dark energy 11, 12, 25, 13, 14, 15, 16. The interactions of light fields with standard model (SM) fields can be written as a sum of effective interaction Lagrangians (so-called “portals”) 6

[TABLE]

where $\mathcal{L}^{\rm PS}$ represents the pseudoscalar (axionic) portal, and $\mathcal{L}^{\rm S^{1}}$ and $\mathcal{L}^{\rm S^{2}}$ are the linear and quadratic scalar portals, respectively. The linear and quadratic scalar portals, as will be demonstrated below, lead to changes in the effective values of certain fundamental constants and thus cause shifts in atomic transition frequencies, and so lend themselves well to searches based on the use of atomic clocks. The axionic portal leads to interactions that could cause spin-dependent shifts due to fictitious magnetic fields, and are thus are well suited to magnetometry searches 29, 30. We note that there are very stringent limits on the interaction strength for the linear scalar interaction coming from astrophysics and gravitational experiments (see, e.g., Refs. 41, 17). However, the constraints on the quadratic portal are substantially weaker 39. For this reason, we will concentrate on the quadratic scalar portal.

While we mainly address topological defect dark matter, and in particular refer the quadratic scalar coupling, it is important to note that our experiment is not limited in scope to this possibility. Any large (on laboratory scales), “clumpy” object (e.g., $Q$-balls 13, 18, 14, strings 19, solitons 16, 20, 21, and other stable objects 18, 19, 22) that interacts with standard model particles in such a way that leads to shifts in atomic transition frequencies is possible to detect using this scheme.

In the assumption of a quadratic scalar coupling between the standard model (SM) and DM fields, the interaction Lagrangian can be parameterized as 6

[TABLE]

[TABLE]

where $\varphi$ is the DM field, $m_{f}$ are the fermion masses, $\psi_{f}$ and $F_{\mu\nu}$ are the SM fermion fields and electromagnetic field tensor, respectively, and there is an implicit sum over $f$ that runs over all SM fermions. The constants $\Gamma$ quantify the strengths of the various DM–SM couplings. From a comparison with the conventional SM Lagrangian,

[TABLE]

it is seen that (to lowest order) the above Lagrangian (S.3) leads to the effective redefinition of certain dimensionless combinations of fundamental constants:

[TABLE]

[TABLE]

[TABLE]

where $\alpha$ and $\Lambda_{\rm QCD}$ are the nominal values of the fine structure constant and the QCD energy scale, respectively, and $m_{e,p}$ are the electron and proton masses. In our notation, the DM field has units of energy $(E)$, thereby $\Gamma_{X}$ is expressed in units of $E^{-2}$. Following Ref. 6, and to aid in the comparison with other works, we will present our results in terms of effective energy scales $\Lambda_{X}=1/\sqrt{|{\Gamma_{X}}|}$ (for $X=\alpha,\,m_{e},\,m_{p},\,m_{q}/\Lambda_{\rm QCD}$). Note that, by the nature of topological defects (TDs), the DM field $\varphi^{2}\to 0$ outside the defect, and $\varphi^{2}\to\varphi_{\rm max}^{2}$ inside the defect; as such the redefined coupling constants are only realised inside the defect. The field amplitude, $\varphi_{\rm max}$, can be linked to $\rho_{\rm inside}$, the energy density inside the defect, as

[TABLE]

where $d$ is the width of the defect, which is set by the Compton wavelength for the field 6:

[TABLE]

Note the ultra-light mass scale for the fields we are interested in: a roughly Earth-sized defect has a mass scale $m_{\varphi}\sim 10^{-14}{\,\rm eV}$ ($m_{\varphi}$ refers to the mass of the field particles, not the mass of the defect itself). We note here that there are other possibilities to search for ultralight DM fields via their non-gravitational interactions; see, e.g., Refs. 29, 30, 23, 32, 24, 25, 26, 23, 27, 28, 29, 31.

From Eqs. (S.4) – (S.6), we may relate the DM-induced frequency shift to a transient variation of fundamental constants. The fractional shift in the clock atom transition frequency, $\omega_{c}$, can be expressed as

[TABLE]

where $K_{X}$ are dimensionless sensitivity coefficients, and $X$ runs over fundamental constants from (S.4)–(S.5),

[TABLE]

and the overall effective coupling constant is defined as

[TABLE]

which depends on the specific clock transition through the sensitivity coefficients $K_{X}$. We also define $\Lambda_{\rm eff}=1/\sqrt{|\Gamma_{\rm eff}|}$, which have a meaning of effective energy scales. The dimensionless sensitivity coefficients $K_{X}$ are known from atomic and nuclear structure calculations 30. For example, considering only the variation in the fine structure constant, $\alpha$, and ignoring relativistic effects, for optical and microwave transitions, the clock frequencies scale as $\omega^{\rm opt}_{c}\propto\alpha^{2}$, and $\omega^{\rm mw}_{c}\propto\alpha^{4}$, respectively. Slight deviations from these dependencies occur due to relativistic corrections to the atomic structure 30, 31.

The clocks on board the GPS satellites work by disciplining the frequency of a quartz oscillator to the resonant frequency of the chosen atomic transition; the quartz oscillator drives the GPS L1 and L2 microwave signals 32. In order to make a meaningful analysis of the data we must analyse the comparison of the output frequencies to the frequency of another reference clock. For thin domain walls, the GPS clock and the reference clock are spatially separated, so only one is affected by the TD at a given time. (When both clocks are affected at the same time we would see no DM signal, as discussed in the next section.) Therefore, for the microwave frequency 87Rb and 133Cs clocks on board the GPS satellites, the sensitivity coefficients are 33, 34

[TABLE]

[TABLE]

respectively, where we have defined the short-hand notation $\Gamma_{e/p}\equiv 2\Gamma_{m_{e}}-\Gamma_{m_{p}}$, and $\Gamma_{q}\equiv\Gamma_{m_{q}/\Lambda_{\rm QCD}}$ (and similarly for $\Lambda_{e/p}$ and $\Lambda_{q}$). Note that the value of $K_{q}$ comes from a combination of a shift in the nuclear magnetic moment, $K_{\mu}$ and a shift in the nuclear size, $K_{\rm hq}$. For Cs, the contributions from these two effects are roughly equal in magnitude and opposite in sign 34 ($K_{q}^{\rm Cs}=K_{\mu}+K_{\rm hq}=0.009-0.007$), and $K_{q}$ in Cs is therefore sensitive to uncertainties in the nuclear calculations. For clocks based on optical transitions, there is only sensitivity to variations in the fine structure constant 30, 31. For clocks that compare the clock frequency to the resonant frequency of an optical cavity, there is an extra factor of $\alpha$ that comes from the variation in the length of the cavity, e.g.,

[TABLE]

as in Ref. 38. In general, the cavity contributes no extra sensitivity to variations in $m_{e}/m_{p}$ or $q/\Lambda_{\rm QCD}$ 28, 27.

The fact that there are a number of different atomic clock types employed within the GPS network means we have sensitivity to several different combinations of the parameters (S.9). The existing limits come from constraints on supernova emission 39: ${\Lambda_{m_{e},\alpha}\gtrsim 3\,{\rm TeV}}$, and $\Lambda_{m_{p}}\gtrsim 15\,{\rm TeV}$, with no existing limit on $\Lambda_{q}$.

S.3 Domain wall signal

One specific example of a TD is a domain wall, a quasi-2D structure that can be characterised by a width, $d$, and DM field amplitude 6, $\varphi_{\rm max}$ (see Fig. 1 of main text). We assume a Gaussian distribution of the field across the wall. As a wall passes a clock, it causes a shift in the clock frequency $\omega_{c}\to\omega_{c}+\delta\omega(t)$, where $\delta\omega(t)$ is proportional to the square of the field value and rapidly goes to zero outside of the wall. Then, the time-difference due to a transient frequency shift with respect to an unperturbed clock is $\Delta t=\int_{-\infty}^{t}\frac{\delta\omega(t^{\prime})}{\omega_{c}}dt^{\prime}$. The accumulated time-reading bias of clock $i$ (measured against some reference clock, $R$) at time $t$ caused by a Gaussian-profile domain-wall that crosses the clock at time $t^{\times}_{i}$ and the reference clock at $t^{\times}_{R}$, is

[TABLE]

[TABLE]

where the (single-clock) crossing duration is defined as

[TABLE]

which is the time scale during which the clock is inside the wall, and $v_{\perp}$ is the component of the wall’s velocity that is perpendicular to the wall.

As mentioned above, to analyse the data, we first perform a single-differencing procedure (S.1). We perform the same differencing procedure for the DM signals. The correspondence between the $S^{(0)}$ and $S^{(1)}$ signals is shown in Fig. S.3.

For thin walls (that is, for walls thin enough such that the interaction time is smaller than the data acquisition interval, $\tau<30\,{\rm s}$), the $S^{(1)}$ spike amplitude would see in the first-order differenced data at the time of the wall crossing (for a system of identical clocks) is

[TABLE]

as shown in Fig. S.3. Note that in the “thin” wall case, the assumed Gaussian profile is unimportant; other profiles (e.g., hard spheres, flat profiles) give similar results. These non-Gaussian profiles can be incorporated through the form factors in Eq. (S.17) arising from the integrals in Eq. (S.15). Our following analysis holds equally in these cases, the difference being only the ratio of the specific form factors. For example, for a flat (hard-edge) profile, $\Gamma_{\rm eff}\to\Gamma_{\rm eff}/\sqrt{\pi}$ in Eq. (S.17).

In order to determine the expected average $S^{(1)}$, we use the average value for the single-clock crossing duration, $\tau_{\rm avg}$, which depends on the object width $d$, orientation, and velocity distribution, as will be discussed in Section S.5. Further, in the assumption that the TDs saturate the DM energy density in the galaxy, one can link $\varphi_{\rm max}$ and $d$ to $\rho_{\rm DM}$:

[TABLE]

where $\mathcal{T}$ is the average time between events (i.e., encounters between the GPS constellation and DM objects). Therefore, the average DM signal amplitude is:

[TABLE]

This equation translates into constraints on $\Gamma_{\rm eff}$, Eq. (5) of the main text.

S.4 Search

Here, we perform a search with the aim of ruling out or detecting large ($\gtrsim 5\sigma$) events (here, $\sigma$ is the standard deviation of $S^{(1)}$ clock data, which is about $0.1\,{\rm ns}$ in the worst-case for the Rb clock solutions when using another Rb clock as a reference, see Table 1). Consider a time-period window, denoted $J_{w}$. If a thin-wall TD-DM object passes through the GPS constellation with a speed

[TABLE]

where $D_{\rm GPS}\approx 50\,000\,{\rm km}$ is the GPS orbital diameter, we can expect that all of the clocks will be affected by the object within this window. However, any clock that is swept within the same 30 s sampling interval as the reference clock will not show any spike in the $S^{(1)}$ data due to clock degeneracy (see Fig. S.3, or satellites 15 and 16 in Fig. 3(a) of the main text). For a DM object with relative velocity $v<700\,{\rm km\,s}^{-1},$ we can safely expect at least 60% of the clocks to show an $S^{(1)}$ spike (this number comes from the orbital geometry of the GPS constellation, noting that by design the satellites are relatively evenly distributed over the sky). For more typical velocities, $v\sim 100-300\,{\rm km\,s}^{-1}$, more than 80% of the satellites are expected to show the $S^{(1)}$ spikes.

To ensure that the potential DM-induced frequency excursions for each clock are of the same magnitude, we only consider sub-networks of identical clocks. For example, for the Rb sub-network, we choose the most recently launched Rb satellite to be the reference clock (as a rule such clocks are the least noisy), and compare all the other Rb clocks to this reference, and do not include the Cs or H-maser clocks in the analysis. To unfold limits on the various coupling constants in Eqs. (S.12) and (S.13), we only consider the Rb and Cs sub-networks, and supplement them with the Sr limits on $\Gamma_{\alpha}$ from Ref. 38. We do not include the ground clocks in our current analysis.

To search for domain wall signals, we analysed the $S^{(1)}$ GPS data streams in two stages. The first stage involved stepping through all the data, one epoch at a time, to identify regions of interest that could potentially be consistent with thin-wall crossings. We call these regions “potential events”. The second stage investigates these regions using a more detailed approach, to determine if we can upgrade any potential events to “candidate events”.

At the first stage, we scanned all the data from May 2000 to October 2016 looking for the most general patterns associated with a domain wall crossing, without taking into account the order in which the satellites were swept. We required at least 60% of the clocks to experience a frequency excursion at the same epoch, as discussed above. This procedure identifies the epoch when the wall crossed the reference clock (vertical blue line in Fig. 3(a) of the main text). If such a frequency excursion is found at epoch $t_{0}$, we form all the possible windows, $J_{w}$, around $t_{0}$, in the range $3\leq J_{w}/30\,{\rm s}\leq 500$. The chosen maximum $J_{w}$ determines the minimum DM velocity we can detect (S.20). Our choice of a maximum window of $500\,{\rm epochs}$, corresponding to sweep durations through the GPS constellation of up to 15,000 s, ensures sensitivity to walls moving as slowly as $\sim 4\,{\rm km\,s}^{-1}$. Less than 0.1% of DM domain walls that cross the GPS network are expected to travel with perpendicular velocities slower than this value.

We define a “potential event” as any time where at least the given fraction of the clocks also experience a frequency excursion of sign opposite to the reference clock excursion anywhere within that window (red tiles in Fig. 3 of main text). Since we consider only the sub-networks of identical clocks (e.g., we only consider the Rb clocks), if an event occurs within the window all the DM-induced frequency excursions must be of the same magnitude (besides the clock noise). Therefore, we only search for frequency excursions within a range of $S^{(1)}$ values, $S^{(1)}_{\rm cut}\to S^{(1)}_{\rm cut}\pm dS$, where $dS=x\sigma$, with $\sigma$ being the calculated standard deviation for that given GPS clock/reference clock combination (see Table 1); for 90% confidence level (C.L.), $x=1.645$. (By using the calculated standard deviation here, we account for both the formal error (see section S.1) and the clock noise directly.) We count the total number of potential events for each window-size and $S^{(1)}_{\rm cut}$ combination.

This method also allows us to place more stringent limits on events for which the average time between events, $\mathcal{T}$, is less than the observation time, $\mathcal{T}_{\rm obs}$. For example, if for a given $S^{(1)}_{\rm cut}$, we saw 10 potential events in a 10 year observation time, we can place a limit at this level for $\mathcal{T}\simeq 1\,{\rm yr}$. Further, we only include days in our analysis for which there is data available for at least 7 clocks; this effectively reduced our observation time, especially for the Cs sub-network, which employs only a small number of clocks in recent times. Note, the choice of 7 simultaneously operational clocks here is conservative; injecting fake DM events demonstrates that this search technique would find practically every event (that otherwise falls within our current sensitivity) with at least 5 available clocks.

The tiled representation of the GPS data stream depends on the chosen signal cut-off $S^{(1)}_{\mathrm{cut}}$ (see Fig. S.2 and Fig. 3 of the main text). We systematically decreased the cut-off values and repeated the above analysis. Wherever the counted number of “potential events” goes to zero gives the overall limit for the total observation time, which we denote $S^{(1)}_{\rm thresh}$. For example, for a window size encompassing sweeps at $v\sim 300\,{\rm km\,s}^{-1}$ ($J_{w}=170\,{\rm s}$), we can exclude events in the Rb network above the $S^{(1)}_{\rm thresh}=0.48\,{\rm ns}$ level (90% C.L.), for a total observation time of $\mathcal{T}_{\rm obs}=16.5\,{\rm yrs}$. Below the thresholds, a number of potential events were identified. An illistrative subset of the results of this analysis is presented in Table 2 for the Rb sub-network, and in Table 3 for the Cs sub-network.

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)}_{\mathrm{thresh}}$, by matching the data streams against the expected patterns; one such example is shown in Fig. 3(a). 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, and we only considered time window lengths of $J_{w}<50\,{\rm epochs}$ ($v_{\rm min}\approx 30\,{\rm km\,s}^{-1}$). As a result of this pattern matching, we did not find any events that were consistent with domain wall DM, thus we have found no candidate events at our current sensitivity. Analyzing numerous potential events well below $S^{(1)}_{\rm thresh}$ (see Tables 2 and 3) has proven to be substantially more computationally demanding, and is beyond the scope of the current work.

While the detailed second stage analysis have improved (lowered) the values of $S^{(1)}_{\rm thresh}$ by rejecting a few hundred potential events, the improvement was minor. Therefore, to be conservative, we have used the $S^{(1)}_{\rm thresh}$ results from the first stage of the analysis to place our constraints on DM walls.

S.5 Crossing duration distribution and domain wall width sensitivity

Before placing limits on $\Gamma_{\rm eff}$, we need to account for the fact that we do not have equal sensitivity to each domain wall width, $d$, or equivalently crossing durations, S.16. This is due in part to aspects of the clock hardware and operation, the time-resolution (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 ${\tau_{\rm min}<\tau<\tau_{\rm max},}$ where $\tau_{\rm min(max)}$ is the minimum (maximum) crossing duration that we are sensitive to.

There are two factors that determine $\tau_{\rm min}$. The first is the so-called “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 known to us at a given epoch, however, it is known 32, 40 to be within $0.01$ and $0.1\,{\rm s}$. As such, we consider the best and worse case scenarios:

[TABLE]

Note, that below $\tau^{\rm Servo}$ we still have sensitivity to DM events, however, the sensitivity in this region is determined by response of the quartz oscillator to the temporary variation in fundamental constants. However, the resulting limits for crossing durations shorter than the servo-loop time are generally weaker than the existing astrophysics limits (see section S.6), so we consider this no further.

The second condition that affects $\tau_{\rm 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 Fig. S.3 and satellites 15 and 16 in Fig. 3(a) of the main text). For less than 60% of the clocks to experience the jump, the (thin-wall) DM object would have to be travelling at over $700\,{\rm 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, less than $0.1\%$ of events are expected to have $v_{\perp}>700\,{\rm km\,s}^{-1}$.) This velocity corresponds to a crossing duration for the entire network of $\sim 70\,{\rm s}$. Transforming this to crossing duration for a single clock, $\tau$, amounts to multiplying by the ratio $d/(D_{\rm GPS})$:

[TABLE]

For the thickest walls we consider ($\sim 10^{4}\,{\rm km}$), this leads to $\tau_{\rm min}^{\rm Degen.}$ of 14 s.

Combining the servo-loop and degeneracy considerations, we arrive at the expression

[TABLE]

For walls thicker than $d\sim 100\,{\rm km}$, $\tau_{\rm min}$ is determined by the ${d}/(D_{\rm GPS})\,70\,{\rm s}$ term.

As to the maximum crossing duration, there are also two factors that affect $\tau_{\rm max}.$ The first, is that 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”:

[TABLE]

If a wall takes longer than 30 s to pass by a clock, the simple single-data–point signals shown in Fig. S.3 would become more complicated, and would require a more-detailed pattern-matching technique. Second, we only consider time “windows” of a certain size in our analysis (see section S.4). If a wall moves so slowly that it does not sweep all the clocks within this window, the event would be missed:

[TABLE]

Therefore, the overall expression for $\tau_{\rm max}$ is:

[TABLE]

Making $J_{w}$ large, however, also tends to increase $S^{(1)}_{\rm thresh}$ (since there is a higher chance that a large window will satisfy the condition for a “potential event”, see section S.4). By performing the analysis for multiple values for $J_{w}$, we can probe the largest portion of the parameter space. In this work, we consider windows of $J_{w}$ up to 500 epochs (15000 s), which corresponds to a minimum velocity of $\sim 4\,{\rm km\,s}^{-1}$, which is roughly the orbit speed of the satellites. This has a negligible effect on our sensitivity, since less than $0.1\%$ of walls are expected to have $v_{\perp}<4\,{\rm km\,s}^{-1}$.

Assuming the standard halo model, the relative scalar velocity distribution of DM objects that cross the GPS network is quasi-Maxwellian

[TABLE]

where $v_{c}=220\,{\rm km\,s}^{-1}$ is the Sun’s velocity in the halo frame, and $C$ is a normalisation constant. The form of Eq. (S.27) is a consequence of the motion of the reference frame. However, the distribution of interest for domain walls is the perpendicular velocity distribution for walls that cross the network

[TABLE]

where $v_{\perp}$ is the component of the wall’s velocity that is perpendicular to the wall, and $C^{\prime}$ is a normalisation constant. Note that this is not the distribution of perpendicular velocities in the galaxy – instead, it is the distribution of perpendicular velocities that are expected to cross paths with the GPS constellation (walls with velocities close to parallel to face of the wall are less likely to encounter the GPS satellites, and objects with higher velocities more likely to). Here, we assumed 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, we also note that our results are fairly independent of the exact form of the velocity distribution. Even if the actual wall velocity distribution is somewhat different, the qualitative feature of a TD “wind” is not expected to change. For example, if a larger proportion of the DM objects move slower than typical galactic velocities, almost nothing changes, since then the relative velocity is essentially given by the Earth velocity in the galactic frame.

Now, define a function $f_{\tau}(d,\tau)$, such that the integral $\int_{\tau_{a}}^{\tau_{b}}f_{\tau}(d,\tau)d\tau$ gives the fraction of events due to walls of width $d$ that have crossing durations between $\tau_{a}$ and $\tau_{b}$. Note, this function must have the following normalization:

[TABLE]

for all $d$, and is given by

[TABLE]

Plots of the velocity and crossing-time distributions are given in Fig. S.4. Then, our sensitivity at a particular wall width is

[TABLE]

Plots of the sensitivity function for a few various cases of parameters are presented in Fig. S.5.

S.6 Placing limits

The most stringent limit on the effective energy scale $\Lambda_{\rm eff}=1/\sqrt{|{\Gamma_{\rm eff}}|}$ is set by assuming one event could have occurred in the observation time, which caused a frequency excursion in the timing data that was equal in magnitude to $S^{(1)}_{\rm thresh}$. In this case:

[TABLE]

[TABLE]

To determine the confidence level for our limits, we must also factor in the uncertainties. The uncertainty from the clock noise is directly built into the $S^{(1)}$ values in Tables 2 and 3. For our current analysis, the $S^{(1)}$ limits are substantially greater than the actual clock-noise, on the order of $5-20\sigma$. The dominating uncertainty in our limits comes from equating the observation time with the average time between events.

To set the limits in the region where $\mathcal{T}\simeq\mathcal{T}_{\rm obs}$, we assume that the frequency of DM–GPS encounters is roughly Poissonian. Suppose we expect to see on average $\lambda$ events in the observation time $\mathcal{T}_{\rm obs}$. The probability for observing at least one event in the time period $\mathcal{T}_{\rm obs}$ is given by

[TABLE]

where $P_{k}(\lambda)={\lambda^{k}e^{-\lambda}}/{k!}$ is the Poisson distribution. For example, to place 90% confidence level limits we require that $P_{k\geq 1}(\lambda)=0.9.$ In this case, solving (S.33) gives $\lambda=2.3$. Therefore, the maximum $\mathcal{T}$ for which we can place 90% C.L. limits is given by

[TABLE]

where for a 90% C.L. limit, $\lambda=2.3$ and $\mathcal{T}_{\rm max}\simeq 7\,{\rm yr}$.

For the region of parameter space where $\mathcal{T}\ll\mathcal{T}_{\rm obs}$, our sensitivity is lower. This is because for small values of the average time between events, the total number density of DM objects must be correspondingly high, leading to smaller DM field values inside a defect. In the assumption that such objects constitute a significant fraction of the DM, their total energy must still add up to the total observed local DM density 35 of ${\rho_{\rm DM}\approx 0.4\,{\rm GeV\,cm}^{-3}}$:

[TABLE]

Therefore, the energy density of each object must be smaller, and as such, the resulting signal would be smaller.

For the applicable region ($\mathcal{T}<\mathcal{T}_{\rm obs}/\lambda$), we combine the $s(d)$ function with (S.19) to get the final constraints:

[TABLE]

Similarly, in the case when one of the specific couplings, $\Gamma_{X}$, dominates over the other coupling strengths in the linear combination in Eqs. (S.12) and (S.13), we have

[TABLE]

Note that the only $d$-dependence in $s(d)$ comes from the integration limits (via $\tau_{\rm min/max}$). Expressing in convenient units:

[TABLE]

Using the relation $d=\hbar/m_{\varphi}c$, one can rewrite the above limit in terms of the field mass $m_{\varphi}$ with the substitution ${(d/{\rm km})\approx{2\times 10^{-10}}{}\,({\rm eV}/m_{\varphi}c^{2}).}$

The resulting 90% C.L. limits on $\Lambda_{\rm eff}$ from combining the $S^{(1)}$ limits from Tables 2 and 3 with Eq. (S.36), for the case when $\mathcal{T}=\mathcal{T}_{\rm obs}/\lambda=7\,{\rm yr}$, are shown in Fig. S.6. For smaller values of $\mathcal{T}$, the limits scale as $\sqrt{\mathcal{T}}$. See also Fig. 4 of the main text, which shows a contour plot of the constraints as a function of $d$ and $\mathcal{T}$.

Recently, the group from Toruń 38 used an optical Sr clock to place limits on the coupling of topological defect DM to atoms. Since this group employed an optical transition in Sr, this experiment is only sensitive to the variation in the fine structure constant (i.e., $\Lambda_{\alpha}$), see Eq. (S.14). Their data covers a period of $\mathcal{T}_{\rm obs}=45700\,{\rm s}\approx 13\,{\rm hrs}$. We combine their derived limits on $\Lambda_{\alpha}$ with our results in the exclusion plots. We note that although tight constraints can be placed using this method38, in order to distinguish a true DM-induced transient event from other external sources (such as electromagnetic interference, or direct physical disturbance of the clocks), a global network is prerequisite . One of the main advantages of our method of employing the GPS constellation is the reliance on such a global network. Another advantage of our approach is the availability of archival data for at least the past 16 years, giving us sensitivity to the region of the parameter space with $\mathcal{T}\gtrsim 1-10\,{\rm yrs}$, which is currently inaccessible by other methods.

If we make assumptions about the relative strengths of the couplings in Eq. (S.9), we can place limits on individual energy scales (S.37). For example, in the assumption that $\Lambda_{\alpha}\ll\Lambda_{e/p},\Lambda_{q}$, we can place limits directly on $\Lambda_{\alpha}$ (and likewise for $\Lambda_{e/p}$ and $\Lambda_{q}$). These resulting limits for $\Lambda_{\alpha}$ (and comparison with the results of Ref. 38) are shown in Fig. 5 of the main text. Note that due to differences in the experimental technique, the optical Sr limits scale with the wall width as $d^{3/4}$, and the sensitivity of the approach of Ref. 38 reduces sharply for widths greater than the Earth radius due to a frequency cut-off used in the analysis 38. In contrast, the GPS limits from this work scale linearly with $d$ (S.37). Both our limits and those of Ref. 38 scale as $\sqrt{\mathcal{T}}$, but have sharp cut-offs above the observation time ($\mathcal{T_{\rm obs}}\sim 16\,{\rm yr}$ for our work, and $\mathcal{T_{\rm obs}}\sim 10^{-3}\,{\rm yr}$ for Ref. 38). Limits on $\Lambda_{e/p}$ and $\Lambda_{q}$ (assuming these respective couplings dominate) are shown in Fig. S.8 (these couplings are unconstrained by the Sr experiment38).

From the three independent limits (one from each of the Rb and Cs sub-networks determined in this work, and one from the optical Sr clock used in Ref. 38), we can derive independent limits on $\Lambda_{e/p}$ and $\Lambda_{q}$, without having to make assumptions about the relative strengths of the individual couplings. This is possible because the three limits (from Cs, Rb, and Sr) depends on a different linear combination of the three available couplings, see Eqs. (S.12)–(S.14). A plot showing the combined allowed region for the coupling parameters $\Gamma_{\alpha}$, $\Gamma_{e/p}$, and $\Gamma_{q}$ is presented in Fig. S.7 (note $\Lambda_{X}=1/\sqrt{|\Gamma_{X}|}$). The resulting limits on $\Lambda_{e/p}$ and $\Lambda_{q}$ are presented in Figs. S.9 and S.10, respectively.

• 1 GPS.gov, http://www.gps.gov/systems/gps/space/.
• 2 International GNSS Service (IGS), http://www.igs.org/network.
• 3 G. Blewitt, J. Geophys. Res. Solid Earth 94, 10187 (1989).
• 4 G. Blewitt, Geophys. Res. Lett. 17, 199 (1990).
• 5 The Jet Propulsion Laboratory (JPL), ftp://sideshow.jpl.nasa.gov/pub/jpligsac/.
• 6 R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
• 7 R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
• 8 M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
• 9 P. Sikivie, Phys. Rev. Lett. 51, 1415 (1983).
• 10 J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
• 11 P. Sikivie, Phys. Rev. Lett. 48, 1156 (1982).
• 12 W. H. Press, B. S. Ryden, and D. N. Spergel, Astrophys. J. 347, 590 (1989).
• 13 R. A. Battye, M. Bucher, and D. N. Spergel (1999), arXiv:astro-ph/9908047 .
• 14 R. Durrer, M. Kunz, and A. Melchiorri, Phys. Rep. 364, 1 (2002).
• 15 A. Friedland, H. Murayama, and M. Perelstein, Phys. Rev. D 67, 043519 (2003).
• 16 P. P. Avelino, C. J. A. P. Martins, J. Menezes, R. Menezes, and J. C. R. E. Oliveira, Phys. Rev. D 78, 103508 (2008).
• 17 B. Bertotti, L. Iess,  and P. Tortora, Nature 425, 374 (2003).
• 18 G. Dvali, A. Kusenko, and M. Shaposhnikov, Phys. Lett. B 417, 99 (1998).
• 19 C. J. Hogan and M. J. Rees, Nature 311, 109 (1984).
• 20 A. Kusenko, Phys. Lett. B 405, 108 (1997).
• 21 T. D. Lee and Y. Pang, Phys. Rep. 221, 251 (1992).
• 22 P. Jetzer, Phys. Rep. 220, 163 (1992).
• 23 D. Budker, P. W. Graham, M. P. Ledbetter, S. Rajendran, and A. O. Sushkov, Phys. Rev. X 4, 021030 (2014).
• 24 A. Arvanitaki, J. Huang, and K. Van Tilburg, Phys. Rev. D 91, 015015 (2015).
• 25 A. Arvanitaki, S. Dimopoulos, and K. Van Tilburg, Phys. Rev. Lett. 116, 031102 (2016a).
• 26 Y. V. Stadnik and V. V. Flambaum, Phys. Rev. Lett. 115, 201301 (2015a).
• 27 Y. V. Stadnik and V. V. Flambaum, Phys. Rev. A 93, 063630 (2016).
• 28 Y. V. Stadnik and V. V. Flambaum, Phys. Rev. Lett. 114, 161301 (2015b).
• 29 A. Arvanitaki, P. W. Graham, J. M. Hogan, S. Rajendran, and K. Van Tilburg,  (2016b), arXiv:1606.04541.
• 30 V. A. Dzuba, V. V. Flambaum, and M. V. Marchenko, Phys. Rev. A 68, 022506 (2003).
• 31 E. J. Angstmann, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. A 70, 014102 (2004).
• 32 R. T. Dupuis, T. J. Lynch, and J. R. Vaccaro, in 2008 IEEE Int. Freq. Control Symp. (IEEE, 2008) pp. 655–660.
• 33 V. V. Flambaum and A. F. Tedesco, Phys. Rev. C 73, 055501 (2006).
• 34 T. H. Dinh, A. Dunning, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. A 79, 054102 (2009).
• 35 F. Nesti and P. Salucci, J. Cosmol. Astropart. Phys. 2013, 016 (2013).