# Detecting the Stochastic Gravitational Wave Background from Massive Gravity
with Pulsar Timing Arrays

Qiuyue Liang<EMAIL_ADDRESS>Center for Particle Cosmology, Department
of Physics and Astronomy, University of Pennsylvania, Philadelphia,
Pennsylvania 19104, USA  Mark Trodden<EMAIL_ADDRESS>Center for
Particle Cosmology, Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia, Pennsylvania 19104, USA

###### Abstract

We explore the potential of Pulsar Timing Arrays (PTAs) such as NANOGrav,
EPTA, and PPTA to detect the Stochastic Gravitational Wave Background (SGWB)
in theories of massive gravity. In General Relativity, the function describing
the dependence of the correlation between the arrival times of signals from
two pulsars on the angle between them is known as the Hellings-Downs curve. We
compute the analogous overlap reduction function for massive gravity,
including the additional polarization states and the correction due to the
mass of the graviton, and compare the result with the Hellings-Downs curve.
The primary result is a complete analytical form for the analog Hellings-Downs
curve, providing a starting point for future numerical studies aimed at a
detailed comparison between PTA data and the predictions of massive gravity.
We study both the massless limit and the stationary limit as checks on our
calculation, and discuss how our formalism also allows us to study the impact
of massive spin-2 dark matter candidates on data from PTAs.

##  I Introduction

One of the most exciting recent developments in astrophysics and cosmology has
been the direct detection of gravitational waves, and the rapid use of this
technique to extend our cosmological observations beyond those made using
electromagnetic radiation or neutrinos. Much attention has rightly been paid
to the information we can glean from observing the gravitational waves
generated by individual compact sources, such as black hole-black hole and
black hole-neutron star mergers. These observations are providing crucial
insights into how these astrophysical bodies behave, and strict constraints on
allowed deviations of the theory of gravity from General Relativity (GR).

Beyond observing individual systems, the totality of gravitational waves
produced within our Hubble volume contributes to an overall stochastic
background of gravitational waves, with the potential to provide complementary
information about astrophysical and cosmological objects, and about the
underlying theory that governs the generation and propagation of gravitational
waves. This stochastic background is the scientific target of a number of
current and upcoming projects, such as the North American NanoHertz
Observatory for Gravitational Waves (NANOGrav), the European Pulsar Timing
Array (EPTA), and the Parkes Pulsar Timing Array (PPTA). These collaborations
use correlations among the precision timings of the arrivals of signals from
tens of millisecond pulsars as a way to discover minute perturbations in
spacetime as gravitational waves permeate the universe. The signal-to-noise
ratio of these observations increases over time, and thus they provide a
particularly interesting probe.

With sensitivities in the tens to hundreds of microHertz, a primary target of
these pulsar timing arrays (PTAs) is supermassive black hole binary systems.
However, they can also provide novel tests of proposed new physics. Some
proposed new sources arise in the matter sector of our theories, including the
background of gravitational waves produced by decaying cosmic string loops,
and that generated during cosmological phase transitions. On the other hand,
new physics can also arise in the gravitational sector of our theories,
affecting how gravitational waves are generated, how they propagate through
spacetime, and opening up the possibility of the production of entirely new
polarization states. Indeed, work has already been carried out exploring the
impact of different polarization modes in PTA gravitational wave searches
Dalang  _et al._ (2021); Qin  _et al._ (2019).

In this paper, we consider the potential of PTAs to constrain or potentially
discover evidence for stochastic gravitational wave background signals in
theories of massive gravity. In particular, we will be interested in comparing
the predictions of ghost-free massive gravity, with 5 polarization modes, to
that of GR.

In pure GR, PTA observations can be used to probe the stochastic gravitational
wave background (SGWB) in the following way. When a gravitational wave passes,
it perturbs the spacetime around the pulsar, and thus changes the frequency of
the signal in a way that can be described as a (non-cosmological) redshift.
Using pulsar observations, one can then measure the correlation function of
the integral of this redshift. Assuming an isotropic SGWB, this correlator can
be decomposed into two parts: the energy density, and the overlap reduction
function. The energy density part encodes the power spectrum, describing the
amplitude as a function of frequency. The overlap reduction function describes
the spatial dependence of the correlator as a function of the angle between
two arbitrary pulsars. In GR, this shape function was calculated in a
particular approximation by Hellings and Downs Hellings and Downs (1983).

In massive gravity both the energy density and the overlap reduction function
receive corrections, partly due to the non-zero mass, but also due to the
presence of additional polarization modes. Since the amplitude and frequency-
dependence of the energy density largely depend on the source of the
gravitational wave, and we are interested in the effects on the SGWB, in this
paper we focus on the correction to the overlap reduction function.

In this paper we first derive an analytical result, which can be directly
compared with the Hellings-Downs curve in GR, for all polarization modes in
massive gravity. In the case of tensor modes, our general result agrees with
that in Lee  _et al._ (2010) up to a normalization factor. We also analyze the
corrections from vector and scalar modes in this theory. These results involve
a particular approximation, the validity of which we numerically verify for
each mode. We then compute this analog Hellings-Downs curve in two opposite
regimes — the massless limit, and the stationary limit. One of the most model-
independent and rigorous constraints on the mass of graviton is that it must
be less than $10^{-23}$ eV de Rham  _et al._ (2017); Shao (2021) 111Note,
however, that in other situations much tighter constraints may hold Gupta and
Desai (2018)., whereas the sensitivity of current PTA measurements is in the
$1-100$nHz frequency range ($\sim 10^{-24}-10^{-22}$ eV). There are two
qualitatively different ways for a graviton with sufficient energy to be
detected in this range. The first is if the graviton mass is sufficiently low
(or even zero, as in GR) that the graviton is relativistic. The massless limit
of our analysis is relevant to this case. The other possibility is that the
graviton mass is in the range $10^{-24}-10^{-23}$eV, in which case the
stationary limit of our analysis is then appropriate. Beyond massive gravity,
it is possible that other spin 2 particles, such as ultra-light dark matter,
which might comprise the galactic halos surrounding pulsars, may also
contribute to any measured PTA signal (see Armaleo  _et al._ (2020) for a
related idea). This is a second situation in which the stationary limit of our
calculations is relevant.

Compared with previous work Lee  _et al._ ; Gair  _et al._ (2015); Qin  _et
al._ (2021); Lee (2013), here we focus particularly on the specific case of
ghost-free massive gravity, we investigate more general settings in which such
massive spin-2 excitations might be relevant, such as in dark matter
applications, and we carry out a detailed analysis of the validity of the
Hellings-Downs approximation.

The structure of this paper is as follows. In Sec.II, we review the
polarization tensors in the theory of ghost-free linearized massive gravity.
In Sec. III, we discuss the corrections arising from massive gravity to the
signal measured in PTA experiments. We calculate the change in the frequency
of the pulses, and obtain the overlap reduction function for all five
polarization modes. In Sec. IV, we then numerically verify that an important
approximation to the overlap reduction function is valid within the frequency
range and for the distances of pulsars relevant to current PTA experiments. We
analytically calculate expressions for this analog Hellings-Downs curve,
compare the result with that in GR, and discuss how the graviton mass and
extra polarization states affect the observed signals. Throughout the paper,
we will follow the $+---$ signature convention. We will write $p^{2}=0$ for
the electromagnetic signals received on earth from the pulsars, and
$k^{2}=m^{2}$ for the massive gravitational wave momentum.

##  II Polarization tensors for massive gravitational waves

The starting point of an analysis of massive gravitational waves is to
construct the polarization tensors describing the different polarization
states of the graviton. To do this, we employ a technique that allows us to
write down a form for such tensors using the simpler and more familiar
polarization vectors for spin-1 fields.

It is well known that a massive spin-1 field contains three modes — two
transverse modes and one longitudinal mode. Its polarization vectors can be
expressed as Gleisberg  _et al._ (2003); Han  _et al._ (1999),

$\begin{array}[]{l}\epsilon_{\mu}^{\pm}(k)=\frac{1}{\sqrt{2}}(0,\cos\theta\cos\varphi\mp
i\sin\varphi,\cos\theta\sin\varphi\pm i\cos\varphi,-\sin\theta)\\\
\epsilon_{\mu}^{0}(k)=\frac{1}{\sqrt{k^{2}}}\left(|\bm{k}|,k_{0}\hat{\bm{\Omega}}\right)=\left(\frac{|\bm{k}|}{m},\frac{k_{0}}{m}\hat{\bm{\Omega}}\right)\
,\quad\hat{\bm{\Omega}}=(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)\
,\end{array}$ (1)

where

$k_{\mu}=\left(k_{0},|\bm{k}|\sin\theta\cos\varphi,|\bm{k}|\sin\theta\sin\varphi,|\bm{k}|\cos\theta\right)=k_{0}\left(1,\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)\
,$ (2)

with $\hat{\bm{\Omega}}$ the unit vector of the gravitational wave’s spatial
direction.

Now, for massive spin 2 fields, the analogous action is the Fierz-Pauli action

$S=\int
d^{4}x\left[\frac{1}{2}\partial_{\lambda}h_{\mu\nu}\partial^{\lambda}h^{\mu\nu}-\partial_{\mu}h_{\nu\lambda}\partial^{\nu}h^{\mu\lambda}+\partial_{\mu}h^{\mu\nu}\partial_{\nu}h-\frac{1}{2}\partial_{\lambda}h\partial^{\lambda}h+\frac{1}{2}m^{2}\left(h_{\mu\nu}h^{\mu\nu}-h^{2}\right)\right]\
,$ (3)

describing a total of five propagating polarization modes. Using the
polarization vectors for a massive spin-1 field, given above, we can construct
the polarization tensors for a massive spin 2 field in the following way
Gleisberg  _et al._ (2003); Han  _et al._ (1999)

$\displaystyle\epsilon_{\mu\nu}^{(i)}$ $\displaystyle=$
$\displaystyle\left\\{\epsilon_{\mu\nu}^{(+2)},\epsilon_{\mu\nu}^{(+1)},\epsilon_{\mu\nu}^{(0)},\epsilon_{\mu\nu}^{(-1)},\epsilon_{\mu\nu}^{(-2)}\right\\}$
(4) $\displaystyle=$
$\displaystyle\left\\{\epsilon_{\mu}^{+}\epsilon_{\nu}^{+},\frac{1}{\sqrt{2}}\left(\epsilon_{\mu}^{+}\epsilon_{\nu}^{0}+\epsilon_{\mu}^{0}\epsilon_{\nu}^{+}\right),\frac{1}{\sqrt{6}}\left(\epsilon_{\mu}^{+}\epsilon_{\nu}^{-}+\epsilon_{\mu}^{-}\epsilon_{\nu}^{+}-2\epsilon_{\mu}^{0}\epsilon_{\nu}^{0}\right),\frac{1}{\sqrt{2}}\left(\epsilon_{\mu}^{-}\epsilon_{\nu}^{0}+\epsilon_{\mu}^{0}\epsilon_{\nu}^{-}\right),\epsilon_{\mu}^{-}\epsilon_{\nu}^{-}\right\\}$

These polarization tensors are transverse $k^{\mu}\epsilon_{\mu\nu}=0$,
traceless $\epsilon^{\mu}_{~{}\mu}=0$ and orthonormal
$\epsilon^{s}_{\mu\nu}\epsilon^{s^{\prime}\mu\nu}=\delta^{ss^{\prime}}$. Since
we are working at the linearized level, we can write the spin 2 piece as a
plane wave (see the discussion in Isi and Stein (2018) and references
therein),

$h_{\mu\nu}(x)=\frac{1}{2\pi}\int
d^{4}k\frac{2\delta(|\bm{k}|^{2}-(k_{0}^{2}-m^{2}))}{|\bm{k}|}e^{ikx}h_{\mu\nu}(k)=\int_{-\infty}^{\infty}df\int_{\text{sky}}d^{2}\hat{\bm{\Omega}}~{}e^{i2\pi
f\left(t-\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\bm{x}\right)}h_{\mu\nu}\left(f,\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)\
,$ (5)

where $2\pi f=k_{0}$ is the frequency of the spin 2 particle. Note that the
$\delta$-function is introduced to impose the on-shell condition $\square
h_{\mu\nu}+m^{2}h_{\mu\nu}=0\Leftrightarrow k_{0}^{2}=|\bm{k}|^{2}+m^{2}$,
therefore in the second equality we have integrated over the magnitude of
$|\bm{k}|$. Note also that, compared with the conventions often used in
quantum field theory, we have defined the transformed quantity with an extra
factor of $|\bm{k}|$ in the denominator, and have also distributed the factors
of $2\pi$ in the definitions of the Fourier transform and its inverse so as to
make the final equality as simple as possible.

Finally we can express $h_{\mu\nu}(k)$ in terms of the polarization tensors,

$h_{\mu\nu}(k)=\sum_{i}h^{(i)}(k)~{}\epsilon_{\mu\nu}^{(i)}\ ,$ (6)

with $i\in\\{0,\pm 1,\pm 2\\}$. We denote the helicity $\pm 2$ polarization
modes as tensor modes, which are related to the usual definitions of cross
($\times$) and plus ($+$) modes in GR by a rotation. Helicity $\pm 1$ modes
are vector modes and the helicity $0$ mode is a scalar mode. We reserve a
comparison with more standard polarization modes in the literature for
Appendix A with a specific choice $\theta=\varphi=0$.

##  III The signal

Having discussed the metric perturbation and its polarization modes, we now
want to know how these quantities affect the observed pulsar signal. Suppose a
pulsar emits a signal of frequency $\nu_{0}$ in flat spacetime. Then, if there
is a gravitational wave passing between the pulsar and our telescopes (e.g.
through the solar system), the measured pulsar frequency $\nu(t)$ will differ
from $\nu_{0}$, leading to an anomalous residual,

$R(t)\equiv\int_{0}^{t}dt^{\prime}\,\left(\frac{\nu_{0}-\nu(t^{\prime})}{\nu_{0}}\right)$
(7)

in the pulse arrival time. The mean square residual $\braket{R^{2}(t)}$,
defined as

$\braket{R^{2}(t)}=\frac{1}{T}\int_{0}^{T}R^{2}(t)dt\ ,$ (8)

where $T$ is the time interval over which the observations are made, is the
crucial quantity that is measured by the PTA system Detweiler (1979).

To calculate this observable, we first note that a pulse follows the null
geodesics of the perturbed spacetime, each of which can be expressed as a null
path described by an affine parameter $\lambda$

$\sigma^{\mu}(\lambda)=s^{\mu}(\lambda)-\frac{1}{2}\eta^{\mu\nu}h_{\nu\gamma}(x(\lambda))s^{\gamma}(\lambda)\
,$ (9)

where $s^{\mu}\equiv\nu(1,-\hat{\bm{p}})\equiv\nu\hat{p}$ is the corresponding
light path in Minkowski spacetime that points from the pulsar to earth, with
frequency $\nu$. We will write $h_{\mu\nu}$ in place of
$h_{\mu\nu}(x(\lambda))$ in the following discussion. The expression (9)
satisfies the geodesic equation222Note that the definition of the Christoffel
symbol used here has a minus sign compared to that adopted in Ref.Anholm  _et
al._ (2009).,

$\frac{d\sigma^{\mu}}{d\lambda}=-\Gamma_{\alpha\beta}^{\mu}\sigma^{\alpha}\sigma^{\beta}=-\frac{1}{2}\eta^{\mu\nu}\left(h_{\nu\alpha,\beta}+h_{\nu\beta,\alpha}-h_{\alpha\beta,\nu}\right)s^{\alpha}s^{\beta}+\mathcal{O}(h_{\mu\nu}^{2})\
,$ (10)

where it is useful to note that, in $(+---)$ signature,
$h^{0}_{~{}\alpha}=h_{0\alpha}$ and
$h_{\alpha\beta}^{~{}~{}~{},0}=h_{\alpha\beta,0}\equiv\dot{h}_{\alpha\beta}$.
Now, since $s^{\mu}$ is a flat space null vector, we have
$s^{\mu}=dx^{\mu}/d\lambda$, and thus the time component of Eq.(10) is

$\displaystyle\frac{d\sigma^{0}}{d\lambda}$ $\displaystyle=$ $\displaystyle-
h_{0\alpha,\beta}s^{\alpha}s^{\beta}+\frac{1}{2}\dot{h}_{\alpha\beta}s^{\alpha}s^{\beta}=-\frac{dh_{0\alpha}}{d\lambda}s^{\alpha}+\frac{1}{2}\dot{h}_{\alpha\beta}s^{\alpha}s^{\beta}\
,$ (11)

where we have used $\frac{dh}{d\lambda}=\frac{\partial h}{\partial
x^{\beta}}\frac{dx^{\beta}}{d\lambda}=\frac{\partial h}{\partial
x^{\beta}}s^{\beta}$.

On the other hand, using the definition of $\sigma^{\alpha}$ in terms of
$s^{\alpha}$ we can eliminate $d\sigma^{0}/d\lambda$ to yield

$\frac{ds^{0}}{d\lambda}=-\frac{1}{2}\frac{dh_{0\alpha}}{d\lambda}s^{\alpha}+\frac{1}{2}\dot{h}_{\alpha\beta}s^{\alpha}s^{\beta}+\mathcal{O}(h^{2})\
,$ (12)

where we have used that the geodesic equation (10) implies that
$ds^{\mu}/d\lambda=\mathcal{O}(h)$ to neglect a term proportional to
$\frac{1}{2}h_{0\alpha}\frac{ds^{\alpha}}{d\lambda}$.

As can be seen from (5), a plane wav metric perturbation propagating in a
specific direction $\hat{\bm{\Omega}}$ is given by

$h_{\mu\nu}\left(t-\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\bm{x}\right)=\int_{-\infty}^{\infty}df~{}e^{i2\pi
f\left(t-\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\bm{x}\right)}h_{\mu\nu}\left(f,\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)\
,$ (13)

We can express $dh_{\mu\nu}/d\lambda$ as

$\frac{dh_{\mu\nu}(t-\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\bm{x})}{d\lambda}=\frac{\partial
h_{\mu\nu}}{\partial x^{0}}\frac{dx^{0}}{d\lambda}+\frac{\partial
h_{\mu\nu}}{\partial\hat{\bm{\Omega}}\cdot\bm{x}}\frac{d\hat{\bm{\Omega}}\cdot\bm{x}}{d\lambda}=\frac{\partial
h_{\mu\nu}}{\partial
x^{0}}\left(\frac{dx^{0}}{d\lambda}-\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\frac{d\bm{x}}{d\lambda}\right)=\dot{h}_{\mu\nu}\nu\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)\
,$ (14)

where $d\bm{x}/d\lambda=-\nu\hat{\bm{p}}$ is the spatial momentum of the
pulsar signal.

Substituting this expression into Eq.(12), and using that $s^{0}=\nu$ we then
obtain,

$\displaystyle\frac{ds^{0}}{d\lambda}=\frac{d\nu}{d\lambda}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{dh_{0\alpha}}{d\lambda}s^{\alpha}+\frac{1}{2\nu\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}\frac{d{h}_{\alpha\beta}}{d\lambda}s^{\alpha}s^{\beta}$
(15) $\displaystyle=$
$\displaystyle\nu\left(-\frac{1}{2}\frac{dh_{0\alpha}}{d\lambda}\hat{p}^{\alpha}+\frac{1}{2\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}\frac{dh_{\alpha\beta}}{d\lambda}\hat{p}^{\alpha}\hat{p}^{\beta}\right)$

Thus we have,

$\frac{d\nu}{\nu}=-\frac{1}{2}dh_{0\alpha}\hat{p}^{\alpha}+\frac{1}{2\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}dh_{\alpha\beta}\hat{p}^{\alpha}\hat{p}^{\beta}\
,$ (16)

and integrating this yields,

$\log\left(\frac{\nu(t)}{\nu_{0}}\right)=-\frac{\hat{p}^{\alpha}}{2}\Delta
h_{0\alpha}+\frac{\hat{p}^{\alpha}\hat{p}^{\beta}}{2\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}\Delta
h_{\alpha\beta}\ ,$ (17)

where

$\Delta h_{\alpha\beta}\equiv
h_{\alpha\beta}\left(t_{\mathrm{p}},\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)-h_{\alpha\beta}\left(t_{\mathrm{e}},\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)$
(18)

is the difference between the metric perturbation at the pulsar, and that
received at the earth after traveling along the direction $\hat{\bm{\Omega}}$.

Exponentiating both sides and expanding at $\mathcal{O}(h)$, we can then
define the redshift as the fractional change in frequency via,

$\displaystyle-z\equiv\frac{\nu(t)-\nu_{0}}{\nu_{0}}$ $\displaystyle=$
$\displaystyle-\frac{\hat{p}^{\alpha}}{2}\Delta
h_{0\alpha}+\frac{\hat{p}^{\alpha}\hat{p}^{\beta}}{2\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}\Delta
h_{\alpha\beta}\ .$ (19)

It is convenient to choose a coordinate system in which the distance between a
pulsar and the earth is denoted by $L$, and we write Anholm  _et al._ (2009)
(note, the photon is massless and travels at the speed of light),

$\bm{x}_{e}=0,\quad\bm{x}_{p}=L\hat{\bm{p}},\quad t_{p}=t_{e}-L\equiv t-L\ .$
(20)

It is also convenient to work in Fourier space, and so we focus on the
quantity $h^{(i)}\left(f,\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)$,
defined in Eq.(13), and use this to express the difference in metric
perturbations Eq.(18) as

$\displaystyle\begin{aligned} \Delta
h_{\mu\nu}=\int_{-\infty}^{\infty}df&e^{i2\pi ft}\left(e^{-i2\pi
fL\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}-1\right)\times\sum_{i}h^{(i)}\left(f,\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)\epsilon_{\mu\nu}^{(i)}(\hat{\bm{\Omega}})\
.\end{aligned}$ (21)

The redshift can then be written as a function of frequency as,

$z(t,\hat{\bm{\Omega}})=\int_{-\infty}^{\infty}dfe^{i2\pi ft}\left(e^{-i2\pi
fL\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}-1\right)\sum_{i}h^{(i)}\left(f,\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)F^{(i)}(\hat{\bm{\Omega}})$
(22)

where we have defined the so-called receiving function
$F^{(i)}(\hat{\bm{\Omega}})$ as,

$F^{(i)}(\hat{\bm{\Omega}})\equiv-\frac{\hat{p}^{\mu}\hat{p}^{\nu}}{2\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}\epsilon^{(i)}_{\mu\nu}+\frac{\hat{p}^{\mu}}{2}\epsilon^{(i)}_{0\mu}\
,$ (23)

which describes how the variation in the frequency depends on the metric
perturbation. These expressions allow us to isolate the Fourier transform of
the redshift as

$z(f,\hat{\bm{\Omega}})=\left(e^{-i2\pi
fL\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}\right)}-1\right)\sum_{i}h^{(i)}\left(f,\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\right)F^{(i)}(\hat{\bm{\Omega}})$
(24)

Since the effect of the stochastic gravitational wave background on the
measured redshift of a given pulsar consists of contributions from
gravitational waves arriving from all directions, a relevant quantity to
calculate is the total redshift

$\tilde{z}(f)\equiv\int_{S^{2}}d^{2}\hat{\bm{\Omega}}~{}z(f,\hat{\bm{\Omega}})\
.$ (25)

Finally, the observable that is relevant to PTA data is the two-point
correlation function of this quantity,
$\braket{\tilde{z}(f)\tilde{z}(f^{\prime})}$. We will assume (see, for example
the NANOGrav 12.5-year results Arzoumanian  _et al._ (2020) or the most recent
results from PPTA Goncharov  _et al._ (2021)) that the frequency-dependence of
this power spectrum is independent of any spatial correlations. Using this so-
called common-process model, we can separate the observable into two distinct
important pieces. The first of these is $\braket{h^{(i)2}}$, which is related
to the fractional energy density of gravitational waves at a given frequency,
$\Omega_{\mathrm{gw}}(|f|)\equiv(3M_{\mathrm{P}}^{2}H_{0}^{2})^{-1}d\rho_{\mathrm{gw}}/d\ln
f$ via Anholm  _et al._ (2009),

$\left\langle
h^{(i)*}(f,\hat{\bm{\Omega}})h^{(i^{\prime})}\left(f^{\prime},\hat{\bm{\Omega}}^{\prime}\right)\right\rangle=\frac{3H_{0}^{2}}{32\pi^{3}}\delta^{2}\left(\hat{\bm{\Omega}},\hat{\bm{\Omega}}^{\prime}\right)\delta_{ii^{\prime}}\delta\left(f-f^{\prime}\right)\times|f|^{-3}\Omega_{\mathrm{gw}}(|f|)\
.$ (26)

The second piece describes the shape of the signal, and is referred to as the
overlap reduction function,

$\Gamma(|f|)=\beta\sum_{i}\int_{S^{2}}d^{2}\hat{\bm{\Omega}}\left(e^{i2\pi
fL_{1}\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}_{1}\right)}-1\right)\times\left(e^{-i2\pi
fL_{2}\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}_{2}\right)}-1\right)F_{1}^{(i)}(\hat{\bm{\Omega}})F_{2}^{(i)}(\hat{\bm{\Omega}})\
,$ (27)

where $\beta$ is a normalization factor introduced to impose $\Gamma(|f|)=1$
for coincident, co-aligned detectors. Together, these yield

$\braket{\tilde{z}(f)\tilde{z}(f^{\prime})}=\frac{3H_{0}^{2}}{32\beta\pi^{3}}\delta^{2}\left(\hat{\bm{\Omega}},\hat{\bm{\Omega}}^{\prime}\right)\delta_{ii^{\prime}}\delta\left(f-f^{\prime}\right)\times|f|^{-3}\Omega_{\mathrm{gw}}(|f|)\Gamma(|f|)\
.$ (28)

In the remainder of this paper, we will almost entirely focus on the overlap
reduction function.

The above analysis is quite general. When the underlying theory is General
Relativity, then the only propagating perturbation of the metric corresponds
to the spin-2 graviton. However, in more general theories, more excitations,
of various spins, may also contribute. In the particular example of massive
gravity, scalar, vector, and tensor contributions decouple from one another,
since their kinetic terms can be diagonalized Hinterbichler (2012), and so we
may define an overlap reduction function for each type of perturbation

$\Gamma_{I}(|f|)=\beta_{I}\sum_{i}\int_{S^{2}}d^{2}\hat{\bm{\Omega}}\left(e^{i2\pi
fL_{1}\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}_{1}\right)}-1\right)\left(e^{-i2\pi
fL_{2}\left(1+\frac{|\bm{k}|}{k_{0}}\hat{\bm{\Omega}}\cdot\hat{\bm{p}}_{2}\right)}-1\right)F_{1}^{(i)}(\hat{\bm{\Omega}})F_{2}^{(i)}(\hat{\bm{\Omega}})\
,$ (29)

where $i\in\text{Modes of type }I$ and $I=T,V,S$, represents tensor, vector or
scalar modes respectively. The full two-point function is a sum over all
contributions and takes the schematic form

$\braket{\tilde{z}^{2}}\propto\left(\frac{\Omega_{T}}{\beta_{T}}\Gamma_{T}+\frac{\Omega_{V}}{\beta_{V}}\Gamma_{V}+\frac{\Omega_{S}}{\beta_{S}}\Gamma_{S}\right)=\frac{\Omega_{T}}{\beta_{T}}\Gamma_{T}\left(1+\frac{\Gamma_{V}}{\Gamma_{T}}\frac{\Omega_{V}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{V}}+\frac{\Gamma_{S}}{\Gamma_{T}}\frac{\Omega_{S}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{S}}\right)\
.$ (30)

If we are unable to discriminate among different polarization modes, then any
signal detected by PTA would be interpreted as a tensor contribution. With
this in mind, it is convenient to define an “effective” overlap reduction
function as,

$\tilde{\Gamma}_{T}=\Gamma_{T}+\Gamma_{V}\frac{\Omega_{V}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{V}}+\Gamma_{S}\frac{\Omega_{S}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{S}}\
.$ (31)

Note, however, that the different polarization modes have different frequency
dependencies, and thus $\Omega_{V}/\Omega_{T}$ is in general a function of
frequency. However, this frequency dependence is limited over the sensitivity
range of current PTA observations, and we will therefore approximate
$\Omega_{V}/\Omega_{T}$ as a frequency-independent quantity.

##  IV The overlap reduction function in massive gravity

In this section, we explore the effective overlap reduction function Eq.(31)
in massive gravity. We will separately discuss the behaviors of tensor,
vector, and scalar modes, before combining their effects to obtain the total
signal.

Our primary goal is to obtain an approximate analytic expression. In the case
of GR this can be done by dropping the exponential terms in the analogous
expression to (29), leading to the Hellings-Downs curve. We seek to make the
same approximation here, neglecting the exponential factors in Eq.(29), and
referring to the simplified quantity as the analog Hellings-Downs curve. This
is appealing, since the exponential factors remove any hope of an analytical
expression for the overlap reduction function, and also make numerical
evaluations of this quantity significantly more challenging. Our strategy will
be to demonstrate numerically that it is safe to neglect the exponential terms
for each polarization mode, and then to pursue an analytic expression for the
simplified quantity. Once we have separately shown that the analog Hellings-
Downs curves for each mode are reasonable approximations, we will then drop
the exponential factor in Eq.(31) to yield the effective analog Hellings-Downs
curve

$\tilde{\Gamma}_{0,T}=\Gamma_{0,T}+\Gamma_{0,V}\frac{\Omega_{V}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{V}}+\Gamma_{0,S}\frac{\Omega_{S}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{S}}~{}\
.$ (32)

In the massless theory, the primary impact of the exponential factor is to
introduce a damped oscillatory behavior. In the massive theory, the situation
is complicated by the appearance of $\frac{|\bm{k}|}{k_{0}}$, which allows for
a new scale in the theory. Therefore, as we will see, there are two different
oscillation modes in the overlap reduction function. One is a slow oscillation
with frequency $2\pi\frac{|\bm{k}|}{k_{0}}$, and the other is a fast
oscillation with frequency $2\pi$. When we integrate over the direction vector
$\hat{\bm{\Omega}}$, it is the slow oscillation that exhibits the damping
behavior seen in GR. As we approach the massless limit, where
$\frac{|\bm{k}|}{k_{0}}\to 1$, the two frequencies coincide, and the overlap
reduction function damps rapidly for $fL>10$ (which is the minimum value
relevant for PTA observations), rapidly agreeing with the value of the analog
Hellings-Downs curve. Therefore we will see that we can safely drop the
exponential factors in this case. In the stationary limit, the slow
oscillation mode vanishes, since $\frac{|\bm{k}|}{k_{0}}=0$, while the fast
oscillation mode is independent of the momentum, and thus independent of the
angle between any two pulsars. Thus this does not contribute to the angular
dependence of the overlap reduction function, and the analog Hellings-Downs
curve is again a reasonable approximation.

We now demonstrate the above claims numerically, and then compute the analog
Hellings-Downs curve for each polarization mode. After dropping the
exponential factor, the integration over sphere in Eq.(27) yields an
analytical expression, which is our main result of the paper. We then show how
each mode contributes to the effective analog Hellings-Downs curve in Eq.(32),
in order to compare with the original Hellings-Downs curve in GR. We discuss
both the massless limit and the stationary limit for tensor modes and vector
modes, and mainly focus on the stationary limit for scalar modes. Finally, we
combine all the polarization modes together to give an effective analog
Hellings-Downs curve, and compare this to the observed data.

###  IV.1 The tensor modes

We first evaluate the overlap reduction function for tensor modes, and compare
it with the original Hellings-Downs curve. To do this we will ignore the
exponential terms in (29), and define the resulting simplified quantity as the
analog Hellings-Downs curve for tensor modes

$\Gamma_{0,T}=\beta_{T}\sum_{i=\pm
2}\int_{S^{2}}d\hat{\bm{\Omega}}F_{1}^{(i)}(\hat{\bm{\Omega}})F_{2}^{(i)}(\hat{\bm{\Omega}})\
,$ (33)

where $F^{(i)}(\hat{\bm{\Omega}})$ is the receiving function defined in
Eq.(23). This is the most direct analogue to the Hellings-Downs curve in
massive gravity. In Fig 1 we choose a representative value of $\xi$ — the
angle between the two pulsars — of $\xi=\pi/8$ and compare a numerical
evaluation of $\Gamma_{T}$ to the simplified quantity (33), for two values of
$\frac{|\bm{k}|}{k_{0}}$. As discussed in Anholm  _et al._ (2009), the
smallest frequency relevant for PTAs is $f_{\text{min}}\sim
0.1~{}\text{yr}^{-1}$, and the closest pulsars used in the observations are at
distances around $L_{\text{min}}\sim 100~{}\text{ly}$, so that $fL\gtrsim 10$.
One can see in Fig 1 that when $fL>10$ the deviation of the full solution
(solid curves) from the simplified one (dashed curves) is less than $\sim
5\%$, and so $\Gamma_{0,T}$ is a reasonable approximation to the full overlap
reduction function. We have also numerically verified that this conclusion
also holds for other values of $\xi$. Thus, for massive tensor modes, and for
$|\bm{k}|/k_{0}$ close to $1$, dropping the exponential factors in (29) and
instead using (33) is a valid approximation.

Figure 1:  Overlap reduction with the exponential factors (solid lines) and
without (dashed lines) for tensor modes. The solid lines are the complete
expression (Eq.(29)) as a function of $fL$, and the dashed lines are results
of Eq.(33) for $\xi=\pi/8$. The numbers next to the dashed lines in the legend
are the values of Eq.(33) for different $\frac{|\bm{k}|}{k_{0}}$ and
$\xi=\pi/8$. Blue lines are those for the massless case, in which
$\frac{|\bm{k}|}{k_{0}}=1$, and red lines represent
$\frac{|\bm{k}|}{k_{0}}=0.9$, as an example of tensor modes in massive
theories. One can see that for $fL>10$, the deviation between the solid lines
and the dashed lines is less than $\sim 5\%$.

An advantage of adopting the simplified form (33) of the overlap reduction
function is that we can derive an analytic form, for general values of $\xi$.
To compute this analog Hellings-Downs curve, we adopt the same coordinate
system as that used in Anholm  _et al._ (2009), in which we take the spatial
part of $\hat{p}_{1}$ to be parallel to the $z$ axis, and the spatial part of
$\hat{p}_{2}$ to lie in the $x-z$ plane, so that we have

$\hat{p}_{1}^{\mu}=(1,0,0,1),\quad\hat{p}_{2}^{\mu}=(1,\sin\xi,0,\cos\xi)\ .$
(34)

The relevant polarization tensors for these tensor modes are then
$\epsilon_{\mu\nu}^{(+2)}=\epsilon_{\mu}^{+}\epsilon_{\nu}^{+}$ and
$\epsilon_{\mu\nu}^{(-2)}=\epsilon_{\mu}^{-}\epsilon_{\nu}^{-}$ where
$\epsilon_{\mu}^{\pm}$ are the spin 1 polarization vectors defined in Eq.(1) .
In fact, these are the exact same polarization tensors as are found in the
usual massless case. Substituting these polarization tensors into the
receiving function (23), we see that the second term, proportional to
$\epsilon^{(i)}_{0\mu}\hat{p}^{\mu}$, vanishes for tensor modes, and so the
numerator in the case of massive tensor modes does not change in comparison
with the massless case. The expression (33) then becomes

$\Gamma_{0,T}=\frac{\beta_{T}}{4}\int
d^{2}\hat{\bm{\Omega}}\frac{\sin^{2}\theta}{2\left(1+\frac{|\bm{k}|}{k_{0}}\cos\theta\right)}\frac{\cos\xi^{2}\sin\theta^{2}-2\cos\theta\sin\theta\cos\xi\sin\xi\cos\varphi+\sin\xi^{2}(\cos\theta^{2}\cos\varphi^{2}-\sin\varphi^{2})}{1+\frac{|\bm{k}|}{k_{0}}\left(\cos\theta\cos\xi+\cos\varphi\sin\theta\sin\xi\right)}\
.$ (35)

Importantly, this differs from the massless case by the presence of the terms
proportional to $|\bm{k}|/k_{0}$ in the denominator. Thus, as expected, in the
massless limit, in which $|\bm{k}|/k_{0}\to 1$, $\Gamma_{0,T}$ reverts to the
original Hellings-Downs curve Hellings and Downs (1983).

The details of performing this integration are not particularly illuminating,
and so we relegate them to Appendix B. The relevant result is

$\displaystyle\Gamma_{0,T}$ $\displaystyle=$
$\displaystyle\frac{-\pi}{6A^{5}}\frac{\beta_{T}}{4}\left(4A\left(-3+\left(-6+5A^{2}\right)\cos\xi\right)+12\left(1+\cos\xi+A^{2}(1-3\cos\xi)\right)\log\frac{1+A}{1-A}\right.$
(36)
$\displaystyle~{}~{}~{}\left.+\frac{3\left(1+2A^{2}(1-2\cos\xi)-A^{4}(1-2\cos\xi^{2})\right)\log
L_{1}}{\sqrt{(1-\cos\xi)\left(2-A^{2}(1+\cos\xi)\right)}}\right)$

where we have defined $A=\frac{|\bm{k}|}{k_{0}}$ for notational convenience,
and where

$L_{1}\equiv\frac{\left(1+2A^{2}(1-2\cos\xi)-A^{4}(1-2\cos\xi^{2})-2A(1-A^{2}\cos\xi)\sqrt{(1-\cos\xi)\left(2-A^{2}(1+\cos\xi)\right)}~{}\right)^{2}}{\left(1-A^{2}\right)^{4}}$
(37)

In its full generality, this is a rather complicated expression. However, two
limits are of phenomenological interest. One of these is the ultra-
relativistic (massless) limit where $A=|\bm{k}|/k_{0}=1-\epsilon$, with
$\epsilon\ll 1$. Expanding to first order in $\epsilon$, we write
$\Gamma_{0,T}=\Gamma_{0,T}^{(0)}+\Gamma_{0,T}^{(1)}+\mathcal{O}(\epsilon^{2})$,
with

$\Gamma_{0,T}^{(0)}=\frac{\beta_{T}}{4}\frac{2\pi}{3}\left(3+\cos\xi+6(1-\cos\xi)\log\frac{1-\cos\xi}{2}\right)\
,$ (38)

the usual Hellings-Downs result, as mentioned earlier, and the leading order
correction given by

$\displaystyle\Gamma_{0,T}^{(1)}=\frac{\beta_{T}}{4}\frac{2\pi\epsilon}{3}\left(15+11\cos\xi+12(2-\cos\xi)\log\frac{\epsilon}{2}-3(2-\cos\xi)\log
L_{1}\right)$ (39)

While this correction may prove useful in future work, we will not focus on it
in this paper, requiring only the leading order result to make our main
points.

A second interesting limit is the stationary one, in which $|\bm{k}|/k_{0}\ll
1$. This limit may be relevant for the proposal of massive gravity in some
allowed parameter regimes, but there are also other settings in which a
massive spin-2 field might contribute to the measurements made by PTAs. One
example is the idea of ultra-light spin-2 dark matter, in which the relevant
particles would comprise the dominant component of galactic halos. The dark
matter mass in such models is such that it may be relevant to the typical
frequency range — $\sim 10^{-9}-10^{-7}$ Hz $=10^{-24}-10^{-22}$ eV — that
PTAs are sensitive to.

In this limit, $0\lesssim|\bm{k}|/k_{0}\ll 1$, it is also possible to show
that it is a reasonable approximation to drop the exponential factors in
$\Gamma_{T}$. In Fig 2 we again choose a representative value the angle
between the two pulsars of $\xi=\pi/8$ and compare a numerical evaluation of
$\Gamma_{T}$ to the simplified quantity (33), for two smaller values of
$|\bm{k}|/k_{0}$. One can see the two oscillation modes, and that as
$|\bm{k}|/k_{0}\to 0$ the fast oscillation becomes more prominent, while the
period of the slow oscillation increases. Since PTA observations take place
over an interval short compared to this period, the slow oscillation can be
safely neglected in our approximation. Indeed, when $fL>10$ the deviation of
the full solution (solid curves) from the simplified one (dashed curves) is
noticeable, but less than a $\sim 25\%$ effect. The fast oscillation, on the
other hand, is independent of the angle, and therefore does not complicate the
integral.

We conclude that $\Gamma_{0,T}$ is a reasonable, but by no means perfect,
approximation to the full overlap reduction function, as long as
$|\bm{k}|/k_{0}$ is not too small. We will find this approximation useful for
carrying out a simplified comparison to data and for providing analytic
insights into the behavior of the relevant quantities. This is in large part
because numerically integrating $\Gamma_{T}$ for many values of $\xi$ is
computationally challenging. However, it is clear that a full numerical
evaluation will be needed for precise predictions. We have also numerically
verified that this conclusion also holds for other values of $\xi$.

Figure 2:  Overlap reduction with the exponential factors (solid lines) and
without (dashed lines) for massive tensor modes for various values of
$0\lesssim|\bm{k}|/k_{0}\ll 1$. The solid lines are the complete expression
(Eq.(29)) as a function of $fL$, and the dashed lines are results of Eq.(33)
for $\xi=\pi/8$. Oscillations are again damped, as in the $|\bm{k}|/k_{0}\sim
1$ case, but here the damping effect is significantly less effective as
$|\bm{k}|/k_{0}$ decreases. The fast oscillation can be seen as wiggles upon
the slow oscillation, which will dominate when taking smaller
$|\bm{k}|/k_{0}$. For $fL>10$, the deviation between the solid lines and the
dashed lines is less than $\sim 25\%$.

It is worth noting that when $\bm{k}/k_{0}=0$, there is no damping effect at
all, since the exponential factors act as pure oscillation terms $2-2\cos(2\pi
fL)$. In this paper, we are particularly focused on the angular dependence, so
we focus on the expression (33), which we have referred to as the analog
Hellings-Downs curve. The analytical expression in the stationary limit is,

$\Gamma_{0,T}=\frac{\beta_{T}}{4}\left(\frac{8\pi}{15}\left(-1+3\cos\xi^{2}\right)+\frac{8\pi}{105}\frac{|\bm{k}|^{2}}{k_{0}^{2}}(-2-3\cos\xi+6\cos\xi^{2}+5\cos\xi^{3})\right)+\mathcal{O}\left(\frac{|\bm{k}|^{4}}{k_{0}^{4}}\right)\
.$ (40)

To obtain more insight into how the graviton mass affects the angular
dependence of tensor modes, in Fig.(3) we plot the analog Hellings-Downs curve
Eq.(36) for different values of $|\bm{k}|/k_{0}$.

Figure 3: Analog Hellings-Downs curve for tensor modes in massive gravity. The
normalization factor $\beta_{T}$ has been chosen to fix $\Gamma_{0,T}=0.5$ at
$\xi=0$ in agreement with the massless case. The blue curve is the massless
original Hellings-Downs curve, while the other lines are the analog Hellings-
Downs curves for different values of $\frac{|\bm{k}|}{k_{0}}$. The red line is
the stationary limit where $\frac{\bm{k}}{k_{0}}=0$. One can see that the
analog Hellings-Downs curve for tensor modes has a suppression at $\xi=0$ and
at $\xi=\pi$.

One can see that even a small deviation from the exact massless limit
corresponding to GR leads to a visible correction to the curve $\Gamma_{0,T}$.
Note that, in Fig.(3), we have chosen the normalization factor
$\beta_{T}=\frac{3}{4\pi}$ for all the different selected values of
$|\bm{k}|/k_{0}$. If we had chosen a different normalization factor so as to
fix $\Gamma_{0,T}=0.5$ at $\xi=0$, then we would have recovered the numerical
result of Lee  _et al._ .

###  IV.2 The vector modes

We now turn to vector modes. As in the tensor case, to simplify our numerical
calculations, we first seek to understand whether we can neglect the
exponential factors in computing the overlap reduction function.

We first note that Eq.(4) implies that the polarization modes for vector modes
are $\epsilon_{\mu\nu}^{(\pm
1)}=\frac{1}{\sqrt{2}}(\epsilon_{\mu}^{\pm}\epsilon_{\nu}^{0}+\epsilon_{\nu}^{\pm}\epsilon_{\mu}^{0})$.
For simplicity, we denote
$\epsilon_{\mu}^{0}=\frac{k_{0}}{m}\left(\frac{|\bm{k}|}{k_{0}},\hat{\bm{\Omega}}\right)$
in the following context. We then use these polarization tensors and the
receiving function Eq.(23) to obtain the analog Hellings-Downs curve for
vector modes as,

$\displaystyle\Gamma_{0,V}$ $\displaystyle=$
$\displaystyle\frac{\beta_{V}}{4}\frac{k_{0}^{2}}{2m^{2}}\int
d^{2}\hat{\bm{\Omega}}\frac{1}{\left(1+A\cos\theta\right)\left(1+A\left(\cos\theta\cos\xi+\cos\varphi\sin\theta\sin\xi\right)\right)}$
(41)
$\displaystyle\times\left\\{\left(A+(2-A)\cos\theta\right)\left(A+(2-A)\cos\theta\cos\xi\right)\cos\xi\sin^{2}\theta\right.$
$\displaystyle~{}-(2-A^{2})\left(A+(2-A)\cos\theta\right)\cos\theta\sin^{2}\theta\sin^{2}\xi~{}\cos^{2}\varphi$
$\displaystyle\left.~{}-\left(\left(A+(2-A)\cos\theta\right)\left(A\cos\theta+(2-A^{2})\cos
2\theta\cos\xi\right)\sin\theta\sin\xi+1\right)\cos\varphi\right\\}$

Finally, we evaluate this using the same method as we used for the tensor
modes, to yield,

$\displaystyle\Gamma_{0,V}=\frac{\beta_{V}}{8}\frac{k_{0}^{2}}{m^{2}}\frac{8\pi}{A^{5}}\left(\frac{A}{3}\left(A^{6}\cos\xi-2A^{4}(5\cos\xi+3)+2A^{2}(11\cos\xi+6)-6(2\cos\xi+1)\right)\right.$
(42) $\displaystyle+$
$\displaystyle\left.\left(A^{2}-1\right)^{2}\left(\left(A^{2}-2\right)\cos\xi-2\right)\log\left(\frac{1-A}{1+A}\right)+\frac{(A^{2}-1)^{2}(1-A^{2}\cos\xi)\log
L_{1}}{2\sqrt{(1-\cos\xi)(2-A^{2}(1+\cos\xi))}}\right)$

where $L_{1}$ is defined in Eq.(37).

In Fig.(4), we plot both the full expression Eq.(29) and the approximate
expression Eq.(41) for two different representative values of the angle $\xi$
between the pulsars. We do this separately for $|\bm{k}|/k_{0}=0.9$ and
$|\bm{k}|/k_{0}=0.2$. In both cases, the overlap reduction function is well-
approximated by the analog Hellings-Downs curve in the region $fL>10$ relevant
for PTA measurements, and we shall therefore adopt this approximation below.

Figure 4: The overlap reduction function and the analog Hellings-Downs curve
for vector modes. The solid lines are the full expression, including
exponential factors, as a function of $fL$, and the dashed lines are the
results of Eq.(42). The top panel is for a relatively large value of
$|\bm{k}|/k_{0}$ and the bottom panel is for a smaller one. In both panels, we
plot two examples: $\xi=\pi/8$ (red) and $\xi=2\pi/3$ (green).

To understand what happens in the massless (GR) limit, we note that the
leading order contribution as $|\bm{k}|/k_{0}\to 1$ is,

$\Gamma_{0,V}=\frac{\beta_{V}}{4}\frac{k_{0}^{2}}{2m^{2}}\frac{8\pi}{3}\cos\xi$
(43)

As one can see, there seems to be a divergence when taking the massless limit,
since $m^{2}/k_{0}^{2}=1-|\bm{k}|^{2}/k_{0}^{2}\to 0$ when $|\bm{k}|/k_{0}\to
1$. However, it is easy to see that this is an artifact of the decomposition
we have chosen. Recall that the actual observable is the time residue of the
two point correlation function of the redshift, which we have separated into
an energy density piece, and the effective overlap reduction function. Once we
drop the exponential factor, we obtain the effective analog Hellings-Downs
curve. Note that observers do not measure $\Gamma_{0,V}$ independently; rather
it is the combination
$\Gamma_{0,V}\frac{\Omega_{V}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{V}}$ that
contributes. In the massless limit, the energy density of the vector mode
vanishes faster than $k_{0}^{2}/m^{2}$. Thus the combination of these two
pieces results in no contribution to the signal from vector modes, as we would
expect. Another way to see this is to use different definitions of the
polarization tensor and $h_{(\pm 1)}$ in Eq.(6), for example following the
convention used in Isi and Stein (2018), where the divergence does not then
appear in what we are calling the analog Hellings-Downs curve, $\Gamma_{0,V}$.

The other interesting limit is the stationary one, where $|\bm{k}|/k_{0}\to
0$. In this case the leading order result is

$\Gamma_{0,V}=\frac{\beta_{V}}{4}\frac{8\pi}{15}(-1+3\cos^{2}\xi)\ .$ (44)

Note therefore that, in contrast with the massless limit, we expect a non-zero
contribution from vector modes in the stationary limit.

To obtain a sense of how the vector mode contributes to the analog Hellings
Downs curve, we choose
$\frac{\Omega_{V}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{V}}=1$ in Eq.(31), and
plot the comparison between $\Gamma_{0,V}$ and the original Hellings-Downs
curve in GR. We show in Fig.(5) that for large $|\bm{k}|/k_{0}$, the vector
modes enhance the curve at $\xi=0$, and suppress it at $\xi=\pi$, whereas for
small $|\bm{k}|/k_{0}$, the vector modes enhance the analog Hellings-Downs
curve at both ends.

Figure 5: Analog Hellings-Down curve for vector modes in massive gravity. The
blue solid curve is the original Hellings-Downs curve in GR, the other curves
are $\Gamma_{0,V}\frac{\Omega_{V}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{V}}$,
for different values of $|\bm{k}|/k_{0}$, describing the contribution of the
vector modes. The plots are made with the parameter
$\frac{\Omega_{V}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{V}}=1$ to manifestly
demonstrate the angular dependence of the vector modes.

###  IV.3 The scalar mode

The story for scalar modes is very similar to that for vector modes. Again, we
define a simplified overlap reduction function, which we refer to as the
analog Hellings-Downs curve for the scalar polarization, via

$\displaystyle\Gamma_{0,S}$ $\displaystyle=$
$\displaystyle\frac{\beta_{S}}{4}\frac{1}{12}\int
d^{2}\hat{\bm{\Omega}}\frac{(1+4A\cos\theta+3\cos
2\theta)}{\left(1+\frac{|\bm{k}|}{k_{0}}\cos\theta\right)\left(1+\frac{|\bm{k}|}{k_{0}}\left(\cos\theta\cos\xi+\cos\varphi\sin\theta\sin\xi\right)\right)}$
(45) $\displaystyle\times\left(\frac{1}{8}\left(-2+6\cos 2\theta+6\cos
2\xi+6\cos 2\theta\cos 2\xi+16A\cos\theta\cos\xi\right)\right.$
$\displaystyle~{}~{}~{}\left.+2(A+3\cos\theta\cos\xi)\sin\theta\sin\xi\cos\varphi+3\sin^{2}\theta\sin^{2}\xi\cos^{2}\varphi\right)\
.$

In Fig.(6) we plot the full scalar overlap reduction function and the analog
Hellings-Downs curve for selected values of the angle $\xi$ between the
pulsars and the value of $|\bm{k}|/k_{0}$. Similar to both tensor and vector
modes, we observe fast and slow oscillations and conclude that the analog
Hellings-Downs curve is a reasonable approximation to the full result for
values of $fL$ of interest in PTA observations. We therefore use this
approximation from now on.

Figure 6: The overlap reduction function and the analog Hellings-Downs curve
for scalar mode. The solid lines are the full expression, including the
exponential factors, as a function of $fL$, and the dashed lines are results
of Eq.(46). One can see that for $fL>10$, it is safe to neglect the
exponential factors.

Evaluating Eq.(45) in the same way yields

$\displaystyle\Gamma_{0,S}$ $\displaystyle=$
$\displaystyle\frac{\beta_{S}}{4}\frac{\pi}{6A^{5}}\left(4A\left(9+4A^{4}+18\cos\xi-3A^{2}(4+5\cos\xi)\right)-12(1-A^{2})(2A^{2}-3-3\cos\xi)\log\frac{1-A}{1+A}\right.$
(46) $\displaystyle\qquad\qquad\left.-\frac{9(1-A^{2})^{2}\log
L_{1}}{\sqrt{(1-\cos\xi)(2-A^{2}(1+\cos\xi))}}\right)\ ,$

where $L_{1}$ is defined in Eq.(37).

We plot the scalar contribution to effective analog Hellings-Downs curve in
Fig.(7), choosing the parameter
$\frac{\Omega_{S}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{S}}=1$. It is well-known
in linearized massive gravity that great care is needed when taking the
massless limit in order to treat the scalar mode correctly. This issue is
sometimes referred to as the vDVZ discontinuity. Since we do not wish to
deviate from the linearized limit, we will not attempt to analyze the massless
limit for the scalar mode in this paper.

The interesting limit we will explore is the stationary limit, where the
analog Hellings-Downs curve for the scalar mode is

$\displaystyle\Gamma_{0,S}$ $\displaystyle=$
$\displaystyle\frac{\beta_{S}}{4}\frac{4\pi}{15}(-1+3\cos^{2}\xi)\ .$ (47)

These results display a similar angular dependence to that of $\Gamma_{0,V}$.
Thus, both vector modes and scalar mode contribute in the stationary limit in
any theory with massive spin-2 excitations.

Figure 7: Analog Hellings-Downs curve for scalar mode in massive gravity. The
blue solid curve is the original Hellings-Downs curve in GR, the other curves
are $\Gamma_{0,S}\frac{\Omega_{S}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{S}}$,
for different values of $|\bm{k}|/k_{0}$, describing the contribution of the
scalar modes. The plots are made choosing the parameter
$\frac{\Omega_{S}}{\Omega_{T}}\frac{\beta_{T}}{\beta_{S}}=1$ to show the
angular dependence of the scalar mode.

###  IV.4 The combined effective analog Hellings-Downs curve

Having separately discussed how each polarization mode contributes to the
effective analog Hellings-Downs curve, we now consider the combined effect.

Recall that in Eq.(30) we separated the correlation function into two pieces —
the power spectrum, which encodes its frequency dependence, and the overlap
reduction function, describing the spatial angular dependence. Recently, two
collaborations — the NANOGrav 12.5-year result Arzoumanian  _et al._ (2020)
and the most recent PPTA result Goncharov  _et al._ (2021) — have claimed
strong evidence for a characteristic power spectrum, but have also claimed
that spatial correlations are not well-described by the Hellings-Downs curve.
Fitting the power spectrum to a power-law model $f^{-\gamma}$, the PPTA
collaboration finds $\gamma\in(1.5,5.5)$ and NANOGrav collaboration finds
$\gamma\in(3.76,6.78)$. There exists a number of different suggestions for
physics that might yield a contribution to the SGWB with a frequency-
dependence in this range. Examples are: supermassive black hole binary systems
($\gamma\sim 13/3$) Burke-Spolaor  _et al._ (2019); primordial gravitational
waves ($\gamma\sim 5$) Grishchuk (2005); and networks of cosmic strings
($\gamma\sim 16/3$) Olmez  _et al._ (2010). It is natural to wonder,
therefore, whether modifications of gravity, such as massive gravity, might
maintain these successful predictions of the frequency-dependence while
modifying the shape of the spatial correlation function, perhaps leading to an
improved agreement with current and upcoming data.

In the massive gravity theory that we have studied here, the vector and scalar
modes can certainly modify the shape of spatial correlation function.
Furthermore, if we ensure that their energy densities are somewhat subdominant
to the tensor one ($\Omega_{V},\Omega_{S}\lesssim\Omega_{T}$), but not so
negligible that they do not contribute to the effective analog Hellings-Downs
curve, then the frequency-dependence should not be strongly modified. In
Fig.(8), we have demonstrated the largest this effect could be by choosing the
parameters $\Omega_{T}=\Omega_{V}=\Omega_{S}$ in the given frequency band of
current PTA observations. One can clearly see that the shape of spatial
correlation function can deviate from the Hellings-Downs curve. For some
allowed values of the parameters in massive gravity, this deviation can be
large enough to be detected as data is accumulated from the current generation
of PTA observing programs. It would be interesting to understand more
comprehensively the extent to which upcoming measurements can probe more
general modifications of gravity, and massive gravity in particular, using
this technique.

Figure 8:  Comparison between the original Hellings-Downs curve and the
effective analog Hellings-Downs curve as defined in Eq.(32). The orange,
dashed line describes the monopole signal. The parameters are chosen to be:
$\frac{|\bm{k}|}{k_{0}}=0.71,\beta_{T}=\beta_{V}=\beta_{S}=3/4\pi,\Omega_{T}=\Omega_{V}=\Omega_{S}$
for the solid red line, and
$\beta_{T}=\beta_{V}=\beta_{S}=3/4\pi,\Omega_{V}/\Omega_{T}=0.8,\Omega_{S}/\Omega_{T}=0.8$
for the dashed red line.

##  V Conclusions

In this paper, we have studied the potential of PTA observations to constrain
or discover new features of gravity, focusing particularly on the case of
massive gravity. PTAs measure the correlation function of arrival times of
pulses from pairs of pulsars as a function of the angle between them. Taking
the sky-averaged and polarization-averaged product of the response of a pair
of Earth-pulsar baselines to a plane wave propagating in a particular
direction, the relevant quantity in GR is the famous Hellings-Downs curve,
which is an approximation to the overlap reduction function part of the
signal’s redshift correlation function. Starting with the Fierz-Pauli action
that describes linearized massive gravity, we have defined the five
polarization tensors of the relevant metric perturbation around flat
spacetime, and have derived how these affect the propagation of signals from
the pulsars to our detectors. We have traced how these changes affect the
overlap reduction function, and have then defined the analog Hellings-Downs
curve for massive gravity.

After numerically justifying the approximations that we make, the main result
of this paper is a full analytical expression for this analog Hellings-Downs
curve for the tensor, vector, and scalar modes. We have analyzed the massless
limit and the stationary limit of these expressions, and have combined the
effects of all the five polarization states together to yield an effective
Hellings-Downs curve. These results indicate that it may be possible to
distinguish massive gravity from GR as future PTA data is accumulated.

Our hope is that the analytic expressions we have derived will be of use to
observers in making detailed comparisons of current and future PTA data with
the predictions of theories with massive spin-2 excitations. Should evidence
for such massive modes be discovered, a natural question is how such a
stochastic gravitational wave background might be generated. In future work we
will explore how such signals might arise in theories with massive spin-2
particles, such as massive gravity and spin-2 dark matter models. For example,
it would be interesting to study how the interactions of spin-2 dark matter
halos can change the effective metric around pulsars, and to understand how
massive gravity might give rise to significant gravitational waves through
supermassive black hole binary systems, or phase transitions in the early
universe.

###### Acknowledgements.

QL thanks Jiewen Chen for helpful discussions, and we thank Kurt Hinterbichler
for useful comments on a draft of this paper. This work is supported in part
by NASA ATP grant 80NSSC18K0694 and by US Department of Energy (HEP) Award DE-
SC0013528.

##  Appendix A Appendix

In this section, we compare the polarization tensors in Eq.(4) to the more
standard definition of polarization tensors in massive gravity theory as has
been discussed in Isi and Stein (2018).

To gain a better understanding of the full five modes in the massive spin-2
theory, we set $\theta=\varphi=0$, so that the polarization vector of a
massive spin-1 particle becomes,

$\epsilon_{\mu}^{\pm}=\frac{1}{\sqrt{2}}(0,1,\mp
i,0),\quad\epsilon_{\mu}^{0}=\frac{1}{m}(|\bm{k}|,0,0,k_{0})\ .$ (48)

Combining these to obtain the polarization tensor basis yields,

$\epsilon_{\mu\nu}^{(+2)}=\frac{1}{2}\left(\begin{array}[]{cccc}0&0&0&0\\\
0&1&-i&0\\\ 0&-i&-1&0\\\
0&0&0&0\end{array}\right),\quad\epsilon_{\mu\nu}^{(-2)}=\frac{1}{2}\left(\begin{array}[]{cccc}0&0&0&0\\\
0&1&i&0\\\ 0&i&-1&0\\\ 0&0&0&0\end{array}\right)$ (49)
$\epsilon_{\mu\nu}^{(+1)}=\frac{1}{2m}\left(\begin{array}[]{cccc}0&|\bm{k}|&-i|\bm{k}|&0\\\
|\bm{k}|&0&0&k_{0}\\\ -i|\bm{k}|&0&0&-ik_{0}\\\
0&k_{0}&-ik_{0}&0\end{array}\right),\quad\epsilon_{\mu\nu}^{(-1)}=\frac{1}{2m}\left(\begin{array}[]{cccc}0&|\bm{k}|&i|\bm{k}|&0\\\
|\bm{k}|&0&0&k_{0}\\\ i|\bm{k}|&0&0&ik_{0}\\\
0&k_{0}&ik_{0}&0\end{array}\right)$

$\epsilon_{\mu\nu}^{(0)}=\frac{1}{\sqrt{6}}\left(\left(\begin{array}[]{cccc}0&0&0&0\\\
0&1&0&0\\\ 0&0&1&0\\\
0&0&0&0\end{array}\right)-\frac{2}{m^{2}}\left(\begin{array}[]{cccc}|\bm{k}|^{2}&0&0&k_{0}|\bm{k}|\\\
0&0&0&0\\\ 0&0&0&0\\\ |\bm{k}|k_{0}&0&0&k_{0}^{2}\end{array}\right)\right)$

Thus, the metric perturbation around flat spacetime can be expressed in matrix
form as,

$h_{\mu\nu}=\left(\begin{array}[]{cccc}-\frac{2}{\sqrt{6}}\frac{|\bm{k}|^{2}}{m^{2}}h^{(0)}&\frac{|\bm{k}|}{2m}\left(h^{(1)}+h^{(-1)}\right)&-i\frac{|\bm{k}|}{2m}\left(h^{(1)}-h^{(-1)}\right)&-\frac{2}{\sqrt{6}}\frac{|\bm{k}|k_{0}}{m^{2}}h^{(0)}\\\
\frac{|\bm{k}|}{2m}\left(h^{(1)}+h^{(-1)}\right)&\left(h^{(2)}+h^{(-2)}\right)+\frac{1}{\sqrt{6}}h^{(0)}&-i\left(h^{(2)}-h^{(-2)}\right)&\frac{k_{0}}{2m}\left(h^{(1)}+h^{(-1)}\right)\\\
-i\frac{|\bm{k}|}{2m}\left(h^{(1)}-h^{(-1)}\right)&-i\left(h^{(2)}-h^{(-2)}\right)&-\left(h^{(2)}+h^{(-2)}\right)+\frac{1}{\sqrt{6}}h^{(0)}&-i\frac{k_{0}}{2m}\left(h^{(1)}-h^{(-1)}\right)\\\
-\frac{2}{\sqrt{6}}\frac{|\bm{k}|k_{0}}{m^{2}}h^{(0)}&\frac{k_{0}}{2m}\left(h^{(1)}+h^{(-1)}\right)&-i\frac{k_{0}}{2m}\left(h^{(1)}-h^{(-1)}\right)&-\frac{2}{\sqrt{6}}\frac{k_{0}^{2}}{m^{2}}h^{(0)}\end{array}\right)\
.$ (50)

In the massless limit, where only two polarization states remain, one can see
that $h_{+}=h^{(2)}+h^{(-2)}$ and $h_{\times}=-i\left(h^{(2)}-h^{(-2)}\right)$
are the normal definitions of the cross mode and the plus mode, respectively,
in GR.

Following the conventions of Isi and Stein (2018), where
$\beta=\frac{|\bm{k}|}{k_{0}}$ and $\alpha\equiv\frac{m}{k_{0}}$, the metric
perturbation is,

$h_{\mu\nu}=\left(\begin{array}[]{cccc}-\beta^{2}h_{l}&\beta
h_{\text{x}}&\beta h_{\text{y}}&-\beta h_{l}\\\ \beta
h_{\text{x}}&h_{+}+\frac{1}{2}\alpha^{2}h_{l}&h_{\times}&h_{\text{x}}\\\ \beta
h_{\text{y}}&h_{\times}&-h_{+}+\frac{1}{2}\alpha^{2}h_{l}&h_{\text{y}}\\\
-\beta h_{l}&h_{\text{x}}&h_{\text{y}}&-h_{l}\end{array}\right)$ (51)

where $\frac{2}{\sqrt{6}}h^{(0)}/\alpha^{2}=h_{l}$,
$(h^{(1)}+h^{(-1)})/2\alpha=h_{\text{x}}$,
$-i(h^{(1)}-h^{(-1)})/2\alpha=h_{\text{y}}$, $(h^{(2)}+h^{(-2)})=h_{+}$ and
$-i(h^{(1)}-h^{(-1)})=h_{\times}$. It is straightforward to show that, up to
an overall factor, the five polarization modes we have defined in (4) are
equivalent to those defined in Isi and Stein (2018).

##  Appendix B Computation of the Hellings-Downs curve

In this appendix, for the convenience of the reader, we describe in detail the
integration of Eq.(35). For simplicity, we denote $\frac{|\bm{k}|}{k_{0}}=A$
in the following equations. Writing $x\equiv\cos\theta$, we have

$\displaystyle\Gamma_{0,T}$ $\displaystyle=$
$\displaystyle\frac{\beta_{T}}{4}\int_{-1}^{1}dx\frac{1-x^{2}}{2(1+Ax)}\times$
(52)
$\displaystyle\int_{0}^{2\pi}d\varphi\frac{\cos\xi^{2}(1-x^{2})-\sin\xi^{2}-2x\sqrt{1-x^{2}}\cos\xi\sin\xi\cos\varphi+\sin\xi^{2}(1+x^{2})\cos\varphi^{2}}{1+Ax\cos\xi+A\sqrt{1-x^{2}}\sin\xi\cos\varphi}$
$\displaystyle=$
$\displaystyle\frac{\beta_{T}}{4}\int_{-1}^{1}dx\frac{1-x^{2}}{2(1+Ax)}\times\int_{0}^{2\pi}d\varphi\left(C_{1}\frac{1}{a+b\cos\varphi}+C_{2}\frac{\cos\varphi}{a+b\cos\varphi}+C_{3}\frac{\cos\varphi^{2}}{a+b\cos\varphi}\right)$

where $C_{i}$ are independent of $\varphi$, and $\sin\theta=+\sqrt{1-x^{2}}$,
since $\sin\theta$ is always positive for $0<\theta<\pi$. We can simplify the
denominator by defining $a\equiv 1+Ax\cos\xi$,$b\equiv
A\sqrt{1-x^{2}}\sin\xi$, and noting that
$a^{2}-b^{2}=1-A^{2}+A^{2}(\cos\xi+x)^{2}$. We can then use standard results
(see, e.g. Gradshteyn and Ryzhik (2007)) to carry out the integration over
$\varphi$, yielding:

$\displaystyle\Gamma_{0,T}=\frac{\beta_{T}}{4}\int_{-1}^{1}dx\frac{1-x^{2}}{2(1+Ax)}\left(\frac{2\pi\left(\cos\xi^{2}(1-x^{2})-\sin\xi^{2}\right)}{\sqrt{1-A^{2}+A^{2}(\cos\xi+x)^{2}}}-2x\cos\xi\frac{2\pi}{A}\left(1-\frac{1+Ax\cos\xi}{\sqrt{1-A^{2}+A^{2}(\cos\xi+x)^{2}}}\right)\right.$
$\displaystyle+\left.\sin\xi^{2}(1+x^{2})2\pi\left(\frac{1}{\sqrt{1-A^{2}+A^{2}(\cos\xi+x)^{2}}}-\frac{1+Ax\cos\xi}{(A\sqrt{1-x^{2}}\sin\xi)^{2}}+\frac{\sqrt{1-A^{2}+A^{2}(\cos\xi+x)^{2}}}{(A\sqrt{1-x^{2}}\sin\xi)^{2}}\right)\right)$
(53)

The integration over $x$ from $-1$ to $1$ is trivial yet tedious. After some
simplification, one recovers the result expressed as in Eq.(36). The
procedures for computing the analogous results for the vector modes, Eq.(41),
and the scalar mode, Eq.(45), are very similar to the above.

## References

  * Dalang  _et al._ (2021) Charles Dalang, Giulia Cusin, and Macarena Lagos, “Polarization distortions of lensed gravitational waves,” (2021), arXiv:2104.10119 [gr-qc] .
  * Qin  _et al._ (2019) Wenzer Qin, Kimberly K. Boddy, Marc Kamionkowski, and Liang Dai, “Pulsar-timing arrays, astrometry, and gravitational waves,” Phys. Rev. D 99, 063002 (2019), arXiv:1810.02369 [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. Lett. 265, L39–L42 (1983).
  * Lee  _et al._ (2010) Kejia Lee, Fredrick A. Jenet, Richard H. Price, Norbert Wex, and Michael Kramer, “Detecting massive gravitons using pulsar timing arrays,” Astrophys. J. 722, 1589–1597 (2010), arXiv:1008.2561 [astro-ph.HE] .
  * de Rham  _et al._ (2017) Claudia de Rham, J. Tate Deskins, Andrew J. Tolley, and Shuang-Yong Zhou, “Graviton Mass Bounds,” Rev. Mod. Phys. 89, 025004 (2017), arXiv:1606.08462 [astro-ph.CO] .
  * Shao (2021) Lijing Shao, “Pulsar tests of the graviton mass,” Astron. Nachr. 342, 300–304 (2021), arXiv:2009.06817 [gr-qc] .
  * Gupta and Desai (2018) Sajal Gupta and Shantanu Desai, “Limit on graviton mass using stacked galaxy cluster catalogs from SPT-SZ, Planck-SZ and SDSS-redMaPPer,” Annals Phys. 399, 85–92 (2018), arXiv:1810.00198 [astro-ph.CO] .
  * Armaleo  _et al._ (2020) Juan Manuel Armaleo, Diana López Nacir, and Federico R. Urban, “Pulsar timing array constraints on spin-2 ULDM,” JCAP 09, 031 (2020), arXiv:2005.03731 [astro-ph.CO] .
  * (9) K. J. Lee, F. A. Jenet, and Richard H. Price, “Pulsar Timing as a Probe of Non-Einsteinian Polarizations of Gravitational Waves,” American Astronomical Society 10.1086/591080.
  * Gair  _et al._ (2015) Jonathan R. Gair, Joseph D. Romano, and Stephen R. Taylor, “Mapping gravitational-wave backgrounds of arbitrary polarisation using pulsar timing arrays,” Phys. Rev. D 92, 102003 (2015), arXiv:1506.08668 [gr-qc] .
  * Qin  _et al._ (2021) Wenzer Qin, Kimberly K. Boddy, and Marc Kamionkowski, “Subluminal stochastic gravitational waves in pulsar-timing arrays and astrometry,” Phys. Rev. D 103, 024045 (2021), arXiv:2007.11009 [gr-qc] .
  * Lee (2013) K. J. Lee, “Pulsar Timing Arrays and Gravity Tests in the Radiative Regime,” Class. Quant. Grav. 30, 224016 (2013), arXiv:1404.2090 [astro-ph.CO] .
  * Gleisberg  _et al._ (2003) Tanju Gleisberg, Frank Krauss, Konstantin T. Matchev, Andreas Schalicke, Steffen Schumann, and Gerhard Soff, “Helicity formalism for spin-2 particles,” JHEP 09, 001 (2003), arXiv:hep-ph/0306182 .
  * Han  _et al._ (1999) Tao Han, Joseph D. Lykken, and Ren-Jie Zhang, “On Kaluza-Klein states from large extra dimensions,” Phys. Rev. D 59, 105006 (1999), arXiv:hep-ph/9811350 .
  * Isi and Stein (2018) Maximiliano Isi and Leo C. Stein, “Measuring stochastic gravitational-wave energy beyond general relativity,” Phys. Rev. D 98, 104025 (2018), arXiv:1807.02123 [gr-qc] .
  * Detweiler (1979) Steven L. Detweiler, “Pulsar timing measurements and the search for gravitational waves,” Astrophys. J. 234, 1100–1104 (1979).
  * Anholm  _et al._ (2009) Melissa Anholm, Stefan Ballmer, Jolien D. E. Creighton, Larry R. Price, and Xavier Siemens, “Optimal strategies for gravitational wave stochastic background searches in pulsar timing data,” Phys. Rev. D 79, 084030 (2009), arXiv:0809.0701 [gr-qc] .
  * Arzoumanian  _et al._ (2020) Zaven Arzoumanian _et al._ (NANOGrav), “The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background,” Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Goncharov  _et al._ (2021) Boris Goncharov _et al._ , “On the evidence for a common-spectrum process in the search for the nanohertz gravitational wave background with the Parkes Pulsar Timing Array,” (2021), arXiv:2107.12112 [astro-ph.HE] .
  * Hinterbichler (2012) Kurt Hinterbichler, “Theoretical Aspects of Massive Gravity,” Rev. Mod. Phys. 84, 671–710 (2012), arXiv:1105.3735 [hep-th] .
  * Burke-Spolaor  _et al._ (2019) Sarah Burke-Spolaor _et al._ , “The Astrophysics of Nanohertz Gravitational Waves,” Astron. Astrophys. Rev. 27, 5 (2019), arXiv:1811.08826 [astro-ph.HE] .
  * Grishchuk (2005) L. P. Grishchuk, “Relic gravitational waves and cosmology,” Phys. Usp. 48, 1235–1247 (2005), arXiv:gr-qc/0504018 .
  * Olmez  _et al._ (2010) S. Olmez, V. Mandic, and X. Siemens, “Gravitational-Wave Stochastic Background from Kinks and Cusps on Cosmic Strings,” Phys. Rev. D 81, 104028 (2010), arXiv:1004.0890 [astro-ph.CO] .
  * Gradshteyn and Ryzhik (2007) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ , seventh ed. (Elsevier/Academic Press, Amsterdam, 2007) pp. xlviii+1171, translated from the Russian, Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger, With one CD-ROM (Windows, Macintosh and UNIX).