# Pulsar timing residual induced by ultralight tensor dark matter

Yu-Mei Wu<EMAIL_ADDRESS>School of Fundamental Physics and Mathematical
Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China
School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China  CAS Key Laboratory of Theoretical
Physics, Institute of Theoretical Physics, Chinese Academy of Sciences,
Beijing 100190, China  Zu-Cheng Chen  Corresponding author:
<EMAIL_ADDRESS>Department of Astronomy, Beijing Normal University,
Beijing 100875, China  Advanced Institute of Natural Sciences, Beijing Normal
University, Zhuhai 519087, China  Department of Physics and Synergistic
Innovation Center for Quantum Effects and Applications, Hunan Normal
University, Changsha, Hunan 410081, China  Qing-Guo Huang  Corresponding
author<EMAIL_ADDRESS>School of Fundamental Physics and Mathematical
Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China
School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China  CAS Key Laboratory of Theoretical
Physics, Institute of Theoretical Physics, Chinese Academy of Sciences,
Beijing 100190, China

###### Abstract

Ultralight boson fields, with a mass around $10^{-23}\text{eV}$, are promising
candidates for the elusive cosmological dark matter. These fields induce a
periodic oscillation of the spacetime metric in the nanohertz frequency band,
which is detectable by pulsar timing arrays. In this paper, we investigate the
gravitational effect of ultralight tensor dark matter on the arrival time of
radio pulses from pulsars. We find that the pulsar timing signal caused by
tensor dark matter exhibits a different angular dependence than that by scalar
and vector dark matter, making it possible to distinguish the ultralight dark
matter signal with different spins. Combining the gravitational effect and the
coupling effect of ultralight tensor dark matter with standard model matter
provides a complementary way to constrain the coupling parameter $\alpha$. We
estimate $\alpha\lesssim 10^{-6}\sim 10^{-5}$ in the mass range $m<5\times
10^{-23}\mathrm{eV}$ with current pulsar timing array.

##  I Introduction

Numerous cosmological observations (Rubin  _et al._ , 1980, 1982; Faber and
Jackson, 1976; Massey  _et al._ , 2010) have provided overwhelming support for
the existence of dark matter, and precise analyses of cosmic microwave
background (CMB) (Aghanim  _et al._ , 2020) in particular have estimated the
fraction of dark matter in the total energy of the Universe. However, the
nature and properties of dark matter, such as its composition and interaction
with the Standard Model (SM), remain mysterious. There are plenty dark matter
candidates (Bertone  _et al._ , 2005; Feng, 2010; Bertone and Tait, 2018),
from tiny axion-like particles with mass around $10^{-22}\,\mathrm{\rm{eV}}$
(Marsh, 2016) to planet-sized massive compact halo objects (MACHOs) of $m\sim
M_{\odot}$ (Alcock  _et al._ , 2000). Among the bunch candidates, the weakly
interactive massive particles (WIMPs) are expected to be highly desirable.
However, despite considerable efforts in searching for them, no signal of
WIMPs has yet been observed (Schumann, 2019).

In recent years, ultralight dark matter (ULDM) with a mass around
$10^{-23}\,\mathrm{\rm{eV}}$ has attracted wide attention (Hu  _et al._ ,
2000; Hui  _et al._ , 2017). The ULDM is assumed to be condensed before
inflation and oscillates coherently when its mass is larger than the Hubble
scale (Fox  _et al._ , 2004). On the one hand, ULDM moves at a non-
relativistic speed, retaining the advantages of traditional cold dark matter
to explain the large-scale structure of galaxies. On the other hand, its large
de Broglie wavelength can suppress the structure on small scales, providing
feasible solutions to the puzzles encountered by the traditional cold dark
matter (Hui  _et al._ , 2017; Niemeyer, 2019), such as the core-cusp problem
Gentile  _et al._ (2004); de Blok (2010) and missing-satellites problem Moore
_et al._ (1999); Klypin  _et al._ (1999).

ULDM has been traditionally represented by scalar (Hu  _et al._ , 2000;
Khmelnitsky and Rubakov, 2014; Hui  _et al._ , 2017) or pseduo-scalar fields
(spin-0) such as axion-like particles (Arvanitaki  _et al._ , 2010; Hlozek
_et al._ , 2015; Payez  _et al._ , 2015; Ivanov  _et al._ , 2019). However,
fields with higher spins, such as the massive vector fields (spin-1), can also
act as alternative candidates for ULDM. Although ULDM candidates with
different spins share many similarities in the cosmological phenomena, each
possesses specific features. For instance, the massive vector particles can
interact with SM particles if they are gauge bosons for $U(1)_{B}$ (“B” refers
to baryon) or $U(1)_{B-L}$ (“L” refers to lepton) symmetry (Pierce  _et al._ ,
2018), and their multiple polarization modes may cause anisotropies in certain
experiments (Nomura  _et al._ , 2020).

Recently, massive tensor fields (spin-2) have been proposed as another viable
candidate for dark matter (Aoki and Mukohyama, 2016; Babichev  _et al._ ,
2016a, b; Aoki and Maeda, 2018). The tensor ULDM can be naturally predicted by
a ghost-free bimetric theory (Hassan and Rosen, 2012) in some mass range
(Marzola  _et al._ , 2018). This theory has two mass eigenstates at the
linearized level: a massless spin-2 field $G_{\mu\nu}$ and a massive spin-2
field $M_{\mu\nu}$. At nonlinear levels, the expansion of the action (Babichev
_et al._ , 2016a, b; Marzola  _et al._ , 2018; Gialamas and Tamvakis, 2023)
reveals that the self-interaction terms of the massless graviton $G_{\mu\nu}$
exactly match the Einstein-Hilbert structure. Meanwhile, the massive mode
gravitates with the same strength as SM matter. But terms that contain one
$M_{\mu\nu}$ and at least one $G_{\mu\nu}$, such as $MG$, $MG^{2}$, and
$MG^{3}$, are missing in the expansion. As a result, the terms with
$G_{\mu\nu}$ can be partially resumed up to an effective background
$g_{\mu\nu}$. In contrast, the massive fluctuation $M_{\mu\nu}$ can be treated
as a propagating massive spin-2 field on the background $g_{\mu\nu}$.

The ULDM can leave imprints on cosmological and astrophysical observables. Its
mass range and properties have been explored by many experiments, such as CMB
(Hlozek  _et al._ , 2015; Poulin  _et al._ , 2018), Lyman-$\alpha$ (Iršič  _et
al._ , 2017, 2017), non-observation of superradiance (Ünal  _et al._ , 2021),
secular variations in binary-pulsar orbital parameters (Armaleo  _et al._ ,
2020a; Blas  _et al._ , 2020), etc. In particular, for dark matter fields with
a mass around $10^{-23}\,\mathrm{\rm{eV}}$, pulsar timing arrays (PTAs)
provide a powerful tool for exploring both its gravitational and non-
gravitational properties (Khmelnitsky and Rubakov, 2014; Porayko and Postnov,
2014; Porayko  _et al._ , 2018; Kato and Soda, 2020; Nomura  _et al._ , 2020;
Xue  _et al._ , 2022; Wu  _et al._ , 2022; Armaleo  _et al._ , 2020b). A PTA
consists of a set of millisecond pulsars for which the arrival time of the
emitted radio pulses are monitored with high precision for decades (Sazhin,
1978; Detweiler, 1979; Foster and Backer, 1990), recording the spacetime
disturbances in the nanohertz range. The coherent motion of ULDM carrying
energy and momentum induces nontrivial oscillations of the metric (Khmelnitsky
and Rubakov, 2014; Nomura  _et al._ , 2020), which can leave a mark on the
timing residuals of PTA, known as the “gravitational effect”. Additionally, if
ULDM interacts with SM particles, it can lead to pulsar spin fluctuations and
reference clock shifts due to the induced variations in fundamental constants
of the SM (Kaplan  _et al._ , 2022) or can exert periodic forces on both the
pulsar and Earth, resulting in relative displacements between them (Xue  _et
al._ , 2022). These effects are called the “coupling effect” and can cause
timing residuals.

The specific form of the pulsar timing signal depends on the spin of ULDM. One
unique aspect of tensor dark matter is that it naturally couples to SM fields
without manually adding a fifth force to mimic the modification of gravity.
The coupling effect on PTA has been has been investigated in previous studies
(Armaleo  _et al._ , 2020b; Sun  _et al._ , 2022; Unal  _et al._ , 2022). Our
present study instead focuses on providing the first comprehensive analysis of
the gravitational effect of tensor dark matter on PTA. This analysis reveals
distinct differences in angular dependence compared to scalar and vector dark
matter due to the presence of more polarization modes. Furthermore, the
gravitational effect also differs from the coupling effect regarding angular
dependence. The gravitational effect influences the arrival time of pulses
through metric fluctuations induced by the energy-momentum tensor, which is
quadratic in the tensor fields. In contrast, the coupling effect depends on
the tensor fields linearly. The paper is organized as follows. In Section II,
we show the time-dependent spacetime perturbations induced by the oscillating
tensor ULDM. In Section III, we derive the pulsar timing residuals induced by
the tensor ULDM from the gravitational effect. We also discuss the constraints
that can be inferred from current data and compare the gravitational and
coupling effects in IV. In Section V, we summarize the results.

##  II Metric perturbations induced by tensor ULDM

The action for a massive spin-2 filed $M_{\mu\nu}$ propagating on a background
$g_{\mu\nu}$ can be formally written as (Marzola  _et al._ , 2018)

$S=\frac{M_{\rm{P}}^{2}}{2}\int d^{4}x\sqrt{-g}R(g)+\int
d^{4}x\sqrt{-g}{\cal{L}}_{\rm{FP}}^{(2)}(M)+{\cal{O}}(M_{\mu\nu}^{3}),$ (1)

where $M_{\rm{P}}$ is the reduced Planck mass, and
${\cal{L}}_{\rm{FP}}^{(2)}(M)$ is the quadratic Fierz-Pauli Lagrangian,

${\cal{L}}_{\rm{FP}}^{(2)}(M)=\frac{1}{8}\left[-\nabla_{\rho}M_{\mu\nu}\nabla^{\rho}M^{\mu\nu}+2\nabla_{\rho}M_{\mu\nu}\nabla^{\nu}M^{\mu\rho}-2\nabla_{\nu}M^{\mu\nu}\nabla_{\mu}M+\nabla^{\mu}M\nabla_{\mu}M+m^{2}(M^{2}-M_{\mu\nu}M^{\mu\nu})\right]$
(2)

with the trace $M=M_{\,\,\,\mu}^{\mu}$. The cosmic expansion is negligible on
galactic scales, so we treat the background to be flat. The equation of motion
of the tensor fields can be derived from the action Eq. (1) by taking the
variation with respect to $M_{\mu\nu}$,

$\square
M^{\mu\nu}-(\partial^{\mu}\partial_{\rho}M^{\nu\rho}+\partial^{\nu}\partial_{\rho}M^{\mu\rho})+\eta^{\mu\nu}\partial_{\rho\sigma}M^{\rho\sigma}+\partial^{\mu\nu}M-\eta^{\mu\nu}\square
M=m^{2}(M^{\mu\nu}-\eta^{\mu\nu}M).$ (3)

Applying $\partial_{\mu}$ to both sides of Eq. (3), the left-hand side gives
zero identically, so we obtain

$m^{2}\partial_{\mu}(M^{\mu\nu}-\eta^{\mu\nu}M)=0.$ (4)

Further, taking the trace of both sides of Eq. (3), we get

$2\partial_{\mu}\partial_{\nu}(M^{\mu\nu}-\eta^{\mu\nu}M)=-3m^{2}M.$ (5)

With Eq. (4) and Eq. (5), we can rewrite Eq. (3) into a set of equations for a
massive tensor field in the standard form

$\displaystyle M$ $\displaystyle=$ $\displaystyle 0,$ (6)
$\displaystyle\partial_{\mu}M^{\mu\nu}$ $\displaystyle=$ $\displaystyle 0,$
(7) $\displaystyle(\square-m^{2})M^{\mu\nu}$ $\displaystyle=$ $\displaystyle
0.$ (8)

Now we consider the scenario where the spin-2 fields with an ultralight mass
act as the dark matter following the discussion in Refs. (Khmelnitsky and
Rubakov, 2014; Nomura  _et al._ , 2020; Armaleo  _et al._ , 2020b). For ULDM
with the typical mass $m\sim 10^{-23}\rm{eV}$ and the typical velocity in the
galaxy $v\sim 10^{-3}$, the occupation number can be estimated to be as high
as $10^{95}$, allowing us to treat the ultralight field as a classical wave.
The wave is characterized by a typical wave number $k=mv$ and frequency
$\omega\approx m+mv^{2}/2\approx m$ due to the non-relativistic velocity.
Additionally, the de Broglie wavelength of the ULDM can be calculated as

$\lambda_{\rm{dB}}=\frac{2\pi}{k}\approx
4{\rm{kpc}}\left(\frac{10^{-23}{\rm{eV}}}{m}\right)\left(\frac{10^{-3}}{v}\right),$
(9)

indicating the spin-2 field oscillates coherently with a monochromatic
frequency determined by its mass on the scale of $\lambda_{\rm{dB}}$.

We return to Eqs. (6-8). Note that the constraint in Eq. (7) indicates that
$M^{00}$ is suppressed by the order of $k/m=v\sim 10^{-3}$ compared to
$M^{0j}$, while $M^{0j}$ is suppressed by the same order compared to $M^{ij}$.
Thus, we can neglect all the time components $M^{0\mu}$. Since the remains six
spatial components $M_{ij}$ are subject to the traceless constraint in Eq.
(6), the solution to Eq. (8) is given by

$M_{ij}(t,\mathbf{x})=A(\mathbf{x})\cos(mt+\alpha(\mathbf{x}))\epsilon_{ij}.$
(10)

We have neglected the spatial derivatives when solving the equation but retain
the spatial dependence of the amplitude $A(\mathbf{x})$ and the phase
$\alpha(\mathbf{x})$. The symmetric, traceless matrix $\epsilon_{ij}$ with a
unit norm encoding the five free degrees of freedom of the massive spin-2
fields is (Armaleo  _et al._ , 2020a),

$\displaystyle\begin{gathered}\epsilon_{ij}=\frac{1}{\sqrt{2}}\begin{pmatrix}\epsilon_{\rm{T}}\cos\chi-\epsilon_{\rm{S}}/\sqrt{3}&\epsilon_{\rm{T}}\sin\chi&\epsilon_{\rm{V}}\cos\xi\\\
\epsilon_{\rm{T}}\sin\chi&-\epsilon_{\rm{T}}\cos\chi-\epsilon_{\rm{S}}/\sqrt{3}&\epsilon_{\rm{V}}\sin\xi\\\
\epsilon_{\rm{V}}\cos\xi&\epsilon_{\rm{V}}\sin\xi&2\epsilon_{\rm{S}}/\sqrt{3}\end{pmatrix}\end{gathered}.$
(12)

The amplitude parameters $\epsilon_{\rm{S}}$, $\epsilon_{\rm{V}}$, and
$\epsilon_{\rm{T}}$, satisfying
$\epsilon_{\rm{S}}^{2}+\epsilon_{\rm{V}}^{2}+\epsilon_{\rm{T}}^{2}=1$,
represent the ratios of the scalar, vector, and tensor sector of the massive
spin-2 fields, respectively. Meanwhile, the angle parameters $\xi$ and $\chi$
determine the azimuthal direction of the helicity-$\pm 1$ and helicity-$\pm 2$
modes, respectively.

Now let’s focus on the metric perturbation induced by the oscillating fields.
Using the definition of energy-momentum from the matter fields,

$T_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta S}{\delta g^{\mu\nu}},$ (13)

the energy-momentum of the massive spin-2 fields in a flat background is given
by,

$T_{\mu\nu}=-\frac{1}{8}\left[\eta_{\mu\nu}\left(\partial_{\lambda}M^{\gamma\delta}\partial^{\lambda}M_{\gamma\delta}+m^{2}M^{\gamma\delta}M_{\gamma\delta}\right)-2\partial_{\mu}M_{\gamma\delta}\partial_{\nu}M^{\gamma\delta}-4\partial_{\lambda}M_{\mu\delta}\partial^{\lambda}M^{\delta}_{\,\,\nu}-4m^{2}M_{\mu\delta}M^{\delta}_{\,\,\nu}\right].$
(14)

with the corresponding time-dependent energy density

$T_{00}=\frac{1}{8}m^{2}A(\mathbf{x})^{2}=\rho_{\rm{TF}}(\mathbf{x}),$ (15)

and the oscillating spatial components

$\displaystyle\begin{gathered}T_{ij}=\frac{1}{24}\\!m^{2}\\!A(\mathbf{x})^{2}\\!\cos(2mt\\!+\\!2\alpha(\mathbf{x}))\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\\
\times\begin{pmatrix}\\!-(\epsilon_{\rm{S}}^{2}\\!+4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}c_{\chi}\\!-3(\epsilon_{\rm{T}}^{2}\\!+\epsilon_{\rm{V}}^{2}c_{2\xi}))&\\!-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}s_{\chi}\\!+3\epsilon_{\rm{V}}^{2}s_{2\xi}&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}c_{\xi}\\!+3\epsilon_{\rm{T}}c_{\chi-\xi})\\\
\\!-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}s_{\chi}\\!+3\epsilon_{\rm{V}}^{2}s_{2\xi}&-(\epsilon_{\rm{S}}^{2}\\!-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}c_{\chi}\\!-3\epsilon_{\rm{T}}^{2}\\!+3\epsilon_{\rm{V}}^{2}c_{2\xi})&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}s_{\xi}\\!+3\epsilon_{\rm{T}}s_{\chi-\xi})\\\
2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}c_{\xi}+3\epsilon_{\rm{T}}c_{\chi-\xi})&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}s_{\xi}+3\epsilon_{\rm{T}}s_{\chi-\xi})&5\epsilon_{\rm{S}}^{2}-3\epsilon_{\rm{T}}^{2}+3\epsilon_{\rm{V}}^{2}\end{pmatrix}.\\\
\quad\\!\\!\\!\\!\\!\\!\\!\\!\\!\end{gathered}$ (19)

Here we have used the notation $c_{x}=\cos x$ and $s_{x}=\sin x$ following
Ref. (Armaleo  _et al._ , 2020a). Note that although the anisotropic pressure
$T_{ij}$ oscillate with the frequency,

$\omega=2m,$ (20)

it averages to zero on the cosmological time scales, much larger than the
oscillation period. Thus, the coherently oscillating spin-2 massive fields can
be considered the pressureless non-relativistic matter. However, as shown
below, the oscillating spatial components induce oscillation in spacetime
metric, leading to observable effects in PTAs.

The oscillating massive spin-2 fields generate disturbance in the galactic-
scale flat metric, which can be expressed as

$ds^{2}=-\left[1+2\Phi(t,\mathbf{x})\right]dt^{2}+\left[1+2\Psi(t,\mathbf{x})\right]\delta_{ij}dx^{i}dx^{j}+h_{ij}(t,\mathbf{x})dx^{i}dx^{j},$
(21)

where the perturbations are described by $\Phi$, $\Psi$ and $h_{ij}$. The
spatial perturbation has been particularly decomposed into a trace part $\Psi$
and a traceless part $h_{ij}$ as the spatial components of the energy-momentum
are anisotropic. Accordingly, we also decompose the symmetric $T_{ij}$ into

$T_{ij}=\frac{1}{3}\delta_{ij}T^{k}_{\,\,k}+\left(T_{ij}-\frac{1}{3}\delta_{ij}T^{k}_{\,\,k}\right),$
(22)

where the first term corresponds the trace part and the second term
corresponds to the traceless part. From Eq. (19), we obtain

$T^{k}_{\,\,k}=\frac{1}{8}m^{2}A(\mathbf{x})^{2}\cos(2mt+2\alpha(\mathbf{x})),$
(23)

and

$\displaystyle\begin{gathered}\qquad
T_{ij}-\frac{1}{3}\delta_{ij}T^{k}_{\,\,k}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\\
=\frac{1}{24}m^{2}\\!A(\mathbf{x})^{2}\cos(2mt\\!+\\!2\alpha(\mathbf{x}))\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\\
\times\begin{pmatrix}-(2\epsilon_{\rm{S}}^{2}\\!-2\epsilon_{\rm{T}}^{2}\\!+\epsilon_{\rm{V}}^{2}\\!+4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}c_{\chi}\\!-3\epsilon_{\rm{V}}^{2}c_{2\xi})&-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}s_{\chi}\\!+3\epsilon_{\rm{V}}^{2}s_{2\xi}&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}c_{\xi}\\!+3\epsilon_{\rm{T}}c_{\chi\\!-\\!\xi})\\\
\\!-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}s_{\chi}\\!+3\epsilon_{\rm{V}}^{2}s_{2\xi}&-(2\epsilon_{\rm{S}}^{2}\\!-2\epsilon_{\rm{T}}^{2}\\!+\epsilon_{\rm{V}}^{2}\\!-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}c_{\chi}\\!+3\epsilon_{\rm{V}}^{2}c_{2\xi})&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}s_{\xi}\\!+3\epsilon_{\rm{T}}s_{\chi\\!-\\!\xi})\\\
2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}c_{\xi}\\!+3\epsilon_{\rm{T}}c_{\chi\\!-\xi})&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}s_{\xi}\\!+3\epsilon_{\rm{T}}s_{\chi\\!-\xi})&2(2\epsilon_{\rm{S}}^{2}-2\epsilon_{\rm{T}}^{2}\\!+\epsilon_{\rm{V}}^{2})\\\
\end{pmatrix}.\\\ \end{gathered}$ (28)

Now we can solve the metric perturbations in Eq. (21). We will start with the
trace part. In the Newtonian gauge, the metric is given by

$ds^{2}=-(1+2\Phi(t,\mathbf{x}))dt^{2}+(1+2\Psi(t,\mathbf{x}))\delta_{ij}dx^{i}dx^{j}.$
(29)

We can split the gravitational potentials $\Phi$ and $\Psi$ into a dominant
time-independent part and a time-dependent part oscillating with the frequency
$\omega$ as follows,

$\displaystyle\Phi(t,\mathbf{x})$ $\displaystyle\approx$
$\displaystyle\Phi_{0}(\mathbf{x})+\Phi_{\rm{osc}}\cos(\omega
t+2\alpha(\mathbf{x})),$ $\displaystyle\Psi(t,\mathbf{x})$
$\displaystyle\approx$
$\displaystyle\Psi_{0}(\mathbf{x})+\Psi_{\rm{osc}}\cos(\omega
t+2\alpha(\mathbf{x})),$ (30)

where $\Phi_{\text{osc}}$ and $\Psi_{\text{osc}}$ are the oscillation
amplitudes. We can determine the time-independent part $\Psi_{0}(\mathbf{x})$
by solving the $00$ component of the linearized Einstein equations,

$\partial_{i}\partial^{i}\Psi=-4\pi GT_{00}=-4\pi
G\rho_{\rm{TF}}(\mathbf{x}).$ (31)

The oscillating part of $\Psi$ can be found by taking the trace of the spatial
components of the linearized Einstein equations,

$-3\ddot{\Psi}+\partial_{i}\partial^{i}(\Phi+\Psi)=4\pi GT^{k}_{\,\,k},$ (32)

of which the time-independent part implies that $\Phi_{0}=-\Psi_{0}$.
Therefore, the oscillation amplitude $\Psi_{\rm{osc}}$ is

$\Psi_{\rm{osc}}=\frac{1}{24}\pi GA^{2}(\mathbf{x})=\frac{\pi
G\rho_{\rm{TF}}(\mathbf{x})}{3m^{2}}.$ (33)

Next, we turn to the traceless part of the metric perturbation, $h_{ij}$. The
metric, in this case, is given by

$ds^{2}=-dt^{2}+\left[\delta_{ij}+h_{ij}(t,\mathbf{x})\right]dx^{i}dx^{j},$
(34)

with $h^{i}_{\,\,i}=0$. After neglecting the suppressed spatial gradients
$\partial_{k}^{2}h_{ij}$, the linearized Einstein equation yields

$\ddot{h}_{ij}=16\pi
G\left(T_{ij}-\frac{1}{3}\delta_{ij}T^{k}_{\,\,k}\right),$ (35)

The traceless part of the metric perturbation can then be obtained as

$\displaystyle\begin{gathered}h_{ij}(t,\mathbf{x})=-\frac{1}{6}\pi
GA^{2}(\mathbf{x})\cos(2mt+2\alpha(\mathbf{x}))\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\\\
\qquad\quad\times\begin{pmatrix}-(2\epsilon_{\rm{S}}^{2}\\!-2\epsilon_{\rm{T}}^{2}\\!+\epsilon_{\rm{V}}^{2}\\!+4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}c_{\chi}\\!-3\epsilon_{\rm{V}}^{2}c_{2\xi})&-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}s_{\chi}\\!+3\epsilon_{\rm{V}}^{2}s_{2\xi}&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}c_{\xi}\\!+3\epsilon_{\rm{T}}c_{\chi\\!-\\!\xi})\\\
\\!-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}s_{\chi}\\!+3\epsilon_{\rm{V}}^{2}s_{2\xi}&-(2\epsilon_{\rm{S}}^{2}\\!-2\epsilon_{\rm{T}}^{2}\\!+\epsilon_{\rm{V}}^{2}\\!-4\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}c_{\chi}\\!+3\epsilon_{\rm{V}}^{2}c_{2\xi})&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}s_{\xi}\\!+3\epsilon_{\rm{T}}s_{\chi\\!-\\!\xi})\\\
2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}c_{\xi}\\!+3\epsilon_{\rm{T}}c_{\chi\\!-\xi})&2\epsilon_{\rm{V}}(\sqrt{3}\epsilon_{\rm{S}}s_{\xi}\\!+3\epsilon_{\rm{T}}s_{\chi\\!-\xi})&2(2\epsilon_{\rm{S}}^{2}-2\epsilon_{\rm{T}}^{2}\\!+\epsilon_{\rm{V}}^{2})\end{pmatrix}.\end{gathered}$
(38)

Similar to the amplitude of the oscillation in the trace part
$\Psi_{\rm{osc}}$, we define the amplitude of the oscillation in the traceless
part by

$\displaystyle h_{\rm{osc}}(\mathbf{x})\equiv\frac{1}{6}\pi
GA^{2}(\mathbf{x})=\frac{4\pi G\rho_{\rm{TF}}(\mathbf{x})}{3m^{2}}.$ (39)

##  III Timing residual induced by tensor ULDM with gravitational effect

The time-dependent metric, Eq. (II) and Eq. (38), induced by the oscillating
fields will cause a redshift and time delay for any propagating signal. One
such signal is the radio pulse emitted by millisecond pulsars, which can be
measured with extremely high timing precision and provide a natural tool for
probing oscillating ULDM. In the following, we will derive the timing residual
induced by tensor ULDM considering the gravitational effect.

Consider a pulse propagating from a pulsar to the Earth, the redshift is
defined as the relative frequency change of the pulse,

$z(t)=\frac{\omega_{0}-\omega_{\rm{obs}}(t)}{\omega_{0}},$ (40)

where $\omega_{0}$ is the angular frequency of the pulse when emitting and
$\omega_{\rm{obs}}(t)$ is the observed angular frequency at the Earth at time
$t$. Integrating the redshift gives the arrival time of the pulse

$R(t)=\int_{0}^{t}dt\,z(t^{\prime}).$ (41)

Solving the photon geodesic equation under the perturbed metric provides the
frequency shift. For the trace part and traceless part, the corresponding
redshifts are (Nomura  _et al._ , 2020),

$\displaystyle z_{\Psi}(t)$
$\displaystyle=\Psi(t,\mathbf{x}_{e})-\Psi(t-|\mathbf{x}_{p}|,\mathbf{x}_{p}),$
(42) $\displaystyle z_{h}(t)$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}[h_{ij}(t,\mathbf{x}_{e})-h_{ij}(t-|\mathbf{x}_{p}|,\mathbf{x}_{p})],$
(43)

where $\hat{\bf{p}}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$
denotes a unit vector pointing from the Earth to the pulsar, and
$\mathbf{x}_{e}$ and $\mathbf{x}_{p}$ represent the location of the Earth and
the pulsar, respectively. Since most monitored pulsars are located relatively
close to the Earth (within 2 kpc), we assume that the tensor ULDM has the same
oscillation amplitude $A$ at both the Earth and the pulsar. The total redshift
is given by

$\displaystyle z_{t}(t)$ $\displaystyle=$ $\displaystyle z_{\Psi}(t)+z_{h}(t)$
(44) $\displaystyle=$ $\displaystyle\frac{1}{24}\pi
GA^{2}\left[\epsilon_{\rm{S}}^{2}(-1-6\cos
2\theta)+\epsilon_{\rm{T}}^{2}(3+6\cos 2\theta)-3\epsilon_{\rm{V}}^{2}\cos
2\theta+8\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}\cos(2\phi-\chi)\sin^{2}\theta\right.$
$\displaystyle\left.\qquad\quad\,\,-2\epsilon_{\rm{V}}\left(2\sqrt{3}\epsilon_{\rm{S}}\cos(\phi-\xi)\sin
2\theta+6\epsilon_{\rm{T}}\cos(\phi-\chi+\xi)\sin
2\theta+3\epsilon_{\rm{V}}\cos 2(\phi-\xi)\sin^{2}\theta\right)\right]$
$\displaystyle\times\left[\cos(2mt+2\alpha(\mathbf{x}_{e}))-\cos\left(2m(t-|\mathbf{x}_{p}|)+2\alpha(\mathbf{x}_{p})\right)\right]$
$\displaystyle=$ $\displaystyle-\frac{2\pi
G\rho_{\rm{TF}}}{3m^{2}}\left[\epsilon_{\rm{S}}^{2}(-1-6\cos
2\theta)+\epsilon_{\rm{T}}^{2}(3+6\cos 2\theta)-3\epsilon_{\rm{V}}^{2}\cos
2\theta+8\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}\cos(2\phi-\chi)\sin^{2}\theta\right.$
$\displaystyle\left.\qquad\quad\,\,-2\epsilon_{\rm{V}}\left(2\sqrt{3}\epsilon_{\rm{S}}\cos(\phi-\xi)\sin
2\theta+6\epsilon_{\rm{T}}\cos(\phi-\chi+\xi)\sin
2\theta+3\epsilon_{\rm{V}}\cos 2(\phi-\xi)\sin^{2}\theta\right)\right]$
$\displaystyle\times\sin(\alpha_{e}-\alpha_{p})\sin\left(2mt+\alpha_{e}+\alpha_{p}\right),$

where we defined the phases in Earth term and pulsar term as
$\alpha_{e}=\alpha(\mathbf{x}_{e})$ and
$\alpha_{p}=\alpha\left(\mathbf{x}_{p}\right)-m\left|\mathbf{x}_{p}\right|$,
respectively. Correspondingly, the timing residuals read,

$\displaystyle R_{t}$ $\displaystyle=$ $\displaystyle\frac{\pi
G\rho_{\rm{TF}}}{3m^{3}}\left[\epsilon_{\rm{S}}^{2}(-1-6\cos
2\theta)+\epsilon_{\rm{T}}^{2}(3+6\cos 2\theta)-3\epsilon_{\rm{V}}^{2}\cos
2\theta+8\sqrt{3}\epsilon_{\rm{S}}\epsilon_{\rm{T}}\cos(2\phi-\chi)\sin^{2}\theta\right.$
(45)
$\displaystyle\left.\qquad\,\,-2\epsilon_{\rm{V}}\left(2\sqrt{3}\epsilon_{\rm{S}}\cos(\phi-\xi)\sin
2\theta+6\epsilon_{\rm{T}}\cos(\phi-\chi+\xi)\sin
2\theta+3\epsilon_{\rm{V}}\cos 2(\phi-\xi)\sin^{2}\theta\right)\right]$
$\displaystyle\times\sin(\alpha_{e}-\alpha_{p})\cos\left(2mt+\alpha_{e}+\alpha_{p}\right).$

Note that the total timing residual of the tensor ULDM depends on the
amplitudes and azimuthal directions of all three sectors (scalar, vector, and
tensor) as well as the pulsar’s spatial orientation, resulting in a
complicated anisotropy. Fig. 1 provides several examples of this complexity:
subfigures (a)-(c) show scenarios where only one of the sectors exists, and
subfigure (d) displays a more representative case containing all three
sectors.

Figure 1:  The redshift induced by the tensor ULDM under four different cases:
(a) only the scalar mode exists:
$\epsilon_{\rm{S}}=1,\epsilon_{\rm{T}}=\epsilon_{\rm{V}}=0$; (b) only the
tensor modes exist:
$\epsilon_{\rm{T}}=1,\epsilon_{\rm{S}}=\epsilon_{\rm{V}}=0$; (c) only the
vector modes exist:
$\epsilon_{\rm{V}}=1,\epsilon_{\rm{S}}=\epsilon_{\rm{T}}=0$; (d) all five
polarization modes exist with the specific values:
$\epsilon_{\rm{V}}=\epsilon_{\rm{S}}=\epsilon_{\rm{T}}=\frac{\sqrt{3}}{3}$ and
$\chi=\xi=\frac{\pi}{4}$. The orange surface and blue surface represent the
contribution of the trace part Eq. (42) and the traceless part Eq. (43),
respectively; their summation, the total redshift Eq. (44), is depicted by the
green surface.

As reported in Kolb  _et al._ (2023), similar to the situation for spin-1
fields as discussed in Graham  _et al._ (2016), the longitudinal scalar mode
of the massive spin-2 field is more abundantly generated than the
longitudinal-transverse vector and transverse tensor polarization modes. Here
and after, we restrict ourselves to the scalar mode of the tensor ULDM, namely
case a) in Fig. 1 where
$\epsilon_{\rm{S}}=1,\epsilon_{\rm{T}}=\epsilon_{\rm{V}}=0$. To highlight the
physical implications behind this case, we will explicitly derive the timing
residual induced by the scalar sector below. The corresponding equation of
motion can be written as

$\displaystyle\begin{gathered}M_{ij}=\frac{1}{\sqrt{6}}A(\mathbf{x})\cos(mt+\alpha(\mathbf{x}))\begin{pmatrix}-1&0&0\\\
0&-1&0\\\ 0&0&2\end{pmatrix},\end{gathered}$ (47)

and the nonvanishing energy-momentum components are

$\displaystyle T_{00}$ $\displaystyle=$
$\displaystyle\frac{1}{8}m^{2}A(\mathbf{x})^{2},$ (48) $\displaystyle T_{xx}$
$\displaystyle=$ $\displaystyle
T_{yy}=-\frac{1}{24}m^{2}A(\mathbf{x})^{2}\cos(2mt+2\alpha(\mathbf{x})),$ (49)
$\displaystyle T_{zz}$ $\displaystyle=$
$\displaystyle\frac{5}{24}m^{2}A(\mathbf{x})^{2}\cos(2mt+2\alpha(\mathbf{x})).$
(50)

In this case, the metric perturbations $\Psi(\mathbf{x})$ and
$h_{ij}(\mathbf{x})$ induced by the trace and traceless parts of the energy-
momentum can be expressed as

$\displaystyle\Psi(\mathbf{x})$ $\displaystyle=\Psi_{0}+\frac{1}{24}\pi
GA(\mathbf{x})^{2}\cos(2mt+2\alpha(\mathbf{x})),$ $\displaystyle
h_{ij}(\mathbf{x})$ $\displaystyle=\frac{1}{3}\pi
GA(\mathbf{x})^{2}\cos(2mt+2\alpha(\mathbf{x}))\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&-2\end{pmatrix}.$ (51)

These perturbations lead to redshifts $z_{\Psi}$ (Eq. (42)) and $z_{h}$ (Eq.
(43)) in the pulsar timing signal, which add up to yield the total redshift

$\displaystyle z_{t}$ $\displaystyle=$ $\displaystyle z_{\Psi}+z_{h}$ (52)
$\displaystyle=$ $\displaystyle-\frac{\pi G\rho_{\rm{TF}}}{3m^{2}}(1+6\cos
2\theta)\times\left[\cos(2mt+2\alpha(\mathbf{x}_{e}))-\cos\left(2m(t-|\mathbf{x}_{p}|)+2\alpha(\mathbf{x}_{p})\right)\right].$

Integrating the redshift results in the final timing residual,

$R_{\rm{tensor}}=-\frac{\pi G\rho_{\rm{TF}}}{3m^{3}}(1+6\cos
2\theta)\sin(\alpha_{e}-\alpha_{p})\cos\left(2mt+\alpha_{e}+\alpha_{p}\right).$
(53)

The above expression can also be derived by simply plugging the conditions
$\epsilon_{\rm{S}}=1$ and $\epsilon_{\rm{T}}=\epsilon_{\rm{V}}=0$ into Eq.
(45).

Recall that the coherent oscillations of scalar and vector ULDM also generate
timing residuals (Khmelnitsky and Rubakov, 2014; Nomura  _et al._ , 2020),

$\displaystyle R_{\rm{scalar}}$ $\displaystyle=$ $\displaystyle\frac{\pi
G\rho_{\rm{SF}}}{m^{3}}\sin(\alpha_{e}-\alpha_{p})\cos\left(2mt+\alpha_{e}+\alpha_{p}\right),$
(54) $\displaystyle R_{\rm{vector}}$ $\displaystyle=$ $\displaystyle-\frac{\pi
G\rho_{\rm{VF}}}{m^{3}}(1+2\cos
2\theta)\sin(\alpha_{e}-\alpha_{p})\cos\left(2mt+\alpha_{e}+\alpha_{p}\right),$
(55)

where $\rho_{\rm{SF}}$ and $\rho_{\rm{VF}}$ denote the energy density of the
scalar fields and vector fields, respectively. As illustrated in Fig. 2, the
amplitude and angular dependence of the signals induced by ULDM of different
spins differ significantly. Notably, the timing residuals of scalar ULDM are
isotropic, while those of vector ULDM and tensor ULDM exhibit nontrivial
direction dependence.

Figure 2:  The redshift induced by oscillating ULDM with different spins. The
red, purple and green lines represent the cases of the scalar, vector and
tensor ULDM, respectively. The direction of oscillation is chosen along the
z-axis for the vector and tensor ULDM.

##  IV Estimating the limits with current PTA

Several PTA data sets have been utilized to search for the scalar or the
vector ULDM (Kato and Soda, 2020; Porayko  _et al._ , 2018; Wu  _et al._ ,
2022). To make an approximate prediction of the tensor ULDM using the existing
results, we take the average of Eq. (45) over the celestial sphere $\hat{p}$,
which gives

$\displaystyle\bar{R}_{t}=\sqrt{\frac{23}{15}}\frac{\pi
G\rho_{\rm{TF}}}{m^{3}}\times\sin(\alpha_{e}-\alpha_{p})\cos\left(2mt+\alpha_{e}+\alpha_{p}\right).$
(56)

Because of the system’s symmetry, the average timing residuals from the
gravitational effect do not depend on the five parameters describing the
quadruple but only on the energy density of the tensor ULDM $\rho_{\rm{TF}}$.
To utilize previous results and facilitate later discussion, we replace
$\rho_{\rm{TF}}$ with the unique quantity $h_{\rm{osc}}$ from the
gravitational effect that is directly constrained by PTA experiments. Using
the relationship Eq. (39), we obtain

$\displaystyle\bar{R}_{t}=\frac{3}{4}\sqrt{\frac{23}{15}}\frac{h_{\rm{osc}}}{m}\times\sin(\alpha_{e}-\alpha_{p})\cos\left(2mt+\alpha_{e}+\alpha_{p}\right).$
(57)

In Fig. 3, we present the constraint on the amplitude $h_{\rm{osc}}$ for
tensor ULDM, which is reproduced from the bounds on the vector ULDM by the
second dataset release of the Parkes PTA collaboration (PPTA DR2) (Wu  _et
al._ , 2022).

Moreover, the spin-2 massive fields are universally coupled to the SM matter
as a requirement of self-consistency of the model, which would cause a time
delay for a pulse signal. The timing residual from the coupling effect is
given by (Armaleo  _et al._ , 2020b),

$\displaystyle\bar{R}_{c}=\sqrt{\frac{1}{15}}\frac{2\alpha\sqrt{\rho_{\rm{TF}}}}{m^{2}M_{\rm{P}}}\times\sin(\frac{\alpha_{e}-\alpha_{p}}{2})\cos\left(mt+\frac{\alpha_{e}+\alpha_{p}}{2}\right),$
(58)

which depends on both the coupling $\alpha$ and the energy density
$\rho_{\rm{TF}}$. Note that the timing residuals induced by the coupling
effect, $\bar{R}_{c}$, and the gravitational effect, $\bar{R}_{t}$, have
frequencies of $m$ and $2m$, respectively. If the spin-2 fields with a certain
mass $m$ make up a significant fraction of the dark matter and the two effects
are comparable, we will observe two signals in the PTA frequency band, with
one having a frequency that is twice that of the other. Theoretically, it is
difficult to determine the values of the two quantities $\alpha$ and
$\rho_{\rm{TF}}$ from the coupling effect, as they are degenerate. However, we
can first use the pure gravitational effect, given by Eq. (56), which depends
solely on $\rho_{\rm{TF}}$, to determine whether massive spin-2 fields make up
all dark matter. If we cannot exclude this possibility, we can estimate an
upper bound for $\alpha$ using Eq. (57) and Eq. (58),

$\displaystyle\alpha<\frac{3\sqrt{23}}{2}\frac{M_{\rm{P}}}{\sqrt{\rho_{\rm{0}}}}mh_{\rm{osc}},$
(59)

where $\rho_{0}=0.4\rm{GeV/cm^{3}}$ is the measured local energy density
(Salucci  _et al._ , 2010). If the spin-2 fields cannot make up all the dark
matter, we need to estimate $\alpha$ using the $\rho_{\rm{TF}}$ obtained from
the gravitational effect. In Fig. 3, we demonstrate several cases with
different values of the coupling $\alpha$. We estimate that the coupling could
get constrained at the level of $10^{-6}\sim 10^{-5}$ in the mass range
$m<5\times 10^{-23}\rm{eV}$ with current PTA data. This result is consistent
with the constraint derived in Ref. Armaleo  _et al._ (2020b) using the PPTA
12-year data set, because both results are obtained under the same assumption
that $\rho_{\rm{TF}}=\rho_{0}$. However, we should emphasize that the green
dashed line representing the coupling effect based on Eq. (59) cannot be
extended to cross the left side of black dot-dashed line representing the
predicted gravitational effect derived from $\rho_{0}$. Moreover, if future
PTA experiments can constrain that $\rho_{\rm{TF}}<\rho_{0}$, then Eq. (59) is
no longer valid and we need to reestimate the upper bound on $\alpha$ using an
updated $\rho_{\rm{TF}}$.

Figure 3:  PPTA bounds on the potential amplitude $h_{\rm{osc}}$ induced by
tensor ULDM versus the ULDM mass $m$. The black dot-dashed line denotes the
model amplitude $h_{\rm{osc}}$ when assuming
$\rho_{\rm{TF}}=0.4\rm{GeV/cm^{3}}$. The green dashed line represents the
deductive coupling effect from $h_{\rm{osc}}$ for several values of the ULDM
coupling $\alpha$.

##  V conclusion

In this work, we explore the gravitational effect of tensor ULDM on pulsar
timing residuals. As a form of energy, the oscillating tensor field perturbs
the spacetime geometry, resulting in the timing residuals that take the form
of a sine function with a frequency equal to twice the particle mass, as shown
in Eq. (45). Although the timing residuals of tensor ULDM depend on all five
parameters describing the spin-2 quadrupole, it is pointed by (Kolb  _et al._
, 2023) that the gravitationally-produced massive spin-2 particles are
predominantly longitudinally polarized, leading to a simpler expression Eq.
(53) that still retains anisotropies. We find that the pulsar timing signal
induced by the tensor ULDM exhibits a different angular dependence than that
by the scalar and vector ULDM, as shown in Fig. 2.

Besides the gravitational effect, the coupling between the tensor ULDM and SM
matter also produces sine-shaped timing residuals, depending both on the
coupling $\alpha$ and the energy density of tensor fields $\rho_{\rm{TF}}$. In
previous studies, the upper bound of $\alpha$ was estimated by assuming that
the tensor ULDM constitutes all dark matter, which is reasonable since the
current sensitivity is insufficient to exclude the assumption. However, as
experimental precision continues to improve, it may become possible to use the
gravitational effects to determine whether massive spin-2 fields can account
for all dark matter. In such cases, it would be necessary to re-evaluate the
coupling parameter using an updated energy density of the tensor ULDM. We
present an estimation using this approach with PPTA DR2 in Fig. 3.

While PTA’s current sensitivity has limitations in efficiently constraining
ULDM, many other cosmological and astrophysical probes have already provided
important constraints on ULDM. For instance, cosmic microwave background put
percent-level constraints on ULDM within the mass range
$10^{-32}{\rm{eV}}<m<10^{-24}{\rm{eV}}$ (Hlozek  _et al._ , 2018; Poulin  _et
al._ , 2018), and constraints from the galaxy ration curves (Bar  _et al._ ,
2018) and Milky Way satellite population (Nadler  _et al._ , 2019) impose a
lower limit on the ULDM mass of $m\sim(10^{-22}{\rm{eV}},10^{-21}{\rm{eV}})$.
The Lyman-alpha forest provides the most stringent constraint, disfavoring
$m>2\times 10^{-20}{\rm{eV}}$ for ultralight fields constituting the entire
dark matter (Rogers and Peiris, 2021). Nevertheless, there is currently no
consensus on ULDM mass constraints, with some exceptions favoring lower ULDM
masses, such as a few times $10^{-23}{\rm{eV}}$ based on the stellar
kinematics of dwarf spheroidal galaxies (González-Morales  _et al._ , 2017;
Kendall and Easther, 2020). On the one hand, many of these constraints are
subject to their own astrophysical uncertainties. For instance, the constraint
from Lyman-alpha forest relies on accurate modeling of fluctuations from
temperature and impact of ionizing background, while the dwarf spheroidal
bound may face challenges due to assumptions made about the sphericity of
halos in these galaxies. On the other hand, these constraints also largely
depend on assumption regarding the structure formation in the ULDM scenario,
and most astrophysical simulations are based on the scalar/axion setting. To
provide more robust constraints, detailed interplay between ULDM physics and
particle physics models, especially baryonic feedback, and simulations on
higher-spin candidates capturing the polarization of vector and tensor ULDM
are essential. Furthermore, the possibility of ultralight fields not
constituting all the dark matter or having multiple masses is also well-
motivated. PTA experiments offer complementary tools to other experiments,
independently constraining the energy density of ULDM. The forthcoming PTA
project, square kilometer arrays (SKA), has the potential to detect a $10\%$
contribution of ULDM to the dark matter for $m<10^{-23}{\rm{eV}}$ within ten
years and place interesting constraints on masses above $10^{-22}{\rm{eV}}$ if
significantly high timing-precision is reached (Porayko  _et al._ , 2018). We
believe cosmological and astrophysical probes across scales and epochs will
offer a consistent understanding of the ULDM.

PTAs also hold the potential to distinguish different sources and even clarify
the nature of the ultralight fields, since the timing residuals induced by
scalar ULDM, vector ULDM, and tensor ULDM, as well as those induced by
continuous gravitational waves emitted by a single supermassive black hole
binary Zhu  _et al._ (2014); Falxa  _et al._ (2023), exhibit distinct angular
dependence. However, just as disentangling multiple astrophysical and
cosmological stochastic gravitational wave backgrounds is challenging (Kaiser
_et al._ , 2022), it may also be difficult to differentiate these
monochromatic signals from each other. Therefore, investigating how they can
be distinguished is an interesting research topic. One potential solution is
to study the correlation pattern of different signals (Unal  _et al._ , 2022).
The underlying concept is that when the same common signal modulates the
arrival time of pulses from various pulsars, the resulting timing residuals of
any two pulsars will be correlated. Just as gravitational waves generate a
Hellings-Downs pattern (Hellings and Downs, 1983) of correlation, the ULDM
signal with different spins should exhibit distinct correlation forms. While
some studies have pointed out the existence of the correlations and attempted
to calculate the form in specific cases (Armaleo  _et al._ , 2020b; Liang and
Trodden, 2021), further in-depth research is still needed.

###### Acknowledgements.

We acknowledge the use of HPC Cluster of ITP-CAS. QGH is supported by the
grants from NSFC (Grant No. 12250010, 11975019, 11991052, 12047503), Key
Research Program of Frontier Sciences, CAS, Grant No. ZDBS-LY-7009, CAS
Project for Young Scientists in Basic Research YSBR-006, the Key Research
Program of the Chinese Academy of Sciences (Grant No. XDPB15). ZCC is
supported by the National Natural Science Foundation of China (Grant No.
12247176 and No. 12247112) and the China Postdoctoral Science Foundation
Fellowship No. 2022M710429.

## References

  * Rubin  _et al._ (1980) V. C. Rubin, Jr. Ford, W. K., and N. Thonnard, “Rotational properties of 21 SC galaxies with a large range of luminosities and radii, from NGC 4605 (R=4kpc) to UGC 2885 (R=122kpc).” Astrophys. J. 238, 471–487 (1980).
  * Rubin  _et al._ (1982) V. C. Rubin, Jr. Ford, W. K., N. Thonnard, and D. Burstein, “Rotational properties of 23Sb galaxies.” Astrophys. J. 261, 439–456 (1982).
  * Faber and Jackson (1976) S. M. Faber and R. E. Jackson, “Velocity dispersions and mass-to-light ratios for elliptical galaxies.” Astrophys. J. 204, 668–683 (1976).
  * Massey  _et al._ (2010) Richard Massey, Thomas Kitching, and Johan Richard, “The dark matter of gravitational lensing,” Rept. Prog. Phys. 73, 086901 (2010), arXiv:1001.1739 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Bertone  _et al._ (2005) Gianfranco Bertone, Dan Hooper, and Joseph Silk, “Particle dark matter: Evidence, candidates and constraints,” Phys. Rept. 405, 279–390 (2005), arXiv:hep-ph/0404175 .
  * Feng (2010) Jonathan L. Feng, “Dark Matter Candidates from Particle Physics and Methods of Detection,” Ann. Rev. Astron. Astrophys. 48, 495–545 (2010), arXiv:1003.0904 [astro-ph.CO] .
  * Bertone and Tait (2018) Gianfranco Bertone and Tim Tait, M. P., “A new era in the search for dark matter,” Nature 562, 51–56 (2018), arXiv:1810.01668 [astro-ph.CO] .
  * Marsh (2016) David J. E. Marsh, “Axion Cosmology,” Phys. Rept. 643, 1–79 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Alcock  _et al._ (2000) C. Alcock _et al._ (MACHO), “The MACHO project: Microlensing results from 5.7 years of LMC observations,” Astrophys. J. 542, 281–307 (2000), arXiv:astro-ph/0001272 .
  * Schumann (2019) Marc Schumann, “Direct Detection of WIMP Dark Matter: Concepts and Status,” J. Phys. G 46, 103003 (2019), arXiv:1903.03026 [astro-ph.CO] .
  * Hu  _et al._ (2000) Wayne Hu, Rennan Barkana, and Andrei Gruzinov, “Cold and fuzzy dark matter,” Phys. Rev. Lett. 85, 1158–1161 (2000), arXiv:astro-ph/0003365 .
  * Hui  _et al._ (2017) Lam Hui, Jeremiah P. Ostriker, Scott Tremaine, and Edward Witten, “Ultralight scalars as cosmological dark matter,” Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Fox  _et al._ (2004) Patrick Fox, Aaron Pierce, and Scott D. Thomas, “Probing a QCD string axion with precision cosmological measurements,” (2004), arXiv:hep-th/0409059 .
  * Niemeyer (2019) Jens C. Niemeyer, “Small-scale structure of fuzzy and axion-like dark matter,” Prog. Part. Nucl. Phys. (2019), 10.1016/j.ppnp.2020.103787, arXiv:1912.07064 [astro-ph.CO] .
  * Gentile  _et al._ (2004) Gianfranco Gentile, P. Salucci, U. Klein, D. Vergani, and P. Kalberla, “The Cored distribution of dark matter in spiral galaxies,” Mon. Not. Roy. Astron. Soc. 351, 903 (2004), arXiv:astro-ph/0403154 .
  * de Blok (2010) W. J. G. de Blok, “The Core-Cusp Problem,” Adv. Astron. 2010, 789293 (2010), arXiv:0910.3538 [astro-ph.CO] .
  * Moore  _et al._ (1999) B. Moore, S. Ghigna, F. Governato, G. Lake, Thomas R. Quinn, J. Stadel, and P. Tozzi, “Dark matter substructure within galactic halos,” Astrophys. J. Lett. 524, L19–L22 (1999), arXiv:astro-ph/9907411 .
  * Klypin  _et al._ (1999) Anatoly A. Klypin, Andrey V. Kravtsov, Octavio Valenzuela, and Francisco Prada, “Where are the missing Galactic satellites?” Astrophys. J. 522, 82–92 (1999), arXiv:astro-ph/9901240 .
  * Khmelnitsky and Rubakov (2014) Andrei Khmelnitsky and Valery Rubakov, “Pulsar timing signal from ultralight scalar dark matter,” JCAP 02, 019 (2014), arXiv:1309.5888 [astro-ph.CO] .
  * Arvanitaki  _et al._ (2010) Asimina Arvanitaki, Savas Dimopoulos, Sergei Dubovsky, Nemanja Kaloper, and John March-Russell, “String Axiverse,” Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Hlozek  _et al._ (2015) Renée Hlozek, Daniel Grin, David J. E. Marsh, and Pedro G. Ferreira, “A search for ultralight axions using precision cosmological data,” Phys. Rev. D 91, 103512 (2015), arXiv:1410.2896 [astro-ph.CO] .
  * Payez  _et al._ (2015) Alexandre Payez, Carmelo Evoli, Tobias Fischer, Maurizio Giannotti, Alessandro Mirizzi, and Andreas Ringwald, “Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles,” JCAP 02, 006 (2015), arXiv:1410.3747 [astro-ph.HE] .
  * Ivanov  _et al._ (2019) M. M. Ivanov, Y. Y. Kovalev, M. L. Lister, A. G. Panin, A. B. Pushkarev, T. Savolainen, and S. V. Troitsky, “Constraining the photon coupling of ultra-light dark-matter axion-like particles by polarization variations of parsec-scale jets in active galaxies,” JCAP 02, 059 (2019), arXiv:1811.10997 [astro-ph.CO] .
  * Pierce  _et al._ (2018) Aaron Pierce, Keith Riles, and Yue Zhao, “Searching for Dark Photon Dark Matter with Gravitational Wave Detectors,” Phys. Rev. Lett. 121, 061102 (2018), arXiv:1801.10161 [hep-ph] .
  * Nomura  _et al._ (2020) Kimihiro Nomura, Asuka Ito, and Jiro Soda, “Pulsar timing residual induced by ultralight vector dark matter,” Eur. Phys. J. C 80, 419 (2020), arXiv:1912.10210 [gr-qc] .
  * Aoki and Mukohyama (2016) Katsuki Aoki and Shinji Mukohyama, “Massive gravitons as dark matter and gravitational waves,” Phys. Rev. D 94, 024001 (2016), arXiv:1604.06704 [hep-th] .
  * Babichev  _et al._ (2016a) Eugeny Babichev, Luca Marzola, Martti Raidal, Angnis Schmidt-May, Federico Urban, Hardi Veermäe, and Mikael von Strauss, “Bigravitational origin of dark matter,” Phys. Rev. D 94, 084055 (2016a), arXiv:1604.08564 [hep-ph] .
  * Babichev  _et al._ (2016b) Eugeny Babichev, Luca Marzola, Martti Raidal, Angnis Schmidt-May, Federico Urban, Hardi Veermäe, and Mikael von Strauss, “Heavy spin-2 Dark Matter,” JCAP 09, 016 (2016b), arXiv:1607.03497 [hep-th] .
  * Aoki and Maeda (2018) Katsuki Aoki and Kei-ichi Maeda, “Condensate of Massive Graviton and Dark Matter,” Phys. Rev. D 97, 044002 (2018), arXiv:1707.05003 [hep-th] .
  * Hassan and Rosen (2012) S. F. Hassan and Rachel A. Rosen, “Bimetric Gravity from Ghost-free Massive Gravity,” JHEP 02, 126 (2012), arXiv:1109.3515 [hep-th] .
  * Marzola  _et al._ (2018) Luca Marzola, Martti Raidal, and Federico R. Urban, “Oscillating Spin-2 Dark Matter,” Phys. Rev. D 97, 024010 (2018), arXiv:1708.04253 [hep-ph] .
  * Gialamas and Tamvakis (2023) Ioannis D. Gialamas and Kyriakos Tamvakis, “Bimetric-affine quadratic gravity,” Phys. Rev. D 107, 104012 (2023), arXiv:2303.11353 [gr-qc] .
  * Poulin  _et al._ (2018) Vivian Poulin, Tristan L. Smith, Daniel Grin, Tanvi Karwal, and Marc Kamionkowski, “Cosmological implications of ultralight axionlike fields,” Phys. Rev. D 98, 083525 (2018), arXiv:1806.10608 [astro-ph.CO] .
  * Iršič  _et al._ (2017) Vid Iršič, Matteo Viel, Martin G. Haehnelt, James S. Bolton, and George D. Becker, “First constraints on fuzzy dark matter from Lyman-$\alpha$ forest data and hydrodynamical simulations,” Phys. Rev. Lett. 119, 031302 (2017), arXiv:1703.04683 [astro-ph.CO] .
  * Ünal  _et al._ (2021) Caner Ünal, Fabio Pacucci, and Abraham Loeb, “Properties of ultralight bosons from heavy quasar spins via superradiance,” JCAP 05, 007 (2021), arXiv:2012.12790 [hep-ph] .
  * Armaleo  _et al._ (2020a) Juan Manuel Armaleo, Diana López Nacir, and Federico R. Urban, “Binary pulsars as probes for spin-2 ultralight dark matter,” JCAP 01, 053 (2020a), arXiv:1909.13814 [astro-ph.HE] .
  * Blas  _et al._ (2020) Diego Blas, Diana López Nacir, and Sergey Sibiryakov, “Secular effects of ultralight dark matter on binary pulsars,” Phys. Rev. D 101, 063016 (2020), arXiv:1910.08544 [gr-qc] .
  * Porayko and Postnov (2014) N. K. Porayko and K. A. Postnov, “Constraints on ultralight scalar dark matter from pulsar timing,” Phys. Rev. D 90, 062008 (2014), arXiv:1408.4670 [astro-ph.CO] .
  * Porayko  _et al._ (2018) Nataliya K. Porayko _et al._ , “Parkes Pulsar Timing Array constraints on ultralight scalar-field dark matter,” Phys. Rev. D 98, 102002 (2018), arXiv:1810.03227 [astro-ph.CO] .
  * Kato and Soda (2020) Ryo Kato and Jiro Soda, “Search for ultralight scalar dark matter with NANOGrav pulsar timing arrays,” JCAP 09, 036 (2020), arXiv:1904.09143 [astro-ph.HE] .
  * Xue  _et al._ (2022) Xiao Xue _et al._ (PPTA), “High-precision search for dark photon dark matter with the Parkes Pulsar Timing Array,” Phys. Rev. Res. 4, L012022 (2022), arXiv:2112.07687 [hep-ph] .
  * Wu  _et al._ (2022) Yu-Mei Wu, Zu-Cheng Chen, Qing-Guo Huang, Xingjiang Zhu, N. D. Ramesh Bhat, Yi Feng, George Hobbs, Richard N. Manchester, Christopher J. Russell, and R. M. Shannon (PPTA), “Constraining ultralight vector dark matter with the Parkes Pulsar Timing Array second data release,” Phys. Rev. D 106, L081101 (2022), arXiv:2210.03880 [astro-ph.CO] .
  * Armaleo  _et al._ (2020b) Juan Manuel Armaleo, Diana López Nacir, and Federico R. Urban, “Pulsar timing array constraints on spin-2 ULDM,” JCAP 09, 031 (2020b), arXiv:2005.03731 [astro-ph.CO] .
  * Sazhin (1978) M. V. Sazhin, “Opportunities for detecting ultralong gravitational waves,” Soviet Ast. 22, 36–38 (1978).
  * Detweiler (1979) Steven L. Detweiler, “Pulsar timing measurements and the search for gravitational waves,” Astrophys. J. 234, 1100–1104 (1979).
  * Foster and Backer (1990) R. S. Foster and D. C. Backer, “Constructing a Pulsar Timing Array,” Astrophys. J. 361, 300 (1990).
  * Kaplan  _et al._ (2022) David E. Kaplan, Andrea Mitridate, and Tanner Trickle, “Constraining fundamental constant variations from ultralight dark matter with pulsar timing arrays,” Phys. Rev. D 106, 035032 (2022), arXiv:2205.06817 [hep-ph] .
  * Sun  _et al._ (2022) Sichun Sun, Xing-Yu Yang, and Yun-Long 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] .
  * Unal  _et al._ (2022) Caner Unal, Federico R. Urban, and Ely D. Kovetz, “Probing ultralight scalar, vector and tensor dark matter with pulsar timing arrays,” (2022), arXiv:2209.02741 [astro-ph.CO] .
  * Kolb  _et al._ (2023) Edward W. Kolb, Siyang Ling, Andrew J. Long, and Rachel A. Rosen, “Cosmological gravitational particle production of massive spin-2 particles,” (2023), arXiv:2302.04390 [astro-ph.CO] .
  * Graham  _et al._ (2016) Peter W. Graham, Jeremy Mardon, and Surjeet Rajendran, “Vector Dark Matter from Inflationary Fluctuations,” Phys. Rev. D 93, 103520 (2016), arXiv:1504.02102 [hep-ph] .
  * Salucci  _et al._ (2010) P. Salucci, F. Nesti, G. Gentile, and C. F. Martins, “The dark matter density at the Sun’s location,” Astron. Astrophys. 523, A83 (2010), arXiv:1003.3101 [astro-ph.GA] .
  * Hlozek  _et al._ (2018) Renée Hlozek, David J. E. Marsh, and Daniel Grin, “Using the Full Power of the Cosmic Microwave Background to Probe Axion Dark Matter,” Mon. Not. Roy. Astron. Soc. 476, 3063–3085 (2018), arXiv:1708.05681 [astro-ph.CO] .
  * Bar  _et al._ (2018) Nitsan Bar, Diego Blas, Kfir Blum, and Sergey Sibiryakov, “Galactic rotation curves versus ultralight dark matter: Implications of the soliton-host halo relation,” Phys. Rev. D 98, 083027 (2018), arXiv:1805.00122 [astro-ph.CO] .
  * Nadler  _et al._ (2019) Ethan O. Nadler, Vera Gluscevic, Kimberly K. Boddy, and Risa H. Wechsler, “Constraints on Dark Matter Microphysics from the Milky Way Satellite Population,” Astrophys. J. Lett. 878, 32 (2019), [Erratum: Astrophys.J.Lett. 897, L46 (2020), Erratum: Astrophys.J. 897, L46 (2020)], arXiv:1904.10000 [astro-ph.CO] .
  * Rogers and Peiris (2021) Keir K. Rogers and Hiranya V. Peiris, “Strong Bound on Canonical Ultralight Axion Dark Matter from the Lyman-Alpha Forest,” Phys. Rev. Lett. 126, 071302 (2021), arXiv:2007.12705 [astro-ph.CO] .
  * González-Morales  _et al._ (2017) Alma X. González-Morales, David J. E. Marsh, Jorge Peñarrubia, and Luis A. Ureña López, “Unbiased constraints on ultralight axion mass from dwarf spheroidal galaxies,” Mon. Not. Roy. Astron. Soc. 472, 1346–1360 (2017), arXiv:1609.05856 [astro-ph.CO] .
  * Kendall and Easther (2020) Emily Kendall and Richard Easther, “The Core-Cusp Problem Revisited: ULDM vs. CDM,” Publ. Astron. Soc. Austral. 37, e009 (2020), arXiv:1908.02508 [astro-ph.CO] .
  * Zhu  _et al._ (2014) X. J. Zhu _et al._ , “An all-sky search for continuous gravitational waves in the Parkes Pulsar Timing Array data set,” Mon. Not. Roy. Astron. Soc. 444, 3709–3720 (2014), arXiv:1408.5129 [astro-ph.GA] .
  * Falxa  _et al._ (2023) M. Falxa _et al._ (IPTA), “Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array,” Mon. Not. Roy. Astron. Soc. 521, 5077–5086 (2023), arXiv:2303.10767 [gr-qc] .
  * Kaiser  _et al._ (2022) Andrew R. Kaiser, Nihan S. Pol, Maura A. McLaughlin, Siyuan Chen, Jeffrey S. Hazboun, Luke Zoltan Kelley, Joseph Simon, Stephen R. Taylor, Sarah J. Vigeland, and Caitlin A. Witt, “Disentangling Multiple Stochastic Gravitational Wave Background Sources in PTA Data Sets,” Astrophys. J. 938, 115 (2022), arXiv:2208.02307 [astro-ph.CO] .
  * Hellings and Downs (1983) R. w. Hellings and G. s. Downs, “UPPER LIMITS ON THE ISOTROPIC GRAVITATIONAL RADIATION BACKGROUND FROM PULSAR TIMING ANALYSIS,” Astrophys. J. 265, L39–L42 (1983).
  * Liang and Trodden (2021) Qiuyue Liang and Mark Trodden, “Detecting the stochastic gravitational wave background from massive gravity with pulsar timing arrays,” Phys. Rev. D 104, 084052 (2021), arXiv:2108.05344 [astro-ph.CO] .