Analysis of atomic-clock data to constrain variations of fundamental constants

We present a new framework to study the time variation of fundamental constants in a model-independent way. Model independence implies more free parameters than assumed in previous studies. Using data from atomic clocks based on $^{87}$Sr, $^{171}$Yb$^+$ and $^{133}$Cs, we set bounds on parameters controlling the variation of the fine-structure constant, $\alpha$, and the electron-to-proton mass ratio, $\mu$. We consider variations on timescales ranging from a minute to almost a day. In addition, we use our results to derive some of the tightest limits to date on the parameter space of models of ultralight dark matter and axion-like particles.


Introduction
The standard model and general relativity are the currently accepted theories of elementary particles and gravity.Their predictions are largely controlled by a set of free parameters known as fundamental constants, which are extracted experimentally and assumed independent of time and spatial position.The underlying origins and potential spacetime variability of the fundamental constants have been rich subjects of investigation, dating back to Dirac's large numbers hypothesis [1,2].
In the years since, dynamical mechanisms, for example from string theory, have been suggested to explain the constants' origins.In some regimes, additional scalar fields imply spacetime variations [3][4][5][6][7].More generally, many realistic models give rise to variations of fundamental constants (comprehensive reviews may be found in, e.g., Refs.[8,9]).Models based on quantum field theory, such as Bekenstein [10] or Barrow [11,12] models can describe fundamental-constant variations in terms of low-energy effective interactions of additional scalar fields coupled to the standard model.More recently, variations of the dimensionless constants due to so-called ultralight fields, which could be related to dark matter [13][14][15], have revived significant attention in this idea (for a recent review see, e.g., [16]).Investigating the variability of fundamental constants is strongly tied to the foundational assumptions and outstanding problems in modern physics.This work reports on new developments in studying temporal variations of the dimensionless fundamental constants.In Sec. 2, a generic approach based on effective field theory, applicable particularly to the analysis of temporal variations, is described.The explicit spacetime dependence of a generic dimensionless constant is represented by a series expansion of a scalar field normalized to the energy scale of new physics responsible for the time variation.A description of scalar-field evolution, including arbitrary damping effects, is given, along with the number of observable parameters.This setup covers a broad range of models describing temporal variations in the literature.It also has significant consequences for the interpretations of measurements, demonstrating that there are more free parameters than typically expected.
Many previous experiments have used the extreme precision of atomic clocks to investigate time-variations of fundamental constants.Examples include constraints placed on linear-in-time variations [17][18][19][20], oscillations [21][22][23][24][25] and transients [26,27].For cases where the experimental data has been interpreted to place bounds on theoretical parameters such as coupling constants, specific models for the scalar fields and their dynamics have always been assumed.Our approach in this paper, however, allows constraints on time-variation to be presented in a model-independent way, without assuming the functional form of the scalar field.The model-independent bounds can be interpreted in terms of specific models due to the framework presented in Sec. 2, where we work with the most generic functional form for the scalar field and effective field theory methods.
In Sec. 3, measurements are presented from atomic clocks ( 87 Sr, 171 Yb + optical and 133 Cs microwave) at the National Physical Laboratory (NPL), assessing the temporal variability (stability) of the fine-structure constant, α, and the electron-to-proton mass ratio, µ = m e /m p , over a period of about two weeks.The frequency ratios between 171 Yb + , 87 Sr and 133 Cs clock transitions place constraints on oscillations in α and µ at and beyond the previous state-of-the-art.
Using these data, in particular the 171 Yb + / 87 Sr and 87 Sr/ 133 Cs ratios, in Sec. 4 modelindependent constraints are placed for the first time on low-dimensional couplings of an ultralight scalar field to matter.In Sec. 5, new constraints are extracted and compared with previous results for the special cases of ultralight scalar and axion-like dark matter.The prospects of future measurement campaigns are discussed in Sec. 6.

Field theory description of varying constants
In this section we introduce our generic framework to describe the spacetime variation of fundamental constants.This framework relies on the concept of effective field theory methods (see e.g.[28]).It is a tool of quantum field theory that has been extensively applied to many areas including particle, nuclear, atomic, condensed matter and gravitational physics.The significance of this approach is that it enables calculations and predictions, which can be tested in experiments, that are generic and universal.
Describing the spacetime variation of a fundamental coupling constant g(t, ⃗ x) can be accomplished by promoting the constant to a series expansion involving a scalar field ϕ(t, ⃗ x): g(t, ⃗ x) = g 0 + 1 Λ ϕ(t, ⃗ x) + . . ., where g 0 is the spacetime-independent contribution and Λ is the high-energy scale relevant to the onset of the physics responsible for the spacetime variation of constants.Coupling ϕ(t, ⃗ x) to conventional matter otherwise described by a coupling g 0 is thus accommodated by replacing g 0 → g(ϕ).
In this paper, we are primarily interested in time-varying fundamental constants and are thus considering only time-dependent scalar fields.The most generic field equation for a scalar field corresponds to a damped harmonic oscillator.From an effective field theory point of view, this is the field equation that corresponds to dimension-four operators: the kinetic term for the scalar field and potential interactions leading to a decay of the scalar field which can be parametrized by operators of dimensions 3 and 4. The equation of motion for a damped harmonic oscillator is given by where ϕ(t) is a time-dependent scalar field, m is the mass and Γ is a damping factor.
The solution to the field equation depends on the boundary conditions.For oscillatory solutions, the classes of behavior are identified by the relation between m and Γ: where ϕ 0 , ϕ 0,1 , and ϕ 0,2 are amplitudes, θ is a phase and Standard effective field theory methods can be used to describe the interactions of ϕ to conventional matter in a general manner.For example, the interactions of the scalar field with the photon field A µ and the electron field ψ e can be described by the Lagrangian: for positive integer n, where κ = √ 4πG = 1/( √ 2M P ) with G being the Newtonian constant of gravitation, M P being the reduced Planck mass, and Note that we normalize the interactions to the strength of the gravitational interaction, which is the weakest interaction in nature known to date.However, the dimensionless couplings d me control, respectively, the strength of interactions between ϕ and the photon and the electron.They may be taken as real numbers, and parameterize the magnitude of the time variation of the fine-structure constant and of the electron mass.The definition used in Eq. ( 6) sets the high-energy scale Λ from Eq. ( 1) such that κ n d are the combinations of parameters that determine the couplings of the scalar field to the different matter components.
The linear scalar-field coupling (n = 1) and quadratic scalar-field coupling (n = 2) are the most relevant as they are the lowest order effective operators and thus their effects are expected to be the strongest.Interactions with n ≥ 3 are suppressed by more powers of the scale of the physics responsible for the time variation and thus very much suppressed if this is a high-energy scale of the order of the Planck scale, for example.
Note that there are a number of free parameters that need to be fitted to data: γ ).We can use this fully generic approach to calculate the shift of the fine-structure constant α due to the scalar field ϕ(t) and obtain in the linear case (n = 1) In the physically relevant underdamped regime, we obtain In the quadratic case (n = 2), we have which leads in the underdamped regime to In all cases, we have five independent parameters: the coupling constants to matter d γ , an amplitude ϕ 0 , a damping factor Γ, an oscillation frequency ω d and a phase θ.Note that they may not all be measurable independently.A measurement of a change of α is only sensitive to the product d Similarly, the coupling of ϕ to gluons is controlled by the coupling d (n) g and that to quarks is controlled by the quark coupling constant d (n) m f .These free parameters control the degree of time variation of the proton mass and of quark masses.We note that the vacuum energy due to gluons, which corresponds to the Quantum Chromodynamics (QCD)-scale Λ QCD , accounts for approximately 90% of the nucleon mass M N .On the other hand, light-quark masses m u , m d and m s account for a mere ∼10% of the nucleon masses [29].The proton mass is thus much more sensitive to a time variation of Λ QCD than it is to a variation of the light quark masses if all of these parameters vary with time.For more details on the QCD couplings, please see Appendix A.
The main motivation for this methodology is that a broad range of models is encompassed by this effective field theory approach.By providing bounds on a time variation of fundamental constants using this formalism, bounds for specific models can easily be obtained and our results can then be interpreted in specific models, for example: • Quintessence-like models, see, e.g., [30][31][32]: in that case Γ = 3H where H is the Hubble parameter, but because today H ∼ 1/t today (where t today is the age of the universe), the damping is irrelevant today.Note that for quintessence models, the mass of the scalar field is of order H ∼ 10 −33 eV and there is some tension with torsion pendulum experiments (Eöt-Wash see, e.g., [33]) as these experiments exclude new bosons with masses lighter than 10 −2 eV for d j ∼ 1, one needs to consider models with d • Ultralight dark matter [13][14][15]: we must assume Γ = 0 as dark matter is stable.As described in detail in Sec. 5, if the scalar field with mass m accounts for all of dark matter, we can relate ϕ 0 to the local density of dark matter ρ DM : ϕ 0 ≈ √ 2ρ DM /m.Note that we could have a multi-component dark matter sector in which case the relation between ϕ 0 and ρ DM doesn't hold.This is further discussed in Sec. 5.
• The scalar field could be a generic scalar from some hidden sector [34].In that case, the amplitude is a free parameter and Γ ̸ = 0 if it can decay today.
• Kaluza-Klein and moduli models: in models with extra dimensions the sizes of compactified extra dimensions can be described by scalar fields (the moduli fields).If the size of extra dimensions changes with cosmological time, we could have a time evolution of these scalar fields.In particular, in string theory, all coupling constants are determined by the expectation values of moduli fields.Coupling constants could thus easily depend on time [3][4][5].
• Dilaton fields, see, e.g., [35][36][37]: they are similar to moduli, but we expect them to couple universally to matter like gravity does.These models include Brans-Dicke type fields and also scalar fields that are coupled non-minimally to the curvature scalar e.g.ϕ 2 R, where R is the Ricci scalar.
• Vacuum evolution models: in some extensions of the Standard Model, the vacuum expectation value of the Higgs boson can evolve with time [38,39].As the vacuum expectation value of the Higgs boson fixes the weak scale, it will lead to a time variation of all fermion masses as well as that of the electroweak bosons masses.This is typical of inflation-type models.
• Test of grand unified theories [3,[46][47][48][49]: in grand unified models, time shifts in α and the strong coupling constant α s (equivalently Λ QCD ) are related.If time variations of both α and µ were observed, predictions of grand unified theories could be directly probed.The same observation applies to shifts in lepton and quark masses.Because the relations between quark and lepton masses are very model dependent, clocks could help to determine the correct unification theory using very low energy data without the need to produce super massive particles.
We thus see that our generic theoretical framework enables us to study a very wide variety of models and also to test different scenarios of ultra-high-energy physics with very lowenergy experiments.This is obviously only a subset of the models that can be studied with these methods.Note that our approach enables us to describe both fundamental variations of constants as originally envisaged by Dirac, but also the effective time evolution of constants where the variation is due to interactions with additional fields, as in the case of, for example, ultralight dark matter.We point out that this framework can also be applied to massive spin-1 and spin-2 fields with only minor modifications in terms of the way these higher spin fields couple to matter.It should also be emphasized that while we looked at a damped oscillator model, there are important other classes of models that could have been considered.Depending on the boundary conditions, other solutions are possible, e.g., soliton models, transient phenomena, cosmic strings, domain wall, kink solution, etc.Great care needs be taken when interpreting data as new physics signals could easily get lost when a specific functional form of the signal is assumed.
Note that for infinitesimal time differences between two measurements of the fundamental constants, we can consider the first-order shifts For the data analysis performed in the next section, it is useful to consider fractional frequency shifts in atomic clocks δν/ν, where ν is the frequency of an atomic transition.Since these fractional frequencies depend on at least one fundamental constant, new physics in the form of g(ϕ) altering spectral widths could be imparted on laboratory measurements.As observables involve comparisons with a reference, in this context fractional frequency ratios are considered.These ratios are related to a linear combination of relative fundamental constant variations, where K g i is a sensitivity coefficient particular to the system of interest [50][51][52].The frequency of optical and microwave transitions can be expressed as where A and B are constants specific to a given transition, F opt (α) and F MW (α) are relativistic corrections and g N is the nuclear g-factor, c is the speed of light and R ∞ is the Rydberg constant.This work focuses on optical-to-optical Yb + /Sr and optical-to-microwave Sr/Cs ratios.The former are sensitive to α variations whereas the latter are sensitive to α, µ and g N variations. 4We now turn out attention to the data analysis.
3 Data and analysis framework

Experimental Setup
In this work we analyze frequency ratio data produced by atomic clocks based on neutral strontium atoms in a lattice trap (Sr) [55], a singly-charged ytterbium ion in a Paul trap (Yb + ) [56], and neutral cesium atoms launched in a fountain configuration (Cs) [57].The properties of the atomic and ionic energy transitions are summarized below in Table 1.
Table 1: Summary of the atomic and ionic energy transitions used to produce high-stability frequency data.K X values are taken from references [51,[58][59][60][61] and given to two decimal places.
The Sr, Yb + , and Cs clock frequencies are all measured relative to an active hydrogen maser (HM).For the optical clock frequencies from Sr and Yb + , this measurement is made via a frequency comb, referenced to the 10 MHz output of the maser.The microwave clock frequency from Cs, however, can be measured directly against the maser.As all the measurements are made during the same observation window, the frequency ratios Yb + /Sr, Yb + /Cs and Sr/Cs can be calculated in post-processing and are independent of the maser frequency.The achieved stabilities in the Yb + /Sr, Yb + /Cs and Sr/Cs frequency ratios are all close to the quantum projection noise limits, determined by the number of atoms interrogated, the clock cycle times and the quality of the local master oscillator [62].Fig. 1 displays the time series of data used in this work, plotted as fractional frequency ratios: with reference ratios, HM taken to be 10 MHz.The Sr and Cs data were available between 1 st -14 th July 2019, with total uptimes of 73% and 93% respectively; Yb + data were available between 1 st -6 th July 2019, with a total uptime of 76%.Fractional frequency ratio of Cs/HM data produced by NPL's cesium fountain, NPL-CsF2, with data pre-processed over 600 s intervals.Data were collected between 1 st -14 th July 2019 (MJD 58665 -58679).

Experimental Results: Frequency Ratio Instability
In the presence of general noise processes, the mean and standard deviation of data sets are not guaranteed to converge as the number of samples increases.It is therefore common practice to estimate the spread of r [i/j] over different averaging times, τ , using the Allan Deviation, σ r (τ ), and its extensions.These are more informative estimators of spread, as they converge for data sets exhibiting the most common kinds of non-stationary statistics [64,65].Specifically in this work, we characterise the instability of our data using the Modified Allan Deviation (MDEV).The MDEV is given by the square root of the Modified Allan Variance, which is defined for a data set of M measurements, y k , averaged over averaging time, τ , as where m is the averaging factor, such that τ = mτ 0 , and τ 0 is the original sampling time interval of the data points [65].As our frequency data were measured on a frequency counter configured in a Λ-counting mode, the MDEV was the appropriate statistic to characterise the instability [66][67][68].To directly compare instability estimates from different types of Allan deviation, the estimators may be converted between each other using the relations documented in [69].
Fig. 2 displays the values of σ r (τ ) calculated at octave intervals of averaging time.For 60 s ≤ τ ≤ 30 000 s, the σ r[Yb + /Sr] curve is well approximated by power-law behavior of σ r[Yb + /Sr] (τ ) ≈ 1.4 × 10 −15 / τ /s.This indicates that the data are well described on these timescales by a white frequency modulation (WFM) noise process with h 0 ≈ 8.3×10 −30 Hz −1 , where h 0 is a constant value of power spectral density [64,65].This is not true for τ < 60 s, because this is shorter than the time constant for steering the Yb + clock laser onto the atomic transition frequency.At these shorter timescales, correlated noise from the clock laser dominates the instability.For averaging times 600 s ≤ τ ≤ 80 000 s, the instabilities of Sr/Cs and Yb + /Cs are well approximated by σ r[Sr/Cs] ≈ σ r[Yb + /Cs] ≈ 1.6 × 10 −13 / τ /s, corresponding to WFM noise with h 0 ≈ 1.0 × 10 −25 Hz −1 .The small difference in instability between the Yb + /Cs and Sr/Cs ratios can be attributed to the fact that the two sets of data span different time periods with different uptimes.While σ r[Sr/Cs] appears to increase at the highest averaging time, we attribute this to error introduced by the routine used to calculate MDEVs, which interpolates gaps in the data record: when downtime is dominated by a small number of large gaps in the data record (as is the case here), this interpolation can introduce spurious drifts and inflate the instability estimate at the largest averaging times [70,71].Therefore we do not consider averaging times above τ = 80 000 s. Some publications [22] have attempted to constrain variations in fundamental constants using frequency data from hydrogen masers (HM), reasoning that ν HM shares the K X sensitivities of the 2 S 1/2 (F=0 − F=1) hyperfine transition in atomic hydrogen to which the maser cavity is tuned.We do not follow this approach in this work, as we cannot confirm over which timescales this reasoning holds true in our commercial maser system; the maser wall shift, the resonant cavity, the voltage-controlled crystal oscillator, etc. all contribute to the instability of ν HM over certain timescales and introduce sensitivity to additional variables,  e.g.temperature [72].Consequently, we do not use optical-to-maser ratio data (Sr/HM or Yb + /HM) directly in this work to place constraints on the variations in fundamental constants.Instead, we use optical-to-cesium ratio data, despite the higher WFM instability in the data sets, as we are confident in the sensitivities of ν Cs to variations in the fundamental constants.
On timescales for which the instability of the frequency ratio data is dominated by the behavior of the atomic transitions, we place constraints on fundamental constants using Eq. ( 12) and K X values taken from Table 1.Due to the negligible sizes of ∆K Yb + /Sr µ and ∆K Yb + /Sr q , we assume that the instability in fractional changes in α to leading order must be less than σ r[Yb + /Sr] scaled by the magnitude of the sensitivity ∆K 01.Thus, we constrain the instability of fractional changes in α to be σ (∆α/α) = σ r[Yb + /Sr] /6.01 ≤ 2.3 × 10 −16 / τ /s on timescales of 60 s ≤ τ ≤ 30 000 s.With σ (∆α/α) constrained to two orders of magnitude below the noise level of Sr/Cs and Yb + /Cs, we make the further assumption that any remaining instability in r [Sr/Cs] would be dominated by ∆µ/µ rather than ∆m q /m q , due to the small size of ∆K .00.Under this assumption, to leading order we may similarly constrain the instability of fractional changes in µ to be no greater than σ r[Sr/Cs] scaled by the magnitude of the sensitivity ∆K Sr/Cs µ = 1.00.Thus, we constrain the instability of fractional changes in µ to be σ (∆µ/µ) = σ r[Sr/Cs] /1.00 ≤ 1.6 × 10 −13 / τ /s on timescales of 600 s ≤ τ ≤ 80 000 s.These constraints on the instability of fractional changes in α and µ as a function of averaging times are shown in Fig. 3 and summarized in Table 2.
These estimates of fractional variations in frequency ratios, and hence α and µ, on different timescales make no assumptions about the functional form of the variations.These results can be translated into model-independent limits, which will be discussed in Sec. 4.

Experimental Results: Sinusoidal Oscillations
If one chooses to focus specifically on oscillatory time variations of fundamental constants, these can be generically described by a damped harmonic oscillator, given in Eq. ( 2).Constraints on damped oscillatory signals could be obtained for a range of oscillation frequencies, f , by finding the best fit amplitude of each oscillation frequency, A(f ), and the best fit to the damping factor, Γ.However, it was decided as a first stage of analysis to fit undamped oscillations (Γ = 0) to reduce the number of parameters that require fitting.If any significant signals or features were detected as a result of fitting to undamped oscillations, reasonable values of Γ could be inferred from the linewidths of any peaks, and further analysis could be performed to fit for Γ and constrain damped oscillations to the data.In the case that significant features were observed, we implemented a routine to fit the data to oscillations weighted by an envelope function, Z(t), which for Z(t) = exp(−Γt/2) would model underdamped oscillations [73] 5 .During testing it was observed that this routine did not perform better than standard periodograms for detecting damped oscillations in the parameter space of interest, so it was decided to use a standard periodogram.
Similar to the approach taken in recent works [22,23,25,74], we constrain the magnitude of undamped oscillations in our data by estimating the power spectral density of the fractional frequency ratios, S r (f ), via the Lomb-Scargle periodogram (LSP) [75,76].The LSP is an estimator of S r (f ) for time series that suffer from irregular sampling or data gaps due to experiment downtime [77].Calculating the LSP is equivalent to performing linear leastsquares fits of a data set to the amplitude of sinusoids at a range of frequencies; it allows algorithms in the style of a fast Fourier transform to be used on time series with incomplete or irregular sampling, without having to account for data gaps by deconvolving the time series with composite window functions [76,78].
Power spectral density estimates for the fractional frequency ratios, S r (f ), were calculated using the implementation of the LSP [75,76] provided in the Astropy Python package [77].The Nyquist frequency, f Ny = (2∆t sample ) −1 , was chosen as the bandwidth upper limit for ease of comparison with previous works, and any future works that may have 100% data uptime.However, S r (f ) is not guaranteed to be free of aliasing for f < f Ny when signals are irregularly sampled, as gaps in data sampling introduce spectral leakage below f Ny [75].The resolution of the frequency grid was not tuned beyond using the default oversampling factor of n 0 = 5 [77].The fidelity of LSP in detecting oscillations in noisy data was validated by injecting sinusoidal signals into data sets with similar noise statistics as those of the observed data.
Fig. 4 shows S r (f ) for each frequency ratio, with black dashed lines indicating the estimated noise levels, h 0 , calculated from values of σ r (τ ) [69].Here it can be seen again that while r [Sr/Cs] and r [Yb + /Cs] are well described by WFM noise, this is not true of r [Yb + /Sr] for f ≥ (60 s) −1 = 1.67 × 10 −2 Hz, where more complex power-law behaviour can be seen: the action of the servos steering the probe light frequency onto the atomic/ionic transition frequencies leads to an approximate S r (f ) ∝ f −1.5 behavior.Limits on oscillations in ∆α/α 5 Since our longest data set had a length of T = 14 days and our minimum sampling time was ∆t = 60 s, the approximate range of Γ detectable with our data was between Γ min ∼ 1/T ≈ 8.3 × 10 −7 Hz and Γ max ∼ 1/∆t ≈ 1.7 × 10 −2 Hz.The range of detectable values for Γ corresponds to a range of lifetimes between 60 s to 1.2 × 10 6 s, using τ * = Γ −1 .It is interesting to compare this range to lifetimes corresponding to known forces of nature.Typical lifetimes in Quantum Electrodynamics are of the order of 10 −20 s to 10 −16 s.For the weak interactions, one finds 10 −13 s to 10 3 s while for the strong force one finds 10 −23 s to 10 −20 s.Thus the range of detectable lifetimes in this work partly overlaps with lifetimes typical of the weak force.The levels of white frequency modulation (WFM) noise, h i 0 , are taken from the best fits to the σ r (τ ) curves presented in Fig. 2. in the frequency range that exhibits WFM noise can be formed by integrating S r[Yb + /Sr] (f ) over the nominal frequency bin width (δf [Yb + /Sr] = 1/T [Yb + /Sr] = 2.5 × 10 −6 Hz) to obtain the total power of each bin, then taking the square-root to obtain the one-sided 6 fractional amplitude spectrum, A r[Yb + /Sr] (f ).
Constraints on oscillations in α and µ can be extracted from the fractional amplitude spectra A α (f ) and A µ (f ) by considering the statistics of the LSP estimator.In the absence of prominent peaks, some authors [23] have placed constraints at the observed power for each frequency bin of the spectra.Following the approach taken in [21,22,25,74] and others, confidence intervals on the constraints at each frequency bin were calculated by simulating a large number of control spectra with equivalent noise statistics to the observed spectra, then using the empirical cumulative distribution function of the simulated power in each  frequency bin to estimate confidence intervals, in this case, at 95% confidence.Whilst we follow this approach to calculate 95% confidence intervals on our power estimates, using these confidence intervals as "limits" has the disadvantage that the bounds placed on a WFM process could differ across neighboring frequency bins by an order of magnitude or more, and would likely fluctuate between repeat experiments even if the noise level in the experiment remained constant.
We believe a more appropriate and reproducible bound for a WFM spectrum is one that is constant across all frequencies: under the null hypothesis that the power differences across frequency bins are merely the result of fluctuations due to a WFM noise process, we produce a global bound based on an estimate of the white noise level, with exclusion limits estimated with false alarm levels, A X (p ≤ p 0 ).These false alarm levels represent the value of A X (f ) that would be exceeded with a probability of no more than p 0 across all frequencies in the case of white noise only.While analytic expressions exist for the distribution of S r (f ) for regularly sampled data (assuming uncorrelated Gaussian errors) [76,79], in this work we employ computational methods, as recommended in [77].
Simulating data sets with the same noise and data gaps as the observed spectra, we use the bootstrap method [80] to estimate bounds at the 68% significance level of A α (p ≤ 0.32) = 5.6 × 10 −18 and A µ (p ≤ 0.32) = 1.3 × 10 −15 .More appropriate to particle physics would be the equivalents of a 5σ significance level, which were estimated using the Baluev method [81], which yield A α (p ≤ 3.5 × 10 −7 ) ≈ 8.9 × 10 −18 and A µ (p ≤ 3.5 × 10 −7 ) ≈ 2.1 × 10 −15 7 .As shown in Fig. 5, all spectral peaks fall well below this fractional amplitude, though it should be admitted that this is a fairly strict detection threshold.Though global bounds across the entire frequency domain of the LSP estimate are more conservative than bounds one could achieve by setting constraints on individual frequency bins, for the reasons outlined above, we believe them to be better motivated and more legitimate for processes that appear to be predominantly WFM noise.
With only two peaks slightly exceeding the 1σ false alarm level for A µ (f ), and no spectral peaks exceeding the 5σ false alarm levels for either A µ (f ) or A α (f ), we have high confidence that all peaks in A α (f ) and A µ (f ) are consistent with instabilities due to WFM noise processes.Therefore, based on the estimated WFM noise levels of fractional frequency ratio data, we place constraints on sinusoidal oscillations in α at the 1σ significance level of A α (f ) ≤ 5.6 × 10 −18 and at the 5σ significance level of A α (f ) ≲ 8.9 × 10 −18 , over the frequency range 2.5 × 10 −6 Hz ≤ f ≤ 1.7 × 10 −2 Hz.We also place constraints on sinusoidal oscillations in µ at the 1σ significance level of A µ (f ) ≤ 1.3 × 10 −15 and at the 5σ significance level of A µ (f ) ≲ 2.1 × 10 −15 , over the frequency range 8.7 × 10 −7 Hz ≤ f ≤ 8.3 × 10 −4 Hz.A summary of these 1σ and 5σ constraints is shown in Table 3.

Amplitude Constraint Amplitude Constraint
Parameter space (1σ significance) (5σ significance) Table 3: Summary of constraints on the powers of fractional oscillations in α and µ produced in this paper, expressed as Fourier spectrum amplitude detection thresholds, A X (f ), at the 1σ and 5σ significance level.
These constraints on sinusoidal oscillations can be further interpreted in the case of specific models such as scalar dark matter, which is the focus of Section 5.

Model-independent constraints
As emphasized before, the main strength of our new theoretical description of a time evolution of constants is that it enables us to set model-independent constraints on the time variation of these constants.Independent of the shape of the function that describes the time variation, we now set some limits which can then trivially be interpreted in specific models.
From Eq. ( 13), a transition frequency ν may be parametrized as [82] where K α , K µ , and K q are sensitivity coefficients characteristic of the transition ν and m q ≡ (m u + m d )/2.The quark coefficient parametrizes changes in the nucleon mass δM N /M N = k M N q (δm q /m q ) and the nuclear magnetic moment δg N /g N = k g N q (δm q /m q ) in terms of m q variations.The ratio of two frequencies r = ν 1 /ν 2 is independent of the dimensionful constants and varying Eq. ( 16) with respect to α, m e /Λ QCD , and m q /Λ QCD noting (11) gives where q is the mass-weighted mean-quark coupling.Using the published values of K Yb + α = −5.95,K Sr α = 0.06, and K Cs α = 2.83 from Ref. [52] and K Cs µ = 1, K Cs q ≈ 0.07 from Ref. [82] (and summarized in Tab. 1) results in the fractional frequency ratios As can be observed in Fig. 2, constraints from Yb + /Sr measurements enable a bound on the coefficient d far below what is possible from Sr/Cs measurements.It is therefore appropriate to neglect d (n) γ in Eq. ( 19), which gives an effective coupling for the Sr/Cs measurements of Comparing Eqs. ( 18) and ( 19) with the data in Fig. 2, model-independent constraints can be placed on the magnitudes of the coupling strengths d Sr/Cs , and the instability of ϕ n (t) over different timescales, σ ϕ n (τ ) given by for timescales 60 s ≤ τ ≤ 30 000 s, and for timescales 600 s ≤ τ ≤ 80 000 s, where σ ϕ n (τ ) is the Modified Allan Deviation of ϕ n , defined in Eq. ( 15), and the contribution from (d g ) is assumed to be subdominant.Note that these constraints do not assume any specific fundamental physics model, nor do they make any assumption about the functional form of ϕ n (t).The constraints in the equations above are only valid for the specified values of τ explored in this work, and cannot constrain fluctuations on timescales outside this range.
The constraints on σ ϕ n (τ ) in Eqs. ( 21) and ( 22) can be roughly interpreted as limits on the average magnitude of fluctuations in ϕ n (t) between any two points in time t and (t + τ ).For example, for two times separated by τ = 1 000 s, the fluctuation [ϕ n (t + τ ) − ϕ n (t)] should roughly satisfy κ n |d Only once the functional form of ϕ is specified is it possible for independent constraints to be placed on the couplings, which we explore in the special case of ultralight dark matter in the following section.

Constraints on ultralight dark matter
One of the strongest cases for positing additional fundamental scalar fields is the problem of dark matter.Particle dark matter in the mass range 10 −22 eV ≲ m ϕ ≲ 1 eV is known as ultralight dark matter (ULDM), and in recent years significant efforts have been focused on detecting ULDM through apparent oscillations of fundamental constants.The upper bound m ϕ ≈ 1 eV occurs when the number density n of bosons in the reduced de Broglie volume (λ/2π) 3 satisfies n(λ/2π) 3 ≫ 1, resulting in a macroscopic phase-space occupation that exhibits Bose-Einstein condensation.Wavelengths λ ∼ O(kpc) span distances comparable to the smallest dwarf galaxy halos and imply a lower bound m ϕ ≈ 10 −22 eV [83].This bound also roughly coincides with the upper limit of dark matter being completely accounted for by m ϕ [84], which is a common assumption in studies seeking to exclude couplings at a given confidence level.
The Standard Halo Model is assumed for the dark matter density and velocity profiles, where for coordinates centered on the galaxy v vir ∼ 10 −3 is the virial speed (in natural units) of dark matter with isotropic distribution ⟨⃗ v vir ⟩ = ⃗ 0. As a result of the solar system's motion through the dark matter halo at speeds comparable to v vir , an Earth-based laboratory experiences a dark matter wind with | ⃗ k| ≈ m ϕ v vir ≪ m ϕ when neglecting subdominant corrections [85].In the dark-matter rest frame, oscillations are controlled by the rest mass m ϕ = 2πf ϕ where f ϕ is the Compton frequency in natural units.Oscillations are coherent in time for τ c ∼ 4π/(m ϕ v 2 vir ) ≳ 10 6 T c ≫ T data where the oscillation timescale T c = 1/f ϕ greatly exceeds the experimental timescale T data .As the coherence length λ ∼ 2π/(m ϕ v vir ) is larger than solar-system scales for all f ϕ of interest, the scalar-field amplitude is approximately constant.Under these conditions, ULDM is described by a macroscopic, nearly constantamplitude waveform oscillating at the underlying particle Compton frequency, up to small velocity corrections ∼ O(v 2 vir ).Measurements from the cosmic microwave background indicate dark matter was present in the early universe and is strongly constrained to be stable on experimental timescales [86].The standard theoretical treatment of ϕ as dark matter assumes cosmological evolution in a flat Friedmann-Lemaître-Robertson-Walker universe where (2) follows with Γ = 0.The field has an amplitude ϕ 0 = 2 ⟨ρ ϕ ⟩/m ϕ , resulting from time-averaging the energy density ρ ϕ = ( φ2 + m 2 ϕ ϕ 2 )/2 in the nonrelativistic limit.As all measurements presented in this paper were taken at a single location, the dynamics of ϕ are described by the solution Typically it is assumed that the scalar field comprises the entirety of the dark-matter density inferred at the solar galactocentric radius R 0 ≃ 8 kpc, or the "local density" ρ DM (R 0 ) ≈ 0.3 GeV/cm 3 [87].This figure should be taken with caution and has a large influence on the constraints, as densities from solar-system objects including planetary ephemerides [88] and more recently asteroid data [89] constrain ρ DM ≲ (10 4 -10 6 ) × ρ DM (R 0 ).The immediate implications are on the sensitivity to the couplings in Eqs. ( 18) and ( 19), since for a generic coefficient (assuming the field saturates the total density) the sensitivity to d and at best linear (quadratic) constraints would scale downward by ∼ 10 3 (∼ 10 6 ).However, they could also be substantially weakened.Constraints would need to be reconsidered in the case of multicomponent dark matter.For example, interactions between the different components could lead to nonzero decay widths (i.e.Γ ̸ = 0), and the density contribution would be spread across the components ρ DM = i ⟨ρ ϕ i ⟩.The extracted limits would be weaker and scaled by ⟨ρ ϕ i ⟩ /ρ DM .

Scalar couplings
Previous experiments monitoring electronic and microwave frequencies have resulted in constraints on the scalar couplings d and d (n) g assuming dark matter [15,[21][22][23][24]74,[90][91][92][93].In the present work, constraints may be extracted by using the amplitude spectra A X (f ) = A r (f )/|∆K| for frequency bin f and relevant sensitivity coefficient ∆K from Fig. 5. Identifying f = f ϕ = m ϕ /2π for linear couplings and f = 2f ϕ = m ϕ /π for quadratic couplings and noting Eq. (17) gives where f ) 2 .Note that the effects of boosting from the dark matter rest frame to the laboratory frame, which introduces a broadening of the oscillation frequency f .As a result, the sum over distinct field modes for measurement times T ≪ τ c reduces the sensitivity to the linear n = 1 coefficients by a stochastic factor ≈ 3 [94]. 8The data encompasses frequencies 10 −6 Hz ≲ f ≲ 2 × 10 −2 Hz, corresponding to masses 4 × 10 −21 eV ≲ m ϕ ≲ 8 × 10 −17 eV.
The results for linear and quadratic scalar couplings discussed in Sec. 2 are presented for the case of ULDM in Figs.6-9.Note that some works assume the slightly larger estimate ρ DM ≈ 0.4 GeV/cm 3 , which for purposes of comparison here amounts to a negligible difference.For all figures, the right axes compare κ n d (n) j = 1/Λ n j for a generic dimensionful scale Λ j .We choose the scale to be identified with Λ j , Λ ′ j for linear and quadratic couplings, respectively.For the coupling d γ in Fig. 6, a new exclusion region up to roughly an order of magnitude improvement over previously published work for 10 −20 eV ≲ m ϕ ≲ 10 −17 eV is observed.This improved sensitivity is mostly explained by the difference in K factors between, e.g., Rb/Cs where |∆K Rb/Cs α | ≈ 0.5 is roughly an order of magnitude smaller than |∆K Yb + /Sr α | ≈ 6.The extracted bounds also essentially surpass Eöt-Wash [95,96] and MI-CROSCOPE [97,98] equivalence-principle (EP) tests assuming a light dilaton.We note that recent experimental results [25], also using Yb + and Sr clocks, have claimed even tighter constraints than extracted here on the linear photon coupling d (1) γ in this mass range.
For the quadratic coupling d γ , a similar trend with respect to clock studies [90,91] using Rb/Cs [21] and 164 Dy/ 162 Dy spectroscopy [74] persists and using EP-test results to extract bounds on quadratic couplings [99] shows greater sensitivity ≳ 4 × 10 −18 eV.We also include constraints from big bang nucleosynthesis (BBN) [91] which surpass the sensitivity of other experiments in the probed mass range.
To the best of our knowledge, the only previously published clock-based studies to account for stochastic degradation factor ≈ 3 were those performed by the BACON collaboration with Al + /(Yb, Hg + ) and Yb/Sr clock comparisons [23], JILA using clock-cavity Sr/Si and H/Si comparisons [22] and NMJI using Yb/Cs clock comparisons [24].Note that EP tests do not rely on assumptions of the amplitude ϕ 0 and thus the contribution of the scalar field to the dark matter abundance.Similarly, though BBN constraints use ⟨ρ ϕ ⟩ = ρ DM , the field is non-oscillating with constant ϕ 0 for m ϕ ≪ 10 −16 eV so coherence considerations are irrelevant.
Limits on linear and quadratic d g .The final scalar constraints are extracted on the parameter A from the scalar-Higgs interaction and are presented in Fig. 9.The simplest (n = 1) linear couplings have garnered theoretical attention since they can emerge from the technically natural operator L ϕH = −AϕH † H for Higgs doublet H [100].The sensitivity coefficient is where b ∼ 0.2 − 0.5 is a dimensionless factor in the Higgs-nucleon Yukawa coupling g hN N = bm N /v, where v ≈ 246 GeV is the Higgs vacuum expectation value and where the nucleon mass m N = (m p + m n )/2 ≈ 0.94 GeV.To easily compare with existing Rb/Cs limits, we choose b = 0.2 and use the relevant Sr/Cs values from Table 1 of Ref. [90].Using the Sr/Cs spectrum and noting κd H ↔ A/m 2 h for Higgs mass m h ≈ 125 GeV in (24) produces the limit.For the ratios considered here, the bulk of the sensitivity comes from (1 − b)K µ as the γ (bottom panel).The best fit from Yb + /Sr (red) along with expected noise level (black dashed line) and 95% confidence level (C.L.) (light red) lines are displayed.Comparisons with constraints on linear couplings include Rb/Cs clocks [21], combined Al + /(Yb, Hg + ), Yb/Sr clocks [23], Yb/Cs clocks [24], Sr/Si and H/Si cavity comparisons [22], Dy/Dy spectroscopy [74] and EP tests [95][96][97][98].Comparisons with constraints on quadratic couplings include clocks [90,91], EP tests [99] and BBN [91].electromagnetic portion is suppressed by an additional factor of α (from radiative corrections) and the sensitivity of quark contributions to m N , g N are suppressed relative to m e /m N by around an order of magnitude.This implies, e.g., optical-optical Yb + /Sr and microwavemicrowave Rb/Cs ratios have weaker sensitivity to A relative to optical-microwave ratios.It is worth noting that b = 0.5 gives almost no sensitivity to Rb/Cs.Accordingly, Sr/Cs comparisons have more robust potential for constraining A, however, as the existing Rb/Cs data set is longer, competitive limits with respect to fifth-force searches exist in the range 10 −24 ≲ m ϕ ≲ 10 −20 .Note again that the current Rb/Cs limits do not include the stochastic degradation factor ≈ 3, which would shift the displayed curve upwards accordingly.g (bottom panel).The best fit from Sr/Cs (green) along with expected noise level (black dashed line) and 95% C.L. (light green) lines are displayed, including comparisons with EP tests and Rb/Cs [21], Yb/Cs [24], and H/Si [22].

Pseudoscalar couplings
Recently Kim and Perez [101] highlighted that atomic clocks can also provide complementary and competitive constraints on an axion-like field a coupled to gluons: where g s and f a are the strong coupling and axion decay constant, respectively.As a result, the square of the pion mass undergoes small oscillations quadratic in the field (bottom panel).The best fit from Sr/Cs (green) along with expected noise level (black dashed line) and 95% C.L. (light green) lines are displayed, including comparisons with EP tests and Rb/Cs [99], and BBN [90,91].
where m u , m d are the up and down quarks and θ = a/f a .The nucleon mass m N (θ) and the nucleon g-factor g N (θ) (and hence nuclear g-factor g(θ)) have an inherent dependence on m π (θ) as parametrized by chiral perturbation theory [61,[102][103][104].Similar to quadratic scalar couplings in the ultralight range, the oscillatory component of θ is given by where here ρ DM = 0.4 GeV/cm 3 , m a ̸ ∝ f −1 a in general and assuming a saturates the local density.The variation of a microwave transition frequency (13) may be expressed as variations in g(θ) and m p using Eqs.( 27) and (28), giving [86] δν Comparing with another microwave or optical standard and identifying m a = πf for signal frequency f from an amplitude spectrum yields the relation where m 15 ≡ m a /(10 −15 eV) and c r is a constant.In [101] microwave-microwave (Rb/Cs, c r ≈ 10 −1 ) and microwave-optical (H/Si, c r ≈ 1) comparisons are studied.One may also consider Sr/Cs, where due to the weak dependence of Sr on nuclear quantities |δr/r| ≈ δν Cs /ν Cs and c r ≈ 2 × 10 −2 , including effects of stochasticity.From Eq. ( 30), we plot the Sr/Cs limits along with Rb/Cs and H/Si from [101] as well as include the constraints from the neutron electric dipole moment (nEDM) in Fig. 10.
Perhaps most compelling are the high-frequency range constraints achievable on the axion-like coupling in Fig. 10 (see blue dashed line).In order to measure and remove density-dependent effects, the cesium fountain interleaves regular measurements with samples recorded with much higher atom densities.For the data presented here, this leads to a Cs clock cycle time of 600 s, corresponding to a maximum (Nyquist) frequency of 833 µHz.However, the cesium fountain could run with lower cycle times, and if the cesium fountain  [21], H/Si [22] and nEDM [105].
were operated without compensating density-dependent effects, cycle times of ≈ 5 s could be achieved, allowing Fourier frequencies up to ≈ 0.1 Hz to be probed, though at the expense of stability at longer times.This would yield constraints over an additional order of magnitude in frequency space outside of the large excluded nEDM region.This contrasts with the projections from future experimental proposals, demonstrating that existing atomic-clock capabilities can provide competitive constraints in axion physics.Taken together, subsequent studies based on future data campaigns are clearly motivated.

Discussion and Conclusion
In this work, we have presented a theoretical framework to describe the time variation of fundamental constants in a model-independent way.This approach demonstrates that a realistic model for the time variation of fundamental constants has many free parameters.
Using data acquired from atomic clocks operating at optical ( 87 Sr, 171 Yb + ) and microwave ( 133 Cs) frequencies, we constrain the instability of fractional changes in α to be σ (∆α/α) (τ ) ≤ 2.3 × 10 −16 / τ /s for averaging times 60 s < τ < 30 000 s, and we constrain the instability of fractional changes in µ to be σ (∆µ/µ) (τ ) ≤ 1.6 × 10 −13 / τ /s for averaging times 600 s < τ < 80 000 s.The theoretical framework then allows us to place constraints on combinations of a new scalar field ϕ n (t) and the different interaction strengths d (n) j coupling that field to other particles: photons, electrons, quarks and gluons.These constraints are independent of the underlying physics and the functional form of ϕ n (t).
As an example of a specific model, we studied ultralight dark matter couplings to matter and presented new constraints on low-dimension dilaton-like operators.The limits on d (1) γ from Yb + /Sr data exclude a new region of parameter space for masses 10 −20 eV ≲ m ϕ ≲ 10 −17 eV, as shown in Fig. 6.We also refer the reader to new experimental results [25] using Yb + and Sr clocks, which claim even tighter constraints on d (1) γ .In the future, the limits on most of the parameters presented in this paper could readily be extended into both higher and lower frequency regions.Higher frequency regions could be accessed by operating the clocks with shorter measurement cycles, and lower frequency regions could be accessed by recording longer data sets.Furthermore, there are clocks being developed that are more sensitive to variations of fundamental constants, such as certain highly-charged ion species [106] or molecular clocks [107].Frequency ratios between these new types of optical clock could place significantly lower bounds on many of the coupling strengths between new scalar fields and particles of the Standard Model [108].
One can disentangle the different parameters ϕ 0 and d j .By considering experiments sensitive to α, µ = m e /m p or α s , one could in principle measure ϕ 0 and some of the d j independently.Furthermore, in general there could be several scalar fields.Some could couple to photons others to gluons.To be very clear, couplings may not be universal.
The mass of the proton m p is mainly sensitive to the time dependence in the QCD coupling constant.Remember that the QCD scale is given by Λ QCD = µ r exp(2π/(α s (µ r ))) 1/b 3 with b 3 = −7 in the Standard Model and where µ r is the energy scale at which α s is measured.Neglecting a possible change in the quark masses, the proton mass m p is proportional to Λ QCD .Using the renormalization group equation for α s we find where in the linear case, we have in the underdamped regime αs α s = − κd for the quadratic case.

Figure 1 :
Figure 1: Top and middle: Fractional frequency ratios of Sr/HM and Yb + /HM data counted by NPL's femtosecond frequency comb, with data averaged over 1 s intervals.Bottom:Fractional frequency ratio of Cs/HM data produced by NPL's cesium fountain, NPL-CsF2, with data pre-processed over 600 s intervals.Data were collected between 1 st -14 th July 2019 (MJD 58665 -58679).

Figure 2 :
Figure 2: Characterisation of each fractional frequency ratio's instability using the Modified Allan Deviation, plotted at octave intervals of averaging time.

Figure 3 :Table 2 :
Figure3: Instability estimates for fractional changes in α and µ, where σ X (τ ) ≡ σ (∆X/X) (τ ) = σ r (τ )/K r X on timescales over which the instability is dominated by the behavior of the atomic transition frequency.

Figure 4 :
Figure 4: Power spectral densities of fractional frequency ratios, estimated with the Lomb-Scargle periodogram.The levels of white frequency modulation (WFM) noise, h i 0 , are taken from the best fits to the σ r (τ ) curves presented in Fig. 2.

Figure 5 :
Figure 5: Fractional amplitude spectra for sinusoidal oscillations in α and µ.The solid horizontal lines indicate the estimated levels of white frequency modulation (WFM) noise of the spectra, and the broken horizontal lines represent the estimated false alarm levels at the equivalent of 1σ (p ≤ 0.32) and 5σ (p ≤ 3.5 × 10 −7 ) significance.
in Fig.7and Fig.8, respectively.The linear constraints are competitive and similar in shape and magnitude with H/Si, Rb/Cs and Yb/Cs comparisons over the range of data, however the sensitivity of Rb/Cs extends up to two orders of magnitude below EP tests around m ϕ ≈ 10 −23 eV.Note that Rb/Cs has no sensitivity to d (n) me − d (n) g .Regarding the quadratic constraints, Sr/Cs data probes a new region for clocks in the d (2) me − d (2) g panel and displays roughly two orders of magnitude more sensitivity than EP tests at the low-mass range.Despite this, BBN still dwarfs the sensitivity by comparison and both EP tests and BBN encompass the range probed for d (2)

Figure 9 :
Figure 9: Constraints on Higgs coupling parameter A. The best fit from Sr/Cs (green) along with expected noise level (black dashed line) and 95% C.L. (light green) lines are displayed, including comparisons with Rb/Cs [90] and fifth-force searches [100].