# Sensitivity of Space-based Gravitational-Wave Interferometers to Ultralight
Bosonic Fields and Dark Matter

Jiang-Chuan Yua,c Two authors contribute equally.  Yue-Hui Yaoa,c Two authors
contribute equally.  Yong Tanga,b,c,d Yue-Liang Wua,b,c,e aUniversity of
Chinese Academy of Sciences (UCAS), Beijing 100049, China
bSchool of Fundamental Physics and Mathematical Sciences,
Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China
cInternational Center for Theoretical Physics Asia-Pacific, Beijing/Hangzhou,
China
dNational Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China
eInstitute of Theoretical Physics, Chinese Academy of Sciences, Beijing
100190, China

###### Abstract

Ultralight bosonic fields (ULBFs) are predicted by various theories beyond the
standard model of particle physics and are viable candidates of cold dark
matter. There have been increasing interests to search for the ULBFs in
physical and astronomical experiments. In this paper, we investigate the
sensitivity of several planned space-based gravitational-wave interferometers
to ultralight scalar and vector fields. Using time-delay interferometry (TDI)
to suppress the overwhelming laser frequency noise, we derive the averaged
transfer functions of different TDI combinations to scalar and vector fields,
and estimate the impacts of bosonic field’s velocities. We obtain the
sensitivity curves for LISA, Taiji and TianQin, and explore their projected
constraints on the couplings between ULBFs and standard model particles,
illustrating with the ULBFs as dark matter.

Dark Matter, Ultralight fields, Gravitational wave, Time-Delay Interferometry

##  I INTRODUCTION

On September 14, 2015 the LIGO detectors detected the first gravitational-wave
signal generated from a binary black hole merger [1], which opens the era of
gravitational astronomy [2]. Since then, over ninety gravitational wave events
have been detected, which not only brings us many valuable information about
the Universe but also shows the incredible ability of interferometry detector.
Precision interferometers can measure the relative displacements of two
mirrors separated several kilometers apart at subatomic level. Besides
gravitational waves, other physical processes may also generate the relative
displacement of mirrors (test masses). Thus interferometer gravitational-wave
detectors can also be used to search for the signals of those physical
processes.

This idea has been recently explored, for example in [3, 4, 5, 6], which
propose to use interferometers to search for ultralight dark matter (ULDM) [7,
8, 9, 10, 11]. Since ULDM particles are very light, typically below $1$ eV,
there are numerous ULDM particles in a de Broglie volume (for $m=10^{-17}$ eV,
the average number of particles $\sim 10^{51}$). The occupancy number is so
large that the quantum fluctuations are negligible and the set of ULDM
particles can be well described by classical bosonic fields that oscillate at
the Compton frequency. When ULDM is coupled to the SM particles, its
oscillation will effectively exert a periodical force on the test masses and
induce spatial displacements of them, which may end up with detectable signal
in the interferometer. Searching for the ULDM signals in the data of LIGO [12,
13, 14] and LISA Pathfinder [15] has been conducted recently. Although
reporting null results, these works show that high precision interferometers
can put stringent constraints on the coupling constants close to or even
comparable with constraints given by equivalence principle tests and fifth-
force experiments [16, 17, 18, 19, 20].

In the above works, it was assumed that the ultralight bosonic fields serve as
dark matter, therefore the amplitudes of fields are determined by the local
dark matter density and the frequency dispersion of fields follows from the
velocity dispersion of DM. However, besides DM, ULBFs are widely predicted in
physical and cosmological theories for various motivations, and have been
actively investigated on different forms [21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 39, 42, 43, 20, 44, 45, 46,
47, 48, 49] since many possible field configurations could exist. In this
work, we take a systematic approach and drop the assumption that ULBF makes up
the local DM, which implies the amplitude of field and the velocity of filed
are generic. We investigate the ability of interferometers to constrain the
combinations of coupling constants and field amplitudes in the more general
setup and illustrate with the DM case. In the prospects of the increasing
sensitivity of operating interferometer detectors [50, 51] and the next
generation ground-based detectors like Cosmic Explorer [52], Einstein
Telescope [53] and space-based detectors like LISA [54], Taiji [55] and
TianQin [56], we show it is possible to detect weaker ULBF signals or put more
stringent constraints on the coupling constants and amplitudes of the fields
in the null case by space-based detectors.

In this work, we focus on two kinds of ULBF, scalar and vector, in the mass
range $10^{-18}\sim 10^{-13}$ eV, which corresponds to the sensitive band of
space-based gravitational-wave interferometers. We investigate their potential
to detect ULBF signals by calculating response functions and sensitivity
curves of several main time-delay interferometry combinations [57, 58], such
as Michelson $X,Y,Z$, Sagnac $\alpha,\beta,\gamma$, fully symmetric Sagnac
$\zeta$ and optimal channel $\eta$. Our results show future space-based
detectors can provide a better probe in some parameter regions.

This paper is organized as follows. In Sec. II, we present the scalar and
vector models and explain how these ultralight bosonic fields affect the
motion of test masses. One of the key result is the single-link transfer
function of interferometer to the bosonic fields. In Sec. III, we construct
the sky and polarization averaged transfer functions of different TDI
combinations from the single-link transfer function and discuss their
asymptotic behaviors in the long wavelength limit and contrast them with the
case for gravitational wave. In particular, we show the effects of field
velocity on the transfer function. Then we calculate the sensitivity curves of
Taiji, LISA and TianQin to ULBF in Sec. IV and illustrate with dark matter
case in Sec. V. Finally, we give our conclusion in Sec. VI.

In this paper, we use natural units ($c=\hbar=1$) unless explicitly stated
otherwise.

##  II Formalism and Conventions

In this section, we introduce our conventions for scalar and vector fields,
and their ways of interacting with ordinary matter. We obtain the equation of
motion for the test masses through the action of the system and give the
spatial displacements of the test masses, which is essential for subsequent
discussions. The physical picture is that the spatial shifts affect the
optical lengths of laser beams, which induces the phase difference or
equivalently the frequency shift that will be detected by the phasemeters.

###  II.1 Scalar field

Ultralight scalar fields for spin-0 particles appear in many popular physical
and cosmological theories, such as axions, dark matter, dilaton and dark
energy models. Here, we do not specify the exact model but only present a
general formalism. We denote $\phi$ as a scalar field which linearly couples
to SM particle mass terms. Then the effective Lagrangian of the system at the
low-energy limit is given by

$\mathcal{L}=\mathcal{L}_{\phi}+\mathcal{L}_{SM}+\mathcal{L}_{\phi-SM},$ (1)

where $\mathcal{L}_{\phi}$ is the Lagrangian of the scalar field and usually
can be expressed by
$\mathcal{L}_{\phi}=-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-\frac{1}{2}m_{\phi}^{2}\phi^{2}$,
$\mathcal{L}_{SM}$ is the Lagrangian of the Standard Model and
$\mathcal{L}_{\phi-SM}$ for the interaction between the scalar field and the
ordinary matter. The details of the interaction will be presented in Sec. V.
At the moment, we can work in the low-energy limit and treat the test mass $m$
as a function of $\phi$ and the action of the test mass is given by [4]

$S=-\int m(\phi)\sqrt{-\eta_{\mu\nu}dx^{\mu}dx^{\nu}}.$ (2)

The effective scalar-matter coupling $\alpha$ is encoded in the mass and
defined by

$\alpha=\frac{1}{\kappa^{2}m(\phi)}\frac{\partial[\kappa
m(\phi)]}{\partial\phi},$ (3)

where $\kappa\equiv\sqrt{4\pi G}=\sqrt{4\pi}/M_{P}$ is the inverse of the
Planck mass and used here for a reference scale.

Without loss of generality, we use a monochromatic plane wave to describe the
scalar field as a transient source:

$\phi\left(t,\vec{x}\right)=\phi_{\vec{k}}e^{i(\omega
t-\vec{k}\cdot\vec{x}+\theta_{0})},$ (4)

where $\phi_{\vec{k}}$ is the amplitude of the wave, $\omega$ the de Broglie
angular frequency, $\vec{k}$ the momentum vector of the scalar field, and
$\theta_{0}$ an initial phase whose exact value does not affect our discussion
and can be set to zero. When the scalar field is referred as DM in our galaxy,
it has a velocity dispersion, $v\sim 10^{-3}$, which implies that the wave is
non-relativistic with $\vec{k}\approx m_{\phi}\vec{v}$ and $\omega\approx
m_{\phi}$. We can also extend the case to the relativistic condition, the
details will be discussed later.

Using the Euler-Lagrangian equation and neglecting the subdominant term
related to the velocity of test mass, we get the additional acceleration and
displacement induced by the oscillation of scalar field:

$a^{i}(t,\vec{x})=i\phi_{\vec{k}}\;\alpha\kappa\;k^{i}e^{i(\omega
t-\vec{k}\cdot\vec{x})}\simeq
i\phi_{\vec{k}}\;\alpha\kappa\;k^{i}e^{im_{\phi}(t-v\hat{k}\cdot\vec{x})}.$
(5)

We can integrate the above equations twice and obtain the spatial
displacement,

$\delta
x^{i}(t,\vec{x})=-i\alpha\kappa\;\phi_{\vec{k}}\frac{k^{i}}{m_{\phi}^{2}}e^{im_{\phi}(t-v\hat{k}\cdot\vec{x})}=-i\mathcal{M}_{s}\hat{k}^{i}e^{im_{\phi}(t-v\hat{k}\cdot\vec{x})},$
(6)

where $\hat{k}$ is the unit vector of $\vec{k}$ and
$\mathcal{M}_{s}=\alpha\kappa\phi_{\vec{k}}|\vec{k}|/m_{\phi}^{2}$.

###  II.2 Vector field

Ultralight vector fields for spin-1 particles are also well motivated
theoretically, such as very light scaling gauge field for DM candidate [59,
60, 61, 32, 62]. Vector field can directly couple with the baryon number (B)
or baryon minus lepton number (B-L) and it will exert an additional force on
the test masses, which leads to the periodic displacements of test masses.

The Lagrangian of the massive vector field coupled with B or B-L current
$J_{D}$ is given by

$\mathcal{L}=-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}+\frac{1}{2}m_{A}^{2}A^{\nu}A_{\nu}-\epsilon_{D}eJ^{\nu}_{D}A_{\nu},$
(7)

where $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, $A^{\mu}$ is
the vector field and $m_{A}$ is the mass of vector field. Here we normalize
the coupling strength of the vector field $e_{D}$ in terms of the
electromagnetic coupling constant $e$ and denote $\epsilon_{D}$ as the ratio
$\epsilon_{D}=e_{D}/e$ [63].

The spatial components of $A^{\mu}$ are described by

$\vec{A}(t,\vec{x})=|\vec{A}|\hat{e}_{A}e^{i(\omega t-\vec{k}\cdot\vec{x})},$
(8)

where $\hat{e}_{A}$ is the unit polarization vector which is parallel to
$\vec{A}$. Here we ignore $A^{0}(t,\vec{x})$ since it is much smaller than
$A^{i}(t,\vec{x})$ under Lorentz gauge in the non-relativistic limit. A object
carrying B or B-L charge can be accelerated by the electric-like field
generated by the vector field. The acceleration acting on the test mass can be
written as [63, 3]

$\displaystyle{a^{i}}(t,\vec{x})=\epsilon_{D}e\frac{q_{D,j}}{M_{j}}\partial_{t}A^{i}(t,\vec{x})\simeq
i\epsilon_{D}e\frac{q_{D,j}}{M_{j}}m_{A}|\vec{A}|\hat{e}_{A}^{i}e^{im_{A}(t-v\hat{k}\cdot\vec{x})}.$
(9)

where $M_{j}$ is the mass of the jth test mass. $q_{D,j}$ is the total baryon
number if the vector field is a $U(1)_{B}$ gauge field, and only counts the
neutrons in the test masses if it is associated with $U(1)_{B-L}$ gauge
symmetry. We can integrate Eq. (9) twice and obtain the spatial displacement
of the test mass:

$\delta
x^{i}(t,\vec{x})=-i\epsilon_{D}e\frac{q_{D,j}}{M_{j}}\frac{|\vec{A}|}{m_{A}}\hat{e}_{A}^{i}e^{im_{A}(t-v\hat{k}\cdot\vec{x})}=-i\mathcal{M}_{v}\hat{e}_{A}^{i}e^{im_{A}(t-v\hat{k}\cdot\vec{x})},$
(10)

where $\mathcal{M}_{v}=\epsilon_{D}eq_{D,j}|\vec{A}|/m_{A}M_{j}$. Comparing it
with the scalar case, Eq. (6), we can find that the difference is that
$\hat{k}^{i}$ in Eq. (6) is replaced by $\hat{e}_{A}^{i}$ in Eq. (10).
Therefore, both transverse and longitudinal modes contribute.

##  III Transfer function

In this section, we discuss how the displacement signals induced by ULBFs
manifest in the physical measurements. We establish the necessary formalism to
describe the transfer function, including the detector response of a single
link and time-delay interferometry.

###  III.1 The detector frame

The detector coordinate system is chosen such that the three arms of the
triangle detector keep fixed in the $x-z$ plane with the first spacecraft (SC)
located at $\vec{x}_{1}=(0,0,0)$, the second located at $\vec{x}_{2}=(0,0,L)$
and the third located at $\vec{x}_{3}=(L\sin\gamma,0,L\cos\gamma)$, where
$\gamma$ is the opening angle between the two arms of the interferometer and
$L$ is the static length of the arm. The wave propagation direction is

$\hat{k}=(\sin\theta_{1}\cos\epsilon_{1},\sin\theta_{1}\sin\epsilon_{1},\cos\theta_{1}),$

and two orthogonal vectors
$\hat{u}=(\cos\theta_{1}\cos\epsilon_{1},\cos\theta_{1}\sin\epsilon_{1},-\sin\theta_{1})$,
$\hat{v}=(-\sin\epsilon_{1},\cos\epsilon_{1},0)$ can be chosen as the
polarization basis vectors. The details of coordinate system is shown in Fig.
1.

Figure 1: Schematic of the detector. 1, 2, 3 represent the labels of
spacecrafts. Three spacecrafts are located at $\vec{x}_{1}=(0,0,0)$,
$\vec{x}_{2}=(0,0,L)$ , and $\vec{x}_{3}=(L\sin\gamma,0,L\cos\gamma)$
respectively. The static distance between two spacecrafts is $L$, and the
angle between two arms is $\gamma=\frac{\pi}{3}$. $\hat{n}_{rs}$ refers to the
unit vector in the direction from $\vec{x}_{r}$ to $\vec{x}_{s}$. The spatial
direction of the wave vector $\vec{k}$ is
$\hat{k}=(\sin\theta_{1}\cos\epsilon_{1},\sin\theta_{1}\sin\epsilon_{1},\cos\theta_{1})$,
and
$\hat{u}=(\cos\theta_{1}\cos\epsilon_{1},\cos\theta_{1}\sin\epsilon_{1},-\sin\theta_{1})$,
$\hat{v}=(-\sin\epsilon_{1}\cos\epsilon_{1},0)$. They form an orthogonal
coordinate system.

###  III.2 Transfer functions of a single link

Now we can calculate the frequency shift of laser caused by the motion of test
masses and we will derive the one-link transfer functions in this section,
which are essential for the discussions of transfer functions of TDI
combinations.

Since the test masses will oscillate under the influence of ULBF, the time
$T_{rs}$ that the laser need to propagate from sender SC ($\vec{x}_{s}$) to
receiver SC ($\vec{x}_{r}$) will vary with time,

$T_{rs}=L_{rs}+\delta t_{rs},$ (11)

where $L_{rs}$ is the static length of the arm in the absence of ULBF and
$\delta t_{rs}$ is perturbation of one-way time due to the ULBF oscillation,

$\delta
t_{rs}=-\hat{n}_{rs}\cdot\left[\delta\vec{x}(t,\vec{x}_{r})-\delta\vec{x}(t-L,\vec{x}_{s})\right],$
(12)

where $\hat{n}_{rs}$ is the unit vector pointing from $\vec{x}_{r}$ to
$\vec{x}_{s}$ and $t$ is the arrival time of laser recorded by the receiver.
The laser phase $\phi_{r}$ received at $t$ equals the laser phase $\phi_{s}$
sent at $t-T_{rs}$,

$\phi_{r}(t)=\phi_{s}(t-T_{rs})=2\pi\nu_{0}(t-L_{rs}-\delta t_{rs}),$ (13)

where $\nu_{0}$ is the fiducial frequency of laser. The instantaneous
frequency $\nu$ measured on the receiving spacecraft is related to $\phi_{r}$
through

$\displaystyle\nu_{rs}(t)=\frac{1}{2\pi}\frac{d\phi_{r}}{dt}.$ (14)

Thus, the relative frequency shift or fractional frequency fluctuation is

$\frac{\delta\nu_{rs}}{\nu_{0}}=\frac{\nu_{rs}-\nu_{0}}{\nu_{0}}=-\frac{d\,\delta
t_{rs}}{dt}.$ (15)

Combining Eq. (12) and Eq. (15), one-way Doppler shift can be expressed as

$y_{rs}(t)\equiv\frac{\delta\nu_{rs}}{\nu_{0}}=\mu_{rs}\;[h(t,\vec{x_{r}})-h(t-L,\vec{x}_{s})],$
(16)

where $h(t,\vec{x})=m\mathcal{M}e^{im(t-v\hat{k}\cdot\vec{x})}$ and $\mu_{rs}$
is the geometric factor given by

$\mu_{rs}=\begin{cases}\hat{k}\cdot\hat{n}_{rs}\;&\;\;\text{for scalar
field},\\\ \hat{e}_{A}\cdot\hat{n}_{rs}\;&\;\;\text{for vector field},\\\
\cfrac{\hat{n}^{i}_{rs}\hat{n}^{j}_{rs}{e}_{ij}(\hat{k},\psi)}{2(1+\hat{n}_{rs}\cdot\hat{k})}\;&\;\;\text{for
gravitational wave},\\\ \end{cases}$ (17)

where we have included the response for gravitational wave for comparison,
$e_{ij}(\hat{k},\psi)$ is the polarization tensor and $\psi$ is the
polarization angle for gravitational wave. As can be seen from Eq. (16), with
the definition of $\mu_{rs}$, the one-way Doppler signals of scalar fields,
vector fields and gravitational waves can be written in a concise and unified
form. Follow from it, the TDI combinations, which are composed of time-shifted
one-way Doppler signals, also have identical forms for scalar fields, vector
fields and gravitational waves.

If the Doppler shift is recorded within a period of time $T$, the
corresponding Fourier transform of $h(t)$ becomes

$h(t)=\frac{\sqrt{T}}{2\pi}\int_{0}^{T}\tilde{h}(\omega)e^{i\omega t}d\omega.$
(18)

Putting Eq.(18) in Eq.(16), we obtain

$\displaystyle
y_{rs}(t)=\mu_{rs}\;\frac{\sqrt{T}}{2\pi}\int_{0}^{T}d\omega\;\tilde{h}(\omega)e^{i\omega
t}\left[e^{-i(\vec{k}\cdot\vec{x}_{r})}-e^{-i(\tau+\vec{k}\cdot\vec{x}_{s})}\right],$
(19)

where $\tau=\omega L\equiv 2\pi fL$ and $\tilde{h}(\omega)$ is the Fourier
amplitude of $h(t)$. Then we can extract the Fourier transformation of
$y_{rs}(t)$ from Eq. (19):

$\tilde{y}_{rs}(\omega)=\mu_{rs}\;\tilde{h}(\omega)\left[e^{-i(\vec{k}\cdot\vec{x}_{r})}-e^{-i(\tau+\vec{k}\cdot\vec{x}_{s})}\right].$
(20)

Here we have neglected the time-dependence of $\hat{n}_{rs}$ in $\mu_{rs}$,
which follows the convention in the gravitational-wave literature and is
compensated by the average of different directions.

The sky and polarization averaged transfer function $R(\omega)$ is defined as

$R(\omega)=\overline{\left|\frac{\tilde{y}_{rs}(\omega)}{\tilde{h}(\omega)}\right|^{2}},$
(21)

where the overline denotes the average over all propagation directions and
polarization angles. Explicitly, for scalar field, we have

$R^{s}(\omega)=\frac{1}{4\pi}\int_{-1}^{1}d\cos\theta_{1}\int_{0}^{2\pi}d\epsilon_{1}\left|\frac{\tilde{y}_{rs}(\omega,\theta_{1},\epsilon_{1})}{\tilde{h}(\omega)}\right|^{2}.$
(22)

For vector field, we have

$R^{v}(\omega)=\frac{1}{16\pi^{2}}\int_{-1}^{1}d\cos\theta_{1}\int_{0}^{2\pi}d\epsilon_{1}\int_{-1}^{1}d\cos\theta_{2}\int_{0}^{2\pi}d\epsilon_{2}\left|\frac{\tilde{y}_{rs}(\omega,\theta_{1},\epsilon_{1},\theta_{2},\epsilon_{2})}{\tilde{h}(\omega)}\right|^{2},$
(23)

where $(\theta_{2},\epsilon_{2})$ is the two-dimensional spherical coordinates
of $\hat{e}_{A}$. For gravitational waves, we have

$R^{GW}(\omega)=\frac{1}{8\pi^{2}}\int_{-1}^{1}d\cos\theta_{1}\int_{0}^{2\pi}d\epsilon_{1}\int_{0}^{2\pi}d\psi\left|\frac{\tilde{y}_{rs}(\omega,\theta_{1},\epsilon_{1},\psi)}{\tilde{h}(\omega)}\right|^{2}.$
(24)

We define the integration operations for later convenience,

$\displaystyle I_{s}$
$\displaystyle\equiv\frac{1}{4\pi}\int_{-1}^{1}d\cos\theta_{1}\int_{0}^{2\pi}d\epsilon_{1}\cdots,$
(25) $\displaystyle I_{v}$
$\displaystyle\equiv\frac{1}{16\pi^{2}}\int_{-1}^{1}d\cos\theta_{1}\int_{0}^{2\pi}d\epsilon_{1}\int_{-1}^{1}d\cos\theta_{2}\int_{0}^{2\pi}d\epsilon_{2}\cdots,$
(26) $\displaystyle I_{GW}$
$\displaystyle\equiv\frac{1}{8\pi^{2}}\int_{-1}^{1}d\cos\theta_{1}\int_{0}^{2\pi}d\epsilon_{1}\int_{0}^{2\pi}d\psi\cdots.$
(27)

###  III.3 Transfer functions of TDI combinations

As can be seen from Eq. (19), the signal in one-link is a modulated Doppler
shift with amplitude proportional to $\tilde{h}(\omega)$ whose typical
magnitudes for astrophysical sources are of order $10^{-21}$. However, the
current baseline for stabilized laser has a power spectral density as high as
$30\;\text{Hz}/\sqrt{\text{Hz}}$ or $10^{-13}/\sqrt{\text{Hz}}$ in relative
frequency fluctuation for laser with wavelength $\lambda=1064$ nm. Thus the
laser frequency noise dominates in the one-link signals. In ground-based
interferometer detectors, we do not need to worry about the laser noise in one
link signals, since it will be suppressed significantly when the beams are
combined at the photon detector due to the equal length of the two arms of
Michelson interferometer. However, for space-borne interferometer detectors,
the distances between pairs of spacecrafts evolve with time and can not keep
equal with each other. Therefore the laser noise will not be cancelled by the
simple Michelson configuration. To tackle this problem, TDI is used [57, 64]
for data postprocessing. The essential idea is to use the data of six links to
synthesize virtual interferometric configurations that have almost equal light
paths. It has been shown that the magnitude of residual laser noise in the
synthesized data meets the requirement of sensitivity of gravitational waves
detection.

In this section, we discuss the response of the TDI-1.5 configurations which
apply to a non-equal arms, rigid but rotating triangle constellation. For the
sake of simplicity and physical intuition, we focus on the equilateral-
triangle configuration and ignore the motion of spacecrafts, since the final
sensitivity curves do not depend on which generation TDI is. We perform the
sky and polarization average of the transfer functions in both semi-analytical
numerical integration and Monte Carlo integration with $39000$ source
positions per frequency bin as a cross check of our results.

We follow the convention used in [65] for TDI discussions. The opposite arm of
spacecraft $i$ is denoted by $L_{i}$ (counter-clockwise) and $L_{i^{\prime}}$
(clockwise), respectively. $y(t)_{rs,ij^{\prime}}=y(t-L_{i}-L_{j^{\prime}})$
is the shorthand for time-delayed series.

####  III.3.1 Michelson interferometry

We start with equal-arm Michelson configuration $M(t)$, which can be
synthesized from the four one-link signals,

$M(t)=[y_{12}+y_{21,3}]-[y_{13}+y_{31,2^{\prime}}].$ (28)

And in the static and equal armlength case, it reduces to

$M(t)=[y_{12}(t)+y_{21}(t-L)]-[y_{13}(t)+y_{31}(t-L)].$ (29)

The corresponding Fourier transform $\tilde{M}(\omega)$ is given by

$\tilde{M}(\omega)=\tilde{h}(\omega)\;\left\\{(\mu_{12}-\mu_{13})(1+e^{-i2\tau})e^{-i\vec{k}\cdot\vec{x}_{1}}-2\left[\mu_{12}e^{-i(\tau+\vec{k}\cdot\vec{x}_{2})}-\mu_{13}e^{-i(\tau+\vec{k}\cdot\vec{x}_{3})}\right]\right\\},$
(30)

where $\tau\equiv 2\pi fL$. After substituting Eq. (30) into Eqs. (22) and
(23), the sky and polarization averaged transfer functions $R^{s}_{M}$,
$R^{v}_{M}$ of Michelson configuration can be expressed as

$\displaystyle R^{s}_{M}(\omega)$ $\displaystyle=$ $\displaystyle
I_{s}\left[A_{1}^{2}(\theta_{1},\epsilon_{1},\omega)+A_{2}^{2}(\theta_{1},\epsilon_{1},\omega)\right],$
(31) $\displaystyle R^{v}_{M}(\omega)$ $\displaystyle=$ $\displaystyle
I_{v}\left[A_{1}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)+A_{2}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)\right],$
(32)

and the expressions of $A_{1}$ and $A_{2}$ can be found in Appendix A.

####  III.3.2 Michelson combinations

The Michelson $X$ combination in TDI-1.5 is given by

$X(t)=\left[y_{13}+y_{31,2^{\prime}}+y_{12,22^{\prime}}+y_{21,322^{\prime}}\right]-\left[y_{12}+y_{21,3}+y_{13,3^{\prime}3}+y_{31,2^{\prime}3^{\prime}3}\right].$
(33)

For the static and equal armlength case,

$\displaystyle X(t)$ $\displaystyle=$
$\displaystyle\left[y_{13}(t)+y_{31}(t-L)+y_{12}(t-2L)+y_{21}(t-3L)\right]$
(34)
$\displaystyle{}-\left[y_{12}(t)+y_{21}(t-L)+y_{13}(t-2L)+y_{31}(t-3L)\right].$

Michelson $Y$ and $Z$ combinations can be obtained by cyclic permutation of
the indices $(1\rightarrow 2\rightarrow 3\rightarrow 1)$. The Fourier
transform $\tilde{X}(\omega)$ is related to $\tilde{M}(\omega)$ through

$\tilde{X}(\omega)=-(1-e^{-i2\tau})\tilde{M}(\omega).$ (35)

This relation is a general consequence of the similarity between structures of
Eq. (29) and Eq. (34). Using Eq. (35), the averaged transfer functions
$R^{s}$, $R^{v}$ of $X$ combination can be expressed as

$\displaystyle R^{s}_{X}(\omega)$ $\displaystyle=$ $\displaystyle
4\sin^{2}\tau\;R^{s}_{M}(\omega),$ (36) $\displaystyle R^{v}_{X}(\omega)$
$\displaystyle=$ $\displaystyle 4\sin^{2}\tau\;R^{v}_{M}(\omega).$ (37)

After sky and polarization averaged, the transfer functions of $Y$, $Z$
combinations are the same as Eq. (36, 37).

In the first plot in Fig. 2, we show the transfer functions for Michelson
combinations. We can see that in the low frequency range $2\times
10^{-5}<f<0.1$ Hz which corresponds to $v<\tau<1$ ($v$ is the velocity of the
field), $R^{s}_{X},R^{v}_{X}\propto f^{6}$. However, when $f<2\times 10^{-5}$
Hz, namely, $\tau<v$, $R^{s}_{X},R^{v}_{X}\propto f^{4}$. This is accord with
the asymptotic behavior of GWs, $R^{GW}_{X}\propto f^{4}$. Thus, $f_{c}=v/2\pi
L$ is the critical frequency at which the asymptotic behaviors change. This
effect is due to the velocity of the scalar and vector fields.

In the long wavelength limit, we can perform the sky average analytically and
Eq. (36) reduces to

$R^{s}_{X}(\omega)\simeq
16\;\left[\frac{4}{3}\sin^{2}\left(\frac{\gamma}{2}\right)\frac{\tau^{6}}{4}+\frac{4}{15}\sin^{2}\gamma\cdot
v^{2}\tau^{4}\right].$ (38)

For non-relativistic fields, the field velocity $v$ is a small quantity. The
first term is dominant when $v<\tau$ and the second term become more important
when $\tau<v$.

####  III.3.3 Sagnac combinations: $\alpha,\beta,\gamma$

The six-link Sagnac combinations $\alpha$, $\beta$, $\gamma$ are the
generators of the space of laser noise free combinations. The expression of
$\alpha$ combination is

$\alpha(t)=\left[y_{13}+y_{32,2^{\prime}}+y_{21,1^{\prime}2^{\prime}}\right]-\left[y_{12}+y_{23,3}+y_{31,13}\right],$
(39)

which in the static and equal armlength case reduces to

$\alpha(t)=\left[y_{13}(t)+y_{32}(t-L)+y_{21}(t-2L)\right]-\left[y_{12}(t)+y_{23}(t-L)+y_{31}(t-2L)\right].$
(40)

$\beta$, $\gamma$ can be obtained by cyclic permutation of the indices
(1,2,3). The transfer function $R^{s}$, $R^{v}$ of $\alpha$ combination can be
expressed as

$\displaystyle R^{s}_{\alpha}(\omega)$ $\displaystyle=$ $\displaystyle
I_{s}[B_{1}^{2}(\theta_{1},\epsilon_{1},\omega)+B_{2}^{2}(\theta_{1},\epsilon_{1},\omega)],$
(41) $\displaystyle R^{v}_{\alpha}(\omega)$ $\displaystyle=$ $\displaystyle
I_{v}[B_{1}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)+B_{2}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)],$
(42)

where the detail expressions of $B_{1}$ and $B_{2}$ are presented in the
Appendix. The averaged transfer functions of $\beta$ and $\gamma$ are the same
as Eq. (41, 42).

The asymptotic behaviors of Sagnac combinations in the low frequency region
are the same as the behaviors of Michelson combinations. In Fig. 2, we can see
that in the low frequency range $2\times 10^{-5}<f<0.1$ Hz,
$R^{s}_{\alpha},R^{v}_{\alpha}\propto f^{6}$. However, when $f<2\times
10^{-5}$ Hz, $R^{s}_{\alpha},R^{v}_{\alpha}\propto f^{4}$, same as the
asymptotic behavior of GWs, $R^{GW}_{\alpha}\propto f^{4}$.

####  III.3.4 Fully symmetric Sagnac combination: $\zeta$

Another important combination is the fully symmetric combination $\zeta$,
given by

$\zeta(t)=y_{21,2^{\prime}}-y_{12,1}+y_{13,1^{\prime}}-y_{31,3}+y_{32,3^{\prime}}-y_{23,2}.$
(43)

For the static and equal armlength case, we have

$\zeta(t)=y_{21}(t-L)-y_{12}(t-L)+y_{13}(t-L)-y_{31}(t-L)+y_{32}(t-L)-y_{23}(t-L).$
(44)

The averaged transfer functions of $\zeta$ combination can be expressed as

$\displaystyle R^{s}_{\zeta}(\omega)$ $\displaystyle=$ $\displaystyle
4\sin^{2}\left(\frac{\tau}{2}\right)\;I_{s}\left[C_{1}^{2}(\theta_{1},\epsilon_{1},\omega)+C_{2}^{2}(\theta_{1},\epsilon_{1},\omega)\right],$
(45) $\displaystyle R^{v}_{\zeta}(\omega)$ $\displaystyle=$ $\displaystyle
4\sin^{2}\left(\frac{\tau}{2}\right)\;I_{v}\left[C_{1}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)+C_{2}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)\right].$
(46)

The expressions of $C_{1}$ and $C_{2}$ can be found in the Appendix.

From Fig. 2, it is clear to see that in the low frequency limit $f<0.1$ Hz,
$R^{s}_{\zeta}$, $R^{GW}_{\zeta}\propto f^{6}$ and $R^{v}_{\zeta}\propto
f^{4}$. It is worth noting that the $f^{6}$ asymptotic behaviors of scalar
field and GWs are true only when the detector array is equilateral triangle.
If the array has unequal armlength, then $R^{s}_{\zeta},R^{GW}_{\zeta}$ will
be proportional to $f^{4}$ in the low frequency limit. At frequencies equal
integers of the inverse of light travel time, the response of the detector is
subtracted to zero due to the overall factor $\sin^{2}\left(\tau/2\right)$ in
Eq. (45, 46).

There are two interesting points here. The first is that $\zeta$ is not
sensitive to the field velocity in the sense that there is only one asymptotic
behavior of transfer function in the low frequency region and no critical
frequency compared to Michelson combinations. The second is that the
unexpected asymptotic behavior of vector field which is proportional to
$f^{4}$ rather than $f^{6}$.

To understand these phenomena, we derive the approximate expression of
$\tilde{\zeta}(\omega)$ in the long wavelength limit,

$\tilde{\zeta}(\omega)\propto\left[\frac{1}{3}(\mu_{12}\kappa_{12}+\mu_{23}\kappa_{23}+\mu_{31}\kappa_{31})+\frac{2}{3}(\mu_{12}\kappa_{23}+\mu_{23}\kappa_{31}+\mu_{31}\kappa_{12})\right](v\tau)^{2}+O((v\tau)^{3}),$
(47)

where $\kappa_{ij}\equiv\hat{k}\cdot\hat{n}_{ij}$. It shows the field velocity
$v$ and $\tau$ appear as a combination in Eq. (47). Furthermore, for scalar
field and the longitudinal polarization of vector field,
$\mu_{ij}=\kappa_{ij}$, and the coefficient of $\tau^{2}$ cancels due to the
geometrical relation $\kappa_{12}+\kappa_{23}+\kappa_{31}=0$. Thus, the
leading-order contribution is $O(\tau^{3})$. However, for the transverse
polarization of vector field, this cancellation doesn’t happen and the leading
order is $O(\tau^{2})$. We further decompose the $R^{v}_{\zeta}$ into
$R^{v,t}_{\zeta}$ for transverse mode and $R^{v,l}_{\zeta}$ for longitudinal
mode. Following Eq. (47), we find that $R^{v,l}_{\zeta}\propto f^{6}$ and
$R^{v,t}_{\zeta}\propto f^{4}$ in the low frequency limit in the equal
armlength case. Thus, transverse polarization dominates the asymptotic
behavior of $R^{v}_{\zeta}$ in the low frequency region.

Besides serving as the noise monitor, $\zeta$ combination may have another
advantage. If a monochromatic signal is detected, we need to distinguish ULBFs
from galactic white dwarf binaries that can also produce nearly monochromatic
gravitational waves. We can utilize the difference of transfer functions
between $X$ and $\zeta$ combinations for ULBFs and GWs. Given a ULBF and a
white dwarf binary which produce the same magnitude signals in the $X$
combination, the signal of ULBF in $\zeta$ combination may be overwhelmed by
the noise due to the smaller value of transfer function while the signal of
GWs in $\zeta$ combination will survive. This suggests that the ratios of
signal-to-noise in two combinations are different in two cases.

Figure 2: The transfer functions of $X$, $\alpha$, $\zeta$, $A$, $E$ and $T$.
We assume an equilateral triangle configuration with $L\simeq 3\times 10^{9}$
m. For $X$, $\alpha$, $A$ and $E$, the transfer functions for scalar and
vector almost overlap each other across the whole band.

####  III.3.5 Optimal combinations

The Michelson combinations $X$, $Y$ and $Z$ can be further assembled into
three optimal combinations, $A,E$ and $T$ [66],

$\displaystyle A(t)$ $\displaystyle=\frac{1}{\sqrt{2}}\left[Z(t)-X(t)\right],$
$\displaystyle E(t)$
$\displaystyle=\frac{1}{\sqrt{6}}\left[X(t)-2Y(t)+Z(t)\right],$ $\displaystyle
T(t)$ $\displaystyle=\frac{1}{\sqrt{3}}\left[X(t)+Y(t)+Z(t)\right].$ (48)

The expressions of the averaged transfer functions of $A,E$, and $T$ can be
found in the Appendix. As can be seen in Fig. 2, $R^{s}_{A}$, $R^{v}_{A}$
$\propto f^{6}$ in $2\times 10^{-5}<f<0.1$ Hz and $R^{s}_{A}$, $R^{v}_{A}$
$\propto f^{4}$ when $f<2\times 10^{-5}$ Hz (the values of $R^{s,v}_{A}$ and
$R^{s,v}_{E}$ are equal across the band). In the equal armlength case,
$R^{s}_{T}$, $R^{GW}_{T}$ $\propto f^{10}$ and $R^{v}_{T}$ $\propto f^{8}$ in
the low frequency region. The different behavior of vector field is also due
to the similar reason that the transverse polarization ($\propto f^{8}$)
dominates the longitudinal polarization ($\propto f^{10}$) as explained in
$\zeta$.

###  III.4 Comparison of transfer functions of ULBFs with GWs

In the previous calculations, we adopted the assumption that $v\sim 10^{-3}$,
which results in a different shape of transfer function compared to the case
of gravitational waves, including the oscillating behavior in the high-
frequency region and the trends in low-frequency region. Now we give a
physical explanation for the difference.

In previous sections, we obtain the displacement caused by bosonic fields in
the non-relativistic case. In fact, we can directly apply the displacement
formula of the detector under the influence of an ultralight scalar field to
the relativistic case. This is because in the calculations we did not impose
any assumption on the velocity of scalar field. Thus, the displacement formula
derived earlier can still be used. But one should notice that, if the scalar
field is relativistic, the approximation of $\omega\sim m$ is no longer
applicable and should be replaced with $\omega\sim m/\sqrt{1-v^{2}}$.

For simplicity, we will first compare the single-link responses and then
discuss the transfer functions of TDI combinations.

The Fourier amplitude of the single-link signal $y_{rs}(t)$ is given by
Eq.(20) and the corresponding transfer function is

$R^{s}_{single}(\omega)=I_{s}\mu_{rs}^{2}\left|e^{-i(\vec{k}\cdot\vec{x}_{r})}-e^{-i(\tau+\vec{k}\cdot\vec{x}_{s})}\right|^{2}=I_{s}\mu_{rs}^{2}\left[2-2\cos(\tau+L\vec{k}\cdot\hat{n}_{rs})\right].$
(49)

For scalar field, $|\vec{k}|=\omega v$, different from $|\vec{k}|=\omega$ for
gravitational waves. Without lose of generality, we use the arm $\hat{n}_{23}$
to calculate the single-link response function in the case of scalar field and
gravitational waves. The result is shown in Fig. 3. In the low-frequency
region, the transfer function of the scalar field differs from that of the
gravitational waves by a velocity-dependent factor. The leading-order term of
the ratio in low frequency approximation is

$\frac{R^{s}_{single}(\omega)}{R^{GW}_{single}(\omega)}\approx\frac{I_{s}\mu_{rs}^{2}(1+v\hat{k}\cdot\hat{n}_{rs})^{2}}{I_{GW}\mu_{rs}^{2}(1+\hat{k}\cdot\hat{n}_{rs})^{2}}=5+3v^{2}.$
(50)

In the high frequency region, the velocity of the field $v$ affects the
oscillation behavior. To be more specific, as $v$ increases, the second term
related to the propagation direction in the trigonometric function in Eq. (49)
becomes non-negligible, which will lead to the change of the oscillation
period. At the same time, averaging the entire sky will result in a decrease
in the oscillation amplitude, since while cancellation appear in some
directions at a fixed frequency, it will not happen in other directions.

Now we consider the TDI combinations. The differences between transfer
functions are similar to the single-link case. We illustrate with the transfer
function of $X$ combination with respect to the frequency and the velocity in
Fig. 4.

Figure 3: The single-link transfer function of scalar field. Here we use
$L\simeq 3\times 10^{9}$ m and compare the transfer function of gravitational
wave with that of scalar field with velocities $10^{-3}$, $0.1$ and $0.9$,
respectively.

Figure 4: The transfer function of $X$ combination for scalar field, using
velocity of scalar field and frequency as independent variables. We use the
arm length $L\simeq 3\times 10^{9}$ m.

From the result, we find that with the increasing velocity, the overall trend
of the transfer function and the oscillation behavior of the high frequency
part are closer to the behavior of gravitational wave. When the velocity of
the scalar field approaches the speed of light, there is only a slight
difference between the two cases, which is caused by $\mu_{rs}$.

##  IV Sensitivity curves

Figure 5: The sensitivity curves of LISA, Taiji and TianQin to scalar field
(top), vector field (middle) and gravitational wave (bottom). We adopt
$L_{\textrm{LISA}}=2.5\times 10^{9}$ m, $L_{\textrm{Taiji}}=3\times 10^{9}$ m
and $L_{\textrm{TianQin}}=1.7\times 10^{8}$ m in the calculation.

The sensitivity of individual TDI combination $O$ is defined by [67]

$S_{O}(f)=\frac{N_{O}(f)}{R_{O}(f)},$ (51)

where $N_{O}(f)$ is the one-sided noise power spectral density (PSD) of
combination $O$ and $R_{O}(f)$ is the averaged transfer function of $O$.

After suppression of the laser noise below the requirement, the two secondary
noises determine the performance of detector. The first one is related to
optical metrology system (oms) noise, which dominates at high frequency. The
second one is the test mass acceleration noise (acc) that dominates at low
frequency. The magnitudes of oms and acc noises for LISA [54], Taiji [55],
TianQin [56] are

$\displaystyle S_{oms}\left(f\right)$ $\displaystyle=$
$\displaystyle\left(s_{oms}\;\frac{2\pi
f}{c}\right)^{2}\left[1+\left(\frac{2\times
10^{-3}~{}\textrm{Hz}}{f}\right)^{4}\right]\;\frac{1}{\textrm{Hz}},$ (52)
$\displaystyle S_{acc}\left(f\right)$ $\displaystyle=$
$\displaystyle\left(\frac{s_{acc}}{2\pi
fc}\right)^{2}\left[1+\left(\frac{0.4\times
10^{-3}~{}\textrm{Hz}}{f}\right)^{2}\right]\;\left[1+\left(\frac{f}{8\times
10^{-3}~{}\textrm{Hz}}\right)^{4}\right]\;\frac{1}{\textrm{Hz}},$ (53)

where we use the following noise parameters,

LISA $\displaystyle:$ $\displaystyle~{}s_{oms}=15\times
10^{-12}~{}\textrm{m},s_{acc}=3\times 10^{-15}~{}\textrm{m}/\textrm{s}^{2},$
Taiji $\displaystyle:$ $\displaystyle~{}s_{oms}=8\times
10^{-12}~{}\textrm{m},s_{acc}=3\times 10^{-15}~{}\textrm{m}/\textrm{s}^{2},$
TianQin $\displaystyle:$ $\displaystyle~{}s_{oms}=1\times
10^{-12}~{}\textrm{m},s_{acc}=1\times 10^{-15}~{}\textrm{m}/\textrm{s}^{2}.$
(54)

We adopt the same frequency dependent factors in Eqs. (52, 53) for a
straightforward comparison of the sensitivity of the three detectors.

Since we perform TDI algorithm with the one-link output, the remaining noises
in the output of TDI are also the linear combinations of the time shifted
noises in one-link output. We derived the PSD for the TDI-1.5 combinations
(PSD of other TDI combinations can be found in Appendix B):

$\displaystyle N_{X}$ $\displaystyle=$ $\displaystyle
16\sin^{2}(\tau)\left\\{\left[3+\cos(2\tau)\right]S_{acc}+S_{oms}\right\\},$
(55) $\displaystyle N_{A}$ $\displaystyle=$ $\displaystyle
N_{E}=8\sin^{2}(\tau)\;\left\\{2\left[3+2\cos\tau+\cos(2\tau)\right]\;S_{acc}+\left(2+\cos\tau\right)\;S_{oms}\right\\},$
(56) $\displaystyle N_{T}$ $\displaystyle=$ $\displaystyle
32\sin^{2}\left(\frac{\tau}{2}\right)\sin^{2}(\tau)\;\left[4\sin^{2}\left(\frac{\tau}{2}\right)\;S_{acc}+S_{oms}\right].$
(57)

With the averaged transfer function and PSD of the TDI combinations at hand,
we can calculate sensitivity curves. In Fig. 5, we plot the sensitivity of
Michelson combination $S_{X}$ of the three detectors to the ULBFs.

The optimal sensitivity $S_{\eta}$ is defined by

$\frac{1}{S_{\eta}}=\frac{1}{S_{A}}+\frac{1}{S_{E}}+\frac{1}{S_{T}}.$ (58)

In Fig. 6 we compare the optimal sensitivity with $X$ combination sensitivity
for Taiji. For gravitational waves, the ratio between $S_{\eta}$ and $S_{X}$
oscillates with an average value $3$ in the high frequency part. However, for
scalar and vector fields, the ratio oscillates above $2$ in the high frequency
part. This difference comes from the relative magnitudes between $S_{A}$,
$S_{E}$ and $S_{T}$ in the high frequency band in the two cases. As shown in
Fig. 6, for gravitational waves, although in the long wavelength region the
sensitivity of $A$ and $E$ are much better than $T$, $T$ will have comparable
sensitivity in the high frequency region. For scalar and vector fields, $T$
has much worse sensitivity than $A$ and $E$ across the whole band, thus will
not contribute to $S_{\eta}$ significantly.

Figure 6: The optimal sensitivity $S_{\eta}$ and the ratios of $S_{X}$ to
$S_{\eta}$ for scalar field, vector field and gravitational wave.

##  V Illustrations of Constraints on Dark Matter

Figure 7: The constraints on the coupling constant $d_{g}$ given by the future
space-based gravitational-wave detectors: LISA, Taiji and TianQin. The
observation time $T=1$ year. For comparison, the constraints given by the
fifth-force experiment [68] and the the tests of the equivalence principle (EP
tests) are also shown. The constraints of EP tests shown here are given by the
experiment from the Eöt-Wash group [69] and MICROSCOPE [70]. Regions above the
lines have been excluded. We consider that the coupling constants other than
$d_{g}$ are zero here.

Figure 8: The constraints on the coupling constants $\epsilon_{B}$ and
$\epsilon_{B-L}$ projected by the future space-based gravitational-wave
detectors: LISA, Taiji and TianQin. The observation time $T=1$ year. For
comparison, the constraints given by the tests of the equivalence principle
(EP tests) are also shown. The constraints of EP tests shown here are given by
the experiment from the Eöt-Wash group [71, 69] and MICROSCOPE [72]. Regions
above the lines are excluded.

In this section, we estimate the future constraints on dark matter (DM) models
based on the results in the previous sections.

For scalar field, we consider the following interactions as an illustration
[73, 74],

$\mathcal{L}_{\phi-
SM}=\kappa\phi\left[-\frac{d_{g}\beta_{3}}{2g_{3}}F^{A}_{\mu\nu}F^{A\mu\nu}-\sum_{i=u,d}(d_{m_{i}}+\gamma_{m_{i}}d_{g})m_{i}\bar{\psi}_{i}\psi_{i}\right],$
(59)

where $\kappa\equiv\sqrt{4\pi}/M_{P}$; $F^{A}_{\mu\nu}$ are the gluon field
strength tensors; $g_{3}$ is the $SU(3)$ gauge coupling; and $\beta_{3}$
refers to the QCD beta functions; $\gamma_{m_{i}}$ are the anomalous
dimensions of the $u$ and $d$ quarks; $d_{g}$ is the couplings to the gluonic
field term, and $d_{m_{u}}$ and $d_{m_{d}}$ are the couplings to the quark
mass terms. It is more convenient to define

$d_{\hat{m}}\equiv\frac{d_{m_{d}}m_{d}+d_{m_{u}}m_{u}}{m_{d}+m_{u}},\qquad
d_{\delta m}\equiv\frac{d_{m_{d}}m_{d}-d_{m_{u}}m_{u}}{m_{d}-m_{u}}.$ (60)

Since the nucleon mass is mostly determined by the QCD energy scale,
$\alpha(\phi)$ in Eq. (3) is approximately given by

$\alpha(\phi)\approx d_{g}^{*}\approx d_{g}+0.093(d_{\hat{m}}-d_{g}).$ (61)

Considering that scalar or vector fields make up all the dark matter
components, we have

$\phi_{\vec{k}}=\frac{\sqrt{2\rho_{\textrm{DM}}}}{m_{\phi}},\qquad|\vec{A}|=\frac{\sqrt{2\rho_{\textrm{DM}}}}{m_{A}},\qquad
v\sim v_{\textrm{DM}}\sim 10^{-3},$ (62)

where $\rho_{\textrm{DM}}=0.3$ GeV$/$cm3 is an estimate of the local density
of DM, and $v_{\textrm{DM}}$ is the velocity dispersion of the dark matter in
our Galaxy.

The minimum detectable Fourier amplitude is

$\left|\tilde{h}(m_{\phi})\right|=m_{\phi}\mathcal{M}_{s}\sqrt{T}=\frac{d_{g}^{*}\kappa\phi_{\vec{k}}|\vec{k}|\sqrt{T}}{m_{\phi}}=\sqrt{S_{O}},$
(63)

where $T$ refers to the duration of observation time. We assume that the
duration is shorter than the coherence time of the field $~{}4\times
10^{-16}\textrm{eV}/m_{\phi}\times 10^{7}$ s, so that the approximation of
monochromatic plane wave is satisfied. Moreover, local DM should be treated as
a superposition of plane waves with different propagation directions and
velocities and a $\mathcal{O}(3)$ factor has been found in the estimation of
laboratory experimental searches [75]. Here, we just use the monochromatic
plane wave for illustration and more dedicated study on the stochastic DM will
be conducted in the coming work.

Using the standard sensitivity curve of $X$ combination, we express the
constrain $d_{g}^{*}$ with respect to $m_{\phi}$,

$d_{g}^{*}=\frac{m_{\phi}}{\kappa
v_{\textrm{DM}}}\sqrt{\frac{S_{X}}{2\rho_{\textrm{DM}}T}}.$ (64)

The constraint curves of coupling constant $d_{g}^{*}(m_{\phi})$ given by
space-based detectors are plotted in Fig. 7, assuming $d_{m}=0$ and
$d_{g}^{*}\approx 0.9d_{g}$. Relevant constraints from other experiments,
including the fifth-force experiment [68] and the tests of equivalence
principle (EP tests) [69, 70], are also shown. All detectors will provide more
stringent constraints than the fifth-force experiments in some parameter
space. The constraint given by TianQin can reach $9.5\times 10^{-6}$, which
approaches the lower limit of the fifth-force experiment. Compared to the EP
tests, LISA and Taiji can reach $5.38\times 10^{-7}$ and $2.07\times 10^{-7}$
at $4.28\times 10^{-17}$ eV and $4.44\times 10^{-17}$ eV, respectively, and
the improvement factors are 6 and 16, respectively. In the low mass region,
Eq. (64) has the approximate expression,
$d_{g}^{*}\propto\tau^{-2}\sqrt{S_{acc}}$, i.e. the sensitivity is inversely
proportional to the square of the arm length. Therefore, LISA and Taiji will
have much better sensitivity of $d_{g}^{*}$ compared to TianQin at low-mass
region.

Similarly, the Fourier amplitude of the vector field $\tilde{h}(\omega)$ is
given by

$|\tilde{h}(m_{A})|=m_{A}\mathcal{M}_{v}\sqrt{T}=\frac{\epsilon_{D}eq_{D,j}|\vec{A}|\sqrt{T}}{M_{j}}=\sqrt{S_{O}}.$
(65)

Using the sensitivity curve of $X$ combination for vector fields, we have

$\epsilon_{D}=\frac{m_{A}M_{j}}{eq_{D,j}}\sqrt{\frac{{S}_{X}}{2\rho_{\textrm{DM}}T}}.$
(66)

For a $U(1)_{B}$ gauge boson, $q_{B}/M=(A/\mu)/m_{n}\approx 1/m_{n}$ is almost
identical among different materials, where $A$ is the mass number, $\mu$ is
the atomic mass in atomic units and $m_{n}$ refers to the mass of neutron. But
for a $U(1)_{B-L}$ gauge boson, $q_{B-L}/M=[(A-Z)/\mu]/m_{n}\approx
0.5/m_{n}$, and the neutron ratio $(A-Z)/\mu$ would be different between
materials. Here we use $A/\mu=1$ for a $U(1)_{B}$ gauge boson and
$(A-Z)/\mu=0.5$ for a $U(1)_{B-L}$ gauge boson.

The constraints on coupling constants $|\epsilon_{B}|(m)$ and
$|\epsilon_{B-L}|(m)$ given by space-based detectors are plotted in Fig. 8. As
a comparison, the results of Eöt-Wash experiment and MICROSCOPE are also shown
111The results of MICROSCOPE shown here are from the first results of the
experiment [70]. According to the latest results [76], the constraints would
be improved by a factor of $\sim 2$. Private communication with P. Fayet..
Compared to EP test, gravitational-wave detectors can provide better
constraints on coupling constants $|\epsilon_{B}|$ and $|\epsilon_{B-L}|$ in
some parameter space. For $|\epsilon_{B}|$, the best constraints projected by
LISA, Taiji and TianQin in the sensitive range are improved by $\sim 12000$,
$\sim 31000$ and $\sim 680$ respectively. And for $|\epsilon_{B-L}|$, the
improvement factors are $\sim 1000$, $\sim 2600$ and $\sim 60$. As we can see,
LISA and Taiji will have much better sensitivity of $|\epsilon_{B}|$ and
$|\epsilon_{B-L}|$ compared to TianQin at low-mass region because of the
longer arm length.

##  VI Summary and Conclusions

Ultralight bosonic fields (ULBFs) are well-motivated and predicted in many
physical and cosmological theories. When having tiny couplings to standard
model particles, they may be probed by precision gravitational-wave
interferometers. We investigate the sensitivity of space-based gravitational-
waves detectors to ULBFs systematically. We obtain the sky and polarization
averaged transfer functions of various time-delay interferometry (TDI)
combinations and discuss their asymptotic behaviors in the long wavelength
region. We find that the velocity of field will affect the behaviors of
transfer function, more significantly at high mass region. For TDI
combinations like Michelson, Sagnac and optimal $A$, $E$, there exists a
critical frequency $f_{c}=v/2\pi L$ determined by the field velocity $v$ and
arm length $L$, below which transfer functions $R^{s,v}\propto f^{4}$ and are
similar to gravitational waves; however, for $f_{c}<f<0.1$ Hz, $R^{s,v}$ are
proportional to $f^{6}$. For combinations like $\zeta$ and $T$, although there
is no critical frequency, the magnitudes of transfer functions for scalar and
vector fields are much smaller than that for gravitational waves. The physical
reason is that, for the same frequency these non-relativistic waves will have
much longer wavelength compared to gravitational waves, and $\zeta$ and $T$
are not sensitive to the common displacements of the spacecrafts due to their
symmetric structure.

A direct application of our result is to identify ULBFs as dark matter (DM).
We estimate the ability of detectors to search DM by calculating the
sensitivity curves of future space-based interferometers, which can improve
the current limits by several orders of magnitude, see Figs. 7 and 8. We
should note that it is more realistic and appropriate to treat the local DM as
a stochastic waves background rather than monochromatic waves [75, 77, 78,
79]. Therefore, the precise sensitivity of detectors to the local DM need a
more dedicated study. We will leave the stochastic effect of the fields to
future work.

## acknowledgement

Y.T is supported by National Key Research and Development Program of China
(Grant No.2021YFC2201901), and National Natural Science Foundation of China
under Grants No.12147103 and 11851302. Y.L.W. is supported by the National Key
Research and Development Program of China under Grant No.2020YFC2201501, and
NSFC under Grants No. 11690022, No. 11747601, No. 12147103, and the Strategic
Priority Research Program of the Chinese Academy of Sciences under Grant No.
XDB23030100.

##  Appendix A Analytic expressions of the transfer functions

Explicit forms of some functions appearing in the transfer functions of
Michelson and Sagnac combinations are listed below.

$\displaystyle A_{1}=$
$\displaystyle(\mu_{12}-\mu_{13})\left[1+\cos(2\tau)\right]-2[\mu_{12}\cos(\tau+\vec{k}\cdot\vec{x}_{2})-\mu_{13}\cos(\tau+\vec{k}\cdot\vec{x}_{3})],$
$\displaystyle A_{2}=$
$\displaystyle-(\mu_{12}-\mu_{13})\sin(2\tau)+2[\mu_{12}\sin(\tau+\vec{k}\cdot\vec{x}_{2})-\mu_{13}\sin(\tau+\vec{k}\cdot\vec{x}_{3})],$
$\displaystyle B_{1}=$ $\displaystyle\mu_{21}[1-\cos
3\tau-\cos(\vec{k}\cdot\vec{x}_{2}+\tau)-\cos(\vec{k}\cdot\vec{x}_{2}+2\tau)]+\mu_{13}[1-\cos
3\tau-\cos(\vec{k}\cdot\vec{x}_{3}+\tau)-$
$\displaystyle\cos(\vec{k}\cdot\vec{x}_{3}+2\tau)]+\mu_{32}[\cos(\vec{k}\cdot\vec{x}_{2}+\tau)+\cos(\vec{k}\cdot\vec{x}_{3}+\tau)-\cos(\vec{k}\cdot\vec{x}_{2}+2\tau)-\cos(\vec{k}\cdot\vec{x}_{3}+2\tau)],$
$\displaystyle B_{2}=$ $\displaystyle\mu_{21}[\sin
3\tau+\sin(\vec{k}\cdot\vec{x}_{2}+\tau)+\sin(\vec{k}\cdot\vec{x}_{2}+2\tau)]+\mu_{13}[\sin
3\tau+\sin(\vec{k}\cdot\vec{x}_{3}+\tau)+\sin(\vec{k}\cdot\vec{x}_{3}+2\tau)]$
$\displaystyle+\mu_{32}[-\sin(\vec{k}\cdot\vec{x}_{2}+\tau)-\sin(\vec{k}\cdot\vec{x}_{3}+\tau)+\sin(\vec{k}\cdot\vec{x}_{2}+2\tau)+\sin(\vec{k}\cdot\vec{x}_{3}+2\tau)],$
$\displaystyle C_{1}=$
$\displaystyle\mu_{21}(1+\cos(\vec{k}\cdot\vec{x}_{2}))+\mu_{13}(1+\cos(\vec{k}\cdot\vec{x}_{3}))+\mu_{32}[\cos(\vec{k}\cdot\vec{x}_{2})+\cos(\vec{k}\cdot\vec{x}_{3})],$
$\displaystyle C_{2}=$
$\displaystyle\mu_{21}\sin(\vec{k}\cdot\vec{x}_{2})+\mu_{13}\sin(\vec{k}\cdot\vec{x}_{3})+\mu_{32}[\sin(\vec{k}\cdot\vec{x}_{2})+\sin(\vec{k}\cdot\vec{x}_{3})].$

The transfer function $R_{s}$, $R_{v}$ for various optimal combinations can be
expressed as

$\displaystyle R_{A}^{s}(\omega)$ $\displaystyle=$ $\displaystyle I_{s}(1-\cos
2\tau)[D_{1}^{2}(\theta_{1},\epsilon,\omega)+D_{2}^{2}(\theta_{1},\epsilon,\omega)],$
(67) $\displaystyle R_{A}^{v}(\omega)$ $\displaystyle=$ $\displaystyle
I_{v}(1-\cos
2\tau)[D_{1}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)+D_{2}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)],$
(68) $\displaystyle R_{E}^{s}(\omega)$ $\displaystyle=$ $\displaystyle
I_{s}\frac{1}{3}(1-\cos
2\tau)[D_{3}^{2}(\theta_{1},\epsilon,\omega)+D_{4}^{2}(\theta_{1},\epsilon,\omega)],$
(69) $\displaystyle R_{E}^{v}(\omega)$ $\displaystyle=$ $\displaystyle
I_{v}\frac{1}{3}(1-\cos
2\tau)[D_{3}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)+D_{4}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)],$
(70) $\displaystyle R_{T}^{s}(\omega)$ $\displaystyle=$ $\displaystyle
I_{s}\frac{2}{3}(1-\cos
2\tau)[D_{5}^{2}(\theta_{1},\epsilon,\omega)+D_{6}^{2}(\theta_{1},\epsilon,\omega)],$
(71) $\displaystyle R_{T}^{v}(\omega)$ $\displaystyle=$ $\displaystyle
I_{v}\frac{2}{3}(1-\cos
2\tau)[D_{5}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)+D_{6}^{2}(\theta_{1},\theta_{2},\epsilon_{1},\epsilon_{2},\omega)],$
(72)

where we have defined the following functions,

$\displaystyle D_{1}=$
$\displaystyle-(\mu_{32}+\mu_{12})\left[1+\cos(2\tau)\right]+2\mu_{32}\cos(\tau+\vec{k}\cdot\vec{x}_{2}+\vec{k}\cdot\vec{x}_{3})+2\mu_{12}\cos(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{2}),$
$\displaystyle D_{2}=$
$\displaystyle(\mu_{32}+\mu_{12})\sin(2\tau)-2\mu_{32}\sin(\tau+\vec{k}\cdot\vec{x}_{2}+\vec{k}\cdot\vec{x}_{3})-2\mu_{12}\sin(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{2}),$
$\displaystyle D_{3}=$
$\displaystyle-(2\mu_{13}+\mu_{23}-\mu_{21})\left[1+\cos(2\tau)\right]+4\mu_{13}\cos(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{3})$
$\displaystyle{}+2\mu_{23}\cos(\tau+\vec{k}\cdot\vec{x}_{2}+\vec{k}\cdot\vec{x}_{3})-2\mu_{21}\cos(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{2}),$
$\displaystyle D_{4}=$
$\displaystyle(2\mu_{13}+\mu_{23}-\mu_{21})\sin(2\tau)-4\mu_{13}\sin(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{3})$
$\displaystyle{}-2\mu_{23}\sin(\tau+\vec{k}\cdot\vec{x}_{2}+\vec{k}\cdot\vec{x}_{3})+2\mu_{21}\sin(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{2}),$
$\displaystyle D_{5}=$ $\displaystyle
2(\mu_{12}-\mu_{13}+\mu_{23})\left[1+\cos(2\tau)\right]-4\mu_{12}\cos(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{2})$
$\displaystyle+4\mu_{13}\cos(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{3})-4\mu_{23}\cos(\tau+\vec{k}\cdot\vec{x}_{2}+\vec{k}\cdot\vec{x}_{3}),$
$\displaystyle D_{6}=$
$\displaystyle-2(\mu_{12}-\mu_{13}+\mu_{23})\sin(2\tau)+4\mu_{12}\sin(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{2})$
$\displaystyle{}-4\mu_{13}\sin(\tau+\vec{k}\cdot\vec{x}_{1}+\vec{k}\cdot\vec{x}_{3})+4\mu_{23}\sin(\tau+\vec{k}\cdot\vec{x}_{2}+\vec{k}\cdot\vec{x}_{3}).$

##  Appendix B Power spectral density of TDI combinations

The noise power spectral density of different TDI channels are calculated and
summarized below. Michelson X, Y, Z:

$N_{X}=16\sin^{2}(\tau)\left\\{\left[3+\cos(2\tau)\right]S_{acc}+S_{oms}\right\\}.$
(73)

Fully symmetric Sagnac $\zeta$:

$N_{\zeta}=24\sin^{2}\left(\frac{\tau}{2}\right)\;S_{acc}+6S_{oms}.$ (74)

Sagnac $\alpha$, $\beta$, $\gamma$:

$N_{\alpha}=N_{\beta}=N_{\gamma}=\left[16\sin^{2}\left(\frac{\tau}{2}\right)+8\sin^{2}\left(\frac{3\tau}{2}\right)\right]S_{acc}+6S_{oms}.$
(75)

Beacon P, Q, R:

$N_{P}=N_{Q}=N_{R}=16\left[2\sin^{4}\left(\frac{\tau}{2}\right)+\sin^{2}(\tau)\right]S_{acc}+8\left[\sin^{2}\left(\frac{\tau}{2}\right)+\sin^{2}(\tau)\right]S_{oms}.$
(76)

Monitor E, F, G:

$N_{E}=N_{F}=N_{G}=16\left[2\sin^{4}\left(\frac{\tau}{2}\right)+\sin^{2}(\tau)\right]S_{acc}+8\left[\sin^{2}\left(\frac{\tau}{2}\right)+\sin^{2}(\tau)\right]S_{oms}.$
(77)

Relay U, V, W:

$\displaystyle N_{U}=N_{V}=N_{W}$ $\displaystyle=$ $\displaystyle
8\left[2\sin^{2}\left(\frac{\tau}{2}\right)+\sin^{2}(\tau)+2\sin^{2}\left(\frac{3\tau}{2}\right)\right]S_{acc}$
(78)
$\displaystyle{}+4\left[\sin^{2}\left(\frac{\tau}{2}\right)+2\sin^{2}(\tau)+\sin^{2}\left(\frac{3\tau}{2}\right)\right]S_{oms}.$

Optimal A, E, T:

$\displaystyle N_{A}$ $\displaystyle=$ $\displaystyle
N_{E}=8\sin^{2}(\tau)\;\left\\{2\left[3+2\cos\tau+\cos(2\tau)\right]\;S_{acc}+\left(2+\cos\tau\right)\;S_{oms}\right\\},$
(79) $\displaystyle N_{T}$ $\displaystyle=$ $\displaystyle
32\sin^{2}\left(\frac{\tau}{2}\right)\sin^{2}(\tau)\;\left[4\sin^{2}\left(\frac{\tau}{2}\right)\;S_{acc}+S_{oms}\right].$
(80)

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Sathyaprakash and Schutz [2009] B. S. Sathyaprakash and B. F. Schutz, Physics, Astrophysics and Cosmology with Gravitational Waves, Living Rev. Rel. 12, 2 (2009), arXiv:0903.0338 [gr-qc] .
  * Pierce  _et al._ [2018] A. Pierce, K. Riles, and Y. Zhao, Searching for dark photon dark matter with gravitational-wave detectors, Phys. Rev. Lett. 121, 061102 (2018).
  * Morisaki and Suyama [2019] S. Morisaki and T. Suyama, Detectability of ultralight scalar field dark matter with gravitational-wave detectors, Phys. Rev. D 100, 123512 (2019).
  * Michimura  _et al._ [2020a] Y. Michimura, T. Fujita, S. Morisaki, H. Nakatsuka, and I. Obata, Ultralight vector dark matter search with auxiliary length channels of gravitational wave detectors, Phys. Rev. D 102, 102001 (2020a), arXiv:2008.02482 [hep-ph] .
  * Morisaki  _et al._ [2021] S. Morisaki, T. Fujita, Y. Michimura, H. Nakatsuka, and I. Obata, Improved sensitivity of interferometric gravitational wave detectors to ultralight vector dark matter from the finite light-traveling time, Phys. Rev. D 103, L051702 (2021), arXiv:2011.03589 [hep-ph] .
  * Kawasaki and Nakayama [2013] M. Kawasaki and K. Nakayama, Axions: Theory and Cosmological Role, Ann. Rev. Nucl. Part. Sci. 63, 69 (2013), arXiv:1301.1123 [hep-ph] .
  * Marsh [2016] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Hui  _et al._ [2017] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Ultralight scalars as cosmological dark matter, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Hui [2021] L. Hui, Wave Dark Matter, Ann. Rev. Astron. Astrophys. 59, 247 (2021), arXiv:2101.11735 [astro-ph.CO] .
  * Ferreira [2021] E. G. M. Ferreira, Ultra-light dark matter, Astron. Astrophys. Rev. 29, 7 (2021), arXiv:2005.03254 [astro-ph.CO] .
  * Abbott  _et al._ [2022] R. Abbott _et al._ (LIGO Scientific, KAGRA, Virgo), Constraints on dark photon dark matter using data from LIGO’s and Virgo’s third observing run, Phys. Rev. D 105, 063030 (2022), arXiv:2105.13085 [astro-ph.CO] .
  * Guo  _et al._ [2019] H.-K. Guo, K. Riles, F.-W. Yang, and Y. Zhao, Searching for Dark Photon Dark Matter in LIGO O1 Data, Commun. Phys. 2, 155 (2019), arXiv:1905.04316 [hep-ph] .
  * Fukusumi  _et al._ [2023] K. Fukusumi, S. Morisaki, and T. Suyama, Upper limit on scalar field dark matter from LIGO-Virgo third observation run, (2023), arXiv:2303.13088 [hep-ph] .
  * Miller and Mendes [2023] A. L. Miller and L. Mendes, First search for ultralight dark matter with a space-based gravitational-wave antenna: Lisa pathfinder, Phys. Rev. D 107, 063015 (2023).
  * Adelberger  _et al._ [2003a] E. G. Adelberger, B. R. Heckel, and A. E. Nelson, Tests of the gravitational inverse square law, Ann. Rev. Nucl. Part. Sci. 53, 77 (2003a), arXiv:hep-ph/0307284 .
  * Schlamminger  _et al._ [2008a] S. Schlamminger, K. Y. Choi, T. A. Wagner, J. H. Gundlach, and E. G. Adelberger, Test of the equivalence principle using a rotating torsion balance, Phys. Rev. Lett. 100, 041101 (2008a), arXiv:0712.0607 [gr-qc] .
  * Touboul  _et al._ [2022a] P. Touboul _et al._ (MICROSCOPE), MICROSCOPE Mission: Final Results of the Test of the Equivalence Principle, Phys. Rev. Lett. 129, 121102 (2022a), arXiv:2209.15487 [gr-qc] .
  * Graham  _et al._ [2016a] P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, and W. A. Terrano, Dark Matter Direct Detection with Accelerometers, Phys. Rev. D 93, 075029 (2016a), arXiv:1512.06165 [hep-ph] .
  * Grote and Stadnik [2019] H. Grote and Y. V. Stadnik, Novel signatures of dark matter in laser-interferometric gravitational-wave detectors, Phys. Rev. Res. 1, 033187 (2019), arXiv:1906.06193 [astro-ph.IM] .
  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, $\mathrm{CP}$ conservation in the presence of pseudoparticles, Phys. Rev. Lett. 38, 1440 (1977).
  * Wilczek [1978] F. Wilczek, Problem of strong $p$ and $t$ invariance in the presence of instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Weinberg [1978] S. Weinberg, A new light boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120, 137 (1983).
  * Nelson and Scholtz [2011] A. E. Nelson and J. Scholtz, Dark Light, Dark Matter and the Misalignment Mechanism, Phys. Rev. D 84, 103501 (2011), arXiv:1105.2812 [hep-ph] .
  * Graham  _et al._ [2016b] P. W. Graham, J. Mardon, and S. Rajendran, Vector Dark Matter from Inflationary Fluctuations, Phys. Rev. D 93, 103520 (2016b), arXiv:1504.02102 [hep-ph] .
  * Ema  _et al._ [2019] Y. Ema, K. Nakayama, and Y. Tang, Production of purely gravitational dark matter: the case of fermion and vector boson, JHEP 07, 060, arXiv:1903.10973 [hep-ph] .
  * Ahmed  _et al._ [2020] A. Ahmed, B. Grzadkowski, and A. Socha, Gravitational production of vector dark matter, JHEP 08, 059, arXiv:2005.01766 [hep-ph] .
  * Agrawal  _et al._ [2020] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, Relic Abundance of Dark Photon Dark Matter, Phys. Lett. B 801, 135136 (2020), arXiv:1810.07188 [hep-ph] .
  * Tang and Wu [2020] Y. Tang and Y.-L. Wu, Weyl scaling invariant $R^{2}$ gravity for inflation and dark matter, Phys. Lett. B 809, 135716 (2020), arXiv:2006.02811 [hep-ph] .
  * Kitajima and Nakayama [2023] N. Kitajima and K. Nakayama, Viable vector coherent oscillation dark matter, JCAP 07, 014, arXiv:2303.04287 [hep-ph] .
  * Kolb and Long [2021] E. W. Kolb and A. J. Long, Completely dark photons from gravitational particle production during the inflationary era, JHEP 03, 283, arXiv:2009.03828 [astro-ph.CO] .
  * An  _et al._ [2021] H. An, F. P. Huang, J. Liu, and W. Xue, Radio-frequency Dark Photon Dark Matter across the Sun, Phys. Rev. Lett. 126, 181102 (2021), arXiv:2010.15836 [hep-ph] .
  * Moroi and Yin [2021] T. Moroi and W. Yin, Light Dark Matter from Inflaton Decay, JHEP 03, 301, arXiv:2011.09475 [hep-ph] .
  * Liang  _et al._ [2023] D. Liang, R. Xu, Z.-F. Mai, and L. Shao, Probing vector hair of black holes with extreme-mass-ratio inspirals, Phys. Rev. D 107, 044053 (2023), arXiv:2212.09346 [gr-qc] .
  * Sun  _et al._ [2022] S. Sun, X.-Y. Yang, and Y.-L. Zhang, Pulsar timing residual induced by wideband ultralight dark matter with spin 0,1,2, Phys. Rev. D 106, 066006 (2022), arXiv:2112.15593 [astro-ph.CO] .
  * An  _et al._ [2023] H. An, S. Ge, W.-Q. Guo, X. Huang, J. Liu, and Z. Lu, Direct Detection of Dark Photon Dark Matter Using Radio Telescopes, Phys. Rev. Lett. 130, 181001 (2023), arXiv:2207.05767 [hep-ph] .
  * Luu  _et al._ [2023] H. N. Luu, T. Liu, J. Ren, T. Broadhurst, R. Yang, J.-S. Wang, and Z. Xie, Detecting Stochastic Wave Dark Matter with Fermi-LAT $\gamma$-ray Pulsar Timing Array, (2023), arXiv:2304.04735 [astro-ph.HE] .
  * Liu  _et al._ [2023] T. Liu, X. Lou, and J. Ren, Pulsar Polarization Arrays, Phys. Rev. Lett. 130, 121401 (2023), arXiv:2111.10615 [astro-ph.HE] .
  * Chen  _et al._ [2020] Y. Chen, J. Shu, X. Xue, Q. Yuan, and Y. Zhao, Probing Axions with Event Horizon Telescope Polarimetric Measurements, Phys. Rev. Lett. 124, 061102 (2020), arXiv:1905.02213 [hep-ph] .
  * Tang  _et al._ [2023] Z. Tang _et al._ , SRF Cavity Searches for Dark Photon Dark Matter: First Scan Results, (2023), arXiv:2305.09711 [hep-ex] .
  * Xia  _et al._ [2021] C. Xia, Y.-H. Xu, and Y.-F. Zhou, Constraining light dark matter upscattered by ultrahigh-energy cosmic rays, Nucl. Phys. B 969, 115470 (2021), arXiv:2009.00353 [hep-ph] .
  * Fayet [1980] P. Fayet, Effects of the Spin 1 Partner of the Goldstino (Gravitino) on Neutral Current Phenomenology, Phys. Lett. B 95, 285 (1980).
  * Fayet [1989] P. Fayet, The Fifth Force Charge as a Linear Combination of Baryonic, Leptonic (Or $B^{-}$l) and Electric Charges, Phys. Lett. B 227, 127 (1989).
  * Fayet [1990] P. Fayet, Extra U(1)’s and New Forces, Nucl. Phys. B 347, 743 (1990).
  * Fayet [2018] P. Fayet, MICROSCOPE limits for new long-range forces and implications for unified theories, Phys. Rev. D 97, 055039 (2018), arXiv:1712.00856 [hep-ph] .
  * Wu  _et al._ [2023] Y.-M. Wu, Z.-C. Chen, and Q.-G. Huang, Pulsar timing residual induced by ultralight tensor dark matter, (2023), arXiv:2305.08091 [hep-ph] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Evans  _et al._ [2021] M. Evans _et al._ , A Horizon Study for Cosmic Explorer: Science, Observatories, and Community, (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The Einstein Telescope: A third-generation gravitational wave observatory, Class. Quant. Grav. 27, 194002 (2010).
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Hu and Wu [2017] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity, Natl. Sci. Rev. 4, 685 (2017).
  * Luo  _et al._ [2016] J. Luo _et al._ (TianQin), TianQin: a space-borne gravitational wave detector, Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Estabrook  _et al._ [2000] F. B. Estabrook, M. Tinto, and J. W. Armstrong, Time delay analysis of LISA gravitational wave data: Elimination of spacecraft motion effects, Phys. Rev. D 62, 042002 (2000).
  * Tinto and Dhurandhar [2005] M. Tinto and S. V. Dhurandhar, TIME DELAY, Living Rev. Rel. 8, 4 (2005), arXiv:gr-qc/0409034 .
  * Wu [2016] Y.-L. Wu, Quantum field theory of gravity with spin and scaling gauge invariance and spacetime dynamics with quantum inflation, Phys. Rev. D 93, 024012 (2016), arXiv:1506.01807 [hep-th] .
  * Wu [2018] Y.-L. Wu, Hyperunified field theory and gravitational gauge–geometry duality, Eur. Phys. J. C 78, 28 (2018), arXiv:1712.04537 [hep-th] .
  * Wu [2023] Y.-L. Wu, Gravidynamics, spinodynamics and electrodynamics within the framework of gravitational quantum field theory, Sci. China Phys. Mech. Astron. 66, 260411 (2023), arXiv:2208.03290 [hep-th] .
  * Wang  _et al._ [2022] Q.-Y. Wang, Y. Tang, and Y.-L. Wu, Dark matter production in Weyl R2 inflation, Phys. Rev. D 106, 023502 (2022), arXiv:2203.15452 [hep-ph] .
  * Michimura  _et al._ [2020b] Y. Michimura, T. Fujita, S. Morisaki, H. Nakatsuka, and I. Obata, Ultralight vector dark matter search with auxiliary length channels of gravitational wave detectors, Phys. Rev. D 102, 102001 (2020b).
  * Tinto and Dhurandhar [2021] M. Tinto and S. V. Dhurandhar, Time-delay interferometry, Living Rev. Rel. 24, 1 (2021).
  * Otto [2015] M. Otto,  _Time-Delay Interferometry Simulations for the Laser Interferometer Space Antenna_ , Ph.D. thesis, Leibniz U., Hannover (2015).
  * Prince  _et al._ [2002] T. A. Prince, M. Tinto, S. L. Larson, and J. W. Armstrong, The LISA optimal sensitivity, Phys. Rev. D 66, 122002 (2002), arXiv:gr-qc/0209039 .
  * Babak  _et al._ [2021] S. Babak, M. Hewitson, and A. Petiteau, Lisa sensitivity and snr calculations (2021), arXiv:2108.01167 [astro-ph.IM] .
  * Adelberger  _et al._ [2003b] E. Adelberger, B. Heckel, and A. Nelson, Tests of the gravitational inverse-square law, Annual Review of Nuclear and Particle Science 53, 77 (2003b), https://doi.org/10.1146/annurev.nucl.53.041002.110503 .
  * Wagner  _et al._ [2012] T. A. Wagner, S. Schlamminger, J. H. Gundlach, and E. G. Adelberger, Torsion-balance tests of the weak equivalence principle, Classical and Quantum Gravity 29, 184002 (2012).
  * Touboul  _et al._ [2017] P. Touboul, G. Métris, M. Rodrigues, Y. André, Q. Baghi, J. Bergé, D. Boulanger, S. Bremer, P. Carle, R. Chhun, B. Christophe, V. Cipolla, T. Damour, P. Danto, H. Dittus, P. Fayet, B. Foulon, C. Gageant, P.-Y. Guidotti, D. Hagedorn, E. Hardy, P.-A. Huynh, H. Inchauspe, P. Kayser, S. Lala, C. Lämmerzahl, V. Lebat, P. Leseur, F. m. c. Liorzou, M. List, F. Löffler, I. Panet, B. Pouilloux, P. Prieur, A. Rebray, S. Reynaud, B. Rievers, A. Robert, H. Selig, L. Serron, T. Sumner, N. Tanguy, and P. Visser, Microscope mission: First results of a space test of the equivalence principle, Phys. Rev. Lett. 119, 231101 (2017).
  * Schlamminger  _et al._ [2008b] S. Schlamminger, K.-Y. Choi, T. A. Wagner, J. H. Gundlach, and E. G. Adelberger, Test of the equivalence principle using a rotating torsion balance, Phys. Rev. Lett. 100, 041101 (2008b).
  * Bergé  _et al._ [2018] J. Bergé, P. Brax, G. Métris, M. Pernot-Borràs, P. Touboul, and J.-P. Uzan, Microscope mission: First constraints on the violation of the weak equivalence principle by a light scalar dilaton, Phys. Rev. Lett. 120, 141101 (2018).
  * Damour and Donoghue [2010] T. Damour and J. F. Donoghue, Equivalence principle violations and couplings of a light dilaton, Phys. Rev. D 82, 084033 (2010).
  * Piazza and Pospelov [2010] F. Piazza and M. Pospelov, Sub-eV scalar dark matter through the super-renormalizable Higgs portal, Phys. Rev. D 82, 043533 (2010), arXiv:1003.2313 [hep-ph] .
  * Centers  _et al._ [2021] G. P. Centers _et al._ , Stochastic fluctuations of bosonic dark matter, Nature Commun. 12, 7321 (2021), arXiv:1905.13650 [astro-ph.CO] .
  * Touboul  _et al._ [2022b] P. Touboul, G. Métris, M. Rodrigues, J. Bergé, A. Robert, Q. Baghi, Y. André, J. Bedouet, D. Boulanger, S. Bremer, P. Carle, R. Chhun, B. Christophe, V. Cipolla, T. Damour, P. Danto, L. Demange, H. Dittus, O. Dhuicque, P. Fayet, B. Foulon, P.-Y. Guidotti, D. Hagedorn, E. Hardy, P.-A. Huynh, P. Kayser, S. Lala, C. Lämmerzahl, V. Lebat, F. m. c. Liorzou, M. List, F. Löffler, I. Panet, M. Pernot-Borràs, L. Perraud, S. Pires, B. Pouilloux, P. Prieur, A. Rebray, S. Reynaud, B. Rievers, H. Selig, L. Serron, T. Sumner, N. Tanguy, P. Torresi, and P. Visser (MICROSCOPE Collaboration), $microscope$ mission: Final results of the test of the equivalence principle, Phys. Rev. Lett. 129, 121102 (2022b).
  * Derevianko [2018] A. Derevianko, Detecting dark-matter waves with a network of precision-measurement tools, Phys. Rev. A 97, 042506 (2018).
  * Foster  _et al._ [2018] J. W. Foster, N. L. Rodd, and B. R. Safdi, Revealing the dark matter halo with axion direct detection, Phys. Rev. D 97, 123006 (2018).
  * Kim [2023] H. Kim, Gravitational Interaction of Ultralight Dark Matter with Interferometers, (2023), arXiv:2306.13348 [hep-ph] .