# Answers to frequently asked questions about the pulsar timing array
Hellings and Downs correlation curve

Joseph D. Romano<EMAIL_ADDRESS>Department of Physics and
Astronomy, Texas Tech University, Box 41051, Lubbock, TX 79409-1051, USA
Department of Physics and Astronomy, University of Texas Rio Grande Valley,
One West University Boulevard, Brownsville, TX 78520, USA  Bruce Allen
<EMAIL_ADDRESS>Max Planck Institute for Gravitational Physics (Albert
Einstein Institute), Leibniz Universität Hannover, Callinstrasse 38, D-30167
Hannover, Germany

###### Abstract

We answer frequently asked questions (FAQs) about the Hellings and Downs
correlation curve—the “smoking-gun” signature that pulsar timing arrays (PTAs)
have detected gravitational waves (GWs). Many of these questions arise because
of intuition based on how ground-based interferometers like LIGO respond to
GWs. These have arms that are short (km scale) compared to the wavelengths of
the GWs that they detect (hundreds to thousands of km). In contrast, PTAs
respond to GWs whose wavelengths (tens of light-years) are much shorter than
their arms (a typical PTA pulsar is thousands of light-years from Earth). To
elucidate this, we calculate the exact response of a “one-arm, one-way”
detector to a passing GW, and compare it in the “short-arm” (LIGO-like) and
“long-arm” (PTA) limits. This provides qualitative and quantitative answers to
many questions about the Hellings and Downs correlation. The resulting “FAQ
sheet” should help in understanding the “evidence for GWs” recently announced
by several PTA collaborations.

##  I Introduction

Four pulsar timing array (PTA) collaborations recently announced the results
of their latest searches for correlated low-frequency ($\sim\\!10^{-9}~{}{\rm
Hz}$) gravitational waves (GWs). Their findings span the range from “weak
evidence” [Reardon et al., 2023] to “some evidence” [Xu et al., 2023] to
“evidence/marginal evidence” [Antoniadis et al., 2023] to “compelling
evidence” [Agazie et al., 2023].

These observations can be challenging to understand, because the “signal” is
unlike the short duration ($\sim\\!0.25~{}{\rm s}$) GW “chirp signal” that
LIGO observed on 14 Sep 2015, from the final inspiral and merger of two
$\sim\\!30M_{\odot}$ black holes. Instead of such short (1 s) audio-frequency
(100 Hz) deterministic GW signals, PTAs search for persistent (years) low-
frequency (nHz) stochastic GW signals. These are too feeble to observe
directly, but should produce correlated perturbations in the arrival times of
pulses from an array of galactic millisecond pulsars.

Although the current PTA observations are not significant enough to allow one
to identify the source of the correlated signal, a credible candidate is the
superposition of GWs from hundreds of thousands of supermassive
($\sim\\!10^{9}M_{\odot}$) black-hole (SMBH) binaries orbiting one another in
the centers of merging galaxies. While any one of those binaries would produce
a deterministic almost-constant-frequency signal, the incoherent sum of many
such sources creates a stochastic hubub of GWs.

PTAs search for correlations which are predicted to follow the so-called
“Hellings and Downs” curve. This correlation curve is named after Ron Hellings
and George Downs, who first derived it in 1983 [Hellings, 1989]. The idea is
that pulsars act as high-quality clocks, producing radio pulses that arrive at
Earth with very predictable and regular arrival times. GWs induce slight
perturbations in these arrival time, and these perturbations (called “timing
residuals”) are correlated between different radio pulsars.

The Hellings and Downs curve is a plot of the expected correlation in the
timing residuals from a pair of pulsars, as a function of their angular
separation on the sky. It is a prediction of general relativity (GR) for a
gravitational-wave background (GWB) that is unpolarized and isotropic. Figure
1 shows how the timing residuals in NANOGrav’s 15-year dataset [Agazie et al.,
2023] appear to match this predicted Hellings and Downs curve.

Figure 1: The spatial correlations observed in the pulsar timing residuals for
the NANOGrav 15-year dataset [Agazie et al., 2023] are shown in blue; the
Hellings and Downs prediction is shown in black. The blue points and error
bars are optimally weighted averages of approximately 150 pulsar-pair
correlations in each angular separation bin, which take into account
covariances between the correlations induced by the GWB itself. (Figure taken
from [Agazie et al., 2023].)

The goal of this paper is to provide qualitative and quantitative answers to
common questions about the properties of the Hellings and Downs correlation
curve. We believe that such questions frequently arise because much of our
intuition regarding the response of a detector to a passing GW was developed
for ground-based GW detectors such as LIGO and Virgo. Not all of this
intuition carries over to the case of PTAs. This is because the distance
between the Earth and a pulsar (one “arm” of a PTA detector) is typically
hundreds to thousands of times longer than the wavelengths (tens of light-
years) of the GWs that PTAs are trying to detect. In contrast, ground-based
interferometers have km-long arms that are much shorter than the wavelengths
(100’s to 1000’s of km) of the GWs that these detectors are sensitive to. In
short, PTA detectors operate in the “long-arm” (or short-wavelength) limit,
whereas ground-based laser interferometers operate in the “short-arm” (or
long-wavelength) limit.

The remainder of this paper is organized as follows: We begin immediately with
the questions and answers in Section II, developing the formal treatment on an
“as needed” basis. The most important concepts are: (i) the exact (timing
residual) response of a “one-arm, one-way” detector to a passing GW, and (ii)
the expected correlation of two different detectors, expressed in terms of
their individual response functions. Limiting expressions for these then give
the response and correlation for “short-arm” LIGO-like detectors, and “long-
arm” PTA pulsar detectors.

For pedagogical purposes, we also suggest exercises for the reader, but don’t
provide solutions. Instead, for the more involved exercises, we point to
specific relevant materials in the published literature.

Although the “Hellings-and-Downs correlation FAQ sheet” is the main content of
the paper, we conclude with a brief summary in Section III.

##  II Frequently asked questions

Q1: Why is the Hellings and Downs correlation normalized at zero angular
separation to $1/3$ in some papers (e.g., [Hellings, 1989; Allen, 2023]) and
to $1/2$ in this FAQ sheet and other papers (e.g., [Anholm et al., 2009;
Agazie et al., 2023])? Does it matter?

A1: This normalization difference is an overall scale: multiply the Hellings
and Downs curve by a factor of $2/3$ to go in one direction, or by $3/2$ to go
in the other. This choice of normalization is completely arbitrary and makes
no difference: any physically observable quantity (for example, the expected
time-averaged product of the timing residuals of two pulsars) is independent
of this normalization.

The reason why $1/3$ and $1/2$ have appeared in the literature is historical.
The Hellings and Downs correlation was first computed (i) as the average
product of the timing residual responses of two pulsars to a unit-amplitude GW
source, averaged over the sphere of GW source directions. This average
introduces a factor of $1/4\pi$ in (23). Later, the Hellings and Downs
correlation was re-interpreted as (ii) a correlation coefficient for the
timing residuals of two pulsars, responding to the same GW source. This
introduces a factor of $3/8\pi$ in (23), which is the normalization choice
used in this paper. In case (i) we get $1/3$ if the pulsars are distinct but
close in sky direction, whereas in case (ii) we get $1/2$ if the pulsars are
distinct but close on the sky, and $1$ (an auto-correlation) for two identical
pulsars.

Q2: Why does the Hellings and Downs correlation have different values at
$0^{\circ}$ and $180^{\circ}$ if the quadrupolar deformation of space produced
by a passing GW affects two test masses $180^{\circ}$ apart in exactly the
same way?

A2: This simple question gets to the heart of the difference between “short-
arm” LIGO-like detectors and “long-arm” PTA detectors. Recall that a short-arm
detector has $fL/c=L/\lambda\ll 1$, where $L$ is the length of the detector
arm and $\lambda\equiv c/f$ is the GW wavelength. A long-arm detector has
$L/\lambda\gg 1$. (Here and throughout, $c$ denotes the speed of light, which
is also the propagation speed of GWs.)

A quick, qualitative answer to this question is that for normal incidence of a
GW on either a short-arm or long-arm detector, the response of two pulsars in
opposite directions will be identical. But for non-normal incidence on a long-
arm detector, the response of the two pulsars will be different from one
another, because their radio pulses encounter different phases of the GW as
they propagate from the pulsars to Earth. Pulses from one of the pulsars will
be (partially) “swimming with the current” of the GW as they travel from the
pulsar to Earth. In contrast, pulses from the other pulsar will “swim against
the current” or “fight upstream” and pass through more cycles of the GW. Since
a GWB is composed of GWs coming from many different directions on the sky, the
correlated response of two pulsars oriented in opposite directions will differ
from that for two pulsars oriented along the same direction. Several
illustrative examples are shown in Figure 2.

Figure 2: A monochromatic plane GW propagating in direction $\hat{k}$ passes
by six “one-arm, one-way” detectors with arm length $L$. These have different
propagation directions $\hat{u}$ for the electromagnetic waves that monitor
the distance between two “test masses” (dots). For PTAs, $\hat{u}$ is the
propagation direction of a radio pulse from pulsar to Earth, and
$\hat{p}=-\hat{u}$ is the direction to the pulsar on the sky. All six
detectors have $L=5\lambda$; the solid horizontal lines indicates the maxima
of the GW at one moment in time. After reading this section, the reader should
be able to show that: for case (a), both detectors have zero response to the
GW because the radio pulses propagate parallel to the propagation direction of
the GW; for case (b), both detectors have zero response to the GW because the
radio pulses pass through an integer number of GW cycles during their transit
from pulsar to Earth; for case (c) with
$\hat{k}\cdot\hat{u}_{1}=-\hat{k}\cdot\hat{u}_{2}=1/2$, the radio pulse moving
in direction $\hat{u}_{1}$ (partially) “surfs” the GW so that this detector
has a $3\times$ larger response to the GW than for the detector with pulse
direction $\hat{u}_{2}$.

In order to give a more complete and quantitative answer to this question, we
need to understand how a generic “one-arm, one-way” detector responds to a
passing GW. We will then combine the resulting response functions to obtain
the expected correlation, obtaining the Hellings and Downs curve. For
presentation purposes, we split the following discussion into a few main
parts.

###  II.1 Calculating the timing residual response

The geometry of a one-arm detector, where the electromagnetic signal (from a
pulsar or laser) travels in just one direction along that arm, is shown in
Figure 3. Such “beam detectors” use electromagnetic radiation to monitor the
separation of two or more freely-falling masses. Laser interferometers (both
ground- and space-based), spacecraft Doppler tracking, and PTAs are all
examples of beam detectors.

Figure 3: Geometry of a one-arm, one-way beam detector. The GW propagates in
the $\hat{k}$ direction; the electromagnetic wave (e.g., a radio pulse from a
pulsar, or laser light traveling between two mirrors) propagates in the
$\hat{u}$ direction, opposite the direction from Earth to the pulsar, or from
the second test mass to the first test mass, $\hat{p}=-\hat{u}$. Here, $L$ is
the “arm” length of the detector.

As a GW passes by the detector, it will “stretch” and “squeeze” space
perpendicular to its direction of propagation $\hat{k}$ (see the upper right
panel of Figure 4). This quadrupolar deformation of space causes the
electromagnetic signal (which we will treat generically here as a photon) to
arrive either slightly earlier or later than it would in the absence of a GW.
The size of the advance or delay in the photon arrival time will depend on the
orientation of the detector (defined by the unit vector $\hat{u}$) relative to
$\hat{k}$, and the length of the arm $L$ relative to the wavelength of the GW,
i.e., $L/\lambda$ or equivalently $fL/c$. We will denote this change in the
arrival time by $\Delta T(t)$, which is a function of the time $t$ that the
observation is made.

An explicit formula for $\Delta T(t)$ involves the metric perturbations
$h_{ij}(t,\vec{x})$, which contain all information about the GWs. In fact,
$\Delta T(t)$ only depends upon the metric perturbations along the spacetime
path followed by the photon as it moves through the GW detector, arriving at
time $t$:

$\Delta T(t)=\frac{1}{2c}u^{i}u^{j}\int_{0}^{L}{\rm
d}s\>h_{ij}(\tau(s),\vec{x}(s))\,.$ (1)

From this equation, you can see that the photon moving through the detector is
“sampling” the GW field $h_{ij}$ along its spacetime path
$(\tau(s),\vec{x}(s))$.

That spacetime path is described by the parametric representation
$\tau(s)=t+(s-L)/c$ and $\vec{x}(s)=\vec{r}_{1}+s\hat{u}$. The parameter $s$
starts at $0$ (photon emission) and increases smoothly to $L$ (photon
reception). This spacetime path is a straight line in 4-dimensional spacetime,
starting at the point $(t-L/c,\vec{r}_{1})$ and ending at the point
$(t,\vec{r}_{2})$. The reader should verify that the magnitude of the velocity
vector ${\rm d}\vec{x}/{\rm d}\tau=({\rm d}\vec{x}/{\rm d}s)({\rm d}s/{\rm
d}\tau)$ is the speed of light.

We will use $t_{1}\equiv t-L/c$ to denote the time that the photon left the
source at $\vec{r}_{1}$, and $t_{2}\equiv t$ to be the time that the photon
was received at the receiver at $\vec{r}_{2}\equiv\vec{r}_{1}+L\hat{u}$. Note
that we do not need to include any corrections to the straight-line path for
the photon given above, as the metric perturbations are already first-order
small and we can ignore all second and higher-order terms in the calculation.

To do the integral we need an explicit expression for the metric
perturbations. Recall that GWs are time-varying perturbations to the geometry
of space-time, which propagate away from the source at the speed of light, see
e.g., [Misner et al., 1973; Hartle, 2002]. In transverse-traceless coordinates
$(t,\vec{x})\equiv(t,x^{i})$, where $i=1,2,3$, the metric perturbations
corresponding to a plane wave have two degrees of freedom, corresponding to
the amplitudes of the plus ($+$) and cross ($\times$) polarizations of the GW
(Figure 4).

Figure 4: Left panel: Coordinate system and unit vectors used in the plane
wave expansion. Right panel (top): two orthogonal polarizations of a GW. A
circular ring of test particles in the plane orthogonal to the propagation
direction $\hat{k}$ of the wave are alternately deformed into ellipses, as
space is “squeezed” and “stretched” by the passing of the wave. Right panel
(bottom): an alternative choice of two orthogonal polarizations may be
obtained by rotating $\hat{l}$ and $\hat{m}$ through angle $\psi$ to obtain
new vectors $\hat{l}^{\prime}$ and $\hat{m}^{\prime}$.

The metric perturbations for the most general GW field can be written as a
superposition of plane waves:

$h_{ij}(t,\vec{x})=\int_{-\infty}^{\infty}{\rm d}f\>\int{\rm
d}\hat{k}\>\sum_{A=+,\times}h_{A}(f,\hat{k})e^{A}_{ij}(\hat{k})e^{i2\pi
f(t-\hat{k}\cdot\vec{x}/c)}\,,$ (2)

where $f$ denotes the frequency of the component waves, $\hat{k}$ their
direction of propagation, and $A=+,\times$ their polarization. Note that this
plane wave expansion can represent any possible pattern of GWs: it does not
have to be from a stochastic background, and might be from a single
deterministic source, or any combination of such (deterministic and/or
stochastic) sources. We only need the correct functions $h_{A}(f,\hat{k})$.

There is only one caveat: this equation does not correctly represent the
metric perturbations close to a GW source. But that doesn’t matter here. The
first significant sources of GWs for these detectors are many of orders of
magnitude farther way than $L$, so in the relevant region of spacetime around
the detectors there are no sources, and (2) holds.

In thinking about (2), it is helpful to imagine that the integral over
directions $\hat{k}$ is really a discrete sum, with one term for each
(distant) GW source. The direction to a particular GW source is given by
$\hat{n}\equiv-\hat{k}$. The quantities $e^{A}_{ij}(\hat{k})$ with
$A=+,\times$ are polarization tensors, defined by

$\displaystyle e_{ij}^{+}(\hat{k})$
$\displaystyle\equiv\hat{l}_{i}\hat{l}_{j}-\hat{m}_{i}\hat{m}_{j}\,,\qquad
e_{ij}^{\times}(\hat{k})$
$\displaystyle\equiv\hat{l}_{i}\hat{m}_{j}+\hat{m}_{i}\hat{l}_{j}\,,$ (3)

where $\hat{l}$, $\hat{m}$ are any two orthogonal unit vectors in the plane
orthogonal to $\hat{k}$. Conventionally, for stochastic background analyses,
we take

$\displaystyle\hat{k}$ $\displaystyle\equiv-$
$\displaystyle\hat{x}\,\sin\theta\cos\phi$
$\displaystyle-\hat{y}\,\sin\theta\sin\phi$
$\displaystyle-\hat{z}\,\cos\theta\,,$ (4) $\displaystyle\hat{l}$
$\displaystyle\equiv$ $\displaystyle\hat{x}\,\sin\phi$
$\displaystyle-\hat{y}\,\cos\phi\,,$ $\displaystyle\hat{m}$
$\displaystyle\equiv-$ $\displaystyle\hat{x}\,\cos\theta\cos\phi$
$\displaystyle-\hat{y}\,\cos\theta\sin\phi$
$\displaystyle+\hat{z}\,\sin\theta\,.$

The reader should verify that these three vectors are orthonormal, meaning
that each has unit length and is perpendicular to the other two. In fact,
$\hat{k}=-\hat{r}$, $\hat{l}=-\hat{\phi}$, $\hat{m}=-\hat{\theta}$, where
$\hat{r}$, $\hat{\theta}$, $\hat{\phi}$ are the standard spherical polar
coordinate unit vectors.

For a given GW propagation direction $\hat{k}$, one can form another
(equivalent) orthonormal set by rotating the vectors $\hat{l}$ and $\hat{m}$
by an angle $\psi$ (called the polarization angle) in the plane orthogonal to
$\hat{k}$, as shown in the bottom right panel of Figure 4. We leave it as an
exercise for the reader to show that such a rotation of $\hat{l}$ and
$\hat{m}$ leads to a new set of polarization tensors

$\displaystyle\begin{bmatrix}e^{\prime}{}^{+}_{ij}\\\
e^{\prime}{}^{\times}_{ij}\end{bmatrix}=\begin{bmatrix}\cos 2\psi&\sin
2\psi\\\ -\sin 2\psi&\cos 2\psi\end{bmatrix}\begin{bmatrix}e^{+}_{ij}\\\
e^{\times}_{ij}\end{bmatrix}\,.$ (5)

This change of polarization basis is a spin-two gauge transformation. It leads
to different values of the expansion coefficients $h_{A}(f,\hat{k})$, but does
not change the metric perturbations $h_{ij}(t,\vec{x})$.

Another useful exercise is to define a complex polarization tensor
$e_{ij}\equiv e^{+}_{ij}-ie^{\times}_{ij}$, and to show that under (5) it
transforms to $e^{\prime}_{ij}=\exp(2i\psi)e_{ij}$. In fact, the “2” that
appears in this transformation is the reason why the field is said to have
“spin two”. As an alternative to the linear polarization basis, the
polarization tensor $e_{ij}$ and its complex conjugate $e^{*}_{ij}$ form a
right- and left-circular polarization basis. (See also Problem 9 from Chapter
16 of [Hartle, 2002].)

Returning to (1), we first need to substitute $\tau(s)$ and $\vec{x}(s)$ for
$t$ and $\vec{x}$ in the plane wave expansion (2) for $h_{ij}(t,\vec{x})$. The
$s$ dependence shows up only in the exponential, so the integration is easy.
Again, we leave it as an exercise for the reader to show that doing the
integration and grouping terms accordingly gives

$\Delta T(t)=\int_{-\infty}^{\infty}{\rm d}f\>\int{\rm
d}\hat{k}\sum_{A=+,\times}h_{A}(f,\hat{k})R^{A}(f,\hat{k})e^{i2\pi
f(t-\hat{k}\cdot\vec{r}_{2}/c)}\,,$ (6)

where

$R^{A}(f,\hat{k})\equiv\frac{1}{2}u^{i}u^{j}e^{A}_{ij}(\hat{k})\frac{1}{i2\pi
f}\,\frac{1}{1-\hat{k}\cdot\hat{u}}\,\left[1-e^{-\frac{i2\pi
fL}{c}(1-\hat{k}\cdot\hat{u})}\right]\,.$ (7)

(For a solution, see equations (6.8)–(6.10) in [Romano, 2019] or equations
(C8)–(C9) in [Allen, 2023] or Appendix B and equations (C1)–(C2) in [Anholm et
al., 2009].) The quantity $R^{A}(f,\hat{k})$ is the detector response function
for a timing residual measurement made by a one-arm, one-way detector of
length $L$ responding to a plane-polarized GW having polarization $A$,
propagation direction $\hat{k}$, and frequency $f$.

Note that we have factored out an exponential term $e^{i2\pi
f(t-\hat{k}\cdot\vec{r}_{2}/c)}$ in the integrand of $\Delta T(t)$. For the
case of pulsar timing, the argument of this exponential is the phase of the GW
at the time $t\equiv t_{2}$ and location $\vec{r}_{2}$ of the reception of the
radio pulse at Earth. It is also easy to show that the product of $e^{i2\pi
f(t-\hat{k}\cdot\vec{r}_{2}/c)}$ and the second term in square brackets in (7)
is $e^{i2\pi f(t_{1}-\hat{k}\cdot\vec{r}_{1}/c)}$. The argument of this
exponential is the phase of the GW at the time $t_{1}$ and location
$\vec{r}_{1}$ of the emission of the radio pulse at the pulsar. Hence, in the
context of pulsar timing observations, the “1” in square brackets in (7)
corresponds to the “Earth term” while the exponential in the square brackets
corresponds to the “pulsar term”. The direction to the pulsar is
$\hat{p}=-\hat{u}$.

The integral leading to (7) is the difference between a function evaluated at
the upper limit (photon received) and lower limit (photon emitted). The
difference has a simple physical explanation first given by Detweiler
[Detweiler, 1979]. Rather than computing the timing residual, he instead gave
an (equivalent) analysis in terms of the effect of the GW on the photon’s
frequency. In this picture, the GW gives the photon source (pulsar) and photon
receiver (Earth) a nonzero velocity. It is the difference between these two
velocities, projected along the line of sight to the pulsar, which creates
(via the Doppler effect) a shift in the photon frequency. If both source and
receiver have the same projected velocity, then the Earth and pulsar term
cancel, and the GW has no net effect on the photon redshift/pulsar timing
residual.

The expression for the timing residual response function (7) may be written as
a product of two terms: (i) a frequency-independent (but polarization-
dependent) geometrical projection factor

$G^{A}(\hat{k})\equiv\frac{1}{2}u^{i}u^{j}\,e_{ij}^{A}(\hat{k})\,,$ (8)

and (ii) a frequency-dependent (but polarization-independent) timing transfer
function

${\cal T}(fL/c,\hat{k}\cdot\hat{u})\equiv\frac{1}{i2\pi
f}\,\frac{1}{1-\hat{k}\cdot\hat{u}}\,\left[1-e^{-\frac{i2\pi
fL}{c}(1-\hat{k}\cdot\hat{u})}\right]\,.$ (9)

The timing transfer function vanishes when the photon encounters an integer
number of complete cycles of GW-induced stretching and squeezing, which
integrate to zero net timing delay. At fixed GW frequency $f$, this occurs for
${\rm int}[2fL/c]$ different incident directions, given by values of
$\hat{k}\cdot\hat{u}$ that satisfy

$1-\hat{k}\cdot\hat{u}=\frac{n}{fL/c}\quad{\rm where}\quad n=1,2,\dots,{\rm
int}[2fL/c]\,.$ (10)

The limit on $n$ arises because the lhs of (10) lies in the interval $[0,2]$.
However, for a given GW direction (fixed $\hat{k}\cdot\hat{u}$) there are an
infinite number of frequencies at which $\cal T$ vanishes. For normal
incidence of the GW on the detector ($\hat{k}\cdot\hat{u}=0$, as in panel (b)
of Figure 2), these zeroes lie at $f=(c/L)n$ for $n=1,2,3,\dots$, as shown in
Figure 5.

Figure 5: Modulus of the timing transfer function ${\cal
T}(fL/c,\hat{k}\cdot\hat{u})$ for a normally incident
($\hat{k}\cdot\hat{u}=0$) GW. It vanishes when the arm length is a positive
integer number of GW wavelengths.

For a short-arm LIGO-like detector, the above expression for the timing
transfer function (9) and hence the response function (7) simplify
considerably. Expanding the exponential in (9) to first order in $fL/c\ll 1$
gives

${\cal T}(fL/c,\hat{k}\cdot\hat{u})\simeq\frac{L}{c}\quad\Rightarrow\quad
R^{A}(f,\hat{k})\simeq\frac{1}{2}\frac{L}{c}u^{i}u^{j}e_{ij}^{A}(\hat{k})\qquad(\text{short-
arm limit})\,.$ (11)

You can convert this to the more common strain response $h(t)\equiv\Delta
L(t)/L$ by simply dividing $\Delta T(t)$ and $R^{A}(f,\hat{k})$ by the nominal
one-way propagation time $T\equiv L/c$ of a photon down the arm of the
detector. Then ${\cal T}\simeq 1\Rightarrow R^{A}(f,k)\simeq G^{A}(\hat{k})$
for the short-arm strain response.

Note that this short-arm response function is symmetric under the interchange
$\hat{u}\rightarrow-\hat{u}$. This means that if the distance from the Earth
to a pulsar was much less than the GW wavelength, then the timing response of
that Earth-pulsar baseline $\hat{p}$ would be exactly the same as that for a
pulsar in the opposite sky direction, $-\hat{p}$. This can be seen graphically
in the leftmost panels of Figures 6 and 7.

For later use, we define another frequency-independent geometrical projection
factor

$F^{A}(\hat{k})\equiv\frac{1}{2}\frac{u^{i}u^{j}}{1-\hat{k}\cdot\hat{u}}\,e^{A}_{ij}(\hat{k})\,,$
(12)

which differs from $G^{A}(\hat{k})$ defined in (8) by a factor of
$1-\hat{k}\cdot\hat{u}$ in the denominator. The full timing response function
(7) may be written as

$R^{A}(f,\hat{k})=\frac{1}{i2\pi f}\,F^{A}(\hat{k})\left[1-e^{-\frac{i2\pi
fL}{c}(1-\hat{k}\cdot\hat{u})}\right]\,.$ (13)

The extra factor of $1-\hat{k}\cdot\hat{u}$ comes from the timing transfer
function (9), and breaks $G$’s $\hat{u}\rightarrow-\hat{u}$ symmetry. The
asymmetry manifests itself graphically in the middle and rightmost panels of
Figure 7. The response function $F^{A}(\hat{k})$ is often called the Earth-
term-only response function, since it doesn’t include the exponential term in
square brackets (i.e., the pulsar term) and is thus independent of frequency.

At first glance it appears that $F^{A}(\hat{k})$ blows up as
$\hat{k}\to\hat{u}$, since the denominator $1-\hat{k}\cdot\hat{u}$ in (12)
approaches zero. But this zero in the denominator is cancelled by a zero from
the geometrical projection factor $u^{i}u^{j}e^{A}_{ij}(\hat{k})$, so that
$F^{A}(\hat{k}=\hat{u})$ is finite. Actually, the full response function
$R^{A}(f,\hat{k}=\hat{u})$ vanishes, because the two terms in square brackets
in (13) cancel. This makes sense: a GW deforms space only in the plane
transverse to its direction of propagation $\hat{k}$. So when
$\hat{k}=\hat{u}$, the radio pulse does not feel these deformations. See
Figure 2, panel (a) for an example.

###  II.2 Visualizing the timing residual response (antenna pattern functions)

We now work out the timing residual response function of a one-arm, one-way
detector to a linearly-polarized GW source at a fixed point in the sky,
radiating at a single frequency, and uniformly averaged over linear
polarization directions. To compute this, start with a linearly-polarized GW
(containing only one of the two possible polarizations (3)) and average over
the choice of directions $\hat{l}$ and $\hat{m}$ that define what “+” and
“$\times$” means. This leads to the polarization-averaged response function

$R(f,\hat{k})\equiv\left[\frac{1}{2\pi}\int_{0}^{2\pi}{\rm
d}\psi\>R^{+}(f,\hat{k},\psi)R^{+*}(f,\hat{k},\psi)\right]^{1/2}\,.$ (14)

Here $R^{+}(f,\hat{k},\psi)$ is given by (7) with the standard polarization
tensors $e^{+}_{ij}(\hat{k})$ replaced by
$e^{\prime}{}^{+}_{ij}(\hat{k},\psi)$ given in (5). (If you are wondering:
defining (14) with only the $\times$ polarization leads to an identical
result.)

Since the $\psi$ dependence of the integrand in (14) involves quadratic
combinations of $\cos 2\psi$ and $\sin 2\psi$, the integral over $\psi$ is an
easy exercise for the reader to do, giving

$R(f,\hat{k})=\sqrt{\frac{1}{2}\Bigl{(}\bigl{|}R^{+}(f,\hat{k})\bigr{|}^{2}+\bigl{|}R^{\times}(f,\hat{k})\bigr{|}^{2}\Bigr{)}}=\frac{1}{\sqrt{2}}\,G(\hat{k})\bigl{|}{\cal
T}(fL/c,\hat{k}\cdot\hat{u})\bigr{|}\,,$ (15)

where

$G(\hat{k})\equiv\sqrt{G^{+}(\hat{k})^{2}+G^{\times}(\hat{k})^{2}}\,.$ (16)

Such quadrature combinations arises whenever we uniformly average over the
polarization direction of sources. For the case of a pulsar located in the
$\hat{p}=-\hat{u}=\hat{z}$ direction

$G(\hat{k})=\frac{1}{2}\sin^{2}\theta\,,\qquad
F(\hat{k})=\frac{1}{2}(1+\cos\theta)\,,$ (17)

where $\hat{k}\cdot\hat{u}=\cos\theta$ is the cosine of the angle of the GW
relative to the one-arm, one-way detector, and $F(\hat{k})$ is defined by a
quadrature combination similar to (16). You will need to use (3), (4), (8),
and (12) to obtain these explicit expressions. Note that, as discussed in the
final paragraph of Section II.1, $F$ is finite as $\theta\to 0$ and vanishes
as $\theta\to\pi$ (see the middle panel of Figure 7).

For fixed frequency $f$ and fixed distance $L$ to a pulsar, the shape of the
polarization-averaged response function $R(f,\hat{k})$ depends only on the
relative orientation of the GW propagation direction $\hat{k}$ and the
direction to the pulsar $\hat{p}=-\hat{u}$. The shape of $R(f,\hat{k})$ is
often called an antenna pattern, because the magnitude of $R$ is the
“strength” of the response. Figure 6 shows plots of the antenna patterns
$R(f,\hat{k})$ for a pulsar located in the $\hat{z}$ direction for different
fixed values of $fL/c$. The left panel shows the antenna pattern in the
“short-arm” LIGO-like limit, $fL/c\ll 1$, which is proportional to
$\sin^{2}\theta$, as shown in (11), (8), and (17). Note that the number of
lobes of the antenna pattern increases as $fL/c$ increases, and that the
overall magnitude of the antenna pattern function decreases like $1/(2\pi f)$.
Both of these behaviors are a consequence of the form of the timing transfer
function given in (9).

Figure 6:  Antenna pattern functions for the polarization-averaged timing
response of a one-arm, one-way detector having $\hat{u}=-\hat{z}$, evaluated
for different fixed values of $fL/c=L/\lambda$ (from left to right, $fL/c=0$,
0.5, 2, 10). The top row of plots are cross-sections: the antenna pattern
functions are axially symmetric about the vertical $z$-axis. For example, the
leftmost pattern resembles a donut with no hole, as shown in the bottom row of
plots, which are the corresponding antenna patterns viewed from slightly above
the $xy$ plane.  Figure 7: Antenna pattern functions for the polarization-
averaged geometrical projection factors $G(\hat{k})$ (left panel) and
$F(\hat{k})$ (middle panel). The right panel compares $F(\hat{k})$ (white
dashed curve) with the full polarization-averaged timing response function
$R(f,\hat{k})$ for $fL/c=100$ rescaled by a factor of $2\pi f$ so that it has
the same scale as $F(\hat{k})$. All of these antenna pattern functions are
axially symmetric around the vertical $z$-axis.

Recall that for most pulsars, $fL/c$ will be of order 100 or more, since the
distance to typical pulsars is of order a kpc (so a few thousand light-years)
or more, while the GW wavelength is only of order tens of light-years. The
rightmost panel of Figure 7 shows a rescaled version of $R(f,\hat{k})$
appropriate for the (semi-realistic) case $fL/c=100$. Also plotted in this
figure are the polarization-averaged antenna pattern functions for the
geometrical projection factors $G(\hat{k})$ (left panel) and $F(\hat{k})$
(middle panel). Note that $F(\hat{k})$ is half of the “envelope” of the full
response function $R(f,\hat{k})$, illustrated by the dashed-white curve in the
rightmost panel plot.

###  II.3 Calculating overlap functions and the Hellings and Downs correlation

Large amplitude GWs would be directly visible in GW detector output. For
example, in the early days of pulsar timing, limits were set on the amplitude
of GWs in this way, by examining the timing residuals from single pulsars. But
the GWs predicted for typical background sources are too weak to allow this:
they cannot be separated from noise in the measurements or other physical
effects. Fortunately, there is a way around this: look for a common GW signal
by correlating the outputs of two or more detectors. With sufficient
integration or averaging, this allows one to “dig below the noise”.

For this reason, we define the overlap function for a pair of detectors, which
is proportional to the expected correlation of their outputs. As we shall show
below, the Hellings and Downs correlation is an example of an overlap
function, calculated in the context of pulsar timing array observations,
evaluated for many different pairs of pulsars.

In our Universe, the GW background is (to very high precision) defined by (2)
for specific functions $h_{+}(f,\hat{k})$ and $h_{\times}(f,\hat{k})$, which
are the Fourier coefficients in the plane wave expansion of the metric
perturbations. Sadly, we do not know what these functions are. So, to
understand what one might see, and how to interpret the data, we are forced to
build imaginary mathematical “model universes”, and ask what might be observed
in those. Typically, these model universes assume that the sources have
particular distances, sky locations, and so on. But, since we don’t know these
values, we have to resort to a statistical approach. We construct many model
universes in which these unknowns take on different values, and then ask,
“what would we typically observe?” In short, our lack of knowledge about the
sources forces us to follow a statistical approach.

The set of model universes is called an “ensemble”. Our Universe itself is not
an ensemble, although, if we have done a very good job in constructing an
ensemble, our Universe might correspond very closely to one of the models in
that ensemble. The different models in the ensemble are called “realizations”;
each one corresponds to a different specific choice of the functions
$h_{A}(f,\hat{k})$. For example, several different ensembles are explicitly
constructed in Section III of [Allen, 2023].

To calculate the expected correlation of a pair of measurements containing
stochastic signals, we need to define the ensemble of possible realizations
associated with that stochastic process. Expectation values $\langle\ \rangle$
are averages over the realizations within the ensemble.

For GWB searches, it is common to consider a Gaussian ensemble. This is a good
choice, because the central limit theorem shows that the incoherent sum of
many weak processes is well described by such an ensemble.

Since the quantities of interest depend upon $h_{ij}(t,\vec{x})$, which in
turn are determined by the Fourier coefficients $h_{A}(f,\hat{k})$, it is
sufficient if we know how to calculate the expected value of products
containing any number of these, for example $\langle
h_{+}(f_{1},\hat{k}_{1})h_{+}(f_{2},\hat{k}_{2})h_{+}(f_{3},\hat{k}_{3})\rangle$.

A beautiful feature of the Gaussian ensemble was proved more than a century
ago [Isserlis, 1918]: all moments, of any order, are expressible in terms of
sums of products of the first two moments, $\langle
h_{A}(f_{1},\hat{k}_{1})\rangle$ and $\langle
h_{A}(f_{1},\hat{k}_{1})h_{A^{\prime}}(f_{2},\hat{k}_{2})\rangle$. We will
assume that these first two moments have the following simple form

$\langle h_{A}(f,\hat{k})\rangle=0\quad\text{and}\quad\langle
h_{A}(f,\hat{k})h^{*}_{A^{\prime}}(f^{\prime},\hat{k}^{\prime})\rangle=\frac{1}{8\pi}H(f)\delta(f-f^{\prime})\delta_{AA^{\prime}}\delta^{2}(\hat{k},\hat{k}^{\prime})\,,$
(18)

where the non-negative real functiom $H(f)=H(-f)$ is the (two-sided) strain
power spectral density of the background, having units of ${\rm
strain}^{2}\,{\rm Hz}^{-1}$. The first moment is zero since for every
realization in the ensemble there is another identical one with its phase
shifted by an arbitrary amount. The $\delta(f-f^{\prime})$ on the rhs of the
second moment arises because the ensemble is stationarity: there is no
preferred origin of time because for every realization in the ensemble, there
is another identical one shifted forwards or backwards in time. The
$\delta_{AA^{\prime}}$ implies that the ensemble is unpolarized: for every
realization, there is another one that is identical except that the
$\hat{l},\hat{m}$ vectors are rotated through an arbitrary angle $\psi$.
Finally, the factor of $\delta^{2}(\hat{k},\hat{k}^{\prime})$ arises because
the ensemble is spatially homogeneous and isotropic: for every realization in
the ensemble and for every point in space, there is another realization that
is related to the first by a rotation about that point by an arbitrary set of
Euler angles.

We are now ready to calculate the expected value of two timing residual
measurements $\Delta T_{a}(t)$ and $\Delta T_{b}(t^{\prime})$, where the two
times $t$ and $t^{\prime}$ may be different. The indices $a$ and $b$ label the
two pulsars; $\hat{p}_{a}$ and $\hat{p}_{b}$ will denote their directions on
the sky, and $L_{a}$ and $L_{b}$ will denote their distances from Earth. We
leave it as another exercise for the reader to show that by using (18) and
(6),

$\displaystyle\langle\Delta T_{a}(t)\Delta T_{b}(t^{\prime})\rangle$
$\displaystyle=\int_{-\infty}^{\infty}{\rm d}f\>e^{i2\pi
f(t-t^{\prime})}\Gamma_{ab}(f)H(f)\,,$ (19)

where

$\Gamma_{ab}(f)\equiv\frac{1}{8\pi}\int{\rm
d}\hat{k}\sum_{A}R_{a}^{A}(f,\hat{k})R_{b}^{A*}(f,\hat{k})e^{-i2\pi
f\hat{k}\cdot(\vec{r}_{a}-\vec{r}_{b})/c}$ (20)

is the overlap function for the two detectors. The concept of an overlap
function is a general construction; it is not restricted to timing residual
measurements for pulsar timing arrays, see e.g., [Flanagan, 1993; Christensen,
1997; Finn et al., 2009]. From (19), one sees that $\Gamma_{ab}(f)$ plays the
role of a transfer function relating the expected correlation in the output of
two detectors to the auto-correlated power.

To proceed further, we use (13) to expand the timing residual response
functions in (20). This yields

$\displaystyle\Gamma_{ab}(f)$ $\displaystyle=\frac{1}{(2\pi
f)^{2}}\frac{1}{8\pi}\int{\rm
d}\hat{k}\sum_{A}F_{a}^{A}(\hat{k})F_{b}^{A}(\hat{k})e^{-i2\pi
f\hat{k}\cdot\cancel{(\vec{r}_{a}-\vec{r}_{b})}/c}$ (21) $\displaystyle\hskip
108.405pt\times\left[1-\cancel{e^{-\frac{i2\pi
fL_{a}}{c}(1+\hat{k}\cdot\hat{p}_{a})}}\right]\left[1-\cancel{e^{+\frac{i2\pi
fL_{b}}{c}(1+\hat{k}\cdot\hat{p}_{b})}}\right]\,,$

where $\hat{p}_{a}=-\hat{u}_{a}$ and $\hat{p}_{b}=-\hat{u}_{b}$, and where we
have crossed off several terms for the following reasons: (i) The argument of
the exponential involving the detector locations $\vec{r}_{a}$ and
$\vec{r}_{b}$ is approximately zero since the separation between two radio
telescopes on Earth is much smaller than the wavelength $c/f$ of the GWs that
pulsar timing arrays are sensitive to. (In fact the observations are generally
referred back to the solar system barycenter (SSB) in which case
$\vec{r}_{a}=\vec{r}_{b}$.) (ii) For distinct pulsars ($a\neq b$), we can
ignore the pulsar terms in the square brackets, since the two Earth-pulsar
product terms and the single pulsar-pulsar product term are rapidly varying
function of $\hat{k}$ and do not contribute significantly when integrated over
the whole sky, see [Hellings, 1989] and [Anholm et al., 2009]. Howevever, for
the auto-correlation of a single pulsar with itself (i.e., $a=b$), the product
of the two pulsar terms equals 1 and hence cannot be ignored. Thus,

$\Gamma_{ab}(f)=\frac{1}{(2\pi f)^{2}}\frac{1}{8\pi}\int{\rm
d}\hat{k}\sum_{A}F_{a}^{A}(\hat{k})F_{b}^{A}(\hat{k})\left[1+\delta_{ab}\right]=\frac{1}{12\pi^{2}f^{2}}\,\Gamma_{ab}\,,$
(22)

where

$\Gamma_{ab}\equiv\frac{3}{8\pi}\int{\rm
d}\hat{k}\sum_{A}F_{a}^{A}(\hat{k})F_{b}^{A}(\hat{k})\left[1+\delta_{ab}\right]\,.$
(23)

The factor of $3/8\pi$ before the integral makes the overlap function approach
$1/2$ for distinct pulsars separated by a small angle on the sky.

As discussed in FAQ Q1, this normalization choice is arbitrary: the physical
quantity is $\langle\Delta T_{a}(t)\Delta T_{b}(t^{\prime})\rangle$ on the lhs
of (19). We are free to rescale the definitions of $\Gamma_{ab}(f)$ and $H(f)$
(given in (22) and (18)) in any way that we want provided that their product
(and hence the physical quantity in (19)) is unchanged.

Note that the final expression for $\Gamma_{ab}$ is independent of frequency.
We have factored out a frequency-dependent term that is independent of the
specifics of the particular pulsars being correlated. This is a consequence of
being able to drop the three exponential terms in (21). (See FAQ Q5 for a
related discussion.)

Simple symmetry arguments show that the integral in (23) can only depend upon
the angle $\gamma_{ab}$ between the lines of sight to the two pulsars. Without
loss of generality, place pulsar $a$ on the $z$-axis and pulsar $b$ in the
$xz$-plane so that $\hat{p}_{a}=\hat{z}$ and
$\hat{p}_{b}=\sin\gamma_{ab}\,\hat{x}+\cos\gamma_{ab}\,\hat{z}$. The
polarisation tensors needed to calculate $F^{A}_{a}(\hat{k})$ and
$F^{A}_{b}(\hat{k})$ are defined by (3), with $\hat{k}$, $\hat{l}$, $\hat{m}$
given in (4). The integral can be approximated numerically or evaluated
analytically using contour integration, see e.g., [Hellings, 1989; Anholm et
al., 2009; Jenet and Romano, 2015; Allen, 2023]:

$\Gamma_{ab}\equiv\frac{1}{2}+\frac{3}{2}\left(\frac{1-\cos\gamma_{ab}}{2}\right)\left[\ln\left(\frac{1-\cos\gamma_{ab}}{2}\right)-\frac{1}{6}\right]+\frac{1}{2}\,\delta_{ab}\,,$
(24)

where $\cos\gamma_{ab}=\hat{p}_{a}\cdot\hat{p}_{b}$. This is an analytic
expression for the Hellings and Downs correlation for pulsar timing, expressed
as a function of the angular separation between the directions to a pair of
pulsars. A plot of this expected correlation is shown in the left panel of
Figure 8, and also as the black dashed curve in Figure 1.

Figure 8: Left panel: Hellings and Downs correlation for pairs of PTA pulsars.
Right panel: The Hellings and Downs correlation for pairs of short-arm LIGO-
like detectors. Both have been normalized to 1/2 at zero angular separation.
Note that the zero angular separation in the left hand panel corresponds to a
pair of pulsars located close on the sky but at significantly different
distances. If they are much closer than a GW wavelength, or if they are the
same pulsar, then the correlation doubles to unity.

We leave it as an exercise for the reader to show that a similar calculation
of the overlap function for an array of “short-arm” LIGO-like detectors yields

$\Gamma_{ab}=\frac{1}{2}P_{2}(\cos\gamma_{ab})\,,$ (25)

where $P_{2}(x)$ is the $l=2$ Legendre polynomial. For this calculation, you
will need to replace the response functions $R^{A}(f,\hat{k})$ in (20) with
the short-arm limiting expressions (11), and assume that the exponential term
involving the detector locations can again be ignored. (But see FAQ Q5
regarding the calculation of overlap functions for real ground-based
interferometers.) This substitution simplifies the integrand in (23) to
quadratics in $\sin\phi$ and $\cos\phi$ and quartics in $\cos\theta$. The
short-arm overlap function is plotted in the right panel of Figure 8. Since
the short-arm overlap function only involves the $l=2$ Legendre polynomial, it
is an example of a purely quadrupolar overlap function.

Finally, to return to the original question about why the Hellings and Down
correlation has different values at $0^{\circ}$ and $180^{\circ}$, we plot the
integrand
$F^{+}_{1}(\hat{k})F^{+}_{2}(\hat{k})+F^{\times}_{1}(\hat{k})F^{\times}_{2}(\hat{k})$
of the Hellings and Downs correlation for these two different cases in the
left panel of Figure 9. It is apparent from this plot that if the GWB
consisted of waves that only came from directions perpendicular to the Earth-
pulsar baselines (so in the $xy$-plane of this plot), then the values of the
Hellings and Downs correlation at $0^{\circ}$ and $180^{\circ}$ would be
exactly the same. But since an isotropic GWB has equal contributions from GWs
coming from all different directions on the sky, the Hellings and Downs
correlation for the case $\hat{p}_{1}=\hat{p}_{2}$ will be larger than that
for $\hat{p}_{1}=-\hat{p}_{2}$. Why these values differ by exactly a factor of
two is the topic of the next question.

Figure 9: Integrands of the Hellings and Downs correlation for two pulsars
separated by $0^{\circ}$ ($\hat{p}_{1}=\hat{p}_{2}=\hat{z}$, solid curve) and
by $180^{\circ}$ ($\hat{p}_{1}=-\hat{p}_{2}=\hat{z}$, dashed curve). The first
two panels are polar plots, which are axially symmetric about the $z$-axis.
The left plot shows
$F_{1}^{+}(\hat{k})F_{2}^{+}(\hat{k})+F_{1}^{\times}(\hat{k})F_{2}^{\times}(\hat{k})$.
The right plot includes an additional factor of $\sin\theta$, which arises
when integrating over $\theta$. The right panel is the same as the middle
plot, but shown as a function of $\cos\theta$. The area under the solid curve
is twice the area under the dashed curve.

Q3: Is there a simple way to see that the value of the Hellings and Downs
correlation for two pulsars separated by $180^{\circ}$ is half that for
$0^{\circ}$ angular separation?

A3: Yes, but to prove this we need to do a few calculations.

We start with (23) for the Hellings and Downs correlation and (12) for the
Earth-term-only response function. Taking the two pulsars to point in the same
direction ($\hat{p}_{1}=\hat{p}_{2}=\hat{z}$), we have

$\sum_{A}F_{1}^{A}(\hat{k})F_{2}^{A}(\hat{k})=\frac{1}{4}(1+\cos\theta)^{2}\,.$
(26)

Having them point in opposite directions ($\hat{p}_{1}=-\hat{p}_{2}=\hat{z}$)
leads to

$\sum_{A}F_{1}^{A}(\hat{k})F_{2}^{A}(\hat{k})=\frac{1}{4}(1+\cos\theta)(1-\cos\theta)=\frac{1}{4}\sin^{2}\theta\,.$
(27)

These functions are plotted in Figure 9, where $\theta$ is the usual polar
angle measured with respect to the $z$-axis. Multiplying by $3/8\pi$ and
integrating over the sphere $\int_{0}^{2\pi}d\phi\int_{0}^{\pi}{\rm
d}\theta\,\sin\theta$ we find

$\Gamma_{12}=\frac{1}{2}\quad({\rm for}\
\hat{p}_{1}=\hat{p}_{2}=\hat{z})\,,\qquad\Gamma_{12}=\frac{1}{4}\quad({\rm
for}\ \hat{p}_{1}=-\hat{p}_{2}=\hat{z})\,.$ (28)

From the three panel plots in Figure 9, we see that when the two pulsars both
point in the $\hat{z}$ direction, the majority of support for the overlap
function comes from sky directions $\hat{n}=-\hat{k}$ having $z>0$. When the
two pulsars point in opposite directions, $\hat{p}_{1}=-\hat{p}_{2}=\hat{z}$,
the majority of support for the overlap function comes from sky directions
having $z<0$, which is a smaller contribution.

Visual inspection of the areas under the two curves shown in the rightmost
panel of Figure 9 shows why the value of the HD curve for two pulsars
separated by $0^{\circ}$ is twice as large as that for two pulsars separated
by $180^{\circ}$.

Q4: Why isn’t the most-negative value of the Hellings and Downs correlation at
$90^{\circ}$?

A4: Again, this is because the Earth-term-only response function is not
symmetric under $\hat{u}\rightarrow-\hat{u}$. If the Hellings and Downs
correlation were purely quadrupolar then its most-negative value would be
exactly at $90^{\circ}$, as it is for the case of an array of short-arm LIGO-
like detectors (see the right panel of Figure 8). If it had only even
multipoles (no odd-$l$ terms) then it would be reflection symmetric about
$90^{\circ}$ but might not have its minimum there.

Although the dominant contribution to the Hellings and Downs correlation is
the quadrupole pattern $P_{2}(\cos\gamma)$, there are non-zero contributions
from all higher order $l$-modes as well. We leave it as an exercise for the
reader to show that

$\Gamma(\gamma)=\frac{3}{2}\sum_{l=2}^{\infty}(2\ell+1)c_{l}P_{l}(\cos\gamma)\,,\quad\text{where}\quad
c_{l}=\frac{1}{(l+2)(l+1)l(l-1)}\,.$ (29)

To obtain this result, you need to use (a) the orthogonality property

$\int_{-1}^{1}{\rm
d}x\>P_{l}(x)P_{l^{\prime}}(x)=\frac{2}{2l+1}\delta_{ll^{\prime}}$ (30)

of the Legendre polynomials, where $x\equiv\cos\gamma$, (b) Rodrigues’ Formula
for Legendre polynomials, and (c) integration by parts. See equations
(97)–(105) in [Gair et al., 2014] for a step-by-step derivation of (29).

Figure 10 compares approximate versions of the Hellings and Downs correlation,
terminating the sum at different maximum values $l_{\rm max}=2,3,4,5$. Since
the $c_{l}$’s in (29) fall off like $1/l^{4}$, only a few values of $l$ are
needed to get a very good approximation.

Figure 10: Approximations to the Hellings and Downs curve (solid black)
obtained by terminating the sum in (29) with two to five Legendre polynomials
(grey curves).

Q5: Why is the Hellings and Downs correlation frequency independent, whereas
overlap functions for ground-based interferometers are frequency dependent?

A5: Actually, the Hellings and Downs correlation for timing residual
measurements is frequency dependent. But this $1/(2\pi f)^{2}$ frequency
dependence is the same for all pairs of pulsars and hence can be “factored
out” of the Hellings and Downs correlation function (22) and absorbed into the
strain spectrum $H(f)$ in (19). (This changes the dimensions of $H(f)$ from
strain${}^{2}\,{\rm Hz}^{-1}$ to s${}^{2}\,{\rm Hz}^{-1}$.) This factorization
is possible because the exponential terms in (21) that depend on the distances
to the pairs of pulsars defining the correlation and the locations of the
radio telescopes making the measurements can be ignored when calculating the
expected correlation, for the reasons discussed in the sentences following
(21).

For ground-based interferometers, the frequency dependence of the correlation
comes from the exponential factor $e^{-i2\pi
f\hat{k}\cdot(\vec{r}_{a}-\vec{r}_{b})/c}$ in (20), since $R^{A}(f,k)$ is
independent of frequency in the short-arm limit, see (11). Here, $\vec{r}_{a}$
and $\vec{r}_{b}$ denote the position vectors of the two interferometers’s
central stations. For two such detectors on Earth, the distance
$D_{ab}\equiv|\vec{r}_{a}-\vec{r}_{b}|$ between the stations is comparable to
the wavelength $\lambda\equiv c/f$ (hundreds to thousands of km) of the GWs
that these detectors are sensitive to. Hence, it cannot be ignored as we were
able to do for the long-arm PTA pulsar calculation in (21). So even though the
individual ground-based interferometers operate in the short-arm limit
$L\ll\lambda$, it is not also the case that $D_{ab}\ll\lambda$, but rather
$\lambda\lesssim D_{ab}$. Since the frequency dependence for ground-based
interferometers is tied to the detector locations, it cannot be factored out
of the correlation $\Gamma_{ab}(f)$ and absorbed into $H(f)$.

Q6: Is the Hellings and Downs correlation only relevant for isotropic GWBs?

A6: No. The Hellings and Downs pattern is a robust prediction for the
correlation produced by any (small or large) set of non-interfering GW
sources, provided that many pairs of pulsars, uniformly distributed on the
sky, are employed. This was first demonstrated for the special case of GWs
produced by a single black hole binary in [Cornish and Sesana, 2013]. It was
then extended to an arbitrary set of non-interferring sources in [Allen,
2023], which further exploited pulsar averaging. Note that for interferring
sources, we do not expect to recover exactly the HD pattern (see FAQ Q7 and
[Allen, 2023]).

The original (and traditional) approach to obtain the Hellings and Down
correlation is to consider just one pair of pulsars, and to average their
correlation over GW propagation directions, assuming that the sources are
unpolarized and isotropic on the sky [Hellings, 1989]. In references [Cornish
and Sesana, 2013] and [Allen, 2023], somewhat involved calculations were done
to show the equivalence of the two approaches.

But, in fact, one can show the equivalence with no calculation. The geometry
of the pulsar averaging is equivalent to that of the isotropic GW-source
averaging after a series of Euler-angle rotations, provided for the GW
averaging calculation one averages over the polarization angle $\psi$ in
addition to the GW propagation direction $\hat{k}=\hat{k}(\theta,\phi)$. Note
that the original [Hellings, 1989] calculation obscures this point by not
giving enough detail. To reproduce the result, one must also assign a
polarization to each source, and average over that. Hence, averages must be
calculated with respect to three angles, not just two.

To see the equivalence purely from geometry, start with the GW averaging
calculation. As often done, without loss of generality, put $\hat{p}_{1}$
along the $z$-axis and $\hat{p}_{2}$ in the $xz$-plane, so that
$\hat{p}_{2}=\sin\gamma\,\hat{x}+\cos\gamma\,\hat{z}$. Pulsar averaging is
carried out with respect to the angular coordinates of $\hat{p}_{1}$ (polar
coordinate $\theta$, azimuthal coordinate $\psi$) and those of $\hat{p}_{2}$
(which lies on a cone about $\hat{p}_{1}$ with fixed half-opening angle
$\gamma$, with its location on the cone parametrized by $\phi\in[0,2\pi)$).
See Figure 11 and compare with Figure 5 of [Allen, 2023].

Figure 11:  Isotropic GW-source averaging (left panel) and pulsar averaging
(right panel) are equivalent, as can be seen by carrying out a series of
rotations.

The sequence of rotations that is needed to go from GW-source averaging to
pulsar averaging is to first rotate by $-\phi$ about the $z$ axis, which moves
the GW propagation direction $\hat{k}$ into the $xz$-plane and $\hat{p}_{2}$
around a cone with central axis $\hat{p}_{1}$ through angle $\phi$. One then
rotates by $-\theta$ around the $y$ axis, which puts $\hat{k}$ along
$-\hat{z}$, and makes $\hat{p}_{1}$ have polar angle $\theta$. Finally, one
rotates around $\hat{k}$ (which now points in the $-\hat{z}$ direction) by
$-\psi$, which becomes the azimuthal angle for $\hat{p}_{1}$ (with respect to
the $-x$ axis).

The fact that one obtains the Hellings and Downs correlation from pulsar
averaging means that the Hellings and Downs correlation is not restricted to
isotropic GWBs. But to exactly recover Hellings and Downs, one must assume
that there are no interference cross-terms that enter the correlation (see the
discussion following Eq. 4.4 of [Allen, 2023] and the discussion in [Allen,
2022]).

In practice, we do not have an infinite number of puslar pairs separated by
the angle $\gamma$. So, to do pulsar averaging, one divides the $N_{\rm
pul}(N_{\rm pul}-1)/2$ inter-pulsar correlations into a smaller number $N_{\rm
bins}$ of angular separation bins. This form of pulsar averaging is always
used to recover the Hellings and Downs curve. For example, for the NANOGrav
15-year analysis shown in Figure 1, $N_{\rm pul}=67$ and $N_{\rm bins}=15$,
which implies that the estimated correlation value in each angular separation
bin involved an average of approximately 150 pulsar pairs.

Q7: In the future, PTA observations will be carried out with many more
pulsars. We also expect that the GWB signal contribution will grow much faster
than pulsar noise as the observation time increases. Does that mean that at
some time in the distant future, we’ll exactly recover the Hellings and Downs
correlation?

A7: This would be the case if the GWB were produced by a set of discrete non-
interferring sources. (See for example Section IV and Equation 4.1 of [Allen,
2023] and the discussion in [Allen, 2022].)

But it is much more likely that our Universe contains many GW sources
radiating at freqencies that are close enough together that they cannot be
resolved individually over the time-span of the PTA experiments (say, tens to
hundreds of years). This means that the GWB is produced by a sum of
interfering sources. An example of such a GWB is the “confusion noise” created
by a uniform distribution of sources emitting monochromatic GWs with random
phases, but occupying the same frequency bin. In this case, pulsar-averaged
noise-free measurements over an infinite number of pulsars pairs with the same
angular separation do not give exactly the Hellings and Downs curve.

If we have an infinite collection of universes, each of which has different
random phases, and we were also able to average over all of those universes,
then we would get exactly the Hellings and Downs curve. But in any given
universe, we would find something different from the Hellings and Downs curve.
The (mean squared) difference between what would be observed in any given
universe, and the average over many universes, is called the “cosmic
variance”.

To compute this cosmic variance $\sigma^{2}_{\rm cosmic}(\gamma)$ for a given
ensemble of model universes, we express the pulsar-averaged correlation
$\Gamma_{i}$ in the $i$th universe as a deviation $\Delta\Gamma_{i}$ away from
the many-universe average $\Gamma$:

$\Gamma_{i}(\gamma)=\Gamma(\gamma)+\Delta\Gamma_{i}(\gamma)\quad(\text{for
universe }i).$ (31)

Then

$\langle\Gamma_{i}(\gamma)\rangle=\Gamma(\gamma)\quad{\rm
but}\quad\left\langle\bigl{(}\Gamma_{i}(\gamma)-\Gamma(\gamma)\bigr{)}^{2}\right\rangle=\langle\Delta\Gamma_{i}^{2}(\gamma)\rangle\equiv\sigma^{2}_{\rm
cosmic}(\gamma)\,,$ (32)

where the ensemble-average angle brackets denote an average over $i$. See
[Allen, 2023, 2022; Allen and Romano, 2023] for explicit calculations of
$\sigma^{2}_{\rm cosmic}(\gamma)$ for several simple but realistic models of
the Universe.

Q8: Assuming noise-free observations, are the fluctuations away from the
expected Hellings and Down curve due solely to the pulsar-term contributions
to the timing residual response?

A8: No. While the pulsar terms do contribute to deviations away from the
expected Hellings and Downs correlation, there are other more fundamental
sources of such deviations. (i) Timing residual measurements from two
different pulsar pairs separated by the same angle $\gamma$ but pointing in
different directions on the sky will have different correlations due to the
different parts of the GWB that they are sensitive to. So if the pulsar pairs
are not uniformly distributed, this gives deviations. (ii) For interfering
sources, even with infinite pulsar averaging as described in FAQ Q7, there
will be fluctuations away from expected Hellings and Downs correlation that
come from the Earth-term-only response functions.

For a given pulsar pair, the pulsar-term contributions to the timing residual
response add to the fluctuations, with the same root-mean-square amplitude as
the Earth terms. But, as mentioned in the previous paragraph, they are not the
sole source of the fluctuations away from the Hellings and Downs correlation,
and as more and more pairs of pulsars are averaged together at a particular
angular separation, this pulsar-term contribution becomes insignificant.

Q9: Why does the Hellings and Downs recovery plot shown in Figure 1) (from the
recent NANOGrav publication [Agazie et al., 2023]) have error bars that are
much larger than one would expect for the reported 3 to $4\sigma$ detection
significance?

A9: While the match at any given angle is not very close, the detection
significance comes from comparing two “waveforms” to compute a single “signal-
to-noise ratio”, in the same way that LIGO compares snippets of its noisy
time-domain data to noise-free chirp templates. Effectively, this signal-to-
noise ratio sums products of the observed pulsar-pair correlations with the
“expected waveform” given by the Hellings and Downs curve. The individual
terms of this sum are inversely weighted by estimates of their expected
squared values in the absence of the Hellings and Downs correlations. So, even
though the match at any given angle is not very close, the overall agreement
is actually quite good. In fact, NANOGrav obtains significance estimates from
two types of analyses:

(i) A 3.8$\sigma$ significance determined using frequentist detection
statistics, such as the noise-weighted inner product of the observed
correlations with the expected Hellings and Downs correlation, as described
above. This is the “stochastic equivalent” of the deterministic matched filter
statistic used by LIGO to search for individual compact binary coalescence
“chirp” signals. The distribution that would be expected in the absence of a
GWB is estimated by artificially introducing phase shifts which remove the
Hellings and Downs correlations (see the right panel of Figure 3 in Ref.
[Agazie et al., 2023]). The observed value of the detection statistic lies in
the “tail” of the null distribution, and thus corresponds to high detection
significance.

(ii) A 3.0$\sigma$ significance obtained via a Bayesian analysis. This
analysis compares the evidence for a data model containing the Hellings and
Downs spatial correlations to a model without these correlations. The ratio of
the evidences is 226, which is the “betting odds” for which of the two models
is correct. The significance of the value 226 is then estimated from the null
distribution of the evidence ratio, again constructed using phase shifts (see
both Figure 2 and the left panel of Figure 3 in Ref. [Agazie et al., 2023]).

What is shown in Figure 1 for the NANOGrav 15-year analysis [Agazie et al.,
2023] is instead a “self-consistency check” of the pulsar-pair cross-
correlations, calculated under the hypothesis that a GWB with Hellings and
Downs correlations is present in the data, with an overall amplitude estimated
from auto-correlation measurements, and assuming that there are no unknown
sources of noise with similar spatial correlations. The estimated correlations
and error bars are optimally weighted combinations of approximately 150
correlation measurements in each angular separation bin, which include the
covariances between the correlations induced by the GWB. If the observed data
are consistent with the Hellings and Downs correlation model then these point
estimates and 1$\sigma$ error bars will cross the expected Hellings and Downs
curve approximately 70% of the time on average. See [Allen and Romano, 2023]
for more details.

##  III Conclusion

We have answered a number of frequently asked questions (FAQs) about the
Hellings and Downs correlation for PTA searches. Many of these questions arise
from inadvertently applying intuition about the effects of GWs on a detector,
which are appropriate for LIGO-Virgo observations, to the case of PTA
observations, where some of that intuition is no longer valid. The key
difference is that a PTA detector has “arm lengths” (Earth-pulsar distances)
much longer than the wavelengths of the GWs that they are sensitive to, while
ground-based laser interferometers (like LIGO, Virgo, KAGRA) operate in the
opposite (i.e., “short-arm”) limit.

By calculating the response of a one-arm, one-way detector to a passing GW,
and then comparing the limiting expressions of this response for “short-arm”
LIGO-like detectors and “long-arm” PTA pulsar detectors, we were able to
obtain qualitative and quantitative answers to these questions. We hope that
this “Hellings-and-Downs correlation FAQ sheet” will help to inform the
scientific community and educate the general public about the current and
ongoing PTA searches for GWs.

###### Acknowledgements.

This paper and its “FAQ sheet” format arose from discussions that JDR had with
his wife and with several students attending the “MaNiTou Summer School on
GWs,” held from 3 July to 8 July 2023, in Nice, France. JDR acknowledges past
conversations with numerous NANOGrav colleagues about various aspects of these
topics, and financial support from the NSF Physics Frontier Center Award
PFC-2020265. He also thanks colleagues Nelson Christensen and Marie-Anne
Bizouard for their support and hospitality while visiting the Artemis Lab,
Observatoire Cote D’Azur in June and July 2023, where the majority of this
paper was written.

## References

  * Reardon et al. (2023) Daniel J. Reardon, Andrew Zic, Ryan M. Shannon, et al., “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,” ApJ 951, L6 (2023).
  * Xu et al. (2023) Heng Xu, Siyuan Chen, Yanjun Guo, et al., “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,” Research in Astronomy and Astrophysics 23, 075024 (2023).
  * Antoniadis et al. (2023) J. Antoniadis, P. Arumugam, S. Arumugam, et al., “The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals,” arXiv e-prints , arXiv:2306.16214 (2023).
  * Agazie et al. (2023) Gabriella Agazie, Akash Anumarlapudi, Anne M. Archibald, et al., “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,” ApJ 951, L8 (2023).
  * Hellings (1989) Ronald W. Hellings, “Pulsar timing and gravitational waves,” in Relativistic Gravitational Experiments in Space: proceedings of NASA Workshop, Annapolis, Maryland, June 28-30, 1988, NASA Conference Publication, Vol. 3046 (1989) pp. 93–97.
  * Allen (2023) B. Allen, “Variance of the Hellings-Downs correlation,” Phys. Rev. D 107, 043018 (2023).
  * 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).
  * Misner et al. (1973) C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation, Misner:1973 (W.H. Freeman and Company, San Francisco, 1973).
  * Hartle (2002) James B. Hartle, Gravity: An Introduction to Einstein’s General Relativity (Addison Wesley, 2002).
  * Romano (2019) Joseph D. Romano, “Searches for stochastic gravitational-wave backgrounds,” arXiv e-prints , arXiv:1909.00269 (2019).
  * Detweiler (1979) S. Detweiler, “Pulsar timing measurements and the search for gravitational waves,” ApJ 234, 1100–1104 (1979).
  * Isserlis (1918) L. Isserlis, “On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables,” Biometrika 12, 134–139 (1918), note that Isserlis’ theorem is often called “Wick’s theorem” in the physics community, although the Wick’s work was three decades later.
  * Flanagan (1993) Éanna É. Flanagan, “Sensitivity of the Laser Interferometer Gravitational Wave Observatory to a stochastic background, and its dependence on the detector orientations,” Phys. Rev. D 48, 2389 (1993).
  * Christensen (1997) N. Christensen, “Optimal detection strategies for measuring the stochastic gravitational radiation background with laser interferometric antennas,” Phys. Rev. D 55, 448–454 (1997).
  * Finn et al. (2009) Lee Samuel Finn, Shane L. Larson, and Joseph D. Romano, “Detecting a stochastic gravitational-wave background: The overlap reduction function,” Phys. Rev. D 79, 062003 (2009).
  * Jenet and Romano (2015) Fredrick A. Jenet and Joseph D. Romano, “Understanding the gravitational-wave Hellings and Downs curve for pulsar timing arrays in terms of sound and electromagnetic waves,” American Journal of Physics 83, 635–645 (2015).
  * Gair et al. (2014) Jonathan Gair, Joseph D. Romano, Stephen Taylor, and Chiara M. F. Mingarelli, “Mapping gravitational-wave backgrounds using methods from CMB analysis: Application to pulsar timing arrays,” Phys. Rev. D 90, 082001 (2014).
  * Cornish and Sesana (2013) Neil J Cornish and A Sesana, “Pulsar timing array analysis for black hole backgrounds,” Classical and Quantum Gravity 30, 224005 (2013).
  * Allen (2022) Bruce Allen, “Will pulsar timing arrays observe the Hellings and Downs correlation curve?” in 18th Vulcano Workshop: Frontier Objects in Astrophysics and Particle Physics, Vol. 74, edited by Antonella Antonelli, Roberto Fusco Femiano, Aldo Morselli, and Gian Carlo Trinchero (2022) pp. 65–80.
  * Allen and Romano (2023) Bruce Allen and Joseph D. Romano, “The Hellings and Downs correlation of an arbitrary set of pulsars,” (2023).