compat=1.0.0

MS-TP-24-13

DESY-24-098

On the overlap reduction function of pulsar timing array
searches for gravitational waves in modified gravity

Nina Cordes1, Andrea Mitridate2, Kai Schmitz1, Tobias Schröder1, Kim Wassner1

$1$ | _University of Münster, Institute for Theoretical Physics, 48149 Münster, Germany_
---|---
$2$ | _Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany_

Pulsar Timing Array (PTA) searches for gravitational waves (GWs) aim to detect
a characteristic correlation pattern in the timing residuals of galactic
millisecond pulsars. This pattern is described by the PTA overlap reduction
function (ORF) $\Gamma_{ab}\left(\xi_{ab}\right)$, which is known as the
Hellings–Downs (HD) curve in general relativity (GR). In theories of modified
gravity, the HD curve often receives corrections. Assuming, e.g., a subluminal
GW phase velocity, one finds a drastically enhanced ORF in the limit of small
angular separations between pulsar $a$ and pulsar $b$ in the sky,
$\xi_{ab}\rightarrow 0$. In particular, working in harmonic space and
performing an approximate resummation of all multipole contributions, the auto
correlation coefficient $\Gamma_{aa}$ seems to diverge. In this paper, we
confirm that this divergence is unphysical and provide an exact and analytical
expression for $\Gamma_{aa}$ in dependence of the pulsar distance $L_{a}$ and
the GW phase velocity $v_{\rm ph}$. In the GR limit and assuming a large
pulsar distance, our expression reduces to $\Gamma_{aa}=1$. In the case of
subluminal phase velocity, we show that the regularization of the naive
divergent result is a finite-distance effect, meaning that $\Gamma_{aa}$
scales linearly with $fL_{a}$, where $f$ is the GW frequency. For superluminal
phase velocity (subluminal group velocity), which is relevant in the case of
massive gravity, we correct an earlier analytical result for $\Gamma_{ab}$.
Our results pave the way for fitting modified-gravity theories with
nonstandard phase velocity to PTA data, which requires a proper understanding
of the auto correlation coefficient $\Gamma_{aa}$.

###### Contents

  1. 1 Introduction
  2. 2 Hellings–Downs curve
  3. 3 Legendre decomposition
  4. 4 Auto correlation coefficient
  5. 5 Massive gravity
  6. 6 Conclusions

##  1 Introduction

In recent years, the field of pulsar timing array (PTA) searches for
gravitational waves (GWs) [1] has entered an exciting phase, with several PTAs
now being on the brink of discovering a stochastic GW background at nHz
frequencies [2, 3, 4, 5]. The central observables in PTA measurements are
timing residuals for a set of galactic millisecond pulsars, $R_{a}$, in which
GWs can leave an imprint in the form of a characteristic correlation pattern.
The GW contribution to the timing residuals can be written as a time integral
over GW-induced redshifts [6, 7],

$R_{a}^{\rm
GW}\left(t\right)=\int_{0}^{t}dt^{\prime}\>z_{a}\left(t^{\prime}\right)\,,$
(1)

which quantify the instantaneous shift in the pulse times of arrival caused by
GWs and whose correlator reads

$\left<z_{a}\left(t\right)z_{b}\left(t\right)\right>=\frac{2}{3}\int_{0}^{\infty}df\>\Gamma_{ab}\left(\xi_{ab},f\right)S_{h}\left(f\right)\,.$
(2)

Here, $S_{h}$ is the power spectrum of GW modes with polarization $A$,
frequency $f$, and propagation direction $\hat{\bm{n}}$,

$\left<\tilde{h}_{A}\left(f,\hat{\bm{n}}\right)\tilde{h}_{A^{\prime}}^{*}\left(f^{\prime},\hat{\bm{n}}^{\prime}\right)\right>=\frac{1}{8\pi}\,\delta_{AA^{\prime}}\,\delta\left(f-f^{\prime}\right)\,\delta^{(2)}\left(\hat{\bm{n}}-\hat{\bm{n}}^{\prime}\right)\,S_{h}\left(f\right)\,,$
(3)

and the function $\Gamma_{ab}$ denotes the PTA overlap reduction function
(ORF) (see Ref. [8] for a pedagogical review),

$\Gamma_{ab}\left(\xi_{ab},f\right)=\frac{3}{2}\int\frac{d^{2}\hat{\bm{n}}}{4\pi}\left[1-e^{-2\pi
ifL_{a}\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}/v_{\rm
ph}\right)}\right]\left[1-e^{2\pi
ifL_{b}\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{b}/v_{\rm
ph}\right)}\right]\sum_{A=+,\times}F_{a}^{A}\left(\hat{\bm{n}}\right)F_{b}^{A}\left(\hat{\bm{n}}\right)\,.$
(4)

In this expression, $L_{a}$ and $L_{b}$ denote the distances to pulsars $a$
and $b$, respectively; the unit vectors $\hat{\bm{p}}_{a}$ and
$\hat{\bm{p}}_{b}$ indicate the positions of pulsars $a$ and $b$ on the
celestial sphere, respectively; $\xi_{ab}$ is the angular separation of the
two pulsars in the sky,
$\xi_{ab}=\arccos\left(\hat{\bm{p}}_{a}\cdot\hat{\bm{p}}_{b}\right)$; and
$F_{a}^{+}$ and $F_{a}^{\times}$ are the PTA detector pattern functions [9],

$F_{a}^{A}\left(\hat{\bm{n}}\right)=\frac{p_{a}^{i}p_{a}^{j}e_{ij}^{A}\left(\hat{\bm{n}}\right)}{2\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}/v_{\rm
ph}\right)}\,,$ (5)

where $e_{ij}^{+}$ and $e_{ij}^{\times}$ are the two GW polarization tensors,

$e_{ij}^{+}\left(\hat{\bm{n}}\right)=\hat{\phi}_{i}\hat{\phi}_{j}-\hat{\theta}_{i}\hat{\theta}_{j}\,,\qquad
e_{ij}^{\times}\left(\hat{\bm{n}}\right)=\hat{\phi}_{i}\hat{\theta}_{j}+\hat{\theta}_{i}\hat{\phi}_{j}\,,$
(6)

which can be constructed in terms of the unit basis vectors $\hat{\bm{\phi}}$
and $\hat{\bm{\theta}}$,

$\hat{\bm{\phi}}=\left(\sin\phi,-\cos\phi,0\right)^{\rm
T}\,,\qquad\hat{\bm{\theta}}=\left(\cos\theta\cos\phi,\cos\theta\sin\phi,-\sin\theta\right)^{\rm
T}\,.$ (7)

In this paper, we are interested in the properties of $\Gamma_{ab}$ in
theories beyond general relativity (GR) [10, 11, 12, 13, 14, 15, 16, 17, 18].
In particular, we shall consider the possibility of a nonstandard GW phase
velocity, $v_{\rm ph}\neq 1$, which is realized in various modified-gravity
models. At nHz frequencies, GWs may, e.g., exhibit a dispersion relation of
the form

$\omega\left(k\right)=v_{\rm ph}\,k\,,\qquad k=\left|\bm{k}\right|\,,\qquad
v_{\rm ph}<1\,.$ (8)

In this case, the GW phase and group velocities, $v_{\rm
ph}=\omega\left(k\right)/k$ and $v_{\rm
gr}=\partial\omega\left(k\right)/\partial k$, are identical, and GWs propagate
at subluminal speed. Another attractive benchmark consists in massive gravity,
where

$\omega\left(k\right)=\sqrt{m_{g}^{2}+k^{2}}\,,$ (9)

with $m_{g}$ denoting the graviton mass. Now, the phase velocity becomes
superluminal, while the group velocity remains subluminal. In fact, in the
case of massive gravity, the two velocities are just the inverse of each
other,

$v_{\rm ph}=\sqrt{1+\left(m_{g}/k\right)^{2}}\,,\qquad v_{\rm
gr}=\frac{1}{\sqrt{1+\left(m_{g}/k\right)^{2}}}\,.$ (10)

In order to be able to fit these two benchmark scenarios of modified gravity
to PTA data [19, 20, 21, 22, 23, 24], it is crucial to have a precise
understanding of the ORF, both for separate pulsars at angular separation
$\xi_{ab}$ as well as for individual pulsars. Earlier work mostly focused on
the cross correlations coefficients $\Gamma_{a\neq b}$. The goal of the
present paper therefore is to provide an in-depth discussion of the auto
correlation coefficient $\Gamma_{aa}$, including new exact and analytical
expressions that are valid both for models with sub- and superluminal phase
velocities.

The rest of the paper is organized as follows. In the next section, we will
briefly review the computation of the ORF in GR, which we will use as an
opportunity to fix the overall normalization of the ORF. In Sec. 3, we will
then turn to the decomposition of the ORF in Legrendre polynomials. In the
context of this discussion, we will notably encounter a spurious divergence at
$\xi_{ab}=0$ in the case of subluminal $v_{\rm ph}$, which serves as a main
motivation for much of the subsequent analysis. Indeed, in Sec. 4, we will
present an exact and analytical expression for the auto correlation
coefficient $\Gamma_{aa}$, which confirms that the ORF does in fact not
diverge. In addition, we will derive a useful and considerably simpler
expression for $\Gamma_{aa}$ that manages to approximate the full result with
excellent precision at $v_{\rm ph}<1$. Similarly, in Sec. 5, we will derive a
useful and simple expression for $\Gamma_{aa}$ that approximates the full
result in the case of massive gravity (i.e., for $v_{\rm ph}>1$) with
excellent precision. Sec. 6, finally, contains our conclusions and an outlook
on remaining open questions and our next steps.

##  2 Hellings–Downs curve

The two square brackets in Eq. (4) account for the so-called earth-term and
pulsar-term contributions to the GW-induced redshifts $z_{a}$ and $z_{b}$,
respectively. In the GR limit, $v_{\rm ph}\rightarrow 1$, the pulsar-term
contributions (i.e., the exponential factors inside the square brackets) yield
only fast-oscillating terms that turn out to be numerically negligible when
integrated over. In the GR case, one therefore typically only keeps the earth-
term contributions (i.e., the factors of $1$ inside the square brackets), such
that $\Gamma_{ab}$ obtains the simple form

$\Gamma_{ab}^{\rm HD}\left(\xi_{ab}\right)\quad\overset{a\neq
b}{=}\quad\frac{3}{2}\int\frac{d^{2}\hat{\bm{n}}}{4\pi}\sum_{A=+,\times}\frac{p_{a}^{i}p_{a}^{j}e_{ij}^{A}\left(\hat{\bm{n}}\right)}{2\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}\right)}\frac{p_{b}^{k}p_{b}^{l}e_{kl}^{A}\left(\hat{\bm{n}}\right)}{2\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{b}\right)}\,.$
(11)

Note that this expression no longer depends on the GW frequency $f$. The
integral in Eq. (11) can be evaluated explicitly and results in the
Hellings–Downs (HD) curve [25], the hallmark signature of GWs in PTA data in
GR,

$\Gamma_{ab}^{\rm HD}\left(\xi_{ab}\right)\quad\overset{a\neq
b}{=}\quad\frac{3}{2}\,x_{ab}\ln x_{ab}-\frac{x_{ab}}{4}+\frac{1}{2}\,,\qquad
x_{ab}=\frac{1}{2}\left(1-\cos\theta_{ab}\right)\,.$ (12)

The expressions in Eqs. (11) and (12) apply to the case of separate pulsars,
$a\neq b$, and hence describe the cross correlation between the timing
residuals $R_{a}$ and $R_{b}$. Starting from Eq. (4), we can, however, also
consider the case of just one pulsar, $a=b$, and compute the auto correlation
coefficient $\Gamma_{aa}$. In this case, the correlation between the two
pulsar terms is exactly of the same size as the correlation between the two
earth terms,

$\left[1-e^{-2\pi
ifL_{a}\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}\right)}\right]\left[1-e^{2\pi
ifL_{a}\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}\right)}\right]=2+\textrm{fast-
oscillating terms}\,,$ (13)

which means that we need to multiply the HD curve by an extra factor of $2$ in
the case of identical pulsars,

$\Gamma_{ab}^{\rm
HD}\left(\xi_{ab}\right)=\left(1+\delta_{ab}\right)\left[\frac{3}{2}\,x_{ab}\ln
x_{ab}-\frac{x_{ab}}{4}+\frac{1}{2}\right]\,.$ (14)

Furthermore, for identical pulsars, the two position vectors
$\hat{\bm{p}}_{a}$ and $\hat{\bm{p}}_{b}$ are of course the same, such that
$\xi_{ab}=0$ and $x_{ab}=0$. In this case, the square brackets in Eq. (14)
yield a factor of $\nicefrac{{1}}{{2}}$, which, together with the prefactor
$1+\delta_{ab}$, results by construction in an auto correlation coefficient of
$\Gamma_{aa}=1$. This convenient normalization of the ORF is a consequence of
the factor of $\nicefrac{{2}}{{3}}$ in Eq. (2) and the corresponding factor of
$\nicefrac{{3}}{{2}}$ in Eq. (4).

Figure 1: Coefficients $a_{\ell}$ in the Legendre decomposition of the ORF for
different values of the GW phase velocity $v_{\rm ph}$, see Eq. (15) (left
panel); and the corresponding approximate ORF based on the first $20$ terms in
this decomposition, see Eq. (17) (right panel).

##  3 Legendre decomposition

In the more general case of $v_{\rm ph}\neq 1$, it is less straightforward to
evaluate the ORF in Eq. (4). In particular, if $v_{\rm ph}<1$, it is crucial
to retain the pulsar terms in the integrand. Otherwise, the factor
$1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}/v_{\rm ph}$ in the denominator of the
detector pattern function can result in a spurious numerical singularity. A
common approach in the literature therefore consists in decomposing the ORF
into Legendre polynomials $P_{\ell}$ [26], which allows for an approximate
numerical evaluation of $\Gamma_{ab}$ order by order in higher multipole
contributions,

$\Gamma_{ab}\left(\xi_{ab},f\right)=\sum_{\ell=2}^{\infty}a_{\ell}\left(f\right)P_{\ell}\left(\cos\xi_{ab}\right)\,,\qquad
a_{\ell}\left(f\right)=\frac{3}{2}\left(2\ell+1\right)\,\frac{\left(\ell-2\right)!}{\left(\ell+2\right)!}\frac{\left|c_{\ell}\left(f\right)\right|^{2}}{16}\,.$
(15)

In GR, the coefficients $c_{\ell}$ evaluate to
$c_{\ell}=\left(-1\right)^{\ell}4$, while in the more general case, we need to
compute

$c_{\ell}\left(f\right)=\int_{-1}^{+1}dx\left[1-e^{-i\,2\pi fL\left(1+x/v_{\rm
ph}\right)}\right]\frac{\left(1-x^{2}\right)^{2}}{1+x/v_{\rm
ph}}\,\frac{d^{2}}{dx^{2}}P_{\ell}\left(x\right)\,,$ (16)

assuming all pulsars to be located at roughly the same distance, $L_{a}\sim
L_{b}\sim L$. Keeping only a finite number of terms in the decomposition in
Eq. (15) will always result in a continuous and smooth approximation of the
ORF. By construction, this approach is therefore incapable of encoding the
correct auto correlation coefficient $\Gamma_{aa}$, which represents a
discontinuous jump in the limit of zero angular separation,
$\Gamma_{aa}\neq\lim_{\xi_{ab}\rightarrow 0^{+}}\Gamma_{ab}$.

For subluminal phase velocity, the coefficients $a_{\ell}$ decrease less
slowly with $\ell$ as in GR (see the left panel of Fig. 1), which results in a
slow convergence of the Legendre decomposition in Eq. (15). To see this,
consider truncating the expansion in Eq. (15) after $\ell=20$, which results
in the following approximate ORF,

$\Gamma_{ab}^{(20)}\left(\xi_{ab},f\right)=\sum_{\ell=2}^{20}a_{\ell}\left(f\right)P_{\ell}\left(\cos\xi_{ab}\right)\,.$
(17)

We plot $\Gamma_{ab}^{(20)}$ for different values of the phase velocity in the
right panel of Fig. 1, which illustrates that, for $v_{\rm ph}<1$, naively
summing the first $\mathcal{O}\left(10\right)$ contributions to the Legendre
decomposition does not yet yield a satisfactory result. Because of the slow
convergence of the coefficients $a_{\ell}$, the approximate ORF
$\Gamma_{ab}^{(20)}$ still features spurious oscillations, which ought to be
absent in the full nonperturbative result.

In order to improve the approximate ORF $\Gamma_{ab}^{(20)}$, it has been
proposed to resum the remaining terms in Eq. (15) [9], making use of the fact
that, for $v_{\rm ph}<1$, the coefficients $a_{\ell}$ remain nearly constant
at large $\ell$,

$\Gamma_{ab}^{(\infty)}\left(\xi_{ab},f\right)=\sum_{\ell=2}^{20}a_{\ell}\left(f\right)P_{\ell}\left(\cos\xi_{ab}\right)+\Theta\left(-\log_{10}v\right)a_{20}\left(f\right)\sum_{\ell=21}^{\infty}P_{\ell}\left(\cos\xi_{ab}\right)\,,$
(18)

which can also be written as

$\Gamma_{ab}^{(\infty)}\left(\xi_{ab},f\right)=\sum_{\ell=2}^{20}a_{\ell}\left(f\right)P_{\ell}\left(\cos\xi_{ab}\right)+\Theta\left(-\log_{10}v\right)a_{20}\left(f\right)\left[\frac{1}{\sqrt{2-2\cos\xi_{ab}}}-\sum_{\ell=0}^{20}P_{\ell}\left(\cos\xi_{ab}\right)\right]\,,$
(19)

where we used the following sum rule for the Legendre polynomials,

$\sum_{\ell=0}^{\infty}P_{\ell}\left(\cos\xi_{ab}\right)=\frac{1}{\sqrt{2-2\cos\xi_{ab}}}\,.$
(20)

The factor $\Theta$ in Eq. (19) denotes the Heaviside theta function, i.e.,
$\Theta\left(x\right)=1$ for $x>1$ and $\Theta\left(x\right)=0$ for $x\leq 0$.

Figure 2: Comparison between the resummed ORF $\Gamma_{ab}^{(\infty)}$ in Eq.
(19) (solid lines) and the approximate ORF $\Gamma_{ab}^{(20)}$ in Eq. (17)
(dashed lines) for different subluminal GW phase velocities $v_{\rm ph}$. The
orange line shows $\Gamma_{ab}^{(20)}$ in the GR limit, which requires no
improvement.

The resummed ORF $\Gamma_{ab}^{(\infty)}$ is shown in Fig. 2, from which it is
evident that it significantly improves over the approximate
$\Gamma_{ab}^{(20)}$. For practical applications, $\Gamma_{ab}^{(\infty)}$
thus provides a reasonable approximation of the exact cross correlation
coefficients $\Gamma_{ab}$ — at least, at sufficiently large angular
separations. In the limit $\xi_{ab}\rightarrow 0$, though, the resummed ORF is
unfortunately not well-behaved, since the
$\left(2-2\cos\xi_{ab}\right)^{-1/2}$ term in Eq. (19) causes a divergence at
$\xi_{ab}=0$. This divergence has been known in the literature for some time
[23, 9, 24]. Still, thus far, it has remained unclear how to regularize the
divergence in the resummed ORF, which poses a severe obstacle to using
$\Gamma_{ab}^{(\infty)}$ in realistic data analyses. Indeed, fitting modified-
gravity models with a nonstandard GW phase velocity to PTA data requires
knowledge of the auto correlation coefficients $\Gamma_{aa}$, which are
typically (much) larger than the cross correlation coefficients at
nonvanishing angular separation, $\Gamma_{a=b}\gtrsim\Gamma_{a\neq b}$. As a
consequence, most of the sensitivity to the GW power spectrum $S_{h}$ in Eq.
(3) derives from measuring the auto correlations of individual pulsars,
$\left<R_{a}R_{a}\right>$; the information contained in the cross correlations
$\left<R_{a}R_{b}\right>$, on the other hand, rather serves as a means to
confirm the GW nature of the signal. In other words, a successful fit of a
modified-gravity model to PTA data must be able to reproduce the expected
hierarchy between the auto and cross correlation terms in the correlation
matrix $\left<R_{a}R_{b}\right>$. Without a proper understanding of the auto
correlation coefficients $\Gamma_{aa}$, fit analyses of this type are not
feasible.

Figure 3: Auto correlation coefficient $\Gamma_{aa}$ as a function of the GW
phase velocity $v_{\rm ph}$ for three choices of $fL_{a}$. The thick black
lines show the exact analytical result in Eq. (23); the red and gray lines
refer to the approximate results in Eq.(34) and Eq. (28), respectively.

In summary, we conclude that the resummed ORF in Eq. (19) is not only
incapable of resolving the discontinuous jump from $\lim_{\xi_{ab}\rightarrow
0^{+}}\Gamma_{ab}$ to $\Gamma_{aa}$ at $\xi_{ab}=0$; it does not even yield a
finite result for $\lim_{\xi_{ab}\rightarrow 0^{+}}\Gamma_{ab}$. The reason
for this breakdown of $\Gamma_{ab}^{(\infty)}$ in the zero-$\xi_{ab}$ limit
clearly lies in the assumption that the coefficients $a_{\ell}$ at large
$\ell$ can all be approximated by $a_{20}$. If we did not make this assumption
and instead summed all exact multipole contributions up to
$\ell\rightarrow\infty$, we would expect to obtain a finite result at
$\xi_{ab}=0$. In Ref. [24], the authors speculate in particular that, when
pushing the Legendre decomposition to $\ell$ values at least as large as
$\ell\sim fL_{a}$, where $L_{a}$ is the pulsar distance and $f$ the GW
frequency, finite-distance effects should become important and regulate the
coefficients $a_{\ell}$ and hence ultimately lead to a finite ORF at
$\xi_{ab}=0$. In practice, however, such an approach appears infeasible (or at
least, less practicable), since $fL_{a}\sim\mathcal{O}\left(100\cdots
1000\right)$ for typical pulsars.

##  4 Auto correlation coefficient

In this paper, we will therefore follow a different approach. We shall abandon
the decomposition in Legendre polynomials and the analysis of the ORF in
harmonic space and return to the expression for the ORF in Eq. (4). As it
turns out, this expression is perfectly suitable to compute the auto
correlation coefficient $\Gamma_{aa}$ in modified-gravity scenarios with
$v_{\rm ph}\neq 1$. To do so, let us consider a single pulsar $a$ located in
the direction of $\hat{\bm{p}}_{a}$ in the sky. Without loss of generality, we
can always choose our coordinate system such that
$\hat{\bm{p}}_{a}=\left(0,0,1\right)^{\rm T}$. Then, together with
$\hat{\bm{n}}=\left(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\right)^{\rm
T}$, the auto correlation coefficient $\Gamma_{aa}$ can be written as

$\displaystyle\Gamma_{aa}\left(f\right)$
$\displaystyle=\frac{3}{2}\int\frac{d^{2}\hat{\bm{n}}}{4\pi}\left[1-e^{-2\pi
ifL_{a}\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}/v_{\rm
ph}\right)}\right]\left[1-e^{2\pi
ifL_{a}\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}/v_{\rm
ph}\right)}\right]\sum_{A=+,\times}\left[\frac{p_{a}^{i}p_{a}^{j}e_{ij}^{A}\left(\hat{\bm{n}}\right)}{2\left(1+\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}/v_{\rm
ph}\right)}\right]^{2}$ (21)
$\displaystyle=\frac{3}{4}\int_{0}^{\pi}d\theta\>\frac{\sin^{5}\theta\,\sin^{2}\left(\pi
fL_{a}\left(1+\nicefrac{{1}}{{v_{\rm
ph}}}\cos\theta\right)\right)}{\left(1+\nicefrac{{1}}{{v_{\rm
ph}}}\cos\theta\right)^{2}}\,.$ (22)

This integral can be solved exactly,

$\displaystyle\Gamma_{aa}\left(f\right)=3\,v_{\rm ph}^{4}-2\,v_{\rm
ph}^{2}+\frac{v_{\rm
ph}^{5}}{64y^{3}}\>\bigg{\\{}\>24\,y\,x_{+}\,x_{-}\left[\textrm{Ci}\left(x_{+}\right)-\textrm{Ci}\left(x_{-}\right)-2\,\textrm{arcoth}\left(v_{\rm
ph}\right)\right]+3\,x_{+}^{2}x_{-}^{2}\left[\textrm{Si}\left(x_{+}\right)-\textrm{Si}\left(x_{-}\right)\right]$
$\displaystyle+3\left(2+x_{+}\,x_{-}\right)\left(x_{-}\cos x_{+}-x_{+}\cos
x_{-}\right)+6\left[1-\left(3+\nicefrac{{1}}{{v_{\rm
ph}}}\right)yx_{-}\right]\sin x_{+}-6\left[1-\left(3-\nicefrac{{1}}{{v_{\rm
ph}}}\right)yx_{+}\right]\sin x_{-}\>\bigg{\\}}\>\,,$ (23)

where we introduced the shorthand symbols $x_{\pm}$ and $y$, and where Ci and
Si are the cosine and sine integrals,

$x_{\pm}=\frac{2\,y\left(v_{\rm ph}\pm 1\right)}{v_{\rm ph}}\,,\qquad y=\pi
fL_{a}\,,\qquad\textrm{Ci}\left(x\right)=-\int_{x}^{\infty}dt\>\frac{\cos
t}{t}\,,\qquad\textrm{Si}\left(x\right)=\int_{0}^{x}dt\>\frac{\sin t}{t}\,.$
(24)

The auto correlation coefficient $\Gamma_{aa}$ in Eq. (23) represents the main
result of this paper; see Fig. 3 for a plot of $\Gamma_{aa}$ as a function of
the GW phase velocity $v_{\rm ph}$ and for three different values of the
product $fL_{a}$.

We shall now discuss some of the properties of the result in Eq. (23). First,
let us consider the GR limit,

$\Gamma_{aa}^{\rm GR}\left(f\right)=\lim_{v_{\rm ph}\rightarrow
1}\Gamma_{aa}\left(f\right)=1-\frac{3\left[1-\textrm{sinc}\left(4\pi
fL_{a}\right)\right]}{8\left(\pi fL_{a}\right)^{2}}\,.$ (25)

Of course, we could have obtained the same result by simply setting $v_{\rm
ph}$ to the speed of light in Eq. (22),

$\Gamma_{aa}^{\rm
GR}\left(f\right)=\frac{3}{4}\int_{0}^{\pi}d\theta\>\frac{\sin^{5}\theta\,\sin^{2}\left(\pi
fL_{a}\left(1+\cos\theta\right)\right)}{\left(1+\cos\theta\right)^{2}}=1-\frac{3\left[1-\textrm{sinc}\left(4\pi
fL_{a}\right)\right]}{8\left(\pi fL_{a}\right)^{2}}\,,$ (26)

Interestingly enough, this result for $\Gamma_{aa}^{\rm GR}$ retains the
dependence on the GW frequency $f$ and pulsar distance $L_{a}$. In this sense,
$\Gamma_{aa}^{\rm GR}$ differs from the HD result, which just assigns an auto
correlation coefficient of $\Gamma_{aa}^{\rm HD}=1$ to each pulsar. This means
that, in principle, the diagonal terms in the correlation matrix
$\left<R_{a}R_{b}\right>$ do not come with exactly the same correlation
coefficients. Instead, the different distances to the pulsars in the PTA
result in minuscule corrections to the HD result. In practice, however, this
effect is safely negligible. Indeed, because PTAs operate in the long-arm
limit, we always have $fL_{a}\gg 1$, such that to excellent approximation,

$\Gamma_{aa}^{\rm GR}\left(f\right)\quad\overset{fL_{a}\gg
1}{\approx}\quad\Gamma_{aa}^{\rm HD}\left(f\right)=1\,.$ (27)

On top, cosmic variance inherent to the HD curve is guaranteed to overshadow
any effect of order $\left(fL_{a}\right)^{-2}$ in the pulsar correlations [27,
28]. The distinction between $\Gamma_{aa}^{\rm GR}$ and $\Gamma_{aa}^{\rm HD}$
is therefore of little relevance in practice. Nonetheless, we mention that the
same result for $\Gamma_{aa}^{\rm GR}$ was already derived in Ref. [15], which
presents a fully analytical expression for the ORF in GR that retains the full
dependence on the pulsar terms in Eq. (4). In the limit of identical pulsars,
$\xi_{ab}\rightarrow 0$ and $a=b$, this analytical expressions reduces exactly
to the result in Eq. (26).

Next, we observe that the full result in Eq. (23) can be drastically
simplified in the large-$fL_{a}$ limit. Expanding in inverse powers of
$fL_{a}$, we find that $\Gamma_{aa}$ simply scales linearly with $fL_{a}$ to
leading order,

$\Gamma_{aa}\left(f\right)=\Gamma_{aa}^{\rm
LO}\left(f\right)+\mathcal{O}\left(\left(\frac{1}{fL_{a}}\right)^{0}\right)\,,\qquad\Gamma_{aa}^{\rm
LO}\left(f\right)=\Theta\left(-\log_{10}v_{\rm
ph}\right)\,\frac{3}{4}\pi^{2}fL_{a}\,v_{\rm ph}\left(v_{\rm
ph}^{2}-1\right)^{2}\,.$ (28)

This result for $\Gamma_{aa}^{\rm LO}$ confirms that the divergence that
appears in $\Gamma_{aa}^{(\infty)}$ at small angular separations $\xi_{ab}$
and in the regime of subluminal phase velocity $v_{\rm ph}$ is indeed
regulated by a finite-distance effect. Moreover, it illustrates that, as soon
as $v_{\rm ph}$ is no longer close to the speed of light, such that $(v_{\rm
ph}^{2}-1)^{2}$ no longer constitutes a limiting suppression factor,
$\Gamma_{aa}$ becomes of the order of $fL_{a}$ and hence parametrically
enhanced compared to the cross correlation coefficients $\Gamma_{ab}$ at
nonzero angular separations. In order to assess the velocity threshold below
which this enhancement becomes relevant, we solve the condition
$\Gamma_{aa}^{\rm LO}\gtrsim 1$ for $v_{\rm ph}$, which results in

$v_{\rm ph}\lesssim\bar{v}_{\rm
ph}=1-\frac{1+\sqrt{3\pi^{2}fL_{a}}}{3\pi^{2}fL_{a}}\qquad\Rightarrow\qquad\Gamma_{a=b}\gtrsim
1\gtrsim\Gamma_{a\neq b}\,.$ (29)

Figure 4: Modifications of the HD curve in models with a subluminal (left
panel) or superluminal (right panel) GW phase velocity.

For typical GW frequencies and pulsar distances, we have
$fL_{a}\sim\mathcal{O}\left(100\cdots 1000\right)$, which means that
deviations of $v_{\rm ph}$ from the speed of light by a few percent are enough
to cause a large hierarchy between the auto and cross correlations
coefficients, as illustrated in the left panel of Fig. 4. One must therefore
be careful not to be misled by the fact that a decrease in the GW phase
velocity by just a few percent leaves much of the ORF at nonzero $\xi_{ab}$
more or less affected. While it is true that, for $1-v_{\rm
ph}\sim\mathcal{O}\left(0.01\right)$, the changes in the cross correlations of
separate pulsars are practically negligible, it is important to account for
the fact that even a small decrease in $v_{\rm ph}$ can drastically increase
the ratio $\Gamma_{a=b}/\Gamma_{a\neq b}$. We therefore conclude that the ORF
is, in fact, more sensitive to a small decrease in $v_{\rm ph}$ than one may
have naively expected; at the smallest possible value of its argument, i.e.,
at $\xi_{ab}=0$ and for $a=b$, small changes in $v_{\rm ph}$ can have a big
effect. The reason for this enhancement consists in the fact that, for $v_{\rm
ph}<1$, the PTA detector pattern functions in Eq. (5) become resonantly
amplified if the GW propagation direction $\hat{\bm{n}}$ and the pulsar
direction $\hat{\bm{p}}_{a}$ satisfy
$\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}=-v_{\rm ph}$. In GR, this is only possible
if $\hat{\bm{n}}=-\hat{\bm{p}}_{a}$, i.e., for a relative angle between
$\hat{\bm{n}}$ and $\hat{\bm{p}}_{a}$ of $\theta=\pi$. At $\theta=\pi$, the
integrand of the integral in Eq. (26), however, vanishes, such that the
$\left(1+\cos\theta\right)^{2}$ term in the denominator has no chance of
enhancing the overall result. The situation is different in the case of a
subluminal phase velocity, where the $\left(1+\nicefrac{{1}}{{v_{\rm
ph}}}\cos\theta\right)^{2}$ term in Eq. (22) can boost the value of the
integrand around $\theta$ values for which $\cos\theta\simeq-v_{\rm ph}$. The
propagation direction of the corresponding GWs is thus not antiparallel to the
pulsar position, but satisfies instead
$\hat{\bm{n}}\cdot\hat{\bm{p}}_{a}=-v_{\rm ph}$. In passing, we also mention
that, despite the fact that the denominator of the PTA detector pattern
functions reaches zero at $0<\theta<\pi$, the integrand in Eq. (22) actually
never diverges. As discussed at the beginning of Sec. 3, the reason for this
is that all pulsar terms are kept in deriving Eq. (22).

The observation that the auto correlation coefficient $\Gamma_{aa}$ becomes
strongly enhanced for $v_{\rm ph}<1$ has important implications for the
analysis of PTA data. That is, for $v_{\rm ph}$ (much) smaller than the upper
limit in Eq. (29), it becomes a reasonable approximation to simply neglect all
cross correlations and work with

$v_{\rm ph}\lesssim\bar{v}_{\rm
ph}\qquad\Rightarrow\qquad\Gamma_{ab}\left(f\right)\approx\delta_{ab}\,\frac{3}{4}\pi^{2}fL_{a}\,v_{\rm
ph}\left(v_{\rm ph}^{2}-1\right)^{2}\,.$ (30)

In the analysis of PTA data, this expression enters the cross power spectrum
for the timing residuals $R_{a}$, which is closely related to the correlator
in Eq. (2) and constructed from the ORF and the strain power spectrum $S_{h}$,

$S_{ab}\left(f\right)=\Gamma_{ab}\left(f\right)\,\frac{S_{h}\left(f\right)}{6\pi^{2}f^{2}}\,.$
(31)

Then, if we express $S_{h}$ in terms of the characteristic strain amplitude,
$S_{h}=h_{c}^{2}/\left(2f\right)$, and choose a power-law ansatz for $h_{c}$,
we find the following timing-residual cross power spectrum for modified
gravity with $v_{\rm ph}\lesssim\bar{v}_{\rm ph}$,

$h_{c}\left(f\right)=A\left(\frac{f}{f_{\rm ref}}\right)^{\alpha}\,,\qquad
S_{ab}\left(f\right)\approx\delta_{ab}\,fL_{a}\,\frac{v_{\rm ph}(v_{\rm
ph}^{2}-1)^{2}A^{2}}{16f_{\rm ref}^{3}}\left(\frac{f}{f_{\rm
ref}}\right)^{-\gamma}\,,\qquad\gamma=3-2\alpha\,.$ (32)

This template for $S_{ab}$ shares some similarity with the usual template for
a common-spectrum spatially uncorrelated red-noise (CURN) process, in the
sense that $S_{ab}\propto\delta_{ab}$. The dependence on $fL_{a}$, however,
represents a departure from a universal “common spectrum” for all pulsars.
Moreover, the large values of the auto correlation coefficient will cause a
suppression in the inferred amplitude $A$ compared to the CURN case.

##  5 Massive gravity

Finally, we turn to superluminal phase velocities (subluminal group
velocities). In this case, the leading-order result in Eq. (28), which simply
vanishes at $v_{\rm ph}\geq 1$, no longer suffices. We therefore need to
include the next-to-leading order contribution to $\Gamma_{aa}$ when expanding
the full result in Eq. (23) in inverse powers of $fL_{a}$,

$\Gamma_{aa}^{\rm NLO}\left(f\right)=\frac{3}{2}\,v_{\rm ph}^{3}\left(v_{\rm
ph}^{2}-1\right)\ln\left(\frac{\left|v_{\rm ph}-1\right|}{v_{\rm
ph}+1}\right)+3\,v_{\rm ph}^{4}-2\,v_{\rm ph}^{2}\,.$ (33)

Together with the leading-order result in Eq. (28), we thus obtain

$\displaystyle\Gamma_{aa}\left(f\right)$ $\displaystyle=\Gamma_{aa}^{\rm
LO}\left(f\right)+\Gamma_{aa}^{\rm
NLO}\left(f\right)+\mathcal{O}\left(\left(\frac{1}{fL_{a}}\right)^{1}\right)\,,$
(34) $\displaystyle\Gamma_{aa}^{\rm LO}\left(f\right)$
$\displaystyle=\Theta\left(-\log_{10}v_{\rm
ph}\right)\,\frac{3}{4}\pi^{2}fL_{a}\,v_{\rm ph}\left(v_{\rm
ph}^{2}-1\right)^{2}\,,$ (35) $\displaystyle\Gamma_{aa}^{\rm
NLO}\left(f\right)$ $\displaystyle=\frac{3}{2}\,v_{\rm ph}^{3}\left(v_{\rm
ph}^{2}-1\right)\ln\left(\frac{\left|v_{\rm ph}-1\right|}{v_{\rm
ph}+1}\right)+3\,v_{\rm ph}^{4}-2\,v_{\rm ph}^{2}\,.$ (36)

In the GR limit, $\Gamma_{aa}^{\rm LO}\rightarrow 0$ and $\Gamma_{aa}^{\rm
NLO}\rightarrow 1$, in accord with the fact that $\Gamma_{aa}^{\rm GR}\approx
1$, up to corrections of order $\left(fL_{a}\right)^{-2}$. The above estimate
of the exact auto correlation coefficient in Eq. (23) is valid at $v_{\rm
ph}\geq 1$ and at $v_{\rm ph}<1$; see Fig. 3, in which we compare our three
results for $\Gamma_{aa}$: the full result in Eq. (23), the leading-order
expression in Eq. (28), and the combination of the leading-order and next-to-
leading-order terms in Eq. (34). In view of this plot, we conclude that
$\Gamma_{aa}^{\rm LO}$ serves as an excellent approximation of $\Gamma_{aa}$
at subluminal $v_{\rm ph}$, while the combination $\Gamma_{aa}^{\rm
LO}+\Gamma_{aa}^{\rm NLO}$ provides an accurate approximation of $\Gamma_{aa}$
at subluminal and superluminal $v_{\rm ph}$.

In massive gravity, one may prefer to express the dependence of the ORF on the
GW speed in terms of the group velocity rather than the phase velocity,
$v_{\rm gr}=1/v_{\rm ph}$. Our result for $\Gamma_{aa}^{\rm NLO}$ in terms of
$v_{\rm gr}$ then reads

$\Gamma_{aa}^{\rm NLO}\left(f\right)=\frac{1}{2\,v_{\rm
gr}^{5}}\left[6\,v_{\rm gr}-4\,v_{\rm gr}^{3}+3\left(v_{\rm
gr}^{2}-1\right)\ln\left(\frac{1+v_{\rm gr}}{1-v_{\rm gr}}\right)\right]\,.$
(37)

This expression does not depend on the GW frequency $f$ explicitly, but only
implicitly via the group velocity,

$v_{\rm
gr}=\frac{1}{\sqrt{1+\left(m_{g}/k\right)^{2}}}=\sqrt{1-\left(\frac{f_{g}}{f}\right)^{2}}\,,$
(38)

where we used that $\omega=2\pi f=\sqrt{m_{g}^{2}+k^{2}}$ in massive gravity,
and where $f_{g}$ is the graviton Compton frequency,

$f_{g}=\frac{m_{g}}{2\pi}\simeq
24.18\,\textrm{nHz}\left(\frac{m_{g}}{10^{-22}\,\textrm{eV}}\right)\,.$ (39)

The result in Eq. (37) can be compared to the known analytical expression for
the ORF in massive gravity [13, 19, 21], which follows from solving the
integral in Eq. (4) after discarding the pulsar terms inside the brackets,

$\displaystyle\Gamma_{\rm a\neq b}^{\rm
mass}\left(\xi_{ab},f\right)=\frac{1}{16\,v_{\rm gr}^{5}}\left[2v_{\rm
gr}\left(3+\left(6-5v_{\rm gr}^{2}\right)\delta\right)-6\left[1+\delta+v_{\rm
gr}^{2}\left(1-3\delta\right)\right]\ln\left(\frac{1+v_{\rm gr}}{1-v_{\rm
gr}}\right)-\frac{3A}{B}\,\ln C\right]\,,$ (40) $\displaystyle A=1+2v_{\rm
gr}^{2}\left(1-2\delta\right)-v_{\rm
gr}^{4}\left(1-2\delta^{2}\right)\,,\qquad
B=\sqrt{\left(1-\delta\right)\left(2-v_{\rm
gr}^{2}\left(1+\delta\right)\right)}\,,\qquad C=\frac{A-2v_{\rm
gr}\left(1-v_{\rm gr}^{2}\delta\right)B}{\left(v_{\rm
gr}^{2}-1\right)^{2}}\,,$ (41)

where the dependence on the angular separation is encoded in
$\delta=\cos\xi_{ab}$. In the limit $\xi_{ab}\rightarrow 0$, we find

$\lim_{\xi_{ab}\rightarrow 0^{+}}\Gamma_{\rm a\neq b}^{\rm
mass}\left(\xi_{ab},f\right)=\frac{1}{4\,v_{\rm gr}^{5}}\left[6\,v_{\rm
gr}-4\,v_{\rm gr}^{3}+3\left(v_{\rm gr}^{2}-1\right)\ln\left(\frac{1+v_{\rm
gr}}{1-v_{\rm gr}}\right)\right]=\frac{1}{2}\,\Gamma_{aa}^{\rm
NLO}\left(f\right)\,,$ (42)

which equals exactly one half of our result for $\Gamma_{aa}^{\rm NLO}$. We
thus conclude that the full expression for the ORF in massive gravity, valid
for both separate and identical pulsars, requires the same $1+\delta_{ab}$
prefactor as in GR,

$\Gamma_{\rm ab}^{\rm
mass}\left(\xi_{ab},f\right)=\frac{1+\delta_{ab}}{16\,v_{\rm
gr}^{5}}\left[2v_{\rm gr}\left(3+\left(6-5v_{\rm
gr}^{2}\right)\delta\right)-6\left[1+\delta+v_{\rm
gr}^{2}\left(1-3\delta\right)\right]\ln\left(\frac{1+v_{\rm gr}}{1-v_{\rm
gr}}\right)-\frac{3A}{B}\,\ln C\right]\,.$ (43)

To the best of our knowledge, the prefactor of $1+\delta_{ab}$ has thus far
been overlooked in earlier studies. Just like in GR, the extra factor of $2$
for $a=b$ originates from the auto correlation of the pulsar term in Eq. (13).

##  6 Conclusions

PTA searches for GWs require a precise understanding of the ORF $\Gamma_{ab}$,
which is sensitive to the presence of new physics in the gravity sector. For
instance, the HD curve in GR receives corrections if the GW phase velocity
$v_{\rm ph}$ does not coincide with the speed of light. In the literature, the
ORF had been previously calculated for such scenarios of modified gravity for
both for $v_{\rm ph}<1$ and $v_{\rm ph}>1$ by various authors. Most existing
analyses, however, only consider the case of separate pulsars ($a\neq b$) and
neglect the case of identical pulsars ($a=b$). In this paper, we therefore
revisited the computation of the ORF in scenarios with a nonstandard phase
velocity, supplementing the known expressions for the cross correlation
coefficients $\Gamma_{a\neq b}$ by a careful analytical evaluation of the auto
correlation coefficients $\Gamma_{a=b}$; see Fig. (4) for a visual
illustration of our main results.

In the GR limit, $v_{\rm ph}\rightarrow 1$, we recover the standard HD result,
$\Gamma_{aa}^{\rm HD}=1$, at leading order, but also find (numerically
negligible) subleading corrections of the order of $(fL_{a})^{-2}$, which
introduce (at least, in principle) a dependence on the pulsar distance $L_{a}$
in the diagonal entries of $\left<R_{a}R_{b}\right>$. For subluminal $v_{\rm
ph}$, we confirm that the divergence of the ORF in the limit of small angular
separations that had been observed in earlier works is unphysical and instead
regulated by the finite distances to the pulsars in the PTA. At leading order,
we notably find that the auto correlation coefficient scales linearly with
$fL_{a}$; see Eq. (28). We therefore conclude that, for a sufficiently small
phase velocity, $\Gamma_{a=b}$ becomes parametrically enhanced compared to
$\Gamma_{a\neq b}$, which means that the cross correlations among pulsars
become insignificant to good approximation in this model.

Meanwhile, we remark that the known expression for $\Gamma_{ab}$ in Eq. (19)
still suffers from a certain degree of uncertainty at very small values of
$\xi_{ab}$. The point is: While we derived the finite auto correlation
coefficient for subluminal $v_{\rm ph}$ from first principles, we did not
study how the approximate resummation of higher multipole terms in the
Legrendre decomposition (or any other way of computing the ORF, for that
matter) needs to be modified such that $\Gamma_{ab}$ does indeed approach
$\Gamma_{aa}$ in the limit $\xi_{ab}\rightarrow 0$ and for $a=b$. In other
words, our calculation tells us the value that the ORF must take at
$\xi_{ab}=0$ and for $a=b$; we, however, still do not have a complete
understanding of how this value is actually reached as $\xi_{ab}\rightarrow 0$
and $L_{a}\rightarrow L_{b}$. For instance, in GR and for superluminal $v_{\rm
ph}$, we know that the ORF first reaches a finite value,
$\lim_{\xi_{ab}\rightarrow 0^{+}}\Gamma_{ab}$, which then needs to be
multiplied by a factor of $1+\delta_{ab}$, in order to produce the final value
of $\Gamma_{aa}$. For subluminal $v_{\rm ph}$, the exact relation between
$\lim_{\xi_{ab}\rightarrow 0^{+}}\Gamma_{ab}$ and $\Gamma_{aa}$, however,
remains unclear at present. For practical applications, this limitation is of
little relevance, though, as long as the angular separations between all
pulsars in the PTA are bounded from below. We checked, e.g., that the
expression in Eq. (19) is insensitive to changes in $fL_{a}$ for angular
separations of at least $\xi_{ab}=\pi/100$. We therefore expect that finite-
distance effects can be neglected for all but the tiniest angular separation,
which, however, are not realized in typical PTAs. In summary, this means that,
in models with $v_{\rm ph}<1$ and for PTAs in which no two pulsars are
extremely close together in the sky, it is justified to use our result for
$\Gamma_{aa}$ to describe the auto correlations of individual pulsars, in
combination with Eq. (19) to compute the cross correlations among pulsars. At
the same time, the exact relation between $\lim_{\xi_{ab}\rightarrow
0^{+}}\Gamma_{ab}$ and $\Gamma_{aa}$ for $v_{\rm ph}<1$ remains an interesting
conceptual question for future work.

Finally, for superluminal $v_{\rm ph}$, which is relevant in the case of
massive gravity, we obtain the known analytical result for $\Gamma_{aa}$, up
to a factor of $2$, which had been overlooked in earlier studies. In
conclusion, the results presented in this paper now pave the way for fitting
modified-gravity models with sub- or superluminal phase velocity to PTA data.
We will present the result of such fits to recent PTA data in an upcoming
publication.

### Acknowledgements

We thank Qiuyue Liang and Mark Trodden for helpful discussions. This work was
supported by the Deutsche Forschungsgemeinschaft under Germany’s Excellence
Strategy – EXC 2121 Quantum Universe – 390833306. The work of K. S. and T. S.
is supported by Deutsche Forschungsgemeinschaft (DFG) through the Research
Training Group (Graduiertenkolleg) 2149: Strong and Weak Interactions – from
Hadrons to Dark Matter.

## References

  * [1] S. R. Taylor, “The Nanohertz Gravitational Wave Astronomer,” arXiv:2105.13270 [astro-ph.HE].
  * [2] NANOGrav Collaboration, G. Agazie et al., “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,” Astrophys. J. Lett. 951 no. 1, (2023) L8, arXiv:2306.16213 [astro-ph.HE].
  * [3] EPTA, InPTA: Collaboration, J. Antoniadis et al., “The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals,” Astron. Astrophys. 678 (2023) A50, arXiv:2306.16214 [astro-ph.HE].
  * [4] D. J. Reardon et al., “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 951 no. 1, (2023) L6, arXiv:2306.16215 [astro-ph.HE].
  * [5] H. Xu et al., “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,” Res. Astron. Astrophys. 23 (2023) 075024, arXiv:2306.16216 [astro-ph.HE].
  * [6] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford University Press, 2007.
  * [7] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology.  Oxford University Press, 3, 2018.
  * [8] J. D. Romano and B. Allen, “Answers to frequently asked questions about the pulsar timing array Hellings and Downs curve,” arXiv:2308.05847 [gr-qc].
  * [9] Q. Liang, M.-X. Lin, and M. Trodden, “A test of gravity with Pulsar Timing Arrays,” JCAP 11 (2023) 042, arXiv:2304.02640 [astro-ph.CO].
  * [10] K. Lee, F. A. Jenet, R. H. Price, N. Wex, and M. Kramer, “Detecting massive gravitons using pulsar timing arrays,” Astrophys. J. 722 (2010) 1589–1597, arXiv:1008.2561 [astro-ph.HE].
  * [11] S. J. Chamberlin and X. Siemens, “Stochastic backgrounds in alternative theories of gravity: overlap reduction functions for pulsar timing arrays,” Phys. Rev. D 85 (2012) 082001, arXiv:1111.5661 [astro-ph.HE].
  * [12] W. Qin, K. K. Boddy, and M. Kamionkowski, “Subluminal stochastic gravitational waves in pulsar-timing arrays and astrometry,” Phys. Rev. D 103 (2021) 024045, arXiv:2007.11009 [gr-qc].
  * [13] Q. Liang and M. Trodden, “Detecting the stochastic gravitational wave background from massive gravity with pulsar timing arrays,” Phys. Rev. D 104 (2021) 084052, arXiv:2108.05344 [astro-ph.CO].
  * [14] J. M. Ezquiaga, W. Hu, M. Lagos, and M.-X. Lin, “Gravitational wave propagation beyond general relativity: waveform distortions and echoes,” JCAP 11 no. 11, (2021) 048, arXiv:2108.10872 [astro-ph.CO].
  * [15] Y. Hu, P.-P. Wang, Y.-J. Tan, and C.-G. Shao, “Full analytic expression of overlap reduction function for gravitational wave background with pulsar timing arrays,” Phys. Rev. D 106 no. 2, (2022) 024005, arXiv:2205.09272 [gr-qc].
  * [16] R. C. Bernardo and K.-W. Ng, “Stochastic gravitational wave background phenomenology in a pulsar timing array,” Phys. Rev. D 107 (2023) 044007, arXiv:2208.12538 [gr-qc].
  * [17] K. Schumacher, N. Yunes, and K. Yagi, “Gravitational wave polarizations with different propagation speeds,” Phys. Rev. D 108 no. 10, (2023) 104038, arXiv:2308.05589 [gr-qc].
  * [18] N. Anil Kumar and M. Kamionkowski, “All the Pretty Overlap Reduction Functions,” arXiv:2311.14159 [astro-ph.CO].
  * [19] Y.-M. Wu, Z.-C. Chen, and Q.-G. Huang, “Search for stochastic gravitational-wave background from massive gravity in the NANOGrav 12.5-year dataset,” Phys. Rev. D 107 no. 4, (2023) 042003, arXiv:2302.00229 [gr-qc].
  * [20] R. C. Bernardo and K.-W. Ng, “Constraining gravitational wave propagation using pulsar timing array correlations,” Phys. Rev. D 107 no. 10, (2023) L101502, arXiv:2302.11796 [gr-qc].
  * [21] Y.-M. Wu, Z.-C. Chen, Y.-C. Bi, and Q.-G. Huang, “Constraining the graviton mass with the NANOGrav 15 year data set,” Class. Quant. Grav. 41 no. 7, (2024) 075002, arXiv:2310.07469 [astro-ph.CO].
  * [22] R. C. Bernardo and K.-W. Ng, “Beyond the Hellings-Downs curve: Non-Einsteinian gravitational waves in pulsar timing array correlations,” arXiv:2310.07537 [gr-qc].
  * [23] Y.-C. Bi, Y.-M. Wu, Z.-C. Chen, and Q.-G. Huang, “Constraints on the velocity of gravitational waves from the NANOGrav 15-year data set,” Phys. Rev. D 109 no. 6, (2024) L061101, arXiv:2310.08366 [astro-ph.CO].
  * [24] Q. Liang, I. Obata, and M. Sasaki, “Testing Gravity with Frequency-Dependent Overlap Reduction Function in Pulsar Timing Array,” arXiv:2405.11755 [astro-ph.CO].
  * [25] R. W. Hellings and G. S. Downs, “Upper limits on the isotropic gravitational radiation background from pulsar timing analysis,” Astrophys. J. Lett. 265 (1983) L39–L42.
  * [26] J. Gair, J. D. Romano, S. Taylor, and C. M. F. Mingarelli, “Mapping gravitational-wave backgrounds using methods from CMB analysis: Application to pulsar timing arrays,” Phys. Rev. D 90 (2014) 082001, arXiv:1406.4664 [gr-qc].
  * [27] B. Allen, “Variance of the Hellings-Downs correlation,” Phys. Rev. D 107 no. 4, (2023) 043018, arXiv:2205.05637 [gr-qc].
  * [28] B. Allen and J. D. Romano, “Hellings and Downs correlation of an arbitrary set of pulsars,” Phys. Rev. D 108 no. 4, (2023) 043026, arXiv:2208.07230 [gr-qc].