# Characterising gravitational wave stochastic background anisotropy with Pulsar Timing Arrays

## Abstract

Detecting a stochastic gravitational wave background, particularly radiation from individually unresolvable supermassive black hole binary systems, is one of the primary targets for Pulsar Timing Arrays. Increasingly more stringent upper limits are being set on these signals under the assumption that the background radiation is isotropic. However, some level of anisotropy may be present and the characterisation of the gravitational wave energy density at different angular scales carries important information. We show that the standard analysis for isotropic backgrounds can be generalised in a conceptually straightforward way to the case of generic anisotropic background radiation by decomposing the angular distribution of the gravitational wave energy density on the sky into multipole moments. We introduce the concept of generalised overlap reduction functions which characterise the effect of the anisotropy multipoles on the correlation of the timing residuals from the pulsars timed by a Pulsar Timing Array. In a search for a signal characterised by a generic anisotropy, the generalised overlap reduction functions play the role of the so-called Hellings and Downs curve used for isotropic radiation. We compute the generalised overlap reduction functions for a generic level of anisotropy and Pulsar Timing Array configuration. We also provide an order of magnitude estimate of the level of anisotropy that can be expected in the background generated by super-massive black hole binary systems.

###### pacs:

04.80.Nn, 04.25.dg, 95.85.Sz, 97.80.-d 97.60.Gb 04.30.-w## I Introduction

The detection of gravitational waves (GWs) is one of the key scientific goals of Pulsar Timing Arrays (PTAs). A PTA uses a network of radio telescopes to regularly monitor stable millisecond pulsars, constituting a galactic-scale GW detector (1); (2); (3); (4). Gravitational radiation affects the propagation of radio pulses between a pulsar and a telescope at the Earth. The difference between the expected and actual time-of-arrival (TOA) of the pulses – the so-called timing residuals – carries information about the GWs (5); (6); (7), which can be extracted by correlating the residuals from different pulsar pairs. This type of GW detector is sensitive to radiation in the Hz frequency band, a portion of the spectrum in which a promising class of sources are super-massive black hole binary (SMBHB) systems with masses in the range of during their slow, adiabatic in-spiral phase (8); (9); (10); (11); (12); (13); (14). Other forms of radiation could be observable by PTAs, such as cosmic strings (21); (15); (16) and/or a background produced by other speculative processes in the early universe, see e.g. (22).

A PTA can be thought of as an all-sky monitor that is sensitive to radiation from the whole cosmic population of SMBHBs radiating in the relevant frequency band. The overwhelming majority of sources are individually unresolvable, but the incoherent superposition of the very weak radiation from the many binaries in the population gives rise to a stochastic background^{1}

In all the searches carried out so far, it has been assumed that the stochastic background, regardless of its origin, is isotropic (25); (26); (27); (28). This is well justified if the background is produced by some physical processes in the early universe or is largely dominated by high-redshift sources. Under the assumption of isotropy, the correlated output from the data from any two pulsars in the array depends only on the angular separation of the pulsars and is known as the Hellings and Downs curve (25). However, a PTA also carries information about the angular distribution of the GW power on the sky. It is therefore important to address how this information is encoded in the data, and the implications for analysis approaches. In fact, if evidence for a signal is found in the data, testing the assumption of isotropy could be one of the methods to confirm its cosmological origin. If, on the other hand, one expects some deviations from isotropy, which may be the case for the SMBHB background created by a finite population (36); (37), it is useful to be able to extract constraints on the underlying physical population.

In this paper we show how the correlated output from pulsar pairs in a PTA is related to the anisotropy of the signal, i.e. the angular distribution of GW power on the sky, and how one can extract this information by measuring the multipole moments that characterise the anisotropy level, following an analogous approach to those applied to the case of ground-based (38) and space-based (39) laser interferometric observations. By doing this, we generalise the Hellings and Downs curve to an arbitrary angular distribution on the sky. We also provide an estimate for the expected level of anisotropy for the background produced by an arbitrary population of sources, and in particular, the population of SMBHB systems.

The paper is organised as follows. In Section II, we review the basic concepts of a GW stochastic background, and we estimate the expected level of anisotropy in a background produced by a population of SMBHB systems. We show that at low frequencies, where the PTA sensitivity is optimal and the number of sources that contribute to the background is very large, the expected level of anisotropy is small, and likely undetectable. However towards the high-frequency end of the sensitivity window, where the actual number of sources decreases sharply, the anisotropy level could be significant, increasing at smaller angular scales. In Section III we show that the present analysis approaches for isotropic backgrounds can be generalised in a conceptually straightforward way to the case of anisotropic signals by decomposing the angular distribution of the GW power on the sky into multipole moments. We introduce the concept of generalised overlap reduction functions, which replace the Hellings and Downs curve. Each one of these characterises the effect of a given anisotropy multipole on the correlation of the timing residuals from a pulsar pair. In Section IV we derive expressions for the generalised overlap reduction functions for an arbitrary stochastic background angular distribution on the sky and PTA configuration. This is essential for future analyses of PTA data which include anisotropy as part of the model. Section V contains our conclusions and suggestions for future work.

For the rest of the paper we will consider geometric units, and therefore set , unless otherwise specified.

## Ii Stochastic backgrounds

Let us consider a plane wave expansion for the metric perturbation produced by a stochastic background:

(1) |

where is the frequency of the GWs, the index labels the two independent polarisations, the spatial indices are , the integral is on the two-sphere , and our sign convention for the Fourier transform of a generic function follows the GW literature convention

(2) |

The unit vector identifies the propagation direction of a single gravitational wave plane, that can be decomposed over the GW polarisation tensors and the two independent polarisation amplitudes, or equivalently (40); (41):

(3a) | ||||

(3b) |

The polarisation tensors are uniquely defined once one specifies the wave principal axes described by the unit vectors and :

(4a) | ||||

(4b) |

For a stationary, Gaussian and unpolarised background the polarisation amplitudes satisfy the following statistical properties:

(5) |

where is the expectation value and is the covariant Dirac delta function on the two-sphere (42). This condition implies that the radiation from different directions are statistically independent. Moreover, we have factorised the power spectrum such that , where the function describes the spectral content of the radiation, and describes the angular distribution on the sky.

The mass-energy density in GWs is (41)

(6) |

and using Eqs. (1) and (5) we have:

(7) |

The GW energy density, Eq. (6) is related to the more frequently used density parameter by:

(8) |

where is the GW energy density in the infinitesimal band from to , is the critical density at the present epoch, and is the value of the Hubble parameter today.

Using Eqs. (1), (5) and (6) we can rewrite Eq. (8) as

(9) |

which shows that the energy-density content of the background is the result of contributions from all the directions . Each direction on the sky need not contribute to the background in the same way, and the function describes the angular dependence (the “hot” and “cold” spots). As in (38), we decompose the angular distribution function on the basis of the spherical harmonic functions,

(10) |

where the sum is over , and . The coefficients are the multipole moments of the radiation which characterise the angular distribution of the background. We adopt the convention that the monopole moment is normalised as

(11) |

which yields

(12) |

Eq. (9) now becomes

(13) |

The frequency spectrum of the background, whether from SMBHBs or other sources or processes in the early Universe, is described by the function . The angular distribution of the radiation is encoded in the values of the radiation multipole moments , which become unknown parameters in the analysis. In Section 3 and 4 we will show how the ’s enter the likelihood function of PTA timing residuals, and how an arbitrary angular distribution affects the correlation of radiation at any two pulsars timed by an array. This provides a way of measuring the multipole moments. In the remainder of this Section we provide an estimate of the expected level of anisotropy in a background generated by the population of SMBHB systems.

In order to gain some insight into this problem, let us consider an idealised situation, constructed as follows. Let us assume that the universe is populated by identical sources with number density . If we want to estimate the level of anisotropy, we need to estimate the expected value of the energy density in GWs coming from sources in a solid angle centred on a direction and compare it to the energy density produced by sources in a cone centred on a different direction . For this example we consider a Euclidean, static universe (or equivalently sufficiently nearby sources, such that we do not take into account effects of expansion and redshift).

In a conical volume within the solid angle and at distance between and , the expected number of sources which contribute to the background is:

(14) |

The actual number of sources is then governed by Poisson statistics, with mean and variance . If the volume is sufficiently small that , then the probability of finding one source is

(15) |

Since the probability of having more than one source within this volume is negligible, the probability of finding no sources is simply .

The expected total number of sources, , present in the whole volume within a solid angle between the minimum and maximum distance, and , respectively (to be discussed later), is given by the sum of the contributions from each slice in the cone. Similarly, the variance is the sum of the variances from each conical slice. We therefore obtain

(16a) | |||||

(16b) |

We now want to compute the expected contribution to the GW energy density per frequency interval and its variance. The GW energy density of each source scales as . If we assume that all the sources are identical – the generalisation to a distribution of masses is straightforward, but is not needed to explain the key points – we can write (with slight abuse of notation) the contribution to the energy density per source simply as

(17) |

where is an appropriate constant factor, equal for all sources.

The expected GW energy density from sources in a small conical volume at distance , again chosen so that it has a vanishingly small probability of having more than one source, , see Eqs. (14) and (15), is

(18) |

The variance of the energy density from sources in this conical volume is

(19) |

where the last equality relies on the consistent application of the condition (which can always be satisfied by choosing a sufficiently small shell thickness ).

We can now compute the expected contribution to the GW energy density and its variance from all sources in a solid angle . The mean energy density and variance are given by the sum of contributions from all slices of thickness ; using Eqs. (18) and (19), this yields:

(20a) | |||||

(20b) | |||||

(20c) |

and, using the fact that the variance of a sum is the sum of variances,

(21a) | |||||

(21b) | |||||

(21c) |

We define the level of anisotropy as the ratio of the standard deviation in the GW power emanating from a given solid angle to the expected power from that angle:

(22) | |||||

We can now return to the choice of the minimal and maximal distance, and . The maximal distance at which sources can be located is set by cosmology and the history of SMBH formation. Meanwhile, the minimal distance of interest to us, , corresponds to the maximal distance at which individual binaries can be resolved. Individually resolvable binaries can be subtracted from the data, and are treated separately from the stochastic background. An individual source can be efficiently searched for with matched filtering techniques (32); (33); (34); (35). Therefore, we expect the power necessary to detect a single SMBH binary to be significantly less than the power necessary to measure a stochastic background. Thus, in order for a stochastic background to be detectable after all individual sources that are presumed to be detectable up to distance are removed, the total power in the background must be significantly greater than the power in the weakest individually resolvable source:

(23) |

Another way to interpret the preceding condition is to consider the idealised situation when the stochastic background provides the dominant noise source: optimal matched filtering would make it possible to individually resolve and subtract coalescing SMBH binaries with signal power far below the noise (background) levels.

We can recast the condition on the detectability of a stochastic background, Eq. (23), as

(24) |

If we define , where , this condition yields

(25) |

where is the total number of sources in the universe, modulo a factor of order unity. We can now rewrite the level of anisotropy (22) in the following form:

(26a) | |||||

(26b) |

where is the total number of sources that contribute to the background and . Note that by virtue of condition (25), the second term in Eq.(26a) is always smaller than unity whenever the stochastic background is detectable, and is actually . The level of anisotropy scales as , and increases by going to small angular scales . However, there is an observational limit on the angular resolution of PTAs which will prevent very small angular scales from being probed. Furthermore, at smaller angular scales, the signal will be progressively dominated by a smaller number of, possibly individually unresolvable, sources. The number of sources in a cone of solid angle is

(27a) | |||||

(27b) | |||||

(27c) |

When this quantity is larger but not much larger than unity, we expect to be in the middle ground between searches for individual sources and standard stochastic-background searches. If this occurs on resolvable angular scales where anisotropy is significant (cf. Eq. (26a) and Eq. (29) below), it will be interesting to check the efficiency of current search pipelines in this regime.

Using the results from e.g. (11) we can provide an order-of-magnitude estimate of the expected level of anisotropy that characterises the SMBHB background. From Figure 4 of Ref. (11) we can see that the total number of sources that contribute in a frequency interval of width , where is the observation time, can be approximated as:

(28) |

where we used the fact that, during a SMBHB inspiral, the time the binary spends in a given frequency band scales as . Substituting Eq. (28) into Eq. (26b) and converting between the average angular scale and the multipole moment index using , we obtain:

(29) | |||||

There will be few SMBHBs beyond redshift , and individual sources are likely to be resolvable up to redshift , so sources that contribute to the stochastic background are within redshift range –, see e.g. (11); (12). Therefore, both and will be factors of order unity. We have confirmed this with a more careful calculation that takes cosmology and the redshifting of gravitational waves into account; however, we note that our simplified treatment relied on a constant density (rate) of coalescing SMBHBs in the Universe, and on a fixed amplitude at a given frequency for all sources, which corresponds to the assumption of a fixed source mass.

As expected, the level of anisotropy at low frequencies and large angular scales is small. However, it can become non-negligible, at the tens of percent level, at frequencies Hz.

## Iii Effect of anisotropy on timing residuals

GWs affect the time of arrival at the telescope of radio pulses from ultra-stable pulsars. Consider a pulsar with frequency whose location in the sky is described by the unit vector . The pulsar is at a distance from the Earth. The effect of a GW source in the direction generating a metric perturbation is to affect the actual frequency at which the radio pulses are received at a telescope, according to

(30) | |||||

where

(31) |

is the difference between the metric perturbation at the Earth , the so-called Earth term, with coordinates , and at the pulsar , the so-called pulsar term, with coordinates .^{2}

(32a) | ||||

(32b) |

where the indices “e” and “p” refer to the Earth and the pulsar. In this frame we can therefore write Eq (31) using Eq (3b)

The fractional frequency shift produced by a stochastic background is simply given by integrating Eq. (30) over all directions. Using Eq (LABEL:e:deltah1), we obtain:

(34) | |||||

where are the antenna beam patterns for each polarisation , defined as

(35) |

The quantity that is actually observed is the time-residual , which is simply the time integral of Eq. (34):

(36) |

The search for a stochastic background contribution in PTA data relies on looking for correlations induced by GWs in the residuals from different pulsars. These correlations in turn depend on the spectrum of the radiation, cf. Eqs (9) and (5), and the antenna beam pattern convolved with the angular distribution of the GW energy density in the sky, cf. Eq. (10).

Regardless of whether the analysis is carried out in a frequentist framework, and therefore one considers a detection statistic, see e.g. Ref (44), or one builds a Bayesian analysis e.g. (45), the key physical quantity that is exploited is the correlation of the timing residuals for every pair of pulsars timed by a PTA. In a frequentist analysis this enters into a suitable detection statistic, whereas in a Bayesian framework it enters the likelihood function. The expected value of the correlation between a residual from a pulsar, say , at time , with that from a different pulsar, say , at time depends on terms of the form:

(37) | |||||

In analogy with Ref. (41), we define the quantity in the previous equation that depends on the relative location of the pulsars in the PTA, and the angular distribution of the GW energy density as the overlap reduction function:

(38) |

where

(39) |

In Eq. (37) contains the information of the spectrum of radiation, and contains information about the angular distribution of GW background power. Under the assumption that the background is isotropic, is a known function that simply depends on the location of the pulsars timed by the array since is constant. In this case, the overlap reduction function (38) becomes:

(40) |

which is the result derived by Hellings and Downs in Ref. (25) and is known (up to a normalisation constant) as the Hellings and Downs curve.

For an anisotropic background, whose angular power spectrum is unknown, is a function of the unknown angular power distribution on the sky. We can generalise the concept of the overlap reduction function by decomposing on the basis of spherical harmonic functions according to Eq. (10). The weight of each of the components is given by an unknown coefficient , which needs to be determined by the analysis. The overlap reduction function (38) becomes therefore:

(41) |

where

(42) |

are the (complex-form) generalised overlap reduction functions. Given an array of pulsars on the sky, the functions are uniquely defined and known.

The generalisation of e.g. the standard Bayesian analysis for an isotropic stochastic background such as the one reported in Ref. (45) to the case in which the assumption of isotropy is relaxed is, at least conceptually, straightforward. The model parameters that describe the stochastic background are not only those that enter the frequency spectrum – for example the overall level and spectral index in the common case of a power-law parametrisation of , appropriate for the background from SMBHBs – but also the coefficients that describe the angular distribution on the sky, that is, how much power is associated to each spherical harmonic decomposition of the overall signal. An initial implementation of this analysis is reported in Ref. (53).

Before we compute the expressions for the generalised overlap reduction functions, it is important to consider the function , defined in Eq. (39) and present in Eqs. (38) and (42), which introduces the frequency dependence of the overlap reduction functions. From a physical point of view encodes the fact that the correlation of the timing residuals carries information about both the Earth and pulsar terms for the two pulsars whose timing residuals are correlated. The relevant scale in the function is

(43) |

which introduces rapid oscillations around unity (44) that depend on the distance and location to the pulsars. For all astrophysically relevant situations , see Eq. (43), and when one computes the integral in Eq. (42) the frequency dependent contributions to the integral rapidly average out to zero as the angle between the pulsar pairs, , increases. The generalised overlap reduction function Eq. (42) is therefore well approximated by

(44) |

where is the Kronecker delta. We will provide some more details in Section IV.3. Here we note that the approximation (44) is equivalent to considering only the correlation of the Earth-term for two distinct pulsars. As we are considering many sources over the whole sky then the pulsar terms will only contribute to the correlation if the distance between two pulsars is of the order of one wavelength or less, and for the frequencies and pulsars being considered this is only true for auto-correlation. The auto-correlation term carries contributions from the Earth and pulsar terms, and therefore the value of of the integral is multiplied by a factor of 2. Note also, that the generalised overlap reduction function (44) does not depend on frequency.

The decompositions (41), (42) and (44) are based on the usual complex-basis spherical harmonic functions , whose definitions are given in Section IV.2. One can alternatively consider a decomposition on a real basis , that are defined in terms of their complex analogs by^{3}

(45) |

Consequently, the real-form generalised overlap reduction functions are:

(46) |

In the next Section we compute the ’s for a generic pulsar pair and discuss their properties.

## Iv Generalised overlap reduction functions

In this Section we compute the generalised overlap-reduction functions, Eq. (44) for a generic pulsar pair and explore their properties. Anholm et al. (44) considered the particular case of the overlap-reduction function between two pulsars for radiation described by dipole anisotropy. Here we go beyond, and consider an arbitrary angular distribution of the background. Our approach is based on decomposing the power of the GW background at different angular scales onto spherical harmonics, cf. Eq. (10) and for the specific case of a dipole distribution we show that our result is equivalent to the one presented in (44).

In the case of an isotropic background, pulsar pairs timed by a
PTA map uniquely into the Hellings and Downs curve. That is to say,
any pulsar pair is uniquely identified by an angular separation, which
in turn corresponds to a value of the overlap reduction
function. This is no longer the case for an anisotropic
distribution. For a given distribution of the GW power on the sky, the
generalised overlap reduction functions depend on the angular
separation between two pulsars and their specific location in
the sky with respect to the background radiation. Equivalently, if one
considers two different pulsar pairs with the same angular separation
but different sky locations, the overlap reduction function that
describes the correlation between the two pulsars will be
different. To illustrate this, we show a selection of the best pulsars currently being timed by the European
Pulsar Timing Array (EPTA) (48) ^{4}

In our analysis we will closely follow the approach considered by Allen and Ottewill (38), who considered the equivalent problem in the case of ground-based laser interferometers.

### iv.1 Choice of coordinate frames

We introduce a “cosmic rest-frame” where the angular dependency of the anisotropy is described, and a “computational frame”, in which some of the key expressions take a particularly simple form, and provide some intuitive clues into the problem. Given any two pulsars, say P and P, we define the computational frame as the frame in which pulsar P is on the z-axis, pulsar P is in the plane, and their angular separation is denoted by . This is the standard frame that is used in e.g. (44) to compute the Hellings and Downs curve for the isotropic case. Therefore, overlap reduction functions in the computational frame only depend on the pulsar pair’s angular separation, . We now outline a method where one can rotate from the cosmic rest-frame to the computational frame, and vice versa, by means of rotation matrices.

Let us consider a generic vector , and let (unprimed) be the component in the cosmic rest-frame and (primed) the component in the computational frame, which will be different for every pulsar pair. The components of the vector in the two different frames are related by:

(47) | |||||

where is the rotation matrix given by:

(48) | |||

Indeed, we must carry out three rotations to go from the cosmic rest-frame to the computational frame. If the pulsars P and P in the cosmic rest-frame have polar coordinates ) and ), respectively, the three angles of the rotations are:

(49a) | ||||

(49b) | ||||

(49c) |

The condition on has two solutions within the range [0,2] and we choose the one that gives a positive coordinate in the computational frame for P.

Having calculated the relevant angles we can apply these to the rotation of spherical harmonics, where we know from Eq. (4.260) in Ref. (50):

(50) |

and

(51) |

where equations (50) and (51) rotate from the computational frame into the cosmic rest-frame, and back to the computational frame, respectively. The matrix is given by Eq. (4.12) in (49)

(52) |

and for

(53) |

where is the hypergeometric Gaussian function. For , can be derived from the unitary property, and yields

(54) |

as in Eq. (4.15) in Ref. (49). We also note that the ’s are real. Since in Eq. (42) is a function of , we can now write the generalised overlap reduction function in the cosmic rest-frame as

(55) |

where (primed) is the generalised overlap reduction function in the computational frame.

### iv.2 Generalised overlap reduction functions in the computational frame

In order to compute the generalised overlap reduction function in the cosmic rest-frame, Eq. (42) or (46), one needs to compute the relevant function in the computational frame then rotate it via Eq. (55) using the matrix (52). Here we compute the generalised overlap reduction functions in the computational frame. For ease of notation, we drop the primes, but it understood that in this section all the analysis is done in the primed, computational frame.

The spherical harmonic function of order and degree , is defined as

(56) | |||||

(57) |

where is the azimuthal angle and is the polar angle and the are the associated Legendre polynomials

(58a) | |||||

(58b) |

and

(59) |

is the normalisation. The Hellings and Downs curve – or equivalently the overlap reduction function for an isotropic background – can be derived (up to a normalisation constant) setting , i.e. .

For each pair of pulsars, the computational frame is defined by the following geometry:

(60a) | ||||

(60b) | ||||

(60c) | ||||

(60d) | ||||

(60e) |

where is the angular separation of the two pulsars, . In this frame , and Eq. (44) reduces to

(61) |

With this choice of frame, the generalised overlap reduction functions can be easily computed. It is worth pointing out that in this frame the ’s are real , and therefore since , where the star here denotes the complex conjugate. One then need only take into account the transformation properties of the associated Legendre polynomials defined in Eq. (58).

In Appendix A we provide comprehensive details of the derivations, whereas here we will just show the main results. For the case , Eq. (61), we obtain the overlap reduction function for the case of an isotropic background (cf. Appendix A.1):

(62) |

is the Hellings and Downs curve up to a multiplicative factor . In fact the Hellings and Downs curve is normalised in such a way that is unity when one considers the auto-correlation of the timing residuals form the same pulsar ( and therefore ). Note that for the isotropic case the rotation from the computational frame into the cosmic frame has no effect.

More generally, it is rather straightforward to compute analytical expressions for the case of a dipole () anisotropy. In this case the generalised overlap reduction functions in the computational frame read (cf. Appendix A.2):

(63a) | ||||

(63b) | ||||