# Ultralight dark matter searches
at the sub-Hz frontier with
atom multi-gradiometry

Leonardo Badurina  Valerie Gibson  Christopher McCabe  and Jeremiah Mitchell

###### Abstract

Single-photon atom gradiometry is a powerful experimental technique that can
be employed to search for the oscillation of atomic transition energies
induced by ultralight scalar dark matter (ULDM). In the sub-Hz regime the
background is expected to be dominated by gravity gradient noise (GGN), which
arises as a result of mass fluctuations around the experiment. In this work we
model the GGN as surface Rayleigh waves and construct a likelihood-based
analysis that consistently folds GGN into the sensitivity estimates of
vertical atom gradiometers in the frequency window between $10^{-3}$ Hz and 1
Hz. We show that in certain geological settings GGN can be significantly
mitigated when operating a multi-gradiometer configuration, which consists of
three or more atom interferometers in the same baseline. Multi-gradiometer
experiments, such as future versions of AION and MAGIS-100, have the potential
to probe regions of scalar ULDM parameter space in the sub-Hz regime that have
not been excluded by existing experiments.

##  1 Introduction

Atom interferometry is a rapidly-developing and powerful experimental
technique that can be employed in different configurations for a wide variety
of precision measurements (see e.g., [1, 2, 3, 4, 5, 6]). For example, two
spatially-separated atom interferometers that are referenced by the same laser
sources can be operated as a gradiometer to shed light on different topics in
fundamental physics. Specifically, owing to their exquisite sensitivity to
timings, accelerations and changes in atomic structures, atom gradiometers
could be used for: the detection of gravitational waves (GW) in the unexplored
‘mid-band’ frequency band [7, 8, 9, 10]; the search of violations of the
universality of free-fall arising from static and time-dependent fifth forces
acting on matter [11]; and the search for oscillating atomic transition
energies induced by scalar ultralight dark matter (ULDM) [12, 13], which will
be the focus of this paper.

It has been shown that several terrestrial long-baseline atom gradiometer and
atom interferometer network projects, including AION [14], MAGIS-100 [15],
MIGA [16], ELGAR [17], ZAIGA [18], and space-borne proposals, such as STE-
QUEST [19] and AEDGE [20], could probe regions of scalar ULDM parameter space
that have not been excluded by other experiments, such as atomic clocks [21],
torsion balances [22] and the MICROSCOPE experiment [23]. For example, compact
atom gradiometers such as AION-10, where the baseline is
$\sim$10\text{\,}\mathrm{m}$$, offer the possibility of detecting scalar ULDM
with a mass of approximately ${10}^{-15}\text{\,}\mathrm{e}\mathrm{V}$ [13],
which corresponds to a signal oscillating at a frequency of approximately
$1\text{\,}\mathrm{H}\mathrm{z}$. Longer baseline gradiometers like AION-km
[14] and MAGIS-km [15] could potentially detect ULDM across an even larger
frequency window: from $\sim${10}^{-3}\text{\,}\mathrm{H}\mathrm{z}$$ to
$\sim${10}^{3}\text{\,}\mathrm{H}\mathrm{z}$$, corresponding to ULDM masses
between $\sim${10}^{-17}\text{\,}\mathrm{e}\mathrm{V}$$ and
$\sim${10}^{-11}\text{\,}\mathrm{e}\mathrm{V}$$ [10].

The projected reach of these experiments is ultimately limited by fundamental
noise sources. Vertical gradiometers such as AION and MAGIS-100 are designed
to reach the atom shot-noise limit above
$\sim$1\text{\,}\mathrm{H}\mathrm{z}$$ [15], but would suffer from gravity
gradient noise (GGN) at lower frequencies [12]. This type of phase noise
arises as a result of mass density fluctuations of the ground and atmosphere
[24], which perturb the local gravitational potential around the atom clouds.
In the absence of a well-modelled GGN signal for vertical atom
interferometers, previous projections were either abruptly interrupted at
$\sim$0.3\text{\,}\mathrm{H}\mathrm{z}$$ [12, 14, 10], or arbitrarily extended
to lower frequencies by assuming atom shot noise only [13].

On the basis of the seminal GGN studies for optical gravitational wave
detectors [25, 26], Refs. [27, 28, 29] provided the first characterisation of
the impact of GGN on a vertical km-long gradiometer consisting of a pair of
atom interferometers. Within the context of gravitational wave searches, Refs.
[27, 28, 29] concluded that a pair of atom interferometers provide a sizable,
though not dramatic, reduction of the GGN noise between 0.1 Hz and 1 Hz
compared to large scale optical interferometers. Subsequently, Ref. [30]
showed how GGN from seismic activity and pressure variations dominates the
sub-Hz background in the MIGA network of atom interferometers. In light of the
fact that GGN is characterised by an exponential profile that decays with
radial distance from the Earth’s surface, while a gravitational wave signal
scales linearly with the length of the baseline, Ref. [31] showed that GGN
could in principle be distinguished from a gravitational wave signal by
placing several interferometers along the baseline. Furthermore, Ref. [31]
demonstrated that increasing the number of interferometers along the baseline
provides better suppression of the GGN signal, in agreement with Ref. [30],
which found a similar suppression for a horizontal configuration. These
previous studies [27, 28, 29, 30, 31] have focused on GGN effects in the
context of gravitational wave searches but have not considered the impact of
GGN on generic physics searches, such as ULDM-dedicated runs, using robust
statistical techniques.

In this paper we develop the formalism to consider the impact of GGN on ULDM
searches with a long-baseline vertical gradiometer. Such a formalism will be
of paramount importance to characterise the sensitivity of gradiometers, such
as long-baseline versions of AION and MAGIS-100, to scalar ULDM across the
frequency range from ${10}^{-3}\text{\,}\mathrm{H}\mathrm{z}$ to
$1\text{\,}\mathrm{H}\mathrm{z}$. In contrast to previous studies, we utilise
a likelihood-based analysis, which is more commonly used in the particle-
physics community, so that the impact of GGN on ULDM searches can be placed on
a more robust statistical footing. In our analysis, we simultaneously account
for an ULDM signal and a background dominated by atom shot noise and
underground GGN. We also extend previous studies by robustly quantifying the
extent to which GGN underground can be mitigated in a multi-gradiometer
experiment that consists of three or more coupled atom interferometers along
the same baseline.

This paper is structured as follows. In section 2 we briefly review how atom
gradiometers can be used for fundamental physics searches, and define the key
features of an atom multi-gradiometer. In section 3 we present the statistical
properties of the ULDM signal and backgrounds in an atom multi-gradiometer to
then construct in section 4 the appropriate likelihood analysis for setting
projected upper limits. In section 5 we discuss the impact of GGN on ULDM
searches, firstly in the context of a single gradiometer, before proceeding to
show how a multi-gradiometer configuration can suppress GGN in different
configurations. Finally, we provide a discussion of these results in section
6. Appendices A-D provide further details and derivations of calculations to
support the results in the main body of the paper.

##  2 Atom gradiometry for precision measurements

An atom interferometer is an experiment that compares the phase between
coherent, spatially delocalised quantum superpositions of atom clouds [32].
Experimentally, this can be achieved by constructing a two-level system,
composed of an excited state $|e\rangle$ and a ground stated $|g\rangle$ with
energy separation $\omega_{A}$.111Unless stated otherwise, we use natural
units and set $\hbar=c=1$. The electronic transitions
$|e\rangle\leftrightarrow|g\rangle$, as well as the degree of quantum
superposition and their physical separation, can be controlled by using a
system of lasers. In this paper, we consider single-photon atom
interferometers, which operate with atoms whose two-level system is controlled
by single-photon transitions, as is the case for the ${}^{87}\mathrm{Sr}$
isotope: the ground state corresponds to
$5\mathrm{s}^{2}\,{}^{1}\mathrm{S}_{0}$, the excited state to
$5\mathrm{s}5\mathrm{p}^{3}\,\mathrm{P}_{0}$, and $\omega_{A}\approx
2.697\times 10^{15}~{}\mathrm{rad}/\mathrm{s}$ (corresponding to a wavelength
$\lambda_{A}\approx$698\text{\,}\mathrm{n}\mathrm{m}$$). The superposition of
ground and excited states is controlled via $\pi/2$-pulses (i.e. beam-splitter
pulses), after which the excited state experiences a momentum kick in the
vertical direction that allows for the two atomic states to spatially separate
during propagation. The atom’s internal state and linear momentum can be
reversed using $\pi$-pulses (i.e. mirror pulses).222A $\pi/2$-pulse is defined
as a pulse of resonant radiation that has a duration of $\pi/(2\Omega)$, where
$\Omega$ is the Rabi frequency. A $\pi$-pulse has a duration of $\pi/\Omega$.

To increase the separation between atom clouds, large momentum transfer (LMT)
mirror pulses have been proposed [33]. With our convention, which follows Ref.
[13], a sequence characterised by $n$ LMT ‘kicks’ features $4n-1$ pulses in
total, of which two are $\pi/2$-pulses and the remaining $4n-3$ are
$\pi$-pulses. These are distributed such that at the beginning and the end of
the sequence, there are ($n-1$) $\pi$-pulses (the ‘beamsplitter sequence’),
while around time $T$, known as the interrogation time, there are ($2n-1$)
$\pi$-pulses (the ‘mirror sequence’). With our labelling convention, the
$\pi$/2-pulses are emitted at times 0 and $2T$. An example of a sequence with
$n=4$ LMT kicks is shown in the left panel of Fig. 1. For clarity, the time
axis in Fig. 1 is not shown to scale: the typical interrogation time $T$ is
$\mathcal{O}(\mathrm{s})$ while the time between individual pulses is
typically $\mathcal{O}(100\,\mathrm{ns})$ [34, 35].

Figure 1: Schematic representation of an atom multi-gradiometer. The left
panel shows the not-to-scale spacetime diagram of the $i^{\mathrm{th}}$ atom
interferometer’s sequence with $n=4$ LMT kicks, where the atom starts in the
ground state at position $z_{i}$. The atom’s excited ($|e\rangle$) and ground
($|g\rangle$) states are shown in blue and cyan, respectively. $\pi/2$\- and
$\pi$-pulses are displayed as wavy lines in fuchsia (dashed) and red (solid),
respectively. Atom-light interactions are indicated with black dots. The left
panel is to be interpreted as a blow-up of each of the right panel’s black
squares, which represent the atom interferometers along the vertical baseline
of length $L$ and are labeled by $\text{AI-}i$, where
$i=\\{1,2,...,\mathcal{N}\\}$. The vertical position of the $i^{\mathrm{th}}$
interferometer is labeled by $z_{i}$. The origin of our coordinate system is
at the Earth’s surface at the centre of the interferometer shaft, so the
smallest vertical position in this underground detector is $z_{1}$. Atom
gradiometers are represented by braces connecting spatially separated pairs of
atom interferometers. The atom gradiometer between the $i^{\mathrm{th}}$ and
$j^{\mathrm{th}}$ atom interferometer is labeled as $\text{AG-}(i,j)$.

In the final stage of the interferometric sequence, each atom is in a
superposition of two velocity states, and thus spatially separate. After an
appropriate drift time, the two paths spatially separate allowing for the
determination of the atom populations via atom fluorescence and point source
atom interferometry [36, 37]. The probability that an atom will be found in a
particular output port depends on the relative phase $\Phi$ acquired along the
two paths of the atom interferometer and can be measured from the number of
atoms in the excited or ground states, $N_{e}$ and $N_{g}$ respectively, as

$\frac{N_{e}-N_{g}}{N_{g}+N_{e}}=C\cos{\Phi}\,,$ (2.1)

where $C\leq 1$, which is known in the literature as constrast, characterises
the amplitude of these oscillations [38, 39]. The contrast is maximised when
the relative displacement between the wave packets vanishes at the end of the
interferometric sequence. Hence, knowing the value of $C$, the total phase
difference $\Phi$ can be determined by using the data collected at the
interferometer’s output ports and Eq. (2.1).

The phase difference $\Phi$ results from both the free-fall evolution of the
quantum state along each path, and the local phase of the laser which is
imprinted onto the atom during each atom-light interaction. Since the
interferometer sequence compares the motion of the atom to the reference frame
defined by the laser pulses, $\Phi$ is highly sensitive to inertial forces
acting on the atom during the interferometer sequence. In this sense, atom
interferometers can be used as accelerometers. Additionally, since the
interferometer sequence is composed of paths during which the atomic state
changes from the ground to the excited state, $\Phi$ is also highly sensitive
to the constants that govern atomic transitions. In this sense, atom
interferometers can be understood as atomic clocks.

Two spatially separated atom interferometers that are referenced by the same
laser pulse via single-photon atomic transitions can be operated as a
gradiometer [8]. The distinctive signal of a coherent oscillating ULDM field
[12] and GW passing through the experiment [11] can be tested with a
gradiometer by finding the difference between the phases measured by two atom
interferometers. Vertical gradiometers like MAGIS-100 and AION, as opposed to
horizontal configurations envisaged for MIGA, ELGAR and ZAIGA, consist of
interferometers that are positioned along the direction of the Earth’s
gravitational acceleration. A key advantage of the vertical gradiometer
experimental set-up over a lone single-photon atom interferometer relies on
the possibility of effectively attenuating laser noise. By operating the same
laser pulse between both interferometers, the laser noise affecting the
experiment’s sensitivity cancels in a differential measurement [8].

###  2.1 Atom multi-gradiometry

We define an experiment consisting of $\mathcal{N}\geq 2$ atom interferometers
that are referenced by the same set of lasers and located within the same
baseline of length $L$ as a an atom multi-gradiometer. This is to be
contrasted with two atom gradiometers operating in different baselines,
irrespective of their synchronisation. A schematic representation of this
setup is shown in the right panel of Fig. 1, where the black squares represent
the $\mathcal{N}$ AIs along the baseline. We will assume that the AIs can be
situated at any point along the baseline and in particular, that they do not
have to be equally spaced.

An atom gradiometer (AG) measurement consists of taking the difference between
the phases collected by two AIs. Figure 1 shows an example of three AG
measurements: using the first and second AIs; the
$(\mathcal{N}-1)^{\mathrm{th}}$ and $\mathcal{N}^{\mathrm{\,th}}$ AIs; and the
second and $(\mathcal{N}-1)^{\mathrm{th}}$ AIs. In total there are
$\mathcal{N}(\mathcal{N}-1)/2$ unique gradiometer measurements that can be
performed. It then follows that a multi-gradiometer experiment employing two
atom interferometers corresponds to a single atom gradiometer, as considered
in previous works (see e.g., [40, 12, 14, 13]).333The set-up envisaged here
differs from the interferometer array planned for the MIGA and ELGAR
experiments, where $\mathcal{N}$ independent gradiometer measurements are
performed on two horizontal baselines, which feature $2\mathcal{N}$
interferometers each [41, 30].  Since the interferometers are located within
the same baseline and are referenced by the same laser, the AIs will be
characterised by the same interrogation time and number of LMT kicks. To
simplify the subsequent analysis, we will assume that the launch velocity and
atom number are identical for all atoms in all of the AIs.

##  3 Characterising signal and noise in an atom gradiometer

Interactions between scalar ULDM and electrons or photons can lead to
oscillations in fundamental constants (see e.g. [42, 43, 44]). These ULDM-
induced oscillations in turn alter atomic transition energies, which can be
searched for with atom interferometer experiments [12, 13]. In this section,
we characterise the properties of the ULDM signal and the backgrounds
considered in this work, namely GGN from seismic effects and atom shot noise,
as recorded by an atom multi-gradiometer operating in a vertical
configuration. While an interferometer’s atom shot noise is an irreducible
noise source and is uncorrelated between detectors, the ULDM signal and the
GGN background will induce non-trivial cross-correlations between pairs of
interferometers that can be exploited to improve the sensitivity to low
frequency ULDM signals.

An atom interferometer measures a time-dependent signal by collecting phase
differences over a campaign of duration $T_{\mathrm{int}}$: we refer to this
as the integration time and it is zero at the time of the first phase
measurement. Measurements are made at a sampling rate $1/\Delta t$, where
$\Delta t$ is the temporal separation between successive measurements. We
assume that $\Delta t$ takes the same value over the measurement campaign. The
integration time satisfies the relation $T_{\mathrm{int}}=(N-1)\Delta t$,
where $N$ is the number of measurements, and we will use the approximation
$N\approx T_{\mathrm{int}}/\Delta t$ since $N\gg 1$.

We define the phase measured by the $i$th interferometer ($\text{AI-}i$) at
the end of a sequence that starts at time $m\Delta t$ as

$\Phi^{(i)}_{m}=\Phi^{(i)}_{\mathrm{Signal},m}+\Phi^{(i)}_{\mathrm{Noise},m}\,,$
(3.1)

where $i\in\\{1,2,...,\mathcal{N}\\}$, $m\in\\{0,1,...,N-1\\}$,
$\Phi^{(i)}_{\mathrm{Signal},m}$ is the signal phase and
$\Phi^{(i)}_{\mathrm{Noise},m}$ is the noise or background contribution. The
gradiometer phase between the $i$th and $j$th AI, which will be referred to as
$\text{AG-}(i,j)$, is defined as

$\Phi^{(i,j)}_{m}=\Phi^{(i)}_{m}-\Phi^{(j)}_{m}=\Phi^{(i,j)}_{\mathrm{Signal},m}+\Phi^{(i,j)}_{\mathrm{Noise},m}\;,$
(3.2)

where $i,j\in\\{1,2,...,\mathcal{N}\\}$, $j>i$, and the gradiometer signal and
noise contributions are respectively defined as

$\displaystyle\Phi^{(i,j)}_{\mathrm{Signal},m}$
$\displaystyle=\Phi^{(i)}_{\mathrm{Signal},m}-\Phi^{(j)}_{\mathrm{Signal},m}\,,$
(3.3) $\displaystyle\Phi^{(i,j)}_{\mathrm{Noise},m}$
$\displaystyle=\Phi^{(i)}_{\mathrm{Noise},m}-\Phi^{(j)}_{\mathrm{Noise},m}\;.$
(3.4)

For ULDM searches, $\Phi^{(i,j)}_{\mathrm{Signal},m}$ would correspond to the
differential ULDM phase recorded by the $i$th and $j$th AI. Instead, for GW
searches, the phase would correspond to the GW gradiometer signal. In either
case, the phase noise $\Phi^{(i,j)}_{\mathrm{Noise},m}$ captures all
background contributions, including static effects, such as the Earth’s
gravitational field experienced by the atoms in free-fall, and time-dependent
effects, such as GGN and atom shot noise. Because of the nature of the signal
of interest, in this paper we shall focus on atom shot noise and GGN, as these
are expected to be the dominant noise sources for long-baseline interferometry
at high and low frequencies, respectively.

###  3.1 Ultralight dark matter signal

ULDM can be modeled as a temporally and spatially oscillating non-relativistic
classical field. This follows from its high occupation number and its small
mean velocity and velocity dispersion, which is characteristic of dark matter
in the Milky Way [45]. We assume that the ULDM field can be decomposed into
Fourier modes with frequencies sampled from the ULDM speed distribution, and
random phases sampled from a uniform distribution between 0 and $2\pi$. We
assume a Maxwellian distribution of ULDM speeds according to the Standard Halo
Model (SHM) [46, 47] given by

$f_{\rm{DM}}(v)=\frac{v}{\sqrt{2\pi}\sigma_{v}v_{\mathrm{obs}}}e^{-\left(v+v_{\mathrm{obs}}\right)^{2}/(2\sigma_{v}^{2})}\left(e^{4vv_{\mathrm{obs}}/(2\sigma_{v}^{2})}-1\right)\,,$
(3.5)

where $\sigma_{v}$ is the velocity dispersion which is set, at the solar
position, by the value of the local standard of rest
$v_{0}=\sqrt{2}\sigma_{v}\approx$238\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$,
and
$v_{\mathrm{obs}}\approx$252\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$
is the average speed of the Earth relative to the halo rest frame [48].
Although the DM speed distribution is characterised by a cut-off at the escape
velocity,
$v_{\rm{esc}}\sim$800\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$ in
the Earth’s frame, and may feature additional, anisotropic components (e.g.,
[49, 50, 51, 52]), the simple form in Eq. (3.5) with no cut-off and no
anisotropies is sufficient for our discussion.

The atom interferometer’s integration time sets the frequency resolution,
$\Delta f=1/T_{\mathrm{int}}$, which in turn sets the ULDM speed $\Delta v$
that can be resolved, since $\Delta f=m_{\phi}v\Delta v/2\pi$, where
$m_{\phi}$ is the mass of the scalar field. The former relation follows from
the properties of the discrete Fourier transform of the data, while the latter
follows from the ULDM’s kinetic energy definition. In an experiment capable of
resolving the ULDM speed resolution to an accuracy of $\Delta v$, the scalar
ULDM field can be written as

$\phi(t)=\frac{\sqrt{\rho_{\mathrm{DM}}}}{m_{\phi}}\sum_{a}{\alpha_{a}\sqrt{F_{\mathrm{DM}}(v_{a})}\cos{\left[\omega_{a}t+\theta_{a}\right]}}\,,$
(3.6)

where $\rho_{\mathrm{DM}}$ is the local dark matter density to which we assign
the value $0.3~{}\mathrm{GeV}/\mathrm{cm}^{3}$, $\omega_{a}\simeq
m_{\phi}(1+v_{a}^{2}/2)$ is the angular frequency of the ULDM wave for a speed
$v_{a}$,

$F_{\mathrm{DM}}(v_{a})=\int_{v_{a}-\Delta v/2}^{v_{a}+\Delta
v/2}dv\,f_{\mathrm{DM}}(v)\;,$ (3.7)

and the sum is over experimentally-resolvable speeds which are indexed by the
variable $a$ [53]. The variable $\theta_{a}\in[0,2\pi)$ is a random phase,
which we assume is uniformly distributed, while $\alpha_{a}$ is Rayleigh
distributed with $\langle\alpha_{a}^{2}\rangle=2$; its probability density
function is given by

$P(\alpha_{a})=\alpha_{a}\exp\left(-\frac{\alpha_{a}^{2}}{2}\right)\,.$ (3.8)

Importantly for the subsequent discussion, the random phase $\theta_{a}$ and
the Rayleigh variable $\alpha_{a}$ are independent [53].

In Eq. (3.6) we neglected the spatial variation of the ULDM wave. We justify
this by noting that on average the wave vector is suppressed by a factor of
$v\sim 10^{-3}$ relative to the angular frequency. Furthermore, the lengths of
the baselines that we consider are significantly smaller than the ULDM’s de
Broglie wavelength, such that any spatially-dependent contribution to the
phase of the DM wave is highly subdominant. Indeed, the frequency range of
interest corresponds to ULDM masses between
${10}^{-17}\text{\,}\mathrm{e}\mathrm{V}$ and
${10}^{-10}\text{\,}\mathrm{e}\mathrm{V}$, which would be associated with de
Broglie wavelengths of length ${10}^{11}\text{\,}\mathrm{k}\mathrm{m}$ and
${10}^{4}\text{\,}\mathrm{k}\mathrm{m}$, respectively.

The phase shift induced by linearly-coupled scalar ULDM for the gradiometer
$\text{AG-}(i,j)$ where the two AIs are separated by a distance $\Delta
z^{(i,j)}=z_{i}-z_{j}$ can be parametrised as

$\Phi^{(i,j)}_{\mathrm{DM},m}=\frac{\Delta
z^{(i,j)}}{L}\sum_{a}\alpha_{a}\sqrt{F_{\mathrm{DM}}(v_{a})}\,A_{a}\cos\phi_{a,m}\;.$
(3.9)

Here, $\phi_{a,m}\supset\omega_{a}m\Delta t+\theta_{a}$ is the phase of the DM
wave at the end of the sequence. The amplitude $A_{a}$ is dependent on the
experimental variables $n$, $T$, and $L$ characterising the interferometric
sequence, and on phenomenological parameters, the most important of which is
the ULDM-Standard Model coupling strength $d_{\phi}$. In the limit $nL\ll T$
and $\omega_{a}nL\ll 1$, which are always satisfied in planned interferometers
[14, 15],

$A_{a}\propto\frac{d_{\phi}\,n\,L}{\omega_{a}}\sin^{2}\left(\frac{\omega_{a}T}{2}\right)\,.$
(3.10)

The complete expressions for $A_{a}$ and $\omega_{a}$ are given in Appendix A.

Since $\alpha_{a}$ and $\theta_{a}$ are independent, over long timescales the
signal is characterised by a vanishing expectation value,

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Big{\rangle}$
$\displaystyle=0\,,$ (3.11)

and a covariance given by

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{DM},m^{\prime}}\Big{\rangle}$
$\displaystyle=\frac{\Delta z^{(i,j)}}{L}\frac{\Delta
z^{(i^{\prime},j^{\prime})}}{L}\sum_{a}F_{\mathrm{DM}}(v_{a})\,A_{a}^{2}\cos(\omega_{a}\Delta
t(m-m^{\prime}))\,.$ (3.12)

Both of these properties follow from the statistical properties of the
amplitude and random phase. A complete derivation of Eqs. (3.11) and (3.12) is
provided in Appendix A.1.

###  3.2 Gravity gradient noise from seismic effects

Time-dependent perturbations to the gravitational potential around the free-
falling atom clouds cause a phase shift to build up between two paths of an
atom interferometer. These perturbations in the gravitational potential
originate from terrestrial density fluctuations close to each atom
interferometer, as illustrated in Fig. 2, and give rise to a noise source
commonly referred to as gravity gradient noise (GGN). For the experimental
configuration and frequency range considered in this paper, the dominant
source of GGN is expected to consist of ground density perturbations induced
by ambient horizontally-propagating seismic waves that are confined near the
Earth’s surface by horizontal geological strata [25, 24]. The modes
responsible for these density fluctuations are referred to as Rayleigh modes,
and they are generated at horizontal discontinuities at strata interfaces,
including the Earth’s surface.

In this work, we employ the homogeneous half space model for the ground around
the interferometer. This model implies that the ground is isotropic around the
detector and that only the fundamental Rayleigh mode is present. The
fundamental mode travels horizontally at the Earth-air interface and is
present irrespective of the geological properties of the site. Following the
GGN treatments for LIGO [25, 24], we represent the seismic field as an
incoherent superposition of monochromatic plane-waves propagating
isotropically at the Earth’s surface. The amplitude of each wave depends on
the vertical displacement $\xi_{a}$ at the Earth’s surface, which depends on
the wave frequency, and is characterised by an exponentially decaying
amplitude while moving away from the surface in either direction. Here, we
bound the upper and lower values of the vertical displacement using Peterson’s
new high noise model (NHNM) and new low noise model (NLNM) [54]. In Fig. 3, we
display the vertical displacement spectrum $\sqrt{S_{\xi,a}}$ according to
Peterson’s model, where
${S_{\xi,a}}\propto\left\langle\xi_{a}^{2}\right\rangle$ is the power spectral
density of the vertical displacement (see section 4 and appendices referenced
therein for more details). We limit our analysis to frequencies greater than
$10^{-3}$ Hz since below this frequency the Earth’s finite size is expected to
significantly alter the properties of the seismic background and render our
analysis based on Rayleigh waves invalid [55, 24, 56].

Figure 2: Schematic representation of the seismic field for a fundamental
Rayleigh mode with angular frequency $\omega_{a}$ travelling horizontally in
the direction $\hat{\mathbf{k}}$. The Rayleigh wave will induce underground
density variations $\delta\rho$ at positions labelled by $\mathbf{r}$. The
density variations induce perturbations in the gravitational potential, which
serve as a noise source for the $i^{\mathrm{th}}$ interferometer AI-$i$.

Before moving on to present the results for the induced phase shift, we
briefly consider other potential sources that could induce time-dependent
perturbations to the gravitational potential and thereby contribute to the
gravity gradient noise. Firstly, we do not include the effect from atmospheric
density fluctuations, which dominate over the Rayleigh wave-induced
contribution below $\sim 0.1$ Hz in horizontal atom interferometer
configurations [30]. Secondly, seismic waves could scatter from the shaft’s
surfaces, but these effects are negligible for cavity radii much smaller than
the wavelength of the seismic waves [24]. For the frequency range between
$10^{-3}$ Hz and 1 Hz, the wavelength of these perturbations falls between
$\sim$40\text{\,}\mathrm{m}$$ and $\sim${10}^{4}\text{\,}\mathrm{m}$$, while
the expected radial size of the access shaft is
$\mathcal{O}($5\text{\,}\mathrm{m}$)$. Thirdly, we also neglect the direct
effect induced by vibrations on the laser system, since we expect that these
vibrations can be mitigated with suspensions and other techniques widely used
in long-baseline interferometry [15]. Finally, anthropic sources such as heavy
vehicles could play a role; owing to their intrinsic site-specific dependence,
they will be considered in a separate study [57].

Figure 3: Peterson’s models of the frequency-dependent spectrum of the
vertical displacement $\xi_{a}$ at the surface of the Earth between
${10}^{-3}\text{\,}\mathrm{H}\mathrm{z}$ and
$10\text{\,}\mathrm{H}\mathrm{z}$. The orange curve refers to the prediction
from the new high noise model (NHNM), while the blue curve is from the new low
noise model (NLNM).

Returning to the assumption that the seismic field arises from an incoherent
superposition of fundamental Rayleigh modes, we obtain the angle-averaged
seismically-induced gravitational potential experienced by a test mass
underground at time $t$ and position $\vec{r_{i}}=(0,0,z_{i})$ by integrating
over all ground density fluctuations around the experiment (recall that
$z_{i}$ is the vertical position from the top of the interferometer shaft). In
Appendix B, we derive the result that the potential due to seismic
fluctuations takes the form

$V\left(z_{i},t\right)\propto\sum_{a}\xi_{a}\,\left[-2\exp\left(-q\frac{\omega_{a}z_{i}}{c_{H}}\right)+\left(1+\sqrt{q/s}\right)\exp\left(-\frac{\omega_{a}z_{i}}{c_{H}}\right)\right]\,\cos\left(\omega_{a}t+\widetilde{\theta}_{a}\right)\,,$
(3.13)

where $q$ and $s$ are $\mathcal{O}(1)$ dimensionless parameters that depend
exclusively on the geological properties of the ground, $c_{H}$ is the
horizontal speed of the Rayleigh wave, $\widetilde{\theta_{a}}$ is a random
phase and the sum is over discrete frequencies $\omega_{a}$ as in the ULDM
case (see Eq. (3.6)). With this potential, we find that the GGN phase shift
recorded by the $i$th AI at the end of a sequence that starts at time $m\Delta
t$ can be expressed as

$\Phi^{(i)}_{\mathrm{GGN},m}=\sum_{a}\xi_{a}\left[\widetilde{A}_{a}\exp\left(-q\frac{\omega_{a}z_{i}}{c_{H}}\right)+\widetilde{B}_{a}\exp\left(-\frac{\omega_{a}z_{i}}{c_{H}}\right)\right]\cos{\widetilde{\phi}_{a,m}}\;,$
(3.14)

where $\widetilde{A}_{a}$ and $\widetilde{B}_{a}$ depend on experimental and
geological parameters, and $\widetilde{\phi}_{a,m}\supset\omega_{a}m\Delta
t+\widetilde{\theta}_{a}$ is the phase of the wave at the end of the sequence.
In Appendix B we provide a full derivation of this result, together with
complete expressions for $\widetilde{A}_{a}$ and $\widetilde{B}_{a}$. Up to a
coefficient,

$\widetilde{A}_{a},\widetilde{B}_{a}\propto\frac{n}{\omega_{a}^{2}}\,\,\sin^{2}\left(\frac{\omega_{a}T}{2}\right)\,,$
(3.15)

which implies that the amplitude of the GGN phase, as for the ULDM phase
shift, is a function of $n$ and $T$.

Unlike ULDM, the GGN phase shift does not depend on the length of the baseline
$L$; rather, it depends on the vertical position of the AIs along the baseline
and on the length scale

$\lambda_{\rm{GGN}}=\frac{c_{H}}{\omega_{a}}\simeq$100\text{\,}\mathrm{m}$~{}\left(\frac{250~{}\mathrm{m\,s^{-1}}}{c_{H}}\right)^{-1}~{}\left(\frac{2.5~{}\mathrm{Hz}}{\omega_{a}}\right)\;.$
(3.16)

The values of $c_{H}$, $s$ and $q$ depend on the precise makeup of the ground
where the interferometer is situated. For example, Ref. [25] quotes a range of
$200\,\mathrm{ms}^{-1}\lesssim c_{H}\lesssim 6000\,\mathrm{ms}^{-1}$ depending
on the rock and mineral composition. As we do not have a specific site in mind
for this work, we focus on geological parameter values that will most clearly
highlight the impact of atom multi-gradiometry and the benefits of the
likelihood formalism to potentially extend the sensitivity in the sub-Hz
frequency range. This corresponds to the lower range of $c_{H}$ where
$\lambda_{\rm{GGN}}\lesssim L$ in the frequency ($f$) range between 0.1 Hz and
1 Hz (corresponding to the angular frequency ($\omega$) range between 0.6 Hz
and 6 Hz). Specifically, our default set of parameters correspond to the
surface materials at the LIGO Livingston site [25], where the ground density
is 1800 kg m-3, $c_{H}=205\,\mathrm{m\,s}^{-1}$, $q=0.88$ and $s=0.36$. In
section 5, we will also briefly consider the scenario where $c_{H}$ is over
$1000\,\mathrm{m\,s}^{-1}$ such that $\lambda_{\rm{GGN}}>L$.

The random variables $\widetilde{\theta}_{a}$ and $\xi_{a}$ that enter Eq.
(3.14) are uncorrelated random variables. As for the ULDM case, this means
that over long timescales the GGN phase measured by a gradiometer will be
characterised by a vanishing expectation value,

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{GGN},m}\Big{\rangle}$
$\displaystyle=0\,,$ (3.17)

and, assuming sufficiently similar geological conditions along the baseline,
covariance given by

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{GGN},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{GGN},m^{\prime}}\Big{\rangle}$
$\displaystyle=\sum_{a}\frac{\langle\xi_{a}^{2}\rangle}{2}F_{\mathrm{GGN},a}^{(i,j)}F_{\mathrm{GGN},a}^{(i^{\prime},j^{\prime})}\cos(\omega_{a}\Delta
t(m-m^{\prime}))\,,$ (3.18)

where

$\displaystyle
F_{\mathrm{GGN},a}^{(i,j)}=F_{\mathrm{GGN},a}^{(i)}-F_{\mathrm{GGN},a}^{(j)}\,,$
(3.19) $\displaystyle
F_{\mathrm{GGN},a}^{(i)}=\widetilde{A}_{a}\exp\left(-q\frac{\omega_{a}z_{i}}{c_{H}}\right)+\widetilde{B}_{a}\exp\left(-\frac{\omega_{a}z_{i}}{c_{H}}\right)\,.$
(3.20)

###  3.3 Atom shot noise

The measurement of the phase recorded by an AI is ultimately limited by atom
shot noise. Indeed, as the phase shift is extracted from the comparison of
atoms in the ground versus excited state at the end of a sequence, the
Poissonian statistics governing atom number affect the interferometer phase.

Formally, the atom shot noise phase of the $i^{\mathrm{th}}$ interferometer
$\Phi_{\mathrm{ASN}}^{(i)}$ is independently sampled at each measurement from
a normal distribution with zero mean and variance
$\sigma_{\mathrm{Atom}}^{2}=(C^{2}N_{\mathrm{Atom}})^{-1}$ for phase
differences close to $\pi/2$, where $N_{\mathrm{Atom}}$ is the number of atoms
in the cloud and $C\leq 1$ is the interferometer contrast [38, 39, 13]. For
simplicity, in our analysis we set $C=1$. Because the variance
$\sigma_{\mathrm{Atom}}^{2}$ and mean are frequency-independent, and given
that the noise phase is independent in each measurement, the phase
contribution from atom shot noise can be modeled as white noise. Hence, the
atom shot noise phase statistics can be summarised as

$\displaystyle\langle\Phi^{(i)}_{\mathrm{ASN},m}\rangle$ $\displaystyle=0\,,$
(3.21)
$\displaystyle\langle\Phi^{(i)}_{\mathrm{ASN},m}\Phi^{(j)}_{\mathrm{ASN},m^{\prime}}\rangle$
$\displaystyle=\delta_{mm^{\prime}}\delta^{ij}\sigma_{\mathrm{Atom}}^{2}\,.$
(3.22)

The atom shot noise contribution to the gradiometer phase will inherit the
statistical properties of the phases measured by each interferometer.
Recalling that the atom shot noise phase contribution for the ($i,j$)th
gradiometer is defined as
$\Phi^{(i,j)}_{\mathrm{ASN}}=\Phi^{(i)}_{\mathrm{ASN}}-\Phi^{(j)}_{\mathrm{ASN}}$
and using Eqs. (3.21)-(3.22), in the limit of long integration time, the atom
shot noise phase will average to zero,

$\displaystyle\langle\Phi^{(i,j)}_{\mathrm{ASN},m}\rangle$ $\displaystyle=0\,$
(3.23)

and is characterised by the covariance

$\displaystyle\langle\Phi^{(i,j)}_{\mathrm{ASN},m}\Phi^{(p,q)}_{\mathrm{ASN},m^{\prime}}\rangle$
$\displaystyle=\delta_{mm^{\prime}}(\delta^{ip}+\delta^{jq}-\delta^{iq}-\delta^{jp})\sigma_{\mathrm{Atom}}^{2}\,.$
(3.24)

Unlike the ULDM and GGN phases, Eq. (3.24) can be exactly zero when comparing
four different interferometers. Hence, in this paper, we expect zero
correlations between some gradiometer pairs to be a unique feature of atom
shot noise.

##  4 Analysis of an atom multi-gradiometer

Having characterised the signal and noise correlations of our model, in this
section we will present a multi-gradiometer analysis scheme for DM searches to
extract additional information available to two or more gradiometers, which
will be of crucial importance to predict and improve an experiment’s
sensitivity at low frequencies. Our approach makes use of the machinery
developed in Refs. [53, 58] for ALP searches.

Because a DM signal is characterised by a frequency largely set by its mass,
and by a frequency spread that is dictated by its speed distribution, we will
present all of the subsequent analysis in the frequency domain. Therefore, the
appropriate tool for analysing the data is the discrete Fourier transform
(DFT) $\widetilde{\Phi}_{k}^{(i,j)}$ of the phase measured by the $(i,j)$th
gradiometer, which is defined as

$\widetilde{\Phi}_{k}^{(i,j)}=\sum_{m=0}^{N-1}\Phi_{m}^{(i,j)}\exp{\left(-\frac{2\pi
ikm}{N}\right)}\;,$ (4.1)

where $k\in\\{0,1,...,N-1\\}$. Following Ref. [58], we decompose the DFT in
terms of real and imaginary parts, namely

$\displaystyle R_{k}^{(i,j)}$ $\displaystyle=\frac{\Delta
t}{\sqrt{T_{\mathrm{int}}}}\operatorname{Re}\left[\widetilde{\Phi}_{k}^{(i,j)}\right]\,,$
(4.2) $\displaystyle I_{k}^{(i,j)}$ $\displaystyle=\frac{\Delta
t}{\sqrt{T_{\mathrm{int}}}}\operatorname{Im}\left[\widetilde{\Phi}_{k}^{(i,j)}\right]\,,$
(4.3)

such that the one-sided power spectral density (PSD) $S_{k}^{(i,j)}$ of the
time-series may be written in terms of this decomposition as

$S_{k}^{(i,j)}\equiv\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\left|\widetilde{\Phi}_{k}^{(i,j)}\right|^{2}=\left|R_{k}^{(i,j)}\right|^{2}+\left|I_{k}^{(i,j)}\right|^{2}\,.$
(4.4)

The information in the $k$-frequency component collected by all gradiometers
can then be organized into a $\mathcal{N}(\mathcal{N}-1)$-dimensional data
vector,

$\mathbf{d}_{k}=\Big{[}R_{k}^{(1,2)},I_{k}^{(1,2)},\ldots,R_{k}^{(\mathcal{N},\,\mathcal{N}-1)},I_{k}^{(\mathcal{N},\,\mathcal{N}-1)}\Big{]}^{\mathrm{\mathbf{T}}}\,.$
(4.5)

Since the phase measured by each interferometer, and hence gradiometer, has an
expectation value of zero (cf. Eqs. (3.11), (3.17) and (3.21)) and is Gaussian
distributed,444We quote this result from Refs. [53, 58], in which the authors
numerically show that white noise and Rayleigh-distributed data is described
by a multi-variable Gaussian distribution. the data vectors will be
characterised by a symmetric
$(\mathcal{N}(\mathcal{N}-1)\times\mathcal{N}(\mathcal{N}-1))$-dimensional
covariance matrix
$\mathbf{\Sigma}_{k}=\left\langle\mathbf{d}_{k}\mathbf{d}^{\mathrm{\mathbf{T}}}_{k}\right\rangle$
whose entries will depend on both the noise and signal covariances presented
in the previous section. Because the noise and signal are uncorrelated, we can
use the statistical properties of the different gradiometer phases to
decompose the covariance matrix $\boldsymbol{\mathrm{\Sigma}}_{k}$ into signal
$\boldsymbol{\mathrm{S}}_{k}$ and noise $\boldsymbol{\mathrm{B}}_{k}$
contributions. In Appendix C we further elucidate this point and provide
detailed analytical calculations for the entries of
$\boldsymbol{\mathrm{\Sigma}}_{k}$ and $\boldsymbol{\mathrm{S}}_{k}$ given the
signal and background models used in this work.

By defining a model $\mathcal{M}$ with parameter vector $\boldsymbol{\theta}$
that has nuisance parameters $\boldsymbol{\theta}_{\mathrm{nuis}}$, which
describe GGN and ASN in the individual gradiometers, and signal parameters
$\boldsymbol{\theta}_{\mathrm{sig}}$, which characterize the ULDM signal
contribution, the likelihood is given by

$\mathcal{L}(d\mid\mathcal{M},\boldsymbol{\theta})=\prod_{k=1}^{N-1}\frac{\exp\left[-\frac{1}{2}\mathbf{d}_{k}^{T}\cdot\boldsymbol{\Sigma}_{k}^{-1}(\boldsymbol{\theta})\cdot\mathbf{d}_{k}\right]}{\sqrt{(2\pi)^{2\mathcal{N}}\left|\boldsymbol{\Sigma}_{k}(\boldsymbol{\theta})\right|}}\,,$
(4.6)

where $d=\\{\mathbf{d}_{k}\\}$ is the set of data vectors at each
$k$-frequency, $|\boldsymbol{\Sigma}_{k}(\boldsymbol{\theta})|$ is the
determinant of the covariance matrix and
$\boldsymbol{\theta}=\\{\boldsymbol{\theta}_{\mathrm{sig}},\boldsymbol{\theta}_{\mathrm{nuis}}\\}$.
In our definition of the likelihood, we omitted the zeroth frequency
component, since the DFT of the data at $k=0$ would be polluted by a tower of
static effects, while the ULDM signals considered here are time-dependent.
Additionally, the mass of a hypothetical ULDM candidate in this frequency bin
would be degenerate with any value below $2\pi/T_{\mathrm{int}}$, thus leading
to poor mass resolution.

With this likelihood, we define the following frequentist tool based on the
log profile likelihood

$\Theta\left(\boldsymbol{\theta}_{\mathrm{sig}}\right)=2\left[\ln\mathcal{L}\left(d\mid\mathcal{M},\left\\{\hat{\boldsymbol{\theta}}_{\text{nuis}},\boldsymbol{\theta}_{\mathrm{sig}}\right\\}\right)\right.\left.-\ln\mathcal{L}\left(d\mid\mathcal{M},\left\\{\hat{\hat{\boldsymbol{\theta}}}_{\text{nuis}},\boldsymbol{\theta}_{\mathrm{sig}}=\mathbf{0}\right\\}\right)\right]\,,$
(4.7)

Here, $\hat{\boldsymbol{\theta}}_{\text{nuis}}$ denotes the values of the
nuisance parameters that maximise the likelihood in the signal plus background
hypothesis, while $\hat{\hat{\boldsymbol{\theta}}}_{\text{nuis}}$ represents
the values of the nuisance parameters that maximise the likelihood in the
background-only hypothesis.

To set upper limits on couplings to scalar ULDM, we set
$\boldsymbol{\theta}_{\mathrm{sig}}=\left(m_{\phi},d_{\phi}\right)$, where we
remind the reader that $m_{\phi}$ is the dark matter mass and $d_{\phi}$ is
the coupling strength of the relevant linear interaction between ULDM and
Standard Model (SM) operators. Then, Eq. (4.7) can be used to define a test
statistic for setting upper limits on the ULDM-SM couplings $d_{\phi}$,

$q\left(m_{\phi},d_{\phi}\right)=\begin{cases}\Theta\left(\\{m_{\phi},d_{\phi}\right\\})-\Theta\left(\\{m_{\phi},\hat{d_{\phi}}\\}\right)&d_{\phi}\geq\hat{d_{\phi}}\\\
0&d_{\phi}<\hat{d_{\phi}}\end{cases}$ (4.8)

where $\hat{d_{\phi}}$ is the value of $d_{\phi}$ that maximises
$\Theta\left(m_{\phi},d_{\phi}\right)$ at fixed $m_{\phi}$.555In light of Eq.
(4.7), this test statistic is independent of the background-only hypothesis,
which is in agreement with the literature [59]. With this definition, the test
statistic $q$ at fixed $m_{\phi}$ is described by a half chi-squared
distribution with one degree of freedom [53, 59], such that for a given
$m_{\phi}$, the 95% limit on $d_{\phi}$ is set when
$q(m_{\phi},d_{\phi,95\%})\approx-2.70$ when $T_{\mathrm{int}}\gg\tau_{c}$ and
$q(m_{\phi},d_{\phi,95\%})\approx-7.55$ when $T_{\mathrm{int}}\ll\tau_{c}$,
where $\tau_{c}$ is the coherence time of the ULDM wave [60].666To make
contact with the signal-to-noise ratio (SNR) defined in Refs. [14, 13, 10],
please note that $\mathrm{SNR}=\sqrt{-q(m_{\phi},d_{\phi})}$. The regime where
$T_{\mathrm{int}}\sim\tau_{c}$ is beyond the scope of this work.

###  4.1 Test Statistics in the Asimov Approach

To infer the expected distribution of the test statistics and understand the
expected sensitivity of the multi-gradiometer experiment, many realisations of
the experiment can be generated via Monte Carlo simulations. Alternatively, it
is also possible to infer the asymptotic properties of the test statistics
analytically by taking the dataset to be equal to the mean predictions of the
model under consideration and neglecting statistical fluctuations, which is
known as the Asimov approach [53, 58]. Since the test statistic for setting
upper limits depends on the modified log profile likelihood, we define the
basic frequentist tool in the Asimov approach as

$\widetilde{\Theta}(\boldsymbol{\theta}_{\mathrm{sig}})=\langle\Theta(\boldsymbol{\theta}_{\mathrm{sig}})\rangle\,.$
(4.9)

Recalling that we have assumed a covariance matrix
$\boldsymbol{\mathrm{\Sigma}}_{k}$ that can be split into signal
$\boldsymbol{\mathrm{S}}_{k}$ and noise $\boldsymbol{\mathrm{B}}_{k}$
contributions, and assuming that the signal is parametrically weaker than the
background contributions, we may express the Asimov test statistic for setting
upper limits as:

$\widetilde{q}\left(m_{\phi},d_{\phi}\right)=\begin{cases}-\frac{1}{2}\sum_{k=1}^{N-1}\mathrm{Tr}\left[\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\right]&d_{\phi}\geq
0\\\ 0&d_{\phi}<0\end{cases}\,,$ (4.10)

which scales with $T_{\mathrm{int}}^{2}$ when $T_{\mathrm{int}}\gg\tau_{c}$,
and with $\tau_{c}\,T_{\mathrm{int}}$ when $T_{\mathrm{int}}\ll\tau_{c}$. In
Appendix D we provide a derivation of Eq. (4.10), and explain the origin of
the behaviour of the test statistic in the two aforementioned regimes. Since
the non-vanishing entries of $\mathbf{S}_{k}$ are proportional to
$d_{\phi}^{2}$ and those of $\mathbf{B}_{k}$, if assuming atom shot noise
only, are proportional to $1/N_{\mathrm{atom}}$, $\widetilde{q}$ is
proportional to $d_{\phi}^{4}\,N_{\mathrm{atom}}^{2}$. This is in agreement
with results in Ref. [10] in which $d_{\phi}\propto
1/\sqrt{N_{\mathrm{atom}}}$.

##  5 ULDM searches in a GGN background

Design |  $L$ [m] |  $T$ [s] | $n$ |  $\Delta z_{\mathrm{max}}$ [m] | $N_{\mathrm{atom}}$ |  $\Delta t$ [s] |  $T_{\mathrm{int}}$ [s]
---|---|---|---|---|---|---|---
Intermediate | $100$ | 1.4 | $1000$ | 85 | $10^{8}$ | 1.5 | $10^{8}$
Advanced | $1000$ | 1.7 | $2500$ | $970$ | $10^{10}$ | 1.0 | $10^{8}$
Table 1: The experimental parameters for the different atom interferometer
scenarios considered in this work. The ‘intermediate’ and ‘advanced’ designs
could be implemented in versions of future vertical gradiometers, such as
AION-100 and MAGIS-100, and AION-km, respectively. $\Delta z_{\mathrm{max}}$
refers to the maximum gradiometer length given the choice of interferometer
parameters. We also consider scenarios where $\Delta z$ is shorter than the
maximum value.

In order to assess the impact of GGN on a gradiometer’s maximum sensitivity,
it is of paramount importance to understand the relative size of GGN with
respect to atom shot noise. Indeed, it will be precisely this comparison that
will guide our analysis into understanding the benefits of a multi-gradiometer
experiment to mitigate and, in an ideal scenario, entirely cancel GGN at
particular frequencies. Additionally, because the GGN background as detected
by an atom interferometer is largely dependent on experimentally tunable
parameters, we will provide a discussion in terms of two vertical gradiometer
configurations: an ‘intermediate’ and an ‘advanced’ design, which feature a
$100\text{\,}\mathrm{m}$ and $1\text{\,}\mathrm{k}\mathrm{m}$ baseline,
respectively. In Table 1 we summarise all of the relevant parameters that will
be necessary for subsequent calculations and discussions. Additionally, we
will restrict our discussion to the impact of GGN on the sensitivity to scalar
ULDM couplings to the electron mass, i.e. $d_{\phi}=d_{m_{e}}$. We will begin
our discussion with a single atom gradiometer before turning to consider the
advantages brought by a multi-atom gradiometer in mitigating the GGN
background.

###  5.1 Single atom gradiometer

In the left panel of Fig. 4 we plot the single atom gradiometer
($\mathcal{N}=2$) noise power spectral density (PSD) as a function of
frequency for the ‘intermediate’ design defined in Table 1. The atom shot
noise (ASN) contribution is shown with a dashed grey line and the gravity
gradient noise (GGN) contribution according to the new high (NHNM) and low
noise models (NLNM) are shown with orange and blue solid lines, respectively.
We see that the PSD of the gradiometer phase induced by GGN dominates the low
frequency range below $\sim 0.5$ Hz for the NHNM. This indicates that in an
environment where the seismic noise is high, GGN can limit an experiment’s
reach in ULDM parameter space for
$L\sim\mathcal{O}{($100\text{\,}\mathrm{m}$)}$.777In the frequency range of
interest, GGN does not limit the sensitivity of $L\sim$10\text{\,}\mathrm{m}$$
baselines, even in the NHNM scenario.

Figure 4:  Impact of GGN on a single atom gradiometer ($\mathcal{N}=2$) for
the ‘intermediate’ design defined in Table 1. Left panel: Contributions to the
gradiometer noise PSD as a function of frequency from atom shot noise (dashed
grey) and GGN according to the NHNM (orange) and NLNM (blue) scenarios. Right
panel: Projected 95% CL exclusion curves on the ULDM-electron coupling when
assuming an atom shot noise-only background (grey dashed), and when GGN
modelled by the NHNM is also included (orange). The NLNM result is degenerate
with the ASN line so is omitted for clarity. The orange shaded region is
excluded by MICROSCOPE [23]. In high seismic-noise scenarios, GGN can limit an
‘intermediate’ design’s reach in ULDM parameter space.

We demonstrate the impact of the NHNM model explicitly in the right panel of
Fig. 4, which shows the sensitivity projections in the ULDM frequency (or
mass) versus ULDM-electron coupling plane for the ‘intermediate’ design.
Typically, the sensitivity projections oscillate as a function of the ULDM
mass. However, for clarity, we follow Refs. [12, 13] and plot only the
envelope of the oscillations by employing the approximation $|\sin
x\,|=\mathrm{min}\\{x,1/\sqrt{2}\\}$ in the phase formulas (e.g., Eqs. (A.5)
and (B.20)). The dashed grey line shows the projected 95% CL exclusion curve
in the presence of a background that only contains atom shot noise. In
contrast, the solid orange line shows the 95% CL exclusion curve when GGN
modelled by the NHNM is included. As anticipated from the left panel, the
sensitivity projections are reduced below $\sim 0.5$ Hz. For clarity the right
panel does not show the projection for the NLNM scenario. As may be
anticipated from the left panel, the sensitivity projection for the NLNM is
degenerate with the grey line as it is sub-dominant to atom shot noise for all
frequencies.

Next we turn to the ‘advanced’ design, in which the baseline is 1 km. The left
panel of Fig. 5 again shows the noise PSD for a single atom gradiometer
($\mathcal{N}=2$). For this design, we see that GGN exceeds atom shot noise in
the sub-Hz range. In the NHNM (orange) it is the dominant noise source below
$\sim 1$ Hz while in the NLNM (blue), it dominates below $\sim 0.5$ Hz. The
right panel of Fig. 5 shows the impact of GGN in the ULDM mass versus ULDM-
electron coupling plane. In the sub-Hz regime, in the NLNM scenario the 95% CL
exclusion curve (blue) is weaker by approximately an order of magnitude with
respect to a naive estimate accounting for atom shot noise only. The loss of
sensitivity in the NHNM is much more dramatic, with a suppression by up to
four orders of magnitude with respect to the atom shot noise only result.
Hence, we conclude that the sensitivity of a single atom gradiometer
($\mathcal{N}=2$) of baseline length
$L\gtrsim\mathcal{O}{($1\text{\,}\mathrm{k}\mathrm{m}$)}$ would be greatly
diminished across the entire low frequency range in both the NHNM and NLNM.

Figure 5:  Impact of GGN on a single atom gradiometer ($\mathcal{N}=2$) for
the ‘advanced’ design defined in Table 1. Right panel: Projected 95% CL
exclusion curve on the ULDM-electron coupling with an atom shot noise-only
background (grey), and when GGN according NHNM (orange) and NLNM (blue) are
also included. The orange shaded region is excluded by MICROSCOPE [23]. GGN
can greatly diminish an ‘advanced’ design’s reach in ULDM parameter space.

Even within the limits of a single atom gradiometer ($\mathcal{N}=2$), there
is scope to regain parts of parameter space by tuning the experimental
parameters to minimise the impact of GGN. We observe that the GGN and the ULDM
signal in an atom interferometer both depend identically on the number of LMT
kicks $n$ and interrogation time $T$ (cf. Eqs. (3.10) and (3.15)), but show
different scaling relationships with respect to the vertical position along
the baseline. Hence, we will focus on the impact of changing the position of
the interferometers along the baseline.

Both panels in Fig. 6 show the 95% CL exclusion curves for the ‘advanced’
design; however, in these panels, we exploit the freedom to locate the two
atom interferometers at different positions along the 1 km baseline. In both
panels the orange line repeats the result from Fig. 5 for the NHNM GGN
scenario in which the distance $\Delta z$ between the two interferometers is
set to the maximum value, $\Delta z_{\rm{max}}=970$ m. The left panel of Fig.
6 additionally shows the 95% CL exclusion curves when the atom interferometers
are separated by $\Delta z=100$ m and are situated closer to the Earth’s
surface (green) or towards the bottom of the shaft (purple). While situating
the atom interferometers towards the Earth’s surface brings no benefit
relative to $\Delta z_{\rm{max}}$, we observe a sensitivity improvement for
the configuration at the bottom of the shaft below $\sim 0.5$ Hz. However,
this comes with the cost of losing sensitivity at frequencies greater than
$\sim 0.5$ Hz. The right panel of Fig. 6 considers the scenario where $\Delta
z=600$ m, and the green and purple lines again show the exclusion curves when
positioning one of the interferometers closer to the Earth’s surface and
towards the bottom of the shaft, respectively. In this case, a similar
behaviour is observed: the configuration with an AI at the bottom of the shaft
is favoured below $\sim 0.7$ Hz but again at the cost of reduced sensitivity
at higher frequencies. Comparing the left and right panels of Fig. 6 we see
that in a narrow frequency range between $\sim 0.2$ Hz and $\sim 0.7$ Hz,
there is a preference for larger $\Delta z$ values at the bottom of the shaft.

Figure 6: Impact on the 95% CL exclusion curves for different values of
$\Delta z$ and different atom interferometer positions for a single atom
gradiometer ($\mathcal{N}=2$) assuming that GGN is modelled by the NHNM. The
orange lines in both panels show the ‘advanced’ design parameters
configuration. We also show exclusion curves with interferometers located
towards the Earth’s surface (green) and towards the bottom of the shaft
(purple) assuming $\Delta z=100$ m (left panel) and $\Delta z=600$ m (right
panel) and with all other ‘advanced’ parameters unchanged. The orange shaded
region is excluded by MICROSCOPE [23]. Placing atom interferometers towards
the bottom of the shaft regains parts of parameter space below $\sim 0.5$ Hz
but at the cost of reduced sensitivity at higher frequencies.

The frequency-dependent behaviour observed in Fig. 6 can be qualitatively
understood by considering the dependence of the GGN and ULDM phase shifts on
the interferometer position along the baseline with respect to the Rayleigh
wave’s profile underground. At fixed $\Delta z$, the GGN background and ULDM
signal scale differently in the high and low frequency limit. At low
frequencies, corresponding to $\omega_{a}\ll 1/T$ and
$\lambda_{\rm{GGN}}\gg\Delta z$, the GGN and ULDM phase shifts both scale like
$\omega_{a}\Delta z$ (cf. Eqs. (3.14) and (3.9)). Therefore, given the
geological parameters values used in this work, little to no sensitivity
enhancement is gained from probing different gradiometer length scales for
angular frequencies $\omega_{a}\ll 10^{-1}$ Hz.

In contrast, in the frequency window between $\sim 0.05$ Hz and $\sim 0.2$ Hz,
where GGN dominates the background, and when $\lambda_{\rm{GGN}}\lesssim\Delta
z_{\mathrm{max}}$, the ULDM and GGN phase shifts scale very differently: as
$\Delta z/\omega_{a}$ for ULDM; and as
$|\exp(-\,\omega_{a}z_{i}/c_{H})-\exp(-\omega_{a}z_{j}/c_{H})|/\omega_{a}^{2}$
for GGN. Writing $z_{j}=z_{i}+\Delta z$, the ratio between the amplitudes of
the ULDM and GGN phase shifts is maximised for $z_{i}=L$ and $\Delta z$ small,
which implies that the GGN background is better mitigated when placing the
interferometers at the bottom of the baseline and choosing small gradiometer
separations. At high frequencies where the background is dominated by atom
shot noise, which is frequency-independent, the maximum reach can be achieved
by choosing the largest possible gradiometer length, $\Delta
z_{\mathrm{max}}$.

Finally, the preference for larger $\Delta z<\Delta z_{\mathrm{max}}$ values
at the bottom of the shaft between $\sim 0.2$ Hz and $\sim 0.7$ Hz arises as a
result of the fact that at these frequencies there exist positions along the
baseline where GGN is subdominant. This can be seen in Fig. 7, where we show
the PSD of GGN over frequencies for long and short values of $\Delta z<\Delta
z_{\mathrm{max}}$ when one of the interferometers is located at 1000 m below
the Earth’s surface. For both $\Delta z=600$ m and $\Delta z=100$ m, the
experiment will be ASN-limited above $\sim 0.2$ Hz, which implies that the
experimental sensitivity will scale like $\Delta z$ above these frequencies.
Hence, between $\sim 0.4$ Hz and $\sim 1$ Hz, the configuration with $\Delta
z=600$ m is favoured, as shown in Fig. 6.

Figure 7: Contributions to the gradiometer noise PSD from atom shot noise
(dashed grey) and GGN according to the NHNM for an ‘advanced’ single atom
gradiometer ($\mathcal{N}=2$) in which the AIs are positioned at: 900 m and
1000 m ($\Delta z=100$ m); and 400 m and 1000 m ($\Delta z=600$ m). We also
show the PSD of GGN according to the NHNM for $\Delta z_{\mathrm{max}}=970$ m,
which is identical to orange curve in the left panel of Fig. 5.

###  5.2 A multi-atom gradiometer experiment

In our discussion of a single atom gradiometer, we observed that GGN can
greatly diminish the reach in ULDM parameter space. Furthermore, we observed
that the impact of GGN could be changed by varying the position of the atom
interferometers and the distance $\Delta z$ between them. For instance, the
configuration with $\Delta z=100$ m where the atom interferometers are
positioned towards the bottom of the 1 km shaft improved the sensitivity in
the sub-Hz regime but at the cost of a loss in sensitivity above 1 Hz compared
to the $\Delta z=970$ m configuration. With these observations in mind, we can
immediately see the benefit of considering a multi-gradiometer experiment
since, by positioning multiple interferometers along the baseline, it will be
possible to probe the long and short $\Delta z$ configurations simultaneously
and thereby maximise the sensitivity across all frequencies.

Although $\mathcal{N}(\mathcal{N}-1)/2$ unique gradiometer measurements can be
performed on $\mathcal{N}$ interferometers, at most $\mathcal{N}-1$
gradiometers can be combined in the multi-gradiometer likelihood-based
analysis presented in section 4. This follows from the non-singularity
condition of the data’s covariance matrix $\boldsymbol{\mathrm{\Sigma}}_{k}$
and can be heuristically understood as the consequence of requiring that
gradiometer measurements in a given set are linearly independent, i.e. that no
gradiometer measurement be written in terms of the sum of other gradiometer
measurements. For example, for three AIs, three gradiometer measurements can
be performed, namely $\\{\Phi^{(1,2)},\Phi^{(1,3)},\Phi^{(2,3)}\\}$; since
$\Phi^{(i,j)}=\Phi^{(i,k)}-\Phi^{(j,k)}$ for
$i,j,k\in\\{1,\ldots\mathcal{N}\\}$, a multi-gradiometer analysis containing
measurement $\Phi^{(1,2)}$ cannot contain both $\Phi^{(2,3)}$ and
$\Phi^{(1,3)}$, but can contain $\Phi^{(1,2)}$ and either $\Phi^{(2,3)}$ or
$\Phi^{(1,3)}$, so two gradiometer measurements in total. Alternatively, this
statement can be recast in the language of graph theory. Mapping an AI to a
vertex and a gradiometer measurement to an edge, a set of gradiometer
measurements can be understood as a graph in which two vertices can be
connected by at most one edge and no cycles are present; the former follows
from the requirement that measurements must not be degenerate, while the
latter follows from the requirement that gradiometer measurements in a set are
linearly independent. Recalling that gradiometer measurements are insensitive
to the ordering of interferometer pairs (i.e. are undirected), the graphs
above are defined as trees (i.e. connected undirected graphs that contains no
cycles). By definition, it then follows that a tree on $\mathcal{N}$ vertices
contains $\mathcal{N}-1$ edges [61], so a set of gradiometer measurements on
$\mathcal{N}$ AIs contains $\mathcal{N}-1$ unique and linearly-independent
elements.

After fixing the interferometers’ positions along the baseline, the number of
sets of cardinality $\mathcal{N}-1$ containing unique and linearly-independent
gradiometer measurements can be efficiently calculated in graph theory.
Indeed, via Caylely’s formula, the number of trees on $\mathcal{N}$ labeled
vertices is $\mathcal{N}^{\mathcal{N}-2}$ [61]; hence, the number of sets each
containing $\mathcal{N}-1$ gradiometer measurements that can be combined in a
multi-gradiometer likelihood analysis is $\mathcal{N}^{\mathcal{N}-2}$.

Any choice of $\mathcal{N}-1$ gradiometer pairs out of
$\mathcal{N}^{\mathcal{N}-2}$ possibilities gives the same sensitivity to
ULDM. Indeed, for any choice of $\mathcal{N}-1$ gradiometers there exists a
map, which preserves the value of the test statistic, from the corresponding
data vector to that of any other choice of $\mathcal{N}-1$ interferometer
pair. For example, let $\mathcal{N}=3$; the set of possible gradiometer
measurements is given by $\\{\Phi^{(1,2)},\Phi^{(1,3)},\Phi^{(2,3)}\\}$.
Suppose that we choose $\\{\Phi^{(1,2)},\Phi^{(1,3)}\\}$ in our likelihood;
then, using Eq. (4.3), the corresponding data vector is
$\mathbf{d}_{k}=(R^{(1,2)},I^{(1,2)},R^{(1,3)},I^{(1,3)})^{\mathrm{T}}\propto(\Phi^{(1,2)},\Phi^{(1,2)},\Phi^{(1,3)},\Phi^{(1,3)})^{\mathrm{T}}$.
Now suppose that we want to perform the multi-gradiometer analysis using a
different set of measurements, namely $\\{\Phi^{(2,3)},\Phi^{(1,3)}\\}$, whose
corresponding data vector is
$\mathbf{d}^{\prime}_{k}=(R^{(2,3)},I^{(2,3)},R^{(1,3)},I^{(1,3)})^{\mathrm{T}}\propto(\Phi^{(2,3)},\Phi^{(2,3)},\Phi^{(1,3)},\Phi^{(1,3)})^{\mathrm{T}}$.
It can be shown that there exists a matrix $\mathbf{M}$ which maps
$\mathbf{d}_{k}$ to $\mathbf{d}^{\prime}_{k}$ and effectively corresponds to
the index relabeling $1\rightarrow 3,3\rightarrow 1$ and $2\rightarrow 2$. In
geometric terms, this can be understood as a relabelling of vertices in a
graph. Recalling that the data’s covariance matrix
$\mathbf{\Sigma}_{k}=\left\langle\mathbf{d}_{k}\mathbf{d}^{\mathrm{\mathbf{T}}}_{k}\right\rangle$
can be decomposed into signal and noise contributions as
$\mathbf{\Sigma}_{k}=\mathbf{S}_{k}+\mathbf{B}_{k}$, it follows that
$\mathbf{S}^{\prime}_{k}=\mathbf{M}\mathbf{S}_{k}\mathbf{M}^{T}$ and
$\mathbf{B}^{\prime}_{k}=\mathbf{M}\mathbf{B}_{k}\mathbf{M}^{T}$. Using the
cyclicity of the trace and Eq. (4.10), the test statistics for the two
different choices of gradiometer measurements are equal at all frequencies,
i.e. $q^{\prime}_{k}=q_{k}$ for all $k$ in the spectrum. Hence, an atom multi-
gradiometer experiment is insensitive to the choice of interferometer pairs.

Having identified the maximum number of gradiometer measurements that can be
combined in a multi-gradiometer analysis and after elucidating the equivalence
in the choice of gradiometer measurements, we are in a position to determine
which configurations maximise the experimental reach across the frequency band
of interest. The optimal choice of configuration, however, is noise model
dependent and so warrants a careful analysis. To gain intuition, we will
explore the effect of different interferometer layouts on the reach of multi-
gradiometer experiments in baselines of length
$L\sim\mathcal{O}($1\text{\,}\mathrm{k}\mathrm{m}$)$.

####  5.2.1 Example configurations

Figure 8: Schematic representation of different interferometer configurations
for different interferometer numbers: $\mathcal{N}=2$, $\mathcal{N}=3$ and
$\mathcal{N}=7$. With purple dots we show the positioning of interferometers
in ‘unequal spacing (ends)’ configurations. We show in red the positioning of
interferometers in ‘equal spacing’ configurations. In green we show the
‘unequal spacing (centre)’ configurations.

To facilitate the comparison between different interferometer configurations,
we will focus our analysis on three distinct multi-gradiometer layouts:
unequally-spaced intereferometers clustered at either end of the baseline;
equally-spaced configurations; and unequally-spaced configurations in which
each end of the baseline is equipped with an AI and the remaining
interferometers are clustered at the centre of the baseline. These three
layouts are shown in Fig. 8 for two, three and seven interferometers.

In detail, in the unequal spacing (ends) configuration, the interferometers
are clustered at either extremity of the baseline with the separation between
nearest neighbours in each cluster determined by $L/2(\mathcal{N}-1)$. The
separation between nearest neighbours in each cluster is shorter than the
shortest separation between interferometers in different clusters. When
$\mathcal{N}$ is odd, one additional interferometer is in the bottom cluster
relative to the cluster near the Earth’s surface. In the equal spacing
configuration, the separation between nearest neighbours is given by
$L/(\mathcal{N}-1)$. Finally, in the unequal spacing (centre) configuration,
each extremity of the baseline is equipped with an AI and the remaining
interferometers are clustered at the centre of the baseline. The separation
between nearest neighbours in the cluster is given by $L/2(\mathcal{N}-1)$,
which again implies that the separation between nearest neighbours in the
central cluster is shorter than the shortest separation between
interferometers in the cluster and at the ends of the baseline. In this case,
the interferometers are positioned symmetrically around the baseline’s
midpoint.

####  5.2.2 Searches below 1 Hz

Figure 9: Projected 95% CL exclusion curves for an atom multi-gradiometer
using the ‘advanced’ parameters with $\mathcal{N}=3$ (left panel) and
$\mathcal{N}=5$ (right panel) assuming GGN is modelled by the NHNM. The grey
and blue lines show the exclusion curve for a single atom gradiometer
($\mathcal{N}=2$) for ASN-only and in an ASN and GGN background, respectively.
The red dot-dashed, purple dashed and green solid lines show the atom multi-
gradiometer exclusion curves for the equally-spaced, unequally-spaced (ends)
and unequally-spaced (centre) configurations defined in Fig. 8, respectively.
The multi-gradiometer exclusion curves regain parts of parameter space between
$\sim 0.1$ Hz and $\sim 1$ Hz and have a similar sensitivity to the
$\mathcal{N}=2$ configuration at higher frequencies. The orange shaded region
is excluded by MICROSCOPE [23].

In Fig. 9 we show the 95% CL exclusion curves in the ULDM frequency (or mass)
versus ULDM-electron coupling plane for different multi-gradiometer
configurations, assuming the ‘advanced’ design and under the assumption of the
NHNM GGN scenario. We begin by discussing the left panel, which compares the
equal and unequal (ends) scenario for $\mathcal{N}=3$ with the single atom
gradiometer ($\mathcal{N}=2$) configuration.888The equal spacing and unequal
spacing (centre) configurations are identical for $\mathcal{N}=3$. We observe
an enhancement in the $\mathcal{N}=3$ exclusion curves at frequencies between
$\sim 0.05$ Hz and $\sim 1$ Hz. The enhancement has a mild dependence on the
position of the interferometers along the baseline since as explained in the
previous section, searches in the $\sim(0.05-0.2)$ Hz frequency range favour
the positioning of an interferometer pair towards the bottom of the baseline,
while above $\sim 0.2$ Hz equally-spaced configurations are preferred.

Below $\sim 0.05$ Hz, we see that the sensitivity is largely insensitive to
the location of the interferometers. This is because at these low frequencies
the GGN length scale is much larger than the baseline ($\lambda_{\rm{GGN}}\gg
L$), so the GGN gradiometer phase is largely insensitive to the location of
the interferometers. Unlike in Fig. 6, where the enhancement in the sub-Hz
regime for a single atom gradiometer ($\mathcal{N}=2$) came at the expense of
sensitivity loss above 1 Hz, for $\mathcal{N}=3$ we find an exclusion curve
above 1 Hz that is comparable to the $\mathcal{N}=2$ configuration. We will
postpone a more thorough discussion of the regime above 1 Hz to the next sub-
section.

The right panel of Fig. 9 compares the equal, unequal (ends) and unequal
(centre) scenarios for $\mathcal{N}=5$ with the single atom gradiometer
($\mathcal{N}=2$) configuration. We see a similar pattern as in the left
panel, with an enhancement at frequencies between $\sim 0.05$ Hz and $\sim 1$
Hz. Again, we find a mild enhancement below $\sim 0.2$ Hz for the
configuration that has the smallest $\Delta z$ towards the bottom of the
baseline, namely, the unequal (ends) configuration. Comparing the left and
right panels, we observe a weak dependence (of at most a factor of a few) on
the number of interferometers. This stems from the fact that $\mathcal{N}\gg
1$ AIs can simultaneously probe a larger array of length scales, which in turn
means that the configurations will asymptotically approach the optimal
gradiometer length and AI vertical position to maximise the reach across all
frequencies. In this frequency range, the advantage of using $\mathcal{N}>3$
instead of $\mathcal{N}=3$ interferometers is minimal.

Therefore, we have shown that since multi-gradiometer configurations can
simultaneously probe different length scales and vertical positions, they can
achieve a frequency-dependent and weakly $\mathcal{N}$-dependent sensitivity
enhancement over the single atom gradiometer ($\mathcal{N}=2$) configuration.

####  5.2.3 Searches above 1 Hz

At high frequencies, the background of all experiments considered in this work
is dominated by atom shot noise. In contrast to the GGN and ULDM signals, atom
shot noise is position-independent and is inversely proportional to the number
of atoms $N_{\rm{atom}}$ employed in each interferometer sequence. As Fig. 9
demonstrates, the sensitivity above 1 Hz is only weakly-dependent on the
position of the interferometers along the baseline. Therefore, in the regime
where GGN can be neglected, we can estimate the sensitivity scaling with the
number of interferometers by making the unphysical assumption that all
interferometers are placed at the extreme ends of the baseline. In this case,
we expect a $\sqrt{\mathcal{N}}$ sensitivity scaling since increasing
$\mathcal{N}$ is equivalent to increasing $N_{\rm{atom}}$ in each cloud, and
we know that the sensitivity scales like $\sqrt{N_{\rm{atom}}}$ in the absence
of GGN (cf. section 4.1).

Considering a more realistic multi-gradiometer configuration with equal
spacing between the atom interferometers, we find analytically that the test
statistic is related to the test statistic for a single gradiometer with
maximum gradiometer length and identical parameters by

$q^{\mathrm{multi}}_{\mathcal{N}}=\frac{1}{9}\left(\frac{\mathcal{N}}{2}\right)^{2}\left(\frac{\mathcal{N}+1}{\mathcal{N}-1}\right)^{2}q_{\mathcal{N}=2}\,.$
(5.1)

Given that the test statistic is proportional to the fourth-power of the ULDM
coupling (cf. section 4.1), we find explicitly the $\sqrt{\mathcal{N}}$
scaling in the limit of large $\mathcal{N}$. Setting $\mathcal{N}=3$ in Eq.
(5.1) shows that there is no sensitivity enhancement when an additional atom
interferometer is positioned at the mid-point along the baseline. This
behaviour is the reason why the unequal spacing (centre) configurations, which
feature an AI at either end of the baseline and equally spaced interferometers
clustered around the baseline’s mid-point, display a smaller sensitivity
increase than equally-spaced interferometers. Conversely, the unequal spacing
(ends) configuration provides the best sensitivity at high frequencies since
the interferometers are furthest from the mid-point.

####  5.2.4 Searches in the large $c_{H}$ regime

Figure 10:  Impact of changing the Rayleigh wave’s horizontal ground speed
$(c_{H})$ on the projected 95% CL exclusion curves of atom gradiometers. Upper
panels: Exclusion curves for a single atom gradiometer ($\mathcal{N}=2$) for
the ‘intermediate’ (left) and ‘advanced’ (right) designs. Solid and dotted
lines show $c_{H}=205\,\mathrm{m\,s}^{-1}$ and
$c_{H}=3232\,\mathrm{m\,s}^{-1}$, respectively. The NHNM and NLNM (upper-right
only) scenarios are shown in orange and blue, respectively. Lower panels:
Comparison of the exclusion curves for $\mathcal{N}=2$ (blue) and
$\mathcal{N}=5$ (red) assuming that the AIs are equally spaced under the NHNM
scenario. When $c_{H}=3232$ m/s, the $\mathcal{N}=5$ configuration provides no
increase in sensitivity relative to $\mathcal{N}=2$. In all panels the orange
shaded region is excluded by MICROSCOPE [23] and the grey lines show the
exclusive curves assuming an atom shot noise-only background.

In the projections shown until this point, we have utilised a set of
geological parameters in which the GGN length scale $\lambda_{\rm{GGN}}$,
defined in Eq. (3.16), is similar in magnitude to the baseline at frequencies
between $\sim 10^{-1}$ Hz and $\sim 1$ Hz. This was achieved by considering a
km-baseline experiment and a ground material in which the Rayleigh wave’s
horizontal speed was $c_{H}\sim 250~{}\mathrm{m\,s^{-1}}$. We now consider the
scenario where $\lambda_{\rm{GGN}}\gtrsim L$. This would occur, for example,
if the ground material consisted of limestone, for which $c_{H}\simeq
3232~{}\mathrm{m\,s^{-1}}$ and
$\lambda_{\rm{GGN}}\simeq$1.3\text{\,}\mathrm{k}\mathrm{m}$~{}\left(\frac{3232~{}\mathrm{m\,s^{-1}}}{c_{H}}\right)^{-1}\,\left(\frac{2.5~{}\mathrm{Hz}}{\omega_{a}}\right)$.

The upper panels of Fig. 10 compare the sensitivity reach between the two
values of the horizontal speed: $205~{}\mathrm{m\,s^{-1}}$, which we show with
solid curves; and $3232~{}\mathrm{m\,s^{-1}}$, which we show with dotted
curves. In the upper-left panel, we focus on a 100 m baseline operating with
$\mathcal{N}=2$ and the ‘intermediate’ parameters from Table 1. Here, we only
display the projected reach in the NHNM scenario, since the NLNM curve is
degenerate with the atom shot noise-limited curve. In the upper-right panel,
we consider the scenario for a 1 km baseline operating with $\mathcal{N}=2$
and the ‘advanced’ parameters from Table 1. Differently from the left panel,
we display the sensitivity curves in both the NHNM and NLNM scenarios, which
we plot in orange and blue respectively. Under both the NHNM and NLNM
scenarios, the $c_{H}=3232~{}\mathrm{m\,s^{-1}}$ sensitivity curves are
significantly lower than the $c_{H}=205~{}\mathrm{m\,s^{-1}}$ curves over all
frequencies in the sub-Hz regime for both baselines, which suggests that
geological materials with large $c_{H}$ values can naturally act as GGN
filters owing to the smaller value of $\lambda_{\mathrm{GGN}}$.999There is an
11% difference in the ground density between the low and high $c_{H}$
scenarios but since the GGN phase scales linearly with the ground density (Eq.
(B.20) and the discussion that follows), this cannot explain the large change
observed in Fig. 10.

The lower panels in Fig. 10 compare the $\mathcal{N}=2$ configuration shown in
blue with the $\mathcal{N}=5$ equal spacing multi-gradiometer configuration,
shown in red dot-dashed and red dashed for $c_{H}=205~{}\mathrm{m\,s^{-1}}$
and $c_{H}=3232~{}\mathrm{m\,s^{-1}}$ respectively. In both lower panels we
focus on the NHNM scenario. In the lower-left panel, where we study the
‘intermediate’ scenario with a 100 m baseline, we observe that when GGN
dominates the background (in the sub-Hz range), only a small sensitivity
enhancement is observed for $\mathcal{N}=5$ compared to $\mathcal{N}=2$ when
$c_{H}=205~{}\mathrm{m\,s^{-1}}$, and there is essentially no improvement when
$c_{H}=3232~{}\mathrm{m\,s^{-1}}$. This follows because the GGN length scale
of the Rayleigh wave’s profile underground is already $\mathcal{O}(100)$ m
when $c_{H}=205~{}\mathrm{m\,s^{-1}}$ and significantly exceeds the length of
the baseline when $c_{H}=3232~{}\mathrm{m\,s^{-1}}$, so the GGN gradiometer
phase is largely insensitive to the location of the interferometers. In the
km-baseline ‘advanced’ scenario, shown in the lower-right panel, a
$\mathcal{N}=5$ configuration significantly improves the sensitivity in the
$\sim 0.1$ Hz to $\sim 1$ Hz frequency range relative to $\mathcal{N}=2$ when
$c_{H}=205~{}\mathrm{m\,s^{-1}}$. However, we again see no such improvement
for $c_{H}=3232~{}\mathrm{m\,s^{-1}}$, since at these frequencies
$\lambda_{\mathrm{GGN}}\gg L$. This highlights the fact that in the
$\lambda_{\mathrm{GGN}}\gg L$ regime, the multi-gradiometer configuration is a
less useful tool to mitigate GGN.

##  6 Discussion and summary

Single-photon atom gradiometry is a powerful experimental technique that can
be employed to search for the oscillation of atomic transition energies
induced by ULDM. By developing a robust statistical analysis based on a
likelihood-based frequentist framework, we have provided a careful treatment
that characterises the impact of GGN on ULDM searches with vertical atom
gradiometers. We have modelled the GGN as arising from an incoherent
superposition of fundamental Rayleigh-wave modes, and used Peterson’s data-
driven models to capture the expected behaviour of low and high seismic-noise
scenarios. We find that GGN has the potential to significantly impact the
sensitivity reach of longer-baseline experiments,
$L\gtrsim\mathcal{O}(100\,~{}\mathrm{m})$, in the sub-Hz regime. However, we
have shown that vertical atom multi-gradiometer experiments, consisting of
three or more atom interferometers in the same baseline, can recover
significant parts of parameter space. Our results demonstrate how robust
sensitivity projections for vertical atom gradiometers, such as future
versions of AION and MAGIS-100, can be obtained for frequencies down to
$10^{-3}$ Hz, corresponding to a ULDM mass of $\sim 10^{-17}$ eV.

There is scope to go beyond the analysis presented in this work. Firstly, the
geological parameters that we used here have been chosen to be a generic
representation, so may not accurately describe the geology of the site where
AION-100, MAGIS-100 or other future vertical atom gradiometers will be
located. Secondly, the homogeneous half space model for the ground surrounding
the experiment may not accurately describe more realistic geological settings.
Indeed, the model does not describe the existence of geological anisotropies
or the presence of different geological strata. For example, it may be the
case that isotropy is broken by a body of water located close to the shaft, as
would occur at the Boulby Underground Laboratory or at CERN owing to their
vicinity to the North Sea and Lake Geneva, respectively. Additionally,
different geological strata would give rise to a much richer spectrum of
Rayleigh modes. This is because the additional strata would create more
interfaces for the reflection and propagation of P- and S-waves, and this
would generate Rayleigh wave overtones beyond the fundamental mode. Thirdly,
in addition to the geological model, our calculations are dependent on the
vertical displacement spectrum in the sub-Hz regime. A specific site may not
follow either of Peterson’s NLNM or NHNM scenarios, which attempts to bracket
the range of possibilities. Therefore, a site-specific spectrum acquired as
part of a dedicated site characterisation campaign would be preferred.
Finally, our model of GGN could be extended to include the contributions
beyond seismic Rayleigh waves, such as atmospheric pressure-induced
perturbations and seismic body waves that propagate through the Earth. In
summary, to further improve the sensitivity predictions in the sub-Hz regime
for future versions of AION or MAGIS-100, the homogeneous half space model
should be extended to account for the possibilities discussed above, and a
detailed site-specific seismic and ground measurement campaign should be
undertaken. We emphasise that a key advantage of the likelihood-based analysis
presented in this work is that the updated geological models and the site-
specific information can be straightforwardly incorporated into the analysis
procedure.

While our focus has been on scalar ULDM searches, we highlight that the
likelihood-based frequentist framework developed here can be extended to
stochastic gravitational wave searches with an atom multi-gradiometer
experiment. A study exploring these signals in the mid-frequency band using
this framework will be the focus of future work.

## Acknowledgments

We are grateful to members of the AION Collaboration for many fruitful
discussions and to John Ellis for comments on the manuscript. LB is a
recipient of an STFC quota studentship. CM acknowledges support from the
Science and Technology Facilities Council (STFC) Grants ST/T00679X/1,
ST/N004663/1. JM acknowledges support from the University of Cambridge Isaac
Newton Trust. We are grateful to the STFC for funding contributions towards
the AION-10 and MAGIS-100 experiments. For the purpose of open access, the
authors have applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising from this submission. No
experimental datasets were generated by this research.

##  Appendix A ULDM signal phase

A linearly-coupled scalar ULDM field $\phi(t)$, defined in Eq. (3.6), which
has interactions governed by the Lagrangian

$\mathcal{L}_{\phi}\supset\phi(t)\sqrt{4\pi
G_{N}}\left[\frac{d_{e}}{4e^{2}}F_{\mu\nu}F^{\mu\nu}-d_{m_{e}}m_{e}\overline{\psi}_{e}\psi_{e}\right],$
(A.1)

will induce oscillations in an atom’s transition frequency (see e.g., [42, 43,
44]):

$\omega_{A}(t)=\omega_{A}+\frac{1}{\sqrt{2}}\,\overline{\Delta\omega_{A}}\,\sum_{a}\alpha_{a}\sqrt{F_{\mathrm{DM}}(v_{a})}\cos(\omega_{a}t+\theta_{a})\;.$
(A.2)

Here, $d_{e}$ and $d_{m_{e}}$ parameterise the coupling strength relative to
Newton’s gravitational constant $G_{N}$, $\alpha_{a}$ is a Rayleigh
distributed variable and $\theta_{a}$ is a random phase. The sum is carried
over DM velocity classes of width $\Delta v$ that are weighted by the DM speed
distribution. For the $5\mathrm{s}^{2}\,{}^{1}\mathrm{S}_{0}\leftrightarrow
5\mathrm{s}5\mathrm{p}\,^{3}\mathrm{P}_{1}$ clock transition in
${}^{87}\mathrm{Sr}$, which we assume throughout this work,
$\omega_{A}=2.697\times 10^{15}~{}\mathrm{rad}/\mathrm{s}$, while

$\frac{\overline{\Delta\omega_{A}}}{\omega_{A}}=d_{\phi}\sqrt{4\pi
G_{N}}\frac{\sqrt{2\rho_{\mathrm{DM}}}}{m_{\phi}}\,,$ (A.3)

where the ULDM-Standard Model coupling strength is

$d_{\phi}=d_{m_{e}}+(2+\xi_{A})d_{e}\;,$ (A.4)

and $\xi_{A}$ is a calculable parameter that takes the value $\xi_{A}\approx
0.06$ for ${}^{87}\mathrm{Sr}$ [62].

Following the procedure described in Ref. [13], the expression for the
gradiometer phase in natural units measured at the end of the $m$th launch
between the $i$th and $j$th AIs, which are separated by the vertical distance
$\Delta z^{(i,j)}$, can be expressed as

$\Phi^{(i,j)}_{\mathrm{DM},m}=\frac{\Delta
z^{(i,j)}}{L}\sum_{a}\alpha_{a}\sqrt{F_{\mathrm{DM}}(v_{a})}\,A_{a}\cos\phi_{a,m}\;,$
(A.5)

where we have defined

$\displaystyle
A_{a}=\frac{8}{\sqrt{2}}\,\frac{\overline{\Delta\omega_{A}}}{\omega_{a}}\sin\left[\frac{\omega_{a}nL}{2}\right]\sin\left[\frac{\omega_{a}(T-(n-1)L)}{2}\right]\sin\left[\frac{\omega_{a}T}{2}\right]\;,$
(A.6) $\displaystyle\phi_{a,m}=\omega_{a}\left(\frac{2T+L}{2}+m\Delta
t\right)+\theta_{a}\;,$ (A.7)

and, as previously stated, $L$ is the baseline, $n$ is the number of LMT
kicks, while $T$ is the interrogation time. For definiteness, the random
variables that enter Eq. (A.5) are $\alpha_{a}$ and $\phi_{a,m}$.

###  A.1 Statistical Properties

The statistical properties of the ULDM signal follow straightforwardly from
the statistical properties of the ULDM wave. The expectation value of the ULDM
phase shift is

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Big{\rangle}=\frac{\Delta
z^{(i,j)}}{L}\sum_{a}\Big{\langle}\alpha_{a}\sqrt{F_{\mathrm{DM}}(v_{a})}\,A_{a}\cos\phi_{a,m}\Big{\rangle}\,.$
(A.8)

Since the random phase $\theta_{a}$ and the Rayleigh distributed random
variable $\alpha_{a}$ are independent, the expectation value will act on the
cosine term and separately on the amplitude, such that

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Big{\rangle}=\frac{\Delta
z^{(i,j)}}{L}\sum_{a}\langle\alpha_{a}\rangle\sqrt{F_{\mathrm{DM}}(v_{a})}\,A_{a}\langle\cos\phi_{a,m}\rangle\,.$
(A.9)

According to the definition of the Rayleigh distribution,
$\langle\alpha_{a}\rangle>0$; however, since $\phi_{a,m}$ depends linearly on
the random phase $\theta_{a}$, $\langle\cos\phi_{a,m}\rangle$ vanishes, such
that

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Big{\rangle}=0\,.$
(A.10)

The covariance of the gradiometer also follows from the statistical properties
of the ULDM wave. Neglecting for brevity the factors depending on the
gradiometer length, using the vanishing expectation value of the ULDM phase
and recalling the fact that $\alpha_{a}$ and $\theta_{a}$ are independent
random variables, the expected covariance is proportional to

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{DM},m^{\prime}}\Big{\rangle}$
$\displaystyle\propto\sum_{a,a^{\prime}}\langle\alpha_{a}\alpha_{a^{\prime}}\rangle\sqrt{F_{\mathrm{DM}}(v_{a})F_{\mathrm{DM}}(v_{a^{\prime}})}\,A_{a}A_{a^{\prime}}\langle\cos\phi_{a,m}\cos\phi_{a^{\prime},m^{\prime}}\rangle\,.$
(A.11)

Since $\phi_{a,m}$ depends on a random phase $\theta_{a}$, the sum over the
velocity classes (labelled by $a$ and $a^{\prime}$) is non-zero only when the
random phases are the same, i.e. $a^{\prime}=a$. Hence, we may write

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{DM},m^{\prime}}\Big{\rangle}$
$\displaystyle\propto\sum_{a}\langle\alpha_{a}^{2}\rangle
F_{\mathrm{DM}}(v_{a})\,A_{a}^{2}\langle\cos\phi_{a,m}\cos\phi_{a,m^{\prime}}\rangle\,.$
(A.12)

To further simplify the expression above, we observe that
$\left\langle\cos\phi_{a,m}\cos\phi_{a,m^{\prime}}\right\rangle$ can be re-
expressed in terms of double angle formulas as

$\begin{split}\left\langle\cos\phi_{a,m}\cos\phi_{a,m^{\prime}}\right\rangle&=\frac{1}{2}\left[\left\langle\cos\left(\phi_{a,m}-\phi_{a,m^{\prime}}\right)+\cos\left(\phi_{a,m}+\phi_{a,m^{\prime}}\right)\right\rangle\right]\\\
&=\frac{1}{2}\cos\left(\omega_{a}\Delta
t(m-m^{\prime})\right)+\frac{1}{2}\left\langle\cos\left(\phi_{a,m}+\phi_{a,m^{\prime}}\right)\right\rangle\\\
&=\frac{1}{2}\cos\left(\omega_{a}\Delta t(m-m^{\prime})\right)\,,\end{split}$
(A.13)

where the second line follows from the definition of $\phi_{a,m}$ and the last
line from the vanishing expectation value of the cosine of a random phase.
Using this result and $\langle\alpha_{a}^{2}\rangle$ = 2, and also restoring
the dependence on the gradiometer length, we find that the covariance takes
the form

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{DM},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{DM},m^{\prime}}\Big{\rangle}$
$\displaystyle=\frac{\Delta r^{(i,j)}}{L}\frac{\Delta
r^{(i^{\prime},j^{\prime})}}{L}\sum_{a}F_{\mathrm{DM}}(v_{a})\,A_{a}^{2}\cos(\omega_{a}\Delta
t(m-m^{\prime}))\,.$ (A.14)

##  Appendix B GGN from the fundamental Rayleigh mode

Our GGN model assumes that the dominant contribution arises from the
fundamental Rayleigh mode. The fundamental Rayleigh mode consists of S-waves
with vertical displacements that are coupled to P-waves. In this appendix, we
begin by deriving an expression for the angle-averaged gravitational potential
caused by the fundamental Rayleigh mode. This potential is used as an input to
calculate the phase shift induced in an atom interferometer. We then consider
the statistical properties of the phase shift. Our notation largely follows
Refs. [25, 63], and we also make extensive use of the discussion in Ref. [24].

We begin with some definitions. The propagation speed of P-waves, $c_{P}$, is
determined by the material’s density $\rho_{0}$, and bulk and shear moduli $K$
and $\mu$, respectively, as

$c_{P}=\sqrt{\frac{K+4\mu/3}{\rho_{0}}}\,.$ (B.1)

The propagation speed of S-waves, $c_{S}$, is related to that of P-waves by

$c_{S}=\sqrt{\frac{1-2\nu}{2-2\nu}}\,c_{P}\;,$ (B.2)

where $\nu$ is the material’s Poisson ratio, which depends on the ground’s
bulk and shear moduli as

$\nu=\frac{3K-2\mu}{2(3K+\mu)}\,.$ (B.3)

We model the geological surroundings of the experimental site as a homogeneous
half space with constant ground density $\rho_{0}$ and Poisson ratio $\nu$
across all depths. This is an idealised scenario which implies that the
Rayleigh waves are isotropic and that only the fundamental Raleigh mode is
present: Rayleigh overtones require geologically stratified structures, which
are not present in our model. The fundamental Raleigh mode has frequency
$f_{a}$, angular frequency $\omega_{a}=2\pi f_{a}$, horizontal wave number
$k_{a}$, horizontal propagation speed $c_{H}=\omega_{a}/k_{a}$, and horizontal
propagation direction $\hat{k}$. The horizontal propagation speed $c_{H}$ is
related to $c_{S}$ by $c_{H}=c_{S}\sqrt{\zeta}$, where $\zeta$ is the real
root of the cubic equation

$\zeta^{3}-8\zeta^{2}+8\left(\frac{2-\nu}{1-\nu}\right)\zeta-\frac{8}{(1-\nu)}=0\;.$
(B.4)

It is convenient to introduce the dimensionless ratios $s$ and $q$,

$\displaystyle s$ $\displaystyle=\sqrt{1-\left(c_{H}/c_{S}\right)^{2}}\;,$
(B.5) $\displaystyle q$
$\displaystyle=\sqrt{1-\left(c_{H}/c_{P}\right)^{2}}\;,$ (B.6)

which will appear as the vertical e-folding rate to horizontal wave number for
S-waves and as the vertical e-folding rate to horizontal wave number for
P-waves, respectively. Equation (B.4) encodes the relation

$1+s^{2}=2\sqrt{qs}\;,$ (B.7)

which we make extensive use of in the manipulations below.

Finally, the definition of our coordinate system is shown in Fig. 2. It is
convenient to express all vectors in cylindrical coordinates,

$\displaystyle\vec{x}$
$\displaystyle=(\varrho\cos\theta,\varrho\sin\theta,0)\,,$ (B.8)
$\displaystyle\vec{r}$
$\displaystyle=(\varrho\cos\theta,\varrho\sin\theta,z)\;,$ (B.9)
$\displaystyle\vec{k}$
$\displaystyle=(k\cos\theta^{\prime},k\sin\theta^{\prime},0)\,.$ (B.10)

In this system, $\vec{x}$ defines a point on the horizontal plane and $z$
defines the depth into the ground. The product of the horizontal wave number
and horizontal distance, $\vec{k}\cdot\vec{x}$, may be written as

$\vec{k}\cdot\vec{x}=k\varrho(\cos\theta^{\prime}\cos\theta+\sin\theta^{\prime}\sin\theta)=k\varrho\cos(\theta-\theta^{\prime})\;.$
(B.11)

Having defined our terms, we now proceed to calculate the surface and
underground density perturbations induced by the Rayleigh waves. Adapting the
expressions from Ref. [24] to our coordinate system, the Rayleigh wave
displacement vector is given by

$\vec{\xi}_{a,R}(\varrho,\theta,z,t)=\left(\xi_{a,H}(z)\hat{k}-\xi_{a}(z)\vec{e}_{z}\right)\exp\left[{i(k_{a}\varrho\cos(\theta-\theta^{\prime})-\omega_{a}t+\widetilde{\theta}_{a})}\right]\,,$
(B.12)

where $\widetilde{\theta}_{a}$ is a random phase, $\vec{e}_{z}$ is the unit
vector pointing downward and

$\displaystyle\xi_{a}(z)$
$\displaystyle=\xi_{a}\left(\frac{1+s^{2}}{1-s^{2}}\right)\left[\exp{(-q\omega_{a}z/c_{H})}-\frac{2}{1+s^{2}}\exp{(-s\omega_{a}z/c_{H})}\right]\,,$
(B.13) $\displaystyle\xi_{a,H}(z)$
$\displaystyle=\frac{i\xi_{a}}{q}\left(\frac{1+s^{2}}{1-s^{2}}\right)\left[\exp{(-q\omega_{a}z/c_{H})}-\frac{2qs}{1+s^{2}}\exp{(-s\omega_{a}z/c_{H})}\right]\;,$
(B.14)

where $\xi_{a}$ is the vertical displacement at the surface and the factor $i$
in $\xi_{a,H}(z)$ ensures that particle motion is retrograde at the surface.
Thus, we see explicitly that the amplitude is characterised by an
exponentially decaying profile while moving away from the surface.

The displacement field induced by the fundamental Rayleigh wave produces a
fractional density perturbation at depth $z$ [25],

$\frac{\delta\rho_{a}(z>0)}{\rho_{0}}=\left[\xi_{a}\delta(z)+\mathcal{R}_{a}(z)\right]\exp\left[{i(k_{a}\varrho\cos(\theta-\theta^{\prime})-\omega_{a}t+\widetilde{\theta}_{a})}\right]\,,$
(B.15)

where the first term is the surface contribution and arises from the
displacement normal to the surface ($\delta(z)$ is the Dirac delta function),
while the second term is the contribution from inside the medium and follows
from the ground-mass continuity equation, such that

$\displaystyle\mathcal{R}_{a}(z)$
$\displaystyle=-i\frac{\omega_{a}}{c_{H}}\xi_{a,H}(z)+\frac{\partial\xi_{a}(z)}{\partial
z}$ (B.16)
$\displaystyle=\frac{\omega_{a}}{c_{H}}\xi_{a}\frac{(q^{2}-1)}{q}\left(\frac{1+s^{2}}{1-s^{2}}\right)\exp{(-q\omega_{a}z/c_{H})}\;.$
(B.17)

Following Ref. [24], the angle-averaged perturbing gravitational potential
induced by a fundamental Rayleigh mode of frequency $\omega_{a}$, and
experienced by a test mass underground at time $t$ and position
$\vec{r}_{0}=(0,0,z_{0})$, is

$V_{a}\left(z_{0},t\right)=-2\pi
G\rho_{0}\,\cos(\omega_{a}t+\widetilde{\theta}_{a})\,\int\varrho\mathrm{d}\varrho\,\mathrm{d}z\,\frac{1}{\sqrt{\varrho^{2}+(z-z_{0})^{2}}}\left[\xi_{a}\delta(z)+\mathcal{R}_{a}(z)\right]J_{0}(\omega_{a}\varrho/c_{H})\,,$
(B.18)

where $J_{0}(\omega_{a}\varrho/c_{H})$ is the Bessel function of first kind.
Performing the integrals, and summing over all frequencies, we arrive at the
expression

$\displaystyle V_{a}\left(z_{0},t\right)=$ $\displaystyle-2\pi
G\rho_{0}\,\sum_{a}\xi_{a}\,\frac{c_{H}}{q\,\omega_{a}}\,\left(\frac{1+s^{2}}{1-s^{2}}\right)$
(B.19)
$\displaystyle\times\left[\left(1+\sqrt{q/s}\right)\exp{(-\omega_{a}z_{0}/c_{H})}-2\exp{(-q\omega_{a}z_{0}/c_{H})}\right]\cos(\omega_{a}t+\widetilde{\theta}_{a})\,.$

###  B.1 Rayleigh wave signal phase

The leading order phase shift induced by the perturbing potential on the atom
interferometer can be calculated using the perturbative approach [64], in
which the GGN potential is treated as a perturbation to Minkowski space. In
flat space time, the atom trajectories follow straight lines and the
interferometer can be geometrically thought of as a kite in spacetime. In this
geometry, the separation, laser and propagation phases are all zero. Using the
standard approach for perturbative atom interferometer calculations [64], the
leading order phase shift from the fundamental Rayleigh wave field with
angular frequency $\omega_{a}$ leaves a non-zero contribution to the
propagation phase. In all experimental configurations considered here, we take
the launch velocity $v_{\mathrm{in}}$ and velocity kick from the $n$ LMT
pulses $v_{\mathrm{kick}}$ to be much smaller than the horizontal speed of the
Rayleigh wave, such that $v_{\mathrm{kick}}/c_{H}$, $v_{\mathrm{in}}/c_{H}\ll
1$. Additionally, for all frequencies considered in this paper, the Rayleigh
wave’s wavelength is significantly larger than the length scales probed by the
atoms during interferometry, i.e. $\omega_{a}Tv_{\mathrm{kick}}/c_{H}\ll 1$,
such that the leading order GGN-induced phase shift measured by the $i$th AI
at the end of the $m$th launch takes the form

$\Phi^{(i)}_{\mathrm{GGN},m}=\sum_{a}\xi_{a}\left(\widetilde{A}_{a}\exp\left(-q\frac{\omega
z_{i}}{c_{H}}\right)+\widetilde{B}_{a}\exp\left(-\frac{\omega
z_{i}}{c_{H}}\right)\right)\cos{\widetilde{\phi}_{a,m}}\;.$ (B.20)

In this expression we have defined

$\displaystyle\widetilde{A}_{a}=nk_{A}\,\frac{16\pi
G\rho_{0}}{\omega_{a}^{2}}\left(\frac{1+s^{2}}{1-s^{2}}\right)\sin^{2}\left(\frac{\omega_{a}T}{2}\right)\,,$
(B.21) $\displaystyle\widetilde{B}_{a}=-nk_{A}\,\frac{8\pi
G\rho_{0}}{q\,\omega_{a}^{2}}\left(\frac{1+s^{2}}{1-s^{2}}\right)\,\left(1+\sqrt{q/s}\right)\sin^{2}\left(\frac{\omega_{a}T}{2}\right)\,,$
(B.22) $\displaystyle\widetilde{\phi}_{a,m}=\omega_{a}\left(T+m\,\Delta
t\right)+\widetilde{\theta}_{a}\,,$ (B.23)

and assumed that the momentum kick from the LMT sequence is $nk_{A}$. For
definiteness, the independent random variables that enter Eq. (B.20) are
$\xi_{a}$ and $\widetilde{\theta}_{a}$.

With our assumption that the ground composition is homogeneous down the
baseline, it follows that the GGN gradiometer phase between the $i$th and
$j$th interferometers is

$\displaystyle\Phi^{(i,j)}_{\mathrm{GGN},m}=\sum_{a}\xi_{a}$
$\displaystyle\Bigg{\\{}\widetilde{A}_{a}\left[\exp\left(-q\frac{\omega
z_{i}}{c_{H}}\right)-\exp\left(-q\frac{\omega z_{j}}{c_{H}}\right)\right]$
(B.24) $\displaystyle+\widetilde{B}_{a}\left[\exp\left(-\frac{\omega
z_{i}}{c_{H}}\right)-\exp\left(-\frac{\omega
z_{j}}{c_{H}}\right)\right]\Bigg{\\}}\cos{\widetilde{\phi}_{a,m}}\;.$

In Figs. 4, 5, 6 and 9 we use the following ground parameters: $\rho_{0}=1800$
kg/m3, $\nu=0.33$, $c_{P}=440$ m/s and $c_{S}=220$ m/s. From these, we find
$s=0.36$, $q=0.88$ and $c_{H}=205$ m/s. In Fig. 10 we also show results with
$\rho_{0}=2000$ kg/m3, $\nu=0.34$, $c_{P}=6964$ m/s and $c_{S}=3464$ m/s,
which lead to $s=0.36$, $q=0.89$ and $c_{H}=3232$ m/s. We stress that these
parameters are representative estimates and that site-specific measurements
would be necessary to deduce these parameters.

###  B.2 Statistical Properties

Since the GGN phase is sourced by a stochastic superposition of waves, it
follows that the statistical properties of the ULDM signal also characterise
the GGN background contribution. Indeed, since the vertical displacement
$\xi_{a}$ is Rayleigh distributed at each angular frequency $\omega_{a}$ and,
importantly, is independent of the random phase $\theta_{a}$, the expectation
value of the GGN gradiometer phase takes the form

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{GGN},m}\Big{\rangle}$
$\displaystyle=\sum_{a}\langle\xi_{a}\rangle\Bigg{\\{}\widetilde{A}_{a}\left[\exp\left(-q\frac{\omega
z_{i}}{c_{H}}\right)-\exp\left(-q\frac{\omega z_{j}}{c_{H}}\right)\right]$
(B.25) $\displaystyle\qquad+\widetilde{B}_{a}\left[\exp\left(-\frac{\omega
z_{i}}{c_{H}}\right)-\exp\left(-\frac{\omega
z_{j}}{c_{H}}\right)\right]\Bigg{\\}}\,\langle\cos{\widetilde{\phi}_{a,m}}\rangle$
$\displaystyle=0\,,$

where the last equality follows from the fact that $\widetilde{\phi}_{a,m}$
depends linearly on a random phase $\widetilde{\theta}_{a}$.

The calculation of the covariance also follows identically from the ULDM case.
The only structural difference arises from the product of the exponential
functions in Eq. (B.24). Hence, using Eqs. (A.12)-(A.13), it can easily be
shown that the covariance is then given by

$\displaystyle\Big{\langle}\Phi^{(i,j)}_{\mathrm{GGN},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{GGN},m^{\prime}}\Big{\rangle}$
$\displaystyle=\sum_{a}\frac{\langle\xi_{a}^{2}\rangle}{2}F_{\mathrm{GGN},a}^{(i,j)}F_{\mathrm{GGN},a}^{(i^{\prime},j^{\prime})}\cos(\omega_{a}\Delta
t(m-m^{\prime}))\,,$ (B.26)

where

$\displaystyle
F_{\mathrm{GGN},a}^{(i,j)}=F_{\mathrm{GGN},a}^{(i)}-F_{\mathrm{GGN},a}^{(j)}\,,$
(B.27) $\displaystyle
F_{\mathrm{GGN},a}^{(i)}=\widetilde{A}_{a}\exp\left(-q\frac{\omega_{a}z_{i}}{c_{H}}\right)+\widetilde{B}_{a}\exp\left(-\frac{\omega_{a}z_{i}}{c_{H}}\right)\,.$
(B.28)

It should be noted that, differently from $\alpha_{a}$,
$\langle\xi_{a}^{2}\rangle\neq 2$; hence, we obtain a term
$\langle\xi_{a}^{2}\rangle/2$ that is absent from Eq. (A.14). This is because
$\alpha_{a}$ is a dimensionless random variable described by Eq. (3.8), such
that $\langle\alpha_{a}^{2}\rangle=2$, while $\xi_{a}$ is not.

##  Appendix C Covariance matrices

The covariance matrix of a multi-dimensional data vector $\mathbf{d}_{k}$ with
zero expectation value is simply expressed as
$\mathbf{\Sigma}_{k}=\langle\mathbf{d}_{k}\mathbf{d}_{k}^{\mathrm{\mathbf{T}}}\rangle$.
We remind the reader that the $\mathcal{N}(\mathcal{N}-1)$-dimensional data
vector is defined as

$\mathbf{d}_{k}=\Big{[}R_{k}^{(1,2)},I_{k}^{(1,2)},\ldots,R_{k}^{(n,n-1)},I_{k}^{(n,n-1)}\Big{]}^{\mathrm{\mathbf{T}}}$
(C.1)

in terms of the real and imaginary parts of the discrete Fourier transform
(DFT), which we can explicitly write out in terms of noise and signal phase
shifts measured by AG-($i,j$) as

$\displaystyle R_{k}^{(i,j)}$ $\displaystyle=\frac{\Delta
t}{\sqrt{T_{\mathrm{int}}}}\sum_{m=0}^{N-1}\left(\Phi^{(i,j)}_{\mathrm{Signal},m}+\Phi^{(i,j)}_{\mathrm{Noise},m}\right)c_{m,k}=R_{k,\mathrm{Signal}}^{(i,j)}+R_{k,\mathrm{Noise}}^{(i,j)}\;,$
(C.2) $\displaystyle I_{k}^{(i,j)}$ $\displaystyle=-\frac{\Delta
t}{\sqrt{T_{\mathrm{int}}}}\sum_{m=0}^{N-1}\left(\Phi^{(i,j)}_{\mathrm{Signal},m}+\Phi^{(i,j)}_{\mathrm{Noise},m}\right)s_{m,k}=I_{k,\mathrm{Signal}}^{(i,j)}+I_{k,\mathrm{Noise}}^{(i,j)}\;.$
(C.3)

For brevity, in Eqs. (C.2)-(C.3) we use the shorthand notation for
trigonometric functions, namely

$c_{m,k}=\cos\left(\frac{2\pi km}{N}\right)\;\text{ and
}\;s_{m,k}=\sin\left(\frac{2\pi km}{N}\right)\,.$ (C.4)

Equipped with Eq. (C.1) the covariance matrix of the $k$-frequency component
collected by all gradiometers can then be organized into a symmetric
$(\mathcal{N}(\mathcal{N}-1)\times\mathcal{N}(\mathcal{N}-1))$-dimensional
matrix taking the form

$\boldsymbol{\Sigma}_{k}=\left[\begin{array}[]{cccc}\left\langle
R^{(1,2)}_{k}R^{(1,2)}_{k}\right\rangle&\left\langle
R^{(1,2)}_{k}I^{(1,2)}_{k}\right\rangle&\ldots&\left\langle
R^{(1,2)}_{k}I^{(m-1,m)}_{k}\right\rangle\\\ \left\langle
I^{(1,2)}_{k}R^{(1,2)}_{k}\right\rangle&\left\langle
I^{(1,2)}_{k}I^{(1,2)}_{k}\right\rangle&\ldots&\left\langle
I^{(1,2)}_{k}I^{(m-1,m)}_{k}\right\rangle\\\ \vdots&&\ddots&\vdots\\\
\left\langle I^{(m-1,m)}_{k}R^{(1,2)}_{k}\right\rangle&\left\langle
I^{(m-1,m)}_{k}I^{(1,2)}_{k}\right\rangle&\ldots&\left\langle
R^{(m-1,m)}_{k}I^{(m-1,m)}_{k}\right\rangle\end{array}\right]\,.$ (C.5)

From Eq. (C.5), it is clear that each entry in the covariance matrix consists
of one of three unique products of imaginary and/or real parts, namely

$\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle,\,\left\langle
R^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle,\,\left\langle
I^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle\,.$ (C.6)

Additionally, from Eqs.(C.2)-(C.3), each of these entries will contain terms
pairing phases originating from identical or different sources. For example,
using Eq. (C.2) and without loss of generality, we find that

$\displaystyle\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\left\langle
R^{(i,j)}_{k,\mathrm{Signal}}R^{(i^{\prime},j^{\prime})}_{k,\mathrm{Signal}}\right\rangle+\left\langle
R^{(i,j)}_{k,\mathrm{Noise}}R^{(i^{\prime},j^{\prime})}_{k,\mathrm{Noise}}\right\rangle$
(C.7) $\displaystyle+\left\langle
R^{(i,j)}_{k,\mathrm{Noise}}R^{(i^{\prime},j^{\prime})}_{k,\mathrm{Signal}}\right\rangle+\left\langle
R^{(i,j)}_{k,\mathrm{Signal}}R^{(i^{\prime},j^{\prime})}_{k,\mathrm{Noise}}\right\rangle\,,$

where the last two terms mix the signal and noise phases. Under our working
assumption of uncorrelated noise and signal contributions, these terms vanish.
Indeed, using the expectation value of the models considered in this work, we
find

$\left\langle
R^{(i,j)}_{k,\mathrm{Noise}}R^{(i^{\prime},j^{\prime})}_{k,\mathrm{Signal}}\right\rangle\propto\left\langle\Phi^{(i,j)}_{\mathrm{Noise},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{Signal},m^{\prime}}\right\rangle=\left\langle\Phi^{(i,j)}_{\mathrm{Noise},m}\right\rangle\left\langle\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{Signal},m^{\prime}}\right\rangle=0\,.$
(C.8)

The fact that the signal and each noise source are uncorrelated across the
frequency spectrum allows us to decompose the covariance matrix
$\boldsymbol{\mathrm{\Sigma}}_{k}$ into a signal $\boldsymbol{\mathrm{S}}_{k}$
and noise or background $\boldsymbol{\mathrm{B}}_{k}$ contribution. Assuming
that each noise contribution is uncorrelated, we can also decompose the noise
covariance matrix $\boldsymbol{\mathrm{B}}_{k}$ in terms of the GGN and atom
shot noise contributions. In the following subsections we will obtain
analytical expressions for these covariance matrices.

###  C.1 ULDM covariance matrix

The covariance matrix for the ULDM matrix can be calculated by recalling the
statistical properties of the ULDM gradiometer phase. Focusing on the
covariance between the real parts of the data’s DFT,

$\begin{split}\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle&=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\sum_{m,m^{\prime}=0}^{N-1}\left\langle\Phi^{(i,j)}_{\mathrm{DM},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{DM},m^{\prime}}\right\rangle
c_{m,k}\,c_{m^{\prime},k}\\\ &=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\frac{\Delta z^{(i,j)}}{L}\frac{\Delta
z^{(i^{\prime},j^{\prime})}}{L}\\\
&\qquad\times\sum_{m,m^{\prime}=0}^{N-1}\sum_{a}F_{\mathrm{DM}}(v_{a})A_{a}^{2}\cos\left(\omega_{a}\Delta
t(m-m^{\prime})\right)c_{m,k}c_{m^{\prime},k}\,,\end{split}$ (C.9)

where the second equality follows from the covariance of the ULDM gradiometer
phase, shown in Eq. (A.14). To make further progress, we recall that the ULDM
speed indices $a$ can be mapped to those labelling the discrete frequencies
probed by the experiment. Hence, we may rewrite the sum over $m$ and
$m^{\prime}$ as

$\begin{split}\sum_{m,m^{\prime}=0}^{N-1}\cos\left(\omega_{a}\Delta
t(m-m^{\prime})\right)c_{m,k}c_{m^{\prime},k}&=\sum_{m,m^{\prime}=0}^{N-1}(c_{m,k^{\prime}}c_{m^{\prime},k^{\prime}}+s_{m,k^{\prime}}s_{m^{\prime},k^{\prime}})c_{m,k}c_{m^{\prime},k}\,,\end{split}$
(C.10)

where we defined the integer $k^{\prime}$ as $k^{\prime}=N\omega_{a}\Delta
t/2\pi=T_{\mathrm{int}}\omega_{a}/2\pi$ and used the trigonometric expansion
of the cosine of two angles. Rearranging the RHS of the equation above,

$\begin{split}\sum_{m,m^{\prime}=0}^{N-1}(c_{m,k^{\prime}}c_{m^{\prime},k^{\prime}}+s_{m,k^{\prime}}s_{m^{\prime},k^{\prime}})c_{m,k}c_{m^{\prime},k}&=\sum_{m^{\prime}=0}^{N-1}c_{m^{\prime},k^{\prime}}c_{m^{\prime},k}\sum_{m=0}^{N-1}c_{m,k^{\prime}}c_{m,k}\\\
&+\sum_{m^{\prime}=0}^{N-1}s_{m^{\prime},k^{\prime}}c_{m^{\prime},k}\sum_{m=0}^{N-1}s_{m,k^{\prime}}c_{m,k}\\\
&=\left(\frac{N}{2}\right)^{2}\delta_{k^{\prime},k}\,,\end{split}$ (C.11)

where the second equality follows from the orthogonality relations of
trigonometric functions. Using this result and expressing the sum in terms of
frequency indices, Eq. (C.9) can be simplified as

$\begin{split}\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle&=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\frac{\Delta z^{(i,j)}}{L}\frac{\Delta
z^{(i^{\prime},j^{\prime})}}{L}\sum_{k^{\prime}}F_{\mathrm{DM}}(v_{k^{\prime}})A_{k^{\prime}}^{2}\left(\frac{N}{2}\right)^{2}\delta_{k^{\prime},k}\,,\\\
&=\frac{\pi}{2}\frac{\Delta z^{(i,j)}}{L}\frac{\Delta
z^{(i^{\prime},j^{\prime})}}{L}\frac{F_{\mathrm{DM}}(v_{k})}{\Delta\omega}A_{k}^{2}\,,\end{split}$
(C.12)

where in going to the second line we used the property of the Kronecker delta
and the definition of the experimental angular frequency resolution
$\Delta\omega=2\pi/T_{\mathrm{int}}$. Finally, writing the integral
$F_{\mathrm{DM}}(v_{k})$ in terms of frequencies, we find

$\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle=\frac{\pi}{2}\frac{\Delta
z^{(i,j)}}{L}\frac{\Delta
z^{(i^{\prime},j^{\prime})}}{L}A_{k}^{2}\frac{1}{\Delta\omega}\int_{\omega_{k}-\Delta\omega/2}^{\omega_{k}+\Delta\omega/2}d\omega\,\frac{f_{\mathrm{DM}}(v_{\omega})}{m_{\phi}\,v_{\omega}}\,,$
(C.13)

where we defined $v_{\omega}=\sqrt{2\omega/m_{\phi}-2}$. For completeness, we
define the PSD of the ULDM signal common to all gradiometer pairs in terms of
all relevant experimental and phenomenological parameters as

$\begin{split}\left\langle
S_{\mathrm{DM},k}\right\rangle&=64\pi\,\frac{\overline{\Delta\omega_{A}}^{2}}{\omega_{k}^{2}}\sin^{2}\left[\frac{\omega_{k}nL}{2}\right]\sin^{2}\left[\frac{\omega_{k}(T-(n-1)L)}{2}\right]\sin^{2}\left[\frac{\omega_{k}T}{2}\right]\,\\\
&\quad\times\frac{1}{\Delta\omega}\int_{\omega_{k}-\Delta\omega/2}^{\omega_{k}+\Delta\omega/2}d\omega\,\frac{f_{\mathrm{DM}}(v_{\omega})}{m_{\phi}\,v_{\omega}}\,,\end{split}$
(C.14)

where the amplitude of the ULDM-induced oscillation in the transition
frequency $\overline{\Delta\omega_{A}}$ is defined in Eq. (A.3). As a result
of these arguments, it follows that Eq. (C.12) can be expressed as

$\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle=\frac{1}{2}\frac{\Delta
z^{(i,j)}}{L}\frac{\Delta z^{(i^{\prime},j^{\prime})}}{L}\left\langle
S_{\mathrm{DM},k}\right\rangle\,.$ (C.15)

Having determined the form of the covariance matrix entry involving only the
real part of the data’s DFT, we proceed with computing the terms involving
only the imaginary components using a simple trick. As pointed out in Ref.
[58], the result for the terms involving only the imaginary components is
equivalent to that involving only the real part under the redefinition
$c_{m,k}\rightarrow s_{m,k}$ and $c_{m^{\prime},k}\rightarrow
s_{m^{\prime},k}$. Since Eq. (C.11) holds under this redefinition, we also
find that

$\left\langle
I^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle\equiv\frac{1}{2}\frac{\Delta
z^{(i,j)}}{L}\frac{\Delta z^{(i^{\prime},j^{\prime})}}{L}\left\langle
S_{\mathrm{DM},k}\right\rangle\,.$ (C.16)

Similarly, we can determine the terms in the covariace matrix involving the
product of the imaginary and real components of the data’s DFT by using the
mapping $c_{m,k}\rightarrow c_{m,k}$ and $c_{m^{\prime},k}\rightarrow
s_{m^{\prime},k}$. Differently from the other products, we now find a sum over
a product of two cosines and one sine. Since the sine function has odd parity
and the cosine is of even parity, the sum will vanish for all values of $k$:

$\left\langle R^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle\equiv
0\,.$ (C.17)

In summary, we find that the entries in the ULDM signal covariance matrix take
the form

$\displaystyle\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{1}{2}\frac{\Delta z^{(i,j)}}{L}\frac{\Delta
z^{(i^{\prime},j^{\prime})}}{L}\left\langle S_{\mathrm{DM},k}\right\rangle\,,$
(C.18) $\displaystyle\left\langle
I^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{1}{2}\frac{\Delta z^{(i,j)}}{L}\frac{\Delta
z^{(i^{\prime},j^{\prime})}}{L}\left\langle S_{\mathrm{DM},k}\right\rangle\,,$
(C.19) $\displaystyle\left\langle
R^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=0\,.$ (C.20)

####  C.1.1 Behaviour under long and short integration times

In light of Eq. (C.13), the non-vanishing entries of the ULDM signal
covariance matrix scale differently with the integration time depending on the
frequency resolution attained by the experimentalist. In the regime where the
integration time $T_{\mathrm{int}}$ is much larger than the coherence time of
the signal $\tau_{c}$, the experimentalist would be able to resolve the full
frequency width of the ULDM signal. If we take the limit
$T_{\mathrm{int}}\rightarrow\infty$, the width of the frequency bins tends to
zero, i.e. $\Delta\omega\rightarrow 0$, such that

$\lim_{\Delta\omega\rightarrow
0}\frac{1}{\Delta\omega}\int_{\omega_{k}-\Delta\omega/2}^{\omega_{k}+\Delta\omega/2}d\omega\,\frac{f_{\mathrm{DM}}(v_{\omega})}{m_{\phi}\,v_{\omega}}=\frac{f_{\mathrm{DM}}(v_{\omega_{k}})}{m_{\phi}\,v_{\omega_{k}}}\,,$
(C.21)

where $\omega_{k}=2\pi k/T_{\mathrm{int}}$. Hence, in agreement with Ref. [53,
13], the amplitude of the PSD would be independent of the integration time.

In the regime where the integration time $T_{\mathrm{int}}$ is shorter than
the coherence time of the signal $\tau_{c}$, the experimentalist would not be
able to resolve the full frequency spread of the ULDM signal. In the regime
where $T_{\mathrm{int}}\ll\tau_{c}$, this further implies that the ULDM signal
would be contained in one frequency bin. In this case, Eq. (C.13) would be
evaluated over the entire dark matter speed distribution. Thus, recalling that
$\Delta\omega=2\pi/T_{\mathrm{int}}$ and that $\left\langle
S_{\mathrm{DM},k}\right\rangle\propto 1/\Delta\omega$, we find that the
expected ULDM signal PSD would scale linearly with the integration time, which
agrees with Refs. [53, 13].

###  C.2 Rayleigh wave covariance matrix

The GGN covariance matrix can be calculated by recalling the statistical
properties of Rayleigh waves. As for the ULDM signal, we shall first focus on
the covariance between the real parts of the data’s DFT,

$\begin{split}\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle&=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\sum_{m,m^{\prime}=0}^{N-1}\left\langle\Phi^{(i,j)}_{\mathrm{GGN},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{GGN},m^{\prime}}\right\rangle
c_{m,k}\,c_{m^{\prime},k}\\\ &=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\sum_{m,m^{\prime}=0}^{N-1}\sum_{a}\frac{\langle\xi_{a}^{2}\rangle}{2}F_{\mathrm{GGN},a}^{(i,j)}F_{\mathrm{GGN},a}^{(i^{\prime},j^{\prime})}\cos(\omega_{a}\Delta
t(m-m^{\prime}))\,c_{m,k}\,c_{m^{\prime},k}\,,\end{split}$ (C.22)

where the second equality follows from the covariance of the GGN gradiometer
phase, shown in Eq. (B.26). We remind the reader that we use the convention

$\displaystyle
F_{\mathrm{GGN},a}^{(i,j)}=F_{\mathrm{GGN},a}^{(i)}-F_{\mathrm{GGN},a}^{(j)}\,,$
(C.23) $\displaystyle
F_{\mathrm{GGN},a}^{(i)}=\widetilde{A}_{a}\exp\left(-q\frac{\omega_{a}z_{i}}{c_{H}}\right)+\widetilde{B}_{a}\exp\left(-\frac{\omega_{a}z_{i}}{c_{H}}\right)\,,$
(C.24)

which was presented in the main body of the paper. Since the statistics of the
ULDM signal are equivalent to those of the GGN gradiometer phase, we follow
the calculations in Appendix C to determine a simplified analytical form of
Eq. (C.22). Specifically, following the steps in Eqs. (C.10)-(C.11), we may
express the RHS of Eq. (C.22) as

$\displaystyle\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{\pi}{2}\frac{1}{\Delta\omega}\frac{\left\langle\xi_{k}^{2}\right\rangle}{2}F_{\mathrm{GGN},k}^{(i,j)}F_{\mathrm{GGN},k}^{(i^{\prime},j^{\prime})}\,.$
(C.25)

To make contact with the power spectral density of vertical displacements
$S_{\xi}(\omega)$ (see Fig. 3), we relate the effective mean-squared vertical
displacement at angular frequency $\omega_{k}$ to $\langle
S_{\xi}(\omega)\rangle$ via

$\displaystyle\left\langle\xi_{k}^{2}\right\rangle$
$\displaystyle=\int_{\omega_{k}-\Delta\omega/2}^{\omega_{k}+\Delta\omega/2}d\omega\,\left\langle
S_{\xi}(\omega)\right\rangle\,.$ (C.26)

In light of the fact that the spectrum of vertical displacements varies
significantly over frequency scales greater than $\sim 10^{-2}$ Hz and that
$\Delta\omega=2\pi/T_{\mathrm{int}}\approx$6\text{\times}{10}^{-6}\text{\,}\mathrm{H}\mathrm{z}$$
in the experiments considered here, we may treat $S_{\xi}(\omega)$ as being
constant within the integral, such that

$\displaystyle\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{\pi}{4}\left\langle S_{\xi}(\omega_{k})\right\rangle
F_{\mathrm{GGN},k}^{(i,j)}F_{\mathrm{GGN},k}^{(i^{\prime},j^{\prime})}\,.$
(C.27)

Following the same mapping arguments presented for the entries of the ULDM
signal covariance matrix, the entries in the GGN covariance matrix take the
form

$\displaystyle\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{\pi}{4}\left\langle S_{\xi}(\omega_{k})\right\rangle
F_{\mathrm{GGN},k}^{(i,j)}F_{\mathrm{GGN},k}^{(i^{\prime},j^{\prime})}\,,$
(C.28) $\displaystyle\left\langle
I^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{\pi}{4}\left\langle S_{\xi}(\omega_{k})\right\rangle
F_{\mathrm{GGN},k}^{(i,j)}F_{\mathrm{GGN},k}^{(i^{\prime},j^{\prime})}\,,$
(C.29) $\displaystyle\left\langle
R^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=0\,.$ (C.30)

###  C.3 Atom shot noise covariance matrix

The atom shot noise covariance matrix can be calculated by recalling the
statistical properties of white noise. Using the statistical properties of
atom shot noise as presented in Eq. (3.24), we find that the covariance
between the real parts of the data’s DFT takes the form

$\begin{split}\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle&=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\sum_{m,m^{\prime}=0}^{N-1}\left\langle\Phi^{(i,j)}_{\mathrm{ASN},m}\Phi^{(i^{\prime},j^{\prime})}_{\mathrm{ASN},m^{\prime}}\right\rangle
c_{m,k}c_{m^{\prime},k}\\\ &=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\sum_{m,m^{\prime}=0}^{N-1}\left(\delta^{ii^{\prime}}+\delta^{jj^{\prime}}-\delta^{ij^{\prime}}-\delta^{ji^{\prime}}\right)\delta_{mm^{\prime}}\sigma_{\mathrm{Atom}}^{2}c_{m,k}c_{m^{\prime},k}\\\
&=\frac{(\Delta
t)^{2}}{T_{\mathrm{int}}}\left(\delta^{ii^{\prime}}+\delta^{jj^{\prime}}-\delta^{ij^{\prime}}-\delta^{ji^{\prime}}\right)\sigma_{\mathrm{Atom}}^{2}\sum_{m=0}^{N-1}c_{m,k}^{2}\,,\end{split}$
(C.31)

where in going to the last line we performed the sum over discrete frequencies
and used the property of the Kronecker delta. Recalling the orthogonality
relation for cosines, we find

$\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle=\frac{\Delta
t}{2}\sigma_{\mathrm{Atom}}^{2}\left(\delta^{ii^{\prime}}+\delta^{jj^{\prime}}-\delta^{ij^{\prime}}-\delta^{ji^{\prime}}\right)\,.$
(C.32)

By defining the PSD of atom shot noise in a single interferometer as $\langle
S_{\mathrm{ASN},k}\rangle=\Delta t\sigma_{\mathrm{Atom}}^{2}$, we can re-
express Eq. (C.32) as

$\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle=\frac{1}{2}\langle
S_{\mathrm{ASN},k}\rangle\left(\delta^{ii^{\prime}}+\delta^{jj^{\prime}}-\delta^{ij^{\prime}}-\delta^{ji^{\prime}}\right)\,.$
(C.33)

Following the same mapping arguments that we presented for the ULDM signal
covariance matrix, we find that the entries in the atom shot noise covariance
matrix take the form

$\displaystyle\left\langle
R^{(i,j)}_{k}R^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{1}{2}\langle
S_{\mathrm{ASN},k}\rangle\left(\delta^{ii^{\prime}}+\delta^{jj^{\prime}}-\delta^{ij^{\prime}}-\delta^{ji^{\prime}}\right)\,,$
(C.34) $\displaystyle\left\langle
I^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=\frac{1}{2}\langle
S_{\mathrm{ASN},k}\rangle\left(\delta^{ii^{\prime}}+\delta^{jj^{\prime}}-\delta^{ij^{\prime}}-\delta^{ji^{\prime}}\right)\,,$
(C.35) $\displaystyle\left\langle
R^{(i,j)}_{k}I^{(i^{\prime},j^{\prime})}_{k}\right\rangle$
$\displaystyle=0\,.$ (C.36)

##  Appendix D Asimov Test Statistics

In section 3 we defined the basic frequentist tool in the Asimov approach as

$\widetilde{\Theta}(\boldsymbol{\theta}_{\mathrm{sig}})=\langle\Theta(\boldsymbol{\theta}_{\mathrm{sig}})\rangle\,,$
(D.1)

where $\Theta(\boldsymbol{\theta}_{\mathrm{sig}})$ is defined in Eq. (4.9) and
the average is taken over data realizations. Recalling that we have assumed a
covariance matrix $\boldsymbol{\mathrm{\Sigma}}_{k}$ that can be split into
signal $\boldsymbol{\mathrm{S}}_{k}$ and noise $\boldsymbol{\mathrm{N}}_{k}$
contributions, the Asimov TS can be expressed as

$\left\langle\Theta\left(\boldsymbol{\theta}_{\mathrm{sig}}\right)\right\rangle=\left\langle\sum_{k=1}^{N-1}\left(\mathbf{d}_{k}^{T}\left[\mathbf{B}^{-1}_{k}-\mathbf{\Sigma}^{-1}_{k}\right]\mathbf{d}_{k}-\mathrm{ln}\left[\frac{|\boldsymbol{\Sigma}_{k}|}{|\mathbf{B}_{k}|}\right]\right)\right\rangle\,,$
(D.2)

where $\mathbf{d}_{k}$ is the multi-dimensional data vector. Because the
entries of each covariance matrix correspond to expectation values, the
average is taken only over the first term in Eq. (D.2). Since the data has
zero mean, we can define
$\mathbf{\Sigma}^{t}_{k}\equiv\left\langle\mathbf{d}_{k}\mathbf{d}_{k}^{T}\right\rangle$
as the true covariance matrix, such that

$\left\langle\mathbf{d}_{k}^{T}\left[\mathbf{B}_{k}^{-1}-\mathbf{\Sigma}_{k}^{-1}\right]\mathbf{d}_{k}\right\rangle=\mathrm{Tr}\left(\mathbf{\Sigma}^{t}_{k}\left[\mathbf{B}_{k}^{-1}-\mathbf{\Sigma}_{k}^{-1}\right]\right)\,.$
(D.3)

Using this result and recalling the trace-determinant identity for symmetric
matrices, Eq. (D.2) takes the form

$\displaystyle\widetilde{\Theta}\left(\boldsymbol{\theta}_{\mathrm{sig}}\right)$
$\displaystyle=\sum_{k=1}^{N-1}\mathrm{Tr}\left(\mathbf{\Sigma}^{t}_{k}\left[\mathbf{B}^{-1}_{k}-\mathbf{\Sigma}^{-1}_{k}\right]+\mathrm{ln}\left[\mathbf{B}_{k}\mathbf{\Sigma}_{k}^{-1}\right]\right)\,.$
(D.4)

Given a fixed background, the true covariance matrix can be defined as
$\mathbf{\Sigma}^{t}_{k}=\mathbf{S}^{t}_{k}+\mathbf{B}_{k}$, where
$\mathbf{S}^{t}_{k}$ is the true signal model. In the limit that the signal is
weaker than the background, i.e.
$\mathbf{\Sigma}_{k}^{-1}\approx\mathbf{B}_{k}^{-1}-\mathbf{B}_{k}^{-1}\mathbf{S}_{k}\mathbf{B}_{k}^{-1}$,
we can expand the natural logarithm of the matrix product close to the
identity matrix such that

$\displaystyle\widetilde{\Theta}\left(\boldsymbol{\theta}_{\mathrm{sig}}\right)$
$\displaystyle\approx\sum_{k=1}^{N-1}\mathrm{Tr}\left[\left(\mathbf{S}^{t}_{k}-\frac{1}{2}\mathbf{S}_{k}\right)\mathbf{B}^{-1}_{k}\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\right]\,.$
(D.5)

For setting upper limits on the ULDM-SM couplings $d_{\phi}$, we use the test
statistic $q\left(m_{\phi},d_{\phi}\right)$. The equivalent test statistic on
the Asimov data set can then be defined as

$\widetilde{q}\left(m_{\phi},d_{\phi}\right)=\begin{cases}\widetilde{\Theta}\left(\\{m_{\phi},d_{\phi}\right\\})-\widetilde{\Theta}\left(\\{m_{\phi},\hat{d_{\phi}}\\}\right)&d_{\phi}\geq\hat{d_{\phi}}\\\
0&d_{\phi}<\hat{d_{\phi}}\end{cases}\;,$ (D.6)

where $\hat{d_{\phi}}$ is the value of $d_{\phi}$ that maximises
$\widetilde{\Theta}\left(m_{\phi},d_{\phi}\right)$ at fixed $m_{\phi}$. In
this case, the appropriate Asimov dataset contains only the background, such
that $\hat{d_{\phi}}=0$, which implies that the true covariance matrix can be
written as $\mathbf{\Sigma}^{t}_{k}=\mathbf{B}_{k}$ \- i.e.
$\mathbf{S}^{t}_{k}=\mathbf{0}$. Therefore, in the weak signal limit, we may
express the test statistic for setting upper limits as:

$\widetilde{q}\left(m_{\phi},d_{\phi}\right)=\begin{cases}-\frac{1}{2}\sum_{k=1}^{N-1}\mathrm{Tr}\left[\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\right]&d_{\phi}\geq
0\\\ 0&d_{\phi}<0\end{cases}\;.$ (D.7)

###  D.1 Behaviour under long and short integration times

In practice, the sum in Eq. (D.5) is only to be performed around a small
frequency window delineated by $\omega_{\mathrm{min}}\lesssim m_{\phi}$ and
$\omega_{\mathrm{max}}\gtrsim m_{\phi}\left(1+v_{\mathrm{esc}}^{2}/2\right)$,
where $v_{\mathrm{esc}}$ is the dark matter escape velocity in the Milky Way.
In the regime where $T_{\mathrm{int}}$ is the longest timescale, i.e.
$T_{\mathrm{int}}\gg\tau_{c}$, the sum over the relevant frequency bins can be
approximated to an integral over frequency indices, such that, for
$d_{\phi}\geq 0$, Eq. (D.7) is approximately

$\displaystyle\widetilde{q}\left(m_{\phi},d_{\phi}\right)$
$\displaystyle\approx-\frac{1}{2}\int
dk\,\mathrm{Tr}\left[\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\right]$
(D.8) $\displaystyle\sim\int dk\left\langle
R_{\mathrm{DM},k}^{2}\right\rangle^{2}/\left\langle
R_{\mathrm{Noise},k}^{2}\right\rangle^{2}\,.$

In addition, since the normalised frequency band over which the sum is
performed is $\sim 10^{-6}$ and GGN is not expected to vary appreciably over
this normalised frequency range, we treat all backgrounds as being constant.
Hence, recalling the relation $2\pi dk/T_{\mathrm{int}}=m_{\phi}vdv$, we find

$\displaystyle\widetilde{q}\left(m_{\phi},d_{\phi}\right)$
$\displaystyle\sim\frac{m_{\phi}T_{\mathrm{int}}}{\langle
R^{2}_{\mathrm{Noise}}\rangle^{2}}\int v\,dv\,\left\langle
R_{\mathrm{DM},\omega}^{2}\right\rangle^{2}$ (D.9)
$\displaystyle\sim\frac{T_{\mathrm{int}}}{\langle
R^{2}_{\mathrm{Noise}}\rangle^{2}}\int\frac{dv}{m_{\phi}v}\,f^{2}_{\mathrm{SHM}}(v)$
$\displaystyle\sim T_{\mathrm{int}}\tau_{c}\,,$

where in going to the second line we used $\langle
R_{\mathrm{DM},\omega}^{2}\rangle\sim f_{\mathrm{SHM}}(v)/m_{\phi}v$, which
follows from the discussion in Appendix C.1.1, and subsequently recalled the
definition of the coherence time $\tau=2\pi/(m_{\phi}v_{0}^{2})$. Hence, as
derived via Bartlett’s method in Ref. [13], the sensitivity to ULDM would
scale like $(T_{\mathrm{int}}\tau_{c})^{-1/4}$.

In the regime where $T_{\mathrm{int}}\ll\tau_{c}$, the experiment is unable to
resolve the DM speed distribution. In this case, the signal is confined to a
single frequency bin, such that, for $d_{\phi}\geq 0$, Eq. (D.7) is
approximately

$\displaystyle\widetilde{q}\left(m_{\phi},d_{\phi}\right)$
$\displaystyle\approx-\frac{1}{2}\mathrm{Tr}\left[\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\mathbf{S}_{k}\mathbf{B}^{-1}_{k}\right]$
(D.10) $\displaystyle\sim\left\langle
R_{\mathrm{DM},k}^{2}\right\rangle^{2}/\left\langle
R_{\mathrm{Noise},k}^{2}\right\rangle^{2}\,.$

Following the analysis presented in Appendix C.1.1, we make use of the fact
that $\langle R_{\mathrm{DM},k}^{2}\rangle\sim 1/\Delta\omega\propto
T_{\mathrm{int}}$ to determine that the test statistic for setting upper
limits scales as $T_{\mathrm{int}}^{2}$. Hence, in the regime where
$T_{\mathrm{int}}\ll\tau_{c}$ the sensitivity to ULDM would scale like
$T_{\mathrm{int}}^{-1/2}$, in agreement with Ref. [13].

## References

  * [1] A. Peters, K.Y. Chung and S. Chu, _High-precision gravity measurements using atom interferometry_ , _Metrologia_ 38 (2001) 25.
  * [2] B. Canuel, F. Leduc, D. Holleville, A. Gauguet, J. Fils, A. Virdis et al., _Six-axis inertial sensor using cold-atom interferometry_ , _Phys. Rev. Lett._ 97 (2006) 010402.
  * [3] G. Rosi, F. Sorrentino, L. Cacciapuoti, M. Prevedelli and G.M. Tino, _Precision measurement of the Newtonian gravitational constant using cold atoms_ , _Nature_ 510 (2014) 518 [1412.7954].
  * [4] R.H. Parker, C. Yu, W. Zhong, B. Estey and H. Müller, _Measurement of the fine-structure constant as a test of the Standard Model_ , _Science_ 360 (2018) 191 [1812.04130].
  * [5] P. Asenbaum, C. Overstreet, M. Kim, J. Curti and M.A. Kasevich, _Atom-interferometric test of the equivalence principle at the ${10}^{-12}$ level_, _Phys. Rev. Lett._ 125 (2020) 191101.
  * [6] B. Stray et al., _Quantum sensing for gravity cartography_ , _Nature_ 602 (2022) 590.
  * [7] S. Dimopoulos, P.W. Graham, J.M. Hogan, M.A. Kasevich and S. Rajendran, _Gravitational Wave Detection with Atom Interferometry_ , _Phys. Lett. B_ 678 (2009) 37 [0712.1250].
  * [8] P.W. Graham, J.M. Hogan, M.A. Kasevich and S. Rajendran, _A New Method for Gravitational Wave Detection with Atomic Sensors_ , _Phys. Rev. Lett._ 110 (2013) 171102 [1206.0818].
  * [9] J. Ellis and V. Vaskonen, _Probes of gravitational waves with atom interferometers_ , _Phys. Rev. D_ 101 (2020) 124013 [2003.13480].