New bounds on dark matter coupling from a global network of optical atomic clocks

See allHide authors and affiliations

Science Advances  07 Dec 2018:
Vol. 4, no. 12, eaau4869
DOI: 10.1126/sciadv.aau4869


We report on the first Earth-scale quantum sensor network based on optical atomic clocks aimed at dark matter (DM) detection. Exploiting differences in the susceptibilities to the fine-structure constant of essential parts of an optical atomic clock, i.e., the cold atoms and the optical reference cavity, we can perform sensitive searches for DM signatures without the need for real-time comparisons of the clocks. We report a two orders of magnitude improvement in constraints on transient variations of the fine-structure constant, which considerably improves the detection limit for the standard model (SM)–DM coupling. We use Yb and Sr optical atomic clocks at four laboratories on three continents to search for both topological defect and massive scalar field candidates. No signal consistent with a DM coupling is identified, leading to considerably improved constraints on the DM-SM couplings.


Evidence for the existence of dark matter (DM) comes entirely from astrophysical observations at large scales (galactic to cosmological). The nature of DM composition, however, will be known only after positive detection of the DM candidates. The viable cold DM candidates, such as axions (13), weakly interacting massive particles (WIMPs) (3), and super-WIMPs (4), require the existence of fields that can be coupled to standard model (SM) fields. Therefore, existing experiments focus on searches for these couplings, although no experimental data have provided a positive detection. In the Large Underground Xenon (LUX) experiment, for example, which studied potential coupling between WIMPs and nucleons, constraints on the scattering cross section per nucleon were reported below 10−45 cm2 (5). New stringent upper limits were recently reported in the 8 to 125 keV/c2 mass range that excluded couplings to electrons with coupling constants of gae > 3 × 10−13 for pseudoscalar and α′ > 2 × 10−28 × α (α is the fine-structure constant, and α′ is its vector boson equivalent) for vector super-WIMPs, respectively, by the XENON100 experiment (6).

Lack of any detected DM in the form of particles yielded alternative theories. For instance, DM can be at least partially explained by the introduction of ultralight (sub–electron volt) massive scalar fields to the general relativistic tensor-scalar theories (7). Moreover, in the early Universe, phase transitions may have produced topological defects in the scalar fields with suitable self-potential (8). If the scalar particles are sufficiently light and weakly interacting, then the scalar DM, in the form of the coherently oscillating (with frequency related to its mass) classical field, could survive to the present day. In addition, depending on the nature of the symmetry breaking during the phase transition, several types of stable topological defects can be formed, including domain walls, strings (vortices), point-like monopoles, or stable hybrid defects such as monopoles connected by strings (9). The standard halo model predicts a quasi-Maxwellian distribution of DM in the galactic frame. The coupling of the scalar field ϕ to the SM fermion and electromagnetic fields may be represented by linear or quadratic interactions (10, 11) and could be detected by observing possible variations of the effective fine-structure constant (8, 1216). The coherently oscillating classical fields may induce linear drifts in time and oscillating variations of the effective fine-structure constant, while topological defects may induce transient variations in time of the effective fine-structure constant. Studies of quasar absorption spectra from independent data observed by the VLT (very large telescope) and Keck telescopes indicate that the fine-structure constant may vary on a cosmological space-time scale (17, 18), which would result in a drift of the fine-structure constant in Earth’s frame of reference at the level of δα/α ~ 10−19 per year. In addition, a topological defect passing through a pulsar may trigger a pulsar glitch, and pulsar glitches have already been observed, although the origin of these glitches is still discussed and in need of further numerical and observational study (17, 19).

Exceptional accuracy and stability of optical atomic clocks and ultranarrow lasers open the door for new tests of fundamental physics that can reach beyond the SM. Optical atomic clocks were used, for example, to constrain the coupling of fundamental constants such as the fine-structure constant α, electron-proton mass ratio μ, and light quark mass to gravity (20). Optical atomic clocks connected by phase-compensated optical fiber links improved a constraint on the Robertson-Mansouri-Sexl parameter quantifying a violation of time dilation (21). Stability of optical clocks and ultranarrow lasers in our network allows for determination of a new constraint on the transient perturbation in the fine-structure constant |δα/α| < 1.6 × 10−16, which is an improvement by two orders of magnitude.

Recently, we have shown that a single optical atomic clock can be used as a detector for DM in the form of stable topological defects (12). This is possible by taking advantage of the differences in the susceptibilities to the fine-structure constant between essential parts of any optical atomic clock, i.e., the atoms and the cavity used to prestabilize the clock laser. Here, we realize the first Earth-scale quantum sensor network based on optical atomic clocks. In our approach, measurements by the distant clocks do not need to be linked in real time. In analogy to the standard radio-astronomical technique, very-long-baseline interferometry (VLBI), clock measurements can be locally recorded (with time stamps that are accurate to the level of 1 ms) and cross-correlated later. Among different types of topological defects, our network is best suited for detecting the domain-wall type, with transverse dimensions that exceed the span of our network regardless of the defect thickness. A concept of archiving data from various types of precise measurements was also proposed in (22). In this letter, we report the results of measurements of Yb and Sr optical clocks at four laboratories on three continents to search for topological defect and massive scalar field candidates.


Our network is composed of optical lattice atomic clocks located at the National Institute of Standards and Technology (NIST), Boulder, CO, USA (23, 24), at LNE-SYRTE (Laboratoire National de Métrologie et d’Essais, Systèmes de Références Temps-Espace), Paris, France (25, 26), at KL FAMO (Krajowe Laboratorium Fizyki Atomowej, Molekularnej i Optycznej), Torun, Poland (27, 28), and at the National Institute of Information and Communications Technology (NICT), Tokyo, Japan (see Fig. 1) (29, 30). Each clock uses an optical cavity–stabilized laser that is frequency tuned to resonantly interrogate the atomic clock transition of cold atoms trapped in an optical lattice. The clock records the local difference ri(t) between the atomic transition frequency, vatom, and the cavity-stabilized laser, vcavity, where i = 1…4 corresponds to the four clocks in the network. Any variation in the electromagnetic fine-structure constant, α, will be imprinted in this difference because of different susceptibilities of these two frequencies to variations in α (vatom/vcavity ∝ α) (12, 31). In particular, we may expect variations in α that can be expressed as δα/α = (ϕ/Λγ,n)n, where ϕ is the DM field and Λγ,n is the energy scale (which inversely parametrizes the strength of the DM-SM coupling) of the nth portal (32).

Fig. 1 Global sensor network.

The participating optical lattice atomic clocks reside at NIST, Boulder, CO, USA (23, 24), at LNE-SYRTE, Paris, France (25, 26), at KL FAMO, Torun, Poland (27, 28), and at NICT, Tokyo, Japan (29, 30).

The general potential of self-interaction of the scalar field can yield an evolution of the scalar field in the form of a coherently oscillating classical field (13), as well as spontaneous symmetry breaking and, as a consequence, to the formation of topological solitons, also called topological defects (9). Both oscillating field and topological defects were recently proposed as possible forms of the DM (7, 8) and may both induce a variation of the fine-structure constant in the form of oscillations or transient changes, respectively. The linear coupling n = 1 is usually considered in the oscillating massive scalar field DM studies (7, 13, 33), while quadratic coupling n = 2 is usually considered for the topological DM (8, 12, 34). One can also translate limits from linear to quadratic couplings (11).

In the case of topological defects, the transient effect can be detected by cross-correlating the pairs of ri(t) recorded in our geographically dispersed clocks. To allow for the comparison between clocks based on different atomic species and different transitions, we define the normalized signals as Embedded Image, where Embedded Image is the sensitivity coefficient [equal to one in the nonrelativistic case (31)] and Embedded Image is the frequency of the clock transition. We applied a smooth high-pass FIR (finite impulse response) filter with frequency cutoff between 0.005 and 0.027 Hz (depending on the size of the sought defect) sampling both past and future points weighted with a Gaussian profile to all Embedded Image to eliminate the low-frequency drifts that originate from the long-term cavities’ instabilities. We model the transient effect by a square perturbation of the fine-structure constant with size d and with an amplitude assumed to be identical for the distant clocks of our network. We take into account rotation and revolution of Earth around the Sun as well as revolution of our Solar System around the center of our galaxy by assuming the topological defects’ speed relative to Earth v = 300 km s−1 (with dispersion of ≈300 km s−1) and considering all possible direction of the topological defects’ relative velocity. We validated this model by showing that, for the most constrained limit, the results calculated with a square profile differ only by a few percent from the results calculated with a simple model of a stable plain domain wall in the xy plane with transverse profile of a scalar field ϕ(z) ~ tanh(z/L), where Ld/2 is the thickness of the wall orthogonal to the direction z of defect propagation (9, 35). Our analysis was performed for the defect size, d, between 3 × 102 and 3 × 104 km. We cross-correlate the pairs of Embedded Image and fit the resulting data with the expected DM signature (see Materials and methods). The left panel in Fig. 2 shows two Embedded Image from the NIST and LNE-SYRTE laboratories. The right panel shows their cross correlation and the fit of the expected shape (triangle) for the 30-s-long event. The inset shows a fragment of the fitted function with triangular shape.

Fig. 2 Cross correlation of two Embedded Image.

Left: Overlapping normalized signals, Embedded Image, from the NIST and LNE-SYRTE laboratories (pink and blue lines, respectively). Right: Cross correlation, A(τ), of the two Embedded Image (brown line) and the fit of the expected cross-correlation function assuming a square perturbation with duration of 30 s (green line). τ is the cross-correlation displacement. The inset depicts the magnification of the region where the fitted triangular function can be seen.

The constraint on the transient variation of the fine-structure constant can be expressed asEmbedded Image(1)where A0 is the fitted amplitude of the cross-correlation peak and ηT is the ratio of the overall defect duration to the duration of the cross-correlated data. The top right panel of Fig. 3 shows the bound on the transient variation in the fine-structure constants, δα/α, derived from our measurement. The black line comes from the direct fit of the expected DM signature to the cross correlation. To give a proper statistical interpretation of our result and to have an ability to distinguish the noise from any exotic physical signal, we also calculated the 5 and 95% confidence levels (CLs). The constraints reported here are the 95% CLs (see the red lines in Fig. 3). A 95% CL indicates that when the signals are assumed to contain only noise, the probability that our determination of α variation is larger than our constraint is 100% − 95% = 5%. To determine the 95% CL, we artificially shift one of the signals by a large period, ensuring that no physical correlation between the two channels should be expected. We repeat this procedure 100 times for different shifts to generate a sufficiently large null signal sample for determination of the CL. For the particular case of a defect size comparable to the size of Earth (d = 10,000 km), we determine a new constraint on the transient perturbation in the fine-structure constantEmbedded Image(2)which is almost two orders of magnitude better than the best limit published so far (12), represented by the orange lines in Fig. 3. For other values of d, the improvement is even larger and approaches three orders of magnitude for d = 900 km. For comparison, we also plotted the actual and possible limits derived from the GPS constellation’s on-board clocks (solid and dashed gray lines, respectively) (16).

Fig. 3 Constraints on the energy scale Λγ,2 and the fine-structure variation δα/α.

Top right: Constraints (95% CL) on the fine-structure variation δα/α (red line). Top left and bottom: Constraints (95% CLs) on the energy scale Λγ,2 under the assumption of the local mean cold DM energy density ρDM = 0.468 GeV cm−3 (red line) (3, 36). The black and green lines are the direct fit and the 5% CL, respectively. The orange lines represent the previous best constraints measured with two colocated 88Sr clocks (12), and the gray solid and dashed lines represent the actual and possible constraints, respectively, estimated from the Global Positioning System (GPS) constellation’s on-board clocks (15, 16). The topological defect (TD) size, d, can be expressed as the TD field mass scale mTD = ℏ/(dc). The constraints in the top left panel have a discontinuity at T ≈ 3000 s because the ri(t) from the most stable clock in our network (NIST, Boulder, CO, USA) was collected over a shorter duration. The full range of the defect size, d, concerns the domain-wall type of TDs, whose transverse dimensions (regardless the d parameter) are much larger than the span of our network.

Equation 1 can be translated into a lower bound for the energy scale Λγ,2 of the quadratic DM-SM couplingEmbedded Image(3)where ρDM is the local mean cold DM energy density in our galactic neighborhood and T is the time between consecutive encounters with topological defects. We assume that the DM field inside a defect is ϕ2 = TvρDMdℏc (8) and zero outside the topological defects. Under the assumption of ρDM = 0.468 GeV cm−3 (3, 36), we derived the constraints on the energy scale Λγ,2, plotted in the top left and bottom panels of Fig. 3 as a function of T and d. The constraints in the top left panel have a visible discontinuity at T ≈ 3000 s. Our sensitivity for longer times between consecutive encounters T is lower because the ri(t) from the most stable clock in our network (NIST, Boulder, CO, USA) was collected over a shorter duration.

In the case of oscillating massive scalar fields, periodic perturbations will occur with the same phase in all ri(t) signals. The strength of the linear coupling is characterized in the literature (10, 13, 14, 33) by the appropriate dimensionless coefficientEmbedded Image(4)where MPc2 = 1.2 × 1019 GeV is the Planck mass energy equivalent (15).

We applied the procedure used in (13, 14) and described in (37) to search for harmonic-oscillation signatures. For each oscillation frequency ω (i.e., for each mass mϕ = ℏω/c2), we fit the normalized signals Embedded Image from all our detectors with an oscillating function R(t) = B + A(ω)cos(ωt + δ) function, where B is a constant offset (the contributions were weighted with inverse corrected sample SD for each Embedded Image). The fitted A2(ω) and the corresponding 5 and 95% CLs as well as the threshold levels are shown in the left panel of Fig. 4. The periodogram allows identification of the noise-type dominant in our data (37). The dependence of the A2(ω) on frequency is consistent with pink (flicker) noise. The noise analysis used to determine the confidence and threshold levels is described in Materials and methods.

Fig. 4 Upper limits on the coupling constant de.

Left: Periodogram of the square of the fitted amplitude A(ω) (light blue). The A2(ω) distribution in the data is consistent with pink (flicker) noise (37). Green and red curves depict the 5 and 95% CLs, respectively. The black curve corresponds to the mean value based on the assumed noise model. The brown curve is the detection threshold. Right: Limit on de derived from our measurement. The blue curve depicts the observed fit. The red curve depicts our limit on de, i.e., the 95% CL. The magenta and orange lines are the 95% confidence limits reported in (13) and (14), respectively. The green line is the result derived from the stability analysis of a single optical atomic clock (15). The brown line represents the equivalence principle tests such as the Eöt-Wash experiment (38) and Lunar Laser Ranging (39). The black line represents the first results of the MICROSCOPE experiment (40).

The corresponding limit on the strength of the DM-SM coupling can be calculated as [see details in (13)]Embedded Image(5)where G is the Newton’s constant. The right panel in Fig. 4 depicts the limit on de derived from our measurements as a function of the scalar field mass mϕ. The red line represents the constraints at the 95% CL. The coherence time of the field is longer than the duration of our observations for masses lower than 10−15.4 eV c−2. In the considered range of mϕ, our results reached the same levels as the limits reported in (13) (magenta line) and improved those of (14) and (15) (orange and green lines, respectively). The reported results for the oscillating massive scalar fields are complementary to the limits derived from the weak equivalence principle tests (3840).


We have demonstrated operation of an Earth-scale network of optical atomic clocks. Different susceptibilities to the fine-structure constant of the atoms and the cavity allowed the intercontinental comparison of different clocks without real-time frequency transfer. Our measurements aim at searching for oscillations and transient variations in the fine-structure constant, which may originate from the presence of DM in the form of topological defects or massive scalar fields. We improved the previous constraints on transient variations in α by two orders of magnitude. The approach demonstrated here has a large potential for further improvements. Longer, coordinated measurement intervals and the incorporation of other optical atomic clocks into the network will allow for further improvements in the detection limits and for applying the determined constraints to larger volumes in the parameter space.


Datasets ri(t)

We used ri(t) from four optical lattice atomic clocks located at NIST, Boulder, CO, USA (23, 24), at LNE-SYRTE, Paris, France (25, 26), at KL FAMO, Torun, Poland (27, 28), and at NICT, Tokyo, Japan (29, 30). During the measurement sessions, the reported fractional clock instabilities at 1 s were equal to 3 × 10−16, 1 × 10−15, 2 × 10−14, and 7 × 10−15, and the reported fractional optical cavity instabilities at 1 s were equal to 2 × 10−16, 8 × 10−16, 2 × 10−14, and 2 × 10−15, in the NIST, LNE-SYRTE, KL FAMO, and NICT laboratories, respectively. The measurements were performed from 9 January 2015 to 17 December 2015. The measurements were collected for 11, 24, 42, and 54 days in the NIST, LNE-SYRTE, KL FAMO, and NICT laboratories, respectively. Together, the laboratories collected data for 114 days. In our topological defect analysis, we used the longest consecutive overlapping records for each pair of clocks. Their length is equal to 3336 and 13,533 s for NIST and LNE-SYRTE and for LNE-SYRTE and KL FAMO, respectively. The constraints reported in Fig. 3 are meaningful only if they are referred to the lengths of the analyzed signal. In oscillating massive scalar field analysis, we used all our data spanning 114 days from all the laboratories. One clock cycle lasts 0.5 s in NIST, 1 s in LNE-SYRTE, 1.5 s in NICT, and 1.3 s in KL FAMO. The clock cycle duration restricts the size of the topological defects that we can detect with the network of unsynchronized clocks. In the case of oscillating scalar fields, the results are limited by the feedback servos’ time constants, which vary from 2 to 20 s depending on the clock.

In our oscillating field analysis, we weighted ri(t) from each clock with their standard deviations. All ri(t) are heavily affected by low-frequency cavity drifts. Before the analysis, we applied a high-pass filter to eliminate the low-frequency cavities’ instabilities. In the case of topological defect analysis, we applied an FIR high-pass Gaussian filter with frequency cutoff at 0.005 to 0.027 Hz to all ri(t). This choice of transmission characteristics increases the sensitivity to the considered range of events. In the oscillating massive scalar field analysis, the main contribution to the signal instabilities are cavity drifts in a much longer time scale. We filtered them by subtracting a second-order polynomial fit from continuous sections of ri(t). To verify the filter response, we added artificial harmonic oscillations with specified parameters to nonfiltered data, and we tested how the filter performs on different frequencies.

Topological defects: Data analysis

For nonzero DM-SM coupling, a topological defect will manifest as a transient perturbation of the fine-structure constant. We denoted its magnitude and duration as δα and T, respectively. The corresponding perturbation of the relative cavity-atom frequency encoded in ri(t) depends on the atom sensitivity Kα and for two different clocks A and B can be written asEmbedded Image(6)

We defined normalized signals as Embedded Image. In the simple case involving two clocks, the goal of our data analysis was to determine the constraint on δα from Embedded Image and Embedded Image. This can be performed by cross-correlating Embedded Image and Embedded ImageEmbedded Image(7)where tBtA is the duration of the cross-correlated signals. For square events, the shape of their cross correlation is triangular with full width at half maximum equal to T and amplitudeEmbedded Image(8)

This gives the constraint on the transient variations in the fine-structure constant (see Eq. 1). For the quadratic coupling, Embedded Image, Eq. 3 directly follows from Eq. 1; note that the DM field inside a defect is Embedded Image (8). The DM field is assumed to be zero outside the topological defects.

The simplest model of DM topological defects requires at least two free parameters whose values have to be arbitrarily chosen; in our case, these parameters are the defect size, d, and the time between consecutive encounters with topological defects, T. We repeated our analysis for every combination of these two parameters within the ranges we accessed with our network. For given values of d and T, we cross-correlated the two Embedded Image and fit the cross correlation with the expected triangular shape. We repeated this fitting procedure for any accessible time and any time delay between the two laboratories. The maximal value of the fitted amplitude, A0, yields the noise level corresponding to the transient variations in α (Eq. 1), and strength of the DM-SM coupling (Eq. 3) (see black lines in Fig. 3). This approach is analogous to the VLBI technique used in radio astronomy.

To give a proper statistical interpretation of our results, we also calculated the 5 and 95% CLs (see the green and red lines, respectively, in Fig. 3). To calculate the CLs, we repeated (100 times) the above algorithm of maximal A0 determination, artificially shifting one of the Embedded Image in time, ensuring that no expected common physical signal is present [for every repetition, we shift Embedded Image by a different time]. Having this ensemble of 100 cross-correlation amplitudes (for a given d), we selected the five smallest and largest elements, which yields the 5 and 95% CLs. We treated the 95% CL as an actual constraint determined in this work; hence, the black line in the top right panel in Fig. 3 should statistically (i.e., without any common component) exceed the red line 95% CL one time per every 20 values in d (the step in d is 300 km). This is consistent with Fig. 3, which indicates that at present level of accuracy, we did not observe any signature of hypothetical DM.

In the dataset considered here, the time of overlap of at least three clocks was considerably smaller than the time of overlap of two clocks. Therefore, without losing the strength of the constraint, we could restrict our analysis to simultaneous comparison of two clocks at a given time. For more than two simultaneous Embedded Image, the approach presented here can be generalized in a straightforward manner. For any value of the model free parameters and for any value of the DM velocity, the expected DM signature can be simultaneously fit to all the cross correlations between the participating clocks. In our analysis, we chose ηT = 1, i.e., the lengths of the cross-correlated Embedded Image segments are equal to the event duration, which is optimal for this analysis. It should be mentioned that the sensitivity of the clocks to DM objects could drop by a factor of 2 for typical servo-loop settings when the DM object duration is comparable to one servo-loop cycle, Tc.

Oscillating massive scalar fields: Data analysis

For each frequency ω, i.e., for each mass mϕ, we performed the linear least-squares analysis of the normalized signal Embedded Image composed of the Embedded Image contributions from all the four laboratories. We fit to the function R(t) = B + A(ω)cos(ωt + δ), where B is a constant offset (the contributions from different laboratories are weighted with inverse variance). In the left panel in Fig. 4, the square of the fitted amplitude is shown as a function of ω. To interpret the obtained result, one has to verify whether the fitted amplitude considerably exceeds the noise level. First, let us consider the simplest case when the noise of Embedded Image is white. The square amplitude of the harmonic component, A2(ω), can be treated as a new random variable. For the case of white noise [N samples of normalized signal Embedded Image with SD σ], its expected value does not depend on ω, 〈A2(ω)〉 = 4σ2/N, and the cumulative probability distribution function (CDF) can be expressed by a simple analytical formula (37)Embedded Image(9)Then, the A2 corresponding to X CL can be expressed asEmbedded Image(10)

The statistical interpretation of Eq. 10 is that the probability that Embedded Image is X. For instance, for X = 95%, the corresponding A2 level is Embedded Image. Equation 10 concerns a single frequency (one of the N/2 frequency channels available in the spectrum). That is, 1 per 20 frequency channels on average should exceed the Embedded Image level. Beyond the 95% CL criterion, we also defined a detection threshold for Nf frequency channels, Embedded Image, in a similar manner as performed in (13) and (37): The random variable, A2, exceeds the Embedded Image level for any frequency channel, on average, only 0.05 times per measurement. That is, if the measurement were repeated 100 times, then the noise would be interpreted as a positive detection, on average, only five times. To make this definition unique, we set an additional condition that the probability of exceeding the detection threshold is the same for all the frequencies. For the case of white noise, the detection threshold defined here can be expressed asEmbedded Image(11)

As an Nf, we took the total considered frequency range (see Fig. 4) divided by Δω = 2π/Ttot, where Ttot is a total time of our measurements (i.e., the difference between the time of the last and first samples in the data combined from all the laboratories). Equations 10 and 11 are strictly valid for white noise. We tested with the Monte Carlo simulations, however, that Eqs. 10 and 11 well approximate the CL and detection threshold for other types of noise; for this, one should replace the white noise parameter 〈A2〉 with the noise model function 〈A2〉(ω). For instance, in the case of pink noise, 〈A2〉(ω) ∝ 1/ω, the A2 corresponding to 5 and 99% CLs calculated as Embedded Image is indistinguishable from accurate Monte Carlo simulations at the scale of Fig. 4. A slight difference occurs for the detection threshold, but it can be easily eliminated with a numerically determined correction factorEmbedded Image(12)

We used the above expressions to determine the CLs and detection threshold. The periodogram in the left panel in Fig. 4 shows that in our measurements, the power distribution in the collected data is consistent with pink noise (37).

Alternative definition of the detection threshold

Another way of defining the detection threshold is that A2 evaluated from white noise does not exceed Embedded Image in any of the Nf frequency channels at the CL X. It can be expressed as Embedded Image. Therefore, the detection threshold for a white noise and for Nf frequency channels with CL X isEmbedded Image(13)which for X close to 1 can be written asEmbedded Image(14)

This is equivalent to Eq. 11.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: The “A next-generation worldwide quantum sensor network with optical atomic clocks” project is carried out within the TEAM IV Programme of the Foundation for Polish Science cofinanced by the European Union under the European Regional Development Fund. Support has been received from the project EMPIR 15SIB03 OC18. This project has received funding from the EMPIR Programme cofinanced by the Participating States and from the European Union’s Horizon 2020 Research and Innovation Programme. This project has received funding from NIST, NASA fundamental physics, and DARPA QuASAR. We acknowledge funding support from the Agence Nationale de la Recherche (Labex First-TF ANR-10-LABX-48-01), Centre National d’Études Spatiales (CNES), Conseil Régional Île-de-France (DIM Nano’K). The contribution of P. Wcisło is supported by the National Science Centre, Poland, project no. 2015/19/D/ST2/02195, and the Polish Ministry of Science and Higher Education “Mobility Plus” Program. P.M. is the JSPS International Research Fellow. The contribution of M.B. is supported by the National Science Centre, Poland, under QuantERA programme no. 2017/25/Z/ST2/03021. Calculations have been carried out using resources provided by the computer cluster founded by the Polish National Science Centre under grant no. 2016/21/D/ST4/00903. The UMK measurements were performed at the National Laboratory FAMO (KL FAMO) in Toruń, Poland, and were supported by a subsidy from the Polish Ministry of Science and Higher Education. This work was performed under the following financial assistance award, 70NANB18H006, from the U.S. Department of Commerce, NIST. Author contributions: P. Wcisło conceived and designed the concept of global network of optical atomic clocks. P. Wcisło and M.Z. developed the concept of data analysis. K.B., S.B., M.B., R.B., R.F., R.C., H.H., T.I., J.L., A.L., W.M., P.M., D.N., M. Schioppo, R.L.T., X.Z., and M.Z. developed the experimental setups and performed the measurements. P. Wcisło, B.Z., and M. Sekido developed the correlation analysis method. P. Wcisło, B.Z., P.A., and M.Z. analyzed all the data. P. Wolf provided the supporting theory and statistical interpretation of the analysis. All authors contributed to the technical discussions and writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. All data supporting the findings of this study are archived on storage drives located at NIST, NICT, LNE-SYRTE, and UMK and may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article