# Fluctuations of atomic energy levels due to axion dark matter

V. V. Flambaum<EMAIL_ADDRESS>I. B. Samsonov
<EMAIL_ADDRESS>School of Physics, University of New South Wales,
Sydney 2052, Australia

###### Abstract

The amplitude of the pseudoscalar (axion) or scalar field fluctuates on a time
scale of order of million field oscillation periods which is a typical
coherence time in the virialized axion galactic dark matter halo model. This
causes fluctuations of frequencies of atomic clocks on the same time scale. We
show that this effect may be employed to search for the axion and scalar field
dark matter with atomic and nuclear clocks. We re-purpose the results of the
atomic clocks experiments comparing the variations of frequencies of hyperfine
transitions in Rb and Cs atoms as well as in hydrogen atom vs cavity frequency
fluctuations, and extract new limits on the axion coupling constant $f_{a}$
for masses in the range $2.4\times 10^{-17}\text{ eV}\lesssim m\lesssim
10^{-13}\text{ eV}$. We also show that similar energy shifts arise in the
second-order perturbation theory with linear in the pseudoscalar field
interaction. These shifts may be potentially measured with nuclear clocks
based on the low-energy transition in 229Th nucleus. We propose a procedure
which could, in principle, help determine the axion mass if the axion dark
matter signal is present in experimental data sets.

##  I Introduction

Scalar and pseudoscalar particles represent promising dark matter candidates
which can fully saturate the local dark matter density
$\rho_{\text{DM}}\approx 0.4\text{ GeV/}\text{cm}^{3}$ [1, 2, 3]. If the mass
of each of these particles is low, $m\ll 1\text{ eV}$, their number per De
Broglie wavelength must be large, and the ensemble of these particles may be
considered as a classical field oscillating harmonically in every particular
point of space, $\phi=\phi_{0}\cos(\omega t)$. The oscillation frequency
$\omega$ is approximately equal to the dark matter particle mass $m$,
$\omega\approx m$, since the kinetic energy is small, $E_{k}\sim 10^{-6}m$ (in
this paper, we assume that the dark matter particles are virialized within the
standard dark matter halo model, see, e.g., Refs. [4, 5, 6]). Interaction of
the standard model particles (electron, photon, quarks, gluons) with this dark
matter field produces oscillating shifts of atomic energy levels which have
been searched for in a number of experiments, see, e.g., Refs. [7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18].

A general problem is that the mass of dark matter particle is unknown;
therefore, one should do Fourier analysis of the data to separate the
oscillating signal. However, if the interaction is proportional to the scalar
field squared,

$V=-g_{f}M_{f}\phi^{2}\bar{\psi}\psi-\frac{g_{\gamma}}{4}\phi^{2}F_{\mu\nu}F^{\mu\nu}+\ldots\,,$
(1)

the energy shift has a non-oscillating contribution owing to the identity

$\phi^{2}=\frac{1}{2}\phi_{0}^{2}[1+\cos(2\omega t)]\,.$ (2)

In Eq. (1), $\psi$ is a Dirac fermion field with mass $M_{f}$, $F_{\mu\nu}$ is
the Maxwell field strength, $g_{f}$ and $g_{\gamma}$ are the corresponding
coupling constants. In general, the scalar field may couple to other Standard
Model fields, denoted by ellipsis in Eq. (1), which are not important in the
present consideration.

Effects of quadratic-in-$\phi$ interaction in Eq. (1) may be described as an
apparent variation of the fine structure constant,
$\alpha^{\prime}=\alpha(1+g_{\gamma}\phi^{2})$, and masses of elementary
particles, $M^{\prime}_{f}=M_{f}(1+g_{f}\phi^{2})$, see, e.g., Refs. [7, 8].
For example, the mass shift immediately follows from comparison of interaction
with the scalar filed $-g_{f}M_{f}\phi^{2}\bar{\psi}\psi$ and fermion mass
term in the Lagrangian $-M_{f}\bar{\psi}\psi$. Dependence of atomic transition
frequencies on $\alpha$, quark masses and $\phi^{2}$ was studied in works [19,
20, 21, 22, 8, 9, 23, 24]. Atomic spectroscopy methods have already allowed
one to improve earlier cosmological limits on the interaction strength of low
mass scalar field $\phi^{2}$ with photons, electrons and quarks by 15 orders
in magnitude [8, 9]. These limits have recently been revisited in Ref. [25]
due to Big Bang Nucleosynthesis considerations. The experimental results were
obtained by the measurements of oscillating frequency ratios of electron
transitions in Dy/Cs [10], Rb/Cs [11], Yb/Cs [12], Sr/H/Si cavity [13],
Cs/cavity [14], Yb/Yb/Sr [15, 16], Rb/Quartz [17] where effects of the
variation of frequencies may be interpreted as variation of $\alpha$ and
fermion masses. In Ref. [26] it was proposed to search for the scalar field
dark matter with interaction (1) by measuring fluctuations of the scalar field
amplitude using magnetometer and optical atomic clock networks. In the case of
linear-in-$\phi$ interaction, the use of a network of precision-measurement
tools for searches of wave-like dark matter was proposed in Ref. [27].

The rest of this paper is organized as follows. In Sec. II we consider
quadratic-in-$\phi$ atomic energy level corrections arising in the first order
of perturbation theory and demonstrate that the mean value of these shifts
should be (approximately) equal to the standard deviation due to the
stochastic nature of the axion field amplitude. This relation allows us to
find new lab-based limits on the axion decay constant $f_{a}$ by re-purposing
the results of the experiments [11] and [13], see Sec. III. Then, in Sec. IV
we propose a procedure which, in principle, would allow one to identify the
axion dark matter signal in experiments measuring fluctuations of energy level
shifts with atomic clocks. This procedure utilizes the fact that the atomic
energy levels fluctuations caused by axion dark matter should have a different
statistical distribution from ordinary noise in the detector. In Sec. V we
show that quadratic-in-$\phi$ contributions to the atomic energy level shifts
appear also in the second order of perturbation theory. We compare these
contributions to the corresponding first-order energy level corrections by
considering the example of low-lying isomeric state in 229Th nucleus. Section
VI is devoted to a summary and discussion of the results of this paper.

We use natural units with $\hbar=c=1$.

##  II First-order perturbation theory energy level corrections due to
quadratic axion-nucleon interaction

In this section, we focus on the quadratic in the QCD axion field $\phi$
interaction with a nucleon as in Eq. (1). In Ref. [28] it was shown that this
interaction originates from the standard QCD $\theta$-term

$\frac{g^{2}\theta}{32\pi^{2}}\tilde{G}^{l\,\mu\nu}G^{l}_{\mu\nu}\,,$ (3)

with $\theta=\phi/f_{a}$, $f_{a}$ is the axion decay constant, $g$ is the
strong interaction coupling constant, $G^{l}_{\mu\nu}$ is the gluon field
strength and $\tilde{G}^{l\,\mu\nu}$ is its Hodge dual. Indeed, this axion-
gluon interaction implies a variation of the pion mass [29],

$\frac{\delta m_{\pi}}{m_{\pi}}\approx-0.05\theta^{2}\,.$ (4)

This causes the corresponding variations of nuclear magnetic moment, nuclear
mass and radius since these quantities depend on the pion mass, see Refs. [30,
22, 31, 32, 33]. As a result, the atomic energy level shift is proportional to
the scalar field squared,

$E\propto\phi^{2}\,.$ (5)

We will study statistical properties of such energy shifts within two models
of the ultralight axion dark matter. In the first one, the axion field at the
observation point is represented by a monochromatic wave with a uniformly
distributed random phase and random amplitude with Rayleigh distributions,
while in the second one we consider the axion field as a wave packet appearing
due to a spread of velocities of dark matter particles. Although the first
model may be considered as a simplified version of the second one, it allows
for a simpler treatment and exact analytical results.

###  II.1 Monochromatic wave model

The axion field dark matter may be modeled by a monochromatic scalar wave with
angular frequency $\omega=m(1+v^{2}/2)$, where $v\sim 10^{-3}$ is the most
probable speed of dark matter particles in the standard dark matter halo model
[4, 5, 6],

$\phi=\phi_{0}\cos(\omega t+\varphi)\,.$ (6)

Here $\varphi\in[0,2\pi)$ is a uniformly distributed random phase and
$\phi_{0}$ is a random amplitude governed by the Rayleigh distribution [34]

$p(\phi_{0})=\frac{2\phi_{0}}{\phi_{\text{DM}}^{2}}\exp{\left(-\frac{\phi_{0}^{2}}{\phi_{\text{DM}}^{2}}\right)}\,,$
(7)

where $\phi_{\text{DM}}=\sqrt{2\rho_{\text{DM}}}/m$ is an average amplitude of
the dark matter field. The variance of the scalar field amplitude stems from
stochastic nature of phases of dark matter particles near the Solar system but
the spread of speeds of these particles is ignored. This spread of speeds of
dark matter particles will be accounted for in the next subsection.

The model (6) ignores also possible nonvirialized dark matter streams [35] and
composite dark matter structures such as boson stars [36] or topological
defects [37], which are to be studied separately.

Substituting Eq. (6) into (5), we have the following axion field dependence of
the atomic energy shift

$E\propto\frac{1}{2}\phi_{0}^{2}[(1+\cos(2\omega t+2\varphi)]\,.$ (8)

Let us assume that the experimental integration time $t_{1}$ significantly
exceeds the oscillation period $T=\pi/\omega\approx\pi/m$, $t_{1}\gg T$. In
this case, the oscillating term in Eq. (8) averages to zero,

$\displaystyle\overline{\cos(2\omega t+2\varphi)}$
$\displaystyle=\frac{1}{t_{1}}\int_{0}^{t_{1}}\cos(2\omega t+2\varphi)dt$ (9)
$\displaystyle=\frac{\sin(2\omega t_{1}+2\varphi)}{2\omega t_{1}}\to 0$

for $\omega t_{1}\gg 1$. As a result, the quadratic interaction (1) implies
the following shift of an atomic energy level in the first order of
perturbation theory

$E\equiv\langle V\rangle=C(\phi_{0})^{2},$ (10)

where $C$ is a time-independent constant.

The problem is that only the time dependence of atomic energy level shifts
produced by new interactions can be measured accurately. Usually, the time-
independent contribution to the energy shift is hidden by uncertainties of
theoretical values of energies in multielectron atoms.

In the case of quadratic interaction (1) with the scalar field (2) this
problem may be addressed as follows. The amplitude of the scalar or
pseudoscalar (axion) field $\phi_{0}$ fluctuates on the time scale $\tau\sim
10^{6}T$ (coherence time), see, e.g., Ref. [34]. This causes fluctuations of
the energy shift (10) of atomic, molecular and nuclear transition energies.
One can set the integration time $t_{1}$ much smaller than the coherence time
$\tau$ but much larger than the oscillation period $T$,

$T\ll t_{1}\ll\tau\,.$ (11)

Repeating these measurements $N$ times such that the total measurement time
$t=Nt_{1}$ exceeds the coherence time $\tau$,

$t\gg\tau\,,$ (12)

one can measure variance of fluctuations of the energy shift $E$

$\sigma_{\phi}^{2}=\overline{(E-\overline{E})^{2}}=\overline{(E^{2})}-(\overline{E})^{2}=(\overline{E})^{2}\,.$
(13)

Here we used the probability distribution function (7) for the scalar field
amplitude. As a result, Eq. (13) allows us to identify a theoretically
calculated atomic energy level shift with experimentally measured standard
deviation of this shift,

$\sigma_{\phi}=\bar{E}\,.$ (14)

Note that the atomic energy level shift (10) due to quadratic in the axion
field interaction (1) is a random variable with exponential distribution

$p_{\phi}(E)=\left\\{\begin{array}[]{ll}\frac{1}{\bar{E}}e^{-E/\bar{E}}&\text{for
}E\geq 0\,,\\\ 0&\text{for }E<0\,,\end{array}\right.$ (15)

since the axion field amplitude follows the Rayleigh distribution. The
relation (14) is simply a property of the exponential distribution. Another
important feature of this distribution is the presence of non-vanishing higher
statistical moments such as skewness $S$ and kurtosis $K$:

$S\equiv\frac{\overline{(E-\bar{E})^{3}}}{\sigma^{3}}=2\,,\qquad
K\equiv\frac{\overline{(E-\bar{E})^{4}}}{\sigma^{4}}=9\,.$ (16)

Thus, the axion signal should manifest itself in the experimental data through
properties of the exponential distribution.

Equation (14) is derived within the assumption that the energy shift
fluctuations are caused solely by the axion field dark matter. In reality, one
has to add the effect of the noise in the detector. Let $p_{\text{noise}}(E)$
be a probability distribution function of this noise with standard deviation
$\sigma_{\text{noise}}$ and vanishing mean. Then the experimentally measured
atomic energy level shift is a random variable which follows a convolution of
these two distributions,

$p(E)=\int_{-\infty}^{\infty}p_{\phi}(E-E^{\prime})p_{\text{noise}}(E^{\prime})dE^{\prime}\,.$
(17)

As a result, the standard deviation for this combined distribution is greater
than the one for exponential distribution (14),

$\sigma>\bar{E}\,.$ (18)

The values of the higher statistical moments (16) are also different for this
distribution. They will be considered in more detail is Sec. IV where a
simulation of experimental data set will be studied.

The relation (18) will be used in the next section for extracting limits on
the axion parameter space from experiments measuring energy level shifts with
atomic clocks. This method is broadband as it does not require performing
Fourier analysis of the data. The principal assumption in this approach is
that the measurement time $t_{1}$ satisfies the conditions (11) and (12).
Given that the energy shift oscillation period $T=\pi/m$, we can convert these
conditions to the constraints on the scalar field mass if the time intervals
$t_{1}$ and $t$ are experimentally fixed,

$\max\\{\pi/t_{1},10^{6}\pi/t\\}\ll m\ll 10^{6}\pi/t_{1}\,.$ (19)

Thus, an experiment measuring fluctuations of frequencies of atomic clocks
during the time $t$ is suitable for searches of the axion and scalar dark
matter with particle mass in the range (19).

This approach may be efficient when axion or scalar mass is not too small. For
instance, assume that the averaging time is $t_{1}>10^{-6}$ s, then the
coherence time obeys $\tau>t_{1}>10^{-6}$ s, and the oscillation period is
$T\sim 10^{-6}\tau>10^{-12}$ s. Assuming also that the total measurement time
is about one day, we have $\tau<t\sim 10^{5}$ s and obtain the range of dark
matter particle masses $10^{-13}$ eV $<m<$ 0.01 eV. The QCD axion with the
mass $m\sim 10^{-5}$ eV falls within this region.

Note that a similar proposal of exploring fluctuations of the dark matter
amplitude has been recently presented in Ref. [26]. The novel feature of the
present work is the idea of using experimental value of the variance and
higher moments to find the limits on the axion decay constant. This idea will
be illustrated by the following two examples.

###  II.2 Wave packet model

Equation (6) may be considered as a toy model for the ultralight axion dark
matter because it does not take into account a spread of frequencies in this
field due to stochastic nature of velocities of dark matter particles in
Galaxy. A more realistic axion (or scalar) field dark matter model was
developed in Ref. [34] where this field in a point of observations is
represented by a wave packet

$\phi(t)=\frac{\sqrt{\rho_{\text{DM}}}}{m}\sum_{i=1}^{n}\alpha_{i}\sqrt{f(v_{i})\Delta
v}\cos\left[m\left(1+\frac{v_{i}^{2}}{2}\right)t+\varphi_{i}\right],$ (20)

where the sum is taken over the speeds $v_{i}$ of dark matter particles with
the distribution from the standard dark matter halo model [4, 5, 6]

$f(v)=\frac{v}{\sqrt{\pi}v_{0}v_{\text{obs}}}e^{-(v+v_{\text{obs}})^{2}/v_{0}^{2}}(e^{4vv_{\text{obs}}/v_{0}^{2}}-1)$
(21)

with $v_{0}\approx 220$ km/s the speed of the local rotation curve and
$v_{\text{obs}}\approx 232$ km/s the speed of the Sun in the galactic rest
frame.111The function (21) is normalized as $\int_{0}^{\infty}f(v)dv=1$. In
reality, the dark matter particle velocities in the local rotation curve vary
from the escape velocity in the Solar system $v_{\text{min}}=72$ km/s to the
galactic escape velocity $v_{\text{max}}\approx 600$ km/s. Restricting the
velocities to this interval slightly changes the overall coefficient in (21),
but we ignore this effect for simplicity.

The sum in Eq. (20) should have sufficiently large number of terms $n$ such
that the variations of speeds of particles in each interval
$[v_{i},v_{i}+\Delta v]$ may be neglected, and these particles may be
described by a monochromatic wave with uniformly distributed random phase
$\varphi_{i}$ and amplitude proportional to a random variable $\alpha_{i}$
with Rayleigh distribution

$p(\alpha_{i})=\alpha_{i}e^{-\alpha_{i}^{2}/2}\,.$ (22)

Recall that we consider experiments measuring atomic or molecular energy level
shifts which are quadratic with respect to the axion field, see Eq. (5).
Assume that such an experiment repeatedly measures atomic energy level shifts
at time instances $t_{i}$ such that $\Delta t=t_{i+1}-t_{i}=t_{1}$ obeys the
conditions (11). Therefore, fast oscillations of the axion field are averaged
out during the integration time,

$\displaystyle\langle\phi^{2}\rangle\equiv\frac{1}{\Delta
t}\int_{t_{i}}^{t_{i+1}}\phi^{2}(t)dt$
$\displaystyle=\frac{\rho_{\text{DM}}}{2m^{2}}\sum_{j,k=1}^{n}\alpha_{j}\alpha_{k}\Delta
v\sqrt{f(v_{j})f(v_{k})}$ (23) $\displaystyle\times\frac{\sin[m\Delta
t(v_{k}^{2}-v_{j}^{2})/4]\cos[m\bar{t}_{i}(v_{k}^{2}-v_{j}^{2})/2+\Delta\varphi_{kj}]}{m\Delta
t(v_{k}^{2}-v_{j}^{2})/4}\,,$

where $\Delta\varphi_{kj}=\varphi_{k}-\varphi_{j}$ and
$\bar{t}_{i}=(t_{i}+t_{i+1})/2=(i+\frac{1}{2})\Delta t$.

Assume now that the integration time $\Delta t$ is such that the following
conditions

$m\Delta t(v_{k}^{2}-v_{j}^{2})/4\ll 1\quad\forall k,j\leq n$ (24)

are satisfied. In this case Eq. (23) reduces to

$\displaystyle\langle\phi^{2}(t_{i})\rangle\approx$
$\displaystyle\frac{\rho_{\text{DM}}}{2m^{2}}\sum_{j,k=1}^{n}\alpha_{j}\alpha_{k}\Delta
v\sqrt{f(v_{j})f(v_{k})}$ (25)
$\displaystyle\times\cos[mt_{i}(v_{k}^{2}-v_{j}^{2})/2+\Delta\varphi_{kj}]\,.$

This equation specifies time fluctuations of atomic energy level shifts due to
ultralight axion dark matter,

$E_{i}\equiv E(t_{i})=C\langle\phi^{2}(t_{i})\rangle\,,$ (26)

where $C$ is a time-independent coefficient. Equations (25) and (26) show that
the time fluctuations of the atomic energy levels due to axion dark matter is
described by a wave packet with a spread of frequencies due to stochastic
nature of speeds of dark matter particles in the point of observation.

Eq. (25) may be conveniently represented as a sum of time-independent and
oscillating terms,

$\displaystyle\langle\phi^{2}(t_{i})\rangle$
$\displaystyle=\frac{\rho_{\text{DM}}}{2m^{2}}\sum_{j=1}^{n}\alpha_{j}^{2}f(v_{j})\Delta
v$ (27)
$\displaystyle+\frac{\rho_{\text{DM}}}{m^{2}}\sum_{j>k}^{n}\alpha_{j}\alpha_{k}\Delta
v\sqrt{f(v_{j})f(v_{k})}$
$\displaystyle\times\cos[mt_{i}(v_{k}^{2}-v_{j}^{2})/2+\Delta\varphi_{kj}]\,.$

Hence, after averaging over a time interval $t=N\Delta t$ significantly
exceeding periods of all oscillating terms, we have a non-vanishing mean
energy level shift

$\bar{E}=\frac{1}{N}\sum_{i=1}^{N}E_{i}=C\frac{\rho_{\text{DM}}}{2m^{2}}\sum_{j=1}^{n}\alpha_{j}^{2}f(v_{j})\Delta
v\,.$ (28)

Note that the coherence time in fluctuations of the atomic energy shifts due
to axion dark matter is usually defined as $\tau\sim 10^{6}\pi/m$. This
definition, however, needs to be formalized because these shifts are described
by a wave packet (25) with no fixed frequency. For this purpose, we consider a
time correlation function of relative energy shifts,

$\displaystyle R(\tilde{t})$ $\displaystyle\equiv\langle\Delta
E(t^{\prime}+\tilde{t})\Delta E(t^{\prime})\rangle_{t^{\prime}}$
$\displaystyle=\frac{1}{t_{w}}\int_{0}^{t_{w}}\Delta
E(t^{\prime}+\tilde{t})\Delta E(t^{\prime})dt^{\prime}\,,$ (29)

where $t_{w}$ is the averaging time window ($\tau\ll t_{w}<t-\tilde{t}$),
$\Delta E(t^{\prime})=E(t^{\prime})-\bar{E}$, $E(t^{\prime})$ is given by Eqs.
(25) and (26), and $\bar{E}$ is defined in Eq. (28). The coherence time $\tau$
may be defined now as a time interval needed for the correlation function (29)
to fall to half its original value,

$R(\tau)=\frac{1}{2}R(0)\,.$ (30)

Given this definition of the coherence time, we stress that the applicability
of this approach is given by the conditions (11) and (12).

The energy level shifts (26) represent a random variable which approximately
follows the exponential distributions because Eq. (25) is a weighted sum of a
product of Rayleigh distributed random amplitudes $\alpha_{i}$. Therefore, the
equation (14) is now satisfied only approximately,

$\sigma_{\phi}\approx\bar{E}\,,$ (31)

where $\sigma_{\phi}$ is the standard deviation of energy level shifts (26).
Higher statistical moments (16) develop their values also only approximately,
$S\approx 2$, $K\approx 9$.

When the noise in the detector is taken into considerations, the standard
deviation increases, and the relation (31) turns into the constraint (18).
Thus, this limit can be used for constraining the parameter space in the
ultralight axion dark matter model from the experiments measuring atomic
energy level shifts with atomic clocks.

##  III Limits from atomic clocks experiments

In the previous section, we have shown that the axion signal may be, in
principle, detected in the experiments comparing atomic energy level shifts
between different atomic clocks and studying their statistical properties. The
limits on the axion parameter space may be found even without full data sets
of these experiments, but just with known leading statistical moments such as
the standard deviation, skewness and kurtosis, because the relations (14) and
(16) may be considered as signatures of stochastic fluctuations of the axion
field amplitude.

In a real experiment, the noise in the detector can hide the axion signal and
spoil the properties (14) and (16). In Sec. II.1, we have shown that this
noise can increase the value of the standard deviation, and the relation (18)
should be rather used for extracting limits from experiments measuring atomic
energy level shifts with atomic clocks. We stress that the relation (18) holds
for both models of the axion dark matter considered in the previous section.

The region of applicability of the constraint (18) is given by conditions (19)
which originate from relations (11) and (12). These relations, in particular,
mean that the coherence time should be much larger than the integration time
in one measurement. In other words, there are many energy shift measurements
per axion dark matter coherence time, which, if the detector noise were
vanishing, would have nearly the same values. These energy shifts all together
would still follow the exponential distribution (at least approximately in the
case of wave packet considered in Sec. II.2) if the total experimental time
$t$ significantly exceeds the coherence time. Indeed, outside the coherence
time the values of energy shifts will differ and if they cover the whole range
of allowed values, the variance and higher moments do not change because of
the repetition of energy shifts within coherence time. This also immediately
follows from the definition of variance (and similar for other moments):

$\sigma^{2}=\frac{1}{N-1}\sum_{i=1}^{N}(E_{i}-\overline{E})^{2}.$ (32)

If we multiply each term in this sum by factor $n$, there will be compensating
increase of $N$ in denominator, $N^{\prime}=Nn$. On the other hand, if the
total duration of the experiment $t$ appears shorter than the coherence time
$\tau$, then the axion-induced energy shifts will have approximately equal
values with a small variance and the limit (18) cannot be applied any more.
Thus, the relation (18) is applicable in the region of axion mass (19).

Below we consider two experiments measuring atomic energy level shifts in
Rb/Cs [11] and H/Si cavity [13] pairs and extract new limits on the axion
coupling constant $f_{a}$ from experimentally measured standard deviation.

###  III.1 Limits from Rb/Cs experiment

Measurements of time dependence of the ratio of frequencies of Rb and Cs
hyperfine transitions were implemented in the work [11]. Using calculations in
Refs. [22, 28], we find this ratio in the form

$\frac{\delta(\nu_{\text{Rb}}/\nu_{\text{Cs}})}{\nu_{\text{Rb}}/\nu_{\text{Cs}}}=10^{-16}\frac{(1+\cos(2mt))}{m_{15}^{2}f_{10}^{2}}\frac{\phi_{0}^{2}}{\phi_{\text{DM}}^{2}}\,,$
(33)

where $m_{15}\equiv m/(10^{-15}$ eV), $f_{10}\equiv f_{a}/(10^{10}$ GeV), $m$
and $f_{a}$ are the axion mass and interaction constant. The reported standard
deviation in measurements of variation of the ratio of frequencies is
$\sigma=3\times 10^{-15}$ [11]. The averaging time in this experiment is
$t_{1}=864$ s, and the number of measurements is $N=100814$. Substituting
these values into the conditions (19), we find the limits for the axion mass
in the range $2.4\times 10^{-17}$ eV $\ll m\ll 2.4\times 10^{-12}$ eV:

$f_{a}>1.8\times 10^{9}\text{GeV}\left(\frac{10^{-15}\text{eV}}{m}\right)\,.$
(34)

The corresponding exclusion region is shown in Fig. 1 by the blue area.
Although this limit is still many orders in magnitude weaker than the QCD
axion line, Eq. (34) gives a new constraint on the axion coupling $f_{a}$ for
axion masses in the range $2.4\times 10^{-17}\text{ eV}\lesssim m\lesssim
10^{-13}\text{ eV}$ which is not covered by other lab-based experiments. Note
that in presenting these constraints in Fig. 1 we assume that the parameters
$f_{a}$ and $m$ are independent while they are related as $f_{a}m\approx
f_{\pi}m_{\pi}$ for the canonical QCD axion.

###  III.2 Limits from H/Si cavity experiment

Similar limit may be obtained from the comparison of the hydrogen hyperfine
transition with the silicon cavity eigenmode performed in Ref. [13].
Dependence of the ratio of corresponding frequencies on the fundamental
constants has been obtained in Refs. [22, 23]

$\frac{\nu_{\text{H}}}{\nu_{\text{Si}}}\propto\alpha^{3}R(Z\alpha)\frac{m_{e}}{m_{p}}g_{p}\,,$
(35)

where $m_{e}$ and $m_{p}$ are electron and proton masses, respectively,
$g_{p}$ is the proton magnetic $g$-factor, $\alpha$ is the fine structure
constant, $Z$ is the nuclear charge and $R(Z\alpha)$ is the relativistic
factor which for hydrogen and silicon is close to 1. Using calculations
presented in Ref. [28] we obtain

$\frac{\delta(\nu_{\text{H}}/\nu_{\text{Si}})}{\nu_{\text{H}}/\nu_{\text{Si}}}=10^{-15}\frac{(1+\cos(2mt))}{m_{15}^{2}f_{10}^{2}}\frac{\phi_{0}^{2}}{\phi_{\text{DM}}^{2}}\,.$
(36)

Equating this frequency variation to the dispersion of the experimental data
in Ref. [13], $\sigma\approx 3\times 10^{-15}$, we find the limit on the axion
decay constant:

$f_{a}>5.8\times 10^{9}\text{ GeV}\left(\frac{10^{-15}\text{eV}}{m}\right)\,.$
(37)

Although this constraint is comparable with that in Eq. (34), it applies to a
slightly different axion mass range $7.3\times 10^{-16}\text{ eV}<m<1.9\times
10^{-10}\text{ eV}$ which corresponds to the integration time $t_{1}=10.7$ s
222This integration time follows from Ref. [13]: The total measurement time in
the H/Si cavity experiment 2826942 s should be divided by 368 data points
presented in Fig. 1b in this work, and by a factor 720 which represents the
decimation of the original data set. Thus, $t_{1}=2826942\text{ s}/(368\times
720)\approx 10.7$ s.. The corresponding exclusion region is shown in Fig. 1 by
a pink area.

Figure 1: New limits on the axion coupling constant $f_{a}$ from re-purposing
the results of the Rb/Cs atomic clocks experiments [11] (blue shaded region)
and H/Si cavity experiment [13] (pink shaded region). These limits are
compared with the earlier found constraints from the same experiments obtained
in Ref. [28] (gray shaded regions labeled as Rb/Cs (old) and H/Si (old)
respectively). Dashed red line represents projected sensitivity of 229Th
nuclear clocks to the axion field through the second-order energy shift in the
perturbation theory, and dashed gray line corresponds to the first-order
perturbation theory correction found in Ref. [28]. The black straight line in
the right bottom corner is the QCD axion benchmark line with $f_{a}m\approx
f_{\pi}m_{\pi}$. For comparison, we included also the constraints from I2
molecular spectroscopy experiment [38], GEO 600 gravitational detector [39],
nuclear spin precession experiment (nEDM) [40], big bang nucleosynthesis (BBN)
[41], supernova explosions (SN1987A) [42], oscillation of nuclear charge
radius in Yb+ experiment [15], Rb/Quartz oscillator experiment [17], and
analysis of black hole spins [43, 44]. The red solid line labeled “Earth”
represents constraints due to possible axion emissions from the Earth [45].

This effect may also be measured in molecules where vibrational and rotational
transitions are sensitive to variation of nucleon mass, see, e.g., Ref. [38,
28]. In fact, variance in the fluctuations of energy levels was measured in
numerous papers searching for the linear drift of the fundamental constants.
This variance is linked to statistical error of the drift measurements.

##  IV Possible signature of the axion signal

In Refs. [27, 26] it was shown that the axion field may manifest itself in
correlated fluctuations of energy shifts in a network of atomic clocks or
magnetometers. Such a correlation is possible if the detectors in a network
are separated by a distance not exceeding the dark matter particle correlation
lengths. Here we will demonstrate that, in principle, it is possible to find a
signature of the ultralight axion dark matter even with a single atomic clock
measuring relative atomic energy level shift within a sufficiently long period
of time. The main idea is that fluctuations of atomic energy levels due to the
axion or scalar field dark matter are governed by a different statistical
distribution as compared with ordinary noise which is often described by a
normal Gaussian distribution, although other types of noise may also be
present.

###  IV.1 Binning of experimental data

Consider an experiment continuously measuring atomic or molecular energy level
shifts with integration time $t_{1}$ and total duration of the experiment
$t=Nt_{1}\gg t_{1}$. As a result, a series of experimental data are collected,
$E_{i}$, $i=1,2,\ldots,N$, with mean $\bar{E}$ and standard deviation
$\sigma$.

Assume that the mass $m$ of the axion (or scalar) field $\phi$ is known, and
it satisfies the conditions (19). As is argued in Sec. II, the amplitude
$\phi_{0}$ of this field fluctuates with a typical coherence time $\tau\sim
10^{6}T=10^{6}\pi/m$. This means that the fluctuations of energy shifts within
one coherence time interval $\tau$ represent the noise in the detector, while
the difference in mean energy shifts in different coherence time intervals is
caused by stochastic fluctuations of the axion (or scalar) field amplitude.

This suggests the following binning procedure of the experimental data set
with the aim to suppress the noise. The total number of measurements $N$ is
divided into $N_{\text{bin}}$ bins with $n_{1}$ data points in each bin such
that $n_{1}t_{1}=\tau$. Thus, the energy shifts may be written as $E_{ik}$
with $i=1,2,\ldots,N_{\text{bin}}$ and $k=1,2,\ldots,n_{1}$. The energy shifts
which fall within one bin may be averaged,
$\bar{E}_{i}=\frac{1}{n_{1}}\sum_{k=1}^{n_{1}}E_{ik}$. As a result, the
detector noise is averaged out, and fluctuations of $\bar{E}_{i}$ are mainly
caused by axion or scalar field dark matter. As is shown in Sec. II,
fluctuations of $\bar{E}_{i}$ should be governed (at least roughly) by the
exponential distribution. This distribution is featured by the properties (14)
and (16). By checking that the energy level shifts $\bar{E}_{i}$ approximately
satisfy these equations, one could conclude that the axion signal is present
in the data.

The problem is, however, that the mass of the axion field and, hence, the
coherence time are not known. Therefore, one has to look for the coherence
time by diving the experimental data set into bins many times with different
bin width $n_{1}$ from 1 to $N/2$. For each such binning one has to find the
set of average energy shifts $\bar{E}_{i}$ and use it for calculation of
standard deviation $\sigma$, skewness $S$ and kurtosis $K$. Thus, these
statistical moments are functions of the bin width $n_{1}$. If the axion or
scalar field signal is present in the experimental data, the values of these
statistical moments should approximately satisfy the equations (14) and (16)
for certain $n_{1}$ such that $n_{1}t_{1}=\tau$. If such number $n_{1}$ is
found, the axion mass is expressed as

$m\approx\kappa\frac{10^{6}\pi}{n_{1}t_{1}}\,.$ (38)

It is usually assumed that the coefficient $\kappa\sim 1$. Numerical
simulations presented below give $\kappa\approx 0.55$. This coefficient
depends on the specific dark matter velocity distribution $f(v)$.

This procedure may help finding an approximate value of the axion mass, which
should serve as a motivation for further experiments to search for the axion
particle with the mass near this value. It is reminiscent of the stacking
procedure proposed in Ref. [46] for optimizing the storage and statistical
analysis of experimental data.

Figure 2: Illustration of binning of measurements of atomic energy level shift
with atomic clocks to determine the axion coherence time. Blue curve
represents contribution to the energy shift from stochastic fluctuations of
the axion or scalar field amplitude. Red dots show random fluctuations of
these energy level shift measurements due to white noise in the detector.
Total measurement time $t=Nt_{1}$ is divided into time intervals equal to the
duration of the scalar field coherence time $\tau\sim 10^{6}\pi/m$ such that
each bin contains results of $n_{1}$ measurements. Energy shift measurements
in each bin are distributed quasi Gaussian with mean $\bar{E}_{i}$ and
dispersion $\sigma_{i}$. The energy shifts $\bar{E}_{i}$ in different bins
follow the exponential distribution if the bin size corresponds to the dark
matter coherence time.

###  IV.2 Numerical simulation

For an illustration of the binning procedure we perform a Monte-Carlo
simulation of experimental data with $N=10^{5}$ energy level shifts generated
as follows. First, we assume that the integration time is $t_{1}=20.7$ s, and
the total duration of the experiment is about 24 days, that is comparable with
the parameters of the experiment [13] considered is Sec. III.2. Then, we
assume that the axion mass is $m=10^{-12}$ eV, so that the conditions (19) are
satisfied. With these parameters set, we generate a function
$E_{\phi}(t)=C\langle\phi^{2}(t)\rangle$, where $\langle\phi^{2}(t)\rangle$ is
given by Eq. (27) with $n=400$ pseudorandom phases $\varphi_{j}$ and Rayleigh
distributed amplitudes $\alpha_{j}$. The coefficient $C$ is chosen, for
simplicity, such that
$\bar{E}_{\phi}\equiv\frac{1}{t}\int_{0}^{t}E_{\phi}(t)dt=1$ in some units.
Given this function, we check that it satisfies the condition (31), as well as
$S\approx 2$ and $K\approx 9$. Thus, it models the contribution to the atomic
energy level shift due to ultralight axion dark matter.

Next, we generate a time series $E_{i}=E_{\phi}(t_{i})+E_{\text{noise}}$,
where $t_{i}=i\times t_{1}$, $i=1,2,\ldots,N$, and $E_{\text{noise}}$
represents a contribution to the energy level shift due to the noise in the
detector. For simplicity, we consider a white Gaussian noise, although, more
generally, other types of noise may be present in a real experiment. We assume
that this noise is characterized by vanishing mean,
$\overline{E_{\text{noise}}}=0$, and unit standard deviation
$\sigma_{\text{noise}}=1$. This choice corresponds to commonly assumed regime
with unit signal-to-noise ratio SNR=1. A part of generated data set (red dots)
and the function $E_{\phi}(t)$ (blue curve) are shown in Fig. 2.

Figure 3: Relative values of the correlation function (39) for a mock data set
considered in Sec. IV.2. This correlation function falls by half its initial
value at $i^{*}=55$ that corresponds to the coherence time $\tau=i^{*}\times
t_{1}\approx 0.55\cdot 10^{6}\pi/m=1150$ sec.

Assume now that the axion dark matter coherence time and the axion mass are
not known, and they should be determined from the given mock data set
$\\{E_{i}\\}$. First, one has to determine the coherence time with the use of
Eq. (30). For discrete time intervals $t_{i}=i\Delta t$, $1\leq i\ll N$, this
correlation function may be represented as a sum:

$R(t_{i})=\frac{1}{N-i}\sum_{j=1}^{N-i}(E_{i+j}-\bar{E})(E_{j}-\bar{E})\,.$
(39)

For the considered data set, the relative values of this correlation function
are plotted in Fig. 3. This graph shows that the correlation function (39)
falls by half its initial value at $i^{*}=55$ corresponding to
$\tau=i^{*}\times t_{1}\approx 0.55\cdot 10^{6}\pi/m=1150$ sec.

Next, we divide the data set $\\{E_{i}\\}$ into $N_{\text{bin}}$ bins with
$n_{1}$ points in each bin, as in Fig. 2. In each bin, the energy shifts are
averaged, and values $\bar{E}_{i}$ are found, $i=1,2,\ldots,N_{\text{bin}}$.
The mean $\bar{E}$ and standard deviation $\sigma$ are calculated using these
$\bar{E}_{i}$. This calculation of $\bar{E}$ and $\sigma$ should be repeated
for different values of the bin width $n_{1}$ in the interval $1\leq n_{1}\leq
N/2$. For the mock data set under consideration, the values of $\bar{E}$ and
$\sigma$ as functions of $n_{1}$ are plotted in Fig. 4. In this figure, the
value of the standard deviation $\sigma$ decreases quasi monotonically with
$n_{1}$ and crosses the line $\sigma=\bar{E}$ near $n_{1}=i^{*}=55$. Given
this value of $n_{1}$, the axion mass is estimated with Eq. (38).

Figure 4: Simulation of the binning procedure for finding the axion coherence
time. Horizontal axis measures the number of points per bin out of total
$10^{5}$ simulated energy shift measurements. Vertical axis corresponds to
mean atomic energy level shift $\bar{E}$ in arbitrary units. Blue curve
represents the standard deviation $\sigma$ calculated for mean energy shifts
$\bar{E}_{i}$ in each bin. This curve intersects the line $\sigma=\bar{E}$
near $n_{1}=55$ (dashed vertical line). This value may be used to estimate the
axion amplitude coherence time $\tau=n_{1}t_{1}$, with $t_{1}$ the integration
time in one measurement. This value of the coherence time is in agreement with
the definition (30). The gray dashed curve shows the values of standard
deviation $\sigma$ with vanishing detector noise. Figure 5: Same a in Fig. 4,
but for the skewness $S$ and kurtosis $K$. Dashed horizontal lines show the
values of the skewness and kurtosis with vanishing detector noise. At
$n_{1}=55$ (vertical dashed line) skewness and kurtosis approximately reach
their expected values (16) specific for exponential distributions. This helps
determining the axion dark matter coherence time $\tau=n_{1}t_{1}$.

Note that for small values of $n_{1}$, the standard deviation $\sigma$ in Fig.
4 grows and reaches the value $\sigma\approx 1.4$ at $n_{1}=1$. This behavior
is natural because for small bin size the detector noise (modelled by Gaussian
distribution with vanishing mean and unit standard deviation in the present
case) dominates. The value $\sigma\approx 1.4$ is simply a standard deviation
of the convolution of the two distributions (17). Thus, our simulation
confirms that for $n_{1}=1$ (no binning), the relation (18) is satisfied, and
experimentally measured standard deviation of the full data set may be used
for extracting limits on the axion parameter space as in Sec. III.

Fig. 4 demonstrates also that when the bin size is $n_{1}\gtrsim 20$, the
detector noise is averaged out, and the standard deviation appears very close
to the gray dashed curve representing the values of the standard deviation
with no detector noise (energy shift is fully produced by the axion signal).
On the other hand, when $n_{1}>\tau/t_{1}$, the bin size becomes so large that
not only detector noise, but also the axion field amplitude fluctuations
average out, and the relation (31) is strongly violated, $\sigma<\bar{E}$.

Additional information may be obtained from the plots of the skewness $S$ and
kurtosis $K$ as functions of bin width $n_{1}$, see Fig. 5. These functions
may be found in a similar way as the standard deviation considered above. Fig.
5 demonstrates that when the bin width is small, the values of these
parameters are much smaller than the predicted ones (16) because the
contribution from the detector noise dominates in this regime. The values (16)
are approximately reached when the bin size becomes $n_{1}=\tau/t_{1}=55$,
that is indicated by vertical dashed line in Fig. 5. As a result, plots of
skewness and kurtosis confirm that the energy level shifts $\bar{E}_{i}$
roughly follow the exponential distribution when the bin width corresponds to
the axion coherence time interval. For real experiments, this could be a
strong indication of presence of the axion signal in the experimental data.

To avoid misunderstanding, we must stress that we do not make ensemble
averaging in our numerical model. The result is obtained with the fixed set of
amplitudes $\alpha_{j}$ and phases $\Delta\varphi_{kj}$ in Eq. (27). This is
an adequate approach to simulate real experimental data. A different sets of
amplitudes $\alpha_{j}$ and phases $\Delta\varphi_{kj}$ may produce slightly
different results but the conclusions will be the same.

The above simulation of the energy level shift qualitatively demonstrates the
procedure for axion mass determination. This procedure may be applied when the
noise in the detector is comparable with the expected energy level shift due
to the axion dark matter; otherwise, if the noise is too high, the axion
signal is totally washed out and the axion mass may not be found. However, if
in a real experimental data a similar behaviour is observed, this would be a
strong indication that variations of frequencies of atomic clocks are caused
by fluctuations of the amplitude of the axion dark matter. It is very tempting
to apply this procedure to the real data sets of experiments reported in Refs.
[10, 12, 11, 13, 14, 15, 17, 16].

##  V Second-order pertubation theory correction to energy levels shift due to
linear pseudoscalar interaction

Standard model spinor fields $\psi$, photon $F_{\mu\nu}$ and gluon
$G^{l}_{\mu\nu}$ fields can have the following interaction vertices with a
pseudoscalar field $\phi$:

$V=\frac{C_{f}}{f_{a}}\partial_{\mu}\phi\bar{\psi}\gamma_{5}\gamma^{\mu}\psi+C_{\gamma}\frac{\phi}{f_{a}}\tilde{F}^{\mu\nu}F_{\mu\nu}+C_{g}\frac{\phi}{f_{a}}\tilde{G}^{l\,\mu\nu}G^{l}_{\mu\nu}\,.$
(40)

Here $C_{f}$, $C_{\gamma}$ and $C_{g}$ are some dimensionless constants which
are of order $O(1)$ for the QCD axion model, but are arbitrary for the general
pseudoscalar (axion-like) particle. In particular, the last term in Eq. (40)
reduces to the QCD $\theta$-term (3) upon the substitution
$C_{g}=g^{2}/(32\pi^{2})$, or

$\theta=\frac{32\pi^{2}C_{g}\phi}{g^{2}f_{a}}\,.$ (41)

In atoms and molecules, the interaction (40) cannot produce energy levels
shifts in the first order of perturbation theory because the pseudoscalar
field mixes the states of opposite parity if one neglects a small axion
momentum corresponding to virialized dark matter particles in the standard
dark matter halo model. Thus, non-trivial corrections to the energy levels
$E_{n}$ start from the second order in the perturbation theory

$E^{(2)}=\sum_{n\neq 0}\frac{\langle 0|V|n\rangle\langle
n|V|0\rangle}{E_{0}-E_{n}}\,.$ (42)

The second-order energy corrections may be useful for studying variance of
fluctuations of the energy shift averaged over the field oscillations. Such
corrections may be significant in the cases of small energy denominators
$E_{0}-E_{n}$ in Eq. (42) which is the case of close metastable states in Dy
atom, polar molecules with rotational doublets, close levels in nuclear clock
based on 229Th. In Dy atom and molecules, the energy interval $E_{0}-E_{n}$
may be reduced to zero by application of magnetic field. Near the level
crossing the interval between the levels becomes linear in the perturbation,
$|E_{n}-E_{0}|=2|\langle 0|V|n\rangle|$. However, level widths and time-
dependent perturbation $V$ make the problem more complicated. We leave this
problem for future study.

The second-order energy correction (42) is quadratic in the pseudoscalar field
amplitude, $E^{(2)}\propto\phi_{0}^{2}$. Therefore, it is interesting to
compare this second-order energy level shift with the effect of $\theta^{2}$
in pion mass (4) discussed in Sec. II. Note that these two contributions to
the energy level shift are produced by two independent mechanisms, although
both originate from the same underlying axion-gluon interaction (3). Below we
estimate the second-order energy level shift in 229Th nucleus caused by
$CP$-violating pion-nucleon interaction $\sim\theta\pi\bar{N}N$ which was
derived in Ref. [47] where neutron EDM due to QCD $\theta$-vacuum was
calculated.

###  V.1 Energy shift of 229Th isomeric state

Consider, for example, the second-order contribution to the energy shift of
the low-lying level $E\equiv\hbar\omega=8.3$ eV of the nuclear clock
transition in 229Th [48]. This shift may be produced by $P,T$-violating
nuclear forces with non-relativistic potential of the general form [49]:

$V=\xi\vec{\sigma}\nabla V_{s}\,.$ (43)

Here $\xi$ is a coupling constant, $\vec{\sigma}$ are the Pauli matrices
corresponding to the spin of the nucleon and $V_{s}$ is the average nucleon-
nucleus potential due to the strong interaction.

In this paper, we consider a simple model where $V_{s}$ is given by an
oscillator-type potential,

$V_{s}=V_{0}(r^{2}/R^{2}-1)\,,$ (44)

with $V_{0}\simeq 50$ MeV and the nuclear radius $R$. This potential vanishes
on the boundary of the nucleus at $r=R$, and is negative inside the nucleus.
Although this potential represents a crude nuclear model, it allows us to
estimate analytically the second-order perturbative correction to the energy
level shift of an isomeric state in 229Th. More accurate and sophisticated
nuclear models would require numerical methods which are beyond the scope of
this paper.

In Ref. [49] the constant $\xi$ was expressed via a dimensionless coupling
$\eta$, $\xi=-2\times 10^{-21}\eta\,$cm, which was related with the QCD vacuum
angle $\theta$ in Ref. [50, 51, 52]: $\eta=4.4\times 10^{5}\theta$. Indeed,
the $P,T$-odd nuclear force (43) is dominated by the $\pi_{0}$ meson exchange
[51], while the coupling constants of $P,T$-odd pion-nucleon interaction were
expressed via $\theta$ in the classic paper [47]. Making use of Eq. (41), we
express $\xi$ via $C_{g}/f_{a}$:

$\xi=-8.8\times 10^{-16}\frac{32\pi^{2}C_{g}\phi}{g^{2}f_{a}}\text{cm}\,.$
(45)

Thus, the $CP$-odd potential (43) is first-order in the axion field
$\phi=\theta f_{a}$, and the corresponding second-order energy correction (42)
may be cast in the form

$E^{(2)}=\langle\delta\psi|V|\psi\rangle\,,$ (46)

where

$|\delta\psi\rangle=\sum_{n\neq 0}\frac{|n\rangle\langle
n|V|\psi\rangle}{E_{0}-E_{n}}$ (47)

is the first-order correction to the wave function $\psi$. In Ref. [49] this
correction was found in the following simple form

$\delta\psi=\xi\vec{\sigma}\nabla\psi\,.$ (48)

Substituting this function into Eq. (46), integrating by parts and using
commutation identities of Pauli matrices we find

$E^{(2)}=-3\frac{\xi^{2}V_{0}}{R^{2}}-4\frac{\xi^{2}V_{0}}{R^{2}}\langle\vec{l}\cdot\vec{s}\rangle\,,$
(49)

where $\vec{l}$ and $\vec{s}$ are the nuclear orbital momentum and spin
operators, respectively.

Remember that the lowest transition frequency in 229Th is given by the
difference between the energies of excited $3/2^{+}$ [633] and the ground
$5/2^{+}$ [631] nuclear states, $\hbar\omega=E_{3/2^{+}}-E_{5/2^{+}}$. Eq.
(49) allows us to find the frequency shift of this transition due to $P,T$-odd
hadronic interaction (43),

$\displaystyle\hbar\delta\omega$
$\displaystyle=-4\frac{V_{0}\xi^{2}}{R^{2}}\left(\langle\vec{l}\cdot\vec{s}\rangle_{3/2^{+}}-\langle\vec{l}\cdot\vec{s}\rangle_{5/2^{+}}\right)$
(50) $\displaystyle=-8\frac{V_{0}\xi^{2}}{R^{2}}\,,$

where we made use of the identity
$\langle\vec{l}\cdot\vec{s}\rangle_{3/2^{+}}-\langle\vec{l}\cdot\vec{s}\rangle_{5/2^{+}}=2$
[32].

The nuclear charge radius of 229Th is $R\approx 7.43$ fm [53]. Substituting
this value into Eq. (50) we find the relative frequency shift

$\delta\omega/\omega\approx-68\theta^{2}\,.$ (51)

This result may be compared with the first-order energy shift due to the pion
mass variation found in Ref. [28] (using the calculation of the dependence of
nuclear energy levels on pion mass from Ref. [32]): $E^{(1)}/E=2\times
10^{5}\delta m_{\pi}^{2}/m_{\pi}^{2}\approx 2\times 10^{4}\theta^{2}$. Thus,
the second-order energy shift (51) is about 300 times smaller than the first-
order contribution from pion mass $\theta$-dependence. This allows us to find
the limits on the axion constant $f_{a}$ by re-scaling the corresponding
limits from Ref. [28]. This limit is represented in Fig. 1 by the dashed line.
We stress that the variation of frequency of nuclear clock considered in this
section originates from the term $\cos(2\omega t)$ in $\phi^{2}$ rather than
from fluctuations of the (pseudo)scalar field amplitude.

##  VI Summary

In this paper, we found two new effects in the QCD axion model which
contribute to variations of fundamental constants.

The first effect appears from quadratic axion-nucleon interaction (1)
originating from the quadratic dependence of the pion mass on the axion (4).
This effect may be observed via variations of frequencies of atomic clocks due
to fluctuations of the (pseudo)scalar field amplitude. We show that these
variations of frequencies may be identified with variance of measured
fluctuations of transition frequencies in the atomic clocks, see Eq. (14).
This allows one to explore the region of QCD axion masses satisfying Eq. (19).
By re-purposing correspondingly the results of the experiments [11] and [13]
we found new laboratory limits on the axion decay constant $f_{a}$ in the
range $2.4\times 10^{-17}\text{ eV}\lesssim m\lesssim 10^{-13}\text{ eV}$.

We propose also a procedure which, in principle, could help finding signatures
of the axion signal in the energy level shift measurements in atomic clock
experiments. This procedure includes binning of the data, averaging of the
data points inside each bin and calculations of standard deviation $\sigma$,
skewness $S$, and kurtosis $K$ for the resulting distribution. The averaging
procedure allows one to suppress the detector noise while keeping the axion
signal if it is present. If the axion signal has magnitude comparable or
bigger than the noise, graphs of the correlation function $R(t)$ and $\sigma$,
$S$ and $K$ as a function of the bin width allows one to approximately find
the value of the coherence time $\tau$, axion mass $mc^{2}\sim 1.7\times
10^{6}\hbar/\tau$ and average energy shift $E$ produced by the axion field,
see Figs. 3, 4 and 5.

The other effect originates from the second-order perturbative correction to
the energy level shift due to linear-in-$\phi$ interaction. Since this effect
is expected to be small, it may manifest itself only in extremely accurate
frequency measurements with future technology based on nuclear clocks. We
estimated this shift for the low-energy nuclear transition in 229Th and found
the projected limits from the expected sensitivity of such nuclear clocks, see
Fig. 1.

Acknowledgements.— We are indebted to the Referee for pointing out that the
Rayleigh distribution for the axion field amplitude results in non-trivial
values of higher statistical moments (16) and for proposing the idea of
binning procedure for searches of possible signatures of the axion dark matter
in atomic clock experiments. We are grateful also to Dmitry Budker and Yevgeny
Stadnik for informing us about Ref. [26] and to Melina Filzinger and Nils
Huntemann for valuable comments. The work was supported by the Australian
Research Council Grants No. DP230101058 and DP200100150.

## References

  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie [1983] L. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Drukier  _et al._ [1986] A. K. Drukier, K. Freese, and D. N. Spergel, Phys. Rev. D 33, 3495 (1986).
  * Pillepich  _et al._ [2014] A. Pillepich, M. Kuhlen, J. Guedes, and P. Madau, Astrophys. J. 784, 161 (2014).
  * Evans  _et al._ [2019] N. W. Evans, C. A. J. O’Hare, and C. McCabe, Phys. Rev. D 99, 023012 (2019).
  * Arvanitaki  _et al._ [2015] A. Arvanitaki, J. Huang, and K. Van Tilburg, Phys. Rev. D 91, 015015 (2015).
  * Stadnik and Flambaum [2015] Y. V. Stadnik and V. V. Flambaum, Phys. Rev. Lett. 115, 201301 (2015).
  * Stadnik and Flambaum [2016] Y. V. Stadnik and V. V. Flambaum, Phys. Rev. A 94, 022111 (2016).
  * Van Tilburg  _et al._ [2015] K. Van Tilburg, N. Leefer, L. Bougas, and D. Budker, Phys. Rev. Lett. 115, 011802 (2015).
  * Hees  _et al._ [2016] A. Hees, J. Guéna, M. Abgrall, S. Bize, and P. Wolf, Phys. Rev. Lett. 117, 061301 (2016).
  * Kobayashi  _et al._ [2022] T. Kobayashi, A. Takamizawa, D. Akamatsu, A. Kawasaki, A. Nishiyama, K. Hosaka, Y. Hisai, M. Wada, H. Inaba, T. Tanabe, and M. Yasuda, Phys. Rev. Lett. 129, 241301 (2022).
  * Kennedy  _et al._ [2020] C. J. Kennedy _et al._ , Phys. Rev. Lett. 125, 201302 (2020).
  * Tretiak  _et al._ [2022] O. Tretiak _et al._ , Phys. Rev. Lett. 129, 031301 (2022).
  * Banerjee  _et al._ [2023] A. Banerjee, D. Budker, M. Filzinger, N. Huntemann, G. Paz, G. Perez, S. Porsev, and M. Safronova, Oscillating nuclear charge radii as sensors for ultralight dark matter (2023), arXiv:2301.10784 [hep-ph] .
  * Filzinger  _et al._ [2023] M. Filzinger, S. Dörscher, R. Lange, J. Klose, M. Steinel, E. Benkler, E. Peik, C. Lisdat, and N. Huntemann, Phys. Rev. Lett. 130, 253001 (2023).
  * Zhang  _et al._ [2023] X. Zhang, A. Banerjee, M. Leyser, G. Perez, S. Schiller, D. Budker, and D. Antypas, Phys. Rev. Lett. 130, 251002 (2023).
  * Sherrill  _et al._ [2023] N. Sherrill _et al._ , Analysis of atomic-clock data to constrain variations of fundamental constants (2023), arXiv:2302.04565 [physics.atom-ph] .
  * Dzuba  _et al._ [1999a] V. A. Dzuba, V. V. Flambaum, and J. K. Webb, Phys. Rev. Lett. 82, 888 (1999a).
  * Dzuba  _et al._ [1999b] V. A. Dzuba, V. V. Flambaum, and J. K. Webb, Phys. Rev. A 59, 230 (1999b).
  * Flambaum and Dzuba [2009] V. V. Flambaum and V. A. Dzuba, Can. J. Phys. 87, 25 (2009).
  * Flambaum and Tedesco [2006] V. V. Flambaum and A. F. Tedesco, Phys. Rev. C 73, 055501 (2006).
  * Pašteka  _et al._ [2019] L. F. Pašteka, Y. Hao, A. Borschevsky, V. V. Flambaum, and P. Schwerdtfeger, Phys. Rev. Lett. 122, 160801 (2019).
  * Flambaum and Munro-Laylim [2023] V. V. Flambaum and P. Munro-Laylim, Phys. Rev. D 107, 015004 (2023).
  * Bouley  _et al._ [2023] T. Bouley, P. Sørensen, and T.-T. Yu, JHEP 03, 104, arXiv:2211.09826 [hep-ph] .
  * Masia-Roig  _et al._ [2023] H. Masia-Roig _et al._ , Phys. Rev. D 108, 015003 (2023).
  * Derevianko [2018] A. Derevianko, Phys. Rev. A 97, 042506 (2018).
  * Kim and Perez [2022] H. Kim and G. Perez, Oscillations of atomic energy levels induced by QCD axion dark matter (2022), arXiv:2205.12988 [hep-ph] .
  * Ubaldi [2010] L. Ubaldi, Phys. Rev. D 81, 025011 (2010).
  * Flambaum  _et al._ [2004] V. V. Flambaum, D. B. Leinweber, A. W. Thomas, and R. D. Young, Phys. Rev. D 69, 115006 (2004).
  * Flambaum and Wiringa [2007] V. V. Flambaum and R. B. Wiringa, Phys. Rev. C 76, 054002 (2007).
  * Flambaum and Wiringa [2009] V. V. Flambaum and R. B. Wiringa, Phys. Rev. C 79, 034302 (2009).
  * Dinh  _et al._ [2009] T. H. Dinh, A. Dunning, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. A 79, 054102 (2009).
  * Centers  _et al._ [2021] G. P. Centers _et al._ , Nature Commun. 12, 7321 (2021).
  * Diemand  _et al._ [2008] J. Diemand, M. Kuhlen, P. Madau, M. Zemp, B. Moore, D. Potter, and J. Stadel, Nature 454, 735 (2008).
  * Eby  _et al._ [2016] J. Eby, C. Kouvaris, N. G. Nielsen, and L. C. R. Wijewardhana, JHEP 02, 028, arXiv:1511.04474 [hep-ph] .
  * Pospelov  _et al._ [2013] M. Pospelov, S. Pustelny, M. P. Ledbetter, D. F. Jackson Kimball, W. Gawlik, and D. Budker, Phys. Rev. Lett. 110, 021803 (2013).
  * Oswald  _et al._ [2022] R. Oswald _et al._ , Phys. Rev. Lett. 129, 031302 (2022).
  * Vermeulen  _et al._ [2021] S. Vermeulen _et al._ , Nature 600, 424–428 (2021).
  * Abel  _et al._ [2017] C. Abel _et al._ , Phys. Rev. X 7, 041034 (2017).
  * Blum  _et al._ [2014] K. Blum, R. T. D’Agnolo, M. Lisanti, and B. R. Safdi, Phys. Lett. B 737, 30 (2014).
  * Lucente  _et al._ [2022] G. Lucente, L. Mastrototaro, P. Carenza, L. Di Luzio, M. Giannotti, and A. Mirizzi, Phys. Rev. D 105, 123020 (2022).
  * Mehta  _et al._ [2020] V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. Mcallister, and M. J. Stott, Superradiance exclusions in the landscape of type IIB string theory (2020), arXiv:2011.08693 [hep-th] .
  * Baryakhtar  _et al._ [2021] M. Baryakhtar, M. Galanis, R. Lasenby, and O. Simon, Phys. Rev. D 103, 095019 (2021).
  * Hook and Huang [2018] A. Hook and J. Huang, JHEP 06, 036, arXiv:1708.08464 [hep-ph] .
  * Foster  _et al._ [2018] J. W. Foster, N. L. Rodd, and B. R. Safdi, Phys. Rev. D 97, 123006 (2018).
  * Crewther  _et al._ [1979] R. Crewther, P. Di Vecchia, G. Veneziano, and E. Witten, Phys. Lett. B 88, 123 (1979).
  * Seiferle  _et al._ [2019] B. Seiferle _et al._ , Nature 573, 243 (2019).
  * Flambaum  _et al._ [1984] V. V. Flambaum, I. B. Khriplovich, and O. P. Sushkov, Sov. Phys. JETP 60, 873 (1984).
  * Pospelov and Ritz [2005] M. Pospelov and A. Ritz, Annals Phys. 318, 119 (2005).
  * Flambaum  _et al._ [2014] V. V. Flambaum, D. DeMille, and M. G. Kozlov, Phys. Rev. Lett. 113, 103003 (2014).
  * de Vries  _et al._ [2015] J. de Vries, E. Mereghetti, and A. Walker-Loud, Phys. Rev. C 92, 045201 (2015).
  * Angeli and Marinova [2013] I. Angeli and K. Marinova, At. Data Nucl. Data Tables 99, 69 (2013).