HARMONIC SPACE ANALYSIS OF PULSAR TIMING ARRAY REDSHIFT MAPS

Pulsar timing arrays (hereafter PTAs) are galactic-scale gravitational wave (GW) detectors based on the precise timing of millisecond pulsars across the sky (Foster & Backer 1990). The nanohertz frequency band of GWs accessible to PTAs has several potential production mechanisms, the most prominent of which is due to the inspiral of subparsec supermassive binary black holes (SMBBHs; see Lommen 2015, and references therein).

SMBBHs with chirp mass ${ \mathcal M }\gt {10}^{8}\,{M}_{\odot }$ at redshifts z ≲ 2 are expected to produce most of the signal (e.g., Sesana et al. 2008). Since there should be many such sources evolving over times much longer than human timescales, the GW signal is expected to form a stochastic background with considerable source confusion. However, individual strong sources may stand out (Sesana et al. 2008; Ravi et al. 2012).

A passing GW induces compression and rarefaction of spacetime along its polarization axes. Periodic signals such as rays of light or pulse trains propagating through this region will be blue- or redshifted according to the strain of the GW. For periodic signals with frequency much higher than that of the GW, the shift will build up, producing a potentially measurable effect. This is the principle on which several models of GW detection are founded, including interferometers such as ligo (Abbott et al. 2016) and lisa (eLISA Consortium 2013) as well as for PTAs (Lommen 2015). There are three PTA consortia: epta (Lentati et al. 2015), nanograv (Arzoumanian et al. 2016), and ppta (Shannon et al. 2015). They combine together to form the ipta (Verbiest et al. 2016).

PTAs search for integrated red- and blueshifts produced by GWs passing the Earth through the careful timing of a network of millisecond pulsars across the sky. Each millisecond pulsar produces an extraordinarily regular train of high-frequency pulses. If this pulse train is redshifted by a GW with typical strain ≲10⁻¹⁴ (e.g., Lommen 2015), no effect will be immediately visible, but after the passage of many pulses, a difference between the expected and actual time of arrival of pulses will become apparent. This timing residual is the basic measurable quantity for a PTA.

A GW of a given polarization will induce red- and blueshifts according to the geometry set by the direction of propagation of the GW and the projection of its polarization axes onto the sky. In order to sample this effect as fully as possible, PTAs time many millisecond pulsars across the sky and search for a correlation in their timing residuals which reflects the redshift pattern induced by GWs.

The expected form of this correlation is the Hellings & Downs curve (Hellings & Downs 1983), which was originally derived for a statistically isotropic unpolarized Gaussian random field of GWs. It also represents the expected correlation pattern for a single SMBBH source of GWs (Cornish & Sesana 2013).

However, the gravitational wave background (GWB) expected to be produced by a population of inspiraling SMBBHs will be neither completely dominated by a single source nor a completely stochastic Gaussian field. In general, it should be somewhere in between (e.g., Sesana et al. 2008).

Although much work has made use of the assumption that a stochastic background would have Gaussian statistics, single sources should not be neglected in the PTA search for GWs (Rosado et al. 2015). This is because the distribution of SMBBH sources is such that the rarest brightest sources dominate the signal in the GWB (Sesana et al. 2008; Kocsis & Sesana 2011; Ravi et al. 2012; Cornish & Sesana 2013; Roebber et al. 2016).

In light of this, it is of interest to search for angular information in the GWB. PTAs can be likened to a collection of GW antennas: their angular resolution is limited but not nonexistent. This has been taken advantage of in the attempt to search for individual sources and hotspots (e.g., Corbin & Cornish 2010; Sesana & Vecchio 2010; Babak & Sesana 2012; Simon et al. 2014). Additionally, recent works have characterized the correlation patterns expected for statistically anisotropic backgrounds made up of a large number of sources (Mingarelli et al. 2013; Taylor & Gair 2013) as well as attempting to map general GWBs (Cornish & van Haasteren 2014; Gair et al. 2014).

Many of these recent works have focused on estimating the distribution of GW signals produced by the source population, either in terms of power or components of the GW tensor. However, the GW strain is not directly measured by PTAs. The large effective beam patterns smear power out across the sky, mixing contributions from different sources. Furthermore, since GWs are tensors and the timing residuals measured by PTAs are scalars, there are components of the strain that cannot be measured (Gair et al. 2014). Both of these complications can be sidestepped by working with the maps of the theoretical timing residuals or equivalently, the redshifts induced by the passing GWs.

In this paper, we consider an alternative analysis of the GWB in the PTA band, inspired by standard cosmic microwave background (CMB) methods. Our primary quantity of interest is the redshift induced in all directions on the sky by GWs passing the Earth. This is related to pulsar timing residuals in the following fashion:

• 1.  
  Sampling the redshift field in a direction $\hat{p}$ gives the amount by which the pulse train of a pulsar at $\hat{p}$ is redshifted or blueshifted due to the influence of GWs passing the Earth.
• 2.  
  Integrating the redshift at $\hat{p}$ gives the shift in the pulsar's timing residuals due to GWs passing the Earth (the "Earth term"). Since we limit our discussion to circular and non-evolving GW sources, the integrals are trivial.

Furthermore, we initially analyze redshift maps in harmonic space, and transform back to real space when considering the implications. This approach may not be practical for experimental analysis and we present it primarily as an alternative framework for understanding the angular information in the GWB.

In Section 2 we review the standard mathematical formalism underlying GWs produced by circular, slowly inspiraling binary systems and their measurement by PTAs, and produce example maps of the redshift patterns produced by various GWBs. In Section 3 we present our harmonic-space analysis of redshift maps and specifically discuss two limiting cases: a single GW source and a statistically isotropic Gaussian random GWB. In Section 4 we discuss the relation between the two-point function in real and harmonic space and present a case where the harmonic analysis provides insight into real-space quantities: how variance in the power spectrum affects the shape of the Hellings & Downs curve. In Section 5 we discuss the degree to which the power spectrum is measurable in an ideal PTA. And finally, in Section 6 we present our conclusions and discuss future directions.

A GW is a transverse plane wave propagating as spatial perturbations in the metric. It has a spin-2 symmetry and two polarizations (Maggiore 2008):

$$\begin{eqnarray}&&{h}_{{ij}}(t,\hat{k})={h}^{+}(t)\,{e}_{{ij}}^{+}(\hat{k})+{h}^{\times }(t)\,{e}_{{ij}}^{\times }(\hat{k}),\end{eqnarray} \tag{ 1 }$$

where h⁺ and h^× are the amplitudes of the two polarizations, ${e}_{{ij}}^{+}$ and ${e}_{{ij}}^{\times }$ are the polarization tensors, and $\hat{k}$ is the direction of propagation of the wave. Sub- and superscripts i, j are written using the Einstein summation notation and denote the tensorial nature of GWs.

The geometry of an incoming GW can be written in terms of a radial vector in the direction of propagation of the GW, and two vectors perpendicular to it which define a basis for the polarization of the wave. Our choice of conventions follows Gair et al. (2014):

$$\begin{eqnarray}\hat{k} & = & \sin \theta \cos \phi \,\hat{x}+\sin \theta \sin \phi \,\hat{y}+\cos \theta \,\hat{z}\\ \hat{l} & = & \cos \theta \cos \phi \,\hat{x}+\cos \theta \sin \phi \,\hat{y}-\sin \theta \,\hat{z}\\ \hat{m} & = & -\,\sin \phi \,\hat{x}+\cos \phi \,\hat{y}.\end{eqnarray} \tag{ 2 }$$

If we consider $\hat{k}$ to be a radial vector along the axis of propagation, the location of the GW source is in the $-\hat{k}$ direction, or equivalently at the angle on the sky $(\pi -\theta ,\phi +\pi )$. The perpendicular vectors $\hat{l}$ and $\hat{m}$ are vectors in the $\hat{\theta }$ and $\hat{\phi }$ directions defining the plus and cross polarizations of the incoming GW:

$$\begin{eqnarray}{e}_{{ij}}^{+}(\hat{k}) & = & {\hat{l}}_{i}{\hat{l}}_{j}-{\hat{m}}_{i}{\hat{m}}_{j}\\ {e}_{{ij}}^{\times }(\hat{k}) & = & {\hat{l}}_{i}{\hat{m}}_{j}+{\hat{m}}_{i}{\hat{l}}_{j}.\end{eqnarray} \tag{ 3 }$$

These are all general properties of GWs, but we are interested in GWs generated by SMBBHs, which can be described more closely. In particular, we restrict ourselves to the case where the binary is circular and very slowly evolving, so that we can ignore its evolution on observational timescales. GWs of this form can be described by four additional parameters: $({ \mathcal A },\iota ,\psi ,{{\rm{\Phi }}}_{0})$, as described in the following paragraphs (e.g., Cutler & Flanagan 1994; Sesana & Vecchio 2010).

The amplitude ${ \mathcal A }$ contains information about the non-angular degrees of freedom of the binary:

$$\begin{eqnarray}&&{ \mathcal A }=2\displaystyle \frac{{(G{ \mathcal M })}^{5/3}}{{c}^{4}D}{(\pi {f}_{\mathrm{emit}})}^{2/3},\end{eqnarray} \tag{ 4 }$$

where ${ \mathcal M }={({m}_{1}{m}_{2})}^{3/5}/{({m}_{1}+{m}_{2})}^{1/5}$ is the chirp mass of the binary system, D is the proper distance, and f_emit is the frequency of the GW in the binary's rest frame.

The inclination ι of the binary tells us the relative contribution of each polarization. A face-on or face-off binary is circularly polarized, and produces equal quantities of the plus and cross polarizations. An edge-on binary only produces plus polarization, and can be considered to be linearly polarized. A general binary is somewhere in-between, and its GW is elliptically polarized. For an inclination ι, the contributions to the plus (a) and cross (b) polarizations can be expressed as

$$\begin{eqnarray}a(\iota ) & = & 1+{\cos }^{2}\iota \\ b(\iota ) & = & -\,2\cos \iota .\end{eqnarray} \tag{ 5 }$$

The angle ψ encodes the transformation between GW polarizations between the source coordinate system and that of the observer. It gives the degree to which the plane of the binary is misaligned with the $(\hat{l},\hat{m})$ basis given above, which leads to mixing between the different polarizations:

$$\begin{eqnarray}{h}_{+}^{\prime } & = & {h}_{+}\cos 2\psi +{h}_{\times }\,\sin 2\psi \\ {h}_{\times }^{\prime } & = & -\,{h}_{+}\,\sin 2\psi +{h}_{\times }\,\cos 2\psi .\end{eqnarray} \tag{ 6 }$$

The mixing takes the form of a rotation by 2ψ since GWs are spin-2: a pure + mode becomes purely × if the $(\hat{l},\hat{m})$ coordinate system is rotated by 45°. The angle ψ and the angles giving the location of the GW source (θ, ϕ) are defined in terms of the coordinate system of the observer, and can be changed by a rotation of the coordinate axes.

The overall temporal phase of the binary is given by

$$\begin{eqnarray}&&{\rm{\Phi }}(t)={\int }_{0}^{t}2\pi f(t^{\prime} ){dt}^{\prime} \approx 2\pi {ft}+{{\rm{\Phi }}}_{0},\end{eqnarray} \tag{ 7 }$$

where Φ₀ is the initial phase of the binary. To make the approximation in Equation (7), we assume non-evolving circular binaries.

Altogether, the components of GWs produced by a non-evolving circular binary can be written as

$$\begin{eqnarray}{h}_{+}(t,\hat{k}) & = & { \mathcal A }[a\cos 2\psi \cos {\rm{\Phi }}(t)+b\sin 2\psi \sin {\rm{\Phi }}(t)]\\ {h}_{\times }(t,\hat{k}) & = & { \mathcal A }[b\cos 2\psi \sin {\rm{\Phi }}(t)-a\sin 2\psi \cos {\rm{\Phi }}(t)].\end{eqnarray} \tag{ 8 }$$

Assuming that the binaries are circular and non-evolving, as in Equation (7), this may be easily written in frequency space as

$$\begin{eqnarray}{\tilde{h}}_{+}(f,\hat{k}) & = & \displaystyle \frac{{ \mathcal A }}{2}[a\cos 2\psi (\cos {{\rm{\Phi }}}_{0}+i\sin {{\rm{\Phi }}}_{0})\\ & & +b\sin 2\psi (\sin {{\rm{\Phi }}}_{0}-i\cos {{\rm{\Phi }}}_{0})]\\ {\tilde{h}}_{\times }(f,\hat{k}) & = & \displaystyle \frac{{ \mathcal A }}{2}[b\cos 2\psi (\sin {{\rm{\Phi }}}_{0}-i\cos {{\rm{\Phi }}}_{0})\\ & & -a\sin 2\psi (\cos {{\rm{\Phi }}}_{0}+i\sin {{\rm{\Phi }}}_{0})],\end{eqnarray} \tag{ 9 }$$

where f is the positive frequency associated with the binary, and $\tilde{h}(f)$ denotes the Fourier transform with respect to time of h(t). Since h₊(t) and h_×(t) are real-valued functions there are also negative frequency terms given by ${\tilde{h}}_{+,\times }(-f)={\tilde{h}}_{+,\times }^{* }(f)$.

As this GW (assumed to originate far outside our galaxy) passes a pulsar and the Earth, the pulse train seen on Earth gains a frequency shift of

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\nu }{\nu }=z(t,\hat{p},\hat{k})=\displaystyle \frac{1}{2}\displaystyle \frac{{p}^{i}{p}^{j}}{1+\hat{k}\cdot \hat{p}}[{h}_{{ij}}({t}_{\mathrm{psr}},\hat{k})-{h}_{{ij}}({t}_{\mathrm{Earth}},\hat{k})],\end{eqnarray} \tag{ 10 }$$

where the direction to the pulsar is written as $\hat{p}$. Frequency shifts will be of the same order of magnitude as h, that is ≲10⁻¹⁵. This is too small to measure. However, over many cycles, the frequency shift will affect the time of arrival of the pulses:

$$\begin{eqnarray}&&r(t,\hat{p})={\int }_{0}^{t}{dt}^{\prime} z(t^{\prime} ,\hat{p}),\end{eqnarray} \tag{ 11 }$$

producing the GW contribution to PTA timing residuals. The amplitude of r will be of order h/f ≲ 100 ns for waves with f ∼ 10 nHz.

Since the GWs will pass through our entire galaxy, $z(t,\hat{p})$ (and equivalently, r) can be split into two terms: the term due to the metric disturbance at the Earth, ${h}_{{ij}}({t}_{\mathrm{Earth}},\hat{k})$, and the term at the pulsar. Earth terms due to the same GW will be correlated between different points on the sky (different pulsars), but pulsar terms will depend on the distance between the Earth and the pulsar. Absent detailed information about pulsar distances, and assuming that the sources do not evolve significantly in frequency between the time that the waves pass the Earth and all pulsars, pulsar terms can be modeled as a term of the same magnitude as the Earth term but with a random additional phase. For simplicity, the rest of the paper will concentrate on the Earth terms, which are correlated on the sky, although they can also be considered as a form of self-noise, which would enter the calculations in Section 5.

Assuming circular binaries, we write

$$\begin{eqnarray}&&{\tilde{z}}_{\mathrm{Earth}}(f,\hat{p},\hat{k})=\displaystyle \frac{1}{2}\displaystyle \frac{{p}^{i}{p}^{j}}{1+\hat{k}\cdot \hat{p}}{\tilde{h}}_{{ij}}(f,\hat{k}),\end{eqnarray} \tag{ 12 }$$

where ${\tilde{h}}_{{ij}}(f,\hat{k})$ is of the form given in Equation (9). This is a complex scalar field. Calculating the total redshift induced in any one direction $\hat{p}$ requires integration over all h_ij coming from all directions $\hat{k}$.

We present two examples of ${\tilde{z}}_{\mathrm{Earth}}$ in a single frequency bin. Figure 1 shows the redshift map produced by two SMBBH sources of the same amplitude but different inclinations and random other parameters. Figure 2 is an example of a likely GWB for a single frequency bin. It is generated from the population models of Roebber et al. (2016). Every binary black hole is assigned a random set of parameters. Although the population contains 150,000 GW sources with ${ \mathcal M }\gt {10}^{7}\,{M}_{\odot }$, relatively few are visible in the maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We consider two toy model GWBs which can be considered as limiting cases for a stochastic GWB produced by a population of SMBBHs with no underlying anisotropy. The first example is for a single source of GWs, which is an idealization of the case where the GW power in a frequency bin is dominated by a single bright source. The second example is the canonical case where the GWB is a stochastic Gaussian random field. This represents the opposite limit of a confusion background, which has no visible individual sources.

3.1. A Single Gravitational Wave Source

For the case of a single source, we will consider a single inspiraling pair of SMBBHs located at the north pole and aligned with our (l, m) coordinate choices, so that ψ = 0 and ϕ = 0. This choice will allow us to do the calculations in a simple form; all other possible single sources can be reproduced by applying a rotation at the end.

In this coordinate system, the direction of the GW propagation is $\hat{k}=-\hat{z}$ and the vectors defining the polarization are $(l,m)=(-\hat{x},\hat{y})$. Plugging these definitions, Equations (3), and (1) into Equation (12) and considering the response in all directions $\hat{r}=(\theta ,\phi )$ produces the redshift induced across the sky by a single source at the north pole:

$$\begin{eqnarray}&&\tilde{z}(\theta ,\phi )=\displaystyle \frac{1}{2}(1+\cos \theta )(\cos 2\phi \,{\tilde{h}}_{+}-\sin 2\phi \,{\tilde{h}}_{\times }).\end{eqnarray} \tag{ 13 }$$

This is a continuous field everywhere except in the direction of the source, where the θ term is constant, but the ϕ term is undefined due to rapid oscillation at small θ. See Figure 3.

Zoom In Zoom Out Reset image size

Since Equation (13) is a scalar field, it can be represented as the sum of spherical harmonics. Doing this expansion (see Appendix) produces

$$\begin{eqnarray}&&{a}_{{lm}}=2\pi \sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l-2)!}{(l+2)!}}({\tilde{h}}_{+}\pm i{\tilde{h}}_{\times }),\ \mathrm{for}\ m=\pm 2.\end{eqnarray} \tag{ 14 }$$

Note that the subscripts l, m here are the usual spherical harmonic labels and not tensor indices. Interestingly, the a_lm only exist for m ± 2. This is a reflection of the four stripes seen in the half-beachball form of the pulsar response function (see Figure 3). Fundamentally, this is due to the spin-2 nature of GWs—a rotation of the $(\hat{l},\hat{m})$ coordinate system by 180° must produce the same result.

Although the underlying form of the redshift pattern is fundamental, the representation in spherical harmonics is a result of our choice to place the GW source at the north pole. A source located elsewhere in the sky can be expressed by a rotation of Equation (14). This will mix between m components, so that a generic source will require a full set of spherical harmonics to reproduce its response function. However, since rotations of spherical harmonics cannot transform one l to another, the scaling of a_lm with l will remain consistent.

A statistical description of the l-scaling of z can be found in the angular power spectrum (Dodelson 2003):

$$\begin{eqnarray}&&{C}_{l}=\displaystyle \frac{1}{2l+1}\displaystyle \sum _{m}{a}_{{lm}}^{* }\,{a}_{{lm}}.\end{eqnarray} \tag{ 15 }$$

For the case of a single source, this becomes

$$\begin{eqnarray}&&{C}_{l}=\displaystyle \frac{2\pi ({\tilde{h}}_{+}^{2}+{\tilde{h}}_{\times }^{2})}{(l+2)(l+1)(l)(l-1)},\ \mathrm{for}\ l\geqslant 2.\end{eqnarray} \tag{ 16 }$$

This is a steeply decreasing function of l. Since it is only a function of l it does not vary under rotations. Therefore, Equation (16) holds for any single SMBBH source of GWs with polarizations ${\tilde{h}}_{+}$ and ${\tilde{h}}_{\times }$.

3.2. A Statistically Isotropic Gaussian Random Field Gravitational Wave Background

The second case that we consider is the case where the GWB is a statistically isotropic Gaussian random field. This represents an idealization of a stochastic background produced by many sources, and similar models have frequently been considered in the PTA literature. Our discussion will follow Burke (1975).

This kind of background is the most similar to the CMB. However, a major difference is that the GWB is stationary but not time-invariant on observational timescales. (This is because it is produced by rotating SMBBHs, which cause the background to rotate through polarizations). When we Fourier transform the data to work with a single frequency bin, the resulting maps will be complex, unlike the real CMB maps.

A GWB produced by SMBBHs produces a redshift field equal to the sum over the redshift field produced by each source. Every source will produce a redshift field of the form of Equation (13), but with an additional random rotation, which sets (θ, ϕ, ψ) to random new values. For a field with maximal source confusion, we consider independent sources of similar amplitude along every line of sight.

As we incoherently add sources, the value of each of the a_lms will change. However, since the l-dependence is left unchanged after a rotation, the new a_lms will be given by a sum of terms with varying complex amplitude, but which all scale with l in the same way as Equation (14). In this way we can see that adding sources preserves the shape of the average power spectrum. An example is shown in Figure 4.

Zoom In Zoom Out Reset image size

The amplitude of the average power spectrum is however strongly affected by the number (and strength) of GW sources. For fixed l, adding randomly located sources can be modeled as a random walk in the amplitudes of each mode. The random walk will have a mean of zero, but the variance will increase proportionally to the number of sources. Since the power spectrum is the variance of the a_lm, a Gaussian random field produced by N_src identical sources should have an underlying power spectrum of the form

$$\begin{eqnarray}&&{C}_{l}\propto \displaystyle \frac{{N}_{\mathrm{src}}}{(l+2)(l+1)(l)(l-1)}.\end{eqnarray} \tag{ 17 }$$

In other words, adding many sources of similar amplitude randomly and incoherently will produce a field whose harmonic decomposition is made up of terms arbitrarily drawn from a distribution given by Equation (17). In real space, this means that as the number of sources increase, points separated by given angle will maintain an average correlation, but actual values will be randomly distributed according to a Gaussian distribution.

A fully Gaussian random field is shown in Figure 5. By comparing Figure 5 and the low-frequency population model GWB shown in Figure 2, we see that the population produces a mostly Gaussian field, but with small artifacts around the brightest sources. While the distribution of source amplitudes in a real population is steeply decreasing, it is still true that a relatively small number of sources produce the majority of the signal.

Zoom In Zoom Out Reset image size

It is important to recall that while Equation (17) gives the expectation value for the power spectrum of the field, any single realization will only approximately reproduce it. We will discuss this further in Section 4.3.

In this section we will discuss the relationship between the power spectrum of the redshift map and the Hellings & Downs curve. We will explore how variance around the fiducial power law power spectrum affects the shape of the two-point correlation function in the GWB models previously discussed. Since redshifts may be readily converted into timing residuals independently of angle, the following analysis applies to both.

4.1. The Power Spectrum is the Harmonic Transform of the Hellings & Downs Curve

In the standard PTA analysis, timing residuals from different pulsars are correlated. This process should average away effects of noise, which is expected to be uncorrelated between different pulsars (e.g., Lommen 2015). Correlations due to passing GWs in the timing residuals of two pulsars separated by an angle θ are expected to take the form of the Hellings & Downs curve (Hellings & Downs 1983):

$$\begin{eqnarray}&&C(\theta )=\displaystyle \frac{1}{2}\left\{1+\displaystyle \frac{3}{2}(1-\cos \theta )\left[\mathrm{ln}\left(\displaystyle \frac{1-\cos \theta }{2}\right)-\displaystyle \frac{1}{6}\right]\right\}.\end{eqnarray} \tag{ 18 }$$

This is a real-space two-point correlation function (also sometimes referred to as an overlap reduction function). It is an especially useful statistic for a Gaussian random field, since the statistics of such a field can be described entirely by the mean (one-point function) and standard deviation (two-point function), e.g., Allen & Romano (1999). For symmetric fields such as the GWB or the CMB, the mean vanishes, ensuring that the two-point function contains the entire statistical information of the field.

If we consider taking the two-point function of a large number of points across the sky (e.g., the pixels in a map such as those shown in Figure 1), it becomes sensible to define a harmonic-space analog of the two-point correlation function. This function is the power spectrum discussed earlier, and the conversion between the two forms may be written (Dodelson 2003) as

$$\begin{eqnarray}&&C(\theta )=\displaystyle \sum _{l=0}^{\infty }{C}_{l}\displaystyle \frac{2l+1}{4\pi }{P}_{l}(\cos \theta ),\end{eqnarray} \tag{ 19 }$$

where C_l is given by Equation (17).

A proof for this relation is shown in Gair et al. (2014), although they use different notation: their C_l is a constant since they are concerned with the power spectrum of the GW point source distribution. Their additional factor of ${N}_{l}^{2}$ is equivalent to the l-scaling in our choice of C_l, and represents the effect of the pulsar response function.

An additional concern is the normalization of Equation (17) required to reproduce the standard form of the Hellings & Downs relation. This can be found from Equation (19), as done by Gair et al. (2014), or by considering Parseval's theorem for spherical harmonics:

$$\begin{eqnarray}\displaystyle \sum _{l=0}^{\infty }\,\displaystyle \sum _{m=-l}^{l}| {a}_{{lm}}{| }^{2} & = & {\displaystyle \int }_{{S}^{2}}d{\rm{\Omega }}\,| \tilde{z}(\theta ,\phi ){| }^{2}\\ \displaystyle \sum _{l=2}^{\infty }(2l+1){C}_{l} & = & 4\pi \,C(0).\end{eqnarray} \tag{ 20 }$$

Doing the sum produces a factor of 1/3, so the normalization of the power spectrum required to satisfy this constraint is

$$\begin{eqnarray}&&{C}_{l}=\displaystyle \frac{6\pi }{(l+2)(l+1)(l)(l-1)}.\end{eqnarray} \tag{ 21 }$$

Since the Hellings & Downs curve is the map space version of the expected form of the angular power spectrum, any effects which modify the power spectrum can be converted into potentially measurable effects on the two-point correlation function.

4.2. Variance in the Hellings & Downs Curve for a Single Gravitational Wave Source

4.3. Variance in the Hellings & Downs Curve for a Gaussian Random Field

For a statistically isotropic Gaussian random GW field, the expectation value of the power spectrum will be of the form of Equation (17). However, we are able to observe only one realization of the GWB.

For a Gaussian random background, all multipole moments a_lm are drawn from a Gaussian distribution with variance C_l. Even if we are able to measure the a_lm perfectly, our ability to correctly estimate the variance of the distribution will be affected by the number of modes for each value of l. Therefore, the observed power spectrum ${\hat{C}}_{l}$ will not be of the same form as the expectation value C_l.

In contrast, for a single source, the choice of each a_lm depends on the sky location and polarization angle of the source, but is otherwise entirely set. The distribution is entirely random for a Gaussian field, and entirely non-random for a single source.

This limitation on the measurability of C_l is the cosmic variance familiar from calculations of the CMB power spectrum:

$$\begin{eqnarray}&&{\rm{\Delta }}{C}_{l}=\displaystyle \frac{{C}_{l}}{\sqrt{2l+1}}.\end{eqnarray} \tag{ 22 }$$

This equation differs from the standard definition (e.g., Dodelson 2003) by a factor of $\sqrt{2}$. Since the Fourier-transformed GWB is a complex field, it contains twice the information of a real field such as the CMB.

The cosmic variance represents the range within which an observed ${\hat{C}}_{l}$ is expected to differ from the true unobservable C_l for each l. When ${\hat{C}}_{l}$ no longer follows Equation (17), $\hat{C}(\theta )$ will no longer follow the Hellings & Downs curve. By changing a single multipole, we are effectively changing the weights of individual Legendre polynomials in Equation (19).

Since C_l, and consequently ΔC_l, is a strong function of l, the effect of cosmic variance will be strongest for the first few multipoles. This is shown in Figure 6. The quadrupole term is by far the most important, but combinations of several other terms can also affect the shape of the two-point correlation function.

Zoom In Zoom Out Reset image size

It will often be the case that the power spectrum of a Gaussian random GWB will have a low or high quadrupole or octopole by chance. It would therefore not be surprising to have a two-point correlation function which does not match the Hellings & Downs curve. Note that in contrast to the work in Mingarelli et al. (2013), Taylor & Gair (2013), and Gair et al. (2014), this change in shape of the two-point correlation function is not due to large-scale anisotropy in the source population, but occurs even in statistically isotropic GWBs.

For two sources which induce a noticeable shape shift on the two-point correlation function, the basic mechanism is the same as for the Gaussian case: specific Legendre polynomials in the expansion of Equation (19) are being up- or downweighted. The primary difference in these two cases is that for a Gaussian field, the amount by which each C_l varies from the expectation value is independent of all the others. For two sources, the specific interference pattern between the two redshift maps leads to oscillations in the power spectrum. These oscillations are set by the relative orientation and distance of the sources and are not random and not independent.

So far we have been discussing a GWB composed of a single frequency bin. However, PTAs are typically sensitive to a range of frequencies. If all the frequency bins under consideration can be described by the same underlying distribution, the effect of cosmic variance can be ameliorated. This is because the separate frequency bins can be considered as independent realizations of the same map—including more frequency bins allows us to sample the distributions more accurately. From Equation (22), ΔC_l ∼ (N_modes)^−1/2, so using n similar bins in the analysis will reduce the cosmic variance by a factor of $\sqrt{n}$. Even for a single frequency bin, an analysis assuming Hellings & Downs behavior may be sufficient to allow an initial detection (Cornish & Sampson 2016).

Our work has focused on the analysis of redshift maps in harmonic space, inspired by CMB analyses. However, unlike CMB experiments which make measurements over large regions of the sky, PTAs are only sensitive to the redshift field in the direction of its pulsars. The observed field is a partial sky map, sampled at M discrete sky locations corresponding to the positions of the pulsars in the array.

The number of pulsars will limit the degree to which a PTA can measure harmonics of the GWB, but the steepness of the power spectrum will turn out to be a more important limitation. We estimate PTA sensitivity to the power spectrum through the following signal-to-noise (S/N) calculation.

For a sparse sampling of the sky, an estimate of a spherical harmonic expansion of a field r(θ, ϕ) sampled at points (θ_i, ϕ_i) can be constructed as

$$\begin{eqnarray}&&{\hat{a}}_{{lm}}=\displaystyle \frac{4\pi }{{\displaystyle \sum }_{i}{w}_{i}}\displaystyle \sum _{i}{w}_{i}{Y}_{{lm}}({\theta }_{i},{\phi }_{i})r({\theta }_{i},{\phi }_{i}),\end{eqnarray} \tag{ 23 }$$

where w_i are weights that can be tuned for each point. The minimum variance estimate will have w_i equal to $1/{\sigma }_{i}^{2}$, where σ_i is the variance at each point. The formally optimal solution would have σ_i only including detector noise and terms intrinsic to the pulsar. However, in practice it would be difficult to separate a given pulsar's noise properties from a GWB.

For simplicity, we assume that all pulsars have equal weight, with rms noise of each pulsar (for GWs plus noise) σ₀. Generalizing to varying noise levels is straightforward. For M pulsars, the estimated angular power spectrum will then have a noise bias

$$\begin{eqnarray}&&{C}_{l}^{N}=\displaystyle \frac{4\pi }{M}{\sigma }_{0}^{2}.\end{eqnarray} \tag{ 24 }$$

We can use Equation (22) to estimate the S/N of the amplitude for a PTA, taking care to realize that the relevant C_ℓ for the noise estimate is the combined signal and noise power spectrum. For the signal power spectrum, we know that the form should follow from Equation (21). Given a variance in residuals from GWs ${\sigma }_{\mathrm{gw}}^{2}$, we write

$$\begin{eqnarray}&&{C}_{l}^{S}={\sigma }_{\mathrm{gw}}^{2}{C}_{l}.\end{eqnarray} \tag{ 25 }$$

The resulting estimate for the S/N for a given multipole l is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\rm{S}}}{{\rm{N}}}\right)}_{l}^{2}={({C}_{l}^{S})}^{2}\displaystyle \frac{2l+1}{{({C}_{l}^{S}+{C}_{l}^{N})}^{2}}.\end{eqnarray} \tag{ 26 }$$

For a first detection, the expectation is that the noise power in the large-scale correlated timing residuals will be much larger than the signal power. We can then simplify the S/N estimate by dropping the signal part of the last term. Explicitly, the S/N in the limit of a weak detection is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\rm{S}}}{{\rm{N}}}\right)}_{l;\mathrm{weak}}^{2}=\displaystyle \frac{9}{4}\displaystyle \frac{2l+1}{{[(l-1)l(l+1)(l+2)]}^{2}}\displaystyle \frac{{M}^{2}{\sigma }_{\mathrm{gw}}^{4}}{{\sigma }_{0}^{4}}.\end{eqnarray} \tag{ 27 }$$

Summing this over all ℓ gives a numerical prefactor of 1/48, with the l = 2 term alone contributing 5/256. The l = 2 term thus contributes 93.75% of the (S/N)². If one only measured the power in the quadrupole anisotropy of the timing residuals, the resulting S/N would be 97% of the total S/N available. This is simply because the quadrupole is contributing such a large fraction of the total power that it is far and away the largest signal to be measured and the S/N adds in quadrature rather than linearly.

To compare with previous work, we can do the similar calculation in map space. As shown in Siemens et al. (2013), the comparable prefactor for this calculation in map space reduces to the total number of pairs times the mean of the square of the Hellings & Downs curve. For a full-sky survey, the mean of the square of the Hellings & Downs curve is 1/48, while the number of unique pulsar pairs is $M(M-1)/2$, very close to M²/2, with the $(M-1)$ instead of M coming from the explicit nulling of autocorrelations in the calculation.

In this work we have introduced an alternative framework for considering spatial variation in GWBs. We primarily work with the all-sky redshift patterns induced by GWs passing the Earth. Using standard techniques from CMB analysis, we do all calculations in harmonic space for computational simplicity, but convert to map space to discuss measurable quantities. Since we assume non-evolving GW sources, all results are also true for the Earth term of the expected pulsar timing residuals, up to a normalization. This assumption breaks down for rare high-mass, high-frequency binaries which evolve on timescales of ∼kyr rather than ∼Myr (Mingarelli et al. 2012).

We explicitly decomposed the redshift pattern produced by a single source of GWs into spherical harmonics, which allowed us to calculate the power spectrum of a single source's redshift map exactly. We showed that the expectation value of the power spectrum for a statistically isotropic Gaussian random GWB has the same form as for a single source. Using the relation between the power spectrum and the real space two-point correlation function, we explored the degree to which variance in the power spectrum changes the shape of the two-point correlation function away from Hellings & Downs. In particular, cosmic variance in the quadrupole moment of the power spectrum for a Gaussian random field can have significant effects on the amplitude of the curve, while also changing its shape. Finally, we showed that the quadrupole term of the power spectrum contributes 97% of the S/N measured by a PTA.

Throughout this work, we have treated the GWB as one of two idealized cases: a single source or a Gaussian random field. A GWB produced by a population of sources will lie somewhere between these two cases, as suggested by Figure 2. It is likely that the degree to which a population of sources resembles one case or another changes as a function of frequency, with shot noise in the SMBBH population becoming more important at higher frequencies.

We have confirmed that GWBs dominated by a single bright source, which are highly anisotropic and non-Gaussian, and those which are isotropic, unpolarized, and Gaussian look very similar from the point of view of a two-point correlation function, as previously reported by Cornish & Sesana (2013). This suggests that two-point correlation functions will be effective for detecting GWBs of all kinds. But they will be ineffective for characterizing GWBs and searching for single sources, despite the clear visual differences between Figures 3 and 5.

This difference should be measurable given a sufficiently high significance measurement of the GWB and some luck in its orientation with respect to low-noise pulsars. A particularly clear example is given in Boyle & Pen (2012): consider the timing residuals for four pulsars, each of which is located in a different stripe near the top of the map in Figure 3. The GW signal in each pulsar will be perfectly correlated, differing only by a phase factor of 180° between adjacent stripes. No such perfect (anti-)correlation is possible for nearby pulsars affected by a Gaussian field such as in Figure 5.

An important difference between these two types of GWBs is that the redshift map produced by a single source is highly non-Gaussian. Although symmetric Gaussian distributions can be statistically completely described by their two-point functions, non-Gaussian distributions may have higher moments. Indeed, the example given by Boyle & Pen (2012) is a kind of four-point function. Future work will explore higher-order correlation functions as a means of characterizing the degree to which a GWB has Gaussian or point-source-like characteristics.

We thank Vicky Kaspi, Chiara Mingarelli, and Pat Scott for useful discussions and comments. We acknowledge the support of NSERC and the Canadian Institute for Advanced Research and resources provided by Calcul Québec.

Software: healpy (Górski et al. 2005), matplotlib (Hunter 2007), numpy, pandas.

From Equation (13), we have the redshift induced in a direction (θ, ϕ) by a source located at the north pole. This is a complex scalar field, and can be expanded in spherical harmonics with coefficients:

$$\begin{eqnarray}&&{a}_{{lm}}=\int d{\rm{\Omega }}\left[\displaystyle \frac{1}{2}(1+\cos \theta )(\cos 2\phi \,{h}_{+}-\sin 2\phi \,{h}_{\times })\right]{Y}_{{lm}}^{* }\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&=\mathop{\underbrace{\displaystyle \frac{1}{2}{(-1)}^{m}\sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l-m)!}{(l+m)!}}}}\limits_{{\rm{A}}}\\ &&\quad \times \mathop{\overbrace{{\displaystyle \int }_{0}^{2\pi }d\phi \,{e}^{-{im}\phi }(\cos 2\phi \,{\tilde{h}}_{+}-\sin 2\phi \,{\tilde{h}}_{\times })}}\limits^{{\rm{B}}}\\ &&\quad \times \mathop{\underbrace{{\displaystyle \int }_{0}^{\pi }d\theta \,\sin \theta \,(1+\cos \theta ){P}_{l}^{m}(\cos \theta )}}\limits_{{\rm{C}}},\end{eqnarray} \tag{ 29 }$$

where ${P}_{l}^{m}(\cos \theta )$ are the associated Legendre polynomials. This factorizes into two integrals (B, C) and one constant term (A). Beginning with the integral over ϕ, we find that it simplifies to

$$\begin{eqnarray}B & = & \displaystyle \frac{1}{2}({\tilde{h}}_{+}+i{\tilde{h}}_{\times }){\displaystyle \int }_{0}^{2\pi }d\phi \,{e}^{-i\phi (m-2)}\\ & & +\displaystyle \frac{1}{2}({\tilde{h}}_{+}-i{\tilde{h}}_{\times }){\displaystyle \int }_{0}^{2\pi }d\phi \,{e}^{-i\phi (m+2)}.\end{eqnarray} \tag{ 30 }$$

Since ${\int }_{0}^{2\pi }d\phi \,\exp (-\alpha i\phi )=2\pi \,\delta (\alpha )$ for real α, only terms with m = ±2 exist. They are given by

$$\begin{eqnarray}B=\left\{\begin{array}{ll}\pi ({\tilde{h}}_{+}\pm i{\tilde{h}}_{\times }), & m=\pm 2\\ 0, & \mathrm{all}\ \mathrm{other}\ m.\end{array}\right.\end{eqnarray} \tag{ 31 }$$

This constraint on m allows us to simplify both our constant term and θ integral:

$$\begin{eqnarray}A\cdot C & = & \displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l\mp 2)!}{(l\pm 2)!}}{\displaystyle \int }_{0}^{\pi }d\theta \,\sin \theta \\ & & \times \,(1+\cos \theta ){P}_{l}^{\pm 2}(\cos \theta )\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l-2)!}{(l+2)!}}{\int }_{0}^{\pi }d\theta \,\sin \theta (1+\cos \theta ){P}_{l}^{2}(\cos \theta ),\end{eqnarray} \tag{ 33 }$$

using the following property of associated Legendre polynomials:

$$\begin{eqnarray}&&{P}_{l}^{-m}={(-1)}^{m}\displaystyle \frac{(l-m)!}{(l+m)!}{P}_{l}^{m}.\end{eqnarray} \tag{ 34 }$$

Writing $\mu =\cos \theta$, we solve the integral over θ:

$$\begin{eqnarray}&&C={\int }_{-1}^{1}d\mu \,(1+\mu ){P}_{l}^{2}(\mu )\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&={\displaystyle \int }_{-1}^{1}d\mu \,(1+\mu )(1-{\mu }^{2})\displaystyle \frac{{d}^{2}}{d{\mu }^{2}}{P}_{l}(\mu )\\ &&\quad (\mathrm{in}\ \mathrm{terms}\ \mathrm{of}\ \mathrm{ordinary}\ \mathrm{Legendre}\ \mathrm{polynomials})\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&={\displaystyle \int }_{-1}^{1}d\mu \,(1+\mu )\left[2\mu \displaystyle \frac{d}{d\mu }{P}_{l}(\mu )-l(l+1){P}_{l}(\mu )\right]\\ &&\quad (\mathrm{using}\ \mathrm{the}\ \mathrm{defining}\ \mathrm{differential}\ \mathrm{equation})\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&={\displaystyle \int }_{-1}^{1}d\mu \,(1+\mu )2\mu \displaystyle \frac{d}{d\mu }{P}_{l}(\mu )\\ &&\quad (\mathrm{second}\ \mathrm{term}\ \mathrm{zero}\ \mathrm{by}\ \mathrm{orthogonality})\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&={[2{\mu }^{2}{P}_{l}(\mu )+2\mu {P}_{l}(\mu )]}_{-1}^{1}-{\int }_{-1}^{1}d\mu \,(4\mu +2){P}_{l}(\mu )\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&=2[{P}_{l}(1)+{P}_{l}(1)+{P}_{l}(-1)-{P}_{l}(-1)]\\ &&\quad (\mathrm{second}\ \mathrm{term}\ \mathrm{zero}\ \mathrm{by}\ \mathrm{orthogonality})\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&=\,4.\end{eqnarray} \tag{ 41 }$$

Putting everything together,

$$\begin{eqnarray}&&{a}_{{lm}}=A\cdot B\cdot C=\left\{\begin{array}{l}2\pi \sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l-2)!}{(l+2)!}}({\tilde{h}}_{+}\pm i{\tilde{h}}_{\times }),\\ \quad m=\pm 2,l\geqslant 2.\\ 0,\quad \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 42 }$$