Searching for Scalar Dark Matter via Coupling to Fundamental Constants with Photonic, Atomic and Mechanical Oscillators

We present a way to search for light scalar dark matter (DM), seeking to exploit putative coupling between dark matter scalar fields and fundamental constants, by searching for frequency modulations in direct comparisons between frequency stable oscillators. Specifically we compare a Cryogenic Sapphire Oscillator (CSO), Hydrogen Maser (HM) atomic oscillator and a bulk acoustic wave quartz oscillator (OCXO). This work includes the first calculation of the dependence of acoustic oscillators on variations of the fundamental constants, and demonstration that they can be a sensitive tool for scalar DM experiments. Results are presented based on 16 days of data in comparisons between the HM and OCXO, and 2 days of comparison between the OCXO and CSO. No evidence of oscillating fundamental constants consistent with a coupling to scalar dark matter is found, and instead limits on the strength of these couplings as a function of the dark matter mass are determined. We constrain the dimensionless coupling constant $d_e$ and combination $|d_{m_e}-d_g|$ across the mass band $4.4\times10^{-19}\lesssim m_\varphi \lesssim 6.8\times10^{-14}\:\text{eV} c^{-2}$, with most sensitive limits $d_e\gtrsim1.59\times10^{-1}$, $|d_{m_e}-dg|\gtrsim6.97\times10^{-1}$. Notably, these limits do not rely on Maximum Reach Analysis (MRA), instead employing the more general coefficient separation technique. This experiment paves the way for future, highly sensitive experiments based on state-of-the-art acoustic oscillators, and we show that these limits can be competitive with the best current MRA-based exclusion limits.

If such DM-SM couplings exist, oscillations in fundamental constants will frequency modulate various types of clocks and oscillators, measurable via clock comparison experiments if the oscillators or clocks exhibit different dependences on the fundamental constants. Clock comparisons present a powerful tool for searching for such variations, owing to the high stability of modern frequency standards [28,[33][34][35].
Confounding these experiments is the fact that the DM particle mass, as well as the strength of its coupling to the SM, is unknown, and only weakly constrained [36]. As a result there is a large parameter space to search, and experiments are required which search for a range of DM particles masses, corresponding to frequency standards which are stable over different time-scales.
In this work we search for variations in the fundamental constants caused by a massive scalar field constituting the local DM halo. We present limits on the possible variation of linear combinations of fundamental constants, thus placing experimental constraints on the coupling of a DM scalar field to the SM. In particular, we contribute new constraints to the relatively unexplored higher mass area of the putative scalar field's parameter space; 4.4 × 10 −19 m ϕ 6.8 × 10 −14 eVc −2 . We achieve this by monitoring frequency fluctuations of a quartz crystal bulk acoustic wave oscillator (OCXO), compared against both cryogenic sapphire oscillator (CSO) and hydrogen maser reference clocks (HM), each of which exhibits a different dependence on the fundamental constants. We consider the model of Darmour and Donoghue [37] as implemented in Ref. [38], where ϕ is a dimensionless, massive scalar field with a quadratic self-interaction potential V (ϕ) = 2 c 2 2 m 2 ϕ ϕ 2 in which the normalisation has been chosen so that m ϕ has dimensions of mass. This model considers ϕ modifications to terms in the effective action that describes the physics of ground state nuclei. At appropriately low energy scales of ∼ 1 GeV this effective action will only contain the electron e, the up quark u and the down quark d as real particles with interactions mediated by electromagnetic (A µ ) and gluonic (A A µ ) fields. Weak interactions and heavy quarks are integrated out at this scale, while it is argued in Ref. [37] that EEP violation effects linked to the strange quark are relatively small, thus it is ignored here.
Each of the five terms described by this effective action can then couple to ϕ. For this work, as is common, we consider a linear couplings for these terms, thus giving a Lagrangian density for scalar field-SM interactions L int as per equation (12) of Ref. [37]. Linear couplings are often considered to be the most "simple" and therefore most "compelling" couplings to the SM, taking the form where F µν is the electromagnetic Faraday tensor, µ 0 the magnetic permeability, F A µν is the gluon strength tensor, g 3 the QCD gauge coupling, β g is the β function for the running of g 3 , m i is the mass of the fermions, γ m i is then the anomalous dimension giving the energy running of the masses of the QCD coupled fermions, and ψ i denotes the fermion spinors.
The constants of interest d j for j = m u , m d , m e , g, e are dimensionless coupling constants that parametrise the scalar field coupling to the SM matter fields, defining the strength of the interaction in a corresponding SM sector, and equivalently, the magnitude of any ϕ dependent oscillations in the corresponding fundamental constants.
The introduction of these coupling constants into the interaction Lagrangian density will modify each term such that corresponding fundamental constants will display dependencies on ϕ of the following forms: Where α is the fine structure constant, m i denotes fermion mass (electron (e), up (u) or down (d) quarks), and Λ QCD represents the QCD mass scale. We also note that the mean quark massm = (m u + m d )/2 displays a similar dependency on ϕ: The periodic evolution of the scalar field is given by the description in Ref. [38] which borrows from the cosmological string-theory dilaton model of Ref. [17]; scalar field and matter. For short time periods (t 1/H) this term can be considered a constant. Following the description presented in Ref. [38], we identify the scalar field as DM with energy density, We see that for typical values for the local DM density of ρ DM = 0.45 GeV/cm 3 [39] the amplitude of scalar field oscillations would range from 1.84 × 10 −13 < ϕ 0 < 1.84 × 10 −17 for 1 mHz < ω ϕ /2π < 10 Hz. It is thus possible to experimentally probe the coupling of a DM scalar field to matter, by searching for ϕ dependent variations in dimensionless ratios of fundamental constants. This can and has been achieved by making comparison measurements between differing frequency standards whose frequency ratio is dependent on some combination of the dimensionless factors: m u /Λ QCD , m d /Λ QCD , m e /Λ QCD and α. An experiment of this nature is thus able to probe scalar field couplings linear in d e and d m i −d g , by considering equations (2) to (7), and the local dark matter parameters.
It is often difficult to extract an exclusion limit on individual coupling parameters, owing to these combinations. There are two common methods for extracting information on individual parameters. One is known as 'Maximum Reach Analysis' (MRA) in which a series of models are considered where it is assumed that the field in question only couples to one SM sector in each model, thus excluding one non-zero coupling parameter at a time. The other method is to take multiple sets of data with differing dependencies on the constants, so that linear combinations of different parameters can be separated. This is known as 'coefficient separation' and is the technique that we employ in this work. By introducing a mechanical resonator, which depends in a different way on the fundamental constants to various photonic and atomic oscillators, we are able to separate coefficients which have not been separated before.
In one of our experiments we analysed fluctuations of the phase difference δφ 21 between a 10 MHz NEL Frequency Controls Inc. ultra low phase noise OCXO and a 10 MHz signal synthesised from a microwave CSO [40]. Fig. 1   The spectral density of the δφ 21 was inferred from the spectrum of PLL correction voltage δu corr (Fig. 1) via the following relationship where F denotes the Fourier frequency, γ is the loop gain, df/du is the frequency-voltage tuning coefficient of the quartz oscillator, S δφ is the spectral density of phase difference fluctuations and S n/f φ is the spectral density of the PLL phase noise floor.
All parameters of the PLL in Eq. (8) can be determined experimentally, thus monitoring the PLL voltage signal δu allows us to measure phase noise variations ϕ 21 that are synchronous with variations in the beat frequency δf 21 = δf 2 − δf 1 of two different frequency standards.
We now consider how a variation in the fundamental constants will affect the beat frequency between each combination of these frequency standards. To do this we write the dependence of each standard's mode frequencies as a combination of the dimensionless constants α, m q /Λ QCD and m e /m p ∝ m e /Λ QCD [22]. Here we are assuming the constants stated above are functions of ϕ as per eqs. (2) to (4), while any other factors contained in the mode frequency dependence are true constants [41].
The derivation of the dependence of the quartz oscillator mode on the fundamental constants is discussed in the appendices and is found to be given by The dependencies of both the CSO and Maser frequencies are given in Appendix A of Ref. [41],as the dependence of the CSO crystal's permittivity on α can be ignored at the frequencies of interest [? ]. We do not consider here the sensitivity to quark mass through the spin g factor of the HM transition induced by small QCD corrections [? ], as its coefficient is more than an order of magnitude smaller than the other coefficients, and inclusion would cause the other coefficients to lose independence. In future work, to perform the technique of coefficient separation whilst taking into account these corrections, we would require an additional atomic oscillator with a different sensitivity to the quark mass.
By normalizing variations in beat frequency with respect to the shared carrier frequency f 0 = 10 MHz, we have, for the quartz against the CSO: Where we have used eqs.
We obtained two initial datasets as per the procedure outlined above. analysable Fourier frequencies is then determined by the sampling rate of the measurement apparatus and the total integration time. In addition to these two initial datasets, further measurements were made at later times over shorter periods in order to provide further complementary results. The CSO-OCXO and HM-OCXO experiments were sampled again, at a higher rate (33 Hz) for 12 hours in order to exclude large noise sources in the initial data, as well as generating fractional frequency noise data at higher frequencies. This data will be subject to the same DM search analysis in the proceeding section, giving less stringent but complementary results. Finally, a further OCXO-OCXO control experiment which displays zero DM sensitivity was run for 12 hours, in order to characterise spurious systematic noise sources in the main data.
The goal of this analysis is to determine a limit corresponding to the weakest possible scalar field-standard model coupling strength that can be confidently excluded in the case of no detection. The general procedure is as follows. The power spectral density of fractional frequency noise, S y , is searched for large deviations from the mean value at a range of Fourier frequencies, which would correspond to signals consistent with dark matter. A threshold "cut" value is chosen, and above this cut value any deviations from the mean are considered dark matter candidates. If all such signals can be excluded as either spurious or systematic noise, they are excluded as dark matter candidates, no detection is reported, and we move to derive exclusion limits. A signal size is determined via simulation, which corresponds to the minimum size of a dark matter signal, and thus DM-SM coupling, which we would expect to pass the chosen threshold cut value with 95% confidence. Given we have excluded all signals above the cut, this simulation-determined signal strength corresponds to our 95% confidence exclusion limit, as we would have expected a signal of this size to remain in the data, survive the cut, and not be excluded, if it were present.
Through the utilised exclusion methods discussed in the appendices, we can exclude all peaks in our data as due to spurious or systematic noise. Furthermore, we are confident that should a signal containing DM characteristics arise; it would fail to be excluded by this analysis, and a strong claim of DM detection could be made.
An example of the 95% confidence exclusion limits in terms of frequency fluctuation signal strength for the CSO-OCXO experiment are as shown in Fig. 2, similar limits were also computed for the HM-OCXO experiment. The excluded amplitude of frequency variation is given by the square root of these limits (when converted back from fractional to absolute frequency), which can then be substituted into Eqs. (12) to give experimental exclusion limits on two different linear combinations of d e , d me and d g (a different set of coupling parameters for each type of oscillator comparison). We then utilise coefficient separation to provide further limits, effectively solving linear equations for |d e | and |d me − d g |. Fig. 3 presents these final exclusion limits derived from both the initial sets of data, and the later higher frequency data, displayed as one combined limit.
We note that previous literature limits in this region [25,43] have been approximated by applying MRA to experimental data from EEP and WEP tests. MRA is performed by effectively assuming experimental sensitivity to only one coupling parameter, setting all other parameters to zero. While this method has produced the most competitive limits to date; we note that it is an idealistic estimate, and underpinned by the inherent assumption that any cancellation between multiple parameters is unlikely [44,45]. In deriving our exclusion limits we have considered no case were the coupling to specific SM sectors is temporarily ignored, making them more general. Although such results are obviously inherently less sensitive, they are of slightly different and complimentary significance to those produced by MRA. as the Red trace. This is a combined limit derived from both the initial and higher frequency data.
The best known limits in the region, as approximated in Ref. [43] [29]. We also further present a projected estimate in Green, for a potential future experiment that would utilise new low phase noise oscillators over a 5 year period.
Recent developments in CSO and OCXO oscillator technologies are giving rise to a new wave of low phase and frequency noise oscillators, far superior to the devices used in this experiment. We have included projected limits in Fig. 3, as well as frequency stability performance in Fig. 1, for a similar hypothetical experiment that searches for variations in CSO-HM and CSO-OCXO frequency differences, using current best-case noise characteristics for such devices from Refs. [28,46,47]. A 5 year experimental run time in this scenario will achieve general sensitivity limits via coefficient separation that are comparable with those produced by the approximations of MRA, however we note that as sensitivity scales with T 1/4 the contribution to the improvement of these limits due to longer run times is inferior to that associated with the use of oscillators with better frequency stability. The process for determining these exclusion limits, along with the assumptions made about oscillator performance, are outlined in the appendices. Further improvements to sensitivity could be made by operating the quartz oscillator in a cryogenic environment, where these oscillators see a boost in quality factor of several orders of magnitude [28,48,49]. However, due to several practical challenges [50,51], such a system is yet to be experimentally realized.
In conclusion, we present exclusion limits on scalar dark matter coupling to the standard model over several orders of magnitude in dark matter particle mass, based on frequency comparisons of stable oscillators. We exclude parameter space for the coupling constant  where v is the speed of sound in the material, L is the relevant length parameter of the resonator, and n is a constant. We will now consider how each of these parameters depends on fundamental constants.
The lengths of solids are (to first order, ignoring small relativistic corrections) proportional to the Bohr radius, a 0 , the characteristic scale for the size of atoms [52,53]. Thus for a solid, such as a BAW resonator, The speed of sound in a medium is given by Here K is the relevant elastic modulus for the type of sound wave and material, and ρ is the density. For the following we consider K to be the Bulk modulus. Given all elastic moduli in solids depend on the balance of the same electrostatic attractive and repulsive forces between atoms, we assume that they depend on the fundamental constants in the same way, and thus sound waves depend in the same way regardless of the modulus chosen.
We can rewrite the velocity as by expressing density in terms of L and m, the mass. Given that the majority of the mass of a BAW is baryons, the variation of m with constants can be considered proportional to variation of m p , assuming that it is some large multiple of the proton mass.
According to [54], for solids with electrostatic inter-atomic bonding, such as quartz where A is a (dimensionful) constant relating to the attractive forces between atoms, and r 0 is the inter-atomic spacing. We assume that, similar to the physical size of the material, the inter-atomic spacing r 0 depends (to first order, ignoring small relativistic corrections)  N]), and the square of the characteristic inter-atomic length scale. F C is given by where z i is the atomic number of element i, and r is the relevant length scale. For the purpose of the fundamental constant dependence it does not matter which length scale we consider, since we assume that they all vary with fundamental constants in the same way.
This ultimately leads us to This is also the dependence of the Coulomb constant on the fundamental constants, which supports this analysis. This, combined with (A2) and (A5), yields Substituting (A4) into (A1) gives us Finally, substituting (A8) and (A2), in terms of fundamental constants (assuming n is a true constant), we arrive at

Appendix B: Data Analysis Considerations
As discussed in the main body of this work, in order to achieve our goal of determining the limit corresponding to the weakest coupling strength that can be excluded in the case of no detection, we employed a modified version of the power bin search analysis method presented by Daw in Ref. [55]. Here we present a discussion of the method and the main considerations that were made when applying it to this work. In performing the power bin search method, data acquired over a broadband frequencybinned range is searched for characteristics that correspond to the effect expected to be induced by a DM signal. In the case of [55] as well as that of our work, the statistic that defines the search is either power excess above the mean level in a single frequency bin, or the total power excess given by the summation of a small number of consecutive bins.
Threshold cuts can then be made to the data to identify candidate bins that contain large power excess potentially due to DM signals. These remaining DM candidates are then systematically eliminated by identifying spurious non-DM related noise signals in the frame of the experiment that co-align at candidate frequencies. If all candidates above the chosen threshold power level can be removed in this way, the threshold can be used to compute a statistical exclusion limit on the DM coupling strength.
The motivation behind choosing between single or multiple bin search channels is derived from consideration of the DM signal's linewidth (if the DM signal is distributed over a large enough band it will deposit power into multiple neighbouring frequency bins), hence the most optimal search channel is defined by our assumptions of the shape of a DM signal. for the oscillators presented in Refs. [28,46,47], that we plan to use for future experiments.
In this work we are considering a scalar field that constituted the entirety of the local DM halo, thus we assume that this field has been thermalised by its motion in the galactic gravitational potential and therefore any signal, when converted to frequency space, will be defined by a Maxwell Boltzmann distribution of linewidth 10 −6 × ω ϕ [39]. This then defines two regimes in the frequency band spanned by our data; the first occurs at low frequencies (DM masses) where the DM linewidth is small and thus its frequency distribution will be completely contained by a single bin. For the data in this regime we conduct a single bin search as described above. The second regime occurs at frequencies where the DM signal's linewidth becomes larger than the width of each frequency bin. For the data in this regime we conduct a 3-bin search, as this provides the most optimal matching of signal linewidth to bin width across this regime of the data set.
We also acknowledge that recent theoretical developments [56] suggest that for integration times less than the coherence time of the DM signal (the first regime in our case), the DM signal exhibits a stochastic fluctuating amplitude, instead of the assumed fixed value ϕ 0 . For this work, as in the other previous works that we have presented, we assume a deterministic approach where the scalar field amplitude is fixed to an RMS value ϕ 0 . In order to provide exclusion limits across the entirety of our dataset, we must make the claim that no spurious DM-like signals were identified above a threshold cut. In order to support this claim we characterised all large spurious signals seen in the data as non-DM noise by one of three different strategies. The first strategy was to utilise the DM insensitive OCXO-OCXO data to find any large signal peaks at coincident frequencies, thus associating these signals as systematic background noise. Similarly, large signals seen in the 0.1-1 Hz frequency band of the initial data could be characterised as noise by looking at the DM sensitive data from the later higher frequency run, in which these large signals do not reoccur. The final method for characterisation is to measure the bandwidth of such signals and compare the measurement to the expected bandwidth of a DM induced effect. As discussed in the preceding section, a DM associated signal would exhibit a bandwidth of 10 −6 ×f m , which for our data corresponds to a width of at most 3 consecutive frequency bins. This allows us to exclude all signals of width > 3 bins as non-DM associated noise.
Using a combination of these three strategies, all spurious signals of large frequency noise excess were able to be characterised as non-DM related noise sources, hence we can claim zero DM detection events and move to determine exclusion limits. In the regions of frequency space where these large signals exist, we simply raise our limits accordingly as to show reduced sensitivity to DM signals in these areas.
Should we see any signal that cannot be exclude via the above methods, re-scans would be acquired to confirm its persistence and further higher resolution scans would follow to confirm the signal spectral shape. We are thus confident that should we see a spurious peak of narrow width and large enough magnitude (such that it cannot be excluded by the presented methods) we would be able to make a strong claim of a dark matter signal. We also note that this is a similar procedure for any clock comparison dark matter experiment, such as those referenced in the main text.
We utilised Monte Carlo (MC) simulations to convert our threshold frequency noise cut into a 95% confidence limit on the magnitude of frequency shifts induced by a DM scalar field. The simulation's algorithm constructs a multiple bin wide 'window' by superimposing a synthesised DM signal injected at some random frequency with a Gaussian-like noise distribution defined by the raw data's statistics. The injected signal's power is then varied until it passes the threshold cut in 95% of simulations. As we interpret our data as containing zero detection events, we define the threshold cut as the power level given by the bin with the highest signal-to-noise ratio (SNR) in each window of the real data, this defines the minimum DM SNR we can exclude. The ensemble of power levels determined by this MC simulation for each window then gives the confidence limit on excluded frequency noise power.
In order to estimate sensitivity for a future, improved experiment, we begin by collecting and constructing models of frequency noise in cutting edge oscillators. We take the results presented in Ref. [28] for an estimate of Quartz BAW phase noise, Ref. [46] for an improved CSO estimate, and the values quoted in the technical manual Ref. [47] for an optimistic HM estimate. Using the gathered frequency noise models we then compute fractional frequency noise power for each oscillator, and add the noise for each set of compared oscillators in quadrature. This gives an estimate of the fractional frequency noise power for a two resonator system, as used in this work.
We consider employing a dynamic bin search method (in which we have a varying integration time) on 5 years worth of frequency noise measurements. In order to estimate the smallest fractional frequency noise power which we could exclude from such data, we employ the well known Dicke Radiometer equation where S y 2 is the spectral density of fractional frequency fluctuations in units of Hz −1 with corresponding spectrum of fractional frequency fluctuations P y , T is the total measurement time, and BW is the spectral bin width. Thus, setting an SN R goal of 1 allows us to estimate the lowest fractional frequency noise power we could exclude (P yexc ) in a two resonator experiment. Following the procedure in the main body, taking the square root of this limit gives the excluded amplitude of frequency variation which can then be substituted into an equation similar to that of Eqs. (13) of the main text to give limits on a linear combination coupling constants. We found that the best limits on |d e | and |d me − d g | would be given by performing coefficient separation on a CSO-HM, CSO-BAW combination, the results of such a projection are presented in the Fig. 4 in the main text.