# Gravitational Interaction of Ultralight Dark Matter with Interferometers

Hyungjin Kim<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

Ultralight dark matter exhibits an order-one density fluctuation over the
spatial scale of its wavelength. These fluctuations gravitationally interact
with gravitational wave interferometers, leading to distinctive signals in
detectors. We investigate the ultralight dark matter-induced effects in the
gravitational wave interferometers. We perform a systematic computation of the
power spectrum of ultralight dark matter in interferometers. We show that the
ultralight dark matter-induced effect is most relevant for the interferometers
with long baseline and that it is only a sub-leading effect compared to the
estimated noise level in the case of Laser Interferometer Space Antenna or
future interferometers with an arm-length comparable to a few astronomical
units. Gravitational wave interferometers can then place upper limits on the
ultralight dark matter density in the solar system. We find that, under
certain assumptions, future interferometers with AU-scale arm-length might
probe the dark matter density a few hundred times the local dark matter
density, which is measured over a much larger spatial scale.

††preprint: DESY-23-085

oneΔ

###### Contents

  1. I Introduction
  2. II Ultralight dark matter
    1. II.1 Density fluctuation
    2. II.2 Metric fluctuation
  3. III Gravitational wave interferometers
    1. III.1 Interferometry response
      1. III.1.1 Michelson interferometer
      2. III.1.2 Time-delay interferometer
    2. III.2 ULDM signal
      1. III.2.1 Deterministic signal
      2. III.2.2 Stochastic signal
    3. III.3 Application
      1. III.3.1 LIGO
      2. III.3.2 Einstein Telescope
      3. III.3.3 LISA
      4. III.3.4 $\mu$Ares
      5. III.3.5 Big Bang Observatory
  4. IV Discussion
  5. V Conclusion
  6. A Power spectrum
  7. B Induced phase
  8. C Density fluctuation
  9. D Coherent signals in pulsar timing array

##  I Introduction

Ultralight dark matter (ULDM) remains an attractive dark matter candidate. Its
mass is smaller than several eV-scale, and it behaves like classical waves.
The QCD axion is one of the most theoretically interesting ultralight dark
matter candidates [1, 2, 3, 4, 5, 6]. Other interesting candidates often arise
from recent developments in dynamical solutions to the electroweak hierarchy
problem [7, 8, 9, 10, 11, 12, 13, 14]. Phenomenologically, ultralight dark
matter of extremely small mass $m\sim 10^{-21}\,{\rm eV}$ has also been
actively investigated as it changes the prediction of cold dark matter at
galactic and sub-galactic scales [15]. Being wave-like, such ultralight dark
matter candidates lead to interesting phenomenology in the early and late
universe. See Ref. [16] for a recent review on wave dark matter.

Recent numerical simulations of ultralight dark matter halo have observed an
order one density fluctuation all over galaxies [17, 18]. The fluctuation is
extended over the spatial scale given by the wavelength of dark matter. This
characteristic density fluctuation can be intuitively understood in terms of
quasiparticles, whose size ($\lambda$) and mass $(m_{\rm eff})$ are given as
its wavelength and the total mass enclosed within the volume of de Broglie
wavelength [19]:

$\displaystyle\lambda$ $\displaystyle=\frac{1}{m\sigma}\simeq 10^{2}\,{\rm
pc}\Big{(}\frac{10^{-22}\,{\rm eV}}{m}\Big{)}\Big{(}\frac{200\,{\rm
km/sec}}{\sigma}\Big{)},$ (1) $\displaystyle m_{\rm eff}$
$\displaystyle=\frac{\pi^{3/2}\rho_{0}}{(m\sigma)^{3}}\simeq 5\times
10^{4}M_{\odot}\Big{(}\frac{\rho_{0}}{{\rm
GeV/cm^{3}}}\Big{)}\Big{(}\frac{10^{-22}\,{\rm
eV}}{m}\Big{)}^{3}\Big{(}\frac{200\,{\rm km/sec}}{\sigma}\Big{)}^{3}.$ (2)

Here $m$ is the mass of ULDM, $\rho_{0}$ is the mean dark matter density, and
$\sigma$ is the velocity dispersion. For $m\ll\,{\rm eV}$, the size and the
mass of the quasiparticles could be astronomical.111The quasiparticle,
however, is not a particle in a traditional sense. The density fluctuation is
a result of wave interference, and it lasts only for a coherence time
$\tau=1/m\sigma^{2}$. In addition, the fluctuation of the mass within a fixed
volume in galaxies follows the exponential distribution in the case of ULDM,
while it is the Poisson distribution for the particle dark matter. These
features make the behavior of ULDM qualitatively different from that of
massive particle dark matter candidates, such as primordial black holes or
massive compact halo objects, which have been previously studied in [20, 21,
22, 23, 24].

These quasiparticles interact continuously with interferometers, designed to
detect gravitational waves (GWs) of astrophysical and cosmological origin.
Since they are engineered to measure an extremely small disturbance of
spacetime due to GWs, any interaction of quasiparticles with such
interferometers leaves some distinctive impacts on detectors. The main goal of
this work is, therefore, to characterize the impacts of ULDM on gravitational
wave interferometers. In particular, we ask two questions: (i) if ULDM signals
can be larger than the projected noise level in each detector and (ii) if the
dark matter density in the solar system can be probed only through the
gravitational interaction of order-one density fluctuations of ULDM with
detectors. We will address these questions by computing the signal power
spectrum of ULDM in detectors and comparing it with the noise power spectrum
of detectors.

Without detailed computation, we can already estimate the expected effects of
ULDM and partially answer one of the above two questions. Consider an
interferometer with two test masses. A useful quantity is a differential
acceleration, $\Delta a=a_{1}-a_{2}$, between two test masses induced by ULDM
density fluctuations. Suppose that a quasiparticle is located right next to
the first test mass. In this case, the differential acceleration is estimated
as

$\Delta a=a_{1}-a_{2}=\frac{Gm_{\rm eff}}{\lambda^{2}}-\frac{Gm_{\rm
eff}}{(L+\lambda)^{2}}\simeq\bar{a}\min\Big{(}1,\frac{2L}{\lambda}\Big{)},$

where $L$ is a typical size of the interferometer and typical acceleration
$\bar{a}$ is defined as

$\displaystyle\bar{a}\equiv\frac{Gm_{\rm
eff}}{\lambda^{2}}=\pi^{3/2}G\rho_{0}\lambda.$ (3)

As the mass $m$ decreases (hence wavelength and effective mass increases), the
differential acceleration increases until it saturates to
$\bar{a}=2G\rho_{0}L$ when $\lambda=2L$; when the size of a quasiparticle
becomes comparable to the size of the interferometer. At this saturation
point, we find $\Delta a\sim 10^{-28}\,{\rm m\,s^{-2}}$ for $L=4\,{\rm km}$
(LIGO/VIRGO), $\Delta a\sim 10^{-22}\,{\rm m\,s^{-2}}$ for $L=2.5{\rm M}\,{\rm
km}$ (LISA), and $\Delta a\sim 10^{-20}$ for $L=400{\rm M}\,{\rm km}$ (e.g.
$\mu$Ares strawman mission concept [25] and Fedderke et al [26]). Here we use
$\rho_{0}=0.4\,{\rm GeV}/{\rm cm^{3}}$.

This can be compared with the sensitivity of detectors. The strain noise power
spectrum can be converted into the root-mean-square fluctuation of
acceleration by $\Delta a\sim[S_{n}(2\pi f)^{4}L^{2}\Delta f]^{1/2}$. We find
$\Delta a=10^{-14}\,{\rm m\,s^{-2}}$ for LIGO with $f\sim\Delta f\sim 50\,{\rm
Hz}$, $\Delta a=10^{-16}\,{\rm m\,s^{-2}}$ for LISA with $f\sim\Delta f\sim
0.1{\rm mHz}$, and $\Delta a=10^{-17}\,{\rm m\,s^{-2}}$ for $\mu$Ares proposal
with $f\sim\Delta f\sim\textrm{a few}\times 10^{-7}{\rm Hz}$. For ground-based
interferometers, the ULDM-induced noise is irrelevant as it is many orders of
magnitude smaller than detector noises. For space-borne interferometers, the
ULDM effect is a few orders of magnitude smaller than the mission-required
detector noise level. We may conclude from this comparison that the ULDM
effect in interferometers through gravitational interaction is unlikely to
interfere with the operation of interferometers for the detection of GWs. For
a more correct comparison, however, we will need to compute the power spectrum
of ultralight dark matter-induced effects since the above estimation on
$\Delta a$ does not carry any spectral information about the ULDM effects in
the detector.

We summarize the main results here. From the detailed computation of the ULDM
power spectrum, we find that the spectrum generally consists of two
distinctive frequency components: one at $\omega=2m$ and the other at
$\omega<m\sigma^{2}$. This has interesting phenomenological implications as it
offers different ways to search for ultralight dark matter signals. In Figure
1, we show the low-frequency part ($\omega<m\sigma^{2}$) of the ULDM signal.
The ULDM noise is sub-dominant even for the interferometers with arm-lengths
larger than a few million km, such as LISA and other proposals with arm-length
comparable to astronomical units, confirming our order-of-magnitude estimation
above. The interferometers then can be used to place an upper limit on the
dark matter content in the solar system. Two distinct frequency components
offer multiple ways for ULDM searches in interferometers. The $\omega=2m$
component behaves similarly to a deterministic signal, which can be searched
with the matched filter. The low-frequency ($\omega<m\sigma^{2}$) component,
on the other hand, behaves similarly to the stochastic background, which can
be searched by cross-correlating detector outputs if more than two detectors
are available. The resulting constraints and projections on $\rho/\rho_{0}$
are shown in Figure 2, suggesting that the density of ULDM in the solar
system, $\rho\sim{\cal O}(10^{2})\rho_{0}$, could be probed with GW
interferometers only through gravitational interaction.

Figure 1: The ultralight dark matter-induced noise for (left) LISA and (right)
$\mu$Ares proposal [25]. The red line is the required noise level for the
mission. The other colored solid lines are the ultralight dark matter-induced
noise for given masses. In both cases, the shown is the strain noise power
spectrum in the Michelson-like TDI-X variable. Note that, for $\mu$Ares, the
low-frequency noise could be dominated by the gravity gradient noise as
discussed in Ref. [27]. See the main text for details. Figure 2: The
constraints and projections on dark matter density in the solar system. Here
$\rho_{0}=0.4\,{\rm GeV}/{\rm cm^{3}}$ is the local dark matter density
measured over $\gtrsim{\cal O}(100)$ pc scale (see e.g. [28, 29] for reviews).
The solid lines are obtained from the matched filter for the coherently
oscillating dark matter signal at $\omega=2m$, while the dashed lines are
obtained by cross-correlating signals in two detectors at
$\omega<m\sigma^{2}$. Note that the current design of LISA does not include
two detectors, but other space-borne interferometer proposals, such as TianQin
or Taiji, would potentially allow the cross-correlation searches [30]. For the
matched filter search of the coherently oscillating signal, we take for
simplicity vanishing dark matter velocity $v_{0}=0$; we check this
approximation does not change the numerical result in a significant way. For
the stochastic signal search with the cross-correlation, we take an average of
the signal over twelve months. See the main text for details.

The work is organized as follows. In section II, we review the basic
statistical properties of ultralight dark matter and compute the spectrum of
density contrast. This will establish a basis for the detailed investigation
of the response of GW interferometers to the ULDM density fluctuations. In
section III, we discuss how GW interferometers respond to the ULDM density
fluctuation, and how the ultralight dark matter signal can be searched with
the matched filter and with cross-correlation of detector outputs in cases
where there are more than two detectors. In section IV, we discuss some of the
assumptions that might affect our results. We conclude in section V. The
natural unit ($c=\hbar=1$) is used in this work.

##  II Ultralight dark matter

We review the basic statistical properties of ultralight dark matter in this
section. We consider a scalar field minimally coupled to gravity without self-
interaction. The ultralight dark matter field operator is given as

$\hat{\phi}(x)=\sum_{i}\frac{1}{\sqrt{2mV}}(a_{i}e^{-ik_{i}\cdot
x}+a_{i}^{\dagger}e^{ik_{i}\cdot x}),$

where $a_{i}$ and $a^{\dagger}_{i}$ are annihilation and creation operator of
the plane wave mode $i$, satisfying $[a_{i},a^{\dagger}_{j}]=\delta_{ij}$. The
continuum expression can be obtained by $\sum_{i}\to V\int d^{3}k/(2\pi)^{3}$,
$a_{i}\to a(\vec{k})/\sqrt{V}$,
$[a(\vec{k}),a^{\dagger}(\vec{q})]=(2\pi)^{3}\delta^{(3)}(\vec{k}-\vec{q})$.
For later convenience, we stick to the above discrete convention. Any
deformation of wave function due to the gravitational potential of the Sun is
ignored since this effect is at the level of ${\cal O}(1)\%$ for the halo dark
matter [31].

The statistical property of ultralight dark matter is defined by the density
operator. The density operator is $\hat{\rho}=\prod_{i}\hat{\rho}_{i}\otimes$
where $\hat{\rho}_{i}$ is the density operator of each mode, given by [31]

$\displaystyle\hat{\rho}_{i}=\int
d^{2}\alpha_{i}\,p(\alpha_{i})|\alpha_{i}\rangle|\alpha_{i}|.$ (4)

Here $|\alpha_{i}\rangle$ is the coherent state, satisfying
$a_{i}|\alpha_{i}\rangle=\alpha_{i}|\alpha_{i}\rangle$
($\alpha_{i}\in{\mathbb{C}}$), and the quasi-probability distribution
$p(\alpha_{i})$ is given as

$\displaystyle p(\alpha_{i})=\frac{1}{\pi
f_{i}}\exp\left[-\frac{|\alpha_{i}|^{2}}{f_{i}}\right].$ (5)

The $f_{i}$ is the mean occupation number for the mode $i$. This quasi-
probability distribution function $p(\alpha_{i})$ corresponds to the Rayleigh
distribution for the amplitude, $|\alpha_{i}|$, and the uniform distribution
for the phase, ${\rm arg}(\alpha_{i})$ [32, 33, 34]. With this description,
any $n$-point correlation function of an operator $\hat{\phi}$ can be easily
computed by the trace operation, $\langle{\cal O}\rangle={\rm Tr}({\cal
O}\hat{\rho})$. In the large occupation number limit, the operators
$(a_{i},a_{i}^{\dagger})$ can be replaced by commuting complex random number
$(\alpha_{i},\alpha_{i}^{*})$, whose probability distribution is given by (5).

It is useful to note the ensemble averages of the creation and annihilation
operators:

$\displaystyle\langle a_{i}a^{\dagger}_{j}\rangle=\langle
a_{i}^{\dagger}a_{j}\rangle=f_{i}\delta_{ij},\quad\textrm{and}\quad\langle
a_{i}a_{j}a_{k}^{\dagger}a_{\ell}^{\dagger}\rangle=f_{i}f_{j}(\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk}).$
(6)

Here we approximate $a_{i}$ and $a_{i}^{\dagger}$ as a commuting random
complex number $\alpha_{i}$ and $\alpha_{i}^{*}$, whose underlying
distribution follows the quasi-probability distribution $p(\alpha_{i})$. The
ultralight dark matter field defined with (5) is a Gaussian random field, and
therefore, the ensemble average of any operators with $N(a)-N(a^{\dagger})\neq
0$ vanishes, where $N(a)$ and $N(a^{\dagger})$ are the number of annihilation
and creation operator, respectively.

The energy density of the ultralight dark matter is

$\displaystyle\rho_{\phi}$
$\displaystyle=\frac{1}{2}\left[\dot{\phi}^{2}+(\nabla\phi)^{2}+m^{2}\phi^{2}\right]$
$\displaystyle=\frac{1}{2V}\sum_{i,j}\frac{1}{\sqrt{2\omega_{i}}\sqrt{2\omega_{j}}}\Big{[}W_{ij}^{-}\Big{(}a_{i}a_{j}e^{-i(k_{i}+k_{j})\cdot
x}+a^{\dagger}_{i}a^{\dagger}_{j}e^{i(k_{i}+k_{j})\cdot
x}\Big{)}+W_{ij}^{+}\Big{(}a_{i}a_{j}^{\dagger}e^{-i(k_{i}-k_{j})\cdot
x}+a^{\dagger}_{i}a_{j}e^{i(k_{i}-k_{j})\cdot x}\Big{)}\Big{]}.$ (7)

For notational simplicity, we have introduced $W^{\pm}_{ij}$, defined as

$\displaystyle
W^{\pm}_{ij}=m^{2}\pm(\omega_{k_{i}}\omega_{k_{j}}+\vec{k}_{i}\cdot\vec{k}_{j})\approx
2m^{2}\times\begin{cases}1&\textrm{for $W_{ij}^{+}$}\\\
-\frac{1}{4}(\vec{v}_{i}+\vec{v}_{j})^{2}&\textrm{for
$W_{ij}^{-}$}\end{cases}$ (8)

where the second expression holds in the non-relativistic limit. The mean
energy density in the non-relativistic limit is

$\displaystyle\langle\rho_{\phi}\rangle\equiv\rho_{0}\approx\frac{m}{V}\sum_{i}f_{i}=m\int\frac{d^{3}k}{(2\pi)^{3}}f(\vec{k})=m\int
d^{3}vf(\vec{v}),$ (9)

By taking the continuum limit, we reproduce the usual expression for the local
dark matter density. In the last line, we change the integration variable to
the velocity. The above expression fixes the normalization of the mean
occupation number $f(\vec{v})$ used in this work.

###  II.1 Density fluctuation

Let us consider the density fluctuation,
$\delta\rho_{\phi}=\rho_{\phi}-\rho_{0}$. The mean value vanishes
$\langle\delta\rho_{\phi}\rangle=0$. The correlator is

$\langle\delta\rho_{\phi}(x)\delta\rho_{\phi}(y)\rangle=\int\frac{d^{4}k}{(2\pi)^{4}}e^{-ik\cdot(x-y)}P_{\delta\rho}(k),$

where the power spectrum is defined as
$\langle\delta\rho_{\phi}(k)\delta\rho_{\phi}^{*}(k^{\prime})\rangle=(2\pi)^{4}\delta^{(4)}(k-k^{\prime})P_{\delta\rho}(k)$
and given by

$\displaystyle P_{\delta\rho}(\omega,\vec{k})=$
$\displaystyle\frac{(2\pi)^{4}}{2}\int
d\Pi_{1}d\Pi_{2}f(\vec{p}_{1})f(\vec{p}_{2})$ $\displaystyle\times$
$\displaystyle\Big{[}(W^{-}_{12})^{2}\big{(}\delta^{(4)}(k-p_{1}-p_{2})+\delta^{(4)}(k+p_{1}+p_{2})\big{)}+(W^{+}_{12})^{2}\big{(}\delta^{(4)}(k-p_{1}+p_{2})+\delta^{(4)}(k+p_{1}-p_{2})\big{)}\Big{]}.$
(10)

This can be straightforwardly computed by using (5) and (7). Here
$d\Pi=[d^{3}p/(2\pi)^{3}](2E)^{-1}$ is the Lorentz invariant phase space
measure. In the frequency space, the correlator is

$\displaystyle\langle\widetilde{\delta\rho}_{\phi}(\omega,\vec{x})\widetilde{\delta\rho}_{\phi}^{*}(\omega^{\prime},\vec{y})\rangle=(2\pi)\delta(\omega-\omega^{\prime})S_{\delta\rho}(\omega,\vec{L}),$
(11)

where $\vec{L}=\vec{x}-\vec{y}$ and the power spectrum
$S_{\delta\rho}(\omega,\vec{L})$ is given by

$\displaystyle
S_{\delta\rho}(\omega,\vec{L})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{L}}P_{\delta\rho}(\omega,\vec{k}).$
(12)

This power spectrum $S_{\delta\rho}(\omega,\vec{L})$ in the frequency space
turns out to be more useful for the study of the response of GW
interferometers to the ULDM fluctuations.

To illustrate the statistical properties of density contrast, we explicitly
compute $S_{\delta}(\omega,\vec{0})$. We consider the normal distribution for
dark matter velocity distribution,

$\displaystyle
f(\vec{v})=\frac{\rho_{0}/m}{(2\pi\sigma^{2})^{3/2}}\exp\left[-\frac{(\vec{v}-\vec{v}_{0})^{2}}{2\sigma^{2}}\right],$
(13)

where $\vec{v}_{0}$ is the mean dark matter velocity, and $\sigma$ is the
velocity dispersion. We work in the rest frame of the GW detectors. In
general, $\vec{v}_{0}\neq 0$ in this frame, but for the sake of simplicity, we
take the isotropic limit $\vec{v}_{0}=0$ for now. Normalization of
$f(\vec{v})$ is fixed by (9). Using (10) and (12), we find

$\displaystyle
S_{\delta}(\omega,\vec{0})=\tau\Big{[}\sigma^{4}A_{\delta}(\omega)+B_{\delta}(\omega)\Big{]},$
(14)

where $\delta=\delta\rho_{\phi}/\rho_{0}$ is density contrast, and
$\tau=1/m\sigma^{2}=1/\omega_{0}$ is the coherence time. In the isotropic
limit ($v_{0}/\sigma\ll 1$), the $A_{\delta}(\omega)$ and $B_{\delta}(\omega)$
are given by

$\displaystyle A_{\delta}(\omega)$
$\displaystyle=\frac{5\pi}{32}\Big{(}\frac{\bar{v}}{\sigma}\Big{)}^{8}\exp\Big{[}-\frac{\bar{v}^{2}}{\sigma^{2}}\Big{]}\theta(\bar{v}^{2})+(\omega\to-\omega),$
(15) $\displaystyle B_{\delta}(\omega)$
$\displaystyle=2|\bar{\omega}|K_{1}(|\bar{\omega}|),$ (16)

where $K_{n}(x)$ is the modified Bessel function. Here we introduce, for
notational simplicity,

$\displaystyle\bar{v}^{2}(\omega)=(\omega/m-2),\quad\textrm{and}\quad\bar{\omega}=\omega/(m\sigma^{2})=\omega/\omega_{0}.$
(17)

The spectrum $S_{\delta}(\omega,\vec{0})$ for $\vec{v}_{0}=0$ and $\sigma=0.1$
is shown in the left panel of Figure 3.

Figure 3: (Left) the density contrast spectrum
$S_{\delta}(\omega)=S_{\delta}(\omega,\vec{0})$ in the isotropic limit
$v_{0}=0$. (Right) the spectrum for the gradient of the metric fluctuation,
$S_{\nabla\Psi}(\omega)=\delta^{ij}S_{\nabla\Psi}^{ij}(\omega,\vec{0})$. Both
spectra exhibit two distinctive frequency components: one at $\omega=2m$ and
another at $\omega<m\sigma^{2}$. For demonstration, we choose unrealistically
large velocity dispersion $\sigma=0.1$. The narrow peak at $\omega=2m$ is
suppressed by $\sigma^{4}$ compared to the smooth low-frequency spectrum. For
$S_{\nabla\Psi}$, there is a logarithmic divergence at low frequencies due to
the long-range nature of gravitational force.

This computation reveals an interesting feature in the ULDM density
fluctuation. The spectrum contains two distinctive frequency components: a
narrow peak at $\omega=2m$ represented by $A_{\delta}(\omega)$ and a broad
low-frequency spectrum at $\omega<m\sigma^{2}$ represented by
$B_{\delta}(\omega)$. Both of them have a width $\Delta\omega\sim
m\sigma^{2}$. They can be understood by investigating the energy density of a
single mode, $\phi(t,\vec{x})=\phi_{0}\cos(\omega t-\vec{k}\cdot\vec{x})$,

$\rho_{\phi}=\frac{1}{2}\Big{[}\dot{\phi}^{2}+(\nabla\phi)^{2}+m^{2}\phi^{2}\Big{]}\approx\rho_{0}\Big{[}1-v^{2}\cos(2mt+\varphi)\Big{]}.$

It consists of dc mode, which is what constitutes the low-frequency part
($\omega<m\sigma^{2}$) for a realistic multi-mode density contrast, and a
coherently oscillating part which constitutes a spectrum at $\omega=2m$ with
its amplitude suppressed by $v^{2}$. This $v^{2}$ suppression is responsible
for a relative suppression of $\sigma^{4}$ in front of $A_{\delta}(\omega)$ in
(14). These two distinctive frequency components in the density contrast
remain in the ULDM-induced signal/noise power spectrum for the gravitational
detectors, offering multiple ways to search for ultralight dark matter at a
different mass range.

###  II.2 Metric fluctuation

ULDM fluctuations source metric perturbations. The metric is given by

$\displaystyle ds^{2}=(1+2\Phi)dt^{2}-(1-2\Psi)dx^{2}.$ (18)

These metric fluctuations perturb the position of test masses in
interferometers and induce a time delay of light. Suppose a test mass at the
position $\vec{x}$. The test mass acceleration and the metric fluctuations
satisfy

$\displaystyle\vec{a}$ $\displaystyle=-\nabla\Phi,$ (19)
$\displaystyle\nabla^{2}\Psi$ $\displaystyle=4\pi G\delta\rho_{\phi},$ (20)
$\displaystyle 6\ddot{\Psi}+2\nabla^{2}(\Phi-\Psi)$ $\displaystyle=24\pi
G\delta P_{\phi}.$ (21)

Here $\delta P_{\phi}$ is the pressure perturbation of the scalar field
$\phi$. The pressure of the field is
$P_{\phi}=\dot{\phi}^{2}/2-(\nabla\phi)^{2}/6-m^{2}\phi^{2}/2$. Statistical
properties of metric fluctuation as well as test mass position can be derived
from those of ULDM by solving the above equations.

The relation between test mass acceleration, metric perturbations, and ULDM
fluctuations becomes transparent in Fourier space:

$\displaystyle\tilde{a}^{i}(k)$ $\displaystyle=-ik^{i}\tilde{\Phi}(k),$ (22)
$\displaystyle\tilde{\Psi}(k)$ $\displaystyle=-\frac{4\pi
G}{k^{2}}\widetilde{\delta\rho}_{\phi}(k),$ (23)
$\displaystyle\tilde{\Phi}(k)$ $\displaystyle=-\frac{4\pi
G}{k^{2}}\left[\Big{(}1-\frac{3\omega^{2}}{k^{2}}\Big{)}\widetilde{\delta\rho}_{\phi}(k)+3\widetilde{\delta
P}_{\phi}(k)\right].$ (24)

From above, one can derive the relation between power spectra, e.g.
$P_{\Psi}(k)=(4\pi G/k^{2})^{2}P_{\delta\rho}(k)$ and
$P_{a}^{ij}(k)=k^{i}k^{j}P_{\Phi}(k)$, where the test mass acceleration power
spectrum is defined as
$\langle\tilde{a}^{i}(k)\tilde{a}^{j*}(k^{\prime})\rangle=(2\pi)^{4}\delta^{(4)}(k-k^{\prime})P^{ij}_{a}(k)$.

Among others, the power spectrum of $\nabla\Psi$ turns out to be particularly
useful for the later discussion. In the frequency space, it is given by

$\displaystyle
S_{\nabla\Psi}^{ij}(\omega,\vec{L})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{L}}k^{i}k^{j}P_{\Psi}(k)=\bar{a}^{2}\tau\big{[}\sigma^{4}A^{ij}+B^{ij}\big{]}$
(25)

where $A^{ij}$ and $B^{ij}$ are

$\displaystyle A^{ij}(\omega,\vec{L})$
$\displaystyle=16\Big{(}\frac{\bar{v}}{\sigma}\Big{)}^{6}\exp\Big{[}-\frac{v_{0}^{2}+\bar{v}^{2}}{\sigma^{2}}\Big{]}\theta(\bar{v}^{2})\left[\delta^{ij}\frac{I_{3}(X)}{X^{3}}+X^{i}X^{j}\frac{I_{4}(X)}{X^{4}}\right]+(\omega\to-\omega),$
(26) $\displaystyle B^{ij}(\omega,\vec{L})$
$\displaystyle=\frac{2}{\sqrt{\pi}}\int_{-\infty}^{\infty}ds\,\frac{e^{-is\bar{\omega}}}{\sqrt{1+s^{2}}}\left[\delta^{ij}\frac{{\rm
erf}(Y)-G(Y)}{Y}+\frac{Y^{i}Y^{j}}{Y^{3}}\big{(}3G(Y)-{\rm
erf}(Y)\big{)}\right],$ (27)

with $G(X)=[1/(2X^{2})]({\rm erf}(X)-2Xe^{-X^{2}}/\sqrt{\pi})$ and

$\displaystyle\vec{X}=2\frac{\bar{v}}{\sigma}\Big{(}\frac{\vec{v}_{0}}{\sigma}+i\vec{L}_{\lambda}\Big{)},\qquad\vec{Y}=\frac{1}{\sqrt{1+s^{2}}}\Big{(}\frac{s\vec{v}_{0}}{\sigma}+\vec{L}_{\lambda}\Big{)},\qquad\vec{L}_{\lambda}=m\sigma\vec{L}.$
(28)

The $\bar{v}$ and $\bar{\omega}$ are defined in (17). Here $I_{n}$ is the
modified Bessel function of the first kind. Recall that $\tau=1/m\sigma^{2}$,
$m_{\rm eff}=\pi^{3/2}\rho_{0}/(m\sigma)^{3}$, and $\lambda=(m\sigma)^{-1}$
and the typical acceleration $\bar{a}=Gm_{\rm eff}/\lambda^{2}$ . We will see
shortly that the response of the GW interferometer depends mostly on
$\nabla\Psi$, and the above power spectrum can be used as a basic building
block to study the response of interferometers to the ultralight dark matter
fluctuations.

The structure of the spectrum is similar to the density contrast spectrum (14)
since $\Psi$ inherits its properties from ULDM fluctuations. Especially,
$A^{ij}(\omega,\vec{L})$ is centered at $\omega=2m$, while
$B^{ij}(\omega,\vec{L})$ is spread over $\omega\lesssim m\sigma^{2}$. In the
isotropic ($v_{0}=0$) and the long-wavelength limit ($m\sigma L\ll 1$), the
$A^{ij}$ and $B^{ij}$ are simplified to

$\displaystyle A_{0}^{ij}(\omega)$
$\displaystyle\approx\frac{1}{3}\Big{(}\frac{\bar{v}}{\sigma}\Big{)}^{6}\exp\Big{[}-\frac{\bar{v}^{2}}{\sigma^{2}}\Big{]}\theta(\bar{v}^{2})\Big{[}\delta^{ij}-\frac{(m\bar{v}L)^{2}}{4}(\delta^{ij}+2\hat{L}^{i}\hat{L}^{j})\Big{]}+(\omega\to-\omega),$
(29) $\displaystyle B_{0}^{ij}(\omega)$
$\displaystyle\approx\frac{16}{3\pi}\Big{[}\delta^{ij}K_{0}(|\bar{\omega}|)-\frac{(m\sigma
L)^{2}}{5}(\delta^{ij}+2\hat{L}^{i}\hat{L}^{j})|\bar{\omega}|K_{1}(|\bar{\omega}|)\Big{]}.$
(30)

The spectrum
$S_{\nabla\Psi}(\omega)=\delta^{ij}S_{\nabla\Psi}^{ij}(\omega,\vec{0})$ is
shown in the right panel of Figure 3. Naively, this can be considered as a
test mass acceleration power spectrum. While the overall structure of the
spectrum is similar, the acceleration spectrum logarithmically diverges at
small $\omega$. More specifically,
$B_{0}^{ij}\propto\delta^{ij}\ln(1/\bar{\omega})$ for $\bar{\omega}\ll 1$, and
this logarithmic divergence is related to the long-range nature of the
gravitational interaction. The derivation of $B^{ij}(\omega,\vec{L})$ is given
in Appendix A.

##  III Gravitational wave interferometers

Figure 4: The laser begins at $t_{0}$ at TM0. It propagates to TM1 and returns
to TM0 at time $t$. During one round-trip, the laser acquires a phase
$e^{2ik_{L}L}$ without ultralight dark matter. With ultralight dark matter,
additional phase $e^{i\Delta\phi}=e^{2i\omega_{L}(\delta t+\delta L)}$ is
induced because of (i) time-delay $\delta t$ of the laser light and (ii) the
fluctuation of the distance between two test masses $\delta L$.

To investigate how each interferometer responds to the ULDM fluctuation, we
first consider a simple setup shown in Figure 4 where a phase-coherent laser
begins from TM0, reaches TM1, and returns to TM0 again at $t$. The phase of
the laser at time $t$ is given by

$e^{-i\omega_{L}t_{0}}=e^{-i\omega_{L}t}e^{2ik_{L}L}e^{i\Delta\phi},$

since the phase along the null trajectory does not change. Here
$(\omega_{L},\vec{k}_{L})$ is the laser frequency and wave vector, and $t_{0}$
is the time at which the laser begins from TM0. Since $t_{0}=t-2L$, the laser
acquires additional phase $e^{2ik_{L}L}$ during this one round-trip.

Ultralight dark matter introduces an additional phase in two different ways:
it introduces the time delay of laser light and it perturbs the test mass
position. The additional phase $e^{i\Delta\phi}$ is introduced as a result,
and it is given by

$\displaystyle\Delta\phi(t)=2\omega_{L}(\delta t+\delta L)$ (31)

where the test mass perturbation $\delta L$ and the time delay $\delta t$ are
given as

$\displaystyle\delta L$
$\displaystyle=+\frac{1}{2}\hat{x}\cdot\Big{[}\big{(}\delta\vec{x}_{1}(t-L)-\delta\vec{x}_{0}(t-2L)\big{)}+\big{(}\delta\vec{x}_{1}(t-L)-\delta\vec{x}_{0}(t)\big{)}\Big{]},$
(32) $\displaystyle\delta t$
$\displaystyle=-\frac{1}{2}\left[\int^{t}_{t-L}dt^{\prime}+\int^{t-L}_{t-2L}dt^{\prime}\right][\Psi(t^{\prime},\vec{x}(t^{\prime}))+\Phi(t^{\prime},\vec{x}(t^{\prime}))],$
(33)

where $\delta\vec{x}_{0}$ and $\delta\vec{x}_{1}$ are the perturbation of the
position of TM0 and TM1 around their nominal positions, $\vec{x}_{0}$ and
$\vec{x}_{1}$. Note that the argument $\vec{x}(t^{\prime})$ in the metric
perturbation must be computed along the null trajectory. The above expression
is obtained by solving null geodesic with the metric perturbation;
specifically, one finds $t_{0}=t-2(L+\delta L)-2\delta t$.

The additional phase $\Delta\phi$ can be expanded in terms of metric
perturbations in the Fourier space. In particular, it can be approximated as

$\displaystyle\widetilde{\Delta\phi}(\omega)\approx\pm\frac{\omega_{L}}{\omega^{2}}e^{i\omega
L}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}\big{[}i\vec{k}\cdot\vec{D}(\hat{n},L)\big{]}\tilde{\Psi}(k),$
(34)

where $\hat{n}=(\vec{x}_{1}-\vec{x}_{0})/L$ and the vector
$\vec{D}(\hat{n},L)=2\hat{n}[e^{i\vec{k}\cdot\hat{n}L}-\cos(\omega L)]$
encodes the response of this simple system to the ULDM fluctuations. The
detailed derivation of the expression is given in Appendix B.

###  III.1 Interferometry response

####  III.1.1 Michelson interferometer

Consider now the Michelson interferometer. The Michelson interferometer
consists of two arms. The laser combines at the asymmetric output, and the
resulting power is $P=P_{0}\sin^{2}(k_{L}\Delta L+\Delta\phi_{\rm Mich})$
where $\Delta\phi_{\rm Mich}=\frac{1}{2}(\Delta\phi_{x}-\Delta\phi_{y})$ with
the additional phase $\Delta\phi_{x,y}$ acquired by the laser that has
propagated along $x-$ and $y-$arm, respectively. The expression for
$\Delta\phi_{x,y}$ is given in (34). In the frequency space, the phase is

$\displaystyle\widetilde{\Delta\phi}_{\rm
Mich}(\omega)=\pm\frac{\omega_{L}}{\omega^{2}}e^{i\omega
L}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}\big{[}i\vec{k}\cdot\vec{D}_{\rm
Mich}(\hat{x},\hat{y},L)\big{]}\tilde{\Psi}(k)$ (35)

where $\vec{D}_{\rm
Mich}(\hat{x},\hat{y};L)=\frac{1}{2}[\vec{D}(\hat{x},L)-\vec{D}(\hat{y},L)]$
is the response vector of the Michelson interferometer, and $\vec{x}_{0}$ is
the beam splitter position. For simplicity, we assume an equal arm-length.

The power spectrum for differential arm-length $\widetilde{\Delta
L}/L=\widetilde{\Delta\phi}_{\rm Mich}/(\omega_{L}L)$ can be directly obtained
from the above expression, and it is given by

$\displaystyle S^{\rm ULDM}_{\Delta
L/L}(\omega)=\frac{1}{\omega^{4}L^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}|\vec{k}\cdot\vec{D}_{\rm
Mich}|^{2}P_{\Psi}(k)=\frac{\bar{a}^{2}}{\omega^{4}L^{2}}\tau\Big{[}\sigma^{4}A_{\rm
Mich}+B_{\rm Mich}\Big{]}.$ (36)

The detector response is encoded in the function $A_{\rm Mich}$ and $B_{\rm
Mich}$. Note again the similar structure of the power spectrum compared to the
density contrast and metric fluctuation $\nabla\Psi$. The $A_{\rm Mich}$ and
$B_{\rm Mich}$ can be decomposed in terms of $A^{ij}$ and $B^{ij}$ in (25) as
follows:

$\displaystyle A_{\rm Mich}(\omega;\hat{x},\hat{y},L)=$ $\displaystyle
2\Big{[}\frac{1+\cos^{2}(\omega
L)}{2}\big{(}A_{R}^{11}(\vec{0})+A_{R}^{22}(\vec{0})\big{)}-\cos(\omega
L)\big{(}A_{R}^{11}(\vec{x}_{1})+A_{R}^{22}(\vec{x}_{2})\big{)}$
$\displaystyle-\big{(}A^{12}_{R}(\vec{x}_{12})+\cos^{2}(\omega
L)A^{12}_{R}(\vec{0})-\cos(\omega
L)(A^{12}_{R}(\vec{x}_{1})+A^{12}_{R}(\vec{x}_{2}))\big{)}\Big{]},$ (37)
$\displaystyle B_{\rm Mich}(\omega;\hat{x},\hat{y},L)=$ $\displaystyle
2\Big{[}\frac{1+\cos^{2}(\omega
L)}{2}\big{(}B_{R}^{11}(\vec{0})+B_{R}^{22}(\vec{0})\big{)}-\cos(\omega
L)\big{(}B_{R}^{11}(\vec{x}_{1})+B_{R}^{22}(\vec{x}_{2})\big{)}$
$\displaystyle-\big{(}B^{12}_{R}(\vec{x}_{12})+\cos^{2}(\omega
L)B^{12}_{R}(\vec{0})-\cos(\omega
L)(B^{12}_{R}(\vec{x}_{1})+B^{12}_{R}(\vec{x}_{2}))\big{)}\Big{]},$ (38)

where $A_{R}^{ab}=\hat{x}_{a}^{i}\hat{x}_{b}^{j}A_{R}^{ij}$,
$B_{R}^{ab}=\hat{x}_{a}^{i}\hat{x}_{b}^{j}B_{R}^{ij}$,

$A_{R}^{ij}(\vec{x})=\frac{1}{2}\Big{[}A^{ij}(\vec{x})+A^{ij}(-\vec{x})\Big{]},\quad\textrm{and}\quad
B_{R}^{ij}(\vec{x})=\frac{1}{2}\Big{[}B^{ij}(\vec{x})+B^{ij}(-\vec{x})\Big{]}.$

Here $(\hat{x}_{1},\hat{x}_{2})=(\hat{x},\hat{y})$ denote the unit vector to
each arm. The dependence of $A_{R}^{ab}$ on $\omega$ is suppressed in the
above expression for notational simplicity. This power spectrum can be
compared with the noise power spectrum in the current and future ground-based
interferometers, such as LIGO and Einstein Telescope. While we only discuss
the Michelson interferometer here, the response of Fabry-Pérot Michelson
interferometer is the same up to the overall multiplication factor.

The response functions can be simplified in the isotropic ($v_{0}/\sigma\ll
1$) and long-wavelength limit ($m\sigma L\ll 1$). In such limits, they become

$\displaystyle A_{\rm Mich}$
$\displaystyle\approx\frac{8}{3}\Big{(}\frac{\bar{v}}{\sigma}\Big{)}^{6}\theta(\bar{v}^{2})e^{-(\bar{v}/\sigma)^{2}}\Big{[}\sin^{4}\Big{(}\frac{\omega
L}{2}\Big{)}+\frac{(m\bar{v}L)^{2}}{4}\Big{(}\\!\cos\omega
L-\sin^{2}\Big{(}\frac{\omega L}{2}\Big{)}\Big{)}\Big{]}+(\omega\to-\omega),$
(39) $\displaystyle B_{\rm Mich}$
$\displaystyle\approx\frac{128}{15\pi}|\bar{\omega}|K_{1}(|\bar{\omega}|)(m\sigma
L)^{2}.$ (40)

In this limit, one sees $B_{\rm Mich}\propto(m\sigma L)^{2}=(L/\lambda)^{2}$,
which signals the tidal limit of differential acceleration as estimated in the
introduction. On the other hand, for $A_{\rm Mich}$, the prefactor strongly
supports $\omega=2m$. The first term $\sin^{4}(\omega L/2)$ dominates the
second term as long as $m\sigma L\gtrsim\sigma^{2}$.

Figure 5: A cartoon for the time-delay interferometer. Each spacecraft is
loaded with two proof masses. The Michelson-like X variable measures the phase
difference of the laser along two paths shown as blue and red dashed lines.

####  III.1.2 Time-delay interferometer

Let us now consider the time-delay interferometer. We take LISA as a benchmark
for the discussion. The constellation consists of three spacecraft (SC1, SC2,
SC3) with three arms, and in each arm, two Doppler shift measurements are
carried out; in this simplified picture, there are six one-way Doppler shift
measurements. For a detailed discussion of the time-delay interferometer, we
refer readers to Refs. [35, 36].

A single one-way Doppler shift measurement can be used for the detection of
gravitational waves, but the sensitivity of such measurement will be limited
by laser frequency noise. Such laser frequency noise can be eliminated in a
certain linear combination of the Doppler shift measurements; such linear
combinations are called time-delay interferometer (TDI) variables [37, 38].

One such TDI variable is called the Michelson-like TDI variable $X$. It is
defined as

$\displaystyle X(t)$
$\displaystyle=[\Delta\phi_{3}(t)+\Delta\phi_{2}(t-2L)]-[\Delta\phi_{2}(t)+\Delta\phi_{3}(t-2L)]$
(41)

where $\Delta\phi_{2}(t)$ [$\Delta\phi_{3}(t)$] is the phase obtained by the
laser along its round trip from SC1 to SC2 (SC3) at time $t$. This TDI-X
variable measures the laser phase difference between two paths shown in Figure
5. Assuming an equal arm-length for simplicity, the $X$ in the Fourier space
is222In each spacecraft, two proof masses are loaded. We denote the position
of proof masses loaded in SCi as $\vec{x}_{i}$. We do not distinguish the
difference in the position of two proof masses loaded in the same spacecraft.
This is equivalent to assuming that the ULDM fluctuation perturbs two proof
masses in the same spacecraft in the same way. This can be justified by noting
that the wavelength of dark matter is always larger than the separation of
proof masses in the same spacecraft.

$\displaystyle\tilde{X}(\omega)=\pm\frac{\omega_{L}}{\omega^{2}}e^{i\omega
L}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{1}}\big{[}i\vec{k}\cdot\vec{D}_{X}(\hat{n}_{3},\hat{n}_{2};L)\big{]}\tilde{\Psi}.$
(42)

where

$\displaystyle\vec{D}_{X}(\hat{n}_{3},\hat{n}_{2};L)=(1-e^{2i\omega
L})[\vec{D}(\hat{n}_{3},L)-\vec{D}(\hat{n}_{2},L)]=2(1-e^{2i\omega
L})\vec{D}_{\rm Mich}(\hat{n}_{3},\hat{n}_{2};L).$ (43)

Here $\hat{n}_{2}$ ($\hat{n}_{3}$) is the unit vector from SC1 and SC2 (SC3).
The response vector $\vec{D}_{X}$ is the same as that of the Michelson
interferometer except for an overall factor of $2(1-e^{2i\omega L})$. In the
computations of the power spectrum, this additional factor leads to a common
multiplication factor $16\sin^{2}(\omega L)$, which will be ignored below.

The power spectrum due to ultralight dark matter fluctuation is

$\displaystyle S^{\rm ULDM}_{\Delta
L/L}(\omega)=\frac{\bar{a}^{2}}{\omega^{4}L^{2}}\tau\Big{[}\sigma^{4}A_{\rm
Mich}(\omega;\hat{n}_{3},\hat{n}_{2},L)+B_{\rm
Mich}(\omega;\hat{n}_{3},\hat{n}_{2},L)\Big{]}$ (44)

which is the same as the simple Michelson interferometer except that the two
arms are not orthogonal to each other, $\hat{n}_{2}\cdot\hat{n}_{3}=1/2\neq
0$. The noise power spectrum is [39]

$\displaystyle S_{\Delta
L/L}(\omega)=\frac{1}{L^{2}}\left[S_{N}(\omega)+2[1+\cos^{2}(\omega
L)]\frac{S_{\rm acc}(\omega)}{\omega^{4}}\right].$ (45)

where $S_{N}(f)$ is the optical metrology noise power spectrum, and $S_{\rm
acc}(\omega)$ is the acceleration noise power spectrum. In this expression,
the optical metrology and acceleration noise in each spacecraft are assumed to
be identical, $S_{N_{ij}}=S_{N}$ and $S_{{\rm acc},ij}=S_{\rm acc}$. The noise
spectrum and the ULDM spectrum will be compared later.

Similar to the Michelson interferometer, the response function can be obtained
analytically in the isotropic and long wavelength limit. Assuming an
equilateral configuration, we find

$\displaystyle A_{X}$
$\displaystyle\approx\frac{4}{3}\Big{(}\frac{\bar{v}}{\sigma}\Big{)}^{6}e^{-(\bar{v}/\sigma)^{2}}\theta(\bar{v}^{2})\Big{[}\sin^{4}\Big{(}\frac{\omega
L}{2}\Big{)}+\frac{3}{8}(m\bar{v}L)^{2}\cos(\omega
L)\Big{]}+(\omega\to-\omega)$ (46) $\displaystyle B_{X}$
$\displaystyle\approx\frac{32}{5\pi}|\bar{\omega}|K_{1}(|\bar{\omega}|)(m\sigma
L)^{2}$ (47)

The behavior of the response function in the time-delay interferometer is the
same as in the Michelson interferometer.

###  III.2 ULDM signal

Two distinct frequency components in the ULDM spectrum allow us to perform two
different sets of analyses for the ultralight dark matter search. In what
follows, we discuss how each frequency component of ultralight dark matter can
be searched with matched filtering and cross-correlation of multiple detector
outputs.

####  III.2.1 Deterministic signal

The $\omega=2m$ part of the spectrum represents the coherently oscillating
ULDM signals. This is a deterministic signal in the sense that the signal form
is known. For such deterministic signals, we can use the matched filter.
Consider the detector output $d(t)=s(t)+n(t)$. By convoluting the detector
output with the filter, $\int dt\,d(t)K(t)$, and choosing the filter $K(t)$
such that it maximally overlaps with the signal $s(t)$, one can add the signal
coherently. The signal-to-noise ratio is

$\displaystyle\frac{S}{N}=\left[\min(T,\sqrt{T\tau})\int_{-\infty}^{\infty}df\,\frac{S_{s}(f)}{S_{n}(f)}\right]^{1/2}$
(48)

where $T$ is the total integration time scale and $\tau=1/m\sigma^{2}$ is the
coherence time scale of the deterministic signal. Here $S_{s}(f)$ is the power
spectrum in differential arm length or strain unit from the ULDM fluctuations.
The factor $\min(T,\sqrt{T\tau})$ denotes the gain in the signal-to-noise
ratio in the coherent search, $S/N\sim\sqrt{T}$ when $T<\tau$ and the gain in
the incoherent search, $S/N\sim T^{1/4}$ when $T>\tau$. This will be used to
place an upper limit on dark matter density in the solar system. Note that the
power spectrum in this expression is defined as a two-sided spectrum.

To proceed further, we note that the power spectrum of the coherently
oscillating mode is narrowly peaked at $\omega=2m$ with a width
$\Delta\omega=m\sigma^{2}$. For the Michelson interferometry, the integral can
be approximated as

$\displaystyle\frac{S}{N}\approx\frac{\bar{a}\sigma^{2}}{2\pi^{2}f_{m}^{2}L}\bigg{[}\frac{\min(T,\sqrt{T\tau})}{S^{\rm
one}_{\Delta
L/L}(f_{m})}\bigg{]}^{\frac{1}{2}}\bigg{[}\tau\int_{0}^{\infty}df\,A_{\rm
Mich}(f)\bigg{]}^{\frac{1}{2}}$ (49)

where $f_{m}=m/\pi$ and $S^{\rm one}_{\Delta L/L}(f)=2S_{\Delta L/L}(f)$ is
the one-sided noise power spectrum. We have evaluated all the other factors in
the integral at $f=f_{m}$ except for the exponential factor strongly peaked at
$\omega=2m$, assuming that they are smoothly distributed over
$\Delta\omega\simeq m\sigma^{2}$. In the isotropic ($v_{0}/\sigma\ll 1$) and
long-wavelength limit $(m\sigma L\ll 1)$, one finds
$\tau\int_{0}^{\infty}df\,{\cal A}_{\rm
Mich}\simeq(8/\pi)\sin^{4}(f_{m}/2f_{*})$ with $f_{*}=(2\pi L)^{-1}$. The
above expression simplified to

$\displaystyle\frac{S}{N}=\frac{2^{5/2}}{\sqrt{\pi}}\sin^{2}(f_{m}/2f_{*})\left[\frac{\bar{a}\sigma^{2}}{(2\pi
f_{m})^{2}L}\right]\bigg{[}\frac{\min(T,\sqrt{T\tau})}{S^{\rm one}_{\Delta
L/L}(f_{m})}\bigg{]}^{\frac{1}{2}}.$ (50)

The quantity in the first squared parenthesis is the typical displacement
$\Delta L/L$ due to ULDM fluctuation during $\Delta t\sim 1/f_{m}$. Note that
$\bar{a}=Gm_{\rm eff}/\lambda^{2}$ and $m_{\rm eff}=\pi^{3/2}\rho\lambda^{3}$.
The above equation can be solved for $\rho/\rho_{0}$ with the local dark
matter density $\rho_{0}=0.4\,{\rm GeV/cm^{3}}$. For the interferometers with
equilateral geometry, $\tau\int_{0}^{\infty}A_{\rm
Mich}\approx(4/\pi)\sin^{4}(f_{m}/2f_{*})$. The additional factor of $1/2$ is
because the two arms are not orthogonal to each other.

####  III.2.2 Stochastic signal

The low-frequency part of the ULDM spectrum behaves similarly to the
stochastic gravitational wave background. To search for this part of the
spectrum, we ideally need more than two detectors, which will allow us to
cross-correlate the outputs in different detectors. If the separation of
detectors is smaller than the wavelength of dark matter, the output at each
detector is expected to be correlated, while the noise is expected to be
uncorrelated. In this way, we can single out the ultralight dark matter signal
in the output data stream of multiple detectors.

Consider two detector output $d_{1,2}(t)=s_{1,2}(t)+n_{1,2}(t)$ with the
signal $s_{1,2}$ and the noise $n_{1,2}$. Then define the cross-correlated
output $Y$ as

$\displaystyle
Y=\int_{-T/2}^{T/2}dt\int_{-T/2}^{T/2}dt^{\prime}\,d_{1}(t)d_{2}(t^{\prime})Q(t-t^{\prime})\simeq\int_{-\infty}^{\infty}df\,\tilde{d}_{1}(f)\tilde{d}_{2}^{*}(f)\tilde{Q}^{*}(f)$
(51)

where $Q(t)$ is a real filter function. Suppose that the cross-correlation of
signals at two detectors is
$\langle\tilde{s}_{1}(f)\tilde{s}_{2}^{*}(f^{\prime})\rangle=\delta(f-f^{\prime})S_{12}(f)$.
Using the optimal filter $\tilde{Q}(f)=S_{12}(f)/S_{n}^{2}(f)$ and assuming
uncorrelated detector noise, one finds the signal-to-noise ratio as [40]

$\displaystyle\frac{S}{N}=\left[T\int_{-\infty}^{\infty}df\,\frac{|S_{12}(f)|^{2}}{S_{n}^{2}(f)}\right]^{1/2},$
(52)

where $T$ is the total integration time scale. In this way, the low-frequency
part of the ultralight dark matter signal can be searched.

For cross-correlation, we are mostly interested in space-borne interferometers
for the following reasons. For the cross-correlation to be relevant, two
conditions need to be met: (i) the ultralight dark matter-induced signal must
lie within the detector bandwidth, and (ii) the ultralight dark matter signals
must be correlated in two detectors. Each interferometer has a range of
frequencies to which the detector is sensitive. For the LIGO, the minimum
frequency that LIGO is sensitive to is about $f_{\rm min}\sim 10\,{\rm Hz}$.
On the other hand, two LIGO detectors, one at Hanford and the other at
Livingston, are separated by $L_{\rm det}\simeq 3000\,{\rm km}$. To guarantee
that the ultralight dark matter signal is correlated over these two detectors,
we require $L_{\rm det}<(1/\Delta p)=(\sigma/\omega)$. This results in
$\omega<5\times 10^{-2}\,{\rm Hz}$, which is more than two orders of magnitude
smaller than the minimum detectable frequency in LIGO detectors. In other
words, for the frequency range where the ultralight dark matter signals would
be correlated over two LIGO detectors, LIGO detectors are not sensitive to
ULDM signals due to a large seismic noise.

The situation in the space-borne interferometer is slightly different. For the
sake of discussion, let us consider two co-planar LISA-like constellations.
The minimum frequency of the LISA detector is $f_{\rm min}\sim 10^{-5}\,{\rm
Hz}$. On the other hand, the detector separation is $L_{\rm det}\sim
L=2.5\times 10^{6}\,{\rm km}$ or less. This translates into the condition
$\omega<\sigma/L_{\rm det}\simeq 10^{-4}\,{\rm Hz}$ for ULDM signals to be
correlated over two detectors. There is a range of frequencies where the
signals are correlated over detectors, while the detector maintains its
sensitivity to the ultralight dark matter signals. The same arguments hold for
other types of proposed space-borne interferometers, such as $\mu$Ares [25], a
proposal using asteroids as GW detectors [26], and Big Bang Observatory (BBO)
[41].

For this reason, we only consider space-borne interferometers. The cross-
correlation of TDI-X variables is

$\displaystyle S_{12}(f)=\frac{\bar{a}^{2}\tau}{(2\pi f)^{4}L^{2}}{\cal
B}_{\rm cross}.$ (53)

The overlap function is

$\displaystyle{\cal B}_{\rm cross}=$ $\displaystyle\cos^{2}(2\pi fL){\cal
B}^{11^{\prime}}(\vec{x}_{11^{\prime}})+{\cal
B}^{33^{\prime}}(\vec{x}_{33^{\prime}})+{\cal
B}^{22^{\prime}}(\vec{x}_{22^{\prime}})-{\cal
B}^{23^{\prime}}(\vec{x}_{23^{\prime}})-{\cal
B}^{32^{\prime}}(\vec{x}_{32^{\prime}})$ $\displaystyle-\cos(2\pi
fL)\Big{[}{\cal B}^{31^{\prime}}(\vec{x}_{31^{\prime}})+{\cal
B}^{13^{\prime}}(\vec{x}_{13^{\prime}})-{\cal
B}^{21^{\prime}}(\vec{x}_{21^{\prime}})-{\cal
B}^{12^{\prime}}(\vec{x}_{12^{\prime}})\Big{]},$ (54)

where $\vec{x}_{i}$ is the position of $i$-th spacecraft in the first
constellation, $\hat{n}_{2,3}=(\vec{x}_{2,3}-\vec{x}_{1})/L$,
$\hat{n}_{1}=(\vec{x}_{3}-\vec{x}_{2})/L$,
$\vec{x}_{ij^{\prime}}\equiv\vec{x}_{i}-\vec{x}_{j^{\prime}}$, and ${\cal
B}^{ab^{\prime}}=\hat{n}^{i}_{a}\hat{n}^{j}_{b^{\prime}}{\cal B}^{ij}$. The
primed quantities are for the second constellation. We take $L=L^{\prime}$ for
simplicity. The signal-to-noise ratio is

$\displaystyle\frac{S}{N}=\frac{\bar{a}^{2}}{(2\pi
f_{0})^{4}L^{2}}\left[\frac{T}{(2\pi)^{2}}\int_{-\infty}^{\infty}df\,\Big{(}\frac{f_{0}}{f}\Big{)}^{8}\frac{|{\cal
B}_{\rm cross}(f)|^{2}}{f_{0}^{2}S_{\Delta L/L}^{2}(f)}\right]^{1/2}.$ (55)

where $f_{0}=\omega_{0}/2\pi=m\sigma^{2}/2\pi$. The prefactor is the squared
value of the typical differential length change, $(\Delta L/L)^{2}$, over
$\Delta t\sim 1/f_{0}$. The integral will be numerically solved.

###  III.3 Application

We apply the discussion in previous sections to the current and proposed
gravitational wave interferometers.

####  III.3.1 LIGO

We take the sensitivity curve of LIGO O3 from Ref. [42]. The strain noise
spectrum from ultralight dark matter can be found by using $S_{\Delta
L/L}(f)={\rm sinc}^{2}(2\pi fL)S_{h}(f)$. For the LIGO with $f\lesssim
10^{4}\,{\rm Hz}$, ${\rm sinc}(2\pi fL)\simeq 1$, we have ignored the
additional ${\rm sinc}^{2}(2\pi fL)$ factor. The arm length is $L_{\rm
LIGO}=4\,{\rm km}$. As we have already estimated in the introduction, the
ULDM-induced noise is by many orders of magnitude smaller than the
instrumental noise. Still, it can provide an upper limit on the dark matter
density in the solar system. For the search with a matched filter, we use (50)
and rewrite it as

$\displaystyle\frac{\rho}{\rho_{0}}\approx\frac{S}{N}\sqrt{\frac{\pi}{2}}\frac{1}{{\rm
sinc}^{2}(f_{m}/2f_{*})}\left[\frac{(2\pi
f_{*})^{2}L}{\bar{a}_{0}\sigma^{2}}\right]\Big{[}\frac{S_{\Delta L/L}^{\rm
one}(f_{m})}{\min(T,\sqrt{T\tau})}\Big{]}^{\frac{1}{2}}$ (56)

where $\bar{a}_{0}=Gm_{\rm eff,0}/\lambda^{2}$, $m_{\rm
eff,0}=\pi^{3/2}\rho_{0}\lambda^{3}$ with $\rho_{0}=0.4\,{\rm GeV/cm^{3}}$,
and $f_{*}=(2\pi L)^{-1}=10^{4}\,{\rm Hz}$ for LIGO. Note that for LIGO, we
can also approximate ${\rm sinc}^{2}(f_{m}/2f_{*})\simeq 1$. For $T=1\,{\rm
yr}$, $\sigma=164\,{\rm km/sec}$, and $S/N=1$, the upper limit on the dark
matter density $\rho/\rho_{0}$ is obtained and shown as a solid red line in
Figure 2. Note that the above expression is valid in the long-wavelength limit
$m\sigma L\ll 1$, which is satisfied not only in LIGO but also in all other
ground-based interferometers.

####  III.3.2 Einstein Telescope

Einstein telescope is a proposed future gravitational wave detector based on a
dual-recycled Michelson interferometer. The arm length is $L_{\rm ET}=10\,{\rm
km}$. We assume an equilateral geometry for the Einstein telescope. The strain
noise power spectrum is taken from [43]. Using the same expression (56) with
$f_{*}=(2\pi L_{\rm ET})^{-1}=5\times 10^{3}\,{\rm Hz}$, we show the upper
limit on the density, $\rho/\rho_{0}$, from Einstein Telescope with the blue
solid line in Figure 2. We choose the same parameters: $T=1\,{\rm yr}$,
$\sigma=164\,{\rm km/sec}$, and $S/N=1$.

####  III.3.3 LISA

The noise source in LISA is conveniently divided into the collective optical
metrology noise $S_{N}(f)$ and the single test-mass acceleration noise $S_{\rm
acc}(f)$. Each of them is given by [44]

$\displaystyle S_{N}^{1/2}(f)$ $\displaystyle=10^{-11}\frac{\rm m}{\sqrt{\rm
Hz}}\bigg{[}1+\left(\frac{2\rm mHz}{f}\right)^{4}\bigg{]}^{1/2},$ (57)
$\displaystyle S_{\rm acc}^{1/2}(f)$ $\displaystyle=3\times
10^{-15}\,\frac{\rm m\,s^{-2}}{\sqrt{\rm Hz}}\left[1+\left(\frac{0.4\,{\rm
mHz}}{f}\right)^{2}\right]^{\frac{1}{2}}\left[1+\left(\frac{f}{8\,{\rm
mHz}}\right)^{4}\right]^{\frac{1}{2}}.$ (58)

LISA Pathfinder has reported differential acceleration noise of two test
masses loaded in a single spacecraft, which is smaller than the LISA proposal
requirement by a factor of few [45, 46]. For the matched filter search, we
reuse (56) with an additional factor of $\sqrt{2}$, arising from the fact that
LISA has an equilateral geometry and two arms are not orthogonal to each
other. When recycling (56), we replace ${\rm sinc}^{2}(f_{m}/2f_{*})\to
1/[1+(f_{m}/2f_{*})^{2}]$ to show the envelope of the projection. The upper
limit on the dark matter density $\rho/\rho_{0}$ with the same parameters is
shown as a green solid line in Figure 2.

For the LISA, we also compare $S_{\Delta L/L}^{\rm ULDM}(f)$ with the mission
requirement $S_{\Delta L/L}(f)$. For the comparison, we compare two spectra in
the GW strain unit. To convert to the strain noise power spectrum, we divided
the power spectrum for the differential arm-length with the sky and
polarization averaged detector response function ${\cal R}(f)$ [47],

${\cal R}(f)=\frac{3}{10(1+0.6(f/f_{*})^{2})}$

with $f_{*}=(2\pi L_{\rm LISA})^{-1}\simeq 19\,{\rm mHz}$. In the left panel
of Figure 1, we show the mission requirement (red solid line) and $S_{n}^{\rm
ULDM}$ with $m=10^{-12}\,{\rm eV}$ (cyan) and $m=10^{-13}\,{\rm eV}$ (orange).
The ULDM-induced noise is several orders of magnitude smaller than the mission
requirement, and therefore, it can safely be ignored for the purpose of
gravitational wave detection.

The subdominant low-frequency ULDM signals might be probed if there are more
than two detectors. Although the current LISA mission concept does not invoke
more than two constellations [44], we assume two LISA-like interferometers and
investigate the reach of such configuration for the ULDM search as an
exercise. Note that other space-borne interferometers, such as Taiji and
TianQin, would allow the cross-correlation searches [30]. For the cross-
correlation, the position of the spacecraft needs to be specified. In the
heliocentric ecliptic coordinate, their position is [48]

$\displaystyle x$
$\displaystyle=a\cos\alpha+ae[\sin\alpha\cos\alpha\sin\beta-(1+\sin^{2}\alpha)\cos\beta]$
(59) $\displaystyle y$
$\displaystyle=a\sin\alpha+ae[\sin\alpha\cos\alpha\cos\beta-(1+\cos^{2}\alpha)\sin\beta]$
(60) $\displaystyle z$ $\displaystyle=\sqrt{3}ae\cos(\alpha-\beta)$ (61)

where $a\approx{\rm AU}$ is the semi-major axis, $e$ is eccentricity,
$\alpha=\omega t+\kappa$, and $\beta=2n\pi/3+\lambda$ with $n=0,1,2$. Here
$\alpha$ is the orbital phase of the guiding center. Note that the arm length
is given by $L=2\sqrt{3}\epsilon R$; for $R={\rm AU}$ and $L_{\rm
LISA}=2.5\times 10^{6}\,{\rm km}$, the eccentricity is $\epsilon=4.8\times
10^{-3}$. The position of spacecraft in the second constellation is given by
the same expression but with a prime ($\kappa^{\prime}$, $\lambda^{\prime}$,
and so on).

Rewriting (55) in terms of $\rho/\rho_{0}$, we find

$\displaystyle\frac{\rho}{\rho_{0}}=\frac{1}{2^{3/4}}\Big{(}\frac{S}{N}\Big{)}^{\frac{1}{2}}\left[\frac{(2\pi
f_{0})^{2}L}{\bar{a}_{0}}\right]\left[\frac{T}{(2\pi)^{2}}\int_{0}^{\infty}df\,\Big{(}\frac{f_{0}}{f}\Big{)}^{8}\frac{|{\cal
B}_{\rm cross}|^{2}}{f_{0}^{2}|S_{\Delta L/L}^{\rm
one}(f)|^{2}}\right]^{-\frac{1}{4}}.$ (62)

Using the halo dark matter parameters, $\sigma=164\,{\rm km/sec}$ and
$\rho_{0}=0.4\,{\rm GeV/cm^{3}}$, and assuming
$\Delta\kappa=\kappa^{\prime}-\kappa=0$ (co-planar configuration),
$\Delta\lambda=\lambda^{\prime}-\lambda=\pi/2$, and $T=5\,{\rm yr}$, we find
the reach of two LISA-like detectors as a dashed green line in Figure 2. For
the frequency range, we use $f_{\rm min}=10^{-5}\,{\rm Hz}$ and $f_{\rm
max}=1\,{\rm Hz}$. We have included a galactic confusion noise, assuming
$T=4\,{\rm yr}$ of operation [49].

####  III.3.4 $\mu$Ares

We also consider the $\mu$Ares proposal [25] for the ultralight dark matter
search. The proposed concept is similar to LISA in the sense that it consists
of three spacecraft, and it operates as a time-delay interferometer. The
configuration of satellites is heliocentric, while the arm-length is
$L=400{\rm M}\,{\rm km}$. For the matched filter and cross-correlation
analysis, we take the strain noise power spectrum in [25], including
astrophysical foreground. The upper limit on the dark matter density
$\rho/\rho_{0}$ with matched filtering is shown as a purple line in Figure 2.

We also compare the ultralight dark matter-induced noise with the required
noise level for the proposed mission concept. The proposal requires $S_{\rm
acc}^{1/2}(f)=10^{-15}\,{\rm m\,s^{-2}/\sqrt{Hz}}$, and
$S^{1/2}_{N}(f)=50\,{\rm pm}/\sqrt{\rm Hz}$. The result is shown on the right
panel of Figure 1 for $m=10^{-12}\textrm{--}10^{-15}\,{\rm eV}$. Note that the
noise curve in this figure is obtained by using (45) with the proposal
requirements, but without any astrophysical foregrounds. Note also that, for
$f\lesssim\textrm{a few}\times 10^{-7}\,{\rm Hz}$, the gravity gradient noise
from asteroids could dominate the total noise spectrum [27].

For the cross-correlation analysis, we specify the position of each spacecraft
in the heliocentric ecliptic coordinate as

$\displaystyle\vec{x}_{n}=R\Big{(}\\!\cos[\phi+2\pi(n-1)/3],\,\sin[\phi+2\pi(n-1)/3],\,0\Big{)}$
(63)

where $n=1,2,3$ and $\phi$ is the orbital phase of the interferometer. We
assume another of such constellations lying on the same plane (co-planar
configuration), with the relative phase with respect to the first
constellation by $\Delta\phi$. Using (62), we find the projected sensitivity
of $\rho/\rho_{0}$ of the cross-correlation search in $\mu$Ares as a purple
dashed line in Figure 2. We use $f_{\rm min}=4\times 10^{-7}\,{\rm Hz}$,
$f_{\rm max}=10^{-1}\,{\rm Hz}$, and $\Delta\phi=\pi/2$.

####  III.3.5 Big Bang Observatory

The configuration of the Big Bang Observatory is similar to the LISA, except
for a smaller arm-length, $L_{\rm BBO}=5\times 10^{7}\,{\rm m}$. The
instrument parameters are [50]

$\displaystyle S_{N}^{1/2}$ $\displaystyle=$ $\displaystyle 1.4\times
10^{-17}\,{\rm m/\sqrt{Hz}},$ (64) $\displaystyle S_{\rm acc}^{1/2}$
$\displaystyle=$ $\displaystyle 3\times 10^{-17}\,{\rm m\,s^{-2}/\sqrt{Hz}}.$
(65)

The cross-correlation search is identical to the one in LISA, except for this
arm-length difference and different instrumental parameters. The matched
filter and cross-correlation search (using two co-planar constellations)
yields the upper limit on $\rho/\rho_{0}$, represented by the (dashed) brown
lines in Figure 2. For the analysis, we use $f_{\rm min}=10^{-4}\,{\rm Hz}$,
$f_{\rm max}=10^{2}\,{\rm Hz}$, and $\Delta\lambda=\pi$. We have included the
same galactic confusion noise as in the case of LISA while assuming the other
astrophysical foreground around $\sim{\cal O}(0.1)\,{\rm Hz}$ from neutron
star binaries can be fully resolved [51].

##  IV Discussion

We discuss now some of the ignored effects in the main discussion and some of
the other implications of our results.

In the main discussion, we have ignored the mean acceleration on the test
masses due to ULDM; we have considered only its variance. The mean
acceleration is nothing but the dynamical friction, and it is present as long
as $\vec{v}_{0}\neq 0$ in the rest frame of the detector. One may directly
compute the dynamical friction from the expression (19). The dynamical
friction for the ultralight dark matter has already been obtained in [52],
where it is given by

$\displaystyle\langle a^{i}\rangle=\hat{v}_{0}^{i}\frac{8\pi G^{2}\rho m_{\rm
eff}\ln\Lambda}{\sigma^{2}}\left[G\Big{(}\frac{v_{0}}{\sigma}\Big{)}+\frac{m_{t}}{2m_{\rm
eff}}G\Big{(}\frac{v_{0}}{\sqrt{2}\sigma}\Big{)}\right]$ (66)

Here $G(X)=1/(2X^{2})[{\rm erf}(X)-(2X/\sqrt{\pi})e^{-X^{2}}]$, $m_{t}$ is the
mass of the test body, and $\ln\Lambda$ is the Coulomb logarithm factor
arising due to the long-range nature of gravitational force. For the halo dark
matter with $\rho_{0}=0.4\,{\rm GeV/cm^{3}}$ and $m=10^{-22}\,{\rm eV}$, this
dynamical friction force is causing $\Delta x\sim{\cal O}(1)\,{\rm meter}$ of
drift of constellation over a year, which is much smaller than the expected
unavoidable change of arm-length of LISA, $\Delta x\sim{\cal
O}(10^{4}\textrm{--}10^{5})\,{\rm km}$, due to the orbital dynamics [44]. In
addition, the dynamical friction force acts on each spacecraft in the same
way, and therefore, we do not expect that it changes the arm length itself.
For these reasons, we do not expect that the dynamical friction is relevant
for the discussion of the ULDM-induced signals in the interferometers.

We have found that some of the future gravitational wave interferometers could
potentially probe $\rho/\rho_{0}\sim{\cal O}(10^{2})$ where
$\rho_{0}=0.4\,{\rm GeV/cm^{3}}$ is the local dark matter density. Since we
have projected the sensitivity of interferometers in units of $\rho_{0}$ and
this local dark matter density is one of the key input parameters for all
terrestrial dark matter detectors, it is worthwhile to discuss how this value
$\rho_{0}\sim{\cal O}(0.1)\,{\rm GeV/cm^{3}}$ is obtained. The current
measurement of local dark matter is usually performed over a much larger
spatial scale. For instance, even the very local measurements of the dark
matter density select stellar tracers within ${\cal O}(10^{2})$ pc around the
solar system, and therefore, the inferred value of dark matter density is an
average value over that spatial volume, $V\sim[{\cal O}(10^{2})\,{\rm
pc}]^{3}$ [28, 29]. When it comes to the measurement within the solar system,
the dark matter density is poorly constrained. Measurement using the planetary
ephemerides places an upper limit on the dark matter density in the solar
system as $\rho/\rho_{0}\lesssim 10^{4}$ [53], while the measurements of lunar
laser ranging and LAGEOS geodetic satellite place an upper limit on the dark
matter density between the Moon and the satellite’s orbital radius as
$\rho/\rho_{0}\lesssim 10^{11}$ [54].

Furthermore, a large density fluctuation might just arise statistically in the
ultralight dark matter halo. For the density of the field, $\rho_{\phi}$, its
underlying distribution is given by the exponential distribution,
$p(\rho_{\phi})d\rho_{\phi}=[\exp(-\rho_{\phi}/\rho_{0})/\rho_{0}]d\rho_{\phi}$
(see Appendix C). Given that the correlation of density field exponentially
drops over $L\gtrsim 1/m\sigma$ [52], we may consider the density fluctuation
$\rho_{\phi}$ over wavelength-sized patches as statistically independent
random variables. The probability of finding a patch with
$\rho_{\phi}>c\rho_{0}$ is then given by

$P(\rho_{\phi}>c\rho_{0})=\int_{c\rho_{0}}^{\infty}p(\rho_{\phi})=e^{-c}.$

On the other hand, the number of statistically independent patches over the
volume $V$ is $N\sim V/\lambda^{3}$. For $m=10^{-15}\,{\rm eV}$, the
wavelength is about AU scale, and the number of patches over $V=(10\,{\rm
kpc})^{3}$ volume is $N=V/\lambda^{3}=(10\,{\rm kpc}/{\rm AU})^{3}\sim
10^{28}$. For $c=60$, one finds $P(\rho_{\phi}>60\rho_{0})=e^{-60}\simeq
10^{-26}$. This would mean that there might be several hundreds of AU-sized
patches within $10\,{\rm kpc}$ volume with $\rho\gtrsim 60\rho_{0}$.

We have only considered cross-correlation searches for space-borne
interferometers. The reason that it was not suitable for the ground-based
interferometers, especially for the current advanced LIGO, was that the
detector separation is too large for the ultralight dark matter signals to be
correlated for the relevant frequency range of the detector. However, in the
initial LIGO phase, two detectors at Hanford, with $4\,$km and $2\,$km arms,
were co-located and co-aligned, and therefore, it is in principle possible to
use the result of two Hanford detectors in the initial LIGO phase for the
cross-correlation search as demonstrated in [55, 56] for the stochastic
gravitational waves. Since the detectors are co-located, environmental noises
are correlated, and this correlated noise must be carefully accounted for in
the analysis.

For the low-frequency ULDM signals, we have considered the cross-correlation,
which requires more than two detectors. One may still be able to search for
these low-frequency ULDM signals with a single detector with the help of a
null data stream. A null data stream is the output of the detector where the
expected signal, either GWs or ULDM, is expected to be suppressed, allowing a
way to calibrate the instrumental noises. For the time-delay interferometers,
like LISA, a symmetrized Sagnac channel is shown to be insensitive to the
gravitational wave signals at low frequencies, $f<f_{*}=(2\pi L)^{-1}$, which
could, at least partially, provide a way to distinguish for stochastic
gravitational wave backgrounds from other noise sources at such low
frequencies [57, 58, 59]. Similarly to the stochastic gravitational wave
backgrounds, ultralight dark matter signal is strongly suppressed in the
symmetrized Sagnac channel as $\propto(m\sigma L)^{4}$ for $m\sigma L<1$, and
therefore, using such null data stream might provide a way to search for ULDM
at sufficiently small masses with the help of additional information such as
annual modulations of the ULDM signal.

##  V Conclusion

We have investigated how the density fluctuations of ultralight dark matter
interact with interferometers, designed for the detection of gravity waves. We
provide a systematic way to compute the ULDM spectrum in any gravitational
wave interferometers. We show that the ultralight dark matter-induced noise is
most significant when the arm-length is large, such as LISA or other proposals
with an astronomical size of arm-length, and that, in all cases, we consider,
the ultralight dark matter effects are subdominant compared to other noise
sources, e.g. interferometric read-out noise and acceleration noise. Then we
consider if such interferometers can be used to place an upper bound on the
dark matter density in the solar system, and find that, under certain
assumptions, $\rho/\rho_{0}\sim{\cal O}(10^{2})$ might be probed with future
gravitational wave interferometers with its arm-length close to astronomical
unit.

We note that the current local dark matter density measurement,
$\rho_{0}=0.4\,{\rm GeV/cm^{3}}$, is performed on a much larger scale, e.g.
$\gtrsim{\cal O}(100)\,{\rm pc}$, and that the constraints on dark matter
abundance in the solar system from planetary ephemerides and laser ranging
experiments are several orders of magnitude larger than $\rho_{0}$. Especially
in the light of a situation where there is no direct and gravitational
detection of dark matter in the solar system, these interferometric approaches
will provide an interesting way to probe ultralight dark matter only through
gravitational interaction. In addition, it has been recently proposed that the
ultralight dark matter density in the solar system might be larger than the
local dark matter density through certain capture processes [60, 61, 62]. Our
interferometric search for ultralight dark matter in the solar system is
expected to provide an interesting probe for such possibilities.

Our analysis can be straightforwardly generalized to pulsar timing array (PTA)
observations. Coherent ULDM signals in PTA observations were already explored
by Khmelnitsky and Rubakov [63], whose results are reproduced in Appendix D
with the formalism discussed in this work. Low-frequency stochastic
fluctuations of ULDM in the context of pulsar timing measurements, however,
have not been explored. As the distance between the Earth and pulsars is on a
kiloparsec scale, and the pulsar timing analyses involve a cross-correlation
of timing residuals of pulsars with different sky locations, PTAs offer
another interesting probe for stochastic ULDM density fluctuations. In future
work, we will investigate how these stochastic ULDM signals affect pulsar
timing measurements and if one can probe ULDM density within the solar system
with PTA observations [64].

###### Acknowledgements.

We would like to thank Julián Rey and Nicholas Rodd for useful conversations.
We especially thank Alessandro Lenoci for the initial collaboration. We also
thank Thomas Konstandin, Germano Nardini, and Mauro Pieroni for their useful
comments and suggestions on the manuscript. This work is supported by the
Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy - EXC 2121
Quantum Universe - 390833306.

##  Appendix A Power spectrum

We provide a detailed computation of the low-frequency part of the power
spectrum $S_{\nabla\Psi}^{ij}(\omega,\vec{L})$. To prepare an explicit
computation, we derive the integral expression of the power spectrum (25) in
the non-relativistic limit:

$\displaystyle S_{\nabla\Psi}^{ij}(\omega,\vec{L})=$ $\displaystyle
32\pi^{3}G^{2}\int d^{3}v_{c}\int
d^{3}v_{d}\,f(\vec{v}_{c}+\vec{v}_{d})f(\vec{v}_{c}-\vec{v}_{d})$ (67)
$\displaystyle\times$
$\displaystyle\Big{[}v_{c}^{i}v_{c}^{j}\Big{(}\delta(\omega-\omega_{1}-\omega_{2})e^{2im\vec{v}_{c}\cdot\vec{L}}+\delta(\omega+\omega_{1}+\omega_{2})e^{-2im\vec{v}_{c}\cdot\vec{L}}\Big{)}$
$\displaystyle+\frac{v_{d}^{i}v_{d}^{j}}{v_{d}^{4}}\Big{(}\delta(\omega-\omega_{1}+\omega_{2})e^{2im\vec{v}_{d}\cdot\vec{L}}+\delta(\omega+\omega_{1}-\omega_{2})e^{-2im\vec{v}_{d}\cdot\vec{L}}\Big{)}\Big{]}$

where we have introduced new integration variables $\vec{v}_{c}$ and
$\vec{v}_{d}$, defined as

$\vec{v}_{c}=\frac{1}{2}(\vec{v}_{1}+\vec{v}_{2}),\qquad\vec{v}_{d}=\frac{1}{2}(\vec{v}_{1}-\vec{v}_{2}).$

The two terms in the first line are related by the complex conjugate with
$\omega\to-\omega$. The two terms in the last line are the same under the
exchange of the integration variable $\vec{v}_{1}\leftrightarrow\vec{v}_{2}$.
The first two terms represent the spectrum at $\omega=\pm 2m$, denoted by
$A^{ij}(\omega,\vec{L})$ in the main text, and the last two terms represent
the low-frequency part, denoted by $B^{ij}(\omega,\vec{L})$. Since the
computation of $\omega=\pm 2m$ mode is straightforward, only the computation
of $B^{ij}$ will be given below.

The relevant part of the power spectrum is

$\displaystyle S_{\nabla\Psi}^{ij}(\omega,\vec{L})=$
$\displaystyle\frac{32\pi^{3}G^{2}}{m}\int d^{3}v_{c}\int
d^{3}v_{d}\,f(\vec{v}_{c}+\vec{v}_{d})f(\vec{v}_{c}-\vec{v}_{d})\frac{v_{d}^{i}v_{d}^{j}}{v_{d}^{4}}\delta\Big{(}\vec{v}_{c}\cdot\vec{v}_{d}-\frac{\omega}{2m}\Big{)}e^{2im\vec{v}_{d}\cdot\vec{L}}$
(68)

where we have used the symmetry of integrand under the exchange of integration
variable $v_{1}\leftrightarrow v_{2}$ in (67) and take the non-relativistic
limit for the $\delta$-function.333In the case of $\vec{L}=0$ and $\omega\to
0$ the above power spectrum is nothing but the second diffusion coefficient
$D[\Delta v^{i}(T)\Delta v^{j}(T)]=\langle\Delta v^{i}(T)\Delta
v^{j}(T)\rangle/T$ of the test mass with the change of velocity $\Delta
v^{i}(T)$ over the time $T$. Compare the expression with (54) in Bar-Or et al
[52]. Using the normal distribution

$f(\vec{v})=\frac{(\rho_{0}/m)}{(2\pi\sigma^{2})^{3}}\exp\left[-\frac{(\vec{v}-\vec{v}_{0})^{2}}{2\sigma^{2}}\right],$

we find

$\displaystyle S_{\nabla\Psi}^{ij}(\omega,\vec{L})=$
$\displaystyle\frac{4(G\rho_{0})^{2}}{m^{3}\sigma^{6}}\int d^{3}v_{c}\int
d^{3}v_{d}\,\exp\left[-\frac{(\vec{v}_{c}-\vec{v}_{0})^{2}}{\sigma^{2}}-\frac{v_{d}^{2}}{\sigma^{2}}\right]\frac{v_{d}^{i}v_{d}^{j}}{v_{d}^{4}}\delta\Big{(}\vec{v}_{c}\cdot\vec{v}_{d}-\frac{\omega}{2m}\Big{)}e^{2im\vec{v}_{d}\cdot\vec{L}}$
$\displaystyle=$ $\displaystyle\frac{4(G\rho_{0})^{2}}{m^{3}\sigma^{4}}\int
d^{3}v_{c}\int
d^{3}v_{d}\,e^{-(\vec{v}_{c}-\vec{v}_{0}/\sigma)^{2}}e^{-v_{d}^{2}}\frac{v_{d}^{i}v_{d}^{j}}{v_{d}^{4}}\delta\Big{(}\vec{v}_{c}\cdot\vec{v}_{d}-\frac{\bar{\omega}}{2}\Big{)}e^{2i\vec{v}_{d}\cdot\vec{L}_{\lambda}}$
$\displaystyle=$
$\displaystyle-\frac{(G\rho_{0})^{2}}{m^{3}\sigma^{4}}\frac{\partial^{2}}{\partial
L_{\lambda}^{i}\partial L_{\lambda}^{j}}\int d^{3}v_{c}\int
d^{3}v_{d}\,e^{-(\vec{v}_{c}-\vec{v}_{0}/\sigma)^{2}}e^{-v_{d}^{2}}\frac{1}{v_{d}^{4}}\delta\Big{(}\vec{v}_{c}\cdot\vec{v}_{d}-\frac{\bar{\omega}}{2}\Big{)}e^{2i\vec{v}_{d}\cdot\vec{L}_{\lambda}}$
(69)

where in the second line we have changed the integration variables
$\vec{v}_{c,d}\to\vec{v}_{c,d}\sigma$, and defined
$\vec{L}_{\lambda}=m\sigma\vec{L}$. In the last line, the vector $v_{d}^{i}$
in the integrand is replaced with the derivative with respect to
$\vec{L}_{\lambda}$.

To proceed further, we use the integral representation of $\delta$-function:
$\delta(x)=\int(ds/2\pi)e^{-ixs}$. We then perform the $\vec{v}_{c}$ integral,
and then the $\hat{v}_{d}$ integral. After some manipulation, we find

$\displaystyle S_{\nabla\Psi}^{ij}(\omega,\vec{L})=$
$\displaystyle-4\pi^{3/2}\frac{(G\rho_{0})^{2}}{m^{3}\sigma^{4}}\frac{\partial^{2}}{\partial
L_{\lambda}^{i}\partial
L_{\lambda}^{j}}\int_{-\infty}^{\infty}ds\,e^{-is\bar{\omega}}\sqrt{1+s^{2}}\int_{0}^{\infty}\frac{e^{-y^{2}}}{y^{2}}{\rm
sinc}(2yY)$ (70)

where $\vec{Y}=(1+s^{2})^{-1/2}(s\vec{v}_{0}/\sigma+\vec{L}_{\lambda})$. The
$y$ integral is divergent. The divergence can be regulated by ${\rm
sinc}(2yY)\to{\rm sinc}(2yY)-1$. The additional $-1$ term would not contribute
to the final result as it does not depend on $L_{\lambda}$. Using this
regularization, we obtain

$\displaystyle S_{\nabla\Psi}^{ij}(\omega,\vec{L})=$
$\displaystyle(\bar{a}^{2}\tau)\left[\frac{2}{\sqrt{\pi}}\int_{-\infty}^{\infty}\frac{e^{-is\bar{\omega}}}{\sqrt{1+s^{2}}}\left(\delta^{ij}\frac{{\rm
erf}(Y)-G(Y)}{Y}+Y^{i}Y^{i}\frac{3G(Y)-{\rm erf}(Y)}{Y^{3}}\right)\right]$
$\displaystyle\equiv$ $\displaystyle(\bar{a}^{2}\tau)B^{ij}(\omega,\vec{L}).$
(71)

This reproduces the result used in the main text.

##  Appendix B Induced phase

In this Appendix, we discuss the approximation for the phase (34). We begin
with

$\displaystyle\Delta\phi(t)=2\omega_{L}(\delta t+\delta L)$ (72)

where the time delay $\delta t$ and the arm-length fluctuation $\delta L$ are

$\displaystyle\delta L$
$\displaystyle=+\frac{1}{2}\hat{x}\cdot\Big{[}\big{(}\delta\vec{x}_{1}(t-L)-\delta\vec{x}_{0}(t-2L)\big{)}+\big{(}\delta\vec{x}_{1}(t-L)-\delta\vec{x}_{0}(t)\big{)}\Big{]},$
(73) $\displaystyle\delta t$
$\displaystyle=-\frac{1}{2}\left[\int^{t}_{t-L}dt^{\prime}+\int^{t-L}_{t-2L}dt^{\prime}\right][\Psi(t^{\prime},\vec{x}(t^{\prime}))+\Phi(t^{\prime},\vec{x}(t^{\prime}))].$
(74)

In the Fourier space, each of them is

$\displaystyle\widetilde{\delta L}(\omega)$ $\displaystyle=e^{i\omega
L}\hat{x}\cdot\big{[}\delta\vec{x}_{1}(\omega)-\cos(\omega
L)\delta\vec{x}_{0}(\omega)\big{]}=\frac{e^{i\omega
L}}{\omega^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}\Big{[}e^{i\vec{k}\cdot\hat{n}L}-\cos(\omega
L)\Big{]}(i\vec{k}\cdot\hat{n})\tilde{\Phi}$ (75)
$\displaystyle\widetilde{\delta t}(\omega)$ $\displaystyle=-e^{i\omega
L}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}\left[\frac{(e^{i\vec{k}\cdot\hat{n}L}-\cos\omega
L)(i\vec{k}\cdot\hat{n})+\omega\sin\omega
L}{\omega^{2}-(\vec{k}\cdot\hat{n})^{2}}\right](\tilde{\Phi}+\tilde{\Psi})$
(76)

where we have used
$\delta\vec{x}_{i}(\omega)=-(1/\omega)^{2}\int[d^{3}k/(2\pi)^{3}]e^{i\vec{k}\cdot\vec{x}_{i}}\tilde{a}(\omega,\vec{k})$
and $\tilde{a}(k)=-i\vec{k}\tilde{\Phi}(k)$. Here
$\hat{n}=(\vec{x}_{1}-\vec{x}_{0})/L$. The phase in the frequency space is

$\displaystyle\widetilde{\Delta\phi}(\omega)$ $\displaystyle=2e^{i\omega
L}\omega_{L}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}\Big{[}\frac{\big{(}e^{i\vec{k}\cdot\hat{n}L}-\cos\omega
L\big{)}(i\vec{k}\cdot\hat{n})}{\omega^{2}}\tilde{\Phi}-\frac{(e^{i\vec{k}\cdot\hat{n}L}-\cos\omega
L)(i\vec{k}\cdot\hat{n})+\omega\sin\omega
L}{\omega^{2}-(\vec{k}\cdot\hat{n})^{2}}(\tilde{\Phi}+\tilde{\Psi})\Big{]}.$
(77)

To simplify this expression, let us consider the power spectrum for the metric
perturbation in the non-relativistic limit. They are given as

$\displaystyle P_{\Psi}(k)$ $\displaystyle=\frac{(4\pi
Gm)^{2}}{2k^{4}}(2\pi)^{4}\int
d^{3}v_{1}d^{3}v_{2}f_{1}f_{2}\Big{[}v_{c}^{4}\delta^{(4)}(k-p_{1}-p_{2})+\delta^{(4)}(k-p_{1}+p_{2})\big{)}\Big{]}+(k\to-k)$
(78) $\displaystyle P_{\Phi}(k)$ $\displaystyle=\frac{(4\pi
Gm)^{2}}{2k^{4}}(2\pi)^{4}\int
d^{3}v_{1}d^{3}v_{2}\,f_{1}f_{2}\Big{[}(v_{c}^{2}+v_{d}^{2}-3(\hat{v}_{c}\cdot\vec{v}_{d})^{2})^{2}\delta^{(4)}(k-k_{1}-k_{2})+\delta^{(4)}(k-k_{1}+k_{2})\Big{]}+(k\to-k)$
(79)

where $f_{i}=f(\vec{v}_{i})$ is the velocity distribution. In both cases, the
$\omega=2m$ mode is velocity suppressed, while $\omega<m\sigma^{2}$ does not
have such velocity suppression.

Let us first consider $\omega<m\sigma^{2}$ modes. In this case,
$\tilde{\Phi}=\tilde{\Psi}$ and $\omega/k\sim\sigma\ll 1$. Therefore, the time
delay contribution is suppressed by $(\omega/k)^{2}\sim\sigma^{2}$ compared to
the test mass acceleration. Hence, the phase is approximated as

$\displaystyle\widetilde{\Delta\phi}(\omega)=e^{i\omega
L}\frac{\omega_{L}}{\omega^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}(i\vec{k}\cdot\vec{D})\tilde{\Psi}$
(80)

with $\vec{D}(\hat{n},L)=\hat{n}(e^{i\vec{k}\cdot\hat{n}L}-\cos\omega L)$.

Now consider $\omega=2m$. In this case, $\omega/k\sim 1/\sigma\gg 1$. The
phase is approximated as

$\displaystyle\widetilde{\Delta\phi}(\omega)$
$\displaystyle\approx-2e^{i\omega
L}\frac{\omega_{L}}{\omega^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}\Big{[}(e^{i\vec{k}\cdot\hat{n}L}-\cos\omega
L)(i\vec{k}\cdot\hat{n})\tilde{\Psi}+\omega\sin\omega
L\Big{(}1+\frac{(\vec{k}\cdot\hat{n})^{2}}{\omega^{2}}\Big{)}(\tilde{\Phi}+\tilde{\Psi})\Big{]}$
$\displaystyle=-2e^{i\omega
L}\frac{\omega_{L}}{\omega^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}\Big{[}(e^{i\vec{k}\cdot\hat{n}L}-\cos\omega
L-i(\vec{k}\cdot\hat{n}L)\,{\rm sinc}\,\omega
L)(i\vec{k}\cdot\hat{n})\tilde{\Psi}+{\rm sinc}\,\omega
L(\vec{k}\cdot\hat{n})^{2}L\tilde{\Phi}\Big{]}$ (81)

In the second line, we neglect the $\omega\sin(\omega
L)(\tilde{\Phi}+\tilde{\Psi})$ term. Interferometers are always sensitive to
the phase difference, and since this term does not have any dependence on the
direction of light propagation, it cancels out in the final interferometer
observables. For $\omega L\gtrsim 1$ and for $kL<(\omega L)^{2}<1$, the terms
proportional to $(\vec{k}\cdot\hat{n})^{2}$ are negligible, and thus, the
phase is approximated as

$\displaystyle\widetilde{\Delta\phi}(\omega)$ $\displaystyle\approx-e^{i\omega
L}\frac{\omega_{L}}{\omega^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}(i\vec{k}\cdot\vec{D})\tilde{\Psi}$
(82)

For $(\omega L)^{2}<kL$, the dominant term is the one proportional to
$\tilde{\Phi}$, so a more correct approximation would be obtained by replacing
$\tilde{\Psi}\to-\tilde{\Phi}$ from (82). Since $\tilde{\Phi}\sim\tilde{\Psi}$
up to an order one numerical factor, the above expression still provides a
reasonable approximation even in this limit. In summary, the phase in both
cases can be approximated as

$\displaystyle\widetilde{\Delta\phi}(\omega)=\pm e^{i\omega
L}\frac{\omega_{L}}{\omega^{2}}\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{x}_{0}}(i\vec{k}\cdot\vec{D})\tilde{\Psi}.$
(83)

##  Appendix C Density fluctuation

We compute the statistical properties of (mass) density fluctuation. In the
non-relativistic limit, the density field becomes

$\displaystyle\hat{\rho}_{\phi}\approx\sum_{i,j}\frac{m}{V}a_{i}^{\dagger}a_{j}e^{i(k_{i}-k_{j})\cdot
x},$ (84)

where we have ignored quantum fluctuations arising from the commutation
relation. From the above expression with the quasi-probability distribution
(5), moments of the density field can be computed as

$\displaystyle\langle\rho_{\phi}^{n}\rangle=n!\rho_{0}^{n}.$ (85)

They are identical to those random variables under the exponential
distribution. Therefore, we find

$\displaystyle
p(\rho_{\phi})d\rho_{\phi}=\frac{e^{-\rho_{\phi}/\rho_{0}}}{\rho_{0}}d\rho_{\phi}$
(86)

with $\rho_{\phi}\in[0,\infty)$. The density fluctuation in the ultralight
dark matter halo is larger than the particle dark matter halo since particle
dark matter density follows the Poisson statistics.

##  Appendix D Coherent signals in pulsar timing array

The pulsar timing array (PTA) measures the time-of-arrival of pulses from
pulsars. Similar to the previous discussion, the time-of-arrival is affected
by ULDM fluctuations through Shapiro time delay, Doppler effect, and Einstein
delay. An arrival time of a pulse is given by

$\displaystyle\tau$ $\displaystyle\simeq\tau_{\rm
em}+L+\delta\tau_{1}+\delta\tau_{2}+\delta\tau_{3}$ (87)

where $\tau$ is the proper time of an observer at the observing time,
$\tau_{\rm em}$ is the proper time at the pulsar at the emission time, $L$ is
the nominal distance between the pulsar and the earth, and

$\displaystyle\delta\tau_{1}(t)$
$\displaystyle=-\hat{n}\cdot\big{[}\delta\vec{x}_{\oplus}(t)-\delta\vec{x}_{\rm
psr}(t-L)\big{]},$ (88) $\displaystyle\delta\tau_{2}(t)$
$\displaystyle=-\int_{t-L}^{t}dt^{\prime}\,\big{[}\Phi(t^{\prime},\vec{x}(t^{\prime}))+\Psi(t^{\prime},\vec{x}(t^{\prime}))\big{]},$
(89) $\displaystyle\delta\tau_{3}(t)$
$\displaystyle=\int^{t}dt^{\prime}\,\Big{[}\Phi(t^{\prime},\vec{x}_{\oplus})-\Phi(t^{\prime}-L,\vec{x}_{\rm
psr})].$ (90)

Here $\delta\tau_{1}$ is due to the change in distance between earth-pulsar
system similar to (73) (Doppler effect), $\delta\tau_{2}$ is the Shapiro time
delay similar to (74), and $\delta\tau_{3}$ is due to the conversion between
the proper time $\tau$ and the coordinate time $t$ (Einstein delay). Note that
$\hat{n}$ is the unit vector pointing to the pulsar.

The next pulse arrives at $\tau^{\prime}$,

$\displaystyle\tau^{\prime}$ $\displaystyle\simeq\tau_{\rm
em}+L+T+\delta\tau_{1}(t+T)+\delta\tau_{2}(t+T)+\delta\tau_{3}(t+T)$ (91)

where $T$ is the rotational period of the pulsar The difference in time
between the two arrival times is

$\displaystyle\tau^{\prime}-\tau=T+\Delta T$ (92)

where the time delay $\Delta T$ is

$\displaystyle\Delta T$
$\displaystyle=\sum_{i=1}^{3}\Big{[}\delta\tau_{i}(t+L)-\delta\tau_{i}(t)\Big{]}.$
(93)

The rotational period of pulsars used in PTA analysis is a few milliseconds,
$T\sim{\cal O}(10^{-3})\,{\rm sec}$, while the coherent signal we are looking
for has a frequency of a few nano Hz, $f\sim{\rm nHz}$. This suggests that the
fluctuations due to ULDM do not significantly change over the rotational time
scale $T$. In the limit $fT\ll 1$, the time decay can be simplified as

$\displaystyle\frac{\Delta T}{T}(t)$
$\displaystyle=-\hat{n}\cdot\big{[}\delta\vec{v}_{\oplus}(t)-\delta\vec{v}_{\rm
psr}(t-L)\big{]}-\Psi(t,\vec{x}_{\oplus})+\Psi(t-L,\vec{x}_{\rm psr})$
$\displaystyle+\int_{t-L}^{t}dt^{\prime}\,\hat{n}\cdot\nabla\Big{[}\Phi(t^{\prime},\vec{x}(t^{\prime}))+\Psi(t^{\prime},\vec{x}(t^{\prime}))\Big{]}.$
(94)

This (approximate) time delay expression coincides with the redshift of the
photon up to the sign [63, 65]. Note that the Doppler shift was ignored in the
work of Khmelnitsky and Rubakov [63].

The timing residual is computed by integrating the above time delay over time:

$\displaystyle R(t)=\int_{0}^{t}dt^{\prime}\frac{\Delta T}{T}(t^{\prime}).$
(95)

For the coherent oscillation at $\omega=2m$, the largest contribution to the
timing residual arises from $-\Psi(t,\vec{x}_{\oplus})+\Psi(t-L,\vec{x}_{\rm
psr})$ part in the time delay. In particular, one finds

$\displaystyle\langle
R(t)^{2}\rangle\approx\int\frac{d^{4}k}{(2\pi)^{4}}\frac{2}{\omega^{4}}P_{\Psi}(k)$
(96)

where we have ignored fast-oscillating terms. The power spectrum is given by
$P_{\Psi}(k)=[(4\pi G)^{2}/k^{4}]P_{\delta\rho}(k)$, where the density
fluctuation power spectrum around $\omega\simeq 2m$ is

$\displaystyle
P_{\delta\rho}(k)=\frac{2\pi^{2}\rho^{2}}{m^{4}\sigma^{5}}\left[\Big{(}\frac{v_{k}^{4}}{16}\Big{(}\frac{\bar{v}^{2}}{\sigma^{2}}-\frac{v_{k}^{2}}{4\sigma^{2}}\Big{)}^{1/2}e^{-\bar{v}^{2}/\sigma^{2}}+(\omega\to-\omega)\Big{)}+\frac{\sigma}{v_{k}}\exp\left(-\frac{v_{k}^{2}}{4\sigma^{2}}-\frac{\sigma^{2}|\bar{\omega}|^{2}}{v_{k}^{2}}\right)\right]$
(97)

The first term in the squared parenthesis represents the coherently
oscillating mode at $\omega=2m$, while the second term represents the
stochastic modes. This expression assumes $v_{0}=0$. Note also that
$\bar{v}^{2}(\omega)=(\omega/m-2)$, $v_{k}=k/m$, and
$\bar{\omega}=\omega/m\sigma^{2}$. An explicit computation with the coherently
oscillating modes reveals

$\displaystyle\langle R(t)^{2}\rangle\approx\frac{\pi^{2}(G\rho)^{2}}{2m^{6}}$
(98)

which agrees with the result in Ref. [63]. While our main discussion focuses
on computing the differential phase induced by ULDM in an interferometric
setup, this exercise shows that the discussion presented in the main text is
consistent with previous works.

## References

  * [1] R.D. Peccei and H.R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * [2] S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * [3] F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * [4] J. Preskill, M.B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127.
  * [5] L.F. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133.
  * [6] M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * [7] P.W. Graham, D.E. Kaplan and S. Rajendran, _Cosmological Relaxation of the Electroweak Scale_ , _Phys. Rev. Lett._ 115 (2015) 221801 [1504.07551].
  * [8] A. Arvanitaki, S. Dimopoulos, V. Gorbenko, J. Huang and K. Van Tilburg, _A small weak scale from a small cosmological constant_ , _JHEP_ 05 (2017) 071 [1609.06320].
  * [9] A. Banerjee, H. Kim and G. Perez, _Coherent relaxion dark matter_ , _Phys. Rev. D_ 100 (2019) 115026 [1810.01889].
  * [10] A. Banerjee, H. Kim, O. Matsedonskyi, G. Perez and M.S. Safronova, _Probing the Relaxed Relaxion at the Luminosity and Precision Frontiers_ , _JHEP_ 07 (2020) 153 [2004.02899].
  * [11] N. Arkani-Hamed, R.T. D’Agnolo and H.D. Kim, _Weak scale as a trigger_ , _Phys. Rev. D_ 104 (2021) 095014 [2012.04652].
  * [12] R. Tito D’Agnolo and D. Teresi, _Sliding Naturalness: New Solution to the Strong- $CP$ and Electroweak-Hierarchy Problems_, _Phys. Rev. Lett._ 128 (2022) 021803 [2106.04591].
  * [13] R. Tito D’Agnolo and D. Teresi, _Sliding naturalness: cosmological selection of the weak scale_ , _JHEP_ 02 (2022) 023 [2109.13249].
  * [14] A. Chatrchyan and G. Servant, _Relaxion Dark Matter from Stochastic Misalignment_ , 2211.15694.
  * [15] W. Hu, R. Barkana and A. Gruzinov, _Cold and fuzzy dark matter_ , _Phys. Rev. Lett._ 85 (2000) 1158 [astro-ph/0003365].
  * [16] L. Hui, _Wave Dark Matter_ , _Ann. Rev. Astron. Astrophys._ 59 (2021) 247 [2101.11735].
  * [17] H.-Y. Schive, T. Chiueh and T. Broadhurst, _Cosmic Structure as the Quantum Interference of a Coherent Dark Wave_ , _Nature Phys._ 10 (2014) 496 [1406.6586].
  * [18] H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst et al., _Understanding the Core-Halo Relation of Quantum Wave Dark Matter from 3D Simulations_ , _Phys. Rev. Lett._ 113 (2014) 261302 [1407.7762].
  * [19] L. Hui, J.P. Ostriker, S. Tremaine and E. Witten, _Ultralight scalars as cosmological dark matter_ , _Phys. Rev. D_ 95 (2017) 043541 [1610.08297].
  * [20] N. Seto and A. Cooray, _Search for small-mass black hole dark matter with space-based gravitational wave detectors_ , _Phys. Rev. D_ 70 (2004) 063512 [astro-ph/0405216].
  * [21] A.W. Adams and J.S. Bloom, _Direct detection of dark matter with space-based laser interferometers_ , astro-ph/0405266.
  * [22] E.D. Hall, R.X. Adhikari, V.V. Frolov, H. Müller, M. Pospelov and R.X. Adhikari, _Laser Interferometers as Dark Matter Detectors_ , _Phys. Rev. D_ 98 (2018) 083019 [1605.01103].
  * [23] J. Jaeckel, S. Schenk and M. Spannowsky, _Probing dark matter clumps, strings and domain walls with gravitational wave detectors_ , _Eur. Phys. J. C_ 81 (2021) 828 [2004.13724].
  * [24] S. Baum, M.A. Fedderke and P.W. Graham, _Searching for dark clumps with gravitational-wave detectors_ , _Phys. Rev. D_ 106 (2022) 063015 [2206.14832].
  * [25] A. Sesana et al., _Unveiling the gravitational universe at $\mu$-Hz frequencies_, _Exper. Astron._ 51 (2021) 1333 [1908.11391].
  * [26] M.A. Fedderke, P.W. Graham and S. Rajendran, _Asteroids for $\mu$Hz gravitational-wave detection_, _Phys. Rev. D_ 105 (2022) 103018 [2112.11431].
  * [27] M.A. Fedderke, P.W. Graham and S. Rajendran, _Gravity Gradient Noise from Asteroids_ , _Phys. Rev. D_ 103 (2021) 103017 [2011.13833].
  * [28] J.I. Read, _The Local Dark Matter Density_ , _J. Phys. G_ 41 (2014) 063101 [1404.1938].
  * [29] P.F. de Salas and A. Widmark, _Dark matter local density determination: recent observations and future prospects_ , _Rept. Prog. Phys._ 84 (2021) 104901 [2012.11477].
  * [30] R.-G. Cai, Z.-K. Guo, B. Hu, C. Liu, Y. Lu, W.-T. Ni et al., _On networks of space-based gravitational-wave detectors_ , 2305.04551.
  * [31] H. Kim and A. Lenoci, _Gravitational focusing of wave dark matter_ , _Phys. Rev. D_ 105 (2022) 063032 [2112.05718].
  * [32] A. Derevianko, _Detecting dark-matter waves with a network of precision-measurement tools_ , _Phys. Rev. A_ 97 (2018) 042506 [1605.09717].
  * [33] J.W. Foster, N.L. Rodd and B.R. Safdi, _Revealing the Dark Matter Halo with Axion Direct Detection_ , _Phys. Rev. D_ 97 (2018) 123006 [1711.10489].
  * [34] G.P. Centers et al., _Stochastic fluctuations of bosonic dark matter_ , _Nature Commun._ 12 (2021) 7321 [1905.13650].
  * [35] M. Otto, _Time-Delay Interferometry Simulations for the Laser Interferometer Space Antenna_ , Ph.D. thesis, Leibniz U., Hannover, 2015.  10.15488/8545.
  * [36] M. Tinto and S.V. Dhurandhar, _Time-delay interferometry_ , _Living Rev. Rel._ 24 (2021) 1.
  * [37] M. Tinto and J.W. Armstrong, _Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation_ , _Phys. Rev. D_ 59 (1999) 102003.
  * [38] J.W. Armstrong, F.B. Estabrook and M. Tinto, _Time-Delay Interferometry for Space-based Gravitational Wave Searches_ , _Astrophys. J._ 527 (1999) 814.
  * [39] S. Babak, A. Petiteau and M. Hewitson, _LISA Sensitivity and SNR Calculations_ , 2108.01167.
  * [40] M. Maggiore, _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford Master Series in Physics, Oxford University Press (2007).
  * [41] S. Phinney et al, _The Big Bang Observer: Direct Detection of Gravitational Waves from the Nirth of the Universe to the Presen_ , _NASA Mission Concept Study_ (2004) .
  * [42] aLIGO collaboration, _Sensitivity and performance of the Advanced LIGO detectors in the third observing run_ , _Phys. Rev. D_ 102 (2020) 062003 [2008.01301].
  * [43] S. Hild et al., _Sensitivity Studies for Third-Generation Gravitational Wave Observatories_ , _Class. Quant. Grav._ 28 (2011) 094013 [1012.0908].
  * [44] LISA collaboration, _Laser Interferometer Space Antenna_ , 1702.00786.
  * [45] M. Armano et al., _Sub-Femto- g Free Fall for Space-Based Gravitational Wave Observatories: LISA Pathfinder Results_ , _Phys. Rev. Lett._ 116 (2016) 231101.
  * [46] M. Armano et al., _Beyond the Required LISA Free-Fall Performance: New LISA Pathfinder Results down to 20 $\mu$Hz_, _Phys. Rev. Lett._ 120 (2018) 061101.
  * [47] T. Robson, N.J. Cornish and C. Liu, _The construction and use of LISA sensitivity curves_ , _Class. Quant. Grav._ 36 (2019) 105011 [1803.01944].
  * [48] N.J. Cornish, _Detecting a stochastic gravitational wave background with the Laser Interferometer Space Antenna_ , _Phys. Rev. D_ 65 (2002) 022004 [gr-qc/0106058].
  * [49] N. Cornish and T. Robson, _Galactic binary science with the new LISA design_ , _J. Phys. Conf. Ser._ 840 (2017) 012024 [1703.09858].
  * [50] J. Crowder and N.J. Cornish, _Beyond LISA: Exploring future gravitational wave missions_ , _Phys. Rev. D_ 72 (2005) 083005 [gr-qc/0506015].
  * [51] K. Yagi and N. Seto, _Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries_ , _Phys. Rev. D_ 83 (2011) 044011 [1101.3940].
  * [52] B. Bar-Or, J.-B. Fouvry and S. Tremaine, _Relaxation in a Fuzzy Dark Matter Halo_ , _Astrophys. J._ 871 (2019) 28 [1809.07673].
  * [53] N.P. Pitjev and E.V. Pitjeva, _Constraints on dark matter in the solar system_ , _Astron. Lett._ 39 (2013) 141 [1306.5534].
  * [54] S.L. Adler, _Placing direct limits on the mass of earth-bound dark matter_ , _J. Phys. A_ 41 (2008) 412002 [0808.0899].
  * [55] LIGO Scientific, VIRGO collaboration, _Improved Upper Limits on the Stochastic Gravitational-Wave Background from 2009–2010 LIGO and Virgo Data_ , _Phys. Rev. Lett._ 113 (2014) 231101 [1406.4556].
  * [56] LIGO Scientific, VIRGO collaboration, _Searching for stochastic gravitational waves using data from the two colocated LIGO Hanford detectors_ , _Phys. Rev. D_ 91 (2015) 022003 [1410.6211].
  * [57] M. Tinto, J.W. Armstrong and F.B. Estabrook, _Discriminating a gravitational wave background from instrumental noise in the LISA detector_ , _Phys. Rev. D_ 63 (2001) 021101.
  * [58] C.J. Hogan and P.L. Bender, _Estimating stochastic gravitational wave backgrounds with Sagnac calibration_ , _Phys. Rev. D_ 64 (2001) 062002 [astro-ph/0104266].
  * [59] J.D. Romano and N.J. Cornish, _Detection methods for stochastic gravitational-wave backgrounds: a unified treatment_ , _Living Rev. Rel._ 20 (2017) 2 [1608.06889].
  * [60] A. Banerjee, D. Budker, J. Eby, H. Kim and G. Perez, _Relaxion Stars and their detection via Atomic Physics_ , _Commun. Phys._ 3 (2020) 1 [1902.08212].
  * [61] A. Banerjee, D. Budker, J. Eby, V.V. Flambaum, H. Kim, O. Matsedonskyi et al., _Searching for Earth/Solar Axion Halos_ , _JHEP_ 09 (2020) 004 [1912.04295].
  * [62] D. Budker, J. Eby, M. Gorghetto, M. Jiang and G. Perez, _A Generic Formation Mechanism of Ultralight Dark Matter Solar Halos_ , 2306.12477.
  * [63] A. Khmelnitsky and V. Rubakov, _Pulsar timing signal from ultralight scalar dark matter_ , _JCAP_ 02 (2014) 019 [1309.5888].
  * [64] H. Kim and A. Mitridate, _to appear_ , .
  * [65] M. Maggiore, _Gravitational Waves. Vol. 2: Astrophysics and Cosmology_ , Oxford University Press (3, 2018).