# Pulsar Timing Array source ensembles

Bruce Allen<EMAIL_ADDRESS>Max Planck Institute for Gravitational
Physics (Albert Einstein Institute), Leibniz Universität Hannover,
Callinstrasse 38, D-30167, Hannover, Germany  Serena Valtolina
<EMAIL_ADDRESS>Max Planck Institute for Gravitational Physics
(Albert Einstein Institute), Leibniz Universität Hannover, Callinstrasse 38,
D-30167, Hannover, Germany

###### Abstract

The stochastic gravitational wave background for pulsar timing arrays is often
modeled by a Gaussian ensemble which is isotropic and unpolarized. However,
the Universe has a discrete set of polarized gravitational wave sources at
specific sky locations. Can we trust that the Gaussian ensemble is an accurate
description? To investigate this, we explicitly construct an ensemble
containing $N$ individual binary sources with circular orbits. The orbital
inclination angles are randomly distributed, hence the individual sources are
elliptically polarized. We then compute the first two moments of the Hellings
and Downs correlation, as well as the pulsar-averaged correlation mean and
(cosmic) variance. The first moments are the same as for a previously studied
ensemble of circularly polarized sources. However, the second moments, and
hence the variances, are different for the two ensembles. While neither
discrete source model is exactly described by a Gaussian ensemble, we show
that in the limit of large $N$, the differences are small.

##  I Introduction

Pulsar timing arrays (PTAs) appear to be on the verge of making the first
detections [1, 2, 3, 4] of nano-Hertz gravitational waves (GWs). Much of the
data analysis and tools [5] assume that the GWs are described by an isotropic
and unpolarized Gaussian ensemble [6, 7, 8]. However, this is an
approximation: a realistic collection of discrete sources is only described by
the Gaussian ensemble in the limit of an infinite density of weak sources [9].
This paper investigates the differences.

The most likely PTA GW sources [10] are gravitationally bound binary systems
consisting of two black holes with masses $10^{4}-10^{10}M_{\odot}$, where
$M_{\odot}$ is the solar mass. In the simplest source models, which we adopt
here, the orbits are assumed to be circular.

The polarization of the GWs produced by a circular binary depend upon the
angle of inclination $\iota$. This is the angle between the line of sight and
the orbital angular momentum. The GWs which are produced are linearly
polarized if $\iota=90^{\circ}$ (orbit seen edge-on). The GWs are circularly
polarized if $\iota=0^{\circ}$ or $\iota=180^{\circ}$ (orbit seen face-on or
face-off). For other values of $\iota$, the GWs are elliptically polarized.
This begs the question, would an ensemble of such elliptically polarized
sources produce a GW background with the same properties as an ensemble of
unpolarized sources? (Note that the language can be misleading. A source is
called “unpolarized” if (a) the time average of $(h^{+})^{2}-(h^{\times})^{2}$
vanishes and (b) the time average of $h^{+}h^{\times}$ vanishes; see the
paragraph following Eq. 24 of [11]. Thus, a circularly polarized source is
“unpolarized”.)

Reference [9] constructs a GW background by summing the GWs produced by a
discrete set of unpolarized sources: the sources are circularly polarized,
corresponding to $\iota=0^{\circ}$ or $\iota=180^{\circ}$. Here, we improve
this simple GW model, assuming instead that $\cos\iota$ is uniformly
distributed on the interval $[-1,1]$. This corresponds to a physically
realistic random distribution of orbital inclination angles. To characterize
the resulting GW background, we compute the mean and the variance of the
Hellings and Downs (HD) correlation.

A brief outline of the paper is as follows. We begin in Sec. II with the GW
amplitude produced by a very distant circular binary system labeled by an
integer $j$. In Sec. III, we show how to treat the parameters as random
variables, and explain the relationship to the simpler ensemble of face-on and
face-off binary systems presented in [9]. In Sec. IV, as a “warm-up exercise”,
we compute the first and second moments and variance of a local quantity: the
time-averaged square strain $s$ defined on the first line of (10). We then
compare this quantity for the ensembles of unpolarized sources from [9] and
our ensemble of polarized sources. In Sec. V, we repeat these calculations of
the mean, second moment, and variance for the HD correlation between two
pulsars. Again, we compare it for the two ensembles. In Sec. VI, we define a
pulsar-averaged correlation at angle $\gamma$, and compute its mean, second
moment, and (cosmic) variance, and compare it for the two ensembles. In Sec.
VII, the results are compared to those of the standard Gaussian ensemble, and
we show that in a certain $N\to\infty$ limit our results approach those of the
Gaussian ensemble. This is followed by a short conclusion.

We use units in which the speed of light is unity, meaning that distances are
expressed in terms of the time that it would take light or GWs to propagate.

##  II GW Source Model

Here, we follow the approach of Sec. 3A of [9]. In a representative universe,
the GW background is produced by a discrete collection of $N$ distant point
sources. The sources are indexed by $j=1,\dots,N$, where $j=1$ is the closest
source, $j=2$ is the second closest, and so on. The GWs arrive at Earth from
the $j$’th source with propagation direction ${\bf\Omega}_{j}$. Hence, the
line of sight to each source is the unit-length vector $-{\bf\Omega}_{j}$,
though we often call ${\bf\Omega}_{j}$ the “source direction”.

In the model of [9], each source has circularly polarized gravitational
waveforms

$\displaystyle h^{+}_{j}(t)$ $\displaystyle={\cal{A}}_{j}\cos(\omega
t+\phi_{j})$ (1) $\displaystyle h^{\times}_{j}(t)$
$\displaystyle={\cal{A}}_{j}\sin(\omega t+\phi_{j}),$

where $h^{+}$ and $h^{\times}$ are the GW amplitudes in the respective
polarizations at Earth at time $t$. Here, $t$ is time at the solar system
barycenter (SSB) 111Since the Earth to SSB distance is about 500 seconds, and
the period of the GWs is orders of magnitude longer than this, the time at the
Earth and SSB may be treated as equivalent for this analysis., and $\omega$ is
an angular frequency. The amplitudes ${\cal{A}}_{j}\equiv{\cal{A}}/j^{1/3}$
fall off with distance, corresponding to a uniform density in 3-dimensional
Euclidean space (out to some maximum radius). The constant ${\cal{A}}$ is the
GW amplitude of the closest source, measured at the SSB.

The orbital frequency of a binary system increases with time as the system
radiates GWs and loses energy, becoming more tightly bound. Here, we assume
that this orbital frequency is evolving slowly on the time scale of PTA
observations. Thus, the GW frequency $\omega$ (nominally, twice the orbital
frequency) may be treated as a constant. We also assume that $\omega$ is
commensurate with the total observation time $T$, meaning that $\omega T/2\pi$
is integer. Note that, in this model, all $N$ sources emit GWs at the same
frequency $\omega$. Since the frequency is commensurate with $T$, they all
fall into a single frequency bin.

In the vicinity of the PTA, far from the emitting binary systems, the GWs are
weak, so the total GW field may be written as the sum of the GWs from the
individual sources:

$h_{ab}(t,{\bf x})\equiv\sum_{j=1}^{N}h_{j,ab}(t,{\bf x})\,.$ (2)

The GW emitted by the $j$’th source is

$h_{j,ab}(t,{\bf x})\equiv h^{+}_{j}(t-{\bf
x}\cdot{\bf\Omega}_{j})e^{+}_{ab}({\bf\Omega}_{j})+h^{\times}_{j}(t-{\bf
x}\cdot{\bf\Omega}_{j})e^{\times}_{ab}({\bf\Omega}_{j})\,,$ (3)

where $+$ and $\times$ label the two GW polarization states, and the traceless
and symmetric polarization tensors $e^{+}_{ab}$ and $e^{\times}_{ab}$ are
defined in Eq. D6 [9]. These polarization states are defined with respect to
the propagation direction ${\bf\Omega}_{j}$ of the GWs arriving from the
$j$’th point source. The definition uses an arbitrary direction orthogonal to
${\bf\Omega}$; our choice is that of Eqs. D6 and D7 of [9]. In the model of
(1), the $+$ and $\times$ waveforms have equal amplitudes, but are shifted in
phase by $90^{\circ}$. This is because [9] assumes, for simplicity, that each
of the individual point sources is circularly polarized. In what follows, for
convenience, we drop the limits from the summation if they are $j=1$ to $N$ as
in (2).

Note that these expressions are only valid in the region of space containing
Earth and the pulsars, a few thousand years in size. In that neighborhood they
are a good approximation: the sources are so distant that the deviations
between their spherical wavefronts and plane waves are much smaller than a
gravitational wavelength.

A binary system in a circular orbit with an inclination angle $\iota$ in the
range $(0^{\circ},180^{\circ})$ does not produce circularly polarized GWs. It
produces an elliptical polarization, with amplitudes [13, 14, 15]

$\displaystyle h^{+}_{j}(t)$
$\displaystyle=A_{j}\left[\frac{1}{2}(1+\cos^{2}\iota_{j})\cos
2\psi_{j}\,\cos(\omega t+\phi_{j})-\cos\iota_{j}\sin 2\psi_{j}\,\sin(\omega
t+\phi_{j})\right]$ (4) $\displaystyle h^{\times}_{j}(t)$
$\displaystyle=A_{j}\left[\frac{1}{2}(1+\cos^{2}\iota_{j})\sin
2\psi_{j}\,\cos(\omega t+\phi_{j})+\cos\iota_{j}\cos 2\psi_{j}\,\sin(\omega
t+\phi_{j})\right]\,,$

which should be substituted into (2). As in (1), the GW amplitudes
$A_{j}\equiv A/j^{1/3}$ correspond to a constant source density within a ball,
where $A$ is the amplitude of the closest source. Later, we will see that the
two different ensembles have the same squared GW strain at Earth if
${\cal{A}}=\sqrt{2/5}\,A$.

Equation (4) is similar to (1), but contains two additional parameters
$\iota_{j}$ and $\psi_{j}$, which describe aspects of the binary orbit. As
previously explained, $\iota\in[0^{\circ},180^{\circ}]$ is the angle between
the line of sight and a vector normal to the orbital plane. For
$\iota=0^{\circ}$ or $\iota=180^{\circ}$, the orbit appears circular when
viewed from Earth (i.e., projected on the plane of the sky). For any other
orbital inclination, the orbit appears to be a noncircular ellipse. The
parameter $\psi_{j}\in[0^{\circ},90^{\circ}]$ describes the orientation of
this ellipse. It is the angle between two vectors, after projection into a
plane orthogonal to ${\bf\Omega}_{j}$. The first is the major (long) axis of
the orbital ellipse, as seen from Earth. The second is $\bf m$, as defined in
Eq. D7 of [9].

Often the parameter $\psi$ does not appear. Indeed, if there is only one
binary system, then $\psi$ can be set to zero in (4). This follows from gauge
invariance: $\psi$ can be eliminated by selecting polarization axes which are
aligned with the ellipse. More precisely, the $+$ and $\times$ polarization
directions are defined by a pair of vectors ${\bf m}$ and ${\bf n}$, which are
orthogonal to the GW propagation direction ${\bf\Omega}$. If there is only a
single source, then picking the direction ${\bf m}$ along the line of the
major axis of the ellipse is equivalent to setting $\psi_{j}=0$. If there is
more than one source, then one could define vectors ${\bf m}^{\prime}$ and
${\bf n}^{\prime}$ on a per-source basis, to eliminate $\psi_{j}$. However, if
we do this, then the calculations that follow would contain new source-
dependent quantities ${\bf m}^{\prime}_{j}$ and ${\bf n}^{\prime}_{j}$; the
number of variables and corresponding complexity is not reduced.

##  III GW Ensemble

To compute statistical properties of the GW background, we consider an
ensemble of different universes, each containing $N$ sources. In each of these
universes (also called “realizations”), (4) defines the GW emission from the
$j$th source. Thus, each realization is defined by particular values of
$\phi_{j}$, ${\bf\Omega}_{j}$, $\iota_{j}$, and $\psi_{j}$, for $j=1,\dots,N$.
But, since we are ignorant of the values that these quantities have in our own
universe, we must resort to statistical inferences about “typical” members of
the ensemble. For this purpose, we treat these variables as random variables,
and use ensemble averages to predict the properties of a typical member of the
ensemble.

For binary systems with random orientations and positions in space, the $5N$
random variables are independent. The orbital phases $\phi_{j}$ are uniformly
distributed over $[0,2\pi]$. The directions to the sources ${\bf\Omega}_{j}$
are uniformly distributed over the two-sphere. The cosines of $\iota_{j}$ are
uniformly distributed on $[-1,1]$, and the polarization angles $\psi_{j}$ are
uniformly distributed over $[0,\pi/2]$.

The correct way to describe this ensemble is as “an ensemble of polarized
sources”, since each realization of the ensemble contains many elliptically
polarized GW sources (4) . However, we frequently employ the shorthand “the
polarized ensemble” or (for the contrasting case) “the unpolarized ensemble”.
[The circularly polarized GW sources in the latter are described by (1).] The
shorthand “polarized ensemble” can be misleading, because the ensemble itself
has no preferred direction or axis or polarization. The reader should bear in
mind that it is not the _ensemble_ which is polarized or unpolarized, but
rather the individual sources in each realization of that ensemble.

##  IV Squared-strain estimator $s$: mean, second moment, and variance

Following [9, 16] we first compute the expected value and variance of the
time-averaged squared strain. The calculation is similar to that used to find
the mean and variance of the HD correlation, but is simpler.

It is helpful to introduce a complex polarization basis

$e_{j,ab}=e_{ab}({\bf\Omega}_{j})\equiv
e^{+}_{ab}({\bf\Omega}_{j})-ie^{\times}_{ab}({\bf\Omega}_{j})\,,$ (5)

where $j$ labels the source and $a$ and $b$ are vector indices.
$e_{ab}({\bf\Omega}_{j})$ contains both the plus and cross components, and
corresponds to a left-circularly-polarized GW. The complex conjugate
$e^{*}_{j,ab}$ is the polarization of a right-circularly-polarized GW.

These complex polarization tensors form a null basis for a two-dimensional
vector space. From the definitions of $e^{+}_{ab}$ and $e^{\times}_{ab}$ given
in Eqs. D6 and D7 of [9], it follows that
$e^{+}_{ab}e^{+\,ab}=e^{\times}_{ab}e^{\times\,ab}=2$. This implies that

$\displaystyle e_{ab}e^{*\,ab}$ $\displaystyle=4$ (6) $\displaystyle
e_{ab}e^{ab}$ $\displaystyle=0$ $\displaystyle e^{*}_{ab}e^{*\,ab}$
$\displaystyle=0\,,$

where the third line follows from the second one.

It is also convenient to introduce several functions that will simplify the
notation. We define

$\displaystyle U_{j}$ $\displaystyle\equiv A_{j}(1+\cos\iota_{j})^{2}$ (7)
$\displaystyle V_{j}$ $\displaystyle\equiv A_{j}(1-\cos\iota_{j})^{2}\,,$

which only depend upon the orbital inclination angle. Using these, we also
define

$\displaystyle{\mathcal{U}}_{j}$ $\displaystyle\equiv U_{j}{\rm
e}^{2i\psi_{j}}$ (8) $\displaystyle{\mathcal{V}}_{j}$ $\displaystyle\equiv
V_{j}{\rm e}^{2i\psi_{j}}\,,$

which depend on both the orbital inclination and on the polarization angle.

We now compute the time average of the squared GW amplitude at Earth (${\bf
x}=0$). First, we substitute (7) and (8) into (4). The resulting expression is
inserted into (2), and (5) and its complex conjugate are used to replace
$e^{+}_{ab}$ and $e^{\times}_{ab}$. The GW strain at Earth becomes

$h_{ab}(t,0)=\frac{1}{8}\sum_{j}\left[\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e_{j,ab}^{*}\right){\rm
e}^{i(\omega
t+\phi_{j})}+\left({\mathcal{U}}^{*}_{j}\,e_{j,ab}^{*}+{\mathcal{V}}_{j}\,e_{j,ab}\right){\rm
e}^{-i(\omega t+\phi_{j})}\right]\,.$ (9)

We now square and contract this quantity, and time-average, obtaining

$\displaystyle s\equiv$ $\displaystyle\,\overline{h_{ab}(t,0)h^{ab}(t,0)}$
(10) $\displaystyle=$
$\displaystyle\,\frac{1}{32}\sum_{j}\Bigl{[}\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{j}\,e^{*\,ab}_{j}+{\mathcal{V}}_{j}\,e^{ab}_{j}\right)\Bigr{]}+$
$\displaystyle\,\frac{1}{64}\sum_{j\neq
k}\Bigl{[}\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{k}\,e^{*\,ab}_{k}+{\mathcal{V}}_{k}\,e^{ab}_{k}\right){\rm
e}^{i(\phi_{j}-\phi_{k})}+$
$\displaystyle\,\quad\quad\quad\left({\mathcal{U}}^{*}_{j}\,e^{*}_{j,ab}+{\mathcal{V}}_{j}\,e_{j,ab}\right)\left({\mathcal{U}}_{k}\,e^{ab}_{k}+{\mathcal{V}}^{*}_{k}\,e^{*\,ab}_{k}\right){\rm
e}^{-i(\phi_{j}-\phi_{k})}\Bigr{]}\,.$

The average over time $\overline{Q}\equiv\frac{1}{T}\int_{-T/2}^{T/2}Q(t)dt$
has eliminated half of the cross terms, and we have grouped the others into
“diagonal” terms with $j=k$ and “off-diagonal” terms with $j\neq k$.

In this paper, we use different indexing and summation conventions for
tangent-space indices such as $a$, $b$, $c$, and $d$, and for source labels
such as $j$, $k$, $\ell$, and $m$. We adopt the Einstein summation convention
for repeated vector indices. For example, in (10) the indices $a$ and $b$ are
implicitly summed from 1 to 3\. In contrast, sums over GW sources are
explicit. Double sums over sources are indicated with a single summation
symbol as $\sum_{j,k}$. If the terms where $j=k$ are dropped, then the double
sums are denoted by $\sum_{j\neq k}$.

Expression (10) for $s$ holds for any realization in the ensemble, and is a
function of the $5N$ variables $\phi_{j}$, ${\bf\Omega}_{j}$, $\iota_{j}$, and
$\psi_{j}$. For a given ensemble, we can calculate the expected value of _any_
function of those $5N$ variables, including $s$ or powers of $s$ or other
functions of $s$.

If $Q$ is a function or functional of the $5N$ variables, then the ensemble
average of $Q$ is defined by

$\langle Q\rangle\equiv\prod_{j=1}^{N}\left(\int_{0}^{2\pi}\frac{{\rm
d}\phi_{j}}{2\pi}\right)\prod_{k=1}^{N}\left(\int_{0}^{2\pi}\frac{{\rm
d}\psi_{k}}{2\pi}\right)\prod_{\ell=1}^{N}\left(\int_{-1}^{1}\frac{{\rm
d}\cos\iota_{\ell}}{2}\right)\prod_{m=1}^{N}\left(\int\frac{{\rm
d}{\bf\Omega}_{m}}{4\pi}\right)Q(\phi_{1},\dots,{\bf\Omega}_{N}).$ (11)

In practice, we carry out these integrals in four stages. If we first carry
out the $N$ integrals over $\phi_{j}$ in (11), the result may be denoted by
$\langle Q\rangle_{\phi}$. Doing the remaining three sets of integrals then
gives the full average. Thus, $\langle
Q\rangle\equiv\langle\langle\langle\langle
Q\rangle_{\phi}\rangle_{\psi}\rangle_{\iota}\rangle_{\bf\Omega}$, where each
of the separate averages refers to one set of the integrals in (11) 222The
physical range over which $\psi$ is averaged is $[0,\pi/2]$. However, since
all of the relevant functions are periodic on this interval, averaging over a
larger interval, as we do in (11), gives the same result..

###  IV.1  First moment of $s$

The first moment of $s$ is easily computed. It follows immediately from (11)
that $\langle{\rm e}^{i(\phi_{j}-\phi_{k})}\rangle_{\phi}$ vanishes if $j\neq
k$ and is unity if $j=k$. Hence, if we take the $\phi_{j}$ average of (10),
only the diagonal term survives, giving

$\langle
s\rangle_{\phi}=\frac{1}{32}\sum_{j}\Bigl{[}\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{j}\,e^{*\,ab}_{j}+{\mathcal{V}}_{j}\,e^{ab}_{j}\right)\Bigr{]}.$
(12)

We next average over $\psi_{j}$, carrying out the next set of integrals in
(11). The rhs of (12) only depends upon $\psi_{j}$ through ${\mathcal{U}}_{j}$
and ${\mathcal{V}}_{j}$, as defined by (8). The
${\mathcal{U}}_{j}{\mathcal{V}}_{j}$ and
${\mathcal{U}}_{j}^{*}{\mathcal{V}}_{j}^{*}$ terms are proportional to
$\exp(\pm 4i\psi_{j})$ and integrate to zero. The surviving
${\mathcal{U}}_{j}{\mathcal{U}}_{j}^{*}$ and
${\mathcal{V}}_{j}^{*}{\mathcal{V}}_{j}$ terms give

$\langle\langle
s\rangle_{\phi}\rangle_{\psi}=\frac{1}{32}\sum_{j}\left(U_{j}^{2}+V_{j}^{2}\right)e^{*}_{j,ab}\,e^{ab}_{j}\,.$
(13)

From (7), one can see that the factor $U_{j}^{2}+V_{j}^{2}$ only depends upon
the inclination angle $\iota_{j}$. Evaluating the third set of integrals in
(11) to average over $\iota_{j}$ gives

$\langle
U_{j}^{2}+V_{j}^{2}\rangle_{\iota}=\frac{1}{2}A_{j}^{2}\int_{-1}^{1}{\rm
d}x\left((1+x)^{4}+(1-x)^{4}\right)=\frac{32}{5}A_{j}^{2},$ (14)

where $x=\cos\iota_{j}$. Thus, the $\iota$-average of (13) is

$\langle\langle\langle
s\rangle_{\phi}\rangle_{\psi}\rangle_{\iota}=\frac{1}{5}\sum_{j}A_{j}^{2}\,e^{*}_{j,ab}\,e^{ab}_{j}\,.$
(15)

The final quantities to average are the polarization tensors, which only
depend upon the direction ${\bf\Omega}_{j}$. In fact, it follows immediately
from their definition (see Eqs. D6 and D7 of [9]) that
$e^{*}_{ab}({\bf\Omega})e^{ab}({\bf\Omega})=4$ is independent of ${\bf\Omega}$
(so the label “$j$” is not needed). Hence, from (15) we obtain the ensemble
average

$\langle s\rangle_{\rm pol}=\frac{4}{5}\sum_{j}A_{j}^{2}.$ (16)

This completes the calculation of the first moment for the polarized ensemble
(4).

###  IV.2  Second moment and variance of $s$

We next compute the second moment of the time-averaged squared strain $s$. We
will use this to determine the variance of $s$. To evaluate $\langle
s^{2}\rangle$, we square $s$ as given in (10) and then use (11) to determine
the ensemble average. This is very similar to the calculation of the first
moment, so we can be more concise.

We begin by computing the average over the phases $\phi_{j}$. First, square
(10) to obtain three terms. These are (i) the square of the diagonal sum (ii)
the square of the off-diagonal sum and (iii) twice their cross-product.
Averaging over the random phases $\phi_{j}$ eliminates (iii), and the
remaining terms are

$\displaystyle\langle
s^{2}\rangle_{\phi}=\Biggl{(}\frac{1}{32}\sum_{j}\Bigl{[}\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{j}\,e^{*\,ab}_{j}+{\mathcal{V}}_{j}\,e^{ab}_{j}\right)\Bigr{]}\Biggr{)}^{2}+\frac{1}{4096}\sum_{j\neq
k}\sum_{\ell\neq m}\Bigl{[}$ (17)
$\displaystyle\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{k}\,e^{*\,ab}_{k}+{\mathcal{V}}_{k}\,e^{ab}_{k}\right)\left({\mathcal{U}}_{\ell}\,e_{\ell,cd}+{\mathcal{V}}^{*}_{\ell}\,e^{*}_{\ell,cd}\right)\left({\mathcal{U}}^{*}_{m}\,e^{*\,cd}_{m}+{\mathcal{V}}_{m}\,e^{cd}_{m}\right)\langle{\rm
e}^{i(\phi_{j}-\phi_{k}+\phi_{\ell}-\phi_{m})}\rangle_{\phi}+$
$\displaystyle\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{k}\,e^{*\,ab}_{k}+{\mathcal{V}}_{k}\,e^{ab}_{k}\right)\left({\mathcal{U}}^{*}_{\ell}\,e^{*}_{\ell,cd}+{\mathcal{V}}_{\ell}\,e_{\ell,cd}\right)\left({\mathcal{U}}_{m}\,e^{cd}_{m}+{\mathcal{V}}^{*}_{m}\,e^{*\,cd}_{m}\right)\langle{\rm
e}^{i(\phi_{j}-\phi_{k}-\phi_{\ell}+\phi_{m})}\rangle_{\phi}+$
$\displaystyle\left({\mathcal{U}}^{*}_{j}\,e^{*}_{j,ab}+{\mathcal{V}}_{j}\,e_{j,ab}\right)\left({\mathcal{U}}_{k}\,e^{ab}_{k}+{\mathcal{V}}^{*}_{k}\,e^{*\,ab}_{k}\right)\left({\mathcal{U}}_{\ell}\,e_{\ell,cd}+{\mathcal{V}}^{*}_{\ell}\,e^{*}_{\ell,cd}\right)\left({\mathcal{U}}^{*}_{m}\,e^{*\,cd}_{m}+{\mathcal{V}}_{m}\,e^{cd}_{m}\right)\langle{\rm
e}^{-i(\phi_{j}-\phi_{k}-\phi_{\ell}+\phi_{m})}\rangle_{\phi}+$
$\displaystyle\left({\mathcal{U}}^{*}_{j}\,e^{*}_{j,ab}+{\mathcal{V}}_{j}\,e_{j,ab}\right)\left({\mathcal{U}}_{k}\,e^{ab}_{k}+{\mathcal{V}}^{*}_{k}\,e^{*\,ab}_{k}\right)\left({\mathcal{U}}^{*}_{\ell}\,e^{*}_{\ell,cd}+{\mathcal{V}}_{\ell}\,e_{\ell,cd}\right)\left({\mathcal{U}}_{m}\,e^{cd}_{m}+{\mathcal{V}}^{*}_{m}\,e^{*\,cd}_{m}\right)\langle{\rm
e}^{-i(\phi_{j}-\phi_{k}+\phi_{\ell}-\phi_{m})}\rangle_{\phi}\Bigr{]}\,.$

The averages of the exponentials are unity if the exponent vanishes, else they
are zero. This implies

$\displaystyle\langle{\rm e}^{\pm
i(\phi_{j}-\phi_{k}+\phi_{\ell}-\phi_{m})}\rangle_{\phi}$
$\displaystyle=\delta_{jm}\delta_{k\ell}$ (18) $\displaystyle\langle{\rm
e}^{\pm i(\phi_{j}-\phi_{k}-\phi_{\ell}+\phi_{m})}\rangle_{\phi}$
$\displaystyle=\delta_{j\ell}\delta_{km}\,,$

provided that $j\neq k$ and $\ell\neq m$. Substituting (18) into (17)
eliminates one of the double sums. The remaining double sum has four terms,
which are all equal. Finally, expressing the square of the single sum as
diagonal and off-diagonal terms gives

$\displaystyle\langle s^{2}\rangle_{\phi}=$
$\displaystyle\frac{1}{1024}\sum_{j}\left[\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{j}\,e^{*\,ab}_{j}+{\mathcal{V}}_{j}\,e^{ab}_{j}\right)\right]^{2}+$
(19) $\displaystyle\frac{1}{1024}\sum_{j\neq
k}\Bigl{[}\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{j}\,e^{*\,ab}_{j}+{\mathcal{V}}_{j}\,e^{ab}_{j}\right)\left({\mathcal{U}}_{k}\,e_{k,cd}+{\mathcal{V}}^{*}_{k}\,e^{*}_{k,cd}\right)\left({\mathcal{U}}^{*}_{k}\,e^{*\,cd}_{k}+{\mathcal{V}}_{k}\,e^{cd}_{k}\right)+$
$\displaystyle\,\quad\quad\quad\quad\left({\mathcal{U}}_{j}\,e_{j,ab}+{\mathcal{V}}^{*}_{j}\,e^{*}_{j,ab}\right)\left({\mathcal{U}}^{*}_{k}\,e^{*\,ab}_{k}+{\mathcal{V}}_{k}\,e^{ab}_{k}\right)\left({\mathcal{U}}_{k}\,e_{k}^{cd}+{\mathcal{V}}^{*}_{k}\,e^{*\,cd}_{k}\right)\left({\mathcal{U}}^{*}_{j}\,e^{*}_{j,cd}+{\mathcal{V}}_{j}\,e_{j,cd}\right)\Bigr{]}\,.$

This completes the averaging over $\phi_{j}$.

Next, we average (19) over $\psi_{j}$. From (8) we see that both
${\mathcal{U}}$ and ${\mathcal{V}}$ depend on $\psi$ only through the factor
$\exp{(2i\psi)}$. Thus, after averaging over $\psi$, terms with different
numbers of conjugated and unconjugated ${\mathcal{U}}$ and ${\mathcal{V}}$
will vanish. For example, ${\mathcal{U}}{\mathcal{U}}^{*}{\mathcal{V}}^{2}$
has one conjugated and three unconjugated; it averages to zero. The terms
inside the single sum of (19) that will survive the averaging are

$\left({\mathcal{U}}^{2}_{j}{\mathcal{U}}_{j}^{*\,2}+{\mathcal{V}}^{2}_{j}{\mathcal{V}}_{j}^{*\,2}\right)\left(e_{j}^{ab}e^{*}_{j,ab}\right)^{2}+2\left({\mathcal{U}}_{j}{\mathcal{U}}_{j}^{*}{\mathcal{V}}_{j}{\mathcal{V}}_{j}^{*}\right)\left(\left(e_{j}^{ab}e^{*}_{j,ab}\right)^{2}+e_{j}^{ab}e_{j,ab}e_{j}^{*\,cd}e^{*}_{j,cd}\right)\,.$
(20)

The terms inside the double sum of (19) that will survive the averaging are

$\displaystyle\left({\mathcal{U}}_{j}{\mathcal{U}}_{j}^{*}+{\mathcal{V}}_{j}{\mathcal{V}}_{j}^{*}\right)e_{j}^{ab}e^{*}_{j,ab}\left({\mathcal{U}}_{k}{\mathcal{U}}_{k}^{*}+{\mathcal{V}}_{k}{\mathcal{V}}_{k}^{*}\right)e_{k}^{cd}e^{*}_{k,cd}\,+\left({\mathcal{U}}_{j}{\mathcal{U}}_{j}^{*}e_{j,ab}e^{*}_{j,cd}+{\mathcal{V}}_{j}{\mathcal{V}}_{j}^{*}e_{j,cd}e^{*}_{j,ab}\right)\left({\mathcal{U}}_{k}{\mathcal{U}}_{k}^{*}e_{k}^{ab}e_{k}^{*\,cd}+{\mathcal{V}}_{k}{\mathcal{V}}_{k}^{*}e_{k}^{cd}e_{k}^{*\,ab}\right)\,.$
(21)

Combining these, using (6) to evaluate the contractions of the polarization
tensors, and making use of (8), the $\psi$-average of (19) becomes

$\displaystyle\langle\langle s^{2}\rangle_{\phi}\rangle_{\psi}=\frac{1}{64}$
$\displaystyle\sum_{j}\Bigl{[}U_{j}^{4}+V_{j}^{4}+2U_{j}^{2}V_{j}^{2}\Bigr{]}$
(22) $\displaystyle+\frac{1}{1024}$ $\displaystyle\sum_{j\neq
k}\Bigl{[}16\left(U_{j}^{2}+V_{j}^{2}\right)\left(U_{k}^{2}+V_{k}^{2}\right)+\left(U_{j}^{2}e_{j,ab}\,e^{*}_{j,cd}+V_{j}^{2}e^{*}_{j,ab}\,e_{j,cd}\right)\left(U_{k}^{2}e^{*\,ab}_{k}\,e^{cd}_{k}+V_{k}^{2}e^{ab}_{k}\,e^{*\,cd}_{k}\right)\Bigr{]}\,.$

Note that the condition $j\neq k$ simplifies the evaluation of the
$\psi$-average of the double sum: the average of the product becomes the
product of the averages.

We next average over $\iota_{j}$, as defined in (11). From (7), the only
$\iota$ dependence is through the functions $U$ and $V$. For example,

$\langle U_{j}^{2}\rangle_{\iota}=\frac{1}{2}\int_{-1}^{1}U_{j}^{2}\,\,{\rm
d}\cos\iota_{j}=\frac{1}{2}A_{j}^{2}\int_{-1}^{1}(1+x)^{4}\,{\rm
d}x=\frac{16}{5}A_{j}^{2}\,,$ (23)

where $x=\cos\iota_{j}$. For this calculation, the following averages,
computed in the same way, are sufficient:

$\displaystyle\langle U_{j}^{2}\rangle_{\iota}=\langle
V_{j}^{2}\rangle_{\iota}$ $\displaystyle=\frac{16}{5}A_{j}^{2}$ (24)
$\displaystyle\langle U_{j}^{4}\rangle_{\iota}=\langle
V_{j}^{4}\rangle_{\iota}$ $\displaystyle=\frac{256}{9}A_{j}^{4}$
$\displaystyle\langle U_{j}^{2}V_{j}^{2}\rangle_{\iota}$
$\displaystyle=\frac{128}{315}A_{j}^{4}\,.$

Using (24) to average (22) over $\iota$ gives

$\displaystyle\langle\langle\langle
s^{2}\rangle_{\phi}\rangle_{\psi}\rangle_{\iota}=$
$\displaystyle\,\,\frac{284}{315}\sum_{j}A_{j}^{4}+\frac{1}{100}\sum_{j\neq
k}A_{j}^{2}A_{k}^{2}\Bigl{[}64+\left(e_{j,ab}e_{j,cd}^{*}+e^{*}_{j,ab}e_{j,cd}\right)\left(e_{k}^{ab}e_{k}^{*\,cd}+e^{*\,ab}_{k}e_{k}^{cd}\right)\Bigr{]}\,.$
(25)

Again, because $j\neq k$, we write the average of the product of $j$\- and
$k$-dependent terms as the product of the averages.

The final step to compute the second moment is to average (25) over the source
directions ${\bf\Omega}_{j}$. For this purpose, we employ the spherical
average of polarization tensors

$\displaystyle\eta_{abcd}$ $\displaystyle\equiv\frac{1}{4\pi}\int
e_{ab}({\bf\Omega})e^{*}_{cd}({\bf\Omega})\,{\rm d}{\bf\Omega}$ (26)
$\displaystyle=-\frac{4}{15}\,\delta_{ab}\delta_{cd}+\frac{2}{5}\bigl{(}\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}\bigr{)}\,,$

which is derived in Eqs. C33 and C34 of [9]. Here, $\delta_{ab}$ denotes the
3-dimensional Kronecker delta. Using (26) to average (25) over the source
directions ${\bf\Omega}_{j}$ gives

$\displaystyle\langle s^{2}\rangle$
$\displaystyle=\frac{284}{315}\sum_{j}A_{j}^{4}+\frac{1}{25}\sum_{j\neq
k}A_{j}^{2}A_{k}^{2}\left(16+\eta_{abcd}\eta^{abcd}\right)$ (27)
$\displaystyle=\frac{284}{315}\sum_{j}A_{j}^{4}+\frac{96}{125}\sum_{j\neq
k}A_{j}^{2}A_{k}^{2}$
$\displaystyle=\frac{1052}{7875}\sum_{j}A_{j}^{4}+\frac{96}{125}\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}^{2}\,.$

Since $j\neq k$, the first line follows from carrying out the final averaging
in (11) _separately_ on the $j$ and the $k$ terms. To obtain the second line,
we use $\eta_{abcd}\eta^{abcd}=16/5$, obtained from (26). The third line
follows from

$\sum_{j\neq
k}A_{j}^{2}A_{k}^{2}=\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}^{2}-\sum_{j}A_{j}^{4}$
(28)

(which holds for _any_ set of amplitudes). This completes the calculation of
the second moment of $s$.

The variance $\sigma^{2}_{s}\equiv\langle s^{2}\rangle-\langle s\rangle^{2}$
is obtained from the first (16) and second (27) moments. The mean and variance
of the polarized ensemble (4) are

$\displaystyle\langle s\rangle_{\rm pol}$
$\displaystyle=\frac{4}{5}\sum_{j}A_{j}^{2}$ (29)
$\displaystyle\sigma^{2}_{s,{\rm pol}}$
$\displaystyle=\frac{1052}{7875}\sum_{j}A_{j}^{4}+\frac{16}{125}\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}^{2}\,,$

where the sum is over the $N$ sources, with amplitudes given after (4). We now
compare this mean and variance to those of the unpolarized ensemble (1).

###  IV.3 Comparison of polarized and unpolarized ensembles

The mean and variance of $s$ for the unpolarized ensemble (1) are given in
Eqs. 3.15 and 3.23 of [9] as

$\displaystyle\langle s\rangle_{\rm unpol}$
$\displaystyle=2\sum_{j}{\cal{A}}_{j}^{2}$ (30)
$\displaystyle\sigma^{2}_{s,{\rm unpol}}$
$\displaystyle=-\frac{4}{5}\sum_{j}{\cal{A}}_{j}^{4}+\frac{4}{5}\Bigl{(}\sum_{j}{\cal{A}}_{j}^{2}\Bigr{)}^{2}\,.$

Both ensembles are characterized by a single number: the amplitude of the
closest source to Earth, denoted respectively by $A$ and ${\cal{A}}$.

In order to compare the two ensembles “on equal footing”, we select these
amplitudes so that both ensembles have the same mean squared GW amplitude
$\langle s\rangle$ at Earth. Comparing (29) and (30), this implies that

${\cal{A}}^{2}=\frac{2}{5}A^{2}\iff{\cal{A}}_{j}=\sqrt{\frac{2}{5}}A_{j}.$
(31)

It follows immediately from (29) and (30) that for a given mean, the variance
of the polarized ensemble is always larger than that of the unpolarized
ensemble.

These expressions include $N$ sources labeled by $j=1,\dots N$. Because
$A_{j}$ and ${\cal{A}}_{j}$ fall off $\propto j^{-1/3}$, the sums
$\sum_{j}A_{j}^{2}$ and $\sum_{j}{\cal{A}}_{j}^{2}$ diverge as the number of
GW sources $N\to\infty$. In contrast, the sums $\sum_{j}A_{j}^{4}$ and
$\sum_{j}{\cal{A}}_{j}^{4}$ converge to $A^{4}\zeta(4/3)\approx 3.6A^{4}$ and
${\cal{A}}^{4}\zeta(4/3)\approx 3.6{\cal{A}}^{4}$ respectively, where $\zeta$
is the Riemann zeta function. Thus, in the limit of many sources, the
fractional difference in the variances approaches zero.

##  V Mean and variance of the HD correlation

The HD correlation $\rho$ between a pair of pulsars is defined in terms of
their redshifts $Z$. Consider a pulsar at spatial location ${\bf
x}=L{\hat{\bf{p}}}$, where $L$ is the distance to the pulsar and
${\hat{\bf{p}}}$ is a unit-length vector. The redshift induced by the $j$’th
source at time $t$ is

$Z_{j}(t)\equiv\frac{1}{2}\frac{{\hat{p}}^{a}{\hat{p}}^{b}}{1+{\bf\Omega}_{j}\cdot{\hat{\bf{p}}}}\Bigl{[}h_{j,ab}(t,0)-\chi
h_{j,ab}(t-L,L{\hat{\bf{p}}})\Bigr{]}\,,$ (32)

where the terms in square brackets are called the “Earth term” and the “pulsar
term” respectively. For a PTA, the correct value of $\chi$ is $\chi=1$, but we
retain this dependence to better understand the effects of the pulsar term.
The GW amplitudes $h_{j,ab}(t,{\bf x})$ are defined in (3) as a linear
combination of the two polarizations $e_{ab}^{+}({\bf\Omega}_{j})$ and
$e_{ab}^{\times}({\bf\Omega}_{j})$.

The GW amplitude of the $j$’th source $h_{j,ab}(t,{\bf x})$ may also be
expressed in other bases, such as the complex polarization basis $e_{j,ab}$
introduced in (5). For what follows, it is convenient to introduce the rotated
basis

${\rm e}^{2i\psi_{j}}e_{j,ab}.$ (33)

We can then write $h_{j,ab}(t,{\bf x})$ as a linear combination of (33) and
its complex conjugate. Defining the complex GW amplitude

$h_{j}(t)\equiv\frac{1}{2}A_{j}\Bigl{[}\left(1+\cos^{2}\iota_{j}\right)\cos(\omega
t+\phi_{j})+i\,2\cos\iota_{j}\sin(\omega t+\phi_{j})\Bigr{]}\,,$ (34)

it follows from (3), (4) and (5) that

$h_{j,ab}(t,{\bf x})=\frac{1}{2}\Bigl{[}h_{j}(t-{\bf
x}\cdot{\bf\Omega}_{j}){\rm e}^{2i\psi_{j}}e_{j,ab}+h_{j}^{*}(t-{\bf
x}\cdot{\bf\Omega}_{j}){\rm e}^{-2i\psi_{j}}e^{*}_{j,ab}\Bigr{]}\,.$ (35)

This can be substituted into (32) to determine the redshift $Z$. In effect,
for a given GW propagation direction, we have replaced the two degrees of
freedom contained in $h_{j,ab}$ (the amplitudes of the cross and plus
polarizations) with a single complex amplitude $h_{j}$.

Since the induced redshift (32) depends upon the difference between Earth
($t,{\bf x}=0$) and pulsar ($t-L,{\bf x}=L{\hat{\bf{p}}}$) terms, we introduce
the differences

$\displaystyle\Delta h_{j}(t)\equiv$ $\displaystyle h_{j}(\tau-{\bf
x}\cdot{\bf\Omega}_{j})\big{|}_{\tau=t,{\bf x}=0}-h_{j}(\tau-{\bf
x}\cdot{\bf\Omega}_{j})\big{|}_{\tau=t-L,{\bf x}=L{\hat{\bf{p}}}}$ (36)
$\displaystyle=$
$\displaystyle\tfrac{1}{2}A_{j}\left(1+\cos^{2}\iota_{j}\right)\Bigl{(}\cos(\omega
t+\phi_{j})-\chi\cos(\omega t+\phi_{j}-\Delta_{j})\Bigr{)}$
$\displaystyle\quad\quad\,\,\,+iA_{j}\cos\iota_{j}\Bigl{(}\sin(\omega
t+\phi_{j})-\chi\sin(\omega t+\phi_{j}-\Delta_{j})\Bigr{)}$ $\displaystyle=$
$\displaystyle\tfrac{1}{4}U_{j}\left[1-\chi{\rm e}^{-i\Delta_{j}}\right]{\rm
e}^{i(\omega t+\phi_{j})}+\tfrac{1}{4}V_{j}\left[1-\chi{\rm
e}^{i\Delta_{j}}\right]{\rm e}^{-i(\omega t+\phi_{j})}\,.$

The second line follows from (34) and the definition

$\Delta_{j}\equiv\omega L(1+{\bf\Omega}_{j}\cdot{\hat{\bf{p}}}),$ (37)

and the third line from the definitions of $U_{j},V_{j}$ given in (7).
Starting with (32), and using (35) and the first line of (36), the induced
redshift may be written as

$Z_{j}(t)=\frac{1}{4}\frac{{\hat{p}}^{a}{\hat{p}}^{b}}{1+{\bf\Omega}_{j}\cdot{\hat{\bf{p}}}}\Bigl{[}\Delta
h_{j}(t){\rm e}^{2i\psi_{j}}e_{j,ab}+\Delta h^{*}_{j}(t){\rm
e}^{-2i\psi_{j}}e^{*}_{j,ab}\Bigr{]}\,.$ (38)

Using the final line of (36), this becomes

$\displaystyle Z_{j}(t)=$
$\displaystyle\frac{1}{16}\frac{{\hat{p}}^{a}{\hat{p}}^{b}}{1+{\bf\Omega}_{j}\cdot{\hat{\bf{p}}}}\Bigl{[}\Bigl{(}U_{j}\left(1-\chi{\rm
e}^{-i\Delta_{j}}\right){\rm e}^{i(\omega t+\phi_{j})}\\!\\!\\!$
$\displaystyle+$ $\displaystyle V_{j}\left(1-\chi{\rm
e}^{i\Delta_{j}}\right){\rm e}^{-i(\omega t+\phi_{j})}\Bigr{)}{\rm
e}^{2i\psi_{j}}e_{j,ab}+$ (39)
$\displaystyle\quad\quad\quad\quad\quad\quad\,\Bigl{(}U_{j}\left(1-\chi{\rm
e}^{i\Delta_{j}}\right){\rm e}^{-i(\omega t+\phi_{j})}\\!\\!\\!$
$\displaystyle+$ $\displaystyle V_{j}\left(1-\chi{\rm
e}^{-i\Delta_{j}}\right){\rm e}^{i(\omega t+\phi_{j})}\Bigr{)}{\rm
e}^{-2i\psi_{j}}e_{j,ab}^{*}\Bigr{]}\,.$

This expression for the redshift can be made more compact by writing it in
terms of “antenna pattern” functions.

For a plane GW with direction ${\bf\Omega}$ and polarization $A=+,\times$, the
(real) antenna pattern function for a given pulsar direction ${\hat{\bf{p}}}$
is defined as

$F^{A}({\bf\Omega})\equiv\frac{1}{2}\frac{{\hat{p}}^{a}{\hat{p}}^{b}}{1+{\bf\Omega}\cdot{\hat{\bf{p}}}}e^{A}_{ab}({\bf\Omega})\,.$
(40)

The two polarizations can be combined into a single complex quantity by
defining a complex antenna pattern function

$F_{j}\equiv
F^{+}({\bf\Omega}_{j})-iF^{\times}({\bf\Omega}_{j})=\frac{1}{2}\frac{{\hat{p}}^{a}{\hat{p}}^{b}}{1+{\bf\Omega}_{j}\cdot{\hat{\bf{p}}}}e_{ab}({\bf\Omega}_{j})\,,$
(41)

where the second equality follows from the definition of the complex
polarization tensor $e_{ab}$ in (5). If we use (41) to rewrite (39), then the
only terms that depend upon the pulsar position ${\hat{\bf{p}}}$ are $F_{j}$
and $(1-\chi{\rm e}^{i\Delta_{j}})$, and their complex conjugates. Thus, we
define

$f_{j}\equiv 1-\chi{\rm e}^{i\Delta_{j}}\,,$ (42)

where $\Delta_{j}$ is defined in (37) and depends upon the source direction
${\bf\Omega}_{j}$. (The term equal to one is the Earth term and the term
proportional to $\chi$ is the pulsar term.) Using (41) and (42), we can write
(39) as a linear combination of antenna pattern functions

$Z_{j}(t)=\frac{1}{8}\left(U_{j}f_{j}^{*}{\rm e}^{i(\omega
t+\phi_{j})}+V_{j}f_{j}{\rm e}^{-i(\omega t+\phi_{j})}\right){\rm
e}^{2i\psi_{j}}F_{j}+\frac{1}{8}\left(U_{j}f_{j}{\rm e}^{-i(\omega
t+\phi_{j})}+V_{j}f^{*}_{j}{\rm e}^{i(\omega t+\phi_{j})}\right){\rm
e}^{-2i\psi_{j}}F_{j}^{*}\,,$ (43)

which is the redshift induced by the $j$’th source.

To compute the HD correlation, we need the total redshift induced by all $N$
GW sources. For weak gravitational fields, the redshifts add linearly. Thus,
the total redshift at the time $t$ is the sum of the contributions (43) from
each source:

$\displaystyle Z_{\hat{\bf{p}}}(t)$ $\displaystyle\equiv\sum_{j}Z_{j}(t)$ (44)
$\displaystyle=\frac{1}{8}\sum_{j}\left[\left(U_{j}{\rm
e}^{2i\psi_{j}}F_{j}+V_{j}{\rm e}^{-2i\psi_{j}}F_{j}^{*}\right)f_{j}^{*}{\rm
e}^{i(\omega t+\phi_{j})}+\left(U_{j}{\rm e}^{-2i\psi_{j}}F_{j}^{*}+V_{j}{\rm
e}^{2i\psi_{j}}F_{j}\right)f_{j}{\rm e}^{-i(\omega t+\phi_{j})}\right]\,.$

We have regrouped the terms to isolate the positive- and negative-frequency
terms. Note that the redshift depends upon the direction to the pulsar via the
functions $f$ and $F$; we indicate this with the subscript ${\hat{\bf{p}}}$.

The correlation between two distinct pulsars with directions
${\hat{\bf{p}}}_{1}$ and ${\hat{\bf{p}}}_{2}$ is defined by

$\rho\equiv\overline{Z_{{\hat{\bf{p}}}_{1}}(t)Z_{{\hat{\bf{p}}}_{2}}(t)},$
(45)

where the overbar [defined immediately after (10)] denotes a time average. We
use (44) in (45) and average over time. Because $\omega$ is assumed to be an
integer multiple of $2\pi/T$, only the cross terms between positive and
negative frequency survive. We obtain

$\displaystyle\rho$
$\displaystyle=\frac{1}{64}\sum_{j,k}\Bigl{[}\left(U_{j}{\rm
e}^{2i\psi_{j}}F_{1,j}+V_{j}{\rm
e}^{-2i\psi_{j}}F^{*}_{1,j}\right)\left(U_{k}{\rm
e}^{-2i\psi_{k}}F_{2,k}^{*}+V_{k}{\rm
e}^{2i\psi_{k}}F_{2,k}\right)f_{1,j}^{*}f_{2,k}{\rm
e}^{i(\phi_{j}-\phi_{k})}+$ (46)
$\displaystyle\quad\quad\quad\quad\,\,\left(U_{j}{\rm
e}^{-2i\psi_{j}}F^{*}_{1,j}+V_{j}{\rm
e}^{2i\psi_{j}}F_{1,j}\right)\left(U_{k}{\rm e}^{2i\psi_{k}}F_{2,k}+V_{k}{\rm
e}^{-2i\psi_{k}}F^{*}_{2,k}\right)f_{1,j}f^{*}_{2,k}{\rm
e}^{-i(\phi_{j}-\phi_{k})}\Bigr{]}\,.$

The dependence of the antenna patterns $f$ and $F$ on pulsar directions
${\hat{\bf{p}}}_{1}$ and ${\hat{\bf{p}}}_{2}$ is denoted via the subscripts
$1$ and $2$.

We next compute the mean and variance of $\rho$ along the same lines as used
for $s$ in Secs. IV.1 and IV.2.

###  V.1  First moment of $\rho$

We proceed as in Sec. IV.1, evaluating the first set of integrals in (11) to
average $\rho$ over the phases $\phi_{j}$. We obtain

$\displaystyle\langle\rho\rangle_{\phi}=\frac{1}{64}\sum_{j}\Bigl{[}$
$\displaystyle\left(U_{j}{\rm e}^{2i\psi_{j}}F_{1,j}+V_{j}{\rm
e}^{-2i\psi_{j}}F^{*}_{1,j}\right)\left(U_{j}{\rm
e}^{-2i\psi_{j}}F_{2,j}^{*}+V_{j}{\rm
e}^{2i\psi_{j}}F_{2,j}\right)f_{1,j}^{*}f_{2,j}+$ (47)
$\displaystyle\left(U_{j}{\rm e}^{-2i\psi_{j}}F^{*}_{1,j}+V_{j}{\rm
e}^{2i\psi_{j}}F_{1,j}\right)\left(U_{j}{\rm e}^{2i\psi_{j}}F_{2,j}+V_{j}{\rm
e}^{-2i\psi_{j}}F^{*}_{2,j}\right)f_{1,j}f^{*}_{2,j}\Bigr{]}\,,$

where only the terms with $j=k$ remain.

We next average over $\psi$. Expanding the rhs of (47) gives terms
proportional to ${\rm e}^{4i\psi_{j}}$, ${\rm e}^{-4i\psi_{j}}$, and ${\rm
e}^{0i\psi_{j}}$. Carrying out the second set of integrals in (11) eliminates
the ${\rm e}^{\pm 4i\psi}$ terms, leaving

$\langle\langle\rho\rangle_{\phi}\rangle_{\psi}=\frac{1}{64}\sum_{j}\Bigl{[}\left(U_{j}^{2}F_{1,j}F_{2,j}^{*}+V_{j}^{2}F_{1,j}^{*}F_{2,j}\right)f_{1,j}^{*}f_{2,j}+\left(U_{j}^{2}F_{1,j}^{*}F_{2,j}+V_{j}^{2}F_{1,j}F_{2,j}^{*}\right)f_{1,j}f_{2,j}^{*}\Bigr{]}\,.$
(48)

The next average is easy because all of the $\iota$ dependence is in $U$ and
$V$.

The average over $\iota$ is defined by the third set of integrals in (11).
Using $\langle U_{j}^{2}\rangle_{\iota}=\langle
V_{j}^{2}\rangle_{\iota}=16A_{j}^{2}/5$ from (24) gives

$\langle\langle\langle\rho\rangle_{\phi}\rangle_{\psi}\rangle_{\iota}=\frac{1}{20}\sum_{j}A_{j}^{2}\left(F_{1,j}F_{2,j}^{*}+F_{1,j}^{*}F_{2,j}\right)\left(f_{1,j}^{*}f_{2,j}+f_{1,j}f_{2,j}^{*}\right)\,.$
(49)

The last average to compute is over the source directions ${\bf\Omega}$.

The average over ${\bf\Omega}$ is defined by the final set of integrals in
(11). We assume that the pulsar distances from Earth, and the distances
between the pulsars, are much greater than the wavelength of the GWs. Thus,
$L_{1}\gg 2\pi/\omega$, $L_{2}\gg 2\pi/\omega$, and
$|L_{1}{\hat{\bf{p}}}_{1}-L_{2}{\hat{\bf{p}}}_{2}|\gg 2\pi/\omega$. Examining
(49), one can see from (37) and (42) that $f_{1,j}^{*}f_{2,j}$ and
$f_{1,j}f_{2,j}^{*}$ each give four terms. Three are rapidly oscillating
functions of ${\bf\Omega}_{j}$. When multiplied in (49) by the slowly varying
antenna pattern functions $F$, they average to zero. The fourth terms are the
only ones that survive, and correspond to setting $f_{1,j}^{*}f_{2,j}$ and
$f_{1,j}f_{2,j}^{*}$ to unity in (49). This gives

$\langle\rho\rangle=\frac{1}{10}\sum_{j}A_{j}^{2}\left\langle
F_{1,j}F_{2,j}^{*}+F_{1,j}^{*}F_{2,j}\right\rangle_{{\bf\Omega}}\,.$ (50)

This spherical average of the antenna pattern functions was first evaluated in
1983 by Hellings and Downs [18].

We briefly digress to contrast the current calculation of $\langle\rho\rangle$
to the earlier calculation of $\langle s\rangle$ in Sec. IV.1. There,
following (15), we average the product
$e^{ab}({\bf\Omega})e^{*}_{ab}({\bf\Omega})$ over ${\bf\Omega}$. Here, those
polarization tensors are replaced by antenna pattern functions, which have a
more elaborate dependence upon ${\bf\Omega}$.

The spherical average in (50) is called the HD curve, and is defined by

$\mu_{u}(\gamma)\equiv\langle F_{1,j}F_{2,j}^{*}\rangle_{{\bf\Omega}}=\langle
F_{1,j}^{*}F_{2,j}\rangle_{{\bf\Omega}}\,.$ (51)

Here, $\gamma$ is the angular separation between the two pulsars, so
$\cos\gamma={\hat{\bf{p}}}_{1}\cdot{\hat{\bf{p}}}_{2}$. The function
$\mu_{u}(\gamma)$ is the mean value of the HD correlation for an unpolarized
unit-amplitude source. It is computed explicitly in [9] (see Eq. 3.35 and App.
D) as

$\mu_{\rm
u}(\gamma)=\frac{1}{4}+\frac{1}{12}\cos\gamma+\frac{1}{2}(1-\cos\gamma)\log\left(\frac{1-\cos\gamma}{2}\right)\,.$
(52)

Thus, we obtain the final expression for the mean value, or first moment,

$\langle\rho\rangle=\frac{1}{5}\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}\mu_{u}(\gamma)\,.$
(53)

Not surprisingly, this is proportional to the standard HD correlation
function.

###  V.2  Second moment and variance of $\rho$

We proceed in analogy with the earlier computation of $\langle s^{2}\rangle$
given in Sec. IV.2. We first define

$\displaystyle P_{1,j}$ $\displaystyle\equiv U_{j}{\rm
e}^{2i\psi_{j}}F_{1,j}+V_{j}{\rm e}^{-2i\psi_{j}}F_{1,j}^{*}$ (54)
$\displaystyle P_{2,j}$ $\displaystyle\equiv U_{j}{\rm
e}^{2i\psi_{j}}F_{2,j}+V_{j}{\rm e}^{-2i\psi_{j}}F_{2,j}^{*}\,,$

where the subscript 1 or 2 denotes the pulsar and $j$ labels the source.
Starting with (46), we can then write the correlation of the two pulsars’
redshifts as

$\displaystyle\rho=$
$\displaystyle\frac{1}{64}\sum_{j}\Bigl{[}P_{1,j}P_{2,j}^{*}f_{1,j}^{*}f_{2,j}+P_{1,j}^{*}P_{2,j}f_{1,j}f_{2,j}^{*}\Bigr{]}+$
(55) $\displaystyle\frac{1}{64}\sum_{j\neq
k}\Bigl{[}P_{1,j}P_{2,k}^{*}f_{1,j}^{*}f_{2,k}{\rm
e}^{i(\phi_{j}-\phi_{k})}+P_{1,j}^{*}P_{2,k}f_{1,j}f_{2,k}^{*}{\rm
e}^{-i(\phi_{j}-\phi_{k})}\Bigr{]}\,.$

The first line is the diagonal terms ($j=k$) and the second line is the off-
diagonal ($j\neq k$) ones.

The first step is to average $\rho^{2}$ over the phases $\phi$. The procedure
is the same one used for (17). The quantity $\rho^{2}$ has three terms, but
averaging eliminates the cross term. Then, the average of exponential factors
(18) introduces Kronecker deltas, eliminating one set of double sums. We are
left with

$\displaystyle\langle\rho^{2}\rangle_{\phi}=$
$\displaystyle\frac{1}{4096}\biggl{(}\sum_{j}\Bigl{[}P_{1,j}P_{2,j}^{*}f_{1,j}^{*}f_{2,j}+P_{1,j}^{*}P_{2,j}f_{1,j}f_{2,j}^{*}\Bigr{]}\biggr{)}^{2}+$
(56) $\displaystyle\,\,\,\,\frac{1}{4096}\sum_{j\neq
k}\Bigl{[}P_{1,j}P_{2,k}^{*}P_{1,k}P_{2,j}^{*}\,f_{1,j}^{*}f_{2,k}f_{1,k}^{*}f_{2,j}+P_{1,j}^{*}P_{2,k}P_{1,k}^{*}P_{2,j}\,f_{1,j}f_{2,k}^{*}f_{1,k}f_{2,j}^{*}$
$\displaystyle\qquad\quad\,\,\,\,\,\,+P_{1,j}P_{2,k}^{*}P_{1,j}^{*}P_{2,k}\,f_{1,j}^{*}f_{2,k}f_{1,j}f_{2,k}^{*}+P_{1,k}P_{2,j}^{*}P_{1,k}^{*}P_{2,j}\,f_{1,k}^{*}f_{2,j}f_{1,k}f_{2,j}^{*}\Bigr{]}\,.$

Note that within the double sum, the final two terms are identical, and just
correspond to swapping $j$ and $k$. If we now explicitly square the single-sum
term, we obtain

$\displaystyle\langle\rho^{2}\rangle_{\phi}=$
$\displaystyle\frac{1}{4096}\sum_{j}\Bigl{[}P_{1,j}P_{2,j}^{*}f_{1,j}^{*}f_{2,j}+P_{1,j}^{*}P_{2,j}f_{1,j}f_{2,j}^{*}\Bigr{]}^{2}+$
(57) $\displaystyle\,\,\,\,\frac{1}{2048}\sum_{j\neq
k}\Bigl{[}P_{1,j}P_{2,j}^{*}P_{1,k}P_{2,k}^{*}\,f_{1,j}^{*}f_{2,j}f_{1,k}^{*}f_{2,k}+P_{1,j}^{*}P_{2,j}P_{1,k}^{*}P_{2,k}\,f_{1,j}f_{2,j}^{*}f_{1,k}f_{2,k}^{*}$
$\displaystyle\qquad\quad\,\,\,\,\,\,+P_{1,j}P_{1,j}^{*}P_{2,k}P_{2,k}^{*}\,f_{1,j}^{*}f_{1,j}f_{2,k}^{*}f_{2,k}+P_{1,j}P_{2,j}^{*}P_{1,k}^{*}P_{2,k}\,f_{1,j}^{*}f_{2,j}f_{1,k}f_{2,k}^{*}\Bigr{]}\,.$

Here, two of the four off-diagonal terms that come from the square of the
single-sum term have doubled the terms in the second line of (56). The
remaining two off-diagonal terms provide the final term that appears in (57).
We have also reordered the terms, putting the $j$-factors in front of the
$k$-factors. This concludes the average over $\phi$.

We next average over $\psi$, carrying out the second set of integrals in (11).
Expanding the square in the single sum of (57) gives terms proportional to
$(P_{1,j}P_{2,j}^{*})^{2}$, $(P^{*}_{1,j}P_{2,j})^{2}$ and
$P_{1,j}P_{1,j}^{*}P_{2,j}P_{2,j}^{*}$. From the definition (54), we see that
$P_{j}$ is a linear combination of ${\rm e}^{2i\psi_{j}}$ and ${\rm
e}^{-2i\psi_{j}}$. Thus, only the $\psi$-independent terms survive the
average. From the single sum these are

$\displaystyle\langle(P_{1,j}P_{2,j}^{*})^{2}\rangle_{\psi}$
$\displaystyle=U_{j}^{4}(F_{1,j}F_{2,j}^{*})^{2}+V_{j}^{4}(F_{1,j}^{*}F_{2,j})^{2}+4\,U^{2}_{j}V^{2}_{j}|F_{1,j}|^{2}|F_{2,j}|^{2}$
(58) $\displaystyle\langle(P_{1,j}^{*}P_{2,j})^{2}\rangle_{\psi}$
$\displaystyle=U_{j}^{4}(F_{1,j}^{*}F_{2,j})^{2}+V_{j}^{4}(F_{1,j}F_{2,j}^{*})^{2}+4\,U^{2}_{j}V^{2}_{j}|F_{1,j}|^{2}|F_{2,j}|^{2}$
$\displaystyle 2\langle P_{1,j}P_{1,j}^{*}P_{2,j}P_{2,j}^{*}\rangle_{\psi}$
$\displaystyle=2(U_{j}^{4}+V_{j}^{4}+2U_{j}^{2}V_{j}^{2})|F_{1,j}|^{2}|F_{2,j}|^{2}+2U_{j}^{2}V_{j}^{2}\left[(F_{1,j}F_{2,j}^{*})^{2}+(F_{1,j}^{*}F_{2,j})^{2}\right]\,.$

Since $j\neq k$, the averages over $\psi_{j}$ and $\psi_{k}$ may be carried
out separately in the double sum of (57). In each case, we always have a
product of two $P$ terms. On the same line of how we computed (48), the ${\rm
e}^{\pm 4i\psi}$ terms average to zero. Thus, we obtain

$\displaystyle\langle\langle\rho^{2}\rangle_{\phi}\rangle_{\psi}=$
$\displaystyle\frac{1}{4096}\sum_{j}\biggl{[}\Bigl{(}U_{j}^{4}(F_{1,j}F_{2,j}^{*})^{2}+V_{j}^{4}(F_{1,j}^{*}F_{2,j})^{2}+4\,U_{j}^{2}V_{j}^{2}|F_{1,j}|^{2}|F_{2,j}|^{2}\Bigr{)}\left(f_{1,j}^{*}f_{2,j}\right)^{2}+$
(59)
$\displaystyle\qquad\qquad\,\Bigl{(}U_{j}^{4}(F_{1,j}^{*}F_{2,j})^{2}+V_{j}^{4}(F_{1,j}F_{2,j}^{*})^{2}+4\,U_{j}^{2}V_{j}^{2}|F_{1,j}|^{2}|F_{2,j}|^{2}\Bigr{)}\left(f_{1,j}f_{2,j}^{*}\right)^{2}+$
$\displaystyle\qquad\qquad\,2\Bigl{(}\left(U_{j}^{4}+V_{j}^{4}+2U_{j}^{2}V_{j}^{2}\right)|F_{1,j}|^{2}|F_{2,j}|^{2}+U_{j}^{2}V_{j}^{2}\left((F_{1,j}F_{2,j}^{*})^{2}+(F_{1,j}^{*}F_{2,j})^{2}\right)\Bigr{)}\bigl{|}f_{1,j}\bigr{|}^{2}\bigl{|}f_{2,j}\bigr{|}^{2}\biggr{]}+$
$\displaystyle\frac{1}{2048}\sum_{j\neq
k}\biggl{[}\left(U_{j}^{2}F_{1,j}F_{2,j}^{*}+V_{j}^{2}F_{1,j}^{*}F_{2,j}\right)\left(U_{k}^{2}F_{1,k}F_{2,k}^{*}+V_{k}^{2}F_{1,k}^{*}F_{2,k}\right)f_{1,j}^{*}f_{2,j}f_{1,k}^{*}f_{2,k}+$
$\displaystyle\qquad\qquad\,\left(U_{j}^{2}F_{1,j}^{*}F_{2,j}+V_{j}^{2}F_{1,j}F_{2,j}^{*}\right)\left(U_{k}^{2}F_{1,k}^{*}F_{2,k}+V_{k}^{2}F_{1,k}F_{2,k}^{*}\right)f_{1,j}f_{2,j}^{*}f_{1,k}f_{2,k}^{*}+$
$\displaystyle\qquad\qquad\,\left(U_{j}^{2}+V_{j}^{2})(U_{k}^{2}+V_{k}^{2}\right)|F_{1,j}|^{2}|F_{2,k}|^{2}\bigl{|}f_{1,j}\bigr{|}^{2}\bigl{|}f_{2,k}\bigr{|}^{2}+$
$\displaystyle\qquad\qquad\,\left(U_{j}^{2}F_{1,j}F_{2,j}^{*}+V_{j}^{2}F_{1,j}^{*}F_{2,j}\right)\left(U_{k}^{2}F_{1,k}^{*}F_{2,k}+V_{k}^{2}F_{1,k}F_{2,k}^{*}\right)f_{1,j}^{*}f_{2,j}f_{1,k}f_{2,k}^{*}\biggr{]}$

This concludes the average over $\psi$.

The average over $\iota$ follows immediately from (24), since the $\iota$
dependence enters via $U$ and $V$. Averaging (59) gives

$\displaystyle\langle\langle\langle\rho^{2}\rangle_{\phi}\rangle_{\psi}\rangle_{\iota}=$
(60)
$\displaystyle\frac{1}{5040}\sum_{j}A_{j}^{4}\biggl{[}\Bigl{(}35(F_{1,j}F_{2,j}^{*})^{2}+35(F_{1,j}^{*}F_{2,j})^{2}+2\,|F_{1,j}|^{2}|F_{2,j}|^{2}\Bigr{)}\Bigl{(}\left(f_{1,j}^{*}f_{2,j}\right)^{2}+\left(f_{1,j}f_{2,j}^{*}\right)^{2}\Bigr{)}+$
$\displaystyle\qquad\quad\quad\quad\quad\Bigl{(}(F_{1,j}F_{2,j}^{*})^{2}+(F_{1,j}^{*}F_{2,j})^{2}+142\,|F_{1,j}|^{2}|F_{2,j}|^{2}\Bigr{)}\bigl{|}f_{1,j}\bigr{|}^{2}\bigl{|}f_{2,j}\bigr{|}^{2}\biggr{]}+$
$\displaystyle\frac{1}{200}\sum_{j\neq
k}A_{j}^{2}A_{k}^{2}\biggl{[}\Bigl{(}F_{1,j}F_{2,j}^{*}+F_{1,j}^{*}F_{2,j}\Bigr{)}\Bigl{(}F_{1,k}F_{2,k}^{*}+F_{1,k}^{*}F_{2,k}\Bigr{)}\Bigl{(}f_{1,j}^{*}f_{2,j}f_{1,k}^{*}f_{2,k}+f_{1,j}f_{2,j}^{*}f_{1,k}f_{2,k}^{*}+f_{1,j}^{*}f_{2,j}f_{1,k}f_{2,k}^{*}\Bigr{)}$
$\displaystyle\qquad\quad\quad\quad\quad+4\,|F_{1,j}|^{2}|F_{2,k}|^{2}\bigl{|}f_{1,j}\bigr{|}^{2}\bigl{|}f_{2,k}\bigr{|}^{2}\biggr{]}\,.$

Some care is needed to average (60) over the source directions.

For the final average over ${\bf\Omega}$, we make the same assumptions about
pulsar and interpulsar distances as given before (50). Hence, the situation in
(60) is similar. The antenna pattern functions $F$ vary slowly with
${\bf\Omega}$, whereas the phase difference terms $f$ vary quickly. Thus, the
same reasoning previously used to obtain (50) can also be used here. It
implies that the average of terms which do not contain $|f|^{2}$ are obtained
by setting $f\to 1$. In contrast, the average of terms that contain $|f|^{2}$
are obtained by setting $|f|^{2}\to 1+\chi^{2}$. See Eqs. 3.34 and 3.41 of [9]
for a more detailed account. So, in the end, all that is needed is averages of
the antenna pattern functions $F$ over source direction.

The necessary antenna pattern averages can be obtained from [9], where they
are given in terms of plus and cross polarizations. We use (41) to obtain the
complex polarization equivalents. These first and second moments are

$\displaystyle\bigl{\langle}F_{1,j}F^{*}_{2,j}\bigr{\rangle}_{{\bf\Omega}_{j}}$
$\displaystyle=$ $\displaystyle\mu_{\rm u}(\gamma)$ (61)
$\displaystyle\Re\,\bigl{\langle}\left(F_{1,j}F^{*}_{2,j}\right)^{2}\bigr{\rangle}_{{\bf\Omega}_{j}}$
$\displaystyle=$ $\displaystyle\mu_{\rm u}^{2}(\gamma)+\sigma_{\rm
u}^{2}(\gamma)-\sigma_{\rm p}^{2}(\gamma)$ (62)
$\displaystyle\bigl{\langle}\left|F_{1,j}\right|^{2}\left|F_{2,j}\right|^{2}\bigr{\rangle}_{{\bf\Omega}_{j}}$
$\displaystyle=$ $\displaystyle\mu_{\rm u}^{2}(\gamma)+\sigma_{\rm
u}^{2}(\gamma)+\sigma_{\rm p}^{2}(\gamma)\,.$ (63)

As before, $\gamma$ is the angle between the pulsar directions
${\hat{\bf{p}}}_{1}$ and ${\hat{\bf{p}}}_{2}$. Only the real part is given in
(62), because (60) does not contain the imaginary part. The HD curve function
$\mu_{\rm u}(\gamma)$ was already given in (51), and the variance functions
$\sigma_{\rm u}^{2}(\gamma)$ and $\sigma_{\rm p}^{2}(\gamma)$ are given in
Eqs. D37 and E8 of [9].

We use relations (61)-(63) to carry out the final average of (60) over source
directions ${\bf\Omega}$. In the single-sum term, both (62) and (63) are
required. In the double-sum term, since $j\neq k$, only (61) is needed. Note
that the average of $|F_{1,j}|^{2}$ (or $|F_{2,k}|^{2}$) is obtained from (61)
by setting $\gamma=0$. We obtain

$\displaystyle\langle\rho^{2}\rangle=$
$\displaystyle\frac{1}{1260}\Bigl{[}36\,\mu_{\rm
u}^{2}(\gamma)+36\,\sigma_{\rm u}^{2}(\gamma)-34\,\sigma_{\rm
p}^{2}(\gamma)+\left(1+\chi^{2}\right)^{2}\left(36\,\mu_{\rm
u}^{2}(\gamma)+36\,\sigma_{\rm u}^{2}(\gamma)+35\,\sigma_{\rm
p}^{2}(\gamma)\right)\Big{]}\,\sum_{j}A_{j}^{4}$ (64)
$\displaystyle+\,\frac{1}{50}\Bigr{[}3\,\mu_{\rm
u}^{2}(\gamma)+\left(1+\chi^{2}\right)^{2}\mu_{\rm
u}^{2}(0)\Bigr{]}\,\sum_{j\neq k}A_{j}^{2}A_{k}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{6300}\Bigl{[}-198\,\mu_{\rm
u}^{2}(\gamma)+180\,\sigma_{\rm u}^{2}(\gamma)-170\,\sigma_{\rm
p}^{2}(\gamma)+\left(1+\chi^{2}\right)^{2}\left(180\,\mu_{\rm
u}^{2}(\gamma)+180\,\sigma_{\rm u}^{2}(\gamma)+175\,\sigma_{\rm
p}^{2}(\gamma)-126\,\mu_{\rm u}^{2}(0)\right)\Big{]}\,\sum_{j}A_{j}^{4}$
$\displaystyle+\,\frac{1}{50}\Bigr{[}3\,\mu_{\rm
u}^{2}(\gamma)+\left(1+\chi^{2}\right)^{2}\mu_{\rm
u}^{2}(0)\Bigr{]}\,\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}^{2}\,,$

where the second equality follows from (28). This concludes the calculation of
the second moment of the correlation $\rho$.

The variance $\sigma^{2}=\langle\rho^{2}\rangle-\langle\rho\rangle^{2}$ is
obtained from (64), and from the first moment (53). Thus, for the polarized
ensemble defined by (4), the mean and variance of the correlation between the
pulsars 1 and 2 are

$\displaystyle\langle\rho\rangle_{\rm pol}=$
$\displaystyle\frac{1}{5}\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}\mu_{\rm u}(\gamma)$
(65) $\displaystyle\sigma^{2}_{\rm pol}=$
$\displaystyle\frac{1}{3150}\Bigl{[}261\,\mu_{\rm
u}^{2}(\gamma)+450\,\sigma_{\rm u}^{2}(\gamma)+265\,\sigma_{\rm
p}^{2}(\gamma)-252\,\mu_{\rm
u}^{2}(0)\Bigr{]}\,\sum_{j}A_{j}^{4}+\frac{1}{50}\Bigr{[}\mu_{\rm
u}^{2}(\gamma)+4\,\mu_{\rm
u}^{2}(0)\Bigr{]}\,\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}^{2}\,,$

where the sums are over the $N$ sources and we set $\chi=1$, which gives the
correct physical description of a PTA. To see the effects of discarding the
pulsar terms, set $\chi=0$, as discussed after (32).

###  V.3  Comparison of polarized and unpolarized ensembles

We now compare these results for the polarized ensemble to the corresponding
results for the unpolarized ensemble.

The mean and variance of the correlation $\rho$ for the unpolarized ensemble
(1) are given in Eqs. 3.36 and 3.48 of [9]. Setting $\chi=1$ gives

$\displaystyle\langle\rho\rangle_{\rm unpol}$
$\displaystyle=\frac{1}{2}\Bigl{(}\sum_{j}{\cal{A}}_{j}^{2}\Bigr{)}\mu_{\rm
u}(\gamma)$ (66) $\displaystyle\sigma^{2}_{\rm unpol}$
$\displaystyle=\frac{1}{8}\Bigl{(}2\,\mu_{\rm
u}^{2}(\gamma)+5\,\sigma^{2}_{\rm u}(\gamma)+3\,\sigma^{2}_{\rm
p}(\gamma)-4\,\mu_{\rm
u}^{2}(0)\Bigr{)}\,\sum_{j}{\cal{A}}_{j}^{4}+\frac{1}{8}\Bigl{(}\mu_{\rm
u}^{2}(\gamma)+4\,\mu_{\rm
u}^{2}(0)\Bigr{)}\Bigl{(}\sum_{j}{\cal{A}}_{j}^{2}\Bigr{)}^{2}\,.$

To obtain this, we have used the relation

$\sigma_{\rm c}^{2}(\gamma)=\sigma^{2}_{\rm p}(\gamma)+\sigma^{2}_{\rm
u}(\gamma)+\mu_{\rm u}^{2}(\gamma),$ (67)

which is given in Eq. F4 of [9].

As in Sec. IV.3, to compare the two ensembles on “equal footing” we assume
that the amplitudes of the individual GW sources are related by (31). With
this, we find

$\displaystyle\langle\rho\rangle_{\rm pol}$
$\displaystyle=\langle\rho\rangle_{\rm unpol}$ (68)
$\displaystyle\sigma^{2}_{\rm pol}(\gamma)$ $\displaystyle=\sigma^{2}_{\rm
unpol}(\gamma)+\left(\frac{3}{70}\mu_{\rm
u}^{2}(\gamma)+\frac{3}{70}\sigma_{\rm
u}^{2}(\gamma)+\frac{38}{1575}\sigma_{\rm
p}^{2}(\gamma)\right)\sum_{j}A_{j}^{4}.$

The additional terms that appear on the rhs are _positive_. This follows from
general principles: introducing additional degrees of freedom (here
$\iota_{j}$ and $\psi_{j}$) increases the variance.

There are two limits of interest. First, for a single polarized source, this
result reduces to that given in Eq. B11 and footnote 43 of [9]. Second, for a
very large number of polarized sources, the
$\left(\sum_{j}A_{j}^{2}\right)^{2}$ terms dominate the $\sum_{j}A_{j}^{4}$
terms, hence in this limit $\sigma^{2}_{\rm pol}/\sigma^{2}_{\rm unpol}\to 1$.

##  VI Cosmic variance

The variance $\sigma_{\rm pol}^{2}$ in (65) quantifies the _total_ range of
variation in the correlation between two pulsars. This variance is the sum of
_pulsar variance_ and _cosmic variance_.

For a given set of GW sources, different pairs of pulsars with the same
angular separation will exhibit different correlations. This variation,
associated with picking different pulsar pairs separated by the same angle, is
the pulsar variance. In contrast, the cosmic variance arises because any given
set of GW sources differs from the mean.

The pulsar variance can be removed from the total variance by averaging the
correlation $\rho$ over all pairs of pulsars separated by the same angle
$\gamma$. This replaces the correlation $\rho$ with a _pulsar-averaged_
correlation $\Gamma(\gamma)$, where $\Gamma(\gamma)$ is the average of $\rho$
over all pairs of pulsars separated by the same angle $\gamma$. (By
definition) the variance of $\Gamma(\gamma)$ is the cosmic variance [9].

We thus define the pulsar-averaged correlation as

$\Gamma(\gamma)\equiv\langle\rho\rangle_{\rm p}\,.$ (69)

Here, the subscript “p” denotes an average over all pulsar pairs
${\hat{\bf{p}}}_{1}$ and ${\hat{\bf{p}}}_{2}$ separated by $\gamma$, where
${\hat{\bf{p}}}_{1}$ and ${\hat{\bf{p}}}_{2}$ are uniformly distributed over
the sphere. The only terms in the correlation $\rho$ (46) which depend upon
the pulsar positions are the antenna pattern functions $F_{j}$ (41) and the
phase difference terms $f_{j}$ (42).

The products of $f$ reduce to unity when computing the pulsar average. This is
because of our assumption that the pulsars 1 and 2 are separated by many
gravitational wavelengths. Each $f$ includes a rapidly oscillating complex
exponential function of the pulsar direction ${\hat{\bf{p}}}$ (42). The
product of this function with the slowly varying $F$ averages to zero. Thus,
provided that pulsars 1 and 2 are separated by many gravitational wavelengths,
only the Earth terms survive. Since the Earth terms are unity, this implies

$\langle F_{1,j}^{*}F_{2,k}\,f_{1,j}f_{2,k}^{*}\rangle_{\rm p}=\langle
F_{1,j}^{*}F_{2,k}\rangle_{\rm p}\,.$ (70)

The same reasoning applies to products where fewer or more of the terms are
complex conjugated.

The only terms left to average to obtain $\Gamma(\gamma)$ from (46) are
products of antenna pattern functions. Those calculations are carried out in
Appendix G of [9] and the results, given in Eq. 4.2 of the same reference, are

$\displaystyle\langle F_{1,j}^{*}F_{2,k}\rangle_{\rm p}=\langle
F_{1,j}F_{2,k}^{*}\rangle_{\rm p}=$ $\displaystyle\mu(\gamma,\beta_{jk})$ (71)
$\displaystyle\langle F_{1,j}F_{2,k}\rangle_{\rm p}=\langle
F_{1,j}^{*}F_{2,k}^{*}\rangle_{\rm p}=$
$\displaystyle\mu(\pi-\gamma,\pi-\beta_{jk})\,.$

Here, $\gamma$ is the angular separation between the directions to the two
pulsars, and $\beta_{jk}$ is the angular separation between the direction to
the two sources $j$ and $k$:
$\cos\beta_{jk}={\bf\Omega}_{j}\cdot{\bf\Omega}_{k}$. The two-point function
$\mu(\gamma,\beta)$ is real and is a generalization of the HD curve $\mu_{\rm
u}(\gamma)$. The two functions coincide in the case where $\beta_{jk}$
vanishes: $\mu(\gamma,0)=\mu_{\rm u}(\gamma)$.

We now calculate the pulsar-averaged correlation $\Gamma(\gamma)$. Using (70)
and (71) to average (46) over all pulsar pairs gives

$\displaystyle\Gamma(\gamma)=\frac{1}{64}\sum_{j,k}\Bigl{[}\Bigl{(}$
$\displaystyle\bigl{(}U_{j}V_{k}{\rm e}^{2i(\psi_{j}+\psi_{k})}+V_{j}U_{k}{\rm
e}^{-2i(\psi_{j}+\psi_{k})}\bigr{)}{\rm e}^{i(\phi_{j}-\phi_{k})}+$ (72)
$\displaystyle\bigl{(}V_{j}U_{k}{\rm e}^{2i(\psi_{j}+\psi_{k})}+U_{j}V_{k}{\rm
e}^{-2i(\psi_{j}+\psi_{k})}\bigr{)}{\rm
e}^{-i(\phi_{j}-\phi_{k})}\Bigr{)}\mu(\pi-\gamma,\pi-\beta_{jk})+$
$\displaystyle\Bigl{(}$ $\displaystyle\bigl{(}U_{j}U_{k}{\rm
e}^{2i(\psi_{j}-\psi_{k})}+V_{j}V_{k}{\rm
e}^{-2i(\psi_{j}-\psi_{k})}\bigr{)}{\rm e}^{i(\phi_{j}-\phi_{k})}+$
$\displaystyle\bigl{(}V_{j}V_{k}{\rm e}^{2i(\psi_{j}-\psi_{k})}+U_{j}U_{k}{\rm
e}^{-2i(\psi_{j}-\psi_{k})}\bigr{)}{\rm
e}^{-i(\phi_{j}-\phi_{k})}\Bigr{)}\mu(\gamma,\beta_{jk})\Bigr{]}\,.$

We now compute the mean and variance of $\Gamma$.

###  VI.1 First moment of $\Gamma$

Calculating the ensemble average of $\Gamma(\gamma)$ is very similar to the
calculation of $\langle\rho\rangle$ in Sec. V.1. We proceed by averaging (72)
over the sources ensemble as described in (11).

We start by calculating the average over the GW source phases $\phi_{j}$,
evaluating the first set of integrals in (11). This gives

$\displaystyle\langle\Gamma(\gamma)\rangle_{\phi}$
$\displaystyle=\frac{1}{32}\sum_{j}\Bigl{[}\left(U_{j}^{2}+V_{j}^{2}\right)\mu(\gamma,\beta_{jj})+U_{j}V_{j}\left({\rm
e}^{4i\psi_{j}}+{\rm
e}^{-4i\psi_{j}}\right)\mu(\pi-\gamma,\pi-\beta_{jj})\Bigr{]}$ (73)
$\displaystyle=\frac{1}{32}\sum_{j}\left(U_{j}^{2}+V_{j}^{2}\right)\mu_{\rm
u}(\gamma)\,,$

where only the diagonal terms ($j=k$) survived the average. Note that since
$\beta_{jj}=0$, the two-point function $\mu(\gamma,\beta_{jj})=\mu_{\rm
u}(\gamma)$ and $\mu(\pi-\gamma,\pi-\beta_{jj})=\mu(\pi-\gamma,\pi)=0$. (For
the final equality, see for example Fig. 12 in Appendix G of [9].) As
expected, the result is proportional to the HD function $\mu_{\rm u}(\gamma)$.

Since the final line of (73) is independent of $\psi_{j}$, the only averages
remaining in (11) are over $\iota_{j}$ and ${\bf\Omega}_{j}$. The quantities
that appear in the final line of (73) are independent of ${\bf\Omega}_{j}$,
because $U_{j}$ and $V_{j}$ only depend upon $\iota_{j}$. Using (24) for the
averaging gives:

$\langle\Gamma(\gamma)\rangle=\frac{1}{5}\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}\mu_{\rm
u}(\gamma)\,.$ (74)

Thus, the expectation value of the pulsar-averaged correlation
$\Gamma(\gamma)$ in (74) is the same as the expectation value of $\rho$ in
(53). This is not surprising: _on average_ an observer would always obtain a
correlation curve proportional to the HD function.

###  VI.2 Second moment and variance of $\Gamma$

We now compute the second moment and the variance of $\Gamma(\gamma)$. The
calculation is very similar to the one carried out in Sec. V.2 for $\rho$, and
corresponds again to evaluating the average given by (11).

We begin by computing $\Gamma^{2}(\gamma)$. As we did for $\rho$ in (55), the
sum in (72) is first divided into diagonal terms ($j=k$) and the off-diagonal
ones ($j\neq k$). $\Gamma^{2}(\gamma)$ then consists of three terms: (i) the
square of the diagonal sum, (ii) the square of the off-diagonal sum, and (iii)
the cross-product terms.

We then average $\Gamma^{2}(\gamma)$ over the phases $\phi_{j}$. Two of the
three terms are easy to average: (i) does not change because it is independent
of $\phi_{j}$, and (iii) vanishes. The square of the off-diagonal sum (ii)
contains the product of two double sums: $\sum_{j\neq k}\sum_{\ell\neq m}$.
This double sum has sixteen terms, each one containing a complex exponential
${\rm e}^{\pm i((\phi_{j}-\phi_{k})\pm(\phi_{\ell}-\phi_{m}))}$. Using (18),
(ii) reduces to a sum over $j\neq k$, and the sixteen terms reduce to four. We
obtain

$\displaystyle\langle\Gamma^{2}(\gamma)\rangle_{\phi}=\frac{1}{1024}\Bigl{(}\sum_{j}\left(U_{j}^{2}+V_{j}^{2}\right)\mu_{\rm
u}(\gamma)\Bigr{)}^{2}+\frac{1}{1024}\sum_{j\neq k}\biggl{[}$ (75)
$\displaystyle\left(U_{j}V_{k}{\rm e}^{2i(\psi_{j}+\psi_{k})}+V_{j}U_{k}{\rm
e}^{-2i(\psi_{j}+\psi_{k})}\right)\left(V_{j}U_{k}{\rm
e}^{2i(\psi_{j}+\psi_{k})}+U_{j}V_{k}{\rm
e}^{-2i(\psi_{j}+\psi_{k})}\right)\mu^{2}(\pi-\gamma,\pi-\beta_{jk})+$
$\displaystyle\left(U_{j}V_{k}{\rm e}^{2i(\psi_{j}+\psi_{k})}+V_{j}U_{k}{\rm
e}^{-2i(\psi_{j}+\psi_{k})}\right)\left(V_{j}V_{k}{\rm
e}^{2i(\psi_{j}-\psi_{k})}+U_{j}U_{k}{\rm
e}^{-2i(\psi_{j}-\psi_{k})}\right)\mu(\pi-\gamma,\pi-\beta_{jk})\mu(\gamma,\beta_{jk})+$
$\displaystyle\left(U_{j}U_{k}{\rm e}^{2i(\psi_{j}-\psi_{k})}+V_{j}V_{k}{\rm
e}^{-2i(\psi_{j}-\psi_{k})}\right)\left(V_{j}U_{k}{\rm
e}^{2i(\psi_{j}+\psi_{k})}+U_{j}V_{k}{\rm
e}^{-2i(\psi_{j}+\psi_{k})}\right)\mu(\pi-\gamma,\pi-\beta_{jk})\mu(\gamma,\beta_{jk})+$
$\displaystyle\left(U_{j}U_{k}{\rm e}^{2i(\psi_{j}-\psi_{k})}+V_{j}V_{k}{\rm
e}^{-2i(\psi_{j}-\psi_{k})}\right)\left(V_{j}V_{k}{\rm
e}^{2i(\psi_{j}-\psi_{k})}+U_{j}U_{k}{\rm
e}^{-2i(\psi_{j}-\psi_{k})}\right)\mu^{2}(\gamma,\beta_{jk})\biggr{]}\,.$

This concludes the average over $\phi_{j}$.

The next step is the average over $\psi_{j}$. Carrying out the second set of
integrals in (11), the second and third terms of the double sum of (75)
vanish. We obtain

$\displaystyle\langle\langle\Gamma^{2}(\gamma)\rangle_{\phi}\rangle_{\psi}=\frac{1}{1024}\Biggl{\\{}$
$\displaystyle\sum_{j\neq
k}\Bigl{[}\left(U_{j}^{2}+V_{j}^{2}\right)\left(U_{k}^{2}+V_{k}^{2}\right)\mu_{\rm
u}^{2}(\gamma)+2U_{j}^{2}V_{k}^{2}\mu^{2}(\pi-\gamma,\pi-\beta_{jk})+\left(U_{j}^{2}U_{k}^{2}+V_{j}^{2}V_{k}^{2}\right)\mu^{2}(\gamma,\beta_{jk})\Bigr{]}$
(76) $\displaystyle+$ $\displaystyle\sum_{j}(U_{j}^{2}+V_{j}^{2})^{2}\mu_{\rm
u}^{2}(\gamma)\Biggr{\\}}\,,$

where we also expanded the square of the single sum in the first line of (75)
and exploited the symmetry of the double sum under the interchange of $j$ and
$k$.

We now average over $\iota_{j}$. The only terms depending upon $\iota_{j}$ are
the functions $U_{j}$ and $V_{j}$. Using (24), we obtain

$\langle\langle\langle\Gamma^{2}(\gamma)\rangle_{\phi}\rangle_{\psi}\rangle_{\iota}=\frac{71}{1260}\mu_{\rm
u}^{2}(\gamma)\sum_{j}A_{j}^{4}+\frac{1}{50}\sum_{j\neq k}\Bigl{(}2\mu_{\rm
u}^{2}(\gamma)+\mu^{2}(\pi-\gamma,\pi-\beta_{jk})+\mu^{2}(\gamma,\beta_{jk})\Bigr{)}A_{j}^{2}A_{k}^{2}\,.$
(77)

This concludes the average over $\iota_{j}$.

The remaining step to evaluate the second moment is to average over the source
directions ${\bf\Omega}_{j}$. The dependence on ${\bf\Omega}_{j}$ is via
$\beta_{jk}$, where $\cos(\beta_{jk})={\bf\Omega}_{j}\cdot{\bf\Omega}_{k}$.
The spherical averages over ${\bf\Omega}_{j}$ and ${\bf\Omega}_{k}$ are
equivalent to evaluating the average over $\cos\beta_{jk}$. To visualize this,
assume that the $j$’th source is aligned with the $\hat{z}$ axis. Thus,
$\beta_{jk}=\beta_{\hat{z}k}\equiv\theta$, where $\theta$ is the traditional
polar angle in spherical polar coordinates. Then, for any function
$g(\beta_{jk})$,

$\int\frac{{\rm d}{\bf\Omega}_{j}}{4\pi}\int\frac{{\rm
d}{\bf\Omega}_{k}}{4\pi}g(\beta_{jk})=\int\frac{{\rm
d}{\bf\Omega}_{k}}{4\pi}g(\beta_{\hat{z}k})=\frac{1}{4\pi}\int_{0}^{2\pi}{\rm
d}\phi\int_{-1}^{1}\\!\\!{\rm
d}(\cos\theta)\,g(\theta)=\frac{1}{2}\int_{-1}^{1}\\!\\!g(\theta)\,{\rm
d}(\cos\theta)\,.$ (78)

We can now use (78) to average (77) over the source directions.

The only terms depending upon $\beta_{jk}$ are
$\mu^{2}(\pi-\gamma,\pi-\beta_{jk})$ and $\mu^{2}(\gamma,\beta_{jk})$. Those
averages are discussed in detail and computed in Appendix G of [9]. Starting
with the second line of Eq. G12 of [9] and making use of (78) we have

$\displaystyle\int\frac{{\rm d}{\bf\Omega}_{j}}{4\pi}\int\frac{{\rm
d}{\bf\Omega}_{k}}{4\pi}\Bigl{(}\mu^{2}(\gamma,\beta_{jk})+\mu^{2}(\pi-\gamma,\pi-\beta_{jk})\Bigr{)}=\frac{1}{2}\int_{-1}^{1}\Bigl{(}\mu^{2}(\gamma,\beta)+\mu^{2}(\pi-\gamma,\pi-\beta)\Bigr{)}{\rm
d}(\cos\beta)=2\tilde{\mu}^{2}(\gamma)\,,$ (79)

where the complete expression for $\tilde{\mu}^{2}(\gamma)$ is given in Eq.
G11 of [9]. Using (79) to average (77) over the source directions then gives

$\displaystyle\langle\Gamma^{2}(\gamma)\rangle$
$\displaystyle=\frac{71}{1260}\mu_{\rm
u}^{2}(\gamma)\sum_{j}A_{j}^{4}+\frac{1}{25}\Bigl{(}\mu^{2}_{\rm
u}(\gamma)+\tilde{\mu}^{2}(\gamma)\Bigr{)}\sum_{j\neq k}A_{j}^{2}A_{k}^{2}$
(80) $\displaystyle=\Bigl{(}\frac{103}{6300}\mu_{\rm
u}^{2}(\gamma)-\frac{1}{25}\tilde{\mu}^{2}(\gamma)\Bigr{)}\sum_{j}A_{j}^{4}+\frac{1}{25}\Bigl{(}\mu^{2}_{\rm
u}(\gamma)+\tilde{\mu}^{2}(\gamma)\Bigr{)}\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}^{2}\,,$

where we used (28) to obtain the second line. This completes the calculation
of the second moment of $\Gamma(\gamma)$.

We can finally write the variance of the pulsar-averaged correlation
$\Gamma(\gamma)$. From (74) and (80), we obtain

$\displaystyle\sigma^{2}_{\begin{subarray}{c}\text{cosmic}\\\ \text{pol
\quad}\end{subarray}}$
$\displaystyle\equiv\langle\Gamma^{2}(\gamma)\rangle-\langle\Gamma(\gamma)\rangle^{2}$
(81) $\displaystyle=\Bigl{(}\frac{103}{6300}\mu_{\rm
u}^{2}(\gamma)-\frac{1}{25}\tilde{\mu}^{2}(\gamma)\Bigr{)}\sum_{j}A_{j}^{4}+\frac{1}{25}\tilde{\mu}^{2}(\gamma)\Bigl{(}\sum_{j}A_{j}^{2}\Bigr{)}^{2}\,,$

which corresponds to the _cosmic variance_ of the correlation $\rho$ (46) for
the polarized ensemble defined by (4).

###  VI.3  Comparison of polarized and unpolarized ensembles

We now compare the cosmic variance (81) of the polarized ensemble (4) with
that obtained for the unpolarized ensemble (1).

The cosmic variance for the unpolarized ensemble (1) is given by Eq. 4.8 of
[9] as

$\sigma^{2}_{\begin{subarray}{c}\text{cosmic}\\\
\text{unpol\,}\end{subarray}}=\frac{1}{4}\tilde{\mu}^{2}(\gamma)\Bigl{(}\Bigl{(}\sum_{j}{\cal{A}}_{j}^{2}\Bigr{)}^{2}-\sum_{j}{\cal{A}}_{j}^{4}\Bigr{)}\,,$
(82)

where $\tilde{\mu}^{2}(\gamma)$ is the same function as the one introduced in
(79).

As in Secs. IV.3 and V.3, to compare the two ensembles on “equal footing”, we
assume that the amplitudes of the individual GW sources are related by (31).
Thus, we obtain

$\sigma^{2}_{\begin{subarray}{c}\text{cosmic}\\\ \text{pol
\quad}\end{subarray}}=\sigma^{2}_{\begin{subarray}{c}\text{cosmic}\\\
\text{unpol\,}\end{subarray}}+\frac{103}{6300}\mu_{\rm
u}^{2}(\gamma)\sum_{j}A_{j}^{4}\,.$ (83)

The additional term of the polarized ensemble variance is _positive_. Thus,
the polarized ensemble has a larger cosmic variance than the unpolarized one.
This is to be expected: see the discussion in Sec. V.3. We note that the
additional degrees of freedom have increased the variance “in quadrature” by a
term proportional to the square of the mean. There is probably a simple
physical explanation for this.

In the many-source limit $N\to\infty$, the ratio of the two cosmic variances
approaches unity,
$\sigma^{2}_{\text{cos,pol}}/\sigma^{2}_{\text{cos,unpol}}\to 1$. The same
behavior was also found for the total variance in Sec. V.3. This follows
immediately from (83), because the sum $\sum_{j}A_{j}^{2}$ diverges as the
number of sources grows, whereas the sum $\sum_{j}A_{j}^{4}$ converges to
$A^{4}\zeta(4/3)$. Thus, in the many-source limit, the two ensembles have the
same behavior.

##  VII Comparison with a Gaussian ensemble

The frequency-domain Gaussian ensemble (introduced in [6, 7]) is often used to
describe PTA sources. In fact, the enterprise software package [5] used to
analyze PTA data assumes that the sources are described by a Gaussian
ensemble.

Here, we compare our results for an ensemble of $N$ discrete sources to the
corresponding results for a Gaussian ensemble. Specifically, we compare the
HD-correlation mean, total variance and cosmic variance of our ensemble of
polarized sources to those of the Gaussian ensemble. These are derived in
Appendix C of [9].

In the weak-signal limit, the gravitational wave amplitude can always be
written as a plane wave expansion

$h_{ab}(t,{\bf x})=\sum_{A}\int{\rm d}f\int{\rm d}{\bf\Omega}\,{\rm e}^{2\pi
if(t-{\bf\Omega}\cdot{\bf x})}h_{A}(f,{\bf\Omega})e^{A}_{ab}({\bf\Omega})\,.$
(84)

Here, $A$ denotes the polarization ($+,\times$), $f\in\Re$ denotes the GW
frequency, ${\bf\Omega}$ denotes a unit vector in the GW propagation
direction, and $e^{A}_{ab}({\bf\Omega})$ is the polarization tensor defined in
Eqs. D6 and D7 of [9]. The complex function $h_{A}(f,{\bf\Omega})$ is the GW
Fourier amplitude, and differs from one realization to the next. For any
realization, the GW amplitude $h_{ab}(t,{\bf x})$ is real, implying that
$h_{A}(f,{\bf\Omega})=h^{*}_{A}(-f,{\bf\Omega})$, where $*$ denotes the
complex conjugate.

Under the assumption that the GWs are generated by some central-limit-theorem
process, we describe their statistics by treating the functions
$h_{A}(f,{\bf\Omega})$ as Gaussian random variables. In this Section, we use
angle brackets $\langle Q\rangle$ to denote the average of $Q$ in this
Gaussian ensemble. Note that this is different than earlier in this paper,
where the ensemble average was defined by (11). In this Section, the ensemble
average refers to a Gaussian ensemble of functions $h_{A}(f,{\bf\Omega})$ (see
Eq. C2 of [9]).

The Gaussian ensemble is completely defined by:

  * (a)

The first moment: the ensemble average of the strain vanishes:

$\langle h_{A}(f,{\bf\Omega})\rangle=0\,.$ (85)

  * (b)

The second moment:

$\langle
h_{A}(f,{\bf\Omega})h^{*}_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})\rangle=\delta_{AA^{\prime}}\delta^{2}({\bf\Omega},{\bf\Omega}^{\prime})\delta(f-f^{\prime})H(f)\,,$
(86)

where $H(f)=H(-f)$ is a real function that defines the measure of the squared
amplitude of perturbations at GW frequency f (see Eqs. C5 and C7 of [9]). The
delta function in frequency implies that the ensemble is second-order
stationary in time, while the delta function on the sphere ensures second-
order stationarity in space. Lastly, the delta function in polarization
implies that the two polarizations are identical, but uncorrelated (the
ensemble is unpolarized).

  * (c)

Isserlis’ theorem [19]: the higher moments (expected value of the product of
$h_{A}(f,{\bf\Omega})h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})\cdots
h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})$) are given
by sums of expected values of all possible terms containing one or two $h$’s.
For example

$\displaystyle\langle
h_{A}(f,{\bf\Omega})h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})\rangle=$
$\displaystyle\,\langle h_{A}(f,{\bf\Omega})\rangle\langle
h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})\rangle+$
(87) $\displaystyle\,\langle
h_{A}(f,{\bf\Omega})h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})\rangle\langle
h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})\rangle+$
$\displaystyle\,\langle
h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})\rangle\langle
h_{A}(f,{\bf\Omega})h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})\rangle=$
$\displaystyle=$ $\displaystyle\,0\,,$

which vanishes as direct consequence of condition (a), and

$\displaystyle\langle
h_{A}(f,{\bf\Omega})h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})h_{A^{\prime\prime\prime}}(f^{\prime\prime\prime},{\bf\Omega}^{\prime\prime\prime})\rangle=$
$\displaystyle\,\langle
h_{A}(f,{\bf\Omega})h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})\rangle\langle
h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})h_{A^{\prime\prime\prime}}(f^{\prime\prime\prime},{\bf\Omega}^{\prime\prime\prime})\rangle+$
(88) $\displaystyle\,\langle
h_{A}(f,{\bf\Omega})h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})\rangle\langle
h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})h_{A^{\prime\prime\prime}}(f^{\prime\prime\prime},{\bf\Omega}^{\prime\prime\prime})\rangle+$
$\displaystyle\,\langle
h_{A}(f,{\bf\Omega})h_{A^{\prime\prime\prime}}(f^{\prime\prime\prime},{\bf\Omega}^{\prime\prime\prime})\rangle\langle
h_{A^{\prime}}(f^{\prime},{\bf\Omega}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},{\bf\Omega}^{\prime\prime})\rangle\,.$

The ensemble average of the products of more than four $h$’s follows the same
pattern, reducing to sums of the expected values of all possible terms
containing one or two $h$’s.

Following the same prescription as the ones reported in Sec. V and VI,
Appendix C of [9] computes the total variance and cosmic variance for a
Gaussian ensemble of GW sources. We compare the results shown in that Appendix
with the ones estimated for our polarized discrete ensemble of sources.

###  VII.1 Correlation mean and total variance: comparison of polarized and
Gaussian ensembles

The mean and total variance of the HD correlation for a Gaussian ensemble
(denoted with subscript “G”) are

$\displaystyle\langle\rho\rangle_{\rm G}$ $\displaystyle=h^{2}\mu_{\rm
u}(\gamma)$ (89) $\displaystyle\sigma^{2}_{\rm G}$
$\displaystyle={\mathcal{h}}^{4}\bigl{(}\mu_{\rm u}^{2}(\gamma)+4\,\mu_{\rm
u}^{2}(0)\bigr{)}\,,$

which correspond to Eqs. C18 and C28 of [9]. The scale factors $h$ and
${\mathcal{h}}$ are measures of the strain, expressed as integrals of the
spectral function $H(f)$. From Eqs. C19 and C26 of [9]:

$\displaystyle h^{2}$ $\displaystyle\equiv 4\pi\int H(f){\rm d}f$ (90)
$\displaystyle{\mathcal{h}}^{4}$ $\displaystyle\equiv(4\pi)^{2}\int{\rm
d}f\int{\rm
d}f^{\prime}{{\rm\,sinc}}^{2}\bigl{(}\pi(f-f^{\prime})T\bigr{)}H(f)H(f^{\prime})\,,$

where $T$ is the total observation time.

We now compare the mean correlation and variance estimated for the Gaussian
ensemble (89) with the results for our polarized ensemble of discrete sources
(65). To compare the two ensembles “on equal footing”, we set the
normalizations (GW-amplitude scale factors) to:

$\displaystyle h^{2}$ $\displaystyle\equiv\frac{1}{5}\sum_{j}A_{j}^{2}$ (91)
$\displaystyle{\mathcal{h}}^{4}$
$\displaystyle\equiv\frac{1}{50}\big{(}\sum_{j}A_{j}^{2}\big{)}^{2}\,,$

which implies that ${\mathcal{h}}^{4}=h^{4}/2$.

With this normalization choice (91), the expected mean correlation and its
variance for the ensemble of polarized sources (65) converge to the Gaussian
ones (89) in the limit of an infinite number of sources. A similar convergence
was found in [9] for an ensemble of unpolarized GW sources.

###  VII.2 Cosmic variance: comparison of polarized and Gaussian ensembles

We now compare the cosmic variance for a polarized and a Gaussian ensemble of
GW sources.

For a polarized ensemble, the expected value of the pulsar-averaged
correlation $\Gamma$ and its variance are given in (74) and (81). For a
Gaussian ensemble of sources, the corresponding quantities are given in
Appendix C of [9]:

$\displaystyle\langle\Gamma(\gamma)\rangle_{\rm G}$
$\displaystyle=h^{2}\mu_{\rm u}(\gamma)$ (92) $\displaystyle\sigma^{2}_{\rm
G,cosmic}$ $\displaystyle=2{\mathcal{h}}^{4}\tilde{\mu}^{2}(\gamma)\,,$

from Eqs. C42 and C45 of that reference.

In the limit of an infinite number of sources and with the normalization
relationship (91), the cosmic variance of the ensemble of polarized sources
(81) converges to that of the Gaussian ensemble of sources (92). A similar
convergence was found in [9] for an ensemble of unpolarized GW sources.

##  VIII Conclusion

The most common model of the gravitational wave stochastic background is the
Gaussian ensemble. However, previous work [9] shows that this does not have
the same statistical properties as the gravitational wave background produced
by a finite set of discrete sources. To study the discrete-source case, that
work constructs an ensemble which is simple, but unrealistic. It is
unrealistic because the discrete sources correspond to circular-orbit binary
systems, viewed face-on. These produce circularly polarized gravitational
waves (1). Physically, we expect the orbits to have all possible orientations
ranging from face-on to face-off, producing elliptically polarized
gravitational waves.

This paper adopts a more realistic model. We employ an ensemble of discrete
sources (4) described in Sec. II, for which the orbital angular momentum
directions are uniformly distributed over the two-sphere of directions. This
is a more realistic gravitational wave ensemble, because it corresponds to a
set of circular-orbit binary systems with randomly distributed orbital planes.
For example, a set of supermassive black-hole binary sources whose circular
orbits have inclination angles $\iota$ which are random variables with
$\cos\iota$ distributed uniformly between $-1$ and $1$. The previous model
assumes $\cos\iota=1$ (face-on).

To construct the ensemble, in Sec. III we use a set of $N$ distant
gravitational waves sources, labeled by an integer $j=1,\dots,N$. Their
amplitudes fall off with distance as would be expected for identical
(intrinsic amplitude) sources, uniformly distributed inside a sphere centered
at Earth. Each source has its own orbital phase $\phi_{j}$, sky direction
$-{\bf\Omega}_{j}$, and orbital orientation parameters: $\cos\iota_{j}$ for
inclination and $\psi_{j}$ for elliptical axis direction. The ensemble is
defined by treating these $5N$ parameters (note: ${\bf\Omega}$ counts as two)
as independent random variables; the ensemble average (11) is a product of
integrals over them. The discrete ensemble of [9] has only $3N$ independent
random variables, because $\iota$ and $\psi$ have fixed values $\iota=\psi=0$.

The mean and variance of the HD correlation are complicated to calculate. So,
before evaluating these, Sec. IV examines a simpler quantity $s$: the time-
averaged squared strain. The mean and variance of $s$ are computed, and then
compared for the polarized and unpolarized ensembles in Sec. IV.3. If the
parameters of the two ensembles are set to give the same mean value of $s$,
then the variance of $s$ is larger for the polarized ensemble. This variance
is “cosmic variance”: it describes how $s$ varies between different
realizations of the universe.

In Sec. V, using identical techniques, we calculate the two quantities most
relevant for pulsar timing arrays. These are the mean and variance of the HD
correlation $\rho(\gamma)$, as a function of the angle $\gamma$ between the
lines of sight to two pulsars. Later, in Sec. VI, we compute the cosmic
variance by first averaging the correlation over many pulsar pairs separated
by angle $\gamma$. As with $s$, if the parameters of the ensembles are picked
so that the mean values of $\rho(\gamma)$ agree, then the polarized ensemble
has larger variance and larger cosmic variance than the unpolarized ensemble.
This is not surprising. The varying orbital inclinations and polarizations add
additional degrees of freedom, which increase the variance between different
realizations drawn from the polarized ensemble.

We also examine the means and variances in the limit of an infinite density of
GW sources with the time-averaged squared strain $s$ at Earth held fixed. In
all cases, these approach the corresponding quantities for the Gaussian
ensemble. This raises an interesting question for future work: is this true
for _all_ observable quantities? Can we quantify how close the ensembles of
discrete sources are to the Gaussian ensemble, and prove that in some limit
they approach the Gaussian ensemble?

In the future, galaxy surveys will provide sky direction, distance, mass and
frequency information for likely GW sources. The methods employed here can be
used to quantify how closely the reconstruction of the HD correlation should
match the (mean value prediction of the) HD curve. This work could also be
extended to take into account the probability distributions of chirp mass and
frequency for the sources, as well as the specific pulsar locations on the
sky, as in [16].

## References

  * Antoniadis  _et al._ [2023] J. Antoniadis _et al._ (EPTA and InPTA Collaborations),  _The second data release from the European Pulsar Timing Array: III. Search for gravitational wave signals_ , Astronomy $\&$ Astrophysics 678, A50 (2023).
  * Agazie  _et al._ [2023] G. Agazie _et al._ (NANOGrav Collaboration),  _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , The Astrophysical Journal Letters 951, L8 (2023).
  * Reardon  _et al._ [2023] D. J. Reardon _et al._ (PPTA Collaboration),  _Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array_ , The Astrophysical Journal Letters 951, L6 (2023).
  * Xu  _et al._ [2023] H. Xu _et al._ (CPTA Collaboration),  _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).
  * Ellis  _et al._ [2020] J. Ellis, M. Vallisneri, S. Taylor, and P. Baker,  _ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE_ , Zenodo 10.5281/zenodo.4059815 (2020).
  * Allen [1997] B. Allen, The stochastic gravity-wave background: sources and detection, in  _Relativistic Gravitation and Gravitational Radiation_ , Cambridge Contemporary Astrophysics, Vol. 26, edited by J.-A. Marck and J.-P. Lasota (Cambridge University Press, 1997) pp. 373–417, arXiv:gr-qc/9604033 [gr-qc] .
  * Allen and Romano [1999] B. Allen and J. D. Romano,  _Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities_ , Phys. Rev. D 59, 102001 (1999).
  * Romano and Cornish [2017] J. D. Romano and N. J. Cornish,  _Detection methods for stochastic gravitational-wave backgrounds: a unified treatment_ , Living Rev. Relativ. 20, 2 (2017).
  * Allen [2023] B. Allen,  _Variance of the Hellings-Downs correlation_ , Phys. Rev. D 107, 043018 (2023).
  * Sesana  _et al._ [2008] A. Sesana, A. Vecchio, and C. N. Colacino,  _The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays_ , Monthly Notices of the Royal Astronomical Society 390, 192 (2008).
  * Conneely  _et al._ [2019] C. Conneely, A. H. Jaffe, and C. M. F. Mingarelli,  _On the amplitude and Stokes parameters of a stochastic gravitational-wave background_ , Monthly Notices of the Royal Astronomical Society 487, 562 (2019).
  * Note [1] Since the Earth to SSB distance is about 500 seconds, and the period of the GWs is orders of magnitude longer than this, the time at the Earth and SSB may be treated as equivalent for this analysis.
  * Peters and Mathews [1963] P. C. Peters and J. Mathews,  _Gravitational Radiation from Point Masses in a Keplerian Orbit_ , Phys. Rev. 131, 435 (1963).
  * Peters [1964] P. C. Peters,  _Gravitational Radiation and the Motion of Two Point Masses_ , Phys. Rev. 136, B1224 (1964).
  * Timpano  _et al._ [2006] S. E. Timpano, L. J. Rubbo, and N. J. Cornish,  _Characterizing the galactic gravitational wave background with LISA_ , Phys. Rev. D 73, 122001 (2006).
  * Allen and Romano [2023] B. Allen and J. D. Romano,  _Hellings and Downs correlation of an arbitrary set of pulsars_ , Phys. Rev. D 108, 043026 (2023).
  * Note [2] The physical range over which $\psi$ is averaged is $[0,\pi/2]$. However, since all of the relevant functions are periodic on this interval, averaging over a larger interval, as we do in (11), gives the same result.
  * Hellings and Downs [1983] R. W. Hellings and G. S. Downs,  _Upper limits on the isotropic gravitational radiation background from pulsar timing analysis._ , Astrophys. J. 265, L39 (1983).
  * 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 (1918), note that Isserlis’ theorem is often called “Wick’s theorem” in the physics community, although the Wick’s work was three decades later.