Pulsar Timing Array Harmonic Analysis and Source Angular Correlations

Gravitational waves (GWs) influence the arrival times of radio signals coming from pulsars. Here, we investigate the harmonic space approach to describing a pulsar's response to GWs. We derive and discuss the"diagonalized form"of the response, which is a sum of spin-2-weighted spherical harmonics of the GW direction multiplied by normal (spin-weight 0) spherical harmonics of the pulsar direction. We show how this allows many useful objects, for example, the Hellings and Downs two-point function, to be easily calculated. The approach also provides a clear description of the gauge dependence. We then employ this harmonic approach to model the effects of angular correlations in the sky locations of GW sources (sometimes called"statistical isotropy"). To do this, we construct rotationally invariant ensembles made up of many Gaussian subensembles, each of which breaks rotational invariance. Using harmonic techniques, we compute the cosmic covariance and the total covariance of the Hellings and Downs correlation in these models. The results may be used to assess the impact of angular source correlations on the Hellings and Downs correlation, and for optimal reconstruction of the Hellings and Downs curve in models where GW sources have correlated sky locations.


I. INTRODUCTION
There is a considerable literature on the topic of pulsar timing arrays (PTAs), which may be on the verge of making five-sigma detections of nHz gravitational waves (GWs) [1][2][3][4].PTAs rely on the effect that GWs have on shifting the arrival times of radio pulses.An introductory discussion of how they work can be found in [5].
As shown in [5] [Sec.2.2, Eq. ( 14)], the shift in arrival times is due to the Sachs-Wolfe effect [6], which also creates temperature fluctuations of the (electromagnetic) cosmic background radiation (CBR).So it is not surprising that there is a considerable literature which applies tools and techniques drawn from CBR to PTAs [7][8][9][10].Here, we use the term "harmonic analysis" for this approach, which describes the response of PTA pulsars in terms of spherical harmonic functions on the two-dimensional sphere.
This paper presents the most important of these tools and results from a physical perspective, and illustrates how they may be used to describe the response of PTAs.While much of this can be found in the corresponding specialist literature [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], and the seminal paper on the topic [8] is a decade old, we hope to offer some fresh insights as well as a few new results.A brief outline of the paper, including links to key equations, follows.
Our analysis assumes that pulsars are perfect clocks: the response of a pulsar to a GW can be described as a "redshift" (or "blueshift") of the clock frequency.In Sec.II, we review the way in which pulsar redshift Z responds to the GW amplitude at Earth.The response FðΩ; Ω p Þ is a function (2.2) of the direction Ω of GW propagation and of the direction Ω p to the pulsar.Here, Ω and Ω p are (coordinates of) points on the unit two-sphere, and Ω and Ωp are the corresponding unit vectors.
The function F is complex, with real and imaginary parts that describe, respectively, the response to plus-and crosspolarized GWs.While the magnitude (2.7) of this response depends only upon the angle between the direction of GW propagation and the direction to the pulsar in the sky, the phase of the response has a more complicated dependence on these two positions.We write F in diagonal form (2.10) as a sum of spin-weight 0 harmonic functions of Ω p and spin-weight 2 harmonic functions of Ω.This reflects the spin-two nature of GWs, correctly incorporating both the magnitude and the phase of the response.
The remainder of the paper exploits this diagonal form.In Sec.III, we use it to derive the harmonic-space form of the Hellings and Downs (HD) curve (3.3), by averaging over source directions Ω, which was the original approach employed in 1983 by Hellings and Downs [26].In Sec.IV, we derive a simple formula that can be used to "pulsar average" any function QðΩ p ; Ω q Þ of the sky directions to pulsars p and q, as illustrated in Fig. 5 of [27].This produces a function QðγÞ of a single variable γ, which is the angle between the directions to two pulsars (4.4).In Sec.V, we again use the harmonic decomposition of F to derive the HD curve (5.1), but this time as the pulsar average of the correlation for a single GW source.This equality between the source direction average for a single pulsar pair and the pulsar average for a single source direction was first demonstrated, highlighted, and discussed in [28].In comparison with the original approach employed by Hellings and Downs [26], we believe that this is a "better way" to define and to think about the HD curve.
In Sec.VI, we turn attention to the HD two-point function μðγ; Ω; Ω 0 Þ, which is defined by (6.1) and was first computed in [27] (Appendix G).This is the average correlation of a pair of pulsars separated by angle γ, for GWs from sources radiating in directions Ω and Ω 0 .We compute it using the pulsar averaging recipe from Sec. IV, showing that the magnitude of this quantity depends only upon γ and upon the angle β between the two sources (6.8).We obtain a beautiful new harmonic form (6.9) for the twopoint function, as a sum of products of Legendre and Jacobi polynomials, which is used later in the paper to study the cosmic covariance and variance.
Starting in Sec.VII, we employ Gaussian statistical ensembles of GW sources.As shown by the central limit theorem, these describe GWs which are produced by the incoherent sum of many weak sources.Working in a circular polarization basis, we first compute the cosmic covariance for the standard ensemble of unpolarized sources [fully defined by first and second moments (7.4)] demonstrating explicitly how the phase of μðγ; Ω; Ω 0 Þ cancels out.Then, in Sec.VIII, we exploit the harmonic form of the two-point function found earlier, to obtain explicit harmonic decompositions of the cosmic covariance (8.4) and variance (8.5).
In Sec.IX, we consider cosmological ensembles where the source sky locations have nonzero angular covariance (they are not a Poisson process [29,30] in the limit of an infinite density of vanishingly weak sources.)This ensemble is sometimes called "statistically isotropic" in contrast with the standard Gaussian ensemble, which is called "purely isotropic."(These names are misleading: while both ensembles are isotropic, a typical representative universe drawn from either ensemble breaks rotational invariance, and thus is not isotropic.) In this paper, a statistical ensemble is a (possibly infinite) collection of representative universes containing GWs.Each of these representative universes is called a "realization": it contains a specific set of GW sources, with specific properties (frequencies, amplitudes, …) at specific sky locations.An ensemble is rotationally invariant if and only if, for every realization, it contains all rotated (about the Earth/origin) versions of that realization.
Our approach is to build statistical ensembles composed of Gaussian subensembles (but note that the resulting full ensemble is not Gaussian [31]).Each individual Gaussian subensemble has preferred directions, defined by a function ψðΩÞ associated with that subensemble, which breaks rotational invariance.Nevertheless, the full ensemble maintains rotational invariance, because for any Gaussian subensemble that it contains with a given ψ, it contains all other Gaussian subensembles described by rotated versions of ψ.This construction is detailed in Sec.IX A. It has been used before, for example in [32][33][34], but without this explicit description [35].
This approach enables the study of GW source models where the (sky) locations of the sources have nontrivial angular covariance.Such correlations arise for any ensemble constructed from a finite number of GW sources at discrete sky locations [36].This explicit construction provides a sound basis for similar "statistical isotropy" calculations which appear in the literature [32][33][34] but whose justification is problematic [37].The correlation function CðΩ; Ω 0 Þ þ 1 ¼ hψðΩÞψðΩ 0 Þi ψ that describes correlations among GW source locations is an average over all subensembles (9.3), and only depends upon the dot product Ω • Ω0 of the directions to the two sources.The coefficients C L of the Legendre-polynomial decomposition (9.4) of Cð Ω • Ω0 Þ characterize the type/degree of the angular correlations.
In Sec.IX B, we compute the cosmic variance and covariance for this correlated-in-angle ensemble.(It would be logical to begin with the total variance/covariance, but that is more complicated, so we do it after.)The averages within a given Gaussian subensemble lead to intermediate results such as (9.6) and (9.8) that are exactly as for the standard case, except that they contain factors of ψðΩÞ.Averaging over the full ensemble then introduces the function C. By employing the harmonic decomposition of the two-point function, a simple result for the cosmic covariance (9.13) is obtained.
In Sec.IX C, we return to the calculation of the total covariance.Again, we compute the mean (9.21) and correlation (9.22) of the HD correlation within a given Gaussian subensemble.Again, these are identical to standard results, apart from containing factors of ψðΩÞ.In Sec.IX D, we then carry out the averages over subensembles, to obtain the covariance C pq;rs (9.29) of the HD correlation for the full ensemble.The latter is determined by a function D pq;rs (9.27) of four pulsar directions Ω p , Ω q , Ω r , and Ω s .(Others have also investigated this quantity, see Appendix E of [8] and citations therein.)We derive a rotationally invariant form for D, given by (9.37) and (9.38).This could be used for optimal reconstructions of the HD correlation curve [38] that take account of correlations among source sky positions.
In Sec.IX E, we compute the total variance for the ensemble of correlated sources, which is the diagonal part C pq;pq of the covariance.From symmetry, this only depends upon the angle γ between the directions to p and q, so it is unaffected by pulsar averaging.By carrying out a pulsar average, the total variance can be expressed as a sum of Legendre polynomials (9.43),where the coefficients for a given C L are explicitly given in (9.47).This is followed by a brief conclusion.
The appendixes contain technical details.In Appendix A we provide key formulas for spin-weighted spherical harmonics.In Appendix B we derive the diagonal form of the HD response, and in Appendix C we compute two-point functions for the four different combinations of linear polarizations (C6).

II. RESPONSE OF A PULSAR TO A GRAVITATIONAL WAVE
It is common to talk about physical effects in terms of fields and particles, for example, the influence of an electric field on an electron.PTAs can be described in similar terms: the influence of GWs on the arrival time of pulses [5,39].The pulsars may be thought of as ideal clocks, and the influence of GWs is to reduce or increase their tick rates.This can be quantified as a time-dependent redshift or blueshift, which is the time derivative of the timing residual.These same quantities, redshift and blueshift, are also used to describe temperature fluctuations in the CBR.
Note that for both PTAs and CBR, the use of "red" and "blue" is a historical misnomer.The clock frequency for PTAs is hundreds of Hz, and the CBR consists of infrared radiation.Both frequencies are well below the visible part of the spectrum.
We begin with the response of a pulsar to a GW traveling in direction Ω, where this symbol indicates a pair of angles θ; ϕ in usual spherical polar coordinates.The corresponding unit-length vector from the origin (Earth) is denoted Note that, as seen from Earth, the sky direction to the very distant GW source is − Ω.
The sky coordinates of the pulsar are Ω p .A unit vector from Earth to the pulsar has the same components as in (2.1) but with coordinates θ p ; ϕ p .Note that we typically label or index pulsars with the subscripts p, q, r, and s.
The (redshift) response of the pulsar consists of an "Earth term" and a "pulsar term."Until the distance to pulsars is known to about light-year precision, the latter cannot be measured, so in this paper we mostly ignore it [40].The Earth term (for a unit-amplitude, circularly polarized GW) is where m and n are a pair of unit length vectors which are (a) perpendicular to the GW direction Ω and (b) perpendicular to each other: This choice of m and n is important for what follows, but note that it is inconsistent with much of the literature.
To describe the GW, we could have made different choices for m and n, rotating them in the m-n plane through an angle which is an arbitrary function of Ω ¼ θ; ϕ.The choice we have made "fixes the gauge," ensuring that FðΩ; Ω p Þ is only a function of θ, θ p , and ϕ − ϕ p .For other choices of gauge, FðΩ; The two GW polarizations are usually labeled "þ" and "×," corresponding to polarization tensors In what follows, we adopt the Einstein summation convention: if a coordinate index appears twice in a given term, then for that term the repeated index should be summed over the three coordinates.Some examples:

Ωb
p e × ab ðΩÞ.Hence, the real part of F is the redshift produced by a plus-polarized GW with unit amplitude at Earth at that moment in time, and the imaginary part of F is the redshift produced by a cross-polarized GW: Thus, F is the instantaneous redshift response of a pulsar to a circularly polarized GW of unit strain amplitude at Earth, where "circularly polarized" means a polarization tensor e þ ab þ ie × ab .It is easy to see that the modulus of the response jFðΩ; Ω p Þj only depends upon the angle between Ω and Ωp .The modulus of the numerator of (2.2) is ð2:6Þ The third equality follows since Ω, m, and n form an orthonormal basis, so δ ab ¼ ma mb þ na nb þ Ωa Ωb is the 3 × 3 identity matrix or the Kronecker delta.The fourth equality is because Ωp has unit length.From (2.2) and (2.6), the modulus of F is Thus, the modulus of F is completely determined by the angle between the GW direction and the pulsar direction.However, the phase of F, which depends upon the polarization of the GW, is not a function of this angle alone.Since GWs at Earth are very weak, to obtain the instantaneous redshift at Earth for a given pulsar, we simply add up the real parts of hðtÞF Ã for each source, where the real and imaginary parts of h are the amplitudes at Earth of the two different polarizations.(If we were using timing residuals rather than redshift to describe the effect of GWs, then an integral over time would be needed rather than simply the instantaneous product.) The antenna pattern function F can be thought of as a propagator or response function which encodes the way that pulsar redshift responds to a GW.Physically, it is enough to specify this response for one source direction (say, Ω ¼ ẑ) and all pulsar directions.Then, the response for any other source and pulsar directions can be obtained by rotation.This embodies a fundamental tenet of the principle of relativity: physical observables are coordinateindependent.
While this is true, there is an important subtlety.To explain it, let us start with the pulsar response for a gravitational wave propagating in the positive z direction, obtained by setting θ ¼ 0 in (2.2).For that case Ωp • ð m þ i nÞ ¼ sin θ p e iðϕ p −ϕÞ , so by inspection one obtains − cos θ p Þe 2iðϕ p −ϕÞ : ð2:8Þ Note that this has a strange feature.Although θ ¼ 0 places the GW source direction at the North Pole for any value of ϕ, the response still depends upon ϕ.The reason has to do with the behavior of the polarization vectors m and n.If we let the GW propagation direction Ω approach the North Pole along different lines of longitude, the limiting values for the polarization vectors m and n depend upon which line of longitude is followed [41].Hence, the response still depends upon ϕ.
At the root of this odd behavior is the following observation.Since (2.8) gives the response of a pulsar at any point on the sky to a GW propagating in the z direction, we should be able to determine the response of a pulsar in any direction, to a GW propagating in any other direction, simply by rotating the z axis to the desired new propagation direction.But the rotation must not only carry the ẑ vector to the new GW propagation direction: it must also carry the pair of vectors m and n to the correct ones at a different point on the sphere.To say it in another way, the response (2.2) only depends upon the dot products of different vectors, which are rotation-invariant.But, if the GW source is carried to a new sky position, then the corresponding vectors m and n must also be carried along in a way that matches their definitions in (2.3).If not, then F rotates by a complex phase, so is not invariant.See the paragraph following (B5) for a precise statement.
There is a simple formula which encodes this complicated invariance in a beautiful way.If we first define a set of numerical coefficients by This expression is extremely useful.In this paper, it plays a central role, similar to that of the spherical harmonic decomposition of the Green's function in electrostatics.
The relationship (2.10) is derived in Appendix B and is a mathematical equality: for any choice of the four arguments θ; ϕ; θ p ; ϕ p , the right-hand sides of (2.2) and (2.10) yield the same complex number.Similar expressions can be found in the literature, for example, [14] [Eq.(18)] or [23] [Eq.(39)], though the latter has an undetermined complex phase.
For convenience, we define A 0 ¼ A 1 ¼ 0, so that sums like the one in (2.10) can be written P lm .It is also helpful to define coefficients ð2:11Þ These simplify the appearance of equations which follow.The functions Y lm ðΩ p Þ which appear on the rhs of (2.10) are the familiar spherical harmonics.These govern the way that the response varies with pulsar direction.In contrast, the 2 Y lm ðΩÞ, through which the GW direction Ω enters the equation, are spin-2 weighted spherical harmonics.These spin-2 weighted harmonics form a complete orthonormal set on the unit sphere and have properties similar to normal spherical harmonics, for example, their ϕ-dependence is e imϕ .While only their general properties are needed for this paper, we give a precise definition in (A5), and full details may be found in [8] (Appendix A).
Note that for all spherical harmonics, we use the sign, phase and normalization conventions of [8] (Appendixes A and B), where a complete set of formulas is given; the most important ones are reproduced in Appendix A, and further details may be found in [42].
The representation of FðΩ; Ω p Þ given in (2.10) is a "factored" or "diagonal" form.In contrast, suppose that we tried to express F in terms of ordinary spherical harmonics.Since it is a square integrable function of Ω and of Ω p , it can be decomposed as a sum of the form FðΩ; Ω p Þ ¼ P lm P l 0 m 0 a lm;l 0 m 0 Y lm ðΩÞY Ã l 0 m 0 ðΩ p Þ for some set of expansion coefficients a lm;l 0 m 0 .As discussed immediately after (2.3), FðΩ; Ω p Þ as defined by (2.2) only depends upon ϕ and ϕ p through the difference ϕ − ϕ p .Thus, the expansion coefficients a lm;l 0 m 0 vanish for m ≠ m 0 .But, unlike the expansion in (2.10), the coefficients a lm;l 0 m 0 are not diagonal in l and l 0 : a lm;l 0 m ≠ 0 for l ≠ l 0 .
The individual plus-and cross-polarization components are easily extracted.Either from (2.2), or from (2.10) and (A3), one can immediately see that Here, the overlines indicate antipodal points on the sphere: This is also how the individual polarization components were extracted in [27].The correlation between pulsars p and q is a function ϱ pq ðΩÞ of their directions, and of the propagation direction Ω of the GW source.This is often called the "HD integrand" and can be written in several equivalent forms: The second equality shows that the real part of FF Ã may be obtained from FF Ã by swapping the locations of the two pulsars.

III. THE HELLINGS AND DOWNS CURVE AS AN AVERAGE OVER SOURCE DIRECTIONS
The HD curve μ u ðγÞ was originally defined [26] as the correlation between two pulsars p and q separated by angle γ in the sky, uniformly averaged over source directions, for a unit amplitude unpolarized source.We use to denote the integral over the unit two-sphere.To average over directions, an additional factor of 1=4π must be included.The angle between the pulsars is where the final equality follows immediately from (2.1).The computation of the HD curve, starting from (2.10), is trivial.We denote the sky locations of the two pulsars by Ω p and Ω q , and let ϱ pq ðΩÞ given in (2.13) denote their correlation.The average of this is a l P l ðcos γÞ; ð3:3Þ where P l ðzÞ is the Legendre polynomial of order l, and γ ¼ cos −1 ð Ωp • Ωq Þ is the angle between the lines of sight to the two pulsars.The second equality follows directly from the definition (2.13), the third from (2.10), the fourth and fifth equalities follow because the spin-2 weighted harmonics form an orthonormal set on the unit sphere, the sixth equality follows from the addition theorem for spherical harmonics, and the final equality follows from the definitions of A l and a l in (2.9) and (2.11).
The final expression in (3.3) is the standard harmonicspace form of the famous HD curve μ u ðγÞ.One can easily carry out the sum [8] (Sec.III.E) to obtain the position space form Note that the same result would have been obtained without taking the real part on the second line of (3.3).This is because the imaginary part of FðΩ; Ω p ÞF Ã ðΩ; Ω q Þ is odd under Ω → Ω, so it integrates to zero.
To reduce clutter, we often omit the summation limits on l and m.In such cases, l ¼ 0; 1; 2; … and m ¼ −l; −l þ 1; …; l − 1; l.Here, the sum is effectively over l ¼ 2; 3; …, because a l and A l vanish for l < 2.

IV. PULSAR AVERAGING
"Pulsar averaging" is a useful calculational method, which was first introduced in [28] and then developed further in [27,38].It is defined as follows.Given a function QðΩ p ; Ω q Þ which depends upon the position of two pulsars p and q, the pulsar average of Q is a function of angle γ ∈ ½0; π, and is defined by where δðxÞ is the ordinary Dirac delta function.Replacing Q by the constant function Q ¼ 1, one can easily verify that this average is correctly normalized, meaning that h1i pq ∈ γ ¼ 1.This definition corresponds to an average over all unit vectors Ωp uniformly distributed on the sphere, and all unit vectors Ωq uniformly distributed in a cone at angle γ around Ωp , as illustrated in Fig. 5 of [27].It has a close analog in experimental practice, for example, when the Hellings and Downs curve is "reconstructed" by binning together measured correlations from large numbers of pulsar pairs [38] with similar separation angles.
For calculational purposes, it is helpful to express the Dirac delta function in (4.1) in terms of spherical harmonics.To do this, begin with the Dirac delta function expressed as a sum of Legendre polynomials P l ðxÞ, as derived in Eq. (4.20) of [38].On the interval x; In (4.2), set x ¼ Ωp • Ωq and x 0 ¼ cos γ on the lhs, and on the rhs replace P l ð Ωp • Ωq Þ using the addition theorem (3.4).This gives , where the lhs is the two-dimensional delta function on the unit two-sphere S 2 .
If we return to the definition (4.1) of the pulsar average, and replace the delta function with (4.3), we obtain This recipe for computing the pulsar average of any function QðΩ p ; Ω q Þ of pulsar positions will be used later for computing the total variance in models with correlated GW source sky locations.
The pulsar average of the function QðΩ p ; Ω q Þ ¼ Y lm ðΩ p ÞY Ã l 0 m 0 ðΩ q Þ will be needed later.This is evaluated starting from the definition (4.1) as X l 00 m 00 P l 00 ðcosγÞδ ll 00 δ mm 00 δ l 0 l 00 δ m 0 m 00 ¼ 1 4π δ ll 0 δ mm 0 P l ðcosγÞ: ð4:5Þ Here, the second equality is obtained using (4.3), the third equality follows from the orthonormality of the spherical harmonics, and the final equality from the definition of the Kronecker delta.We now use this to carry out some additional harmonic-space computations.

V. THE HELLINGS AND DOWNS CURVE AS A PULSAR AVERAGE FOR A SINGLE GW SOURCE
An alternative definition of the HD curve is as the pulsar average of the cross-correlation (2.13) for one fixed circularly polarized GW point source.This approach was first investigated in [28] and then further developed in [27].Here, we compute this pulsar average, starting from the harmonic expansion (2.10) of the response function F.
For this computation, we fix Ω, and compute the pulsar average of the correlation (2.13) The second equality is obtained by substituting the diagonal form (2.10) for F, the third by substituting the pulsar average of two spherical harmonics given by (4.5), the fourth from the sum of spin-2 weighted harmonics the fifth equality follows from the definition of a l in (2.11) and the final equality from comparison with the average over source directions computed in (3.3).This equality, between (a) the pulsar average for a single source and (b) the average response of a single pair of pulsars to an isotropically distributed set of (noninterfering) sources, was first demonstrated in [28].
In the next section, we will discuss the addition theorem for spin-2 weighted spherical harmonics, from which (5.2) may be obtained as a special case by setting β ¼ χ ¼ 0 in (6.4).

VI. THE HELLINGS AND DOWNS TWO-POINT FUNCTION
Previous work [27,38,43,44] on the variance of the HD correlation exploited a two-point function.Here, this is defined in analogy with Eq. (G1) of [27] as This is averaging the complex redshift response of a pulsar with sky direction Ωp to a distant unit-amplitude GW point source with sky direction − Ω, with the corresponding response for a second pulsar Ωq to a second unit-amplitude point source with sky direction − Ω0 .[The minus signs are explained after (2.1).]As before, γ is the angular separation on the sky of the two pulsars.
The original definition given in [27] is slightly different: it is a real quantity μðγ; βÞ whose square is the squared modulus of the complex quantity μðγ; Ω; Ω 0 Þ defined here.We will see that the modulus depends upon the directions to the two GW point sources only via the angle β ∈ ½0; π between their lines of sight, where The magnitude μ 2 ðγ; βÞ ¼ jμðγ; Ω; Ω 0 Þj 2 is what matters: the phase of μðγ; Ω; Ω 0 Þ is a "gauge artifact" that drops out of observable quantities.
To evaluate the two-point function (6.1), we substitute F from (2.10) into the definition and use (4.5) to compute the pulsar average of Y lm ðΩ p ÞY Ã l 0 m 0 ðΩ q Þ.We obtain The final sum over m is the spin-2 equivalent of the traditional addition theorem (3.4) for scalar harmonics.
The addition theorem for spin-weighted harmonics is given in [8] [(A9)-(A11)].Using [8] [(A6)] and the relation between the Wigner "big D" and "small d" matrices D j m 0 m ðϕ; θ; ψÞ ¼ e −im 0 ϕ d j m 0 m ðθÞe −imψ , the sum appearing in (6.3) may be written where we have expressed the Wigner small d matrix in terms of Jacobi polynomials.These are polynomials in sin 2 ðβ=2Þ and cos 2 ðβ=2Þ, and are illustrated in Fig. 1.
In contrast with the corresponding sum of scalar harmonics, the sum on the final line of (6.3) does not just depend upon the angular separation β between Ω and Ω 0 .While the magnitude of (6.4) is only a function of β, its phase has a complicated dependence upon the positions of the two GW sources.This dependence is via the real angle χ defined by In the notation of [8], The angle χðΩ; Ω 0 Þ ∈ ½−π; π may be defined by inverting (6.5), with arctan in the range ½−π=2; π=2 or in the range ½0; π.Alternatively, χ may be defined in the range ½0; 4π as the argument of the complex number whose imaginary and real parts are (respectively) the numerator and denominator in (6.5).Because χ only enters (6.4) via e 2iχ , these different choices are equivalent.
An important property of χ is that it is an antisymmetric function of its two arguments: This proves that χ cannot be written as a function of β, since β is a symmetric function of Ω and Ω 0 .Another important property of χ, which also follows directly from its definition (6.5), is that χ changes sign if both arguments are sent to their antipodal points.Using the notation introduced in (2.12), this is written But χ is mostly a nuisance: as discussed in [27] (Appendix G), we will see that e 2iχ is a gauge artifact that drops out of physically observable quantities.Making use of the addition theorem for spin-weighted spherical harmonics provides an elegant harmonic decomposition of the two-point function.Substituting the sum over m in (6.4)Note that in these equations, the real quantity μðγ; βÞ may have either sign, so it cannot be interpreted as a radius in the complex plane.As it must, the two-point function (6.8) reduces to the normal HD curve in the limit of coincident GW sources Ω 0 → Ω, where β → 0 and χ → 0. Since the Jacobi polynomials are normalized to P ðα;βÞ l ð1Þ ¼ 1, (6.8) and (6.9) immediately give μðγ; Ω; ΩÞ ¼ μðγ; 0Þ ¼ μ u ðγÞ, in agreement with the HD curve of (3.3).

VII. COSMIC VARIANCE AND COVARIANCE AND GAUGE INDEPENDENCE
We now investigate the complex phase e 2iχðΩ;Ω 0 Þ which appears in (6.8), and show that it drops out of the cosmic covariance, which is a physical observable.This also establishes the gauge independence of the cosmic variance, since it is the covariance restricted to the diagonal.While the results are more general, here we demonstrate them for the specific case of the Gaussian ensemble [38].
The GW metric perturbations in any representative universe may be defined via a plane-wave expansion [ where the reader should keep in mind that h þ ðf; ΩÞ ¼ h Ã þ ð−f; ΩÞ and h × ðf; ΩÞ ¼ h Ã × ð−f; ΩÞ are complex quantities.An ensemble is defined by a set of Fourier amplitude functions hðf; ΩÞ.Each specific function corresponds to a particular universe within the ensemble.
An ensemble may equivalently be defined by specifying all moments of hðf; ΩÞ.Letting angle brackets hi denote FIG. 1.The Jacobi polynomials of (6.4) are shown for l ¼ 2; …; 5.
averages over that ensemble, the Gaussian ensemble is fully defined by the first and second moments taken together with Isserlis's theorem [45].Here, HðfÞ ¼ Hð−fÞ is a real spectral function, and the factor of two is to maintain notational consistency with [27] [Eq.(C4)] and [38].Isserlis's theorem defines the higher-order moments hhðf; ΩÞhðf 0 ; Ω 0 Þ…hðf 00 ; Ω 00 Þi, where any of the functions might also be complex-conjugated, in terms of the first and second moments given by (7.4).
The relations in (7.4) are usually stated for plus-and cross-polarization components.At first glance it appears that (7.4) provides two second moments, whereas the conventional expressions have only one.That is misleading: the conventional expressions have four second moments, for the four combinations of plus and cross, two of which vanish.If we had used a right-and leftcircular polarization basis, then the last two lines of (7.4) could be combined into a single equation, with a Kronecker delta for the two polarization states on the rhs.The apparent extra factor of two arises because the hhh Ã i term is the sum of the linear polarization plus-plus and crosscross terms.
We define the correlation ρ pq between pulsars p and q following Eq.(C15) of [27].For any representative universe in the ensemble, the correlation between pulsars is ρ pq ≡ Z p ðtÞZ q ðtÞ; ð7:5Þ where Z p ðtÞ is the (real, physical) redshift of pulsar p as a function of time, and overline denotes a time average.In what follows, the averaging-time interval is denoted T, which may equivalently be taken as the total observation time.
The pulsar-averaged correlation ΓðγÞ is defined following Eq.(C41) of [27], as ΓðγÞ ≡ hρ pq i pq ∈ γ : ð7:6Þ Here, the angle brackets denote the average over all pulsar pairs p and q separated by angle γ on the sky, as defined in Sec.IV.The pulsar-averaged correlation in any representative universe in the ensemble may be computed in the same way as [27] [Eq.(C41)].It is where γ ≡ π − γ is the angular sky separation between Ω p and Ωq or between Ω q and Ωp , and sincx ≡ sinðxÞ=x.The first equality follows from the definition of ΓðγÞ as the pulsar-averaged correlation at angle γ, with the factors of 1=2 arising from taking the real part as given in (7.1).(As shown in [27], only Earth terms survive the pulsar averaging, so pulsar terms have been dropped.)The second equality follows from the definition (6.1) of the two-point function μðγ; Ω; Ω 0 Þ and the use of to generate complex conjugates of F, as previously employed in (2.12).
The ensemble average of Γ may be computed from inspection of (7.7), using the second moments (7.4) for the Gaussian ensemble.The second and third terms vanish, and the first and fourth terms give To obtain the final equality, we have used the fact that μðγ; Ω; ΩÞ ¼ μ u ðγÞ is the HD curve, and independent of source direction Ω.The squared GW strain at Earth is defined using notation compatible with [27] and [38].
To compute the covariance and variance, we need the deviation of the correlation away from the mean, for any representative of the ensemble.This is ΔΓðγÞ ≡ ΓðγÞ − hΓðγÞi; ð7:10Þ where, as before, angle brackets with no trailing subscript denote an ensemble average.It follows immediately from the definition above that hΔΓðγÞi vanishes.The cosmic covariance is the ensemble average Note that the cosmic covariance may have either sign, whereas the cosmic variance (the value of the covariance along the diagonal γ ¼ γ 0 ) is non-negative.Notationally, they are easily distinguished, because the cosmic variance has one argument, whereas the cosmic covariance has two.The cosmic covariance σ 2 cos ðγ; γ 0 Þ can be computed directly from (7.7).The expression for ΓðγÞΓðγ 0 Þ contains 16 terms, whose average over the Gaussian ensemble can be evaluated using Isserlis's theorem.This implies that hðf;ΩÞh Ã ðf 0 ;Ω 0 Þhðf 00 ; Ω 00 Þh Ã ðf 000 ; and that the ensemble average of terms with unequal numbers of h and h Ã vanish.
To evaluate the cosmic covariance σ 2 cos ðγ; γ 0 Þ, we start by noting that among the 16 terms of ΓðγÞΓðγ 0 Þ are 10 terms containing unequal numbers of h and h Ã ; their ensemble averages vanish.Each of the remaining six terms contains two delta functions in frequency and two delta functions on the sphere.Integrating those out gives where γ0 ≡ π − γ 0 and we have defined (see Appendixes A and B of [38]) The key point is that the complex phase expð2iχÞ cancels out of the cosmic covariance, as can be seen by inspection of (7.13).For example, the first two terms are The first equality is obtained using the harmonic form (6.8) for the two-point function, and the second equality follows from the antisymmetry of χ under interchange of the arguments (6.6) or antipodal reflection (6.7).Carrying out similar phase cancellations for the remaining terms in (7.13) yields a simple expression for the cosmic covariance.If we let which should be compared with [27] [Eq.(G12)] and is identical to [38] [Eq.(4.32)].

VIII. HARMONIC FORM OF THE COSMIC VARIANCE AND COVARIANCE
Starting from the harmonic decomposition (6.9) of the two-point function, it is straightforward to obtain a harmonic form for the cosmic variance and covariance.
The Jacobi polynomials satisfy the orthogonality condition ð8:2Þ This allows the integrals appearing in (7.17) to be evaluated by inspection.
The harmonic form of the cosmic covariance is obtained from (7.17) by replacing the two-point functions with the harmonic sums given in (6.9), and integrating using (8.2).The first of the two integrals is Z π 0 dβ sin βμðγ; βÞμðγ 0 ; βÞ where the coefficients a l are positive quantities defined in (2.11).With (8.3), it is easy to see that the second of the two integrals in (7.17) is equal to the first, since P l ðcos γÞ ¼ ð−1Þ l P l ðcos γÞ.Hence, from (7.17) and ( 8.3) we obtain a beautiful harmonic expansion of the cosmic covariance In the same way, we can evaluate cosmic variance function e μ 2 ðγÞ defined by (7.18).This function encodes the angular (γ) dependence of the cosmic variance σ 2 cos ðγÞ for the Gaussian ensemble.From (7.18)  where P 2 l denotes the square of a Legendre polynomial and not the associated Legendre function with m ¼ 2. This harmonic form of the cosmic variance was first given in [27] [Eq.(C53)] and was found independently in [10].

IX. MEAN AND VARIANCE OF HD CORRELATION IN MODELS WITH CORRELATED SOURCE SKY LOCATIONS
On the large scale, the universe appears to be fairly isotropic.However, the most likely PTA sources (pairs of supermassive black holes at the centers of merging galaxies) are discrete point sources at specific sky locations.Even if they are distributed via a discrete Poisson process, their apparent angular locations or intensities exhibit correlations.These can occur at the largest angular scales (for example dipole anisotropies [7,46] due to our motion with respect to the average Hubble flow), or they may be at much smaller angular scales.
If PTA sources are distributed in the same way as galaxies, then the relevant quantity is the power spectrum of matter density perturbations, which peaks at a distance scale of about 70 Mpc.However, it appears that even at the peak of the spectrum, these correlations are overwhelmed by "shot noise" arising from the discreteness of the individual PTA sources [47] (Fig. 3).The shot noise produces angular covariance coefficients C l which are l-independent for l ≲ l max .Here, l max ≈ 100 Mpc=10 pc ≈ 10 7 is the ratio of the distance to the closest PTA sources to the characteristic GW wavelength.
To model and understand the effects of these correlations, we construct ensembles of cosmological models in which each realization breaks rotational invariance, but for which the full ensemble is rotationally invariant and thus has no preferred directions.

A. Modeling angular correlations among sources:
A collection of Gaussian subensembles One way to do this is to create an ensemble of Gaussian ensembles.To avoid confusion, we will say that the full ensemble is made up of Gaussian subensembles.In this construction, each of the Gaussian subensembles breaks rotational invariance, but the full ensemble contains all rotated versions of each subensemble, and thus is rotationally invariant.While the full ensemble is no longer Gaussian [31], the key calculational methods can still be used.To compute ensemble averages, we first average over a given Gaussian subensemble, and then average over all subensembles.
Each Gaussian subensemble is constructed as in Sec.VII, but replacing the second moments given in (7.4) with hðf; Here, ψðΩÞ is a real non-negative dimensionless function of the GW source direction, which describes the anisotropic distribution of GW sources within any particular Gaussian subensemble.
Averages within a given subensemble are computed as in previous sections of this paper, using Isserlis's theorem if and as needed.The angle brackets without subscripts hi in (9.1) refer to an average only over the subensemble labeled by ψ.If ψ is not a constant function, then this subensemble breaks rotational invariance [48].In contrast, averages over the full ensemble, which includes many different choices of ψ, will be written with a subscript as hi ψ .
We assume that the full ensemble is described by a set of functions ψ whose first and second moments are given by Here, the angle brackets hi ψ denote a full ensemble average, but in practice we only use this to carry out the final average over all subensembles.We note that knowledge of the first and second moments alone is not enough information for us to compute the ensemble average of any functional.Since only the first and second moments of ψ are known, we can only compute ensemble averages of quantities that are linear or quadratic in ψ.Since this includes the mean and variance of the HD correlation, it is sufficient for our purposes.However, we have no equivalent of Isserlis's theorem to compute higher moments-although each subensemble is Gaussian, our full ensemble is not Gaussian [31].
Because the first moment (9.2) is independent of direction and the second moment (9.3)only depends upon the angle between Ω and Ω 0 , the full ensemble has no preferred directions [49].Because the first moment of ψ is normalized to unity, any quantity linear in H has the same expectation value as previously calculated.Thus, the normalization and interpretation of the spectral function HðfÞ is unchanged.
The function Cð Ω • Ω0 Þ describes the power spectrum of angular fluctuations in the GW background energy density.Using a standard normalization convention (see final paragraph of this section) it can be written as a sum of Legendre polynomials C L P L ðcos βÞ; ð9:4Þ where, as before, cos β ¼ Ω • Ω0 .The expansion coefficients C L are constrained by (9.2), so they cannot have arbitrary values.For example, the sum of ð9:5Þ The first inequality holds because the mean value of a nonnegative quantity is non-negative, the second from completing the square, the third from (9.3), the fourth from Ω • Ω ¼ 1, and the fifth follows from (9.4) and P l ð1Þ ¼ 1.
Our ensemble definition is quite general, so it can be used to model different effects.For example, suppose we want to construct an ensemble of universes which is like the standard Gaussian ensemble, but for which the power spectrum HðfÞ varies in overall amplitude from one subensemble to the next.The first moment normalization (9.2) implies that HðfÞ is the average power spectrum of the complete ensemble (taking into account GW sources in all directions).If each subensemble has exactly that power spectrum, then Cð Ω • Ω0 Þ ¼ 0, so the C l vanish for all l, and we recover the standard Gaussian ensemble.Alternatively, we can construct an ensemble with the same angular properties as the standard Gaussian ensemble, but for which the power spectrum varies in overall amplitude by a factor of 1 about HðfÞ.To obtain this, set Cð Ω • Ω0 Þ ¼ 1, corresponding to C 0 ¼ 4π and C l ¼ 0 for l > 0. We return to this after (9.10).
Note that there is an alternative approach, which we did not follow.For that, each function ψ is decomposed into spherical harmonics, ψðΩÞ ¼ P lm ψ lm Y lm ðΩÞ.The set of ψ used to define the ensemble is then specified via a set of complex coefficients ψ lm ¼ ð−1Þ m ψ Ã l;−m .The properties in (9.2) are equivalent to hψ 00 i ψ ¼ ffiffiffiffiffi ffi 4π p and hψ lm i ψ ¼ 0 for l > 0, and those of (9.3) are equivalent to This approach is equivalent to ours, and provides another way to understand the normalization conventions.Those employing it should beware that the ψ lm cannot be a set of Gaussian random variables with the above first and second moments, because those would not satisfy ψðΩÞ ≥ 0 for every representative function ψðΩÞ in the ensemble.

B. Cosmic (co)variance for the ensemble with correlated source locations
Using this computational framework, we now compute the cosmic variance and covariance for the ensemble with correlated source locations.We do these quantities first, because they are considerably easier to obtain than the total variance and covariance.Those are computed later in this paper, starting in Sec.IX C.
For an extended discussion of the differences between total and cosmic (co)variance, please see [27].
Our starting point is the pulsar-averaged redshift correlation ΓðγÞ given by (7.7) for any realization of the universe.The average of ΓðγÞ over a Gaussian subensemble follows immediately from computing the expected value of (7.7) using (9.1).This simply inserts ψðΩÞ into (7.8)To obtain the second equality, we have removed the twopoint function μðγ; Ω; ΩÞ from the integral, since when the two points are coincident, it is independent of Ω and equal to the HD curve μ u ðγÞ.
To obtain the expected value of ΓðγÞ for the full ensemble, we average (9.6) over ψ, using the first moment (9.2).This gives hΓðγÞi ψ ¼ h 2 μ u ðγÞ; ð9:7Þ which is in agreement with the isotropic result.
To find the cosmic variance and covariance, we need to compute the second moment of Γ.For a given subensemble, we carry out the same calculation which led to (7.13) in Sec.VII.We obtain The first equality follows by repeating the calculation leading to (7.13) (the only change is that two factors of ψ appear), and the second equality follows from (9.6) and the cancellation of the phase of the two-point function (6.8).
We now average (9.8) over the full ensemble using the second moment (9.3).This gives To find the covariance of the full ensemble, we subtract hΓðγÞi ψ hΓðγ 0 Þi ψ , which is obtained from (9.7).Using (7.17) and (9.4), this gives the cosmic covariance for an ensemble with correlated sources, as As C → 0, only the first term on the rhs survives, consistent with (7.16).Below, we simplify the C ≠ 0 case by evaluating the integral and writing the result in terms of the coefficients C l .Before doing so, we examine the special case discussed after (9.5): setting Cðcos βÞ ¼ constant > 0. This corresponds to introducing variations in the overall normalization of the GW background power spectrum HðfÞ.With this choice, different realizations of the universe have the same shape for HðfÞ, but with differing overall scales.Setting Cðcos βÞ to a positive constant corresponds to C 0 > 0 and C l>0 ¼ 0. Using (7.18), this shifts the covariance (9.10) by an amount C 0 ½h 4 μ u ðγÞμ u ðγ 0 Þ þ 2h 4 e μ 2 ðγ; γ 0 Þ=4π.Setting γ ¼ γ 0 , one can see that the variance increases, compared to the standard Gaussian ensemble with C 0 ¼ 0.
Returning to the general case, the integrals over β can be evaluated by using the harmonic decomposition (6.9) of the two-point function μðγ; βÞ.For this, we need the integral of three Jacobi polynomials, which can be written in terms of Clebsch-Gordon coefficients or Wigner 3j symbols [50].Since the first Jacobi polynomial is a normal Legendre polynomial, the integral that we need is Note that the rhs is symmetric in l and l 0 .Letting z ¼ cos β and using (6.9), this implies that Z π 0 sin β dβ P L ðcos βÞμðγ; βÞμðγ 0 ; βÞ To exploit this, we return to (9.10), replacing Cðcos βÞ with its harmonic sum (9.4).
The final form of the cosmic covariance for the ensemble of correlated sky location sources follows immediately, and is The power of −1 arises from the μðγ; βÞμðγ 0 ; βÞ term of (9.10), because P L ðcos βÞ ¼ ð−1Þ L P L ðcos βÞ; P l ðcos γÞ ¼ ð−1Þ l P l ðcos γÞ; P l 0 ðcos γ0 Þ ¼ ð−1Þ l 0 P l 0 ðcos γ 0 Þ: ð9:14Þ In the limit Cðcos βÞ → 0, which is equivalent to C L → 0, we recover the cosmic covariance (8.4) which was computed for uncorrelated sources.Note that the L ¼ 0 term in the sum, which is proportional to e μ 2 ðγ; γ 0 Þ, is easily recovered from (9.13).This is because, for m ≥ 0 which immediately leads to (8.4).

C. Mean and second moment of the HD correlation for a subensemble
To apply the "ensemble of Gaussian subensembles" computational framework, we need the first and second moments of the (total, not pulsar averaged) HD correlation in a given subensemble.The calculations are similar to the ones carried out in Sec.VII for the cosmic variance, but differ in one important way: the total covariance is affected by the pulsar terms, which drop out of the cosmic variance.
The pulsar term changes the frequency-independent redshift response FðΩ; Ω p Þ given in (2.2) to a frequencydependent response, given by Here, the frequency-dependent term is where T p is the light travel time from pulsar p to Earth.In (9.17), the first term is the Earth term, and the second term is the pulsar term, which can add constructively or destructively to the Earth term.Detailed derivations may be found in [5] [Eq.(32)] or [27] [Eq.(C17)].
The correlation ρ pq between pulsars p and q for any universe in any subensemble is defined by (7.5), and given by an expression similar to (7.7): The subsensemble average of (9.18), for a given anisotropy ψðΩÞ, is obtained by using (9.1).Only the first and last terms of (9.18) survive, introducing a delta function of frequency and a delta function on the sphere.Integrating out those delta functions gives the subensemble average By making reasonable assumptions about the pulsar distances and sky locations, this expression can be further simplified.
Consider the integration over Ω in (9.19).The products T Ã T consist of a slowly varying part and a rapidly varying part.We make two assumptions.(A) That ψðΩÞ only varies on angular scales greater than ≈1=fT pulsar ≈ T=T pulsar ≈ 10 −2 radians, where T pulsar is a typical Earth-pulsar or pulsarpulsar light propagation time, and T is the total observation time.(B) That T pulsar is larger than the characteristic coherence time of the GW background (see Eq. (C13) in [27]).Then, as discussed before Eq. ( 45) of [5], the rapidly varying terms integrate to zero.The slowly varying terms in the products T Ã T are 1 if p and q denote distinct pulsars, and are 2 if p and q denote the same pulsar, so ð9:20Þ This doubling arises through the autocorrelation of the pulsar term, and is discussed in more detail below.Assumption (A) about the angular scale of variation of ψ, which corresponds to l values of a few hundred, is sufficient but not necessary.The pulsar terms will also integrate to zero for p ≠ q if the variations in ψ are uncorrelated with the pulsar locations.The average over the full ensemble (which we construct to be rotationally invariant) should guarantee that this condition is satisfied for nonpathological choices of ψ.
With these assumptions, the correlation between pulsars in any Gaussian subensemble is given by ψðΩÞϱ pq ðΩÞ; ð9:21Þ where the squared strain h 2 is defined by (7.9), and ϱ pq ðΩÞ is defined by (2.13).The effect of the 1 þ δ pq factor is to double the correlation of a pulsar with itself, in comparison with the correlation between two pulsars that lie along the same line of sight separated by a distance larger than a GW wavelength; a detailed discussion may be found between Eqs. (C11) and (C13) of [27] and around Eq. ( 45) of [5].
For later use, we now compute the second moment of ρ pq for a particular Gaussian subensemble, starting from (9.18).The calculation is very similar to the one which leads to (7.13).The result is where h 4 is given in (7.14).
It is not surprising that (9.22) takes exactly the same form as [38] [Eq.(2.11)], since the average is over a single Gaussian subensemble.Since (9.21) is linear in ψ, the rhs of (9.22) is quadratic in ψ.Hence, the average of (9.22) over different subensembles can be computed by employing (9.3).

D. The mean and total covariance (of the HD correlation) for the full ensemble
To obtain the mean and total covariance of the HD correlation, we average over the different Gaussian subensembles.This average over different ψðΩÞ yields the full ensemble average.
The mean correlation between pulsars p and q is obtained by averaging (9.21) over the different subensembles using (9.2).This gives the full ensemble average where γ pq is the sky separation angle between pulsars p and q.The first equality in (9.23) indicates the average over all subensembles, the second equality follows from (9.2), which sets ψ → 1 in (9.21), and from (3.3), which shows that the F Ã F and FF Ã terms both average to the same real quantity.The third equality follows from (3.3), which defines the unpolarized HD curve μ u ðγÞ.In the final line, we have adopted the notation of [38] [Eqs.(2.3) and (2.4)] and defined the HD correlation matrix The entries in this matrix of HD correlations are obtained from the values of the HD curve at the angular separations of pulsars p and q, apart from the entries along the diagonal.Those have been doubled by the pulsar term contribution to the autocorrelation.Thus, (9.23) implies that the mean of HD correlation for the full ensemble is identical to that found in the standard Gaussian ensemble: both are proportional to the HD curve.The total covariance of ρ pq is easily obtained for the full ensemble.To compute this, it is helpful to first compute the ensemble average over ψ of hρ pq ihρ rs i.We replace hρ pq i and hρ rs i with (9.21), and then use (9.3) and (9.23) to compute the ensemble average over ψ.Since where μ pq is the HD correlation matrix given in (9.24), and we have defined In all of these equations, p, q, r and s label pulsars, any or all of which could be distinct or identical.The full ensemble average hρ pq ρ rs i ψ ≡ ⟪ρ pq ρ rs ⟫ ψ is obtained from (9.22) by using (9.25) to average the three terms over ψ.This gives The covariance of the full ensemble is obtained by evaluating hρ pq i ψ hρ rs i ψ with (9.23), and then subtracting it from (9.28).This eliminates the first term on the rhs of (9.28), giving the total covariance This is one of our paper's main results, since we will now derive explicit formulas for the different terms.
If the GW source locations in the ensemble are uncorrelated and all Gaussian subensembles have the same overall GW intensity, then D pq;rs → 0. The covariance in (9.29) then reduces to the terms on the first line of the final equality.This is the covariance of the Gaussian ensemble, as given in [38] [Eq.(2.10)].It only depends upon h and not upon h.
If the GW source locations in the ensemble are uncorrelated but the GW intensity varies between subensembles (meaning C l ¼ 0 for l > 0 and C 0 ≠ 0) then the covariance also depends upon h.In this case D pq;rs → ðC 0 =4πÞμ pq μ rs , so that If the GW background is broadband (in the sense of [27] [Eq.(C31)], where the ratio h 4 =h 4 ∝ 1=T → 0 as the observation time T → ∞) then the h 4 term dominates, reflecting the overall variation in GW intensity among different subensembles.
The covariance is useful in several contexts.For example, to reconstruct the HD correlation from experimental data, the relative weighting of correlations within a given angular bin (in γ) is determined from the covariance [38] [Eq.(3.10)].Here, we evaluate the covariance and variance in closed form.
From the Legendre polynomial expansion (9.4) of Cð Ω • Ω0 Þ and the addition theorem (3.4), it follows from (9.27) and (9.31) that The harmonic amplitudes P lm ðΩ p ; Ω q Þ were first studied in [7].There, they were computed in "position space" for small l ≤ 2.Then, they were obtained for all l in Appendix E of [8].In both cases, special pulsar positions were used, with p at the North Pole of the sphere and q along the line of longitude ϕ ¼ 0. Here, we provide an (infinite harmonic sum) expression which is valid for any pulsar pair.We evaluate A lm ðΩ p ; Ω q Þ by substituting the diagonal form (2.10) for F into (9.33)twice, and integrating over Ω.The integrand is a product of three spherical harmonics, with spin weights 0, 2 and −2 [the −2 arises from complex conjugation, see (A2)].Using (A7), this may be written in terms of Wigner 3j symbols, as The final equality of (A4) implies that the summand vanishes unless m 2 ¼ m þ m 1 , so the double sum over m 1 and m 2 may be rewritten as a single sum.
Combining the A lm according to (9.32) gives the spherical harmonic coefficients and flipped the sign of m 2 , to obtain an expression that is explicitly symmetric under interchange of pulsars p and q.This is because changing the sign of the second row of either of the Wigner 3j symbols introduces a factor of ð−1Þ lþl 1 þl 2 , see (A4).
Note that if l ¼ m ¼ 0, then by virtue of (9.15),only the diagonal terms l 1 ¼ l 2 and m 1 ¼ m 2 survive in expressions (9.35) and (9.36).These reduce to the sum in (3.3), giving P 00 ¼ A 00 ¼ ffiffiffiffiffi ffi 4π p μ u ðγÞ.Combining these results provides an explicit expression for D pq;rs , from which the covariance matrix for any sky positions may be obtained.Substituting (9.36) into (9.34)gives where the rotationally invariant function of the pulsar sky positions is The total variance σ 2 tot is obtained from the covariance C pq;rs in (9.29) by setting r → p and s → q.From rotational invariance, σ 2 tot is only a function of the angle γ between pulsars p and q, with cos γ ¼ Ωp • Ωq .So, from (9.26) and (9.29), we obtain For the third equality, we have used (9.24) and δ 2 pq ¼ δ pq , and explicitly indicated the dependence of D on the angle γ.
To determine this completely, we return to (9.37) and evaluate D pq;pq ðγÞ and D pp;qq ðγÞ as sums of Legendre polynomials in γ.The results may be found in (9.43), (9.47), and (9.52).
Since D pq;pq is only a function of γ, it can be pulsar averaged without changing its value.The required quantity is the pulsar average hG Ll 1 l 2 l 3 l 4 ðΩ p ; Ω q ; Ω p ; Ω q Þi pq ∈ γ .The pulsar average of the four spherical harmonics which appear in (9.38) is ð9:40Þ The first equality follows from the recipe (4.4) for pulsar averaging, and the second equality from the standard formula (A7) for the integral of three spin-weighted spherical harmonics.This immediately gives the total variance D pq;pq .Inserting (9.40) into (9.38) and then inserting (9.38) into (9.37)yields ð9:41Þ For this, we have defined constants to simplify the appearance of (9.41) and subsequent equations.Note that the sign disappears, since the summand of (9.41) vanishes unless To obtain (9.41), we made several simplifications.First, we replaced ð−1Þ m 1 þm 2 with ð−1Þ M , because nonvanishing terms must have the bottom row of each Wigner 3j symbol sum to zero (A4), implying that M ¼ −m 1 − m 2 .Second, we used since the only nonzero terms in the sum have . This is because a Wigner 3j symbol vanishes if the bottom row vanishes and the sum of the top row is odd.
Expression (9.41) provides a convenient decomposition of the total variance into a sum of Legendre polynomials of cos γ, where γ is the angle between the directions to pulsars p and q.For this purpose, define a matrix of coefficients d Ll via D pq;pq ðγÞ ¼ There are two alternative approaches which allow further simplifications in the formula for the coefficients d Ll .
In the first approach, note that the Wigner 3j symbols in (9.41) vanish if the sum of the bottom row is nonzero (A4).This means that the summation over M; m; m 1 ; m 2 ; m 3 ; m 4 can be replaced by a summation over M; m; m 1 , with ð9:44Þ which has four infinite sums over l 1 , l 2 , l 3 , l 4 and three finite sums over M, m, and m 1 .
A simpler and more symmetric expression can be obtained by returning to (9.41) and using the Wigner 6j symbol to carry out the sums over M; m; m 1 ; m 2 ; m 3 ; m 4 .The Wigner 6j symbol satisfies the equation & j 1 j 2 j 3 j 4 j 5 j 6 We make the following substitutions into (9.45): Then, we exploit properties of the Wigner 3j symbol (A4).Swapping any pair of columns or inverting the signs of the bottom row multiplies the Wigner 3j symbol by ð−1Þ S , where S denote the sum of the top row.Using these, we arrive at This formula for the numerical coefficients is very pretty, and it may be possible to some terms in this sum by exploiting further symmetries of the Wigner 3j and 6j symbols.
To complete the evaluation of the total variance (9.39), we also need to evaluate evaluate D pp;qq ðγÞ.Return to the definition (9.27),where from (2.13) and (2.7) the HD integrand is and the final equality follows from (2.7).If we let z ¼ Ω • Ωp , then Y 00 ðΩ p ÞY Ã 00 ðΩÞ: ð9:49Þ For the second equality, we have expressed the quadratic polynomial in terms of Legendre polynomials, and for the third equality, we have used the addition theorem (3.4) for l ¼ 0, 1, and 2. Thus, the expansion coefficients P lm ðΩ p ; Ω p Þ given in (9.32) are Corresponding expressions for P lm ðΩ q ; Ω q Þ are obtained by setting Ω p → Ω q in (9.50).We now complete the evaluation of D pp;qq ðγÞ.From (9.34), the quantity required is where the values are taken from (9.50) with corresponding expressions with Ω p → Ω q , and the sum over M is done via the addition theorem (3.4)

X. CONCLUSION
We have shown how harmonic analysis, based on the diagonal decomposition (2.10), makes it straightforward to calculate the most important quantities of interest for pulsar timing arrays.We then use these methods to model universes whose GW source sky positions have nontrivial angular correlations.To do this modeling, we build "statistically isotropic ensembles" from anisotropic Gaussian subensembles.This leads to simple equations for the cosmic variance/covariance, and for the total variance/covariance.Investigations for realistic cosmological models are underway [51], though for large l these effects may be too small to be observable in the near future.
Note added.Recently, we learned that Agarwal and Romano had independently carried out the calculation of the cosmic variance for the ensemble with nontrivial angular correlations [52].Their results are consistent with those obtained in Sec.IX B of this paper; in fact we have unified our notation so that the results may be easily compared.Subsequently, we applied these methods to estimate the impact of galaxy discreteness and clustering on the HD correlations [47].

ACKNOWLEDGMENTS
I would like to thank Joe Romano and Neil Cornish for helpful discussions, Nastassia Grimm for persuading me that it makes sense to construct an ensemble of Gaussian subensembles, Serena Valtolina for checking some of the calculations and discovering that (9.44) could be written in terms of the Wigner 6j symbol (9.47),Daniel Pook-Kolb for numerically checking (2.10) and computing a table of d Ll values, Marc Favata for confirming that (A7) can be violated if s 1 þ s 2 þ s 3 ≠ 0, Neha Anil Kumar for pointing out an equation equivalent to (2.10) in Sec.V.B of [14], and Joe Romano and Deepali Agarwal for comparing results and converging notation, and Gabi Richardson for assistance in manuscript preparation.I also thank the anonymous referee for useful suggestions and for identifying several oversights and omissions.
Hence, the symbol is invariant under (a) any even permutation of columns or (b) any odd permutation of the columns accompanied by a sign flip of the bottom row.

Spin-2 harmonics used in this paper
These vanish for l < 2 and may be obtained for l ≥ 2 by taking derivatives of the normal (spin-weight 0) spherical harmonics: where the "edth" spin-raising operators are 6. Integral of three spherical harmonics A7) may not hold: the lhs may be nonzero, but the rhs vanishes.For example, R dΩ 0 Y 00 ðΩÞ 0 Y 22 ðΩÞ 2 Y 2;−2 ðΩÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 3=32π p , whereas the rhs of (A7) vanishes for s 1 ¼ s 2 ¼ 0 and s 3 ¼ 2. In such cases, the integral may be evaluated using the method of [53] (Appendix A).

APPENDIX B: DERIVATION OF THE DIAGONAL FORM OF
Here, we derive the diagonal form of FðΩ; Ω p Þ given in (2.10), following an approach inspired by [8] (Sec.III.D).We also explain how (2.10) can be checked/verified directly, by explicitly carrying out the sums.Lastly, we perform two simple sanity checks.
To verify the diagonal form in (2.10) directly, use the addition theorem for spin-weighted harmonics [8] [Eqs.(A9)-(A11) with s ¼ 2; s 0 ¼ 0] to carry out the sum over m.Then, use and [derived in [9] Eq. ( 42)] to complete the sum over l.Some algebra with trigonometric identities leads directly to (2.2).
[In (B1) and (B2), P m l ðzÞ denotes an associated Legendre function, not the square of a Legendre polynomial.] To derive the diagonal form in (2.10), begin with Fðẑ; Ω p Þ as given in (2.8).This the response of a pulsar at an arbitrary sky direction Ωp to a GW with direction ẑ.We rotate this pattern to obtain the response to a GW with arbitrary direction Ω. (Here, and in what follows, it is often helpful to write the arguments of F and spherical harmonics as unit vectors rather than as coordinates on the sphere.) There are many different rotations that will bring ẑ to Ω.For the reasons explained in Sec.II, we select the unique rotation that consistently maintains the directions of the polarization vectors m and n, as defined by (2. This matrix acts from the left, on column vectors whose three entries are the x, ŷ and ẑ components.(In this appendix, γ and β denote rotation angles.Elsewhere in the paper, they denote the angles between pairs of pulsars or between pairs of GW sources.)By inspection, the rotation in (B3) acting on ẑ gives Thus, to obtain the GW direction Ω ¼ Rðα; β; γÞẑ as given in (2.1), we must set β ¼ −θ and γ ¼ −ϕ.Note that α can take any value.This is also obvious from inspection of (B3), since the rightmost matrix leaves ẑ invariant.However, there is only a single value of α which yields the correct polarization vectors m and n, as given in (2.3).To see this, act on nðẑÞ with the rotation Rðα; −θ; −ϕÞ.The ẑ component of Rðα; −θ; −ϕÞ n is sin θ sinðϕ − αÞ.Since nð ΩÞ has no ẑ component, we must have α ¼ ϕ þ Nπ for N integer.Only even N, equivalent to N ¼ 0, maintains the orientation of n.Thus, the only acceptable rotation which carries ẑ to Ω and which carries nðẑÞ to nð ΩÞ is This rotation matrix also carries mðẑÞ to mð ΩÞ.
We emphasize this point one last time.For an arbitrary rotation R, FðR Ω; R Ωp Þ ≠ Fð Ω; Ωp Þ. Equality is only obtained for rotations that satisfy R mð ΩÞ ¼ mðR ΩÞ and R nð ΩÞ ¼ nðR ΩÞ.In words: the pulsar response F is only invariant under simultaneous rotations of the GW source and pulsar directions which also preserve the polarization vectors m and n.This consistency requirement has also been noted by others [14] [ Footnote 4].
From here, it is straightforward.We first express the unrotated response function (2.8) as a sum of spherical harmonics with expansion coefficients q l .Because the ϕ p dependence in the first equality is e 2iϕ p , the sum only includes spherical harmonics with m ¼ −2, which implies that the q l vanish if l < 2. For l ≥ 2 they are which follows immediately from (B2).Next, we rotate the response function, by rotating the spherical harmonics.Since the rotation matrix (B5) preserves the polarization directions, rotational invariance implies that where the final equality comes from (B6).Setting Rẑ ¼ Ω in (B8), and then noting that, since the equation holds for all Ωp , we can send Ωp → R −1 Ωp , we obtain where D l mm 0 is the Wigner D-matrix.(For fixed l, the Y lm form a 2l þ 1-dimensional vector space representation of the group SOð3Þ.Thus, the rotated Y l;−2 is a sum of harmonics with the same l and all allowed m values [42] [Pg. 51]. The second equality of Eq. (B9) is obtained using D l m 0 m ðRÞY lm 0 ðΩÞ ðB10Þ [Eq. (16.52)], which is consistent with our choice of Euler angles in (B3) and with [46] [Eqs.(7.3)-(7.7)].Note that the corresponding relationship in [42] [Eqs.(2.43) and (2.45)] replaces D l m;−2 ðR −1 Þ in (B9) with D l m;−2 ðRÞ.This is equivalent: since [42] uses active rather than passive rotations, the signs of the Euler angles and their ordering are inverted, swapping R and R −1 , see [42] [Eq.(1.54)] and [55] [Eq.(6.39)].
The inverse of the rotation matrix (B5) can be found by inspection of (B3), and is R −1 ¼ Rðϕ; θ; −ϕÞ.This rotation carries Ω to ẑ, while also preserving the polarization vectors.
APPENDIX C: LINEAR POLARIZATION COMPONENTS OF THE TWO-POINT FUNCTION Some calculations (see [27] for examples) are best carried out using two-point functions for linear polarization components, written μ þþ , μ ×× , μ ×þ , and μ þ× .Here, we extract these from the complex two-point function μðγ; Ω; Ω p Þ.
The linear polarization components of the two-point function (C6) satisfy a relationship which is useful for some calculations [47].If all four terms on the lhs of (C6) are squared and summed, one obtains the "gauge-independent" relationship While the summands on the lhs depend upon the phases χ and χ, on the rhs these have canceled out: the sum only depends upon γ and β.

0 sin β dβ cos β 2 8 P
, if we send l → l − 2 and l 0 → l 0 − 2 in (8.1), the orthogonality condition takes the form Z π :38Þ It should be possible to express G Ll 1 l 2 l 3 l 4 as a function of Legendre polynomials of the dot products Ωp • Ωr , Ω p • Ωs , Ωq • Ωr , and Ωq • Ωs .E. The total Hellings and Downs variance for the full ensemble