## Abstract

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.

## INTRODUCTION

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 (*1*–*3*), 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} cm^{2} (*5*). New stringent upper limits were recently reported in the 8 to 125 keV/c^{2} mass range that excluded couplings to electrons with coupling constants of *g*_{ae} > 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*, *12*–*16*). 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.

## RESULTS

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 *r*_{i}(*t*) between the atomic transition frequency, *v*_{atom}, and the cavity-stabilized laser, *v*_{cavity}, 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 α (*v*_{atom}/*v*_{cavity} ∝ α) (*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 *n*th portal (*32*).

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 *r*_{i}(*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 , where is the sensitivity coefficient [equal to one in the nonrelativistic case (*31*)] and 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 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 *L* ≈ *d*/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 × 10^{2} and 3 × 10^{4} km. We cross-correlate the pairs of and fit the resulting data with the expected DM signature (see Materials and methods). The left panel in Fig. 2 shows two 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.

The constraint on the transient variation of the fine-structure constant can be expressed as(1)where *A*_{0} 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 constant(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*).

Equation 1 can be translated into a lower bound for the energy scale Λ_{γ,2} of the quadratic DM-SM coupling(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} = T*v*ρ_{DM}*dℏ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 *r*_{i}(*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 *r*_{i}(*t*) signals. The strength of the linear coupling is characterized in the literature (*10*, *13*, *14*, *33*) by the appropriate dimensionless coefficient(4)where *M*_{P}*c*^{2} = 1.2 × 10^{19} 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*_{ϕ} = ℏω/*c*^{2}), we fit the normalized signals 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 ). The fitted *A*^{2}(ω) 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 *A*^{2}(ω) 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.

The corresponding limit on the strength of the DM-SM coupling can be calculated as [see details in (*13*)](5)where *G* is the Newton’s constant. The right panel in Fig. 4 depicts the limit on *d*_{e} 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 (*38*–*40*).

## DISCUSSION

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.

## MATERIALS AND METHODS

### Datasets *r*_{i}(*t*)

We used *r*_{i}(*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 *r*_{i}(*t*) from each clock with their standard deviations. All *r*_{i}(*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 *r*_{i}(*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 *r*_{i}(*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 *r*_{i}(*t*) depends on the atom sensitivity *K*_{α} and for two different clocks *A* and *B* can be written as(6)

We defined normalized signals as . In the simple case involving two clocks, the goal of our data analysis was to determine the constraint on δα from and . This can be performed by cross-correlating and (7)where *t*_{B} − *t*_{A} 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 amplitude(8)

This gives the constraint on the transient variations in the fine-structure constant (see Eq. 1). For the quadratic coupling, , Eq. 3 directly follows from Eq. 1; note that the DM field inside a defect is (*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 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, *A*_{0}, 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 *A*_{0} determination, artificially shifting one of the in time, ensuring that no expected common physical signal is present [for every repetition, we shift 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 , 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 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, *T*_{c}.

### 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 composed of the 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 is white. The square amplitude of the harmonic component, *A*^{2}(ω), can be treated as a new random variable. For the case of white noise [*N* samples of normalized signal with SD σ], its expected value does not depend on ω, 〈*A*^{2}(ω)〉 = 4*σ*^{2}/*N*, and the cumulative probability distribution function (CDF) can be expressed by a simple analytical formula (*37*)(9)Then, the *A*^{2} corresponding to *X* CL can be expressed as(10)

The statistical interpretation of Eq. 10 is that the probability that is *X*. For instance, for *X* = 95%, the corresponding *A*^{2} level is . 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 level. Beyond the 95% CL criterion, we also defined a detection threshold for *N*_{f} frequency channels, , in a similar manner as performed in (*13*) and (*37*): The random variable, *A*^{2}, exceeds the 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 as(11)

As an *N*_{f}, we took the total considered frequency range (see Fig. 4) divided by Δω = 2π/*T*_{tot}, where *T*_{tot} 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 〈*A*^{2}〉 with the noise model function 〈*A*^{2}〉(ω). For instance, in the case of pink noise, 〈*A*^{2}〉(ω) ∝ 1/ω, the *A*^{2} corresponding to 5 and 99% CLs calculated as 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 factor(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 *A*^{2} evaluated from white noise does not exceed in any of the *N*_{f} frequency channels at the CL *X*. It can be expressed as . Therefore, the detection threshold for a white noise and for *N*_{f} frequency channels with CL *X* is(13)which for *X* close to 1 can be written as(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.

## REFERENCES AND NOTES

**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.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).