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

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–7]. More generally, many realistic models give rise to variations of fundamental constants (comprehensive reviews may be found in, e.g. [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–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 section 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–20], oscillations [21–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 section 2, where we work with the most generic functional form for the scalar field and effective field theory methods.

In section 3, measurements are presented from atomic clocks (⁸⁷Sr, ¹⁷¹Yb⁺ optical and ¹³³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, $\mu = m_e/m_p$, over a period of about two weeks. The frequency ratios between ¹⁷¹Yb⁺, ⁸⁷Sr and ¹³³Cs clock transitions place constraints on oscillations in α and µ at and beyond the previous state-of-the-art.

Using these data, in particular the ¹⁷¹Yb⁺/⁸⁷Sr and ⁸⁷Sr/¹³³Cs ratios, in section 4 model-independent constraints are placed for the first time on low-dimensional couplings of an ultralight scalar field to matter. In section 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 section 6.

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,\vec{x})$ can be accomplished by promoting the constant to a series expansion involving a scalar field $\phi(t,\vec{x})$:

$$\begin{align} g(t,\vec{x}) = g_0+\frac{1}{\Lambda}\phi(t,\vec x)+\ldots \ , \end{align} \tag{ 1 }$$

where g₀ 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 $\phi(t,\vec{x})$ to conventional matter otherwise described by a coupling g₀ is thus accommodated by replacing $g_0 \rightarrow g(\phi)$.

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 parameterized by operators of dimensions 3 and 4. The equation of motion for a damped harmonic oscillator is given by

$$\begin{align} \ddot \phi + \Gamma \dot \phi+m^2 \phi = 0, \end{align} \tag{ 2 }$$

where $\phi(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 Γ:

$$\begin{align} \phi(t) = \begin{cases}\phi_{0,1} \exp{\left(r_1 t \right)} +\phi_{0,2} \exp{\left(r_2 t \right)}, \quad \Gamma^2 - 4 m^2 > 0\; \text{(overdamped)}, \\ \exp{\left(-\frac{\Gamma t}{2}\right)} (\phi_{0,1} + \phi_{0,2} t), \quad \Gamma^2 - 4 m^2 = 0\; \text{(critically damped)}, \\ \phi_0 \exp{\left(-\frac{\Gamma t}{2}\right)} \cos \left(\omega_d t - \theta \right), \quad \Gamma^2 - 4 m^2 < 0\; \text{(underdamped)}, \end{cases} \end{align} \tag{ 3 }$$

where φ₀, $\phi_{0,1}$, and $\phi_{0,2}$ are amplitudes, θ is a phase and

$$\begin{align} &r_{1,2} = -\frac{1}{2} \Gamma \pm \frac{1}{2} \sqrt{\Gamma^2- 4 m^2}, \end{align} \tag{ 4 }$$

$$\begin{align} &\omega_d = \frac{1}{2} \sqrt{| \Gamma^2- 4 m^2|}. \end{align} \tag{ 5 }$$

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:

$$\begin{align} \mathcal{L} = \left(\kappa\phi\right)^n\left(d^{(n)}_\gamma\tfrac{1}{4} F_{\mu\nu}F^{\mu\nu} - d^{(n)}_{m_e} m_e \bar \psi_e \psi_e \right), \end{align} \tag{ 6 }$$

for positive integer n, where $\kappa = \sqrt{4{\pi}G} = 1/(\sqrt{2}M_P)$ with G being the Newtonian constant of gravitation, M_P being the reduced Planck mass, and $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$. 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_\gamma^{(n)}$ and $d_{m_e}^{(n)}$ 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 equation (6) sets the high-energy scale Λ from equation (1) such that $\kappa^n d_j^{(n)} \equiv g_{0,j}/\Lambda_j^n$. In that sense, $\kappa^n d_{j}^{(n)}$ 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\geqslant 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:

$(\phi_0,\omega_d, \Gamma,\theta,d^{(n)}_{m_e}, d^{(n)}_\gamma)$. We can use this fully generic approach to calculate the shift of the fine-structure constant α due to the scalar field $\phi(t)$ and obtain in the linear case (n = 1)

$$\begin{align} \frac{\Delta \alpha}{\alpha}& = \kappa d^{(1)}_\gamma (\phi(t_2)-\phi(t_1)). \end{align} \tag{ 7 }$$

In the physically relevant underdamped regime, we obtain

$$\begin{align} \frac{\dot \alpha}{\alpha} = -\frac{\kappa d^{(1)}_\gamma \phi_0}{2} \exp{\left(-\frac{\Gamma t}{2}\right)} \bigg( \Gamma \cos\left( \theta - t \omega_{d}\right) + 2 \omega_{d} \sin\left( \theta - t \omega_{d}\right) \bigg). \end{align} \tag{ 8 }$$

In the quadratic case (n = 2), we have

$$\begin{align} \frac{\Delta \alpha}{\alpha}& = \kappa^2 d^{(2)}_\gamma (\phi(t_2)^2-\phi(t_1)^2), \end{align} \tag{ 9 }$$

which leads in the underdamped regime to

$$\begin{align} \frac{\dot \alpha}{\alpha} = -\frac{\kappa^2 d^{(2)}_\gamma\phi_0^2}{2} \exp{\left(-\frac{\Gamma t}{2}\right)} \bigg( \Gamma \cos\left( \theta - t \omega_{d}\right) + 2 \omega_{d} \sin\left( \theta - t \omega_{d}\right) \bigg). \end{align} \tag{ 10 }$$

In all cases, we have five independent parameters: the coupling constants to matter $d^{(n)}_\gamma$, an amplitude φ₀, 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^{(n)}_\gamma \phi^{n}_0$.

Similarly, the coupling of φ to gluons is controlled by the coupling $d_g^{(n)}$ and that to quarks is controlled by the quark coupling constant $d_{m_f}^{(n)}$. 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 $\Lambda_\textrm{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 $\Lambda_\textrm{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.

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–32]: in that case $\Gamma = 3H$ where H is the Hubble parameter, but because today $H\sim 1/t_\textrm{today}$ (where $t_\textrm{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\sim 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⁻² eV for $d_j^{(i)}\sim 1$, one needs to consider models with $d_j^{(i)}\ll1$.
• Ultralight dark matter (ULDM) [13–15]: we must assume $\Gamma = 0$ as dark matter is stable. As described in detail in section 5, if the scalar field with mass m accounts for all of dark matter, we can relate φ₀ to the local density of dark matter ρ_DM: $\phi_0 \approx \sqrt{2\rho_\textrm{DM}}/m$. Note that we could have a multi-component dark matter sector in which case the relation between φ₀ and ρ_DM does not hold. This is further discussed in section 5.
• The scalar field could be a generic scalar from some hidden sector [34]. In that case, the amplitude is a free parameter and $\Gamma \neq 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–5].
• Dilaton fields, see, e.g. [35–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. $\phi^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.
• Quantum gravity: quantum gravity predicts that $d_j^{(i)}\sim1$ [40–45]. Clocks have the capacity of probing the gravitationally generated interactions between very light scalar fields and photons or electrons (i.e. with $d_j^{(i)}\sim 1$) and in that sense they can probe quantum gravity.
• Test of grand unified theories [3, 46–49]: in grand unified models, time shifts in α and the strong coupling constant α_s (equivalently $\Lambda_\textrm{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 low-energy 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, ULDM. 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

$$\begin{align} \frac{\delta\alpha}{\alpha} &\equiv d_\gamma^{(n)}(\kappa\phi)^n, \\ \frac{\delta m_f}{m_f} &\equiv d_{m_f}^{(n)}(\kappa\phi)^n, \nonumber\\ \nonumber \frac{\delta \Lambda_{\text{QCD}}}{\Lambda_{\text{QCD}}} &\equiv d_g^{(n)}(\kappa\phi)^n, \end{align} \tag{ 11 }$$

with $f = e, u, d$.

For the data analysis performed in the next section, it is useful to consider fractional frequency shifts in atomic clocks $\delta\nu/\nu$, 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(\phi)$ 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,

$$\begin{align} \frac{\delta (\nu_1/\nu_2)}{\nu_1/\nu_2} = \sum_{g} (K_{g_1} - K_{g_2})\frac{\delta g}{g}, \end{align} \tag{ 12 }$$

where $K_{g_i}$ is a sensitivity coefficient particular to the system of interest [50–52]. The frequency of optical and microwave transitions can be expressed as

$$\begin{align} &\nu_\textrm{opt} = A\, F_\textrm{opt}(\alpha)\, cR_{\infty}, \nonumber\\ &\nu_\textrm{MW} = B\,\alpha^2 F_\textrm{MW}(\alpha)\, g_{N}\, (m_e/m_p) \, cR_{\infty}, \end{align} \tag{ 13 }$$

where A and B are constants specific to a given transition, $F_\textrm{opt}(\alpha)$ and $F_\textrm{MW}(\alpha)$ are relativistic corrections, g_N is the nuclear g-factor, c is the speed of light and $R_{\infty}$ 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 $\alpha, \mu$ and g_N variations ⁴ . We now turn out attention to the data analysis.

3.1. 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 [51, 58–61] and given to two decimal places.

Species	Transition	$K_{cR_{\infty}}+K_{\alpha}$	Kµ	Kq
87Sr	1S$_{0}- ^{3}$P0	$2+0.06$	0.00	0.00
171Yb+	2S$_{1/2}- ^{2}$F$_{7/2}$	$2-5.95$	0.00	0.00
133Cs	2S$_{1/2}(F = 3 \rightarrow F = 4)$	$2+2.83$	1.00	0.07

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]. Figure 1 displays the time series of data used in this work, plotted as fractional frequency ratios:

$$\begin{align} r_{[i/j]} = \frac{(\nu_{i}/\nu_{j} - R^{*}_{ij})}{R^{*}_{ij}} \, , \end{align} \tag{ 14 }$$

with reference ratios, $R^{*}_{ij} = \nu^{*}_{i}/\nu^{*}_{j}$, defined using the reference frequency values: $\nu^{*}_{\text{Sr}} = 429,228,004,229,$ $872.99$ Hz, $\nu^{*}_{\text{Yb}^{+}} = 642,121,496,772,645.12$ Hz, $\nu^{*}_{\text{Cs}} = 9,192,631,770$ Hz [63], and with the reference value $\nu^{*}_{\text{HM}}$ taken to be 10 MHz. The Sr and Cs data were available between $1{\textrm{st}}$ and $14{\textrm{th}}$ July 2019, with total uptimes of 73% and 93% respectively; Yb⁺ data were available between $1{\textrm{st}}$ and $6{\textrm{th}}$ July 2019, with a total uptime of 76%.

Zoom In Zoom Out Reset image size

3.2. 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, $\sigma_{r}(\tau)$, 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 characterize 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

$$\begin{align} \sigma^{2}_{y}(\tau) = \frac{\tau^{4}_{0}}{2\tau^4(M-3m+2)} \sum_{j = 1}^{M-3m+2} \left\{\sum_{i = j}^{j+m-1} \left( \sum_{k = i}^{i+m-1} [\, y_{k+m} - y_{k}] \right) \right\}^2 \, , \end{align} \tag{ 15 }$$

where m is the averaging factor, such that $\tau = m \tau_{0}$, and τ₀ 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 characterize the instability [66–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].

Figure 2 displays the values of $\sigma_{r}(\tau)$ calculated at octave intervals of averaging time. For 60 s $\leqslant\tau\leqslant$ 30 000 s, the $\sigma_{r[\text{Yb}^{+}\text{/Sr}]}$ curve is well approximated by power-law behavior of $\sigma_{r[\text{Yb}^{+}\text{/Sr}]}(\tau)\approx1.4\times 10^{-15}/$ $\sqrt{\tau /\textrm{s}}$. This indicates that the data are well described on these timescales by a white frequency modulation (WFM) noise process with $h_0 \approx 8.3 \times 10^{-30}$ Hz⁻¹, where h₀ 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 $\leqslant \tau \leqslant 80~000$ s, the instabilities of Sr/Cs and Yb⁺/Cs are well approximated by $\sigma_{r[\text{Sr/Cs}]} \approx \sigma_{r[\text{Yb}^{+}\text{/Cs}]} \approx 1.6\times 10^{-13}/\sqrt{\tau / \textrm{s}}$, corresponding to WFM noise with $h_0\approx 1.0 \times 10^{-25}$ Hz⁻¹. 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 $\sigma_{r[\text{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 $\tau = 80~000$ s.

Zoom In Zoom Out Reset image size

Some publications [22] have attempted to constrain variations in fundamental constants using frequency data from hydrogen masers (HMs), reasoning that $\nu_{\text{HM}}$ shares the K_X sensitivities of the ²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 $\nu_{\text{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 $\nu_{\text{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 equation (12) and K_X values taken from table 1. Due to the negligible sizes of $\Delta K_{\mu}^{\text{Yb}^{+}\text{/Sr}}$ and $\Delta K_{q}^{\text{Yb}^{+}\text{/Sr}}$, we assume that the instability in fractional changes in α to leading order must be less than $\sigma_{r[\text{Yb}^{+}\text{/Sr}]}$ scaled by the magnitude of the sensitivity $\Delta K_{\alpha}^{\text{Yb}^{+}\text{/Sr}} = | K_{\alpha}^{\text{Yb}^{+}}~-~K_{\alpha}^{\text{Sr}}| = 6.01$. Thus, we constrain the instability of fractional changes in α to be $\sigma_{(\Delta\alpha/\alpha)} = \sigma_{r[\text{Yb}^{+}\text{/Sr}]} / 6.01 \leqslant 2.3\times 10^{-16}/\sqrt{\tau / \textrm{s}}$ on timescales of 60 s $\leqslant \tau \leqslant 30~000$ s. With $\sigma_{(\Delta\alpha/\alpha)}$ 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_{[\mathrm{Sr/Cs}]}$ would be dominated by $\Delta \mu / \mu$ rather than $\Delta m_{q} /m_{q}$, due to the small size of $\Delta K_{q}^{\text{Sr/Cs}} = |K_{q}^{\text{Sr}}~-~K_{q}^{\text{Cs}}| \approx 0.07$ compared with $\Delta K_{\mu}^{\text{Sr/Cs}} = |K_{\mu}^{\text{Sr}}~-~K_{\mu}^{\text{Cs}}| = 1.00$. Under this assumption, to leading order we may similarly constrain the instability of fractional changes in µ to be no greater than $\sigma_{r[\text{Sr/Cs}]}$ scaled by the magnitude of the sensitivity $\Delta K_{\mu}^{\text{Sr/Cs}} = 1.00$. Thus, we constrain the instability of fractional changes in µ to be $\sigma_{(\Delta\mu/\mu)} =$ $\sigma_{r[\text{Sr/Cs}]} / 1.00$ $\leqslant 1.6\times 10^{-13}/\sqrt{\tau / \textrm{s}}$ on timescales of 600 s $\leqslant \tau \leqslant 80~000$ s. These constraints on the instability of fractional changes in α and µ as a function of averaging times are shown in figure 3 and summarized in table 2.

Zoom In Zoom Out Reset image size

Table 2. Summary of constraints on the instability of fractional changes in α and µ produced in this paper, expressed as Modified Allan Deviations, $\sigma_{(\Delta X /X)}(\tau)$.

Instability constraint	Parameter space
$\sigma_{(\Delta \alpha / \alpha)}(\tau) \leqslant 2.3\times 10^{-16}/\sqrt{\tau /\textrm{s}}$	60 s $\leqslant \tau \leqslant 30~000$ s
$\sigma_{(\Delta \mu / \mu)}(\tau) \leqslant 1.6 \times 10^{-13}/\sqrt{\tau /\textrm{s}}$	600 s $\leqslant \tau \leqslant 80~000$ s

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 section 4.

3.3. 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 equation (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 $(\Gamma = 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(-\Gamma t / 2)$ would model underdamped oscillations [73]. ⁵ 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 least-squares 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_{\text{Ny}} = (2\Delta t_{\text{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\lt f_{\text{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.

Figure 4 shows $S_{r}(\,f\,)$ for each frequency ratio, with black dashed lines indicating the estimated noise levels, h₀, calculated from values of $\sigma_{r}(\tau)$ [69]. Here it can be seen again that while $r_{[\text{Sr/Cs}]}$ and $r_{[\text{Yb}^{+}\text{/Cs}]}$ are well described by WFM noise, this is not true of $r_{[\text{Yb}^{+}\text{/Sr}]}$ for $f \geqslant (60\text{s})^{-1} = 1.67 \times 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\,)\propto f^{-1.5}$ behavior. Limits on oscillations in $\Delta\alpha/\alpha$ in the frequency range that exhibits WFM noise can be formed by integrating $S_{r[\text{Yb}^{+}\text{/Sr}]}(\,f\,)$ over the nominal frequency bin width ($\delta f_{[\text{Yb}^{+}\text{/Sr}]} = 1/T_{[\text{Yb}^{+}\text{/Sr}]} = 2.5 \times 10^{-6}$ Hz) to obtain the total power of each bin, then taking the square-root to obtain the one-sided ⁶ fractional amplitude spectrum, $A_{r[\text{Yb}^{+}\text{/Sr}]}(\,f\,)$.

Zoom In Zoom Out Reset image size

Dividing this spectrum by the magnitude of $\Delta K_{\alpha}^{\text{Yb}^{+}\text{/Sr}}$ yields the fractional amplitude spectrum of $\Delta \alpha / \alpha$ oscillations: $A_{\alpha}(\,f\,) = A_{r[\text{Yb}^{+}\text{/Sr}]}(\,f\,) / 6.01 \approx 7.6 \times 10^{-19}$. Similarly, we constrain oscillations in $\Delta\mu/\mu$ using $S_{r[\text{Sr/Cs}]}(\,f\,)$ ($\delta f_{[\text{Sr/Cs}]} = 1/T_{[\text{Sr/Cs}]} = 8.7\times 10^{-7}$ Hz) resulting in a fractional amplitude spectrum of µ as $A_{\mu}(\,f\,) \approx A_{[\text{Sr/Cs}]}(\,f\,) / 1.00$ $\approx 2.7 \times 10^{-16}$. Both fractional amplitude spectra $A_{\alpha}(\,f\,)$ and $A_{\mu}(\,f\,)$ are shown in figure 5 for the frequency ranges over which they may be used to constrain oscillations in α and µ.

Zoom In Zoom Out Reset image size

Constraints on oscillations in α and µ can be extracted from the fractional amplitude spectra $A_{\alpha}(\,f\,)$ and $A_{\mu}(\,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 by false alarm levels, $A_{X}(p \leqslant 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₀ 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_{\alpha}(p \leqslant 0.32) = 5.6\times 10^{-18}$ and $A_{\mu}(p \leqslant 0.32) = 1.3\times 10^{-15}$. More appropriate to particle physics would be the equivalents of a 5σ significance level, which were estimated using the Baluev method [81], yielding $A_{\alpha}(p \leqslant 3.5 \times 10^{-7}) \approx$ $8.9\times 10^{-18}$ and $A_{\mu}(p \leqslant 3.5\times 10^{-7}) \approx 2.1\times 10^{-15}$. ⁷ As shown in figure 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_{\mu}(\,f\,)$, and no spectral peaks exceeding the 5σ false alarm levels for either $A_{\mu}(\,f\,)$ or $A_{\alpha}(\,f\,)$, we have high confidence that all peaks in $A_{\alpha}(\,f\,)$ and $A_{\mu}(\,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_{\alpha}(\,f\,) \leqslant 5.6\times 10^{-18}$ and at the 5σ significance level of $A_{\alpha}(\,f\,) \lesssim 8.9\times10^{-18}$, over the frequency range $2.5 \times 10^{-6}$ Hz $\leqslant f \leqslant 1.7\times 10^{-2}$ Hz. We also place constraints on sinusoidal oscillations in µ at the 1σ significance level of $A_{\mu}(\,f\,) \leqslant 1.3\times 10^{-15}$ and at the 5σ significance level of $A_{\mu}(\,f\,) \lesssim 2.1\times 10^{-15}$, over the frequency range $8.7\times 10^{-7}$ Hz $\leqslant f \leqslant 8.3\times 10^{-4}$ Hz. A summary of these 1σ and 5σ constraints is shown in table 3.

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.

Amplitude constraint (1σ significance)	Amplitude constraint (5σ significance)	Parameter space
$A_{\alpha}(\,f\,) \leqslant 5.6\times 10^{-18}$	$A_{\alpha}(\,f\,) \lesssim 8.9\times 10^{-18}$	$2.5\times 10^{-6}$ Hz $\leqslant f \leqslant 1.7\times 10^{-2}$ Hz
$A_{\mu}(\,f\,) \leqslant 1.3\times 10^{-15}$	$A_{\mu}(\,f\,) \lesssim 2.1\times 10^{-15}$	$8.7\times 10^{-7}$ Hz $\leqslant f \leqslant 8.3\times 10^{-4}$ Hz

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.

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 equation (13), a transition frequency ν may be parameterized as [82]

$$\begin{align} \nu = (\text{const.})\;(cR_\infty)\;\alpha^{K_{\alpha}}(m_e/\Lambda_\textrm{QCD})^{K_\mu} (m_{q}/\Lambda_\textrm{QCD})^{K_q}, \end{align} \tag{ 16 }$$

where $K_\alpha, K_\mu$, and K_q are sensitivity coefficients characteristic of the transition ν and $m_q \equiv (m_u+m_d)/2$. The quark coefficient $K_q = k_q^{M_N} + k_q^{g_N}$ parameterizes changes in the nucleon mass $\delta M_N/M_N = k_q^{M_N}(\delta m_q/m_q)$ and the nuclear magnetic moment $\delta g_{N}/g_{N} = k_q^{g_N}(\delta m_q/m_q)$ in terms of m_q variations. The ratio of two frequencies $r = \nu_1/\nu_2$ is independent of the dimensionful constants, and varying equation (16) with respect to $\alpha, m_e/\Lambda_{\text{QCD}}$, and $m_q/\Lambda_{\text{QCD}}$ noting (11) gives

$$\begin{align} \frac{\delta r}{r} = \left(\Delta K_{\alpha}d_\gamma^{(n)} + \Delta K_\mu\left(d_{m_e}^{(n)} - d_g^{(n)}\right) + \Delta K_q\left(d_{q}^{(n)} - d_g^{(n)}\right)\right)(\kappa\phi)^{n}, \end{align} \tag{ 17 }$$

where $d^{(n)}_{q}$ is the mass-weighted mean-quark coupling. Using the published values of $K_\alpha^{\text{Yb}^+} = -5.95, \ K_\alpha^{\text{Sr}} = 0.06$, and $K_\alpha^{\text{Cs}} = 2.83$ from [52] and $K_\mu^{\text{Cs}} = 1, K_q^{\text{Cs}} \approx 0.07$ from [82] (and summarized in table 1) results in the fractional frequency ratios

$$\begin{align} &\left(\frac{\delta r}{r}\right)_{\text{Yb}^+/\text{Sr}} = -6.01d_\gamma^{(n)}(\kappa \phi)^n, \end{align} \tag{ 18 }$$

$$\begin{align} &\left(\frac{\delta r}{r}\right)_{\text{Sr}/\text{Cs}} = -\left(2.77d_\gamma^{(n)} + d_{m_e}^{(n)}-d_g^{(n)} + 0.07\left(d_{q}^{(n)} - d_g^{(n)}\right)\right)(\kappa \phi)^n. \end{align} \tag{ 19 }$$

As can be observed in figure 2, constraints from Yb⁺/Sr measurements enable a bound on the coefficient $d_\gamma^{(n)}$ far below what is possible from Sr/Cs measurements. It is therefore appropriate to neglect $d_\gamma^{(n)}$ in equation (19), which gives an effective coupling for the Sr/Cs measurements of

$$\begin{align} d_{\text{Sr/Cs}}^{(n)} \approx d_{m_e}^{(n)}-d_g^{(n)} + 0.07\left(d_{q}^{(n)} - d_g^{(n)}\right). \end{align} \tag{ 20 }$$

Comparing equations (18) and (19) with the data in figure 2, model-independent constraints can be placed on the magnitudes of the coupling strengths $d_\gamma^{(n)}$, $d_{\text{Sr/Cs}}^{(n)}$, and the instability of $\phi^n(t)$ over different timescales, $\sigma_{\phi^n}(\tau)$ given by

$$\begin{align} \kappa^n |d_\gamma^{(n)}| \sigma_{\phi^n}(\tau) \lesssim 2.3\times 10^{-16}/\sqrt{\tau /\textrm{s}}\;, \end{align} \tag{ 21 }$$

for timescales 60 s $\leqslant \tau \leqslant$ 30 000 s, and

$$\begin{align} \kappa^n |d_{\text{Sr/Cs}}^{(n)}| \sigma_{\phi^n}(\tau) \lesssim 1.6\times 10^{-13}/\sqrt{\tau /\textrm{s}}\;, \end{align} \tag{ 22 }$$

for timescales 600 s $\leqslant \tau \leqslant$ 80 000 s, where $\sigma_{\phi^n}(\tau)$ is the MDEV of φⁿ , defined in equation (15), and the contribution from $(d_{q}^{(n)} - d_g^{(n)})$ 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 $\phi^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 $\sigma_{\phi^n}(\tau)$ in equations (21) and (22) can be roughly interpreted as limits on the average magnitude of fluctuations in $\phi^n(t)$ between any two points in time t and $(t+\tau)$. For example, for two times separated by $\tau = 1~000$ s, the fluctuation $[\phi^n(t+\tau) - \phi^n(t)]$ should roughly satisfy $\kappa^n |d_\gamma^{(n)}| [\phi^n(t+\tau) -$ $\phi^n(t)] \lt7\times 10^{-18}.$

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 ULDM in the following section.

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}\;\text{eV} \lesssim m_{\phi} \lesssim 1\;\text{eV}$ is known as ULDM, and in recent years significant efforts have been focused on detecting ULDM through apparent oscillations of fundamental constants. The upper bound $m_\phi \approx 1$  eV occurs when the number density n of bosons in the reduced de Broglie volume $(\lambda/2\pi)^3$ satisfies $n(\lambda/2\pi)^3 \gg 1$, resulting in a macroscopic phase-space occupation that exhibits Bose–Einstein condensation. Wavelengths $\lambda \sim \mathcal{O}(\text{kpc})$ span distances comparable to the smallest dwarf galaxy halos and imply a lower bound $m_\phi \approx 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_\textrm{vir}\sim 10^{-3}$ is the virial speed (in natural units) of dark matter with isotropic distribution $\langle \vec{v}_{\text{vir}}\rangle = \vec{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 $|\vec{k}| \approx m_\phi v_{\text{vir}} \ll$ $m_\phi$ when neglecting subdominant corrections [85]. In the dark-matter rest frame, oscillations are controlled by the rest mass $m_\phi = 2{\pi}f_\phi$ where f_φ is the Compton frequency in natural units. Oscillations are coherent in time for $\tau_c \sim 4\pi/(m_\phi v_\textrm{vir}^2) \gtrsim 10^6 \ T_{c} \gg T_\textrm{data}$ where the oscillation timescale $T_c = 1/f_\phi$ greatly exceeds the experimental timescale $T_\textrm{data}$. As the coherence length $\lambda \sim 2\pi/(m_\phi v_\textrm{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 constant-amplitude waveform oscillating at the underlying particle Compton frequency, up to small velocity corrections ${\sim}\mathcal{O}(v^2_{\text{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 $\Gamma = 0$. The field has an amplitude $\phi_0 = \sqrt{2\langle \rho_{\phi}\rangle }/m_\phi$, resulting from time-averaging the energy density $\rho_{\phi} = (\dot\phi^2 + m_\phi^2\phi^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

$$\begin{align} \phi(t) \approx \frac{\sqrt{2\langle \rho_{\phi}\rangle }}{m_\phi}\cos(m_\phi t). \end{align} \tag{ 23 }$$

Typically it is assumed that the scalar field comprises the entirety of the dark-matter density inferred at the solar galactocentric radius $R_0\simeq 8$  kpc, or the 'local density' $\rho_\textrm{DM}(R_0)\approx 0.3\;\text{GeV}\,\text{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 $\rho_{\text{DM}}\lesssim (10^4$-$10^6)\times \rho_{\text{DM}}(R_0)$. The immediate implications are on the sensitivity to the couplings in equations (18) and (19), since for a generic coefficient (assuming the field saturates the total density) the sensitivity to $d_j^{(n)} \propto 1/\rho_{\text{DM}}^{n/2}$ and at best linear (quadratic) constraints would scale downward by ${\sim}10^3$ (${\sim}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. $\Gamma\neq 0$), and the density contribution would be spread across the components $\rho_{\text{DM}} = \sum_i \langle \rho_{\phi_i}\rangle$. The extracted limits would be weaker and scaled by $\sqrt{\langle \rho_{\phi_i}\rangle /\rho_\textrm{DM}}$.

5.1. Scalar couplings

Previous experiments monitoring electronic and microwave frequencies have resulted in constraints on the scalar couplings $d_\gamma^{(n)}, d_{m_f}^{(n)}$, and $d_g^{(n)}$ assuming dark matter [15, 21–24, 74, 90–93]. In the present work, constraints may be extracted by using the amplitude spectra $A_X(\,f\,) = A_r(\,f\,)/|\Delta K|$ for frequency bin f and relevant sensitivity coefficient $\Delta K$ from figure 5. Identifying $f = f_\phi = m_\phi/2\pi$ for linear couplings and $f = 2f_\phi = m_\phi/\pi$ for quadratic couplings and noting equation (17) gives

$$\begin{align} d_j^{(n)} = \zeta_{f}^{(n)} A_X(\,f\,), \end{align} \tag{ 24 }$$

where $\zeta_{f}^{(1)} \approx 2.0\times 10^{16}f$/Hz, $\zeta_{f}^{(2)} \approx (0.2\zeta_{f}^{(1)})^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_\phi\rightarrow f_\phi +$ $m_\phi v_\textrm{vir.}^2(4\pi)^{-1}$, have been included in $\zeta_{f}^{(1)}$. As a result, the sum over distinct field modes for measurement times $T\ll \tau_c$ reduces the sensitivity to the linear n = 1 coefficients by a stochastic factor ≈ 3 [94]. ⁸ The data encompasses frequencies $10^{-6}\;\text{Hz} \lesssim f \lesssim 2\times 10^{-2}\;\text{Hz}$, corresponding to masses $4\times 10^{-21}\;\text{eV} \lesssim m_\phi \lesssim$ $8\times 10^{-17}\;\text{eV}$.

The results for linear and quadratic scalar couplings discussed in section 2 are presented for the case of ULDM in figures 6–9. Note that some works assume the slightly larger estimate $\rho_{\text{DM}}\approx 0.4\;\text{GeV}\,\text{cm}^{-3}$, which for purposes of comparison here amounts to a negligible difference. For all figures, the right axes compare $\kappa^n d_j^{(n)} = 1/\Lambda_j^n$ for a generic dimensionful scale $\Lambda_j$. We choose the scale to be identified with $\Lambda_j, \Lambda_j^{^{\prime}}$ for linear and quadratic couplings, respectively. For the coupling $d_\gamma^{(1)}$ in figure 6, a new exclusion region up to roughly an order of magnitude improvement over previously published work for $10^{-20}\;\text{eV} \lesssim m_\phi \lesssim 10^{-17}\;\text{eV}$ is observed. This improved sensitivity is mostly explained by the difference in K factors between, e.g. Rb/Cs where $|\Delta K^{\text{Rb/Cs}}_\alpha| \approx 0.5$ is roughly an order of magnitude smaller than $|\Delta K^{\text{Yb}^+\text{/Sr}}_\alpha| \approx 6$. The extracted bounds also essentially surpass Eöt-Wash [95, 96] and MICROSCOPE [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_\gamma^{(1)}$ in this mass range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the quadratic coupling $d_\gamma^{(2)}$, a similar trend with respect to clock studies [90, 91] using Rb/Cs [21] and ¹⁶⁴Dy/¹⁶²Dy spectroscopy [74] persists and using EP-test results to extract bounds on quadratic couplings [99] shows greater sensitivity ${\gtrsim }4\times 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 ${\approx}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 NMIJ using Yb/Cs clock comparisons [24]. Note that EP tests do not rely on assumptions of the amplitude φ₀ and thus the contribution of the scalar field to the dark matter abundance. Similarly, though BBN constraints use $\langle \rho_\phi\rangle = \rho_{\text{DM}}$, the field is non-oscillating with constant φ₀ for $m_\phi \ll 10^{-16}$  eV so coherence considerations are irrelevant.

Limits on linear and quadratic $d_{m_e}^{(n)} - d_g^{(n)}$ and $d_{q}^{(n)} - d_g^{(n)}$ couplings are presented in figures 7 and 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_\phi \approx 10^{-23}$  eV. Note that Rb/Cs has no sensitivity to $d_{m_e}^{(n)} - d_g^{(n)}$. Regarding the quadratic constraints, Sr/Cs data probes a new region for clocks in the $d_{m_e}^{(2)} - d_g^{(2)}$ 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_{q}^{(2)} - d_g^{(2)}$.

The final scalar constraints are extracted on the parameter A from the scalar-Higgs interaction and are presented in figure 9. The simplest (n = 1) linear couplings have garnered theoretical attention since they can emerge from the technically natural operator $\mathcal{L}_{\phi H} = -A\phi H^\dagger H$ for Higgs doublet H [100]. The sensitivity coefficient is

$$\begin{align} K_H = \frac{\alpha}{2\pi}K_\alpha - (1-b)K_\mu - 1.05(1-b)K_q, \end{align} \tag{ 25 }$$

where $b\sim 0.2-0.5$ is a dimensionless factor in the Higgs-nucleon Yukawa coupling $g_{hNN} = 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 \approx 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 [90]. Using the Sr/Cs spectrum and noting $\kappa d_H \leftrightarrow A/m_h^2$ for Higgs mass $m_h \approx 125$  GeV in (24) produces the limit. For the ratios considered here, the bulk of the sensitivity comes from $(1-b)K_\mu$ as the 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 microwave-microwave 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} \lesssim m_\phi \lesssim 10^{-20}$. Note again that the current Rb/Cs limits do not include the stochastic degradation factor ${\approx}3$, which would shift the displayed curve upwards accordingly.

5.2. 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:

$$\begin{align} \mathcal{L}_a = \frac{g_s^2}{32{\pi}^2}\frac{a}{f_a} G_{\mu\nu}^b\widetilde{G}^{b\mu\nu}, \end{align} \tag{ 26 }$$

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

$$\begin{align} \frac{\delta m_{\pi}^2 }{m_{\pi}^2} = - \frac{m_u m_d}{2(m_u + m_d)^2} \theta^2, \end{align} \tag{ 27 }$$

where $m_u, m_d$ are the up and down quarks and $\theta = a/f_a$. The nucleon mass $m_N(\theta)$ and the nucleon g-factor $g_N(\theta)$ (and hence nuclear g-factor $g(\theta)$) have an inherent dependence on $m_{\pi}(\theta)$ as parameterized by chiral perturbation theory [61, 102–104]. Similar to quadratic scalar couplings in the ultralight range, the oscillatory component of θ is given by

$$\begin{align} \theta^2 (t) = \frac{\rho_{\text{DM}}}{m_a^2 \,f_a^{\;2}}\cos(2m_a t), \end{align} \tag{ 28 }$$

where here $\rho_\textrm{DM} = 0.4$  GeV/cm³, $m_a \not\propto f_a^{-1}$ in general and assuming a saturates the local density. The variation of a microwave transition frequency (13) may be expressed as variations in $g(\theta)$ and m_p using equations (27) and (28), giving [86]

$$\begin{align} \frac{\delta \nu_\textrm{MW}}{\nu_\textrm{MW}} \approx -0.11\left[\frac{\partial \ln g(\theta)}{\partial \ln m_{\pi}^2} - \frac{\partial \ln m_p}{\partial \ln m_{\pi}^2} \right]\frac{\rho_{\text{DM}}}{m_a^2 \,f_a^{\;2}}\cos(2m_a t). \end{align} \tag{ 29 }$$

Comparing with another microwave or optical standard and identifying $m_a = {\pi}f$ for signal frequency f from an amplitude spectrum yields the relation

$$\begin{align} \frac{1}{f_a\cdot\text{GeV}^{-1}} = 10^{-10}\sqrt{\frac{m_{15}^2}{c_r\cdot 10^{-15}}\bigg|\frac{\delta r}{r}\bigg|}, \end{align} \tag{ 30 }$$

where $m_{15} \equiv m_a/(10^{-15}\;\text{eV})$ and c_r is a constant. In [101] microwave-microwave (Rb/Cs, $c_r\approx 10^{-1}$) and microwave-optical (H/Si, $c_r \approx 1$) comparisons are studied. One may also consider Sr/Cs, where due to the weak dependence of Sr on nuclear quantities $|\delta r/ r| \approx \delta \nu_\text{Cs}/\nu_\text{Cs}$ and $c_r\approx 2\times 10^{-2}$, including effects of stochasticity. From equation (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 figure 10.

Zoom In Zoom Out Reset image size

Perhaps most compelling are the high-frequency range constraints achievable on the axion-like coupling in figure 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 were operated without compensating density-dependent effects, cycle times of ${\approx}5$ s could be achieved, allowing Fourier frequencies up to ${\approx}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.

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 (⁸⁷Sr, ¹⁷¹Yb⁺) and microwave (¹³³Cs) frequencies, we constrain the instability of fractional changes in α to be $\sigma_{(\Delta\alpha/\alpha)}(\tau)~\leqslant 2.3\times 10^{-16}/\sqrt{\tau/\textrm{s}}$ for averaging times 60 s $\lt \tau \lt$ 30 000 s, and we constrain the instability of fractional changes in µ to be $\sigma_{(\Delta\mu/\mu)}(\tau) \leqslant 1.6\times 10^{-13}/\sqrt{\tau/\textrm{s}}$ for averaging times 600 s $\lt \tau \lt 80~000$ s. The theoretical framework then allows us to place constraints on combinations of a new scalar field $\phi^n(t)$ and the different interaction strengths $d_j^{(n)}$ coupling that field to other particles: photons, electrons, quarks and gluons. These constraints are independent of the underlying physics and the functional form of $\phi^n(t)$.

As an example of a specific model, we studied ULDM couplings to matter and presented new constraints on low-dimension dilaton-like operators. The limits on $d_\gamma^{(1)}$ from Yb⁺/Sr data exclude a new region of parameter space for masses $10^{-20}\textrm{eV} \lesssim m_\phi \lesssim 10^{-17}\textrm{eV}$, as shown in figure 6. We also refer the reader to new experimental results [25] using Yb⁺ and Sr clocks, which claim even tighter constraints on $d_\gamma^{(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].

This work was carried out as part of the QSNET project within the Quantum Technologies for Fundamental Physics Programme, with support from the Science and Technology Facilities Council (Grant Numbers ST/T00102X/1, ST/T00598X/1 and ST/T006048/1). We also acknowledge the support of the UK government Department for Science, Innovation and Technology through the UK National Quantum Technologies Programme. We thank Abhishek Banerjee and Tejas Deshpande for providing exclusion-region results from previous experiments in figures 6 and 7, respectively, and Takumi Kobayashi for providing Yb/Cs data. X C would like to dedicate this work to the memory of Professor Harald Fritzsch who passed away on the 16th of August 2022.

The data that support the findings of this study are available upon request from the authors.

The scalar field φ can couple to quarks or gluons linearly (n = 1)

$$\begin{align} {\cal L} = \kappa \phi \left (\frac{d^{(1)}_g}{4} G_{\mu\nu}G^{\mu\nu} -d^{(1)}_{m_f} m_f \,\bar \psi_f \,\psi_f \right ), \end{align} \tag{ 31 }$$

with $\kappa = \sqrt{4 {\pi}G}$, but we could also have quadratic couplings (n = 2)

$$\begin{align} {\cal L} = \kappa^2 \phi^2 \left (\frac{d^{(2)}_g}{4} G_{\mu\nu}G^{\mu\nu} -d^{(2)}_{m_f} m_f \,\bar \psi_f \,\psi_f \right). \end{align} \tag{ 32 }$$

One can disentangle the different parameters φ₀ and d_j . By considering experiments sensitive to α, $\mu = m_e/m_p$ or α_s , one could in principle measure φ₀ 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 $\Lambda_\textrm{QCD} = \mu_r \exp(2 \pi/(\alpha_s(\mu_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 $\Lambda_\textrm{QCD}$. Using the renormalization group equation for α_s we find

$$\begin{align} \frac{\dot m_p}{m_p} \approx\frac{\dot \Lambda_\textrm{QCD}}{\Lambda_\textrm{QCD}} = -\frac{2\pi}{\beta} \frac{\dot \alpha_s}{\alpha_s^2}, \end{align} \tag{ 33 }$$

where in the linear case, we have in the underdamped regime

$$\begin{align} \frac{\dot \alpha_s}{\alpha_s} = -\frac{\kappa d^{(1)}_g \phi_0}{2} \exp{\left(-\frac{\Gamma t}{2}\right)} \bigg( \Gamma \cos\left( \theta - t \omega_{d}\right) + 2 \omega_{d} \sin\left( \theta - t \omega_{d} \right) \bigg) \end{align} \tag{ 34 }$$

and similarly again in the underdamped regime

$$\begin{align} \frac{\dot \alpha_s}{\alpha_s} = -\frac{\kappa^2 d^{(2)}_g \phi_0^2}{2} \exp{\left(-\frac{\Gamma t}{2}\right)} \bigg( \Gamma \cos\left( \theta - t \omega_{d}\right) + 2 \omega_{d} \sin\left( \theta - t \omega_{d}\right) \bigg) \end{align} \tag{ 35 }$$

for the quadratic case.