Here we demonstrate the proposed Bayesian nonparametric approach based on the
corrected multivariate likelihood on a bivariate time series that has been
extensively studied in the literature and a six-dimensional times series of
windspeed measurements. Section Se5.1 first introduces a method to choose the
order of the parametric VAR($p$) working model of the VNPC approach. Both
studies use the same prior settings as the simulation study (see Section
Se4.1).

Elbow criterion Se5.1

Under some weak assumptions, the VAR($p$) model with sufficiently large order
$p$ can capture the dependence structure of a Gaussian linear process to any
degree of accuracy. Accordingly, using a VAR model with a sufficiently large
order can be regarded as a nonparametric procedure. This result has been
employed by the VAR-sieve-bootstrap method proposed by meyer2015. In addition,
that result also implies that the misspecified VAR model can still be used for
the spectral density matrix estimation if the order is chosen to be
sufficiently large. This explains why standard model selection techniques such
as AIC tend to return a large order when misspecification occurs. This can
lead to increased computational time, while also potentially causing the
procedure to exhibit excessive confidence in the parametric model. Hence, we
want to use a relatively small order for the parametric working model in our
corrected likelihood procedure that is capable of describing the main features
but not minute detail of the data. With this objective in mind, we follow the
order selection approach used by the univariate corrected likelihood procedure
in kirch2019, which applies a rule-of-thumb rule to check for a bend or elbow
in a plot of model order versus negative maximum log-likelihood evaluated at
the least squares estimates. To elaborate, in most cases we will see an
obvious bend in this scree-like plot as the order increases. The penalty term
of standard penalisation techniques is tuned to find the best fit but too
small to find the elbow at which only the most obvious features are explained.
In analogy to Principal Components Analysis (PCA), which uses the scree plot
to reduce the dimension, we can also use the order where we see a clear bend
in the negative maximum log-likelihood plot. With the corresponding
autoregressive order, the main features of the data are captured and
increasing the number of parameters further will contribute very little to the
corrected likelihood procedure except computational costs.

Southern oscillation index

The Southern Oscillation Index (SOI) and the associated recruitment series are
two simultaneously recorded series used to explore the El Niño cycle
(shumway2011). The SOI is a standardised index of sea level air pressure
differences between Tahiti and Darwin, Australia. The index is related to the
El Niño Southern Oscillation (ENSO), the movement of equatorial sea surface
temperature across the Pacific Ocean and its atmospheric response. ENSO is
known to be periodically occurring every 2-7 years and last 6-18 months, and
it contains three phases: neutral, El Niño (the warm phase) and La Niña (the
cool phase) (noaa2022; statsnz2020). ENSO is an extremely important
meteorological movement in the Pacific Ocean as it has an impact on the
weather through changes in air pressure, sea temperature and wind direction.
Due to the periodic behaviour of ENSO, the observations in the SOI series are
expected to have a time-dependent structure. For our analysis, we make the
same simplifying assumption as rosen2007 that all effects (e.g. the seasonal
effect) can be represented by the autocovariance structure. The recruitment
series is composed of the number of newly spawned fish, re-scaled to the
interval $[0,100]$ as shown in Figure Fig5.1. It is believed that the fish
spawn is associated with the water temperature so we can expect a cross-
dependence between the SOI and the recruitment series. Both series consist of
a period of 453 monthly values over the years 1950-1987, which are shown in
Figure Fig5.1. The data is available in the R package astsa as datasets soi
and rec (stoffer2022).

[width=1height=0.18]graph/SOI_Rec.pdf

Figure : The SOI and recruitment series

[width=.9height=0.4]graph/SOI_REC_elbow_criterion.pdf

Figure : Negative maximum log likelihood for different VAR(p) models applied
to the SOI and Recruitment series.

Here, the standardised bivariate time series is analysed. We employ the VNPC
procedure with four orders ($p=0,2,5$ and $15$) for the parametric working
model and the VAR procedure with three orders ($p=2,5$ and $15$), and we
compare the performance of those two procedures. In fact, when the order of
the parametric working model is $0$, the VNPC procedure is the same as the VNP
procedure. The order of $5$ is chosen by the elbow criterion (in Figure ,
there is a significant bend at that order). The order of $15$ is chosen by
minimising the AIC. The arrangement of the bivariate time series is the SOI to
be the first series and the recruitment to be the second. To make the
comparison for the three procedures with different orders, we demonstrate the
results of the two series separately, the SOI in Figure  and the recruitment
series in Figure . It should be noted that since both series are monthly data,
the frequency corresponding to a year period within $[0,2\pi]$ is
$\omega_{\text{yearly}}:=2\pi/12\approx 0.52$.

[width=.9height=0.35]graph/SOI_Fit.pdf

Figure : Logarithmic spectral estimates for the SOI series by the VNPC(p)
procedure (red line) and the VAR(p) procedure (black dashed line). The
periodogram is given in grey.

[width=.9height=0.35]graph/Rec_Fit.pdf

Figure : Logarithmic spectral estimates for the recruitment series by the
VNPC(p) procedure (red line) and the VAR(p) procedure (black dashed line). The
periodogram is given in grey.

From the periodograms of SOI in Figure , it can be seen that there is an
annual peak. The peak is detected by the VNP procedure, while the spectral
density estimate is overly smooth. However, when $p=2$, the VAR procedure
performs even worse as it is too smooth to pick up any spectral density peaks.
In this case, the VNPC procedure with $p=2$ does not improve much compared to
VNP, i.e., VNPC with $p=0$. When the order is increased to $5$, which is the
relatively small order chosen by the elbow criterion, the VNPC procedure has
successfully estimated the main annual peak, in contrast to the VAR procedure.
When using the large order given by AIC, i.e., $p=15$, the performances of the
VNPC and the VAR procedures are very similar. Both estimators pick up the
annual peak and its harmonics as well as a broader low-frequency that can be
interpreted as the El Niño effect which occurs irregularly every 2–7 years.
The findings for the recruitment series in Figure  are similar. The VNPC
procedure can make use of the relatively small VAR order, $p=5$, to capture
the main annual peak in that series, while the estimates obtained by the VAR
procedure with that order are too smooth. Figure  provides the $90\%$
pointwise region for VNPC with $p=5$. The top-right plot diplays the real part
of the cross-spectrum, which is called the cospectrum; and the bottom-left
plot displays the imaginary part of the cross-spectrum, which is called the
quadrature spectrum.

[width=1height=0.35]graph/SOI_Rec_VNPC_05.pdf

Figure : Posterior credible regions for the SOI and recruitment series for the
VNPC procedure with a parametric working model with order $5$. The diagonal
pair shows the logarithmic spectral estimate for the two series and the off-
diagonal pair displays the real and imaginary parts of the cross-spectrum.
Pointwise $90\%$ regions are visualised in shaded pink and uniform $90\%$
regions are in shaded blue. The posterior median is given by the black solid
line and the periodogram is shown in grey.

The squared coherency is more adequate for investigating the correlation
between two series than the cospectrum and the quadrature spectrum. Recall
that the coherency $\kappa$ between two series at frequency $\omega$ is
defined as κ(ω—):=f12(ω)(f11(ω)f22(ω))1/2, 0≤ω≤π. This statistic can be viewed
as the counterpart of the cross-correlation function in the time domain and it
holds $|\kappa(\omega|\mf)|^{2}\leq 1$ for all $0\leq\omega\leq\pi$ (see
Section 11.6 in [brockwell1991]). We calculate the squared coherency for all
5000 spectral density estimates and then consider the median and the $90\%$
pointwise region. The result in Figure  shows that there is an obvious peak at
$\omega_{\text{yearly}}$, and this agrees with the findings by the cospectrum
and the quadrature spectrum in Figure . There are also a few minor peaks at
the higher frequencies, which are the harmonics of $\omega_{\text{yearly}}$,
and the two series are correlated at those frequencies as well. In addition,
the strong coherency shown at the low frequencies close to $0$ might be
because of the 2–7 years El Niño period mentioned previously. In fact, all the
results in this Section are in line with the findings in [rosen2007].

[width=1height=0.4]graph/SOI_Rec_Squared_Coherency_VNPC5_chol.pdf

Figure : Estimated squared coherency for the SOI and recruitment series by the
VNPC procedure with a parametric working model with order $5$. The posterior
median is given by the black line and the pointwise $90\%$ credible region is
in shaded pink.

##  Wind speed profiles

.3 [width=0.95height=0.4]graph/map_six_copy.png .55
[width=0.75height=0.43]graph/wind_speed_paper.pdf

Figure : The left graph shows the geographical locations for the six airports.
The right graph shows the standardised first differences of median wind speed
measurements every two hours from 01-06-2019 12:00 am to 31-08-2019 11:59 pm.

In this case study, we apply the VNPC method to the California wind speed data
from the Iowa State University Environmental Mesonet (IEM) Automated Surface
System (ASOS) database ([iowa2022]; [mannarano1998]). This example has been
studied by [hu2023] using a stochastic gradient variational Bayes approach for
fast Bayesian inference. In their work, they consider the wind speed in knots
at six airports in California from 01-06-2019 12:00 am to 31-08-2019 11:59 pm:
EDU (Davis), SAC and SMF (Sacramento), MRY (Monterey), SNS (Salinas), and WVI
(Watsonville). That particular time interval is used to avoid extreme values
and non-stationarities caused by rainfall and storms, which occurred in other
months. In particular, they apply their method to the standardised first
difference of median wind speed every two hours for each location. In our
study, we use the same processed data, which is accessible from the
supplementary material in [hu2023]. Figure  shows the geographical locations
for those six airports and the standardised detrended data. All six series are
treated as a six-dimensional time series and analysed simultaneously. The
order for the parametric working model suggested by the elbow criterion is
$10$ (see Figure ).

[width=1height=0.4]graph/wind_speed_all_six_elbow_criterion.pdf

Figure : Negative maximum log likelihood for different VAR(p) models applied
for the standardised detrended wind speed data.

The results in Figure  show that all six spectral densities have a prominent
peak at around $2\pi/12\approx 0.524$, which indicates a strong daily period
of wind profile patterns at all locations. This is consistent with the finding
in [hu2023]. However, in contrast to spectral estimates obtained by
variational inference in [hu2023], peaks in the VPNC estimates are much more
pronounced and sharper, for the main daily peak as well as its harmonics. The
spectral estimates at $2\pi/12\approx 0.524$ for the coastal locations WVI,
SNS and MRY have slightly larger power than for the inland airports as also
noted by [hu2023]. Furthermore, we are interested in the squared coherencies
between the windspeeds at the various pairs of airports. The estimates are
given in Figure . Not surprisingly, all pairs have a strong daily periodic
coherency, which is in accordance with the daily period found in all series in
Figure . Moreover, all pairs except for the ones with SMF have a strong second
harmonic.

[width=1height=0.5]graph/wind_speed_all_six_chol_10.pdf

Figure : Logarithmic spectral estimates for the standardised detrended wind
speed data by the VNPC procedure with a parametric working model with order
$10$. The periodogram is given by the grey line and the posterior median by
the black line. The pointwise $90\%$ posterior credible region is visualised
in shaded pink and the uniform $90\%$ region is in shaded blue.

In [hu2023] correlations tend to appear locally, between locations close by,
and no correlation between locations separated by 100-200 miles, i.e., between
those pairs where one wind profile is taken from the group (EDU, SAC, SMF) and
one from (WVI, SNS, MRY). In contrast, our results in Figure  show high
squared coherencies also between locations in different groups, most
pronounced at the daily frequency. Coherencies at higher frequencies drop off
more quickly for between-group pairs than within-group pairs, e.g. WVI vs. SNS
as compared to WVI vs. EDU. Even over geographical distances of 100–200 miles
such as between the airports in the Monterey and Sacramento region, one would
expect a daily correlation between the wind profiles, however hourly or
minutely patterns of the wind profiles could vary greatly between distant
airports but be similar for neighbouring airports. In fact, wind speed is
affected by multiple factors including wind direction, topography and local
climate. [wang2020] analysed large-scale wind patterns in California using
model simulations from the variable-resolution Community Earth System Model
(VR-CESM). They described ten wind patterns by directions and strengths in the
northern California and patterns are associated with geopotential height
pattern, geography and climate. Wind pattern is not variable with distance
necessarily and wind speed correlation may not be strictly related to distance
between locations.

[width=.9height=0.7]graph/wind_speed_all_six_squared_coherency_chol_10.pdf

Figure : Estimated squared coherency for the standardised detrended wind speed
data by the VNPC procedure with a parametric working model with order $10$.
The posterior median is given by the black line and the pointwise $90\%$
credible region is in shaded pink.


Supplemental.pdf [figure]style=plain,subcapbesideposition=center

#  Search for oscillations of fundamental constants using molecular
spectroscopy

R. Oswald, A. Nevsky, V. Vogt, S. Schiller<EMAIL_ADDRESS>kezhang@uni-
mainz.de Institut für Experimentalphysik, Heinrich-Heine-Universität
Düsseldorf, 40225 Düsseldorf, Germany  N. L. Figueroa, K. Zhang, O. Tretiak,
D. Antypas  Johannes Gutenberg-Universität Mainz, 55128 Mainz, Germany
Helmholtz-Institut, GSI Helmholtzzentrum für Schwerionenforschung, 55128
Mainz, Germany  D. Budker  Johannes Gutenberg-Universität Mainz, 55128 Mainz,
Germany  Helmholtz-Institut, GSI Helmholtzzentrum für Schwerionenforschung,
55128 Mainz, Germany  Department of Physics, University of California,
Berkeley, California 94720, USA  A. Banerjee, G. Perez  Department of Particle
Physics and Astrophysics, Weizmann Institute of Science, Rehovot, Israel
7610001

###### Abstract

A possible implication of an ultralight dark matter (UDM) field interacting
wibeginth the Standard Model (SM) degrees of freedom is oscillations of
fundamental constants. Here, we establish direct experimental bounds on the
coupling of an oscillating UDM field to the up, down, and strange quarks and
to the gluons, for oscillation frequencies between 10 Hz and $10^{8}\,$Hz. We
employ spectroscopic experiments that take advantage of the dependence of
molecular transition frequencies on the nuclear masses. Our results apply to
previously unexplored frequency bands, and improve on existing bounds at
frequencies $>5\,$MHz. We identify a sector of UDM - SM coupling space where
the bounds from Equivalence Principle tests may be challenged by next-
generation experiments of the present kind.

## Introduction

There are strong theoretical reasons to assume that fundamental constants (FC)
are, in fact, dynamical and can be effectively described as expectation values
of scalar fields (see Uzan (2011) for a review). Temporal evolution of these
fields results in a time variation of the ‘constants’ that can be searched for
at the precision frontier (see, for example, review Safronova  _et al._
(2018a)). If a scalar field constitutes ultralight dark matter (UDM)
Arvanitaki  _et al._ (2015a); Graham  _et al._ (2016) with sub-eV mass, then
its amplitude oscillates at its Compton frequency, $f_{\phi}=m_{\phi}c^{2}/h$,
where $m_{\phi}$ is the scalar-particle mass, $c$ is the speed of light in
vacuum, and $h$ is Planck’s constant.

Constructing a natural theoretical model of an UDM is challenging. However,
two concrete proposals relevant to this study have been put forward. In the
first, the UDM mass is protected by an approximate scale-invariance symmetry
Arvanitaki  _et al._ (2015a). In the second, UDM is an axion-like particle,
whose mass is protected by an approximate shift symmetry according to the
Goldstone theorem Chadha-Day  _et al._ (2021) that is broken, together with
the combined charge-parity (CP) invariance Flacke  _et al._ (2017); Choi and
Im (2016), by two independent sectors Banerjee  _et al._ (2019). This model is
inspired by the relaxion paradigm Graham  _et al._ (2015). The two mod.els are
qualitatively different, yet, in both frameworks, DM couples to the standard-
model (SM) fields either due to the fact that the couplings break scale
invariance Goldberger  _et al._ (2008) or via mixing with the Higgs Flacke
_et al._ (2017), resulting in time-varying FC. An additional theoretical
approach, that also leads to time-varying FC, is based on discrete symmetries
Hook (2018); Brzeminski  _et al._ (2020).

As neither observations nor theoretical arguments can constrain the DM-
particle mass Bertone and Tait (2018), broadband searches, such as present
here, are particularly motivated. Note that the preferred region of the model
of Refs. Banerjee  _et al._ (2019, 2021) is $m_{\phi}\gtrsim
10^{-11}{\rm\,eV}\sim{\rm kHz}\,,$ the frequency range studied here.

In the UDM-particle mass range above roughly $10^{-18}\,$eV, the most
stringent constraints on time-varying UDM have been provided by equivalence-
principle (EP) tests of gravity (see Tino  _et al._ (2020) and references
therein). Here, we argue that there is a sense in which the bounds arising
from direct DM searches are independent from any single EP test. Using this
insight we shobeginw how the EP-bounds on UDM models can be challenged by
atomic and molecular experiments in the near future. Prior to discussing
direct searches for scalar UDM, we introduce the phenomenology of EP tests.

EP tests are conveniently expressed in terms of the Eötvös parameter,
$\eta^{\rm Exp}_{\rm EP}\equiv 2|\vec{a}_{\rm A}-\vec{a}_{\rm
B})/|\vec{a}_{\rm A}+\vec{a}_{\rm B}|\,$, that is sensitive to the
differential acceleration ($\vec{a}$) of two test bodies, $\rm A$ and $\rm B$
(Damour and Donoghue (2010), for examaple). The parameter can be expressed in
terms of the relevant DM couplings to the SM fields $d_{i}$ [see Eq. (1)]. One
defines the ‘dilatonic charge’ of a body, $Q^{\rm X}_{i}=\partial\ln m^{\rm
X}/\partial\ln g_{i}\,,$ $m^{\rm X}$ being the mass of the body X and $g_{i}$
a FC. Then,

$\eta^{\rm Exp}_{\rm EP}\propto\sum_{i,j}{(\Delta Q)^{\rm Exp}_{i}d_{i}\times
Q^{\rm source}_{j}d_{j}}\,,$

with the dilatonic charge difference $(\Delta Q)^{\rm Exp}_{i}\equiv Q^{\rm
A}_{i}-Q^{\rm B}_{i}\,$.

In contrast to EP-violating acceleration searches, direct scalar-UDM
experiments probe observables in either quantum or macroscopic systems arising
due to the dependence of atomic transition energies, the length of solid
objects, or the refractive indices of materials on the FCs. For a review, see,
for example, Safronova (2019), for proposals, see Arvanitaki  _et al._
(2015b); Safronova  _et al._ (2018b); Dzuba  _et al._ (2018); Banerjee  _et
al._ (2020a); Peik  _et al._ (2021); Manley  _et al._ (2020); Geraci  _et al._
(2019); Stadnik and Flambaum (2015, 2016); Grote and Stadnik (2019), and for
experiments providing bounds on FC oscillations see Van Tilburg  _et al._
(2015); Aharony  _et al._ (2021); Antypas  _et al._ (2019); Hees  _et al._
(2016); Kennedy  _et al._ (2020); Wcisło et al. (2018); Vermeulen et al.
(2021); Savalle  _et al._ (2021); Campbell  _et al._ (2021); Aiello  _et al._
(2021). The sensitivities of different kinds of experiments are intrinsically
different. Indeed, atomic experiments are sensitive to variation in the
electron mass but are SSalmost insensitive to changes in nuclear masses,
molecular experiments probe for variation of both electron and nuclear masses,
whereas EP tests probe nuclear masses and are largely insensitive to electron
mass. Thus, EP tests and oscillating FC experiments are complementary to each
other. We further quantify this statement below.

Moreover, the level of FC oscillations might be enhanced at $\sim$kHz
frequencies and higher due to the presence of UDM halos around the Earth and
the Sun Banerjee  _et al._ (2020b). In such cases, the DM density and the
coherence time are increased, leading to increased sensitivity of given
experimental setups to FC oscillations. Such enhancement would not, however,
apply to fifth-force experiments, as in these, the test masses exchange
virtual DM particles, a process independent of the background DM density. In
recent experiments searching for FC oscillations, the investigated parameter
space was extended to frequencies higher than 1 Hz, to cover the audio and
radiofrequency (RF) range Aharony  _et al._ (2021); Antypas  _et al._ (2019);
Savalle  _et al._ (2021); Vermeulen et al. (2021).

While oscillations of the fine-structure constant $\alpha$ and $m_{\rm e}$
have received substantial attention, here, we focus on ‘nuclear’ FCs: the
quantum chromodynamics (QCD) energy scale $\Lambda_{{\rm QCD}}\simeq
0.33\,$GeV, and the masses of the light quarks. These constants determine the
nuclear mass. We show that molecular spectroscopy can be used to search for
oscillations of these FCs with fractional $10^{-14}-10^{-15}$ sensitivity over
a seven-orders-wide frequency band, from 10 Hz to 100 Hz.

## Theoretical model

To illustrate the interaction of a sub-eV scalar field $\phi$ with SM fields,
we write the low-energy effective Lagrangian as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\cal L}_{\rm
eff}\\!\\!\supset\\!\frac{\phi}{M_{\rm
Pl}}\left(\\!\sum\limits_{X}d_{m_{X}}m_{X}\bar{X}X+\frac{d_{\alpha}}{4}F^{2}+\frac{d_{g_{\rm
s}}\beta(g_{\rm s})}{2g_{\rm s}}G^{2}\\!\\!\right)\\!,$ (1)

where, $X={\rm e,\,u,\,d,\,s}$ are the fermions with mass $m_{X}$,
$F^{2}=F^{\mu\nu}F_{\mu\nu}$, $G^{2}=\frac{1}{2}{\rm
Tr}(G^{\mu\nu}G_{\mu\nu})$, $F_{\mu\nu}$, $G_{\mu\nu}$ are the electromagnetic
field and gluon field strength, respectively, $d_{j}$ are dimensionless
coupling constants, and $M_{\rm Pl}=\sqrt{\hbar c/({8\pi G_{\rm
N}})}=2.4\times 10^{18}$ GeV is the Planck mass. The parameter $g_{\rm s}$ is
the strong-interaction coupling constant, $\alpha_{\rm s}\equiv g_{\rm
s}^{2}/4\pi$. The function $\beta(g_{\rm s})$ describes the evolution
("running") of the coupling constant with energy, via the renormalization-
group equation (RGE) $\beta(g_{\rm s})/(2g_{\rm s})=-(11-2n_{f}/3)\alpha_{\rm
s}/8\pi$, with $n_{f}$ being the number of dynamical quarks.

As a consequence of the UDM-SM couplings in Eq. (1), the SM constants
effectively acquire a dependence on the scalar field,

$\displaystyle m_{X}(\phi)$ $\displaystyle=$ $\displaystyle
m_{X}\left(1+d_{m_{X}}\frac{\phi}{M_{\rm Pl}}\right),$ (2)
$\displaystyle\alpha(\phi)$ $\displaystyle\simeq$
$\displaystyle\alpha\left(1-d_{\alpha}\frac{\phi}{M_{\rm Pl}}\right),$ (3)
$\displaystyle\alpha_{\rm s}(\phi)$ $\displaystyle\simeq$
$\displaystyle\alpha_{\rm s}\left(1-\frac{2d_{g_{\rm s}}\beta(g_{\rm
s})}{g_{\rm s}}\frac{\phi}{M_{\rm Pl}}\right).$ (4)

The QCD scale $\Lambda_{\rm QCD}$ depends on $g_{\rm s}$ through the RGE and
dimensional transmutation (see, for example, Peskin and Schroeder (1995)).
Thus, the variation of $\Lambda_{\rm QCD}$ can be written in terms of the
variation of $\alpha_{\rm s}$ as

$\frac{\partial\ln\Lambda_{\rm QCD}}{\partial\phi}=-\frac{g_{\rm
s}}{2\beta(g_{\rm s})}\frac{\partial\ln\alpha_{\rm
s}}{\partial\phi}=\frac{d_{g_{\rm s}}}{M_{\rm Pl}}\,.$ (5)

The mass of a nucleus $m_{\rm N}$ is the sum of the nucleon masses, strong and
electromagnetic binding energy. Neglecting the small electromagnetic binding
energy proportional to $\alpha$, the nucleon mass depends on the QCD scale
$\Lambda_{\rm QCD}$ and the light quark masses $m_{\rm u,d,s}$ Shifman  _et
al._ (1978); Damour and Donoghue (2010). The variation of $m_{\rm N}$ can be
related to variation of FCs as Junnarkar and Walker-Loud (2013)

$\frac{\delta m_{\rm N}}{m_{\rm N}}=0.909\frac{\Delta\Lambda_{\rm
QCD}}{\Lambda_{\rm QCD}}+0.084\frac{\Delta\hat{m}}{\hat{m}}+3\times
10^{-4}\frac{\Delta\delta{m}}{\delta{m}}+0.043\frac{\Delta m_{\rm s}}{m_{\rm
s}}\,,$

where $\hat{m}\equiv(m_{u}+m_{d})/2$ is the mean mass of the up and down
quarks and $\delta m\equiv m_{u}-m_{d}$ is the mass difference. Note that, for
the contribution of $m_{s}$ to the nucleon mass, we have used the lattice QCD
result Junnarkar and Walker-Loud (2013). Combining this result with Eqs. (2)
and (5), we find

$\frac{\delta m_{\rm N}}{m_{\rm N}}=0.910\,\,\hat{Q}_{\rm
N}\cdot\vec{d}\,\,\frac{\phi(t)}{M_{\rm Pl}}\,,$ (6)

with $\vec{d}\equiv\left(d_{\alpha},d_{m_{e}},d_{g_{\rm
s}},d_{\hat{m}},d_{\delta m},d_{m_{\rm s}}\right)$ and $\hat{Q}_{\rm
N}\approx\left(0,0,0.999,0.092,3\times 10^{-4},0.047\right)$ defined to be a
unitSS-length vector. Assuming that $\phi$ is a viable UDM candidate, it can
be treated as a classical oscillating field,

$\phi({\vec{x}},t)\approx m_{\phi}^{-1}\sqrt{2\rho_{\rm
DM}^{\oplus}}\sin\left[m_{\phi}\left(t+\vec{\beta}_{\oplus}\vec{x}\right)\right]\,,$
(7)

with $\rho_{\rm DM}^{\oplus}$ and $\vec{\beta}_{\oplus}$ being the UDM density
and its typical velocity on the surface of the Earth, respectively. Gravity-
based measurements yield a weak direct bound on $\rho_{\rm DM}^{\oplus}$ (see,
for instance, Hogan and Rees (1988); Banerjee  _et al._ (2020b); Anderson  _et
al._ (2020)).

Below, we consider various scenarios for the properties of the DM around the
Earth. In the standard scenario where UDM constitutes a galactic halo Marsh
(2016), with $\rho_{\rm DM}^{\oplus}\equiv\rho_{\rm DM}^{\rm G}\simeq
0.3\,{\rm GeV/cm}^{3}$ and $\beta_{\oplus}c\simeq 220$ km/s, it is reasonable
to assume that during the UDM virialization process around the galaxy,
different patches or quasiparticles obtain random phases. This results in the
UDM field-amplitude admitting stochastic fluctuations around its commonly
assumed value Foster  _et al._ (2018); Centers  _et al._ (2020); Lisanti  _et
al._ (2021); Gramolin  _et al._ (2021). In addition, we shall consider a
scenario where the UDM forms a solar halo Anderson  _et al._ (2020) with
$\rho_{\rm DM}^{\oplus}=10^{5}\rho_{\rm DM}^{\rm G}$ and
$\beta_{\oplus}c=20\,$km/s and as long as $m_{\phi}\gtrsim 10^{-13}\,$eV the
number of patches in the halo is large and we still assume that the field
amplitude is stochastic. Finally, we discuss an Earth halo phenomenological
model, with $\rho_{\rm DM}^{\oplus}$, taking its maximally allowed value, that
depends on the SSUDM-particle mass Banerjee  _et al._ (2020b), and the UDM has
a negligible SSvelocity dispersion velocity (see Supplemental Material [Supp.
Mat.]).

Experimental approach.

## Experimental approach

Atomic clocks  can be used to search for oscillations of  the proton mass and
the nuclear g-factor. However, the accessible frequency range is
$f_{\phi}\lesssim 1\,$Hz due to the operation mode of the clocks.

A recently suggested alternative approach in this context is spectroscopy of
molecules Hanneke  _et al._ (2021); Antypas  _et al._ (2021). Their transition
frequencies contain contributions stemming from changes in rotational and
vibrational energy. Here we focus on the latter. The vibrational energy
$\hbar\omega_{{\rm vib}}$ of a diatomic molecule containing two nuclei
$N_{1},\,N_{2}$ scales approximately as $E_{{\rm Ryd}}\sqrt{m_{{\rm e}}/\mu}$,
where $\mu=m_{N_{1}}m_{N_{2}}/(m_{N_{1}}+m_{N_{2}})$ is the reduced nuclear
mass. Thus, molecular transitions with a change of vibrational energy are
sensitive to the nuclear mass. Furthermore, the electron-mass dependence is
enhanced, beyond the scaling contained in $E_{\rm Ryd}\propto m_{\rm e}$.

In a detection instrument based on spectroscopy, a reference quantum system
having a resonance frequency $\nu^{(1)}$ is interrogated by the wave of
frequency $\nu^{(2)}$ emitted from an oscillator. $\nu^{(2)}$ is tuned to the
proximity of $\nu^{(1)}$. In practice, the oscillator is often stabilized to
another reference (atomic ensemble or cavity). Both frequencies may depend on
more than one FC. The fluctuation spectrum of the frequency deviation
$\Delta\nu(t)/\nu=[\nu^{(1)}(t)-\nu^{(2)}(t)]/\nu$ is measured. The dependence
of a frequency $\nu^{(i)}$ on a particular FC $g$ may be characterized by the
fractional derivative $R^{(i)}_{g}=d\ln\nu^{(i)}/d\ln g$. A hypothetical
modulation $\delta g/g$ of a constant $g$ causes a modulation of the frequency
deviation $\delta\nu/\nu=(R_{g}^{(1)}-R_{g}^{(2)})\delta g/g$. One key
parameter of a given experiment is, therefore, the differential sensitivity
$\Delta R_{g}=R_{g}^{(1)}-R_{g}^{(2)}$, determined by the choice of reference
and oscillator.

## Apparatus and Operation

In our experiments, we use an electronic transition of molecular iodine (I2)
between the ground electronic state X and the excited electronic state B
Gerstenkorn and Luc (1985). The concepts of the experiments are shown in Fig.
1. Details are presented in the Supp. Mat.  We have performed two experiments,
A and B. In the apparatus A the reference is the well-known transition
R(56)32-0 at 532 nm ($\upsilon=0$, $\upsilon^{\prime}=32$), with sensitivity
$R_{\rm N}^{(1,{\rm A})}\simeq-0.06$. In apparatus B, the transition is
R(122)2-10 at 725 nm ($\upsilon=10$, $\upsilon^{\prime}=2$) Rakowsky  _et al._
(1989), with $R_{\rm N}^{(1,{\rm B})}\simeq 0.07$.

In both experiments, the interrogating oscillator is a laser (frequency
$\nu^{(2)}$). Apparatus A uses a monolithic Nd:YAG laser whose frequency is
actively stabilized to a reference cavity.  For this configuration, and the
considered frequency range, $R_{\alpha}^{({\rm 2})}=1$, $R_{\rm e}^{({\rm
2})}=1$ (see Supp. Mat.).  In apparatus B, the laser is a Ti:Sapphire laser
and frequencies in the range 100 kHz–100 MHz are considered. As this range is
above the acoustic cutoff frequency of the laser $f_{2}^{({\rm B})}\simeq
50\,$kHz Antypas  _et al._ (2019), the frequency $\nu^{(2,{\rm B})}$ is
essentially independent of the FCs Antypas  _et al._ (2020).

Summarizing, experiments A and B provide sensitivity to $\alpha$, $m_{\rm e}$,
and $m_{\rm N}$. For experiment A, $\Delta R_{\alpha}^{({\rm A})}\simeq
2-1=1$, $\Delta R_{\rm e}^{({\rm A})}=(1-R_{\rm N}^{(1,{\rm
A})})-1\simeq+0.06$, $\Delta R_{\rm N}^{({\rm A})}=R_{\rm N}^{(1,{\rm
A})}-0\simeq-0.06$, while for experiment B, $\Delta R_{\alpha}^{({\rm
B})}\simeq 2$, $\Delta R_{\rm e}^{({\rm B})}=(1-R_{\rm N}^{(1,{\rm
B})})-0\simeq+0.93$ $\Delta R_{\rm N}^{({\rm B})}=R_{\rm N}^{(1,{\rm
B})}-0\simeq 0.07$. In both experiments, the instantaneous frequency deviation
$\Delta\nu$ is converted into a voltage signal
$V^{(k)}(t)=D^{(k)}\Delta\nu^{(k)}(t)$, with the discriminators $D^{(k)}$
being system parameters, and $k={\rm A,B}$. This allows us to obtain the
spectrum of the fractional frequency variation $\delta\nu^{(k)}/\nu^{(k)}$.
The time-varying FC ($\alpha$, $m_{\rm{e}}$ and $m_{\rm{N}}$) contribute to
the variation according to:

$\frac{\delta\nu^{(k)}}{\nu^{(k)}}=\Delta
R_{\alpha}^{({k})}\frac{\delta\alpha}{\alpha}+\Delta
R_{\rm{e}}^{({k})}\frac{\delta m_{\rm e}}{m_{\rm e}}+\Delta
R_{\rm{N}}^{({k})}\frac{\delta m_{\rm N}}{m_{\rm N}}\ .$ (8)

The Eq. (8) is used to probe oscillations of the FC. The cavity vs. molecular
transition comparison can be used to constrain a combination of several DM -
SM coupling constants.

Figure 1: The two molecular iodine experiments to search for FC oscillations.
(a) Experiment A; (b) Experiment B.

## Search for oscillating fundamental constants

In the experiments, the lasers are tuned to the respective iodine transitions
and the signals $V^{(k)}(t)$ are recorded.  In experiment A, the voltage
$V(t)$ was recorded continuously with a 16-bit data acquisition (DAQ) system
and a sampling rate of $250\,$kSa/s. Here, $D\simeq 1\,$V/MHz. We analyzed a
set of $N=2^{34}$ samples spanning $T\simeq 19\,$h. From the data
$V_{i}=V^{\rm(A)}(t_{i})$ the normalized periodogram
$P_{k}=|\tilde{V_{k}}|^{2}/N^{2}$ was calculated. This is the same as the
discretized power spectral density (PSD), multiplied by  $1/T$. Various peaks
in the periodogram were investigated and identified as being of technical
origin, in part by shifting the interrogation-laser frequency away from the
resonance and repeating the measurement. This left no obvious candidate UDM
signals in the spectrum. SSA number of frequency intervals exhibiting spectral
peaks of technical origin are listed in the Supp. Mat. We do not give limits
for these excluded intervals, SSthat have widths of 5 Hz or smaller. From the
periodogram, the upper limit of the coupling parameters $d_{g}$ was determined
using the analysis of Ref. Derevianko (2018).  The spectral amplitude of the
recorded signal is shown  in Fig. 2a. For interpretation see Supp. Mat.

(a)
Figure 2:  (a): Experiment A. The spectral amplitude ${\cal A}=$ $\sqrt{{\rm
PSD_{F}}/T}$ of the scaled discriminator signal $\Delta\nu^{(\rm
A)}(t)=V^{(\rm A)}(t)/D^{(\rm A)}$. ${\rm PSD_{F}}$ is the optimally filtered
power spectral density (PSD), with a filter chosen appropriately for signals
having the same linewidth $f/Q_{0}$, $Q_{0}\simeq 1\times 10^{6}$, as standard
galatic halo UDM. The width of the orange band corresponds to the mean of
${\cal A}\pm\sigma(f)$. (b): Experiment B. The bound (95% confidence level) on
fluctuations of the signal $\Delta\nu^{(B)}(t)=V^{(\rm B)}(t)/D^{(\rm B)}$.

In experiment B, the voltage $V(t)$ was measured with the laser frequency
tuned either on the slope of the I2 resonance, or off-resonance, alternating
between these UDM-sensitive and insensitive acquisition modes to account for
spurious signals due to sources other than UDM. A 12-bit DAQ system sampled
$V(t)$ in successive $T=0.1$ s-long intervals at a rate of 250 MSa/s and the
corresponding periodograms were computed and continuously averaged. The
periodogram difference between the on- and off-resonance acquisition modes was
also computed and averaged over a 60-hr-long run. This spectrum will contain
power in excess of statistical noise in the presence of FC oscillations, and
it is subsequently investigated for UDM detection. A number of candidate peaks
were identified with power in excess of a 95% detection threshold, which was
computed with consideration of the ‘look elsewhere’ effect Scargle (1982). All
spurious signals were checked in auxiliary experiments and eventually
attributed to technical noise. The post-inspected spectrum is used to obtain a
constraint for $\delta\nu^{(\rm B)}/\nu^{(B)}$ that is shown in Fig. 2 b. See
Supp. Mat. for details.

We analyse the experimental data in the framework of three mentioned UDM
models that differ in terms of the amplitude of the UDM field and its
coherence time $\tau_{\rm coh}=1/(\omega_{\phi}(v_{\rm vir}/c)^{2})$:

$\displaystyle\tau_{\rm coh}=\begin{cases}5.9\times 10^{5}f_{\phi}^{-1},\
Q=1.1\times 10^{6},\,\,{\rm galactic\,halo}\\\ 7.1\times 10^{7}f_{\phi}^{-1},\
Q=9.0\times 10^{7},\,\,{\rm solar\,halo}\\\ \infty,\,\,{\rm
Earth\,halo}\,.\end{cases}$ (9)

In order to derive bounds to the UDM couplings, we assume that only one of the
constants $m_{e}$, $\alpha$, or $m_{\rm N}$ in Eq. (8) oscillates and analyze
the three cases separately. See Supp. Mat. for details. In Fig. 3, we present
our constraints together with existing EP constraints (turquoise line) on the
combined quark and gluons couplings $\hat{Q}_{\rm N}\cdot\vec{d}$ as a
function of the UDM mass.

Constraints on the variation of $\alpha$ and $m_{e}$ and $\alpha$ are
presented in Fig. 4a and Fig. 4b, respectively, alongside previous
constraints. We show only the strongest existing constraints on the relevant
parameter space. In principle, astrophysical bounds on our scenario could also
apply, however, these are typically weaker than those discussed here and are
less robust (see Budnik  _et al._ (2021); DeRocco  _et al._ (2020); Bar  _et
al._ (2020) for recent discussions). For the sake of clarity, these
constraints are presented for the standard galactic UDM halo only. Our results
cover the previously unexplored bands $10-50$ Hz and $5-10$ kHz, and improve
on existing bounds in the range $5-100$ MHz.

More generally, our experiments are sensitive to the following linear
combinations of the full set of couplings, defined in Eq. (6),(8), that can be
written as $\vec{Q}\cdot\vec{d}\,=|\vec{Q}\,|\hat{Q}\cdot\vec{d}$ with
$|\vec{Q}^{\rm A}|=1\,$, $|\vec{Q}^{\rm B}|=2.21\,$, and

$\displaystyle\hat{Q}^{\rm A}\simeq$
$\displaystyle\,\,(1.0,0.06,0.05,0.005,1.5\times 10^{-5},0.0024)\ ,$
$\displaystyle\hat{Q}^{\rm B}\simeq$
$\displaystyle\,\,(0.90,0.42,-0.027,-0.0025,-8\times 10^{-6},-0.0013)\,.$

We now discuss the complementarity between direct UDM searches and the bounds
arising from EP tests.  Since UDM is allowed to couple to all fields, the
relevant parameter space is of dimension five. One can find a direction
$\hat{Q}^{\perp}_{\rm Full}(m_{\phi})$ in this space that is orthogonal to the
best four EP-test bounds for a given mass. For example, in the mass range of
$2\times 10^{-12}\lesssim m_{\phi}/{\rm eV}\lesssim 5\times 10^{-9}$, these
are the Be-Al Wagner  _et al._ (2012), Be-Ti Schlamminger  _et al._ (2008),
Cu-Pb Smith  _et al._ (1999a) and Be-Cu Su  _et al._ (1994) experiments. and
we find $\hat{Q}^{\perp}_{\rm
Full}(m_{\phi})\simeq\big{(}0.003\,,\,-0.987\,,\,0.002\,,\,-0.001\,,\,-0.162\,\big{)}\,.$
For masses above $5\times 10^{-9}\,$eV, $\hat{Q}^{\perp}_{\rm Full}$ is
perpendicular to the Be-Al Wagner  _et al._ (2012), Be-Ti Schlamminger  _et
al._ (2008), Cu-Pb Smith  _et al._ (1999a) and Cu-Pb-alloy Nelson  _et al._
(1990) SS sensitivity vectors, with corresponding $\hat{Q}^{\perp}_{\rm
Full}$. Models of light scalar UDM with coupling direction defined according
to $\hat{Q}^{\perp}_{\rm Full}(m_{\phi})\cdot\vec{d}$ would not be constrained
by these four leading EP bounds. Note that throughout the whole the mass range
$\hat{Q}^{\perp}_{\rm Full}(m_{\phi})$ has a substantial overlap with the
$d_{m_{e}}$ direction (the second entry). Thus, experiments that are
particularly sensitive to time-variation of $m_{\rm e}$, such as the ones
being discussed here, test a sector of coupling space that the first four-best
EP experiments are insensitive to. In Fig. 5, we present our bounds on
$\hat{Q}^{\perp}_{\rm Full}\cdot\vec{d}$, projected SS(for clarity) in the
$d_{m_{e}}$ direction, as dotted red and blue lines. The fifth-best EP bound
projected onto $\hat{Q}^{\perp}_{\rm Full}(m_{\phi})$ and further on
$d_{m_{\rm e}}$, is shown by a brown dotted line. Note that we could only
calculate the projection of $\hat{Q}^{\perp}_{\rm Full}$ into the remaining
5th-best EP bound to an accuracy of 1:103 SS due to the limited precision of
the published test mass composition data. We find that in this sector of
coupling space the bounds related to our direct UDM experiments are SSonly two
to three orders weaker than the bounds from the EP tests.

Figure 3: Exclusion plot of the combination of couplings to the QCD sector,
$\hat{Q}_{\rm N}\cdot\vec{d}=0.999\,d_{g_{\rm s}}+0.092\,d_{\hat{m}}+3\times
10^{-4}\,d_{\delta m}+0.047\,d_{m_{\rm s}}$.  Turquoise line:  fifth-force/EP-
violation experiments  Smith  _et al._ (1999b); Schlamminger  _et al._ (2008);
Wagner  _et al._ (2012); Bergé  _et al._ (2018); Hees  _et al._ (2018).
Figure 4: Exclusion plot of the coupling to (a) $\alpha$ and (b) $m_{\rm e}$.
Existing constraints:  shaded regions in orange Vermeulen et al. (2021),
yellow Savalle  _et al._ (2021), pink Aharony  _et al._ (2021), green Antypas
_et al._ (2020), Magenta Aiello  _et al._ (2021), and purple Branca  _et al._
(2017).  Fifth-force/EP tests: turquoise Smith  _et al._ (1999b); Schlamminger
_et al._ (2008); Wagner  _et al._ (2012); Bergé  _et al._ (2018); Hees  _et
al._ (2018).  Figure 5:  Exclusion plot for $d_{m_{e}}$; the solid lines
assume a model where only $d_{m_{e}}\neq 0$. The dotted lines depict the
bounds for a model defined by a vector of sensitivities, $\hat{Q}_{\rm
Full}^{\perp}(m_{\phi})$, that is orthogonal to the sensitivities of four
leading EP test experiments (see Supp. Mat. for details). The bounds from our
experiments are shown in red and blue, whereas bounds from other published
experiments, shown in Fig. 4, are not shown again here, for simplicity. The
bound from the fifth-best EP test experiment in a given mass range, projected
onto the $\hat{Q}^{\perp}_{\rm Full}(m_{\phi})$ direction, and further on the
$d_{m_{\rm e}}$ direction, is shown as dotted brown line.

## Conclusion

Our molecular-spectroscopy experiments have resulted in the first bounds on
the coupling of an oscillating UDM field to the gluon and quark fields, in a
broad frequency range that spans seven decades (10 Hz–100 MHz). Within this
range, improvements on previous limits for the coupling to the electromagnetic
field and the electron field were also obtained within a small frequency
window. Our experiments are not entirely free of technical noises. A new
generation of similar experiments, with minimization of all noise sources,
long acquisition times, high sampling rates, and, possibly, multiple setups
enabling reduction of the noise level by statistical averaging, could further
improve the present limits by several orders of magnitude. Furthermore, we
have argued that there is a special class of dark matter couplings where the
bounds from equivalence principle tests are significantly SS less stringent
than expected. Consequently, in the near future experiments of the kind
described here SSmay be able to probe this class of UDM models SSwith
sensitivity competitive to EP tests.

###### Acknowledgements.

_Acknowledgments._ We thank D. Kanta for help with the project, and D.
Iwaschko and R. Gusek for electronics development. A.B. thanks the Galileo
Galilei Institute for Theoretical Physics for the hospitality and the INFN for
partial support during the completion of this work. This work was supported by
the Cluster of Excellence “Precision Physics, Fundamental Interactions, and
Structure of Matter” (PRISMA+ EXC 2118/1) funded by the German Research
Foundation (DFG) within the German Excellence Strategy (Project ID 39083149),
by the European Research Council (ERC) under the European Union Horizon 2020
research and innovation program (project Dark-OST, grant agreement No 695405,
project YbFUN, grant agreement No 947696), by the DFG Reinhart Koselleck
project. The work of A.B. was supported by the Azrieli Foundation. The work of
G.P. is supported by grants from BSF-NSF (No. 2019760), Friedrich Wilhelm
Bessel research award, GIF, the ISF (grant No. 718/18), Minerva, SABRA - Yeda-
Sela - WRC Program, the Estate of Emile Mimran, and The Maurice and Vivienne
Wohl Endowment.

## References

  * Uzan (2011) J.-P. Uzan, Living Rev. Relativ. 14, 2 (2011).
  * Safronova  _et al._ (2018a) M. S. Safronova, D. Budker, D. DeMille, D. F. J. Kimball, A. Derevianko, and C. W. Clark, Rev. Mod. Phys. 90, 025008 (2018a).
  * Arvanitaki  _et al._ (2015a) A. Arvanitaki, J. Huang, and K. Van Tilburg, Phys. Rev. D91, 015015 (2015a).
  * Graham  _et al._ (2016) P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, and W. A. Terrano, Phys. Rev. D93, 075029 (2016), 1512.06165 .
  * Chadha-Day  _et al._ (2021) F. Chadha-Day, J. Ellis, and D. J. E. Marsh, (2021).
  * Flacke  _et al._ (2017) T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta, and G. Perez, J. High Energy Phys. 2017, 50 (2017).
  * Choi and Im (2016) K. Choi and S. H. Im, JHEP 12, 093 (2016), 1610.00680 .
  * Banerjee  _et al._ (2019) A. Banerjee, H. Kim, and G. Perez, Phys. Rev. D100, 115026 (2019), arXiv:1810.01889 [hep-ph] .
  * Graham  _et al._ (2015) P. W. Graham, D. E. Kaplan, and S. Rajendran, Phys. Rev. Lett. 115, 221801 (2015), 1504.07551 .
  * Goldberger  _et al._ (2008) W. D. Goldberger, B. Grinstein, and W. Skiba, Phys. Rev. Lett. 100, 111802 (2008).
  * Hook (2018) A. Hook, Phys. Rev. Lett. 120, 261802 (2018).
  * Brzeminski  _et al._ (2020) D. Brzeminski, Z. Chacko, A. Dev, and A. Hook, (2020), arXiv:2012.02787 [hep-ph] .
  * Bertone and Tait (2018) G. Bertone and M. P. Tait, Tim, Nature 562, 51 (2018).
  * Banerjee  _et al._ (2021) A. Banerjee, E. Madge, G. Perez, W. Ratzinger, and P. Schwaller, (2021), arXiv:2105.12135 [hep-ph] .
  * Tino  _et al._ (2020) G. M. Tino, L. Cacciapuoti, S. Capozziello, G. Lambiase, and F. Sorrentino, Prog. Part. Nucl. Phys. 112, 103772 (2020).
  * Damour and Donoghue (2010) T. Damour and J. F. Donoghue, Phys. Rev. D 82, 084033 (2010).
  * Safronova (2019) M. S. Safronova, Ann. Phys. 531, 1800364 (2019).
  * Arvanitaki  _et al._ (2015b) A. Arvanitaki, J. Huang, and K. Van Tilburg, Phys. Rev. D 91, 015015 (2015b).
  * Safronova  _et al._ (2018b) M. S. Safronova, S. G. Porsev, C. Sanner, and J. Ye, Phys. Rev. Lett. 120, 173001 (2018b).
  * Dzuba  _et al._ (2018) V. A. Dzuba, V. V. Flambaum, and S. Schiller, Phys. Rev. A 98, 022501 (2018).
  * Banerjee  _et al._ (2020a) A. Banerjee, H. Kim, O. Matsedonskyi, G. Perez, and M. S. Safronova, J. High Energy Phys. 2020, 153 (2020a).
  * Peik  _et al._ (2021) E. Peik, T. Schumm, M. S. Safronova, A. Pálffy, J. Weitenberg, and P. G. Thirolf, Quantum Science and Technology 6, 034002 (2021).
  * Manley  _et al._ (2020) J. Manley, D. J. Wilson, R. Stump, D. Grin, and S. Singh, Phys. Rev. Lett. 124, 151301 (2020).
  * Geraci  _et al._ (2019) A. A. Geraci, C. Bradley, D. Gao, J. Weinstein, and A. Derevianko, Phys. Rev. Lett. 123, 031304 (2019).
  * Stadnik and Flambaum (2015) Y. V. Stadnik and V. V. Flambaum, Phys. Rev. Lett. 114, 161301 (2015).
  * Stadnik and Flambaum (2016) Y. V. Stadnik and V. V. Flambaum, Phys. Rev. A 93, 063630 (2016).
  * Grote and Stadnik (2019) H. Grote and Y. V. Stadnik, Phys. Rev. Research 1, 033187 (2019).
  * Van Tilburg  _et al._ (2015) K. Van Tilburg, N. Leefer, L. Bougas, and D. Budker, Phys. Rev. Lett. 115, 011802 (2015).
  * Aharony  _et al._ (2021) S. Aharony, N. Akerman, R. Ozeri, G. Perez, I. Savoray, and R. Shaniv, Phys. Rev. D 103, 075017 (2021).
  * Antypas  _et al._ (2019) D. Antypas, O. Tretiak, A. Garcon, R. Ozeri, G. Perez, and D. Budker, Phys. Rev. Lett. 123, 141102 (2019).
  * Hees  _et al._ (2016) A. Hees, J. Guéna, M. Abgrall, S. Bize, and P. Wolf, Phys. Rev. Lett. 117, 061301 (2016).
  * Kennedy  _et al._ (2020) C. J. Kennedy, E. Oelker, J. M. Robinson, T. Bothwell, D. Kedar, W. R. Milner, G. E. Marti, A. Derevianko, and J. Ye, Phys. Rev. Lett. 125, 201302 (2020).
  * Wcisło et al. (2018) P. Wcisło et al., Sci. Adv. 4 (2018), 10.1126/sciadv.aau4869.
  * Vermeulen et al. (2021) S. M. Vermeulen et al.,  (2021), 10.21203/rs.3.rs-322391/v1.
  * Savalle  _et al._ (2021) E. Savalle, A. Hees, F. Frank, E. Cantin, P.-E. Pottie, B. M. Roberts, L. Cros, B. T. Mcallister, and P. Wolf, Phys. Rev. Lett. 126, 051301 (2021).
  * Campbell  _et al._ (2021) W. M. Campbell, B. T. McAllister, M. Goryachev, E. N. Ivanov, and M. E. Tobar, Phys. Rev. Lett. 126, 071301 (2021).
  * Aiello  _et al._ (2021) L. Aiello, J. W. Richardson, S. M. Vermeulen, H. Grote, C. Hogan, O. Kwon, and C. Stoughton, (2021), arXiv:2108.04746 [gr-qc] .
  * Banerjee  _et al._ (2020b) A. Banerjee, D. Budker, J. Eby, H. Kim, and G. Perez, Commun. Phys. 3, 1 (2020b).
  * Peskin and Schroeder (1995) M. E. Peskin and D. V. Schroeder,  _An Introduction to quantum field theory_ (Addison-Wesley, Reading, USA, 1995).
  * Shifman  _et al._ (1978) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Phys. Lett. B 78, 443 (1978).
  * Junnarkar and Walker-Loud (2013) P. M. Junnarkar and A. Walker-Loud, Phys. Rev. D 87, 114510 (2013).
  * Hogan and Rees (1988) C. J. Hogan and M. J. Rees, Phys. Lett. B 205, 228 (1988).
  * Anderson  _et al._ (2020) N. B. Anderson, A. Partenheimer, and T. D. Wiser, (2020), arXiv:2007.11016 [hep-ph] .
  * Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Foster  _et al._ (2018) J. W. Foster, N. L. Rodd, and B. R. Safdi, Physical Review D 97 (2018), 10.1103/physrevd.97.123006.
  * Centers  _et al._ (2020) G. P. Centers, J. W. Blanchard, J. Conrad, N. L. Figueroa, A. Garcon, A. V. Gramolin, D. F. J. Kimball, M. Lawson, B. Pelssers, J. A. Smiga, A. O. Sushkov, A. Wickenbrock, D. Budker, and A. Derevianko, (2020), arXiv:1905.13650 [astro-ph.CO] .
  * Lisanti  _et al._ (2021) M. Lisanti, M. Moschella, and W. Terrano, Phys. Rev. D 104, 055037 (2021), arXiv:2107.10260 [astro-ph.CO] .
  * Gramolin  _et al._ (2021) A. V. Gramolin, A. Wickenbrock, D. Aybas, H. Bekker, D. Budker, G. P. Centers, N. L. Figueroa, D. F. J. Kimball, and A. O. Sushkov, (2021), arXiv:2107.11948 [hep-ph] .
  * Hanneke  _et al._ (2021) D. Hanneke, B. Kuzhan, and A. Lunstad, Quantum Science and Technology 6, 014005 (2021).
  * Antypas  _et al._ (2021) D. Antypas, O. Tretiak, K. Zhang, A. Garcon, G. Perez, M. G. Kozlov, S. Schiller, and D. Budker, Quantum Science and Technology 6, 034001 (2021).
  * Gerstenkorn and Luc (1985) S. Gerstenkorn and P. Luc, J. Phys. France 46, 867 (1985).
  * Rakowsky  _et al._ (1989) S. Rakowsky, D. Zimmermann, and W. E. Ernst, Appl. Phys. B 48, 463 (1989).
  * Antypas  _et al._ (2020) D. Antypas, D. Budker, V. V. Flambaum, M. G. Kozlov, G. Perez, and J. Ye, Annalen der Physik 532, 1900566 (2020).
  * Derevianko (2018) A. Derevianko, Phys. Rev. A 97, 042506 (2018).
  * Scargle (1982) J. D. Scargle, The Astrophys. J. 263, 835 (1982).
  * Budnik  _et al._ (2021) R. Budnik, H. Kim, O. Matsedonskyi, G. Perez, and Y. Soreq, Phys. Rev. D 104, 015012 (2021).
  * DeRocco  _et al._ (2020) W. DeRocco, P. W. Graham, and S. Rajendran, Phys. Rev. D 102, 075015 (2020).
  * Bar  _et al._ (2020) N. Bar, K. Blum, and G. D’Amico, Phys. Rev. D 101, 123025 (2020).
  * Wagner  _et al._ (2012) T. A. Wagner, S. Schlamminger, J. H. Gundlach, and E. G. Adelberger, Class. Quant. Grav. 29, 184002 (2012).
  * Schlamminger  _et al._ (2008) S. Schlamminger, K. Y. Choi, T. A. Wagner, J. H. Gundlach, and E. G. Adelberger, Phys. Rev. Lett. 100, 041101 (2008).
  * Smith  _et al._ (1999a) G. L. Smith, C. D. Hoyle, J. H. Gundlach, E. G. Adelberger, B. R. Heckel, and H. E. Swanson, Phys. Rev. D 61, 022001 (1999a).
  * Su  _et al._ (1994) Y. Su, B. R. Heckel, E. G. Adelberger, J. H. Gundlach, M. Harris, G. L. Smith, and H. E. Swanson, Phys. Rev. D 50, 3614 (1994).
  * Nelson  _et al._ (1990) P. G. Nelson, D. M. Graham, and R. D. Newman, Phys. Rev. D 42, 963 (1990).
  * Smith  _et al._ (1999b) G. L. Smith, C. D. Hoyle, J. H. Gundlach, E. G. Adelberger, B. R. Heckel, and H. E. Swanson, Phys. Rev. D 61, 022001 (1999b).
  * Bergé  _et al._ (2018) J. Bergé, P. Brax, G. Métris, M. Pernot-Borràs, P. Touboul, and J.-P. Uzan, Phys. Rev. Lett. 120, 141101 (2018).
  * Hees  _et al._ (2018) A. Hees, O. Minazzoli, E. Savalle, Y. V. Stadnik, and P. Wolf, Phys. Rev. D 98, 064051 (2018).
  * Branca  _et al._ (2017) A. Branca, M. Bonaldi, M. Cerdonio, L. Conti, P. Falferi, F. Marin, R. Mezzena, A. Ortolan, G. A. Prodi, L. Taffarello, G. Vedovato, A. Vinante, S. Vitale, and J.-P. Zendri, Phys. Rev. Lett. 118, 021302 (2017).

in 1,…,0 See pages \x of Supplemental.pdf