0809/0809.0701.md

Optimal strategies for gravitational wave stochastic background searches in pulsar timing data

Melissa Anholm,^1,* Stefan Ballmer,^2,† Jolien D. E. Creighton,^1,‡ Larry R. Price,^1,§ and Xavier Siemens^1,¶

¹Center for Gravitation and Cosmology, Department of Physics,

University of Wisconsin–Milwaukee, P.O. Box 413, Milwaukee, Wisconsin 53201, USA

²LIGO Laboratory, California Institute of Technology, MS 18-34, Pasadena, California 91125, USA

(Dated: October 29, 2018)

A low frequency stochastic background of gravitational waves may be detected by pulsar timing experiments in the next five to ten years. Using methods developed to analyze interferometric gravitational wave data, in this paper we lay out the optimal techniques to detect a background of gravitational waves using a pulsar timing array. We show that for pulsar distances and gravitational wave frequencies typical of pulsar timing experiments, neglecting the effect of the metric perturbation at the pulsar does not result in a significant deviation from optimality. We discuss methods for setting upper limits using the optimal statistic, show how to construct skymaps using the pulsar timing array, and consider several issues associated with realistic analysis of pulsar timing data.

PACS numbers:

I. INTRODUCTION

The search for gravitational waves is at the forefront of current fundamental physics research. The direct detection of gravitational waves will usher in a new era in astronomy and astrophysics. Gravitational waves will reveal information about black holes, supernovae, and neutron stars that cannot be gleaned from electromagnetic observations. Furthermore, the detection of a gravitational wave background will open an observational window onto a time in the early universe before recombination, prior to which the universe is opaque to electromagnetic waves. The scientific rewards for such a detection would be truly exceptional. Several international efforts are underway to detect gravitational waves and two of these efforts are expected to result in the detection of gravitational waves in the next 5 to 10 years: Interferometric ground-based gravitational wave detectors and pulsar timing observations.

Neutron stars radiate powerful beams of radio waves from their magnetic poles. If a neutron star's magnetic poles are not aligned with its rotational axis, the beams sweep through space like the beacon on a lighthouse and the neutron star is said to be a pulsar. If the Earth lies within the sweep of a pulsar's beams, the star is observed as a point source in space emitting short, rapid bursts of radio waves [1]. Due to their enormous mass, neutron stars have a very large moment of inertia and the radio pulses we observe arrive at a very constant rate. Pulsar timing experiments exploit this regularity [2, 3, 4, 5]. Fluctuations in the time of arrival of radio pulses, af-

ter all known effects have been subtracted out, could be due to the presence of gravitational waves. Since the 1970s, when these ideas were first conceived, pulsar timing precision has improved dramatically. Several known pulsars can now be timed with a precision of about 1 micro-second, and a handful can be timed with a precision around several hundred nanoseconds [6]. Recent work [7] has shown that the presence of nanohertz gravitational waves could be detected by observing 20 pulsars with timing precisions of 100 nanoseconds over a period of 5 to 10 years. Non-detection would still improve current bounds on the low frequency stochastic gravitational wave background [8].

Gravitational waves from supermassive black hole binary systems could be detected via pulsar timing observations [9, 10, 11, 12, 13]. In addition, pulsar timing has the potential to measure the polarization properties of gravitational waves which could confirm (or even change) the current theory of gravity [14, 15]. Gravitational wave observations in the nanohertz band could also yield information about the early universe [16]. Cosmic strings, line-like topological defects, could produce gravitational waves in the nanohertz band. Cosmic strings can form during phase transitions in the early universe due to the rapid cooling that takes place after the Big Bang [17, 18, 19]. Cosmic string production is generic in supersymmetric grand unified theories [20]. Additionally, in string theory motivated cosmological models cosmic strings may also form (dubbed cosmic superstrings to differentiate them from field theoretic cosmic strings) [21, 22, 23, 24, 25, 26]. Cosmic strings and superstrings are expected to produce a background of stochastic gravitational waves and bursts of gravitational waves [27, 28, 29, 30, 31, 32] that could be detected using pulsar timing observations. Pulsar timing observations are already producing some of the most interesting constraints on cosmic string models and a detection would have profound implications [31, 32].

Recently Lommen [33] produced an upper limit on the

*Electronic address: anholm@gravity.phys.uwm.edu

†Electronic address: sballmer@caltech.edu

‡Electronic address: jolien@gravity.phys.uwm.edu

§Electronic address: larry@gravity.phys.uwm.edu

¶Electronic address: siemens@gravity.phys.uwm.edustochastic gravitational wave background using observations of three millisecond pulsars spanning 17 years. The methods used by Lommen were based on those developed by Kaspi and collaborators [34] and have been the subject of some criticism in the literature [8, 30]. More recently Jenet and collaborators [7, 8] developed a new technique for gravitational wave stochastic background searches in pulsar timing data and applied it to Parkes Pulsar Timing Array data [35, 36].

In this paper we consider optimal strategies for extraction of a gravitational wave stochastic background signal using data from a pulsar timing array. Our methods are based on those developed for and used in ground based interferometric gravitational-wave detectors such as LIGO and Virgo [37, 38, 39, 40, 41, 42, 43, 44, 45], and improve on existing methods in several ways. In Section II we write the redshift of pulsar signals induced by passing gravitational waves, first derived by Detweiler [3], in a coordinate-independent way more suitable for our analysis, and discuss its form in the frequency domain including the long-wavelength limit. In Section III we construct the optimal cross-correlation filter for a pulsar pair by maximizing the signal to noise. We find that the overlap reduction function is well approximated by a constant, or equivalently, that the metric perturbation at the pulsar can be neglected for values of pulsar distances and gravitational wave frequencies typical of pulsar timing experiments, without significant losses in sensitivity. In Section IV we also show how to construct the optimal combination of cross-correlations of pulsar pairs in a pulsar timing array and include a more sophisticated derivation of the optimal detection statistic based on the likelihood ratio. In Section V we discuss upper limit and detection methods. In Section VI we show how to construct sky maps using pulsar timing data—a pulsar timing radiometer. In Section VII we discuss several important issues relating to the realistic analysis of pulsar timing data, including the Lomb-Scargle periodogram for power spectrum estimation of unevenly sampled data, and optimal procedures for computing Fourier transforms. We conclude in Section VIII. Lommen, Romano and Woan [46] will extend our work using a likelihood based approach developed in [47], and consider the case of stochastic backgrounds that are loud compared to the noise, closely examine time-domain implementations of the optimal statistic, and provide a detailed comparison of the optimal statistic described here with the methods of Jenet and collaborators [7, 8].

II. THE SIGNAL

Gravitational waves affect pulsar timing measurements by creating perturbations in the null geodesics that the radio signals emitted from the pulsar travel on [3]. In this section we will describe the relationship between the metric perturbation and the signal measured in pulsar timing experiments.

A metric perturbation in a spatial, transverse, traceless gauge has a plane wave expansion given by [40]

$$h_{ij}(t,\vec{x})=\sum _A{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{\int }_{S^2}d\widehat{\mathrm{\Omega }}{e}^{i2\pi f(t-\widehat{\mathrm{\Omega }}\cdot \vec{x})}{h}_A(f,\widehat{\mathrm{\Omega }})e_{ij}^A(\widehat{\mathrm{\Omega }}),(1)h\_\{ij\}(t,\; \backslash vec\{x\})\; =\; \backslash sum\_A\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; e^\{i2\backslash pi\; f(t\; -\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash vec\{x\})\}\; h\_A(f,\; \backslash hat\{\backslash Omega\})\; e\_\{ij\}^A(\backslash hat\{\backslash Omega\}),\; \backslash quad\; (1)$$

where $f$ is the frequency of the gravitational waves, $\vec{k} = 2\pi f \hat{\Omega}$ is the wave vector, $\hat{\Omega}$ is a unit vector that points along the direction of travel of the waves, $i, j = x, y, z$ are spatial indices, and the index $A = +, \times$ labels polarizations. The polarization tensors $e_{ij}^A(\hat{\Omega})$ are

$$e_{ij}^{+}(\widehat{\mathrm{\Omega }})={\widehat{m}}_i{\widehat{m}}_j-{\widehat{n}}_i{\widehat{n}}_j,(2a)e\_\{ij\}^+(\backslash hat\{\backslash Omega\})\; =\; \backslash hat\{m\}\_i\; \backslash hat\{m\}\_j\; -\; \backslash hat\{n\}\_i\; \backslash hat\{n\}\_j,\; \backslash quad\; (2a)$$

$$e_{ij}^{\times }(\widehat{\mathrm{\Omega }})={\widehat{m}}_i{\widehat{n}}_j+{\widehat{n}}_i{\widehat{m}}_j,(2b)e\_\{ij\}^\backslash times(\backslash hat\{\backslash Omega\})\; =\; \backslash hat\{m\}\_i\; \backslash hat\{n\}\_j\; +\; \backslash hat\{n\}\_i\; \backslash hat\{m\}\_j,\; \backslash quad\; (2b)$$

where

$$\widehat{\mathrm{\Omega }}=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta ),(3a)\backslash hat\{\backslash Omega\}\; =\; (\backslash sin\; \backslash theta\; \backslash cos\; \backslash phi,\; \backslash sin\; \backslash theta\; \backslash sin\; \backslash phi,\; \backslash cos\; \backslash theta),\; \backslash quad\; (3a)$$

$$\widehat{m}=(\sin \varphi ,-\cos \varphi ,0),(3b)\backslash hat\{m\}\; =\; (\backslash sin\; \backslash phi,\; -\backslash cos\; \backslash phi,\; 0),\; \backslash quad\; (3b)$$

$$\widehat{n}=(\cos \theta \cos \varphi ,\cos \theta \sin \varphi ,-\sin \theta )\mathrm{.}(3c)\backslash hat\{n\}\; =\; (\backslash cos\; \backslash theta\; \backslash cos\; \backslash phi,\; \backslash cos\; \backslash theta\; \backslash sin\; \backslash phi,\; -\backslash sin\; \backslash theta).\; \backslash quad\; (3c)$$

Now consider the metric perturbation from a single gravitational wave traveling along the $z$ -axis so that $\hat{\Omega} = \hat{z}$ . The metric perturbation is given explicitly by

The physical metric due to the perturbation is given by

where $\eta_{ab} = \text{diag}{-1, 1, 1, 1}$ is the Minkowski metric, $a, b$ are spacetime indices, and

$$h_{+,\times }=h_{+,\times }(t-z)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{e}^{i2\pi f(t-z)}{h}_{+,\times }(f,\widehat{z})\mathrm{.}(6)h\_\{+,\; \backslash times\}\; =\; h\_\{+,\; \backslash times\}(t-z)\; =\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; e^\{i2\backslash pi\; f(t-z)\}\; h\_\{+,\; \backslash times\}(f,\; \backslash hat\{z\}).\; \backslash quad\; (6)$$

In this background, a pulsar emitting pulses at frequency $\nu_0$ and direction cosines $\alpha$ , $\beta$ , and $\gamma$ , with respect to the $x$ -, $y$ -, and $z$ -axes, respectively, will be observed to change its frequency in the solar system reference frame according to [3]

where $h_{+, \times}^e$ , $h_{+, \times}^p$ are the gravitational wave strains at the solar-system barycenter and the pulsar, respectively. This central result was obtained by Detweiler [3], who generalized a result of Estabrook and Wahlquist [2] toinclude both gravitational wave polarizations and for pulsars at arbitrary locations. They in turn based their calculation on an earlier one by Kaufmann [48]. A detailed derivation of this result is provided in Appendix A.

Looking at Eq. (7) (and as shown in Appendix B) we can write the redshift $z(t, \hat{\Omega})$ of signals from a pulsar in the direction of the unit vector $\hat{p}$ produced by a gravitational wave coming from the direction $\hat{\Omega}$ as,

$$z(t,\widehat{\mathrm{\Omega }})=\frac{1}{2}\frac{{\widehat{p}}^i{\widehat{p}}^j}{1+\widehat{\mathrm{\Omega }}\cdot \widehat{p}}\mathrm{\Delta }{h}_{ij},(8)z(t,\; \backslash hat\{\backslash Omega\})\; =\; \backslash frac\{1\}\{2\}\; \backslash frac\{\backslash hat\{p\}^i\; \backslash hat\{p\}^j\}\{1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\}\; \backslash Delta\; h\_\{ij\},\; \backslash quad\; (8)$$

where

$$\mathrm{\Delta }{h}_{ij}\equiv h_{ij}(t_p,\widehat{\mathrm{\Omega }})-h_{ij}(t_e,\widehat{\mathrm{\Omega }}),(9)\backslash Delta\; h\_\{ij\}\; \backslash equiv\; h\_\{ij\}(t\_p,\; \backslash hat\{\backslash Omega\})\; -\; h\_\{ij\}(t\_e,\; \backslash hat\{\backslash Omega\}),\; \backslash quad\; (9)$$

is the difference in the metric perturbation traveling along the direction $\hat{\Omega}$ at the pulsar and at the center of the solar system. The vectors $(t_e, \vec{x}_e)$ and $(t_p, \vec{x}_p)$ give the spacetime coordinates of the solar-system barycenter and the pulsar, respectively. The metric perturbation at each location takes the form,

$$h_{ij}(t,\widehat{\mathrm{\Omega }})=\sum _A{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{e}^{i2\pi f(t-\widehat{\mathrm{\Omega }}\cdot {\vec{x}}_0)}{h}_A(f,\widehat{\mathrm{\Omega }})e_{ij}^A(\widehat{\mathrm{\Omega }}),(10)h\_\{ij\}(t,\; \backslash hat\{\backslash Omega\})\; =\; \backslash sum\_A\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; e^\{i2\backslash pi\; f(t\; -\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash vec\{x\}\_0)\}\; h\_A(f,\; \backslash hat\{\backslash Omega\})\; e\_\{ij\}^A(\backslash hat\{\backslash Omega\}),\; \backslash quad\; (10)$$

for a fixed $\vec{x}_0$ .

We choose a particular coordinate system by placing the solar-system barycenter at the origin and the pulsar some distance $L$ away. With these conventions

$$t_p=t_e-L\equiv t-L,(11)t\_p\; =\; t\_e\; -\; L\; \backslash equiv\; t\; -\; L,\; \backslash quad\; (11)$$

$${\vec{x}}_e=0,(12)\backslash vec\{x\}\_e\; =\; 0,\; \backslash quad\; (12)$$

$${\vec{x}}_p=L\widehat{p}\mathrm{.}(13)\backslash vec\{x\}\_p\; =\; L\backslash hat\{p\}.\; \backslash quad\; (13)$$

If assume that the amplitude of the metric perturbation is the same at the solar-system barycenter and the pulsar then we can use Eq. (10) to write out $\Delta h_{ij}$ in our coordinate system as

Ultimately, we will be interested in the Fourier transform of this quantity which is simply

$$\mathrm{\Delta }{\tilde{h}}_{ij}(f,\widehat{\mathrm{\Omega }})=(e^{-i2\pi fL(1+\widehat{\mathrm{\Omega }}\cdot \widehat{p})}-1)\sum _A{h}_A(f,\widehat{\mathrm{\Omega }})e_{ij}^A(\widehat{\mathrm{\Omega }})\mathrm{.}(15)\backslash Delta\; \backslash tilde\{h\}\_\{ij\}(f,\; \backslash hat\{\backslash Omega\})\; =\; \backslash left(\; e^\{-i2\backslash pi\; fL(1+\backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\})\}\; -\; 1\; \backslash right)\; \backslash sum\_A\; h\_A(f,\; \backslash hat\{\backslash Omega\})\; e\_\{ij\}^A(\backslash hat\{\backslash Omega\}).\; \backslash quad\; (15)$$

We can then write the Fourier transform of Eq. (8) as

$$\tilde{z}(f,\widehat{\mathrm{\Omega }})=(e^{-i2\pi fL(1+\widehat{\mathrm{\Omega }}\cdot \widehat{p})}-1)\sum _A{h}_A(f,\widehat{\mathrm{\Omega }})F^A(\widehat{\mathrm{\Omega }}),(16)\backslash tilde\{z\}(f,\; \backslash hat\{\backslash Omega\})\; =\; \backslash left(\; e^\{-i2\backslash pi\; fL(1+\backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\})\}\; -\; 1\; \backslash right)\; \backslash sum\_A\; h\_A(f,\; \backslash hat\{\backslash Omega\})\; F^A(\backslash hat\{\backslash Omega\}),\; \backslash quad\; (16)$$

where we have defined

$$F^A(\widehat{\mathrm{\Omega }})\equiv e_{ij}^A(\widehat{\mathrm{\Omega }})\frac{1}{2}\frac{{\widehat{p}}^i{\widehat{p}}^j}{1+\widehat{\mathrm{\Omega }}\cdot \widehat{p}}\mathrm{.}(17)F^A(\backslash hat\{\backslash Omega\})\; \backslash equiv\; e\_\{ij\}^A(\backslash hat\{\backslash Omega\})\; \backslash frac\{1\}\{2\}\; \backslash frac\{\backslash hat\{p\}^i\; \backslash hat\{p\}^j\}\{1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\}.\; \backslash quad\; (17)$$

As shown in Appendix B the total redshift is given by summing over the contributions coming from gravitational waves in every direction:

$$\tilde{z}(f)={\int }_{S^2}d\widehat{\mathrm{\Omega }}\tilde{z}(f,\widehat{\mathrm{\Omega }}),(18)\backslash tilde\{z\}(f)\; =\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; \backslash tilde\{z\}(f,\; \backslash hat\{\backslash Omega\}),\; \backslash quad\; (18)$$

and similarly for $z(t)$ .

In fact, it is not the redshift, but a related quantity called the residual that gets reported in pulsar timing measurements. The residual, $R(t)$ , is defined as the integral of the redshift:

$$R(t)\equiv {\int }_0^{t}d{t}^{\mathrm{\prime }}z(t^{\mathrm{\prime }})\mathrm{.}(19)R(t)\; \backslash equiv\; \backslash int\_0^t\; dt\text{'}\; z(t\text{'}).\; \backslash quad\; (19)$$

This simple relationship gives us the freedom to develop the data analysis for either variable and we henceforth limit our attention to the redshift, but the results here can be phrased in terms of the residual with minimal effort.

In the literature, searches for gravitational waves using pulsar timing data are typically performed in the time domain. The (unknown) metric perturbation at the pulsar in, say, Eqs. (7) or (8) is neglected because one can treat it as another noise term which averages to zero when performing correlations between measurements of different pulsars. In the frequency domain this is unnecessary. Eq. (16) does not depend explicitly on the metric perturbation at the pulsar, rather the dependence is all in a distance and frequency dependent phase factor. It is then conceivable that if we could determine the distance to a pulsar $L$ with sufficient accuracy we could use the metric perturbation at the pulsar to improve the sensitivity of our searches. Unfortunately, such measurements of pulsar distances are unavailable. We will show in Section III A, however, that for the case of a stochastic background search in pulsar timing data the phase factor can be neglected without any significant loss in sensitivity. It is unclear whether this is true for other types of gravitational wave searches.

From the ground-based interferometer perspective Eqs. (7) or (8) are somewhat counter-intuitive. This difficulty arises from the factor of $1 + \hat{\Omega} \cdot \hat{p}$ in the denominator; in Appendices A and B we show explicitly how this factor enters the expression. When $\hat{\Omega} \cdot \hat{p} = \pm 1$ , i.e. the gravitational wave and the pulsar directions are parallel or anti-parallel, Eqs. (16) and (17) lead to no redshifting of the pulsar signal for completely different reasons. When they are parallel the reason is the transverse nature of gravitational waves, and when they are anti-parallel it is because the pulsar signals “surf” the gravitational waves. Our surprise is a result of our long-wavelength limit intuition.Eq. (16) has an obvious long-wavelength limit. We can use this limit to compare the form of our results with those of ground-based interferometers such as LIGO. When $2\pi fL \ll 1$ we can Taylor expand the exponential and to first order Eq. (16) becomes

$$\tilde{z}(f,\widehat{\mathrm{\Omega }})\approx -i\pi fL{\widehat{p}}^i{\widehat{p}}^j\sum _A{h}_A(f,\widehat{\mathrm{\Omega }})e_{ij}^A(\widehat{\mathrm{\Omega }})\mathrm{.}(20)\backslash tilde\{z\}(f,\; \backslash hat\{\backslash Omega\})\; \backslash approx\; -i\backslash pi\; fL\; \backslash hat\{p\}^i\; \backslash hat\{p\}^j\; \backslash sum\_A\; h\_A(f,\; \backslash hat\{\backslash Omega\})\; e\_\{ij\}^A(\backslash hat\{\backslash Omega\}).\; \backslash quad\; (20)$$

Typical values of $f$ are in the range $1/10 \text{ yr}^{-1}$ to $10 \text{ yr}^{-1}$ . Typical values of the Earth pulsar distance $L$ are in the range $100 \text{ ly}$ to $10^4 \text{ ly}$ . This means $fL$ is in the range $10$ to $10^5$ and pulsar timing experiments are never in the long wavelength limit. However, the Taylor expansion can also be done for large $fL$ when the angle between $\hat{\Omega}$ and $\hat{p}$ is sufficiently close to $\pi$ . In this case the pulsar signals can “surf” the gravitational waves and not undergo redshifting. Writing that angle as $\pi - \epsilon$ with $\epsilon \ll 1$ , then the Taylor expansion is also valid when $\epsilon \ll (\pi fL)^{-1/2}$ .

Taking the inverse Fourier transform of Eq. (20) yields

$$\tilde{z}(t)\approx -\frac{L}{2}{\widehat{p}}^i{\widehat{p}}^j{\dot{h}}_{ij}(t,{\vec{x}}_e),(21)\backslash tilde\{z\}(t)\; \backslash approx\; -\backslash frac\{L\}\{2\}\; \backslash hat\{p\}^i\; \backslash hat\{p\}^j\; \backslash dot\{h\}\_\{ij\}(t,\; \backslash vec\{x\}\_e),\; \backslash quad\; (21)$$

which is the projection of the time derivative of the metric perturbation at the solar-system barycenter onto the unit vector that points to the pulsar. Note that unlike Eq. (8), this equation no longer depends on the direction of the gravitational wave and can be expressed in terms of the full metric perturbation (derivative). For the case of ground based interferometers the signal, the so-called strain, is proportional to the difference in length of the two arms because the signal at the dark port of the interferometer depends on that difference. If the arms point in the directions of the unit vectors $\hat{X}$ and $\hat{Y}$ the strain is given by

$$h(t)\equiv h_{ij}(t,\vec{x})\frac{1}{2}({\widehat{X}}^i{\widehat{X}}^j-{\widehat{Y}}^i{\widehat{Y}}^j),(22)h(t)\; \backslash equiv\; h\_\{ij\}(t,\; \backslash vec\{x\})\; \backslash frac\{1\}\{2\}\; (\backslash hat\{X\}^i\; \backslash hat\{X\}^j\; -\; \backslash hat\{Y\}^i\; \backslash hat\{Y\}^j),\; \backslash quad\; (22)$$

which is the metric perturbation $h_{ij}(t, \vec{x})$ projected onto the difference of the arms.

III. DETECTION STATISTIC

With an understanding of the signal in hand we now turn our attention to developing an optimal detection strategy. In this section we will first derive the optimal cross-correlation statistic for a single pulsar pair using arguments based on maximizing signal to noise ratio. We will then determine the best way to combine measurements from multiple pulsar pairs to obtain the most constraining upper limit. This section will conclude with a more sophisticated derivation of the optimal detection statistic based on the likelihood ratio.

A. The optimal filter

In this section we will derive the optimal filter for detecting a stochastic background of gravitational waves

from the cross-correlation of redshift measurements of two different pulsars. This problem was addressed in detail by Allen and Romano [40], for the case of interferometric gravitational wave detectors and our analysis follows theirs closely.

Consider the signals from two pulsars

$$s_1(t)=z_1(t)+n_1(t),(23)s\_1(t)\; =\; z\_1(t)\; +\; n\_1(t),\; \backslash quad\; (23)$$

$$s_2(t)=z_2(t)+n_2(t),(24)s\_2(t)\; =\; z\_2(t)\; +\; n\_2(t),\; \backslash quad\; (24)$$

where $z_i(t)$ is the redshift and $n_i(t)$ is the noise intrinsic in the measurement. Throughout this work we will assume that each $n_i(t)$ is stationary and Gaussian, and is greater in magnitude than the redshift. Additionally we assume that

for all $i$ and $j$ , where the angle brackets denote an expectation value.

A stochastic background will show up in the data as correlated noise between measurements with different detectors. Our goal is to find a filter, $Q(t - t')$ , that optimizes the cross-correlation statistic

$$S\equiv {\int }_{-T\mathrm{/}2}^{T\mathrm{/}2}dt{\int }_{-T\mathrm{/}2}^{T\mathrm{/}2}d{t}^{\mathrm{\prime }}{s}_1(t)s_2(t^{\mathrm{\prime }})Q(t-t^{\mathrm{\prime }}),(26)S\; \backslash equiv\; \backslash int\_\{-T/2\}^\{T/2\}\; dt\; \backslash int\_\{-T/2\}^\{T/2\}\; dt\text{'}\; s\_1(t)\; s\_2(t\text{'})\; Q(t\; -\; t\text{'}),\; \backslash quad\; (26)$$

where $T$ is the observation time. We will define the optimal filter to be the $Q(t - t')$ that maximizes the signal to noise ratio (SNR)

$$\text{SNR}\equiv \frac{\mu }{\sigma },(27)\backslash text\{SNR\}\; \backslash equiv\; \backslash frac\{\backslash mu\}\{\backslash sigma\},\; \backslash quad\; (27)$$

where $\mu$ and $\sigma$ are the mean and square root of the variance, respectively, associated with the cross-correlation signal defined in Eq. (26).

We start by assuming that the observation time is much greater than the separation of the two detectors and extend the limits of the integral over $dt'$ to $\pm\infty$ . Technically our assumption is not correct because pulsars are typically separated by distances far greater than the observation time. Later we will see that neglecting the phase terms that correspond to the metric perturbation at the pulsar location is an excellent approximation. In effect this makes our detectors co-located though not co-aligned, and our assumption about the observation time is appropriate. We work in the frequency domain so that Eq. (26) becomes

$$S={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}d{f}^{\mathrm{\prime }}{\delta }_T(f-f^{\mathrm{\prime }}){\tilde{s}}_1^{\ast }(f){\tilde{s}}_2(f^{\mathrm{\prime }})\tilde{Q}(f^{\mathrm{\prime }}),(28)S\; =\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\text{'}\; \backslash delta\_T(f\; -\; f\text{'})\; \backslash tilde\{s\}\_1^*(f)\; \backslash tilde\{s\}\_2(f\text{'})\; \backslash tilde\{Q\}(f\text{'}),\; \backslash quad\; (28)$$

where $\delta_T(f - f')$ is the finite-time approximation to the delta function

$${\delta }_T(f)=\frac{\sin (\pi fT)}{\pi f}\mathrm{.}(29)\backslash delta\_T(f)\; =\; \backslash frac\{\backslash sin(\backslash pi\; f\; T)\}\{\backslash pi\; f\}.\; \backslash quad\; (29)$$The mean of the cross-correlation is

$$\mu \equiv ⟨S⟩={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}d{f}^{\mathrm{\prime }}{\delta }_T(f-f^{\mathrm{\prime }})⟨{\tilde{z}}_1^{\ast }(f){\tilde{z}}_2(f^{\mathrm{\prime }})⟩\tilde{Q}(f^{\mathrm{\prime }})\mathrm{.}(30)\backslash mu\; \backslash equiv\; \backslash langle\; S\; \backslash rangle\; =\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\text{'}\; \backslash delta\_T(f\; -\; f\text{'})\; \backslash langle\; \backslash tilde\{z\}\_1^*(f)\; \backslash tilde\{z\}\_2(f\text{'})\; \backslash rangle\; \backslash tilde\{Q\}(f\text{'}).\; \backslash quad\; (30)$$

Because of Eqs. (16), (18), and (25), taking the expectation value above requires us to evaluate $\langle h_A^*(f, \hat{\Omega}) h_{A'}(f', \hat{\Omega}') \rangle$ . The assumptions that the stochastic background is stationary, unpolarized and isotropic lead us to take

$$⟨h_A^{\ast }(f,\widehat{\mathrm{\Omega }})h_{A^{\mathrm{\prime }}}(f^{\mathrm{\prime }},{\widehat{\mathrm{\Omega }}}^{\mathrm{\prime }})⟩={\delta }^2(\widehat{\mathrm{\Omega }},{\widehat{\mathrm{\Omega }}}^{\mathrm{\prime }}){\delta }_{A{A}^{\mathrm{\prime }}}\delta (f-f^{\mathrm{\prime }})H(f),(31)\backslash langle\; h\_A^*(f,\; \backslash hat\{\backslash Omega\})\; h\_\{A\text{'}\}(f\text{'},\; \backslash hat\{\backslash Omega\}\text{'})\; \backslash rangle\; =\; \backslash delta^2(\backslash hat\{\backslash Omega\},\; \backslash hat\{\backslash Omega\}\text{'})\; \backslash delta\_\{AA\text{'}\}\; \backslash delta(f\; -\; f\text{'})\; H(f),\; \backslash quad\; (31)$$

where $H(f) = H(-f)$ is the gravitational wave spectrum. $H(f)$ is related to $\Omega_{\text{gw}}(f)$ through

$${\mathrm{\Omega }}_{\text{gw}}(f)\equiv \frac{1}{{\rho }_{\text{crit}}}\frac{d{\rho }_{\text{gw}}}{d\ln f},(32)\backslash Omega\_\{\backslash text\{gw\}\}(f)\; \backslash equiv\; \backslash frac\{1\}\{\backslash rho\_\{\backslash text\{crit\}\}\}\; \backslash frac\{d\backslash rho\_\{\backslash text\{gw\}\}\}\{d\; \backslash ln\; f\},\; \backslash quad\; (32)$$

where $\rho_{\text{crit}} = 8\pi/3H_0^2$ and

$${\rho }_{\text{gw}}=\frac{1}{32\pi }⟨{\dot{h}}_{ab}(t,\vec{x}){\dot{h}}^{ab}(t,\vec{x})⟩,(33)\backslash rho\_\{\backslash text\{gw\}\}\; =\; \backslash frac\{1\}\{32\backslash pi\}\; \backslash langle\; \backslash dot\{h\}\_\{ab\}(t,\; \backslash vec\{x\})\; \backslash dot\{h\}^\{ab\}(t,\; \backslash vec\{x\})\; \backslash rangle,\; \backslash quad\; (33)$$

is the energy density in gravitational waves. It follows from the plane wave expansion Eq. (1) along with Eqs. (31) and (32) in Eq. (33) that

$$H(f)=\frac{3H_0^2}{32{\pi }^3}\mathrm{\mid }f{\mathrm{\mid }}^{-3}{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid }),(34)H(f)\; =\; \backslash frac\{3H\_0^2\}\{32\backslash pi^3\}\; |f|^\{-3\}\; \backslash Omega\_\{\backslash text\{gw\}\}(|f|),\; \backslash quad\; (34)$$

and therefore

$$⟨h_A^{\ast }(f,\widehat{\mathrm{\Omega }})h_{A^{\mathrm{\prime }}}(f^{\mathrm{\prime }},{\widehat{\mathrm{\Omega }}}^{\mathrm{\prime }})⟩=\frac{3H_0^2}{32{\pi }^3}{\delta }^2(\widehat{\mathrm{\Omega }},{\widehat{\mathrm{\Omega }}}^{\mathrm{\prime }}){\delta }_{A{A}^{\mathrm{\prime }}}\delta (f-f^{\mathrm{\prime }})\times \mathrm{\mid }f{\mathrm{\mid }}^{-3}{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid }),(35)\backslash langle\; h\_A^*(f,\; \backslash hat\{\backslash Omega\})\; h\_\{A\text{'}\}(f\text{'},\; \backslash hat\{\backslash Omega\}\text{'})\; \backslash rangle\; =\; \backslash frac\{3H\_0^2\}\{32\backslash pi^3\}\; \backslash delta^2(\backslash hat\{\backslash Omega\},\; \backslash hat\{\backslash Omega\}\text{'})\; \backslash delta\_\{AA\text{'}\}\; \backslash delta(f\; -\; f\text{'})\; \backslash times\; |f|^\{-3\}\; \backslash Omega\_\{\backslash text\{gw\}\}(|f|),\; \backslash quad\; (35)$$

which is sometimes written in terms of the characteristic strain

$$h_c^2(f)=\frac{3H_0^2}{2{\pi }^2}\frac{1}{f^2}{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid })\mathrm{.}(36)h\_c^2(f)\; =\; \backslash frac\{3H\_0^2\}\{2\backslash pi^2\}\; \backslash frac\{1\}\{f^2\}\; \backslash Omega\_\{\backslash text\{gw\}\}(|f|).\; \backslash quad\; (36)$$

The expectation value we set out to evaluate is then

$$⟨{\tilde{z}}_1^{\ast }(f){\tilde{z}}_2(f^{\mathrm{\prime }})⟩=\frac{3H_0^2}{32{\pi }^3}\frac{1}{\beta }\delta (f-f^{\mathrm{\prime }})\mathrm{\mid }f{\mathrm{\mid }}^{-3}{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid })\mathrm{\Gamma }(\mathrm{\mid }f\mathrm{\mid }),(37)\backslash langle\; \backslash tilde\{z\}\_1^*(f)\; \backslash tilde\{z\}\_2(f\text{'})\; \backslash rangle\; =\; \backslash frac\{3H\_0^2\}\{32\backslash pi^3\}\; \backslash frac\{1\}\{\backslash beta\}\; \backslash delta(f\; -\; f\text{'})\; |f|^\{-3\}\; \backslash Omega\_\{\backslash text\{gw\}\}(|f|)\; \backslash Gamma(|f|),\; \backslash quad\; (37)$$

where we defined

$$\mathrm{\Gamma }(\mathrm{\mid }f\mathrm{\mid })=\beta \sum _A{\int }_{S^2}d\widehat{\mathrm{\Omega }}(e^{i2\pi f{L}_1(1+\widehat{\mathrm{\Omega }}\cdot {\widehat{p}}_1)}-1)\times (e^{-i2\pi f{L}_2(1+\widehat{\mathrm{\Omega }}\cdot {\widehat{p}}_2)}-1)F_1^A(\widehat{\mathrm{\Omega }})F_2^A(\widehat{\mathrm{\Omega }}),(38)\backslash Gamma(|f|)\; =\; \backslash beta\; \backslash sum\_A\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; \backslash left(\; e^\{i2\backslash pi\; f\; L\_1\; (1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\_1)\}\; -\; 1\; \backslash right)\; \backslash times\; \backslash left(\; e^\{-i2\backslash pi\; f\; L\_2\; (1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\_2)\}\; -\; 1\; \backslash right)\; F\_1^A(\backslash hat\{\backslash Omega\})\; F\_2^A(\backslash hat\{\backslash Omega\}),\; \backslash quad\; (38)$$

the pulsar timing analogue of the overlap reduction function [40], which has a normalization factor $\beta$ . The normalization is chosen so that $\Gamma(|f|) = 1$ for coincident, co-aligned detectors. As we show below, pulsar timing experiments are in a regime where the exponential factors

FIG. 1: Plot of the full overlap reduction, Eq. (38), along with the approximation Eq. (39) for two pulsars a distance $L$ from the solar-system barycenter. The overlap reduction function is a function of $fL$ . The top two (blue) curves show Eq. (38) with $\beta = 3/4\pi$ (solid line) and Eq. (39) (dashed line) for two pulsars at an angle $\xi = \pi/8$ as a function of $fL$ . The middle (red) and bottom (green) curves show the same quantities for two pulsars at $\xi \approx 0.86$ and $\xi = \pi/2$ respectively. The smallest value of the frequency $f_{\min} \sim 0.1 \text{ yr}^{-1}$ and the closest pulsars used in timing experiments are at a distance of $L_{\min} \sim 100 \text{ ly}$ so that $fL \gtrsim 10$ . This range of $fL$ puts pulsar timing experiments in the regime where Eq. (39) is an excellent approximation to Eq. (38) and we are justified in throwing out the pulsar term while remaining close to optimal.

in Eq. (38) can be neglected. In this situation, which we will assume henceforth, the normalization factor is easy to determine and we have that

where $\xi = \cos^{-1}(\hat{p}_1 \cdot \hat{p}_2)$ is the angle between the two pulsars. This quantity is proportional to the Hellings and Downs curve [49]. A detailed derivation of this result is provided for completeness in Appendix C 1.

The rationale given in the literature for throwing out the pulsar term in Eq. (7), or equivalently Eqs. (8) and (9), is that the unknown metric perturbation at the pulsars can be thought of as a kind of noise term which averages to zero when performing a correlation between different pulsars. The equivalent procedure in the frequency domain is to neglect the phase factors in Eq. (16), or in terms of our optimal filter, approximating Eq. (38) with Eq. (39). The regime where the approximate Eq. (39) isvalid, is helpful in quantifying the accuracy of the rationale. Figure 1 shows the overlap reduction function for two pulsars a distance $L$ from the solar-system barycenter. Since the distance to both pulsars is the same, the overlap reduction function can be written as just a function of $fL$ . The top two (blue) curves show Eq. (38) with $\beta = 3/4\pi$ (solid line) and Eq. (39) (dashed line) for two pulsars at an angle $\xi = \pi/8$ as a function of $fL$ . The middle (red) and bottom (green) curves show the same quantities for two pulsars at $\xi \approx 0.86$ and $\xi = \pi/2$ respectively. As discussed in the last section the smallest value of the frequency $f_{\min} \sim 0.1 \text{ yr}^{-1}$ and the closest pulsars used in timing experiments are at a distance of $L_{\min} \sim 100 \text{ ly}$ so that $fL \gtrsim 10$ . As shown in Fig. 1 this range of $fL$ puts pulsar timing experiments in the regime where Eq. (39) is an excellent approximation to Eq. (38), and we can neglect the pulsar term while remaining close to optimal.

Returning to Eq. (30), we now have

With the assumption that the noise is much greater than the signal, the variance, $\sigma^2$ , depends only on the statistical properties of the noise in each detector. We have

where

$$⟨{\tilde{n}}_i^{\ast }(f){\tilde{n}}_i(f^{\mathrm{\prime }})⟩=\frac{1}{2}\delta (f-f^{\mathrm{\prime }})P_i(\mathrm{\mid }f\mathrm{\mid }),(42)\backslash langle\; \backslash tilde\{n\}\_i^*(f)\; \backslash tilde\{n\}\_i(f\text{'})\; \backslash rangle\; =\; \backslash frac\{1\}\{2\}\; \backslash delta(f\; -\; f\text{'})\; P\_i(|f|),\; \backslash quad\; (42)$$

is the (one-sided) noise power spectrum.

With $\mu$ and $\sigma^2$ in hand, we next define a positive-definite inner product using the noise power spectra of the two detectors

$$(A,B)\equiv {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{A}^{\ast }(f)B(f)P_1(\mathrm{\mid }f\mathrm{\mid })P_2(\mathrm{\mid }f\mathrm{\mid })\mathrm{.}(43)(A,\; B)\; \backslash equiv\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; A^*(f)\; B(f)\; P\_1(|f|)\; P\_2(|f|).\; \backslash quad\; (43)$$

With this definition it is easy to see that

$$\mu \approx \frac{H_0^2}{8{\pi }^2}T(\tilde{Q},\frac{{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid }){\mathrm{\Gamma }}_0}{\mathrm{\mid }f{\mathrm{\mid }}^3{P}_1(\mathrm{\mid }f\mathrm{\mid })P_2(\mathrm{\mid }f\mathrm{\mid })}),(44)\backslash mu\; \backslash approx\; \backslash frac\{H\_0^2\}\{8\backslash pi^2\}\; T\; \backslash left(\; \backslash tilde\{Q\},\; \backslash frac\{\backslash Omega\_\{\backslash text\{gw\}\}(|f|)\; \backslash Gamma\_0\}\{|f|^3\; P\_1(|f|)\; P\_2(|f|)\}\; \backslash right),\; \backslash quad\; (44)$$

$${\sigma }^2\approx \frac{T}{4}(\tilde{Q},\tilde{Q}),(45)\backslash sigma^2\; \backslash approx\; \backslash frac\{T\}\{4\}\; (\backslash tilde\{Q\},\; \backslash tilde\{Q\}),\; \backslash quad\; (45)$$

from which it follows from the definition of SNR in Eq. (27) and Schwartz's inequality that the optimal filter is given by

$$\tilde{Q}(f)=\chi \frac{{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid }){\mathrm{\Gamma }}_0}{\mathrm{\mid }f{\mathrm{\mid }}^3{P}_1(\mathrm{\mid }f\mathrm{\mid })P_2(\mathrm{\mid }f\mathrm{\mid })},(46)\backslash tilde\{Q\}(f)\; =\; \backslash chi\; \backslash frac\{\backslash Omega\_\{\backslash text\{gw\}\}(|f|)\; \backslash Gamma\_0\}\{|f|^3\; P\_1(|f|)\; P\_2(|f|)\},\; \backslash quad\; (46)$$

for some normalization constant, $\chi$ . Our primary interest will be in stochastic backgrounds with power law spectra, $\Omega_{\text{gw}}(f) = \Omega_\alpha f^\alpha$ (for constant $\Omega_\alpha$ ). In that case the normalization constant for the optimal filter, $\tilde{Q}_\alpha(f)$ , is chosen so that

$$\mu ={\mathrm{\Omega }}_{\alpha }{T}_0,(47)\backslash mu\; =\; \backslash Omega\_\backslash alpha\; T\_0,\; \backslash quad\; (47)$$

where $T_0$ is some arbitrary constant with dimensions of time. From Eq. (44) it follows that

$$\chi ={\mathrm{\Omega }}_{\alpha }\frac{T_0}{T}\frac{8{\pi }^2}{H_0^2}{\left[{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df\frac{{\mathrm{\Omega }}_{\text{gw}}^2(\mathrm{\mid }f\mathrm{\mid }){\mathrm{\Gamma }}_0^2}{f^6{P}_1(\mathrm{\mid }f\mathrm{\mid })P_2(\mathrm{\mid }f\mathrm{\mid })}\right]}^{-1}\mathrm{.}(48)\backslash chi\; =\; \backslash Omega\_\backslash alpha\; \backslash frac\{T\_0\}\{T\}\; \backslash frac\{8\backslash pi^2\}\{H\_0^2\}\; \backslash left[\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash frac\{\backslash Omega\_\{\backslash text\{gw\}\}^2(|f|)\; \backslash Gamma\_0^2\}\{f^6\; P\_1(|f|)\; P\_2(|f|)\}\; \backslash right]^\{-1\}.\; \backslash quad\; (48)$$

Finally, we can compute

$$\text{SNR}\approx \frac{H_0^2}{4{\pi }^2}{T}^{1\mathrm{/}2}{\left[{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df\frac{{\mathrm{\Omega }}_{\text{gw}}^2(\mathrm{\mid }f\mathrm{\mid }){\mathrm{\Gamma }}_0^2}{f^6{P}_1(\mathrm{\mid }f\mathrm{\mid })P_2(\mathrm{\mid }f\mathrm{\mid })}\right]}^{1\mathrm{/}2}\mathrm{.}(49)\backslash text\{SNR\}\; \backslash approx\; \backslash frac\{H\_0^2\}\{4\backslash pi^2\}\; T^\{1/2\}\; \backslash left[\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash frac\{\backslash Omega\_\{\backslash text\{gw\}\}^2(|f|)\; \backslash Gamma\_0^2\}\{f^6\; P\_1(|f|)\; P\_2(|f|)\}\; \backslash right]^\{1/2\}.\; \backslash quad\; (49)$$

The differences between these results and those for interferometers can all be traced to the differing overlap reduction function $\Gamma(f) \approx \Gamma_0$ . The normalization of $\Gamma_0$ means that the maximal SNR (for co-incident, co-aligned detectors) is only 5/6 of that obtainable from interferometers, assuming the noise power spectra are the same in each case.

To construct the optimal filter, Eq. (46), the noise power spectra for the two pulsars $P_1(|f|)$ and $P_2(|f|)$ must be determined. These can either be modeled, or measured with the methods described in Section VI. Once constructed the optimal filter can be applied in the frequency domain. Section VI gives a prescription for taking Fourier transforms of unevenly sampled data. The optimal filter can also be inverse Fourier transformed and the correlation performed in the time domain. It is unclear which of these two methods is more robust and the authors of [46] will explore the time-domain approach in detail.

B. The pulsar timing array

The question we would like to address in this section is: Given redshift measurements from $N$ different pulsars (which each have a different noise profile), what is the best way to combine those measurements to produce the most constraining upper limit? One can consider the cross-correlations between any even number of detectors, but it has been shown [40, 50] that the optimal choice is the combination of pairwise cross-correlations. As it turns out, the solution to this problem also solves the problem of non-stationarity in the noise power spectra over periods longer than the typical observation time, $T$ .

First let

$${}^{(ij)}{S}_1{,}^{(ij)}{S}_2,\dots {,}^{(ij)}{S}_{n_{ij}},(50)^\{(ij)\}\; S\_1,\; ^\{(ij)\}\; S\_2,\; \backslash dots,\; ^\{(ij)\}\; S\_\{n\_\{ij\}\},\; \backslash quad\; (50)$$

be $n_{ij}$ measurements of the cross-correlation between the $i$ -th and $j$ -th pulsar. We will assume that eachmeasurement is taken with an optimal filter normalized so that while searching for a background of the form $\Omega_{\text{gw}}(f) = \Omega_\alpha f^\alpha$ ,

$${⟨}^{(ij)}{S}_k⟩={\mathrm{\Omega }}_{\alpha }{T}_0\equiv \mu ,(51)\backslash langle^\{(ij)\}S\_k\backslash rangle\; =\; \backslash Omega\_\backslash alpha\; T\_0\; \backslash equiv\; \backslash mu,\; \backslash quad\; (51)$$

where $T_0$ is an arbitrary constant introduced for dimensional reasons. Each measurement therefore has the form

$${}^{(ij)}{S}_k={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}d{f}^{\mathrm{\prime }}{\delta }_T(f-f^{\mathrm{\prime }}){\tilde{s}}_{i,k}^{\ast }(f){\tilde{s}}_{j,k}(f^{\mathrm{\prime }}{)}^{(ij)}{Q}_k(f)(52)^\{(ij)\}S\_k\; =\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\text{'}\; \backslash delta\_T(f-f\text{'})\; \backslash tilde\{s\}\_\{i,k\}^*(f)\; \backslash tilde\{s\}\_\{j,k\}(f\text{'})^\{(ij)\}\; Q\_k(f)\; \backslash quad\; (52)$$

with

$${}^{(ij)}{Q}_k(f){=}^{(ij)}{\chi }_k\frac{{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid }{)}^{(ij)}{\mathrm{\Gamma }}_0}{\mathrm{\mid }f{\mathrm{\mid }}^3{P}_{i,k}(\mathrm{\mid }f\mathrm{\mid })P_{j,k}(\mathrm{\mid }f\mathrm{\mid })}\mathrm{.}(53)^\{(ij)\}Q\_k(f)\; =^\{(ij)\}\backslash chi\_k\; \backslash frac\{\backslash Omega\_\{\backslash text\{gw\}\}(|f|)^\{(ij)\}\backslash Gamma\_0\}\{|f|^3\; P\_\{i,k\}(|f|)\; P\_\{j,k\}(|f|)\}.\; \backslash quad\; (53)$$

where $^{(ij)}\Gamma_0$ is the overlap reduction function of the $(ij)$ pulsar pair, $s_{i,k}$ is the $k$ -th measurement of the signal from the $i$ -th pulsar, and $P_{i,k}$ is the associated noise power spectrum. Additionally,

$${}^{(ij)}{\chi }_k={\mathrm{\Omega }}_{\alpha }\frac{T_0}{{}^{(ij)}{T}_k}\frac{8{\pi }^2}{H_0^2}{\left[{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df\frac{{\mathrm{\Omega }}_{\text{gw}}^2(\mathrm{\mid }f\mathrm{\mid }{)}^{(ij)}{\mathrm{\Gamma }}_0^2}{f^6{P}_{i,k}(\mathrm{\mid }f\mathrm{\mid })P_{j,k}(\mathrm{\mid }f\mathrm{\mid })}\right]}^{-1},(54)^\{(ij)\}\backslash chi\_k\; =\; \backslash Omega\_\backslash alpha\; \backslash frac\{T\_0\}\{^\{(ij)\}T\_k\}\; \backslash frac\{8\backslash pi^2\}\{H\_0^2\}\; \backslash left[\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash frac\{\backslash Omega\_\{\backslash text\{gw\}\}^2(|f|)^\{(ij)\}\backslash Gamma\_0^2\}\{f^6\; P\_\{i,k\}(|f|)\; P\_\{j,k\}(|f|)\}\; \backslash right]^\{-1\},\; \backslash quad\; (54)$$

where $^{(ij)}T_k$ is the observation time of the $k$ -th measurement of the $(ij)$ pulsar pair. Our task is to combine the $^{(ij)}S_k$ in a way that optimizes SNR. The first step is to form the sample mean for each set of measurements

^{(ij)}\hat{\mu} \equiv \frac{1}{n_{ij}} \sum_{k=1}^{n_{ij}} ^{(ij)}S_k, \quad (55)

which is both an unbiased estimator and random variable. It therefore has a mean

$${\mu }_{ij}\equiv {⟨}^{(ij)}\widehat{\mu }⟩=\mu ,(56)\backslash mu\_\{ij\}\; \backslash equiv\; \backslash langle^\{(ij)\}\backslash hat\{\backslash mu\}\backslash rangle\; =\; \backslash mu,\; \backslash quad\; (56)$$

and a variance

$${\sigma }_{ij}^2\equiv {⟨}^{(ij)}{\widehat{\mu }}^2⟩-{⟨}^{(ij)}\widehat{\mu }{⟩}^2=\frac{{}^{(ij)}{\sigma }^2}{n_{ij}},(57)\backslash sigma\_\{ij\}^2\; \backslash equiv\; \backslash langle^\{(ij)\}\backslash hat\{\backslash mu\}^2\backslash rangle\; -\; \backslash langle^\{(ij)\}\backslash hat\{\backslash mu\}\backslash rangle^2\; =\; \backslash frac\{^\{(ij)\}\backslash sigma^2\}\{n\_\{ij\}\},\; \backslash quad\; (57)$$

where

$${}^{(ij)}{\sigma }^2=\sum _{k=1}^{n_{ij}}\frac{{}^{(ij)}{T}_k}{4}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df{}^{(ij)}{\chi }_k^2\frac{{\mathrm{\Omega }}_{\text{gw}}^2(\mathrm{\mid }f\mathrm{\mid }{)}^{(ij)}{\mathrm{\Gamma }}_0^2}{f^6{P}_{i,k}(\mathrm{\mid }f\mathrm{\mid })P_{j,k}(\mathrm{\mid }f\mathrm{\mid })}\mathrm{.}(58)^\{(ij)\}\backslash sigma^2\; =\; \backslash sum\_\{k=1\}^\{n\_\{ij\}\}\; \backslash frac\{^\{(ij)\}T\_k\}\{4\}\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash \; ^\{(ij)\}\backslash chi\_k^2\; \backslash frac\{\backslash Omega\_\{\backslash text\{gw\}\}^2(|f|)^\{(ij)\}\backslash Gamma\_0^2\}\{f^6\; P\_\{i,k\}(|f|)\; P\_\{j,k\}(|f|)\}.\; \backslash quad\; (58)$$

The next step is combine the sample mean for each set of measurements into a single estimator we can use to determine an upper bound on $\Omega_\alpha$ and hence $\Omega_{\text{gw}} = \Omega_\alpha f^\alpha$ . We do so by introducing an unbiased estimator consisting of a weighted average of the sample means

for some constants, $\lambda_{ij}$ , which has mean

$${\mu }_{\widehat{\mu }}\equiv ⟨\widehat{\mu }⟩=\mu ,(60)\backslash mu\_\{\backslash hat\{\backslash mu\}\}\; \backslash equiv\; \backslash langle\backslash hat\{\backslash mu\}\backslash rangle\; =\; \backslash mu,\; \backslash quad\; (60)$$

and variance

The object is to now determine the $\lambda_{ij}$ that maximizes the SNR of $\hat{\mu}$ . The (squared) SNR of $\hat{\mu}$ is

To find the $\lambda_{ij}$ that maximize the SNR, we exploit the same trick that led us to the optimal filter. Namely, we introduce an inner product

which allows us to write

$${\text{SNR}}_{\widehat{\mu }}^2\equiv {\mu }^2\frac{(\lambda ,{\sigma }^{-2})}{(\lambda ,\lambda )},(64)\backslash text\{SNR\}\_\{\backslash hat\{\backslash mu\}\}^2\; \backslash equiv\; \backslash mu^2\; \backslash frac\{(\backslash lambda,\; \backslash sigma^\{-2\})\}\{(\backslash lambda,\; \backslash lambda)\},\; \backslash quad\; (64)$$

from which it follows that choosing

$${\lambda }_{ij}\propto {\sigma }_{ij}^{-2},(65)\backslash lambda\_\{ij\}\; \backslash propto\; \backslash sigma\_\{ij\}^\{-2\},\; \backslash quad\; (65)$$

maximizes the SNR. The optimal statistic, choosing $\lambda_{ij} = \sigma_{ij}^{-2}$ , is then given by

\begin{aligned} S_{\text{opt}} &= \frac{\sum_{i=1}^l \sum_{j<i}^l \sigma_{ij}^{-2} \hat{\mu}}{\sum_{i=1}^l \sum_{j<i}^l \sigma_{ij}^{-2}} \\ &= \frac{\sum_{i=1}^l \sum_{j<i}^l \sigma_{ij}^{-2} \sum_{k=1}^{n_{ij}} ^{(ij)}S_k}{\sum_{i=1}^l \sum_{j<i}^l n_{ij} \sigma_{ij}^{-2}}. \end{aligned} \quad (66)

Because the estimator defined in Eq. (59) is unbiased and defined so that $\mu = \langle S_{\text{opt}} \rangle = \Omega_\alpha T_0$ , the estimate of $\hat{\Omega}_\alpha$ is found using

$${\widehat{\mathrm{\Omega }}}_{\alpha }=\frac{{\widehat{S}}_{\text{opt}}}{T_0},(67)\backslash hat\{\backslash Omega\}\_\backslash alpha\; =\; \backslash frac\{\backslash hat\{S\}\_\{\backslash text\{opt\}\}\}\{T\_0\},\; \backslash quad\; (67)$$where $\hat{S}{\text{opt}}$ is the measured value of the optimal statistic. The expected variance of $\hat{S}{\text{opt}}$ follows from Eq. (61),

Aside from maximizing the SNR, the linear combination of sample means that forms the optimal statistic in Eq. (66) serves two important and related purposes. First of all, as mentioned at the beginning of this section, weighing each $^{(ij)}\hat{\mu}$ by the inverse of the squared variance means that less noisy measurements (those with smaller variances) contribute more to the sum, which helps minimize the effect of long term non-stationarity. This is augmented by the normalization convention we chose in Eq. (51) for the mean of each measurement. Using Eq. (45) and Eq. (46) with the $\lambda$ that follows from Eq. (51), we see that $^{(ij)}\sigma^{-2} \propto ^{(ij)}T$ , and so longer observation times are also favored in the sum.

C. Computational procedure

In this subsection we describe how the quantities necessary for a stochastic background search are computed. The goal is to produce a measurement of the optimal statistic, $\hat{S}_{\text{opt}}$ , using Eq. (66). The optimal statistic can then be used to make detection or upper limit statements (see Section IV).

First the power spectra spectra for each pulsar (and each stretch), $P_{i,k}(|f|)$ , must be determined. The spectra can either be modeled or measured with the methods described in Section VI. Then the overlap reduction functions, $^{(ij)}\Gamma_0$ , need to be computed for each pulsar pair. To optimize the statistic for particular spectra the value of $\alpha$ (in $\Omega_{\text{gw}}(f) = \Omega_{\alpha} f^{\alpha}$ ) needs to be chosen. The normalizations, $^{(ij)}\chi_k$ , can then be computed using Eq. (54). The normalizations allow us to compute the variances, $^{(ij)}\sigma^{-2}$ , given by Eq. (58), in the numerator and denominator of Eq. (66), as well as the filters, $^{(ij)}Q_k(f)$ , through Eq. (53). Note that the unknown factors of $\Omega_{\alpha}$ cancel everywhere: From Eq. (54) it is easy to see the normalization $^{(ij)}\chi_k \propto \Omega_{\alpha}^{-1}$ , so there is a cancellation a factor of $\Omega_{\alpha}$ in Eq. (53), and a factor of $\Omega_{\alpha}^2$ in Eq. (58). With these quantities in hand the cross-correlations, $^{(ij)}S_k$ , in Eq. (52) can be computed by taking Fourier transforms of the data (see Section VI). Alternatively, a set of time-domain filters, $^{(ij)}Q_k(t)$ , can be created by taking inverse Fourier transforms of Eq. (53) and applied to the data in the time domain using Eq. (26).

Note that there is no dependence on the arbitrary constant $T_0$ introduced in Eq. (51) for dimensional reasons. The $^{(ij)}\chi_k$ are linear in $T_0$ and enter the variances quadratically (see Eq. (58)). The dependence cancels in Eq. (66) because it is present in both numerator and denominator. $T_0$ also enters $S_{\text{opt}}$ linearly through $^{(ij)}\chi_k$ in

$^{(ij)}Q_k$ but cancels in Eq. (67) so that the point estimate of $\hat{\Omega}_{\alpha}$ is independent of $T_0$ .

D. Likelihood approach

The detection statistic that has been derived is also an optimal statistic in the sense that it is the logarithm of the likelihood ratio, at least in the limit where the expected signal is smaller than the noise, and therefore it is the optimal statistic in both the Bayesian sense and by the Neyman-Pearson criterion. This section is based on the likelihood analysis of [47], generalized to consider multiple detector pairs.

As we did previously, we assume that the noise is stationary and Gaussian, as is the stochastic background. For any given pulsar $i$ we assume that there are discrete samples of data which forms a vector $\mathbf{s}_i$ . Although the discussion below does not place requirements on the data sampling, we will assume that the observations of the pulsars all involve the same number of points $N$ at the same evenly spaced sampling interval so that sample $j$ of pulsar $i$ is $\mathbf{s}_i[j] = s_i(j\Delta t)$ where $\Delta t$ is the sampling interval. This signal vector is the sum of a noise vector $\mathbf{n}_i$ and the redshift vector $\mathbf{z}_i$ , $\mathbf{s}_i = \mathbf{z}_i + \mathbf{n}_i$ . The data is a combination of two random processes: the instrumental noise and the contribution from the stochastic background. The auto-correlation matrix $\mathbf{R}_i = \langle \mathbf{s}_i^{\dagger} \otimes \mathbf{s}i \rangle$ is an $N \times N$ matrix which contains both of these contributions and, since we assume Gaussian noise and stochastic background, this matrix completely characterizes the distribution of the data. As we did previously, we assume that the measurement noise in a pulsar observation is independent of the noise in the observations of other pulsars; the stochastic background, however, is correlated amongst the pulsar signals. This correlation is characterized by the stochastic background correlation matrix $\epsilon^2 \mathbf{S}{ij} = \langle \mathbf{z}_i^{\dagger} \otimes \mathbf{z}_j \rangle$ . Here $\epsilon$ is an order parameter which we will use to expand the probability distribution in powers of the small stochastic background signal. It can also be interpreted as an overall amplitude parameter of the stochastic background. The probability distribution for the collection of all pulsar observations is given by a multidimensional Gaussian distribution

$$p(\mathbf{x}\mathrm{\mid }ϵ)=\frac{1}{\sqrt{\det (2\pi \mathbf{\Sigma })}}\exp (-\frac{1}{2}{\mathbf{x}}^{†}\cdot {\mathbf{\Sigma }}^{-1}\cdot \mathbf{x})(69)p(\backslash mathbf\{x\}|\backslash epsilon)\; =\; \backslash frac\{1\}\{\backslash sqrt\{\backslash det(2\backslash pi\backslash mathbf\{\backslash Sigma\})\}\}\; \backslash exp\backslash left(-\backslash frac\{1\}\{2\}\backslash mathbf\{x\}^\{\backslash dagger\}\; \backslash cdot\; \backslash mathbf\{\backslash Sigma\}^\{-1\}\; \backslash cdot\; \backslash mathbf\{x\}\backslash right)\; \backslash quad\; (69)$$

where

$$\mathbf{x}=\left[\begin{array}{c}{\mathbf{s}}_1\\ {\mathbf{s}}_2\\ \mathrm{⋮}\\ {\mathbf{s}}_l\end{array}\right](70)\backslash mathbf\{x\}\; =\; \backslash begin\{bmatrix\}\; \backslash mathbf\{s\}\_1\; \backslash \backslash \; \backslash mathbf\{s\}\_2\; \backslash \backslash \; \backslash vdots\; \backslash \backslash \; \backslash mathbf\{s\}\_l\; \backslash end\{bmatrix\}\; \backslash quad\; (70)$$

is a column vector formed from all of the data vectorsand

is the correlation matrix for the collective observation vector $\mathbf{x}$ . In this weak signal limit we find

and

The logarithm of the likelihood ratio is

This is the optimal detection statistic for a weak stochastic background. We have identified $\mathcal{S}$ as the $O(\epsilon^2)$ term and $-2\mathcal{N}^2$ as the $O(\epsilon^4)$ term of the log-likelihood ratio.

The locally optimal detection statistic is obtained in the $\epsilon \rightarrow 0$ limit; it is the leading $O(\epsilon^2)$ term:

Although this presentation has been described in terms of observational vectors in the time-domain, the derivation of the likelihood ratio has not explicitly required this choice of basis. It is convenient to perform a unitary transformation that diagonalizes the various correlation matrices. This transformation is called a Karhunen-Loeve transformation; for a stationary process with a correlation time much shorter than the time spanned by the $l$ samples, the linear combinations of the time seriesthat diagonalize the correlation matrices asymptotically approach the discrete Fourier transform. Therefore we can approximately express our result in the frequency domain where the $\mathbf{R}i$ and $\mathbf{S}{ij}$ matrices can be understood in terms of the power spectrum and the expectation value of the redshift cross-correlation respectively [cf. Eq. (42) and Eq. (37)]. The locally-optimal detection statistic is therefore

where $\Omega_{\text{gw}}(f) = \epsilon^2 \hat{\Omega}{\text{gw}}(f)$ and $\Omega_\alpha = \epsilon^2 \hat{\Omega}\alpha$ . This is the same optimal detection statistic $S{\text{opt}}$ of Eq. (66) (with the simplification of $n{ij} = 1$ ) up to a normalization constant.

The locally optimal statistic is optimal in the limit of weak signals. However, the likelihood ratio is only determined by this statistic up to a unknown factor which depends on the (unknown) strength of the signal. It is important now to distinguish between the assumed amplitude of the stochastic background, $\epsilon$ , and the true amplitude, $\epsilon_{\text{true}}$ . The true gravitational wave spectrum $\Omega_{\text{gw}}(f)$ is now related to the template spectrum $\hat{\Omega}{\text{gw}}(f)$ via $\Omega{\text{gw}}(f) = \epsilon_{\text{true}}^2 \hat{\Omega}{\text{gw}}(f)$ . To measure the strength of the stochastic background given a set of pulsar observations, we can use the maximum likelihood estimator: the value of $\epsilon$ , $\epsilon{\text{MLE}}$ , for which the likelihood ratio is a maximum. That is, we wish to find the value of $\epsilon_{\text{MLE}}$ for which $d \ln \Lambda / d\epsilon^2|{\epsilon{\text{MLE}}} = 0$ . From Eq. (74) we see that this estimate is

$${ϵ}_{\text{MLE}}^2={\mathcal{N}}^{-2}\mathcal{S}(77)\backslash epsilon\_\{\backslash text\{MLE\}\}^2\; =\; \backslash mathcal\{N\}^\{-2\}\; \backslash mathcal\{S\}\; \backslash quad\; (77)$$

where $\mathcal{N}^2$ , from the $O(\epsilon^4)$ term of the log-likelihood ratio, is a normalizing factor which also includes the data. By substituting Eq. (77) into Eq. (74) we obtain the maximum likelihood detection statistic

$$\underset{ϵ}{\max }\ln \mathrm{\Lambda }\simeq \frac{1}{2}\frac{{\mathcal{S}}^2}{{\mathcal{N}}^2}(78)\backslash max\_\{\backslash epsilon\}\; \backslash ln\; \backslash Lambda\; \backslash simeq\; \backslash frac\{1\}\{2\}\; \backslash frac\{\backslash mathcal\{S\}^2\}\{\backslash mathcal\{N\}^2\}\; \backslash quad\; (78)$$

where the terms of $O(\epsilon^6)$ have been discarded. Notice that this statistic is not simply the square of the cross correlation statistic. The data also appears in the factor $\mathcal{N}^{-2}$ . This factor effectively suppresses elements of the pulsar network where the data measured greatly exceeds the normal noise level.

Some insight into the maximum likelihood detection statistic and the maximum likelihood amplitude estimate can be obtained by computing the expectation value of the log-likelihood ratio, Eq. (74). We find, to leading order in $\epsilon$ ,

Therefore

If we ignore the $O(\epsilon^6)$ terms, this is maximized when $\epsilon_{\text{MLE}} = \epsilon_{\text{true}}$ , in which case

This gives a scale of the value of the likelihood ratio we would expect to achieve.

IV. UPPER LIMITS AND DETECTION

Several methods exist in the LIGO literature that are appropriate for upper limit computation and detection using pulsar timing data [40, 41, 42, 43, 44, 45]. These methods can be divided into two classes: Frequentist and Bayesian.

We expect that the optimal statistic Eq. (66) will be formed from a large number of pulsar pairs. For example, the Parkes Pulsar Timing Array [35, 36] consists of 20 pulsars and the optimal statistic could be constructed from up to 190 cross-correlation pairs. In this case we can make use of the central limit theorem: The distribution of $\hat{S}{\text{opt}}$ should be well approximated by a Gaussian with a mean $\mu = \langle S{\text{opt}} \rangle = \Omega_\alpha T_0$ and variance $\sigma_{\hat{\mu}}^2$ given by Eq. (68), namely,

$$p({\widehat{S}}_{\text{opt}}\mathrm{\mid }\mu {\sigma }_{\widehat{\mu }})=\frac{1}{{\sigma }_{\widehat{\mu }}\sqrt{2\pi }}\exp \left(\frac{-({\widehat{S}}_{\text{opt}}-\mu {)}^2}{2{\sigma }_{\widehat{\mu }}^2}\right)\mathrm{.}(82)p(\backslash hat\{S\}\_\{\backslash text\{opt\}\}\; |\; \backslash mu\; \backslash sigma\_\{\backslash hat\{\backslash mu\}\})\; =\; \backslash frac\{1\}\{\backslash sigma\_\{\backslash hat\{\backslash mu\}\}\; \backslash sqrt\{2\backslash pi\}\}\; \backslash exp\; \backslash left(\; \backslash frac\{-(\backslash hat\{S\}\_\{\backslash text\{opt\}\}\; -\; \backslash mu)^2\}\{2\backslash sigma\_\{\backslash hat\{\backslash mu\}\}^2\}\; \backslash right).\; \backslash quad\; (82)$$

A straightforward frequentist upper limit can then be set by finding the value of $\mu_{\text{ul}}$ such that in some pre-determined fraction $C$ (called the confidence) of hypothetical experiments, the value of the optimal statistic exceeds the actual value $\hat{S}{\text{opt}}$ found in the search. In other words we would like to find the value $\mu{\text{ul}}$ such that

$${\int }_{{\widehat{S}}_{\text{opt}}}^{\mathrm{\infty }}d{S}_{\text{opt}}p(S_{\text{opt}}\mathrm{\mid }{\mu }_{\text{ul}}{\sigma }_{\widehat{\mu }})=C\mathrm{.}(83)\backslash int\_\{\backslash hat\{S\}\_\{\backslash text\{opt\}\}\}^\{\backslash infty\}\; dS\_\{\backslash text\{opt\}\}\; p(S\_\{\backslash text\{opt\}\}\; |\; \backslash mu\_\{\backslash text\{ul\}\}\; \backslash sigma\_\{\backslash hat\{\backslash mu\}\})\; =\; C.\; \backslash quad\; (83)$$

The solution to this is

$${\mu }_{\text{ul}}={\widehat{S}}_{\text{opt}}+\sqrt{2}{\sigma }_{\widehat{\mu }}{\text{erfc}}^{-1}(2(1-C))\mathrm{.}(84)\backslash mu\_\{\backslash text\{ul\}\}\; =\; \backslash hat\{S\}\_\{\backslash text\{opt\}\}\; +\; \backslash sqrt\{2\}\; \backslash sigma\_\{\backslash hat\{\backslash mu\}\}\; \backslash text\{erfc\}^\{-1\}(2(1-C)).\; \backslash quad\; (84)$$The assertion is that the real value of $\mu$ is less than $\mu_{\text{ul}}$ with confidence $C$ , because if $\mu = \mu_{\text{ul}}$ , a fraction $C$ of the time we would have observed a value of $S_{\text{opt}}$ greater than $\hat{S}{\text{opt}}$ . An equivalent, though potentially more robust, frequentist method to set upper limits involves performing simulated signal injections in the timing data set. Multiple injections are performed to determine the value of $\mu{\text{ul}}$ such that a fraction $C$ of the time the value of the optimal statistic measured in the data sets with injections exceeds the value found in the search. Frequentist detection methods such as Neyman-Pearson or maximum-likelihood are well described in the literature (see, for example, [40] and references therein) and we will not discuss them here. Additionally Feldman and Cousins [51] provide a means to smoothly transition between upper limits and detection.

Bayesian upper limits can be computed by constructing a posterior distribution using the value of the optimal statistic found in the search, and variance along with priors. We begin by applying the product rule to the probability density of $\mu$ along with the measured value $\hat{S}{\text{opt}}$ given $\sigma{\hat{\mu}}$ to write,

then solve for $p(\mu|\hat{S}{\text{opt}}\sigma{\hat{\mu}})$ to obtain Bayes' theorem,

$$p(\mu \mathrm{\mid }{\widehat{S}}_{\text{opt}}{\sigma }_{\widehat{\mu }})=p({\widehat{S}}_{\text{opt}}\mathrm{\mid }\mu {\sigma }_{\widehat{\mu }})\frac{p(\mu \mathrm{\mid }{\sigma }_{\widehat{\mu }})}{p({\widehat{S}}_{\text{opt}}\mathrm{\mid }{\sigma }_{\widehat{\mu }})},(86)p(\backslash mu|\backslash hat\{S\}\_\{\backslash text\{opt\}\}\backslash sigma\_\{\backslash hat\{\backslash mu\}\})\; =\; p(\backslash hat\{S\}\_\{\backslash text\{opt\}\}|\backslash mu\backslash sigma\_\{\backslash hat\{\backslash mu\}\})\; \backslash frac\{p(\backslash mu|\backslash sigma\_\{\backslash hat\{\backslash mu\}\})\}\{p(\backslash hat\{S\}\_\{\backslash text\{opt\}\}|\backslash sigma\_\{\backslash hat\{\backslash mu\}\})\},\; \backslash quad\; (86)$$

the posterior probability density for $\mu$ , or equivalently $\Omega_{\alpha}$ . One can then choose a prior $p(\mu|\sigma_{\hat{\mu}})$ (for example requiring $\mu > 0$ ) and normalize the probability distribution (the probability $p(\hat{S}{\text{opt}}|\sigma{\hat{\mu}})$ does not depend on $\mu$ so it is a prior dependent normalization constant), and find the $\mu_{\text{ul}}$ such that

$$\mathcal{M}{\int }_{-\mathrm{\infty }}^{{\mu }_{\text{ul}}}d\mu p(\mu \mathrm{\mid }{\widehat{S}}_{\text{opt}}{\sigma }_{\widehat{\mu }})=C,(87)\backslash mathcal\{M\}\; \backslash int\_\{-\backslash infty\}^\{\backslash mu\_\{\backslash text\{ul\}\}\}\; d\backslash mu\; p(\backslash mu|\backslash hat\{S\}\_\{\backslash text\{opt\}\}\backslash sigma\_\{\backslash hat\{\backslash mu\}\})\; =\; C,\; \backslash quad\; (87)$$

where $\mathcal{M}$ is the normalization constant. For sufficiently simple choices of the prior distribution $p(\mu|\sigma_{\hat{\mu}})$ the integral Eq. (87) can be performed analytically to obtain the Bayesian analog of Eq. (84). As with frequentist methods [40, 51], Bayesian detection methods involve selecting thresholds, in this case on the odds ratio, which is the ratio of the posteriors, suitably integrated over, say, different ranges of $\mu$ . For more details see Refs. [52, 53].

V. A PULSAR TIMING RADIOMETER: CONSTRUCTING SKYMAPS OF THE STOCHASTIC BACKGROUND

A sky map may be created by computing $\Omega_{\text{gw}}(f)$ for a collection of pixels in the sky. We do this by assuming that the only signal present comes from a single location

FIG. 2: Plots of $|^{(ij)}\Gamma_{\hat{\Omega}}|$ from Eq. (92) for two pulsars with $\xi = \pi/2$ (top panel) and $\xi = \pi$ degrees (bottom panel)

on the sky. We begin by relaxing the assumption that the stochastic background is isotropic. That is, we take

where $P(\hat{\Omega})$ is the strength or brightness [54] of gravitational waves from the direction $\hat{\Omega}$ .

In this case, the overlap reduction function takes the modified form,

$${}^{(ij)}{\mathrm{\Gamma }}_P=\frac{3}{4\pi }\sum _A{\int }_{S^2}d\widehat{\mathrm{\Omega }}P(\widehat{\mathrm{\Omega }})F_i^A(\widehat{\mathrm{\Omega }})F_j^A(\widehat{\mathrm{\Omega }})\mathrm{.}(89)^\{(ij)\}\backslash Gamma\_P\; =\; \backslash frac\{3\}\{4\backslash pi\}\; \backslash sum\_A\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; P(\backslash hat\{\backslash Omega\})\; F\_i^A(\backslash hat\{\backslash Omega\})\; F\_j^A(\backslash hat\{\backslash Omega\}).\; \backslash quad\; (89)$$

where we've ignored the pulsar phase factors, and the optimal filter is given by

$${}^{(ij)}{\tilde{Q}}_P(f){=}^{(ij)}\chi \frac{{}^{(ij)}{\mathrm{\Gamma }}_P{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid })}{\mathrm{\mid }f{\mathrm{\mid }}^3{P}_i(\mathrm{\mid }f\mathrm{\mid })P_j(\mathrm{\mid }f\mathrm{\mid })},(90)^\{(ij)\}\backslash tilde\{Q\}\_P(f)\; =\; ^\{(ij)\}\backslash chi\; \backslash frac\{^\{(ij)\}\backslash Gamma\_P\; \backslash Omega\_\{\backslash text\{gw\}\}(|f|)\}\{|f|^3\; P\_i(|f|)\; P\_j(|f|)\},\; \backslash quad\; (90)$$

where we have suppressed the $k$ index which specifies the particular measurement of the $(ij)$ pulsar pair. We can further optimize for point sources by taking $P(\hat{\Omega}) = \delta^2(\hat{\Omega} - \hat{\Omega}')$ . The optimal filter then becomes,

$${}^{(ij)}{\tilde{Q}}_{\widehat{\mathrm{\Omega }}}(f){=}^{(ij)}\chi \frac{{}^{(ij)}{\mathrm{\Gamma }}_{\widehat{\mathrm{\Omega }}}{\mathrm{\Omega }}_{\text{gw}}(\mathrm{\mid }f\mathrm{\mid })}{\mathrm{\mid }f{\mathrm{\mid }}^3{P}_i(\mathrm{\mid }f\mathrm{\mid })P_j(\mathrm{\mid }f\mathrm{\mid })},(91)^\{(ij)\}\backslash tilde\{Q\}\_\{\backslash hat\{\backslash Omega\}\}(f)\; =\; ^\{(ij)\}\backslash chi\; \backslash frac\{^\{(ij)\}\backslash Gamma\_\{\backslash hat\{\backslash Omega\}\}\; \backslash Omega\_\{\backslash text\{gw\}\}(|f|)\}\{|f|^3\; P\_i(|f|)\; P\_j(|f|)\},\; \backslash quad\; (91)$$with,

$${}^{(ij)}{\mathrm{\Gamma }}_{\widehat{\mathrm{\Omega }}}=\frac{3}{4\pi }\sum _A{F}_i^A(\widehat{\mathrm{\Omega }})F_j^A(\widehat{\mathrm{\Omega }})\mathrm{.}(92)\{\}^\{(ij)\}\backslash Gamma\_\{\backslash hat\{\backslash Omega\}\}\; =\; \backslash frac\{3\}\{4\backslash pi\}\; \backslash sum\_A\; F\_i^A(\backslash hat\{\backslash Omega\})\; F\_j^A(\backslash hat\{\backslash Omega\}).\; \backslash quad\; (92)$$

Figure 2 shows two examples of the sky location dependent overlap reduction function. The top panel shows $|\Gamma_{\hat{\Omega}}|$ from Eq. (92) for two pulsars with $\xi = \pi/2$ located at $0^\circ$ Dec and 9h and 15h RA respectively. The bottom panel shows the same quantity for two pulsars with $\xi = \pi$ located at $0^\circ$ Dec and 6h and 18h RA respectively.

One could also imagine computing the overlap reduction function for each term in a multipole expansion of $P(\hat{\Omega})$ . The overlap reduction function for the monopole term in the expansion (appropriate for an isotropic stochastic gravitational wave search) is the Hellings-Downs curve given by Eq. (39) in Section III. Surprisingly, the dipole overlap reduction function is given by a similarly simple equation. We find

$${\mathrm{\Gamma }}_{\text{dip}}=(\cos {\alpha }_1+\cos {\alpha }_2)\times (2-\frac{3}{2}\cos \xi +6{\tan }^2\frac{\xi }{2}\ln (\sin \frac{\xi }{2})),(93)\backslash Gamma\_\{\backslash text\{dip\}\}\; =\; (\backslash cos\; \backslash alpha\_1\; +\; \backslash cos\; \backslash alpha\_2)\; \backslash times\; \backslash left(\; 2\; -\; \backslash frac\{3\}\{2\}\; \backslash cos\; \backslash xi\; +\; 6\; \backslash tan^2\; \backslash frac\{\backslash xi\}\{2\}\; \backslash ln\; \backslash left(\; \backslash sin\; \backslash frac\{\backslash xi\}\{2\}\; \backslash right)\; \backslash right),\; \backslash quad\; (93)$$

where as before $\xi$ is the angle between the two pulsars, and $\alpha_1$ and $\alpha_2$ are the angles each of two pulsars make to the direction of the dipole. A detailed derivation of this result is given in Appendix C2. This result is relevant to searches for a dipole anisotropy in the gravitational wave sky using pulsar timing data.

The sky-dependence of the sensitivity of a pulsar network can be estimated by computing the signal to noise for sources at the sky locations of interest. We start by taking the expectation value of the optimal statistic, Eq. (66), using the optimal filter for a sky location $\hat{\Omega}$ assuming the redshift data contain a stochastic signal from that location. We then divide by the square root of the variance given in Eq. (68). The result is proportional to

where we have assumed (for illustrative purposes) that the noise spectra of all pulsars is the same, the observation times for all pairs is the same, and $n_{ij} = 1$ for all pulsar pairs. Figure 3 shows the Eq. (94), for the 20 pulsars of the Parkes Pulsar Timing Array [35, 36]. Since most of the pulsars are in the Southern hemisphere the Parkes Pulsar Timing Array is most sensitive in that region.

Another quantity of interest is the point spread function, which measures the intrinsic spatial correlation of the skymap, or equivalently, the ability of a pulsar network to locate a stochastic source of gravitational waves. We construct the point spread function by computing the signal to noise for a source at some sky location that we search for using the optimal filter for some other location. In particular, we take the expectation value of

FIG. 3: Skymap of the sensitivity, Eq. (94), for the Parkes Pulsar Timing Array.

the optimal statistic, Eq. (66), using the optimal filter for a sky location $\hat{\Omega}$ assuming the redshift data contain a signal from another location $\hat{\Omega}'$ , then we divide by the square root of the variance given in Eq. (68). The result is proportional to

where we have again assumed that the noise spectra of all pulsars is the same, the observation times for all pairs is the same, and $n_{ij} = 1$ for all pulsar pairs. Figure 4 shows the point spread function, Eq. (95), for the 20 pulsars of the Parkes Pulsar Timing Array [35, 36] for a source at 6h RA $45^\circ$ Dec (top panel) and another at 18h RA $-45^\circ$ Dec (bottom panel).

The point spread function can be understood in terms of the likelihood ratio of Sec. IIID: Suppose that the likelihood ratio is computed using the overlap reduction function $\Gamma_{\hat{\Omega}}$ appropriate for a stochastic signal coming from direction $\Gamma_{\hat{\Omega}}$ when the true signal is in fact coming from direction $\Gamma_{\hat{\Omega}'}$ . The the expectation value of the log-likelihood ratio is [cf. Eq. (80)]

with

$${}^{(ij)}C={\left(\frac{3H_0}{32{\pi }^3}\right)}^{2}T{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}df\frac{{\widehat{\mathrm{\Omega }}}_{\text{gw}}^2(\mathrm{\mid }f\mathrm{\mid })}{f^6{P}_i(\mathrm{\mid }f\mathrm{\mid })P_j(\mathrm{\mid }f\mathrm{\mid })}\mathrm{.}(97)\{\}^\{(ij)\}C\; =\; \backslash left(\; \backslash frac\{3H\_0\}\{32\backslash pi^3\}\; \backslash right)^2\; T\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; df\; \backslash frac\{\backslash hat\{\backslash Omega\}\_\{\backslash text\{gw\}\}^2(|f|)\}\{f^6\; P\_i(|f|)\; P\_j(|f|)\}.\; \backslash quad\; (97)$$

If ${}^{(ij)}C$ is approximately the same for all pulsar pairsFIG. 4: Plot of the point spread function Eq. (95) for the Parkes pulsar timing array for a source at 6h RA 45° Dec (top panel) and 18h RA -45° Dec (bottom panel).

then

$$\underset{ϵ}{\max }⟨\ln \mathrm{\Lambda }⟩\propto A^2(\widehat{\mathrm{\Omega }},{\widehat{\mathrm{\Omega }}}^{\mathrm{\prime }})\mathrm{.}(98)\backslash max\_\{\backslash epsilon\}\; \backslash langle\; \backslash ln\; \backslash Lambda\; \backslash rangle\; \backslash propto\; A^2(\backslash hat\{\backslash Omega\},\; \backslash hat\{\backslash Omega\}\text{'}).\; \backslash quad\; (98)$$

In this sense the point spread function describes the degree to which the position of a point source of stochastic gravitational waves can be located in terms of the likelihood ratio.

VI. ISSUES WITH PULSAR TIMING DATA

We have derived the optimal statistic for detecting a stochastic background. In this section we would like to discuss some issues associated with departures from the idealizations made to arrive at the optimal statistic.

A. Colored noise and non-stationarity

In contrast to previous methods [7, 8] the techniques presented here do not rely on the data being white. The power spectra $P_i(|f|)$ in the optimal statistic account for colored noise. However, the methods assume the data is stationary.

If the data is non-stationary over long timescales it can be divided into short stationary (or almost stationary) stretches and the power spectrum can be estimated

using the Lomb-Scargle periodogram described below, or modeled for each stretch. The optimally filtered data stretches can then be combined along the lines discussed in Section III B. One concern associated with breaking the data up into small stretches is loss of low frequency information: Gravitational waves with periods larger than the length of the short stretches will be lost in this procedure. The problem can be avoided by first computing the quantities $s_i(f)/P_i(|f|)$ for each of the short stretches and then combine them using the Dirichlet kernel to construct full time baseline versions of these quantities.

If the spectrum is measured it can be smoothed by performing a running average over a small frequency window, which if the data are stationary in the stretch the spectrum is estimated, is equivalent to ensemble averaging.

B. Unevenly sampled data

The fact that pulsar timing measurements are not taken continuously leads to a data set that is unevenly sampled in time. This poses a problem for frequency-domain analyses not present in their time-domain counterparts. The authors of [46] will explore the time domain approach in detail. It is unclear which of these two methods will turn out to be more robust. In what follows we address the specific issues of computing periodograms and Fourier transforms for unevenly sampled data sets which we think is useful in any case.

1. The Lomb-Scargle Periodogram

The problem of constructing periodograms from unevenly sampled data comes up in the data analysis of variable stars. It was in precisely this context that Lomb [55] and Scargle [56] proposed a least-squares solution to the problem. The basic idea is as follows: Let $x(t_i)$ be a time series with zero mean sampled at $i = 1 \dots N$ unevenly spaced times. Now fit the time series by finding the coefficients $a_{\min}$ and $b_{\min}$ that minimize the square of the residual

$$r^2(f)\equiv \sum _{i=1}^N\{x(t_i)-a\cos [2\pi f(t_i-\tau )]-b\sin [2\pi f(t_i-\tau )]{\}}^2,(99)r^2(f)\; \backslash equiv\; \backslash sum\_\{i=1\}^N\; \backslash \{x(t\_i)\; -\; a\; \backslash cos[2\backslash pi\; f(t\_i\; -\; \backslash tau)]\; -\; b\; \backslash sin[2\backslash pi\; f(t\_i\; -\; \backslash tau)]\backslash \}^2,\; \backslash quad\; (99)$$

where

$$\tan (4\pi f\tau )=\frac{\sum _{i=1}^N\sin (4\pi f{t}_i)}{\sum _{i=1}^N\cos (4\pi f{t}_i)}\mathrm{.}(100)\backslash tan(4\backslash pi\; f\backslash tau)\; =\; \backslash frac\{\backslash sum\_\{i=1\}^N\; \backslash sin(4\backslash pi\; f\; t\_i)\}\{\backslash sum\_\{i=1\}^N\; \backslash cos(4\backslash pi\; f\; t\_i)\}.\; \backslash quad\; (100)$$Then the periodogram is defined up to normalization by the difference

where $r_{\min}^2(f)$ is the quantity in Eq. (99) with $a = a_{\min}$ and $b = b_{\min}$ . After normalization [57] and generalization to data with nonzero mean, we have the Lomb-Scargle periodogram

$$P_X^{\text{LS}}(f)=\frac{1}{2{\sigma }_x^2}\{\frac{{(\sum _{i=1}^N[x(t_i)-{\mu }_x]\cos [2\pi f(t_i-\tau )])}^2}{\sum _{i=1}^N{\cos }^2[2\pi f(t_i-\tau )]}+\frac{{(\sum _{i=1}^N[x(t_i)-{\mu }_x]\sin [2\pi f(t_i-\tau )])}^2}{\sum _{i=1}^N{\sin }^2[2\pi f(t_i-\tau )]}\},(102)P\_X^\{\backslash text\{LS\}\}(f)\; =\; \backslash frac\{1\}\{2\backslash sigma\_x^2\}\; \backslash left\backslash \{\; \backslash frac\{\backslash left(\backslash sum\_\{i=1\}^N\; [x(t\_i)\; -\; \backslash mu\_x]\; \backslash cos[2\backslash pi\; f(t\_i\; -\; \backslash tau)]\backslash right)^2\}\{\backslash sum\_\{i=1\}^N\; \backslash cos^2[2\backslash pi\; f(t\_i\; -\; \backslash tau)]\}\; +\; \backslash frac\{\backslash left(\backslash sum\_\{i=1\}^N\; [x(t\_i)\; -\; \backslash mu\_x]\; \backslash sin[2\backslash pi\; f(t\_i\; -\; \backslash tau)]\backslash right)^2\}\{\backslash sum\_\{i=1\}^N\; \backslash sin^2[2\backslash pi\; f(t\_i\; -\; \backslash tau)]\}\; \backslash right\backslash \},\; \backslash quad\; (102)$$

where $\mu_x$ and $\sigma_x^2$ are the mean and variance, respectively of $x(t_i)$ . Note that the definition of $\tau$ in Eq. (100) ensures that the resulting periodogram is independent of where $t = 0$ .

2. Fourier transforms

The idea of using a least-squares minimization is also useful for constructing Fourier transforms. To do so, we borrow an idea from the radar community [58]. Suppose we have a timeseries, $x(t_i)$ , non-uniformly sampled at times $t_0 \dots t_N$ and we wish to construct its Fourier transform,

$$\tilde{x}(f_m)=\sum _{j=0}^{N}x(t_j)e^{-i2\pi f_m{t}_j}(103)\backslash tilde\{x\}(f\_m)\; =\; \backslash sum\_\{j=0\}^N\; x(t\_j)\; e^\{-i2\backslash pi\; f\_m\; t\_j\}\; \backslash quad\; (103)$$

over $M$ evenly spaced frequencies, $f_0 \dots f_M$ . The strategy we will employ is to use a least-squares procedure to find the best fit to the original timeseries after an inverse Fourier transform. That is, the (squared) residual to be minimized is given by

$$r^2=\sum _{j=0}^N{\mid x(t_j)-\sum _{k=0}^M\tilde{\xi }(f_k)e^{i2\pi f_k{t}_j}\mid }^2,(104)r^2\; =\; \backslash sum\_\{j=0\}^N\; \backslash left|\; x(t\_j)\; -\; \backslash sum\_\{k=0\}^M\; \backslash tilde\{\backslash xi\}(f\_k)\; e^\{i2\backslash pi\; f\_k\; t\_j\}\; \backslash right|^2,\; \backslash quad\; (104)$$

where the $\tilde{\xi}(f_k)$ are to be determined. Defining

$$A_{kj}=e^{i2\pi f_k{t}_j},(105)A\_\{kj\}\; =\; e^\{i2\backslash pi\; f\_k\; t\_j\},\; \backslash quad\; (105)$$

we can write

$$r^2={\parallel \vec{x}-A\overrightarrow{\stackrel{~}{\xi }}\parallel }^2\mathrm{.}(106)r^2\; =\; \backslash left\backslash |\; \backslash vec\{x\}\; -\; A\; \backslash vec\{\backslash tilde\{\backslash xi\}\}\; \backslash right\backslash |^2.\; \backslash quad\; (106)$$

The least-squares solution to this problem is given by

$$\overrightarrow{\stackrel{~}{\xi }}=(A^{†}A{)}^{-1}{A}^{†}\vec{x}\mathrm{.}(107)\backslash vec\{\backslash tilde\{\backslash xi\}\}\; =\; (A^\backslash dagger\; A)^\{-1\}\; A^\backslash dagger\; \backslash vec\{x\}.\; \backslash quad\; (107)$$

The problem is then purely a computational one, which, because of the limited amount of pulsar timing data available, is completely tractable on a modern computer, regardless of the efficiency of the algorithm. On a final note, one can actually improve upon this procedure [58] by weighing the residual in Eq. (104) by the square root of the variance, $\sigma_{x(t_j)}$ , associated with each data point

$$r^2=\sum _{j=0}^N\frac{1}{{\sigma }_{x(t_j)}}{\mid x(t_j)-\sum _{k=0}^M{\tilde{\xi }}_k{e}^{i2\pi f_k{t}_j}\mid }^2,(108)r^2\; =\; \backslash sum\_\{j=0\}^N\; \backslash frac\{1\}\{\backslash sigma\_\{x(t\_j)\}\}\; \backslash left|\; x(t\_j)\; -\; \backslash sum\_\{k=0\}^M\; \backslash tilde\{\backslash xi\}\_k\; e^\{i2\backslash pi\; f\_k\; t\_j\}\; \backslash right|^2,\; \backslash quad\; (108)$$

which has the advantage of automatically including the error bars associated with individual pulsar timing data points.

VII. SUMMARY AND CONCLUSIONS

A stochastic background of gravitational waves could be detected via pulsar timing observations in the next 5 to 10 years. This background may be astrophysical, such as that produced by supermassive black holes, or cosmological, such as that produced by a network of cosmic (super)strings. In the latter case a detection would open a window onto a time in the early universe prior to recombination and could have profound consequences. Leveraging techniques developed for ground-based instruments such as LIGO and Virgo, in this paper we have shown how to optimally extract the signal produced by a stochastic background of gravitational waves using cross-correlations of timing data from a pulsar timing array.We started by considering the redshift induced by a gravitational wave on the frequency of arrival of radio pulses from a pulsar first derived by Detweiler [3]. The redshift is proportional to the difference in the metric perturbation at the pulsar (when a pulse is emitted) and at the Earth (when that pulse is received). Using a convenient coordinate independent description of the signal we examined the form of the signal in the frequency domain. The term involving the metric perturbation at the pulsar is typically neglected because it can be treated as a sort of noise term which averages to zero in correlations of timing measurements of different pulsars. In the frequency domain the dependence on the metric perturbation at the pulsar is in a phase factor that depends on the distance to the pulsar. It is possible that if we could determine the distance to pulsars with sufficient accuracy we could use the metric perturbation at the pulsar to improve the sensitivity of our searches. Unfortunately, accurate measurements of pulsar distances are unavailable. By first finding the optimal cross-correlation filter, we have shown that for pulsar distances and gravitational wave frequencies typical of pulsar timing experiments, the metric perturbation at the pulsar can be neglected without a significant deviation from optimality. It is unclear whether this is true for other types of gravitational wave searches. We have also determined the optimal way to combine pulsar timing data from a pulsar timing array, which is constructed from pairs of optimally filtered cross-correlations.

We have discussed and illustrated frequentist and Bayesian methods for setting upper limits using the distribution of the optimal statistic. We have shown how to construct a pulsar timing radiometer: A map of the sky created by optimizing the cross-correlation statistic for particular sky directions. We have also shown how to determine the intrinsic spatial correlation of such maps, which in turn determines the ability of a pulsar timing array to locate a source of stochastic gravitational waves.

We have ended with a discussion of some problems related to realistic analysis of pulsar timing data, particularly the issues of non-stationarity and uneven sampling. The optimal filter is constructed from power spectra of the pulsar timing data, which can be modeled or measured, and accounts for the effects of colored noise. We have described a technique, the Lomb-Scargle periodogram, for robust spectrum estimation that can be used to construct the optimal filter. The optimal filter can then be applied in the frequency domain and we have described a procedure for taking Fourier transforms of unevenly sampled data that accounts for error bars in the individual pulsar timing data points. The optimal filter can also be inverse Fourier transformed and applied in the time domain where uneven sampling is not an issue. Regardless of which method turns out to be more useful and robust for stochastic background searches, we believe the development of Fourier techniques for unevenly sampled data will be beneficial. Lommen, Romano and Woan [46] will examine time-domain methods in detail.

Acknowledgments

We would like to thank Luis Anchordoqui, Steven Detweiler, Nick Fotopoulos, Rick Jenet, Adam Mercer, and Joe Romano for many useful discussions. Additionally we would like to thank Nick Fotopoulos and Adam Mercer for carefully reading the manuscript. JC is supported in part by NSF Grant No. PHY-0701817. LP is supported by NSF Grant No. PHY-0503366 and the Research Growth Initiative at the University of Wisconsin-Milwaukee. XS is supported in part by NSF Grant No. PHY-0758155. SB gratefully acknowledges the support of the California Institute of Technology. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-0107417.

APPENDIX A: DERIVATION OF DETWEILER'S FORMULA USING THE GEODESIC EQUATION

For completeness of presentation, we include a derivation of Detweiler's formula. We consider, as we did in Section II, the metric perturbation due to a single gravitational wave traveling in the $\hat{\Omega} = \hat{z}$ direction, so Eqs. (4-6) hold. Then, if a vector $s^a$ is null in Minkowski space, the corresponding null vector, $\sigma^a$ , in the perturbed spacetime $g_{ab} = \eta_{ab} + h_{ab}$ is given by [59],

$${\sigma }^a=s^a-\frac{1}{2}{\eta }^{ab}{h}_{bc}{s}^c\mathrm{.}(A1)\backslash sigma^a\; =\; s^a\; -\; \backslash frac\{1\}\{2\}\backslash eta^\{ab\}h\_\{bc\}s^c.\; \backslash quad\; (A1)$$

The null vector in Minkowski space that points from the pulsar to the solar system is $s^a = \nu(1, -\alpha, -\beta, -\gamma)$ , where, as before, $\alpha$ , $\beta$ and $\gamma$ are the direction cosines with the $x$ -, $y$ - and $z$ -directions, respectively. The corresponding perturbed vector is readily computed from Eq. (A1)

$${\sigma }^a=\nu \left(\begin{array}{c}1\\ -\alpha (1-\frac{1}{2}{h}_{+})+\frac{1}{2}\beta h_{\times }\\ -\beta (1+\frac{1}{2}{h}_{+})+\frac{1}{2}\alpha h_{\times }\\ -\gamma \end{array}\right)\mathrm{.}(A2)\backslash sigma^a\; =\; \backslash nu\; \backslash begin\{pmatrix\}\; 1\; \backslash \backslash \; -\backslash alpha(1\; -\; \backslash frac\{1\}\{2\}h\_+)\; +\; \backslash frac\{1\}\{2\}\backslash beta\; h\_\backslash times\; \backslash \backslash \; -\backslash beta(1\; +\; \backslash frac\{1\}\{2\}h\_+)\; +\; \backslash frac\{1\}\{2\}\backslash alpha\; h\_\backslash times\; \backslash \backslash \; -\backslash gamma\; \backslash end\{pmatrix\}.\; \backslash quad\; (A2)$$

The geodesic equation tells us that the $t$ -component of $\sigma^a$ satisfies

$$\frac{d{\sigma }^t}{d\lambda }=-{\mathrm{\Gamma }}_{ab}^t{\sigma }^a{\sigma }^b\mathrm{.}(A3)\backslash frac\{d\backslash sigma^t\}\{d\backslash lambda\}\; =\; -\backslash Gamma\_\{ab\}^t\; \backslash sigma^a\; \backslash sigma^b.\; \backslash quad\; (A3)$$

It follows from the form of the metric perturbation in Eq. (5) that

where the overdot denotes a derivative with respect to $t$ . The geodesic equation then reads

After a little algebra Eq. (A2) leads to

$$({\sigma }^x{)}^2-({\sigma }^y{)}^2={\nu }^2({\alpha }^2-{\beta }^2)+\mathcal{O}(h),(\text{A6})(\backslash sigma^x)^2\; -\; (\backslash sigma^y)^2\; =\; \backslash nu^2(\backslash alpha^2\; -\; \backslash beta^2)\; +\; \backslash mathcal\{O\}(h),\; \backslash quad\; (\backslash text\{A6\})$$

as well as

$${\sigma }^x{\sigma }^y={\nu }^2\alpha \beta +\mathcal{O}(h),(\text{A7})\backslash sigma^x\backslash sigma^y\; =\; \backslash nu^2\backslash alpha\backslash beta\; +\; \backslash mathcal\{O\}(h),\; \backslash quad\; (\backslash text\{A7\})$$

so that

$$-\frac{d\nu }{d\lambda }=\frac{1}{2}{\dot{h}}_{+}{\nu }^2({\alpha }^2-{\beta }^2)+{\dot{h}}_{\times }{\nu }^2\alpha \beta \mathrm{.}(\text{A8})-\backslash frac\{d\backslash nu\}\{d\backslash lambda\}\; =\; \backslash frac\{1\}\{2\}\backslash dot\{h\}\_+\backslash nu^2(\backslash alpha^2\; -\; \backslash beta^2)\; +\; \backslash dot\{h\}\_\backslash times\backslash nu^2\backslash alpha\backslash beta.\; \backslash quad\; (\backslash text\{A8\})$$

We proceed by writing the time derivatives in Eq. (A8) as derivatives with respect to the affine parameter $\lambda$ . In particular, since $h_{+, \times} = h_{+, \times}(t - z)$ we have

$$\frac{d{h}_{+,\times }}{d\lambda }=\frac{\mathrm{\partial }{h}_{+,\times }}{\mathrm{\partial }t}\frac{dt}{d\lambda }+\frac{\mathrm{\partial }{h}_{+,\times }}{\mathrm{\partial }z}\frac{dz}{d\lambda }\mathrm{.}(\text{A9})\backslash frac\{dh\_\{+,\; \backslash times\}\}\{d\backslash lambda\}\; =\; \backslash frac\{\backslash partial\; h\_\{+,\; \backslash times\}\}\{\backslash partial\; t\}\; \backslash frac\{dt\}\{d\backslash lambda\}\; +\; \backslash frac\{\backslash partial\; h\_\{+,\; \backslash times\}\}\{\backslash partial\; z\}\; \backslash frac\{dz\}\{d\backslash lambda\}.\; \backslash quad\; (\backslash text\{A9\})$$

Moreover, we also have that the frequency $\nu = dt/d\lambda$ , in addition to $\partial h_{+, \times}/\partial z = -\partial h_{+, \times}/\partial t$ , and $dz/d\lambda = -\nu\gamma$ . Therefore we can write Eq. (A9) as

$${\dot{h}}_{+,\times }=\frac{1}{\nu (1+\gamma )}\frac{d{h}_{+,\times }}{d\lambda }\mathrm{.}(\text{A10})\backslash dot\{h\}\_\{+,\; \backslash times\}\; =\; \backslash frac\{1\}\{\backslash nu(1\; +\; \backslash gamma)\}\; \backslash frac\{dh\_\{+,\; \backslash times\}\}\{d\backslash lambda\}.\; \backslash quad\; (\backslash text\{A10\})$$

Then Eq. (A8), the geodesic equation, becomes

$$-\frac{1}{\nu }\frac{d\nu }{d\lambda }=\frac{1}{2}\frac{{\alpha }^2-{\beta }^2}{1+\gamma }\frac{d{h}_{+}}{d\lambda }+\frac{\alpha \beta }{1+\gamma }\frac{d{h}_{\times }}{d\lambda },(\text{A11})-\backslash frac\{1\}\{\backslash nu\}\; \backslash frac\{d\backslash nu\}\{d\backslash lambda\}\; =\; \backslash frac\{1\}\{2\}\; \backslash frac\{\backslash alpha^2\; -\; \backslash beta^2\}\{1\; +\; \backslash gamma\}\; \backslash frac\{dh\_+\}\{d\backslash lambda\}\; +\; \backslash frac\{\backslash alpha\backslash beta\}\{1\; +\; \backslash gamma\}\; \backslash frac\{dh\_\backslash times\}\{d\backslash lambda\},\; \backslash quad\; (\backslash text\{A11\})$$

which we integrate to find

$$\frac{\nu (t)}{{\nu }_0}=\exp [-\frac{1}{2}\frac{{\alpha }^2-{\beta }^2}{1+\gamma }\mathrm{\Delta }{h}_{+}-\frac{\alpha \beta }{1+\gamma }\mathrm{\Delta }{h}_{\times }]\mathrm{.}(\text{A12})\backslash frac\{\backslash nu(t)\}\{\backslash nu\_0\}\; =\; \backslash exp\; \backslash left[\; -\backslash frac\{1\}\{2\}\; \backslash frac\{\backslash alpha^2\; -\; \backslash beta^2\}\{1\; +\; \backslash gamma\}\; \backslash Delta\; h\_+\; -\; \backslash frac\{\backslash alpha\backslash beta\}\{1\; +\; \backslash gamma\}\; \backslash Delta\; h\_\backslash times\; \backslash right].\; \backslash quad\; (\backslash text\{A12\})$$

It is worth pointing out that the direction cosines are functions of the affine parameter. The dependence is in terms of $\mathcal{O}(h)$ , and we have neglected this dependence in going from Eq. (A11) to Eq. (A12). The final result is obtained by expanding this expression to first order in $h_{+, \times}$

$$\frac{{\nu }_0-\nu (t)}{{\nu }_0}=\frac{1}{2}\frac{{\alpha }^2-{\beta }^2}{1+\gamma }\mathrm{\Delta }{h}_{+}+\frac{\alpha \beta }{1+\gamma }\mathrm{\Delta }{h}_{\times },(\text{A13})\backslash frac\{\backslash nu\_0\; -\; \backslash nu(t)\}\{\backslash nu\_0\}\; =\; \backslash frac\{1\}\{2\}\; \backslash frac\{\backslash alpha^2\; -\; \backslash beta^2\}\{1\; +\; \backslash gamma\}\; \backslash Delta\; h\_+\; +\; \backslash frac\{\backslash alpha\backslash beta\}\{1\; +\; \backslash gamma\}\; \backslash Delta\; h\_\backslash times,\; \backslash quad\; (\backslash text\{A13\})$$

where $\Delta h_{+, \times} = h_{+, \times}^{\text{P}} - h_{+, \times}^{\text{e}}$ is the difference between the metric perturbation at the pulsar and the detector. This expression is precisely Eq. (7).

APPENDIX B: LINEARITY OF THE REDSHIFT

In this appendix we will explicitly demonstrate that the total redshift is the sum of the contributions from gravitational waves in every direction, as written in Eq. (18). The derivation is merely a generalization of the results derived in the previous appendix. We begin by considering a metric perturbation, $h_{ab}$ , in a spatial transverse-traceless gauge comprised of the sum of metric perturbations, $h_{ab}^{(i)}$ , from gravitational waves in $N$ different directions, $\hat{\Omega}_{(i)}$ . Namely,

$$h_{ab}=\sum _{(i)}^N{h}_{ab}^{(i)}(t-{\widehat{\mathrm{\Omega }}}_{(i)}\cdot \vec{x}),(\text{B1})h\_\{ab\}\; =\; \backslash sum\_\{(i)\}^N\; h\_\{ab\}^\{(i)\}(t\; -\; \backslash hat\{\backslash Omega\}\_\{(i)\}\; \backslash cdot\; \backslash vec\{x\}),\; \backslash quad\; (\backslash text\{B1\})$$

where $t$ and $\vec{x}$ form a four vector, $x^a$ , in the background (Minkowski) geometry. Adjusting the notation from Appendix A,

as a null vector in Minkowski space. The null geodesic is perturbed by Eq. (B1), resulting in a

$${\sigma }^a=s^a+\delta s^a\mathrm{.}(\text{B3})\backslash sigma^a\; =\; s^a\; +\; \backslash delta\; s^a.\; \backslash quad\; (\backslash text\{B3\})$$

As before, our interest is in the quantity

$$\frac{d{\sigma }^t}{d\lambda }=-{\mathrm{\Gamma }}_{ab}^t{\sigma }^a{\sigma }^b\mathrm{.}(\text{B4})\backslash frac\{d\backslash sigma^t\}\{d\backslash lambda\}\; =\; -\backslash Gamma\_\{ab\}^t\; \backslash sigma^a\; \backslash sigma^b.\; \backslash quad\; (\backslash text\{B4\})$$

The spatial nature of the gauge we've chosen ensures that

$${\mathrm{\Gamma }}_{ab}^t=\frac{1}{2}{\dot{g}}_{ab}(\text{B5})\backslash Gamma\_\{ab\}^t\; =\; \backslash frac\{1\}\{2\}\; \backslash dot\{g\}\_\{ab\}\; \backslash quad\; (\backslash text\{B5\})$$

$$=\frac{1}{2}{\dot{h}}_{ab},(\text{B6})=\; \backslash frac\{1\}\{2\}\; \backslash dot\{h\}\_\{ab\},\; \backslash quad\; (\backslash text\{B6\})$$

which is evident from the first line of Eq. (A4). It follows from these definitions that

where $i$ and $j$ are spatial indices. As before, we want to write the expression above in terms of the affine parameter along $s^a$ . We begin by noting that for term inEq. (B1)

where we have used the $(t - \hat{\Omega}_{(i)} \cdot \vec{x})$ dependence of the metric perturbation to write the spatial derivatives as time derivatives along with Eq. (B2). Putting this together with Eq. (B7), the result is that

$$\begin{array}{cccc}& -\frac{1}{\nu }\frac{d\nu }{d\lambda }=\sum _{(i)}\frac{1}{2}\frac{{\widehat{p}}^i{\widehat{p}}^j}{1+{\widehat{\mathrm{\Omega }}}_{(i)}\cdot \widehat{p}}{h}_{ij}^{(i)}(t-{\widehat{\mathrm{\Omega }}}_{(i)}\cdot \vec{x}),& & \text{(B9)}\end{array}-\backslash frac\{1\}\{\backslash nu\}\; \backslash frac\{d\backslash nu\}\{d\backslash lambda\}\; =\; \backslash sum\_\{(i)\}\; \backslash frac\{1\}\{2\}\; \backslash frac\{\backslash hat\{p\}^i\; \backslash hat\{p\}^j\}\{1\; +\; \backslash hat\{\backslash Omega\}\_\{(i)\}\; \backslash cdot\; \backslash hat\{p\}\}\; h\_\{ij\}^\{(i)\}(t\; -\; \backslash hat\{\backslash Omega\}\_\{(i)\}\; \backslash cdot\; \backslash vec\{x\}),\; \backslash tag\{B9\}$$

which can be integrated to give the redshift

$$\begin{array}{cccc}& z\equiv \frac{{\nu }_0-\nu (t)}{{\nu }_0}=\sum _{(i)}\frac{1}{2}\frac{{\widehat{p}}^i{\widehat{p}}^j}{1+{\widehat{\mathrm{\Omega }}}_{(i)}\cdot \widehat{p}}\mathrm{\Delta }{h}_{ij}^{(i)}(t-{\widehat{\mathrm{\Omega }}}_{(i)}\cdot \vec{x}),& & \text{(B10)}\end{array}z\; \backslash equiv\; \backslash frac\{\backslash nu\_0\; -\; \backslash nu(t)\}\{\backslash nu\_0\}\; =\; \backslash sum\_\{(i)\}\; \backslash frac\{1\}\{2\}\; \backslash frac\{\backslash hat\{p\}^i\; \backslash hat\{p\}^j\}\{1\; +\; \backslash hat\{\backslash Omega\}\_\{(i)\}\; \backslash cdot\; \backslash hat\{p\}\}\; \backslash Delta\; h\_\{ij\}^\{(i)\}(t\; -\; \backslash hat\{\backslash Omega\}\_\{(i)\}\; \backslash cdot\; \backslash vec\{x\}),\; \backslash tag\{B10\}$$

which is the discrete version of Eq. (18).

APPENDIX C: THE OVERLAP REDUCTION FUNCTION IN THE HIGH-FREQUENCY LIMIT

1. Derivation of the Hellings-Downs curve

In this section we derive the Hellings and Downs curve given by Eq. (39) [49]. We begin with the definition of

the overlap reduction function, Eq. (38), and we ignore the exponential factors. Thus we wish to evaluate

$$\begin{array}{cccc}& {\mathrm{\Gamma }}_0=\beta \sum _{A=+,\times }{\int }_{S^2}d\widehat{\mathrm{\Omega }}{F}_1^A(\widehat{\mathrm{\Omega }})F_2^A(\widehat{\mathrm{\Omega }})& & \text{(C1)}\end{array}\backslash Gamma\_0\; =\; \backslash beta\; \backslash sum\_\{A=+,\; \backslash times\}\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; F\_1^A(\backslash hat\{\backslash Omega\})\; F\_2^A(\backslash hat\{\backslash Omega\})\; \backslash tag\{C1\}$$

and using the definition of $F^A(\hat{\Omega})$ given by Eq. (17) we find

$$\begin{array}{cccc}& {\mathrm{\Gamma }}_0=\frac{1}{4}\beta \sum _{A=+,\times }{\int }_{S^2}d\widehat{\mathrm{\Omega }}\frac{{\widehat{p}}_1^i{\widehat{p}}_1^j}{1+\widehat{\mathrm{\Omega }}\cdot {\widehat{p}}_1}\frac{{\widehat{p}}_2^k{\widehat{p}}_2^l}{1+\widehat{\mathrm{\Omega }}\cdot {\widehat{p}}_2}{e}_{ij}^A(\widehat{\mathrm{\Omega }})e_{kl}^A(\widehat{\mathrm{\Omega }})\mathrm{.}& & \text{(C2)}\end{array}\backslash Gamma\_0\; =\; \backslash frac\{1\}\{4\}\; \backslash beta\; \backslash sum\_\{A=+,\; \backslash times\}\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; \backslash frac\{\backslash hat\{p\}\_1^i\; \backslash hat\{p\}\_1^j\}\{1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\_1\}\; \backslash frac\{\backslash hat\{p\}\_2^k\; \backslash hat\{p\}\_2^l\}\{1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\_2\}\; e\_\{ij\}^A(\backslash hat\{\backslash Omega\})\; e\_\{kl\}^A(\backslash hat\{\backslash Omega\}).\; \backslash tag\{C2\}$$

The two unit vectors $\hat{p}1$ and $\hat{p}2$ are those pointing from the Earth toward the first and second pulsar respectively and the polarization tensors $e{ij}^+(\hat{\Omega})$ and $e{ij}^\times(\hat{\Omega})$ for a gravitational wave traveling in direction $\hat{\Omega}$ are given by Eqs. (2a) and (2b) respectively. To evaluate the integral we choose a coordinate system in which $\hat{p}_1$ is parallel to the $z$ -axis and $\hat{p}_2$ is in the $x$ - $z$ plane. Then

$$\begin{array}{cccc}& {\widehat{p}}_1=(0,0,1)& & \text{(C3a)}\end{array}\backslash hat\{p\}\_1\; =\; (0,\; 0,\; 1)\; \backslash tag\{C3a\}$$

$$\begin{array}{cccc}& {\widehat{p}}_2=(\sin \xi ,0,\cos \xi )& & \text{(C3b)}\end{array}\backslash hat\{p\}\_2\; =\; (\backslash sin\; \backslash xi,\; 0,\; \backslash cos\; \backslash xi)\; \backslash tag\{C3b\}$$

where $\xi$ is the angular separation between the two pulsars. Because we have chosen coordinates in which $\hat{p} \cdot \hat{m} = 0$ [cf. Eq. (3b)], the $\times$ -polarization terms vanish and our expression for $\Gamma_0$ becomes

$$\begin{array}{cccc}& {\mathrm{\Gamma }}_0=-\frac{1}{4}\beta {\int }_{S^2}d\widehat{\mathrm{\Omega }}\frac{{\sin }^2\theta ({\sin }^2\xi {\sin }^2\varphi -{\sin }^2\xi {\cos }^2\theta {\cos }^2\varphi -{\cos }^2\xi {\sin }^2\theta +2\sin \xi \cos \xi \sin \theta \cos \theta \cos \varphi )}{(1+\cos \theta )(1+\cos \xi \cos \theta +\sin \xi \sin \theta \cos \varphi )}& & \text{(C4)}\end{array}\backslash Gamma\_0\; =\; -\backslash frac\{1\}\{4\}\; \backslash beta\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; \backslash frac\{\backslash sin^2\; \backslash theta\; (\backslash sin^2\; \backslash xi\; \backslash sin^2\; \backslash phi\; -\; \backslash sin^2\; \backslash xi\; \backslash cos^2\; \backslash theta\; \backslash cos^2\; \backslash phi\; -\; \backslash cos^2\; \backslash xi\; \backslash sin^2\; \backslash theta\; +\; 2\; \backslash sin\; \backslash xi\; \backslash cos\; \backslash xi\; \backslash sin\; \backslash theta\; \backslash cos\; \backslash theta\; \backslash cos\; \backslash phi)\}\{(1\; +\; \backslash cos\; \backslash theta)(1\; +\; \backslash cos\; \backslash xi\; \backslash cos\; \backslash theta\; +\; \backslash sin\; \backslash xi\; \backslash sin\; \backslash theta\; \backslash cos\; \backslash phi)\}\; \backslash tag\{C4\}$$

Straightforward manipulation shows that this integral becomes

$$\begin{array}{cccc}& {\mathrm{\Gamma }}_0=\frac{1}{4}\beta (I+J)& & \text{(C5)}\end{array}\backslash Gamma\_0\; =\; \backslash frac\{1\}\{4\}\; \backslash beta\; (I\; +\; J)\; \backslash tag\{C5\}$$

with

and

$$J=-2{\sin }^2\xi {\int }_0^{\pi }d\theta \sin \theta (1-\cos \theta )K(C7)J\; =\; -2\; \backslash sin^2\; \backslash xi\; \backslash int\_0^\backslash pi\; d\backslash theta\; \backslash sin\; \backslash theta\; (1\; -\; \backslash cos\; \backslash theta)\; K\; \backslash quad\; (C7)$$

where we have defined

$$K\equiv {\int }_0^{2\pi }d\varphi \frac{{\sin }^2\varphi }{1+\cos \xi \cos \theta +\sin \xi \sin \theta \cos \varphi }\mathrm{.}(C8)K\; \backslash equiv\; \backslash int\_0^\{2\backslash pi\}\; d\backslash phi\; \backslash frac\{\backslash sin^2\; \backslash phi\}\{1\; +\; \backslash cos\; \backslash xi\; \backslash cos\; \backslash theta\; +\; \backslash sin\; \backslash xi\; \backslash sin\; \backslash theta\; \backslash cos\; \backslash phi\}.\; \backslash quad\; (C8)$$

$K$ may be trivially evaluated by contour integration in the complex plane. The result is

where the negative sign applies when $0 < \theta < \pi - \xi$ and the positive sign applies when $\pi - \xi < \theta < \pi$ . Hence we find that

Combining Eqs. (C5), (C6), and (C10), we obtain

$${\mathrm{\Gamma }}_{\text{dip}}=-\frac{1}{4}\beta {\int }_{S^2}d\widehat{\mathrm{\Omega }}\cos \chi \frac{{\sin }^2\theta ({\sin }^2\xi {\sin }^2\varphi -{\sin }^2\xi {\cos }^2\theta {\cos }^2\varphi -{\cos }^2\xi {\sin }^2\theta +2\sin \xi \cos \xi \sin \theta \cos \theta \cos \varphi )}{(1+\cos \theta )(1+\cos \xi \cos \theta +\sin \xi \sin \theta \cos \varphi )}\mathrm{.}(C16)\backslash Gamma\_\{\backslash text\{dip\}\}\; =\; -\backslash frac\{1\}\{4\}\; \backslash beta\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; \backslash cos\; \backslash chi\; \backslash frac\{\backslash sin^2\; \backslash theta\; (\backslash sin^2\; \backslash xi\; \backslash sin^2\; \backslash phi\; -\; \backslash sin^2\; \backslash xi\; \backslash cos^2\; \backslash theta\; \backslash cos^2\; \backslash phi\; -\; \backslash cos^2\; \backslash xi\; \backslash sin^2\; \backslash theta\; +\; 2\; \backslash sin\; \backslash xi\; \backslash cos\; \backslash xi\; \backslash sin\; \backslash theta\; \backslash cos\; \backslash theta\; \backslash cos\; \backslash phi)\}\{(1\; +\; \backslash cos\; \backslash theta)(1\; +\; \backslash cos\; \backslash xi\; \backslash cos\; \backslash theta\; +\; \backslash sin\; \backslash xi\; \backslash sin\; \backslash theta\; \backslash cos\; \backslash phi)\}.\; \backslash quad\; (C16)$$

As in the previous section, we write

$${\mathrm{\Gamma }}_{\text{dip}}=\frac{1}{4}\beta (I+J)(C17)\backslash Gamma\_\{\backslash text\{dip\}\}\; =\; \backslash frac\{1\}\{4\}\; \backslash beta\; (I\; +\; J)\; \backslash quad\; (C17)$$

where the first term is now given by

and $J$ is as in Eq. (C7), but $K$ is now given by

$$K={\int }_0^{2\pi }d\varphi \frac{\sin \xi \cos {\alpha }_1\cos \theta {\sin }^2\varphi +(\cos {\alpha }_2-\cos {\alpha }_1\cos \xi )\sin \theta (\cos \varphi -\sin \varphi ){\sin }^2\varphi }{\sin \xi (1+\cos \xi \cos \theta )+{\sin }^2\xi \sin \theta \cos \varphi }(C19)K\; =\; \backslash int\_0^\{2\backslash pi\}\; d\backslash phi\; \backslash frac\{\backslash sin\; \backslash xi\; \backslash cos\; \backslash alpha\_1\; \backslash cos\; \backslash theta\; \backslash sin^2\; \backslash phi\; +\; (\backslash cos\; \backslash alpha\_2\; -\; \backslash cos\; \backslash alpha\_1\; \backslash cos\; \backslash xi)\; \backslash sin\; \backslash theta\; (\backslash cos\; \backslash phi\; -\; \backslash sin\; \backslash phi)\; \backslash sin^2\; \backslash phi\}\{\backslash sin\; \backslash xi\; (1\; +\; \backslash cos\; \backslash xi\; \backslash cos\; \backslash theta)\; +\; \backslash sin^2\; \backslash xi\; \backslash sin\; \backslash theta\; \backslash cos\; \backslash phi\}\; \backslash quad\; (C19)$$

and may be evaluated by the same methods. The result is

$$K=\frac{2\pi }{{\sin }^2\xi }(a_{\pm }\cot \theta \csc \theta +b_{\pm }{\cot }^2\theta +c_{\pm }{\csc }^2\theta )(C20)K\; =\; \backslash frac\{2\backslash pi\}\{\backslash sin^2\; \backslash xi\}\; (a\_\backslash pm\; \backslash cot\; \backslash theta\; \backslash csc\; \backslash theta\; +\; b\_\backslash pm\; \backslash cot^2\; \backslash theta\; +\; c\_\backslash pm\; \backslash csc^2\; \backslash theta)\; \backslash quad\; (C20)$$

The expression in braces achieves a maximum value of unity when $\xi = 0$ , so the correct normalization constant is $\beta = 3/4\pi$ . With this normalization we recover Eq. (39).

2. Generalization to a Dipole Stochastic Background

We now generalize the Hellings-Downs curve to the case of a stochastic background with a dipole moment in the direction $\hat{D}$ . We will start by defining the following quantities:

$$\widehat{D}=(\sin {\alpha }_1\cos \eta ,\sin {\alpha }_1\sin \eta ,\cos {\alpha }_1)(C12)\backslash hat\{D\}\; =\; (\backslash sin\; \backslash alpha\_1\; \backslash cos\; \backslash eta,\; \backslash sin\; \backslash alpha\_1\; \backslash sin\; \backslash eta,\; \backslash cos\; \backslash alpha\_1)\; \backslash quad\; (C12)$$

$$\widehat{D}\cdot {\widehat{p}}_1\equiv \cos {\alpha }_1(C14a)\backslash hat\{D\}\; \backslash cdot\; \backslash hat\{p\}\_1\; \backslash equiv\; \backslash cos\; \backslash alpha\_1\; \backslash quad\; (C14a)$$

This derivation differs from the derivation of the Hellings-Downs curve only in that a factor $\hat{D} \cdot \hat{\Omega}$ must be included in the integral.

$${\mathrm{\Gamma }}_{\text{dip}}=\frac{1}{4}\beta \sum _{A=+,\times }{\int }_{S^2}d\widehat{\mathrm{\Omega }}\frac{{\widehat{p}}_1^i{\widehat{p}}_1^j{\widehat{p}}_2^k{\widehat{p}}_2^l{e}_{ij}^A(\widehat{\mathrm{\Omega }})e_{kl}^A(\widehat{\mathrm{\Omega }})}{(1+\widehat{\mathrm{\Omega }}\cdot {\widehat{p}}_1)(1+\widehat{\mathrm{\Omega }}\cdot {\widehat{p}}_2)}\widehat{D}\cdot \widehat{\mathrm{\Omega }}\mathrm{.}(C15)\backslash Gamma\_\{\backslash text\{dip\}\}\; =\; \backslash frac\{1\}\{4\}\; \backslash beta\; \backslash sum\_\{A=+,\; \backslash times\}\; \backslash int\_\{S^2\}\; d\backslash hat\{\backslash Omega\}\; \backslash frac\{\backslash hat\{p\}\_1^i\; \backslash hat\{p\}\_1^j\; \backslash hat\{p\}\_2^k\; \backslash hat\{p\}\_2^l\; e\_\{ij\}^A(\backslash hat\{\backslash Omega\})\; e\_\{kl\}^A(\backslash hat\{\backslash Omega\})\}\{(1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\_1)(1\; +\; \backslash hat\{\backslash Omega\}\; \backslash cdot\; \backslash hat\{p\}\_2)\}\; \backslash hat\{D\}\; \backslash cdot\; \backslash hat\{\backslash Omega\}.\; \backslash quad\; (C15)$$

This integral can be written aswhere the following constant terms have been defined:

$$a_{\pm }=\cos {\alpha }_1(1\mp \cos \xi )\pm \frac{\cos {\alpha }_2-\cos {\alpha }_1\cos \xi }{{\sin }^2\xi }(1\mp \cos \xi {)}^2(\text{C21a})a\_\{\backslash pm\}\; =\; \backslash cos\; \backslash alpha\_1\; (1\; \backslash mp\; \backslash cos\; \backslash xi)\; \backslash pm\; \backslash frac\{\backslash cos\; \backslash alpha\_2\; -\; \backslash cos\; \backslash alpha\_1\; \backslash cos\; \backslash xi\}\{\backslash sin^2\; \backslash xi\}\; (1\; \backslash mp\; \backslash cos\; \backslash xi)^2\; \backslash quad\; (\backslash text\{C21a\})$$

$$b_{\pm }=\mp \cos {\alpha }_1(1\mp \cos \xi )-\frac{\cos {\alpha }_2-\cos {\alpha }_1\cos \xi }{2{\sin }^2\xi }(1\mp \cos \xi {)}^2(\text{C21b})b\_\{\backslash pm\}\; =\; \backslash mp\; \backslash cos\; \backslash alpha\_1\; (1\; \backslash mp\; \backslash cos\; \backslash xi)\; -\; \backslash frac\{\backslash cos\; \backslash alpha\_2\; -\; \backslash cos\; \backslash alpha\_1\; \backslash cos\; \backslash xi\}\{2\; \backslash sin^2\; \backslash xi\}\; (1\; \backslash mp\; \backslash cos\; \backslash xi)^2\; \backslash quad\; (\backslash text\{C21b\})$$

$$c_{\pm }=-\frac{\cos {\alpha }_2-\cos {\alpha }_1\cos \xi }{2{\sin }^2\xi }(1\mp \cos \xi {)}^2(\text{C21c})c\_\{\backslash pm\}\; =\; -\backslash frac\{\backslash cos\; \backslash alpha\_2\; -\; \backslash cos\; \backslash alpha\_1\; \backslash cos\; \backslash xi\}\{2\; \backslash sin^2\; \backslash xi\}\; (1\; \backslash mp\; \backslash cos\; \backslash xi)^2\; \backslash quad\; (\backslash text\{C21c\})$$

and $a_+$ , $b_+$ , and $c_+$ are to be used in the case where the inequality $0 < \theta < \pi - \xi$ holds, and $a_-$ , $b_-$ , and $c_-$ are to be used otherwise. Thus, the integral $J$ must again be split into two sections, and the result of the integration is

$$J=4\pi (\cos \xi -1)(\cos {\alpha }_1+\cos {\alpha }_2)-16\pi (\cos {\alpha }_1+\cos {\alpha }_2){\tan }^2\frac{\xi }{2}\ln (\sin \frac{\xi }{2})(\text{C22})J\; =\; 4\backslash pi\; (\backslash cos\; \backslash xi\; -\; 1)\; (\backslash cos\; \backslash alpha\_1\; +\; \backslash cos\; \backslash alpha\_2)\; -\; 16\backslash pi\; (\backslash cos\; \backslash alpha\_1\; +\; \backslash cos\; \backslash alpha\_2)\; \backslash tan^2\; \backslash frac\{\backslash xi\}\{2\}\; \backslash ln\; \backslash left(\; \backslash sin\; \backslash frac\{\backslash xi\}\{2\}\; \backslash right)\; \backslash quad\; (\backslash text\{C22\})$$

Thus, we see that

$${\mathrm{\Gamma }}_{\text{dip}}=\pi \beta (\cos {\alpha }_1+\cos {\alpha }_2)(\cos \xi -\frac{4}{3}-4{\tan }^2\frac{\xi }{2}\ln (\sin \frac{\xi }{2}))(\text{C23})\backslash Gamma\_\{\backslash text\{dip\}\}\; =\; \backslash pi\backslash beta\; (\backslash cos\; \backslash alpha\_1\; +\; \backslash cos\; \backslash alpha\_2)\; \backslash left(\; \backslash cos\; \backslash xi\; -\; \backslash frac\{4\}\{3\}\; -\; 4\; \backslash tan^2\; \backslash frac\{\backslash xi\}\{2\}\; \backslash ln\; \backslash left(\; \backslash sin\; \backslash frac\{\backslash xi\}\{2\}\; \backslash right)\; \backslash right)\; \backslash quad\; (\backslash text\{C23\})$$

Because we wish for $\Gamma_{\text{dip}}$ to have maximal value of unity at $\xi = 0$ and $\xi = \pi$ (where it is clear that $\alpha_1 = \alpha_2 = 0$ ), we must select a normalization constant of $\beta = -3/2\pi$ .

---

[1] R.A Hulse and J.H. Taylor, Astrophys. J. 195 (1975) L51.

[2] F.B. Estabrook and H.D. Wahlquist, General Relativity and Gravitation, 6 (1975) 439.

[3] S. Detweiler, Astrophys. J. 234 (1979) 1100.

[4] M. V. Sazhin, Sov. Astron. 22 36, 1978.

[5] D. Lorimer and M. Kramer, Handbook of Pulsar Astronomy (Cambridge University Press, 2005).

[6] G. Hobbs, D. R. Lorimer, A. G. Lyne, and M. Kramer, A statistical study of 233 pulsar proper motions. Mon. Not. Roy. Astron. Soc. 360 963, 2005.

[7] F. A. Jenet, G. B. Hobbs, K. J. Lee, and R. N. Manchester, Astrophys. J. 625 L123, 2005.

[8] F. A. Jenet et al., Astrophys. J. 653 1571, 2006.

[9] J. S. B. Wyithe and A. Loeb, Astrophys. J. 590 691, 2003.

[10] A. H. Jaffe and D. C. Backer, Astrophys. J. 583 616, 2003.

[11] M. Enoki, K. T. Inoue, M. Nagashima, and N. Sugiyama, Astrophys. J. 615 19, 2004.

[12] M. Kramer, D. C. Backer, J. M. Cordes, T. J. W. Lazio, B. W. Stappers, and S. Johnston, New Astron. Rev. 48 993, 2004.

[13] A. Sesana, A. Vecchio, and C.N. Colacino, arXiv:0804.4476v2 [astro-ph], 2008.

[14] D. M. Eardley, D. L. Lee, A. P. Lightman, R. V. Wagoner, and C. M. Will, Phys. Rev. Lett. 30(18) 884, Apr 1973.

[15] D. M. Eardley, D. L. Lee, and A. P. Lightman, Phys. Rev. D 8(10) 3308, Nov 1973.

[16] L. A. Boyle and A. Buonanno, Phys. Rev. D 78 (2008) 043531. arXiv:0708.2279.

[17] T. W. B. Kibble, J. Phys. A9 1387, 1976.

[18] M. B. Hindmarsh and T. W. B. Kibble, Rept. Prog. Phys. 58 477, 1995.

[19] A. Vilenkin and E. Shellard, Cosmic strings and other Topological Defects (Cambridge University Press, 2000).

[20] R. Jeannerot, J. Rocher, and M. Sakellariadou, Phys. Rev. D68 103514, 2003.

[21] N. T. Jones, H. Stoica, and S. H. H. Tye, JHEP 07 051, 2002.

[22] S. Sarangi and S. H. H. Tye, Phys. Lett. B536 185, 2002.

[23] G. Dvali and A. Vilenkin, JCAP 0403 010, 2004.

[24] N. T. Jones, H. Stoica, and S. H. H. Tye, Phys. Lett. B563 6, 2003.

[25] E. J. Copeland, R. C. Myers, and J. Polchinski, JHEP 06 013, 2004.

[26] M. G. Jackson, N. T. Jones, and J. Polchinski, JHEP 10 013, 2005.

[27] V. Berezinsky, B. Hnatyk, and A. Vilenkin, 2000, astro-ph/0001213.

[28] T. Damour and A. Vilenkin, Phys. Rev. Lett. 85 3761, 2000.

[29] T. Damour and A. Vilenkin, Phys. Rev. D 64 064008, 2001.

[30] T. Damour and A. Vilenkin, Phys. Rev. D 71 063510, 2005.

[31] X. Siemens et al., Phys. Rev. D 73 105001, 2006.

[32] X. Siemens, V. Mandic, and J. Creighton, Phys. Rev. Lett. 98 111101, 2007.

[33] A. N. Lommen, PhD Thesis, UC Berkeley, 2001.

[34] V. M. Kaspi, J. H. Taylor and M. F. Ryba, Astrophys. J. 428, 713 (1994).

[35] R. N. Manchester, AIP Conf. Proc. 983, 584 (2008) [arXiv:0710.5026 [astro-ph]].

[36] http://www.atnf.csiro.au/research/pulsar/ppta/index.php?n=Main.PPTA.- [37] Mon. N. R. Astron. Soc., 227 (1987) 933.

• [38] N. Christensen, Phys. Rev. D 46, 5250 (1992).
• [39] E. E. Flanagan, Phys. Rev. D 48, 2389 (1993) [arXiv:astro-ph/9305029].
• [40] B. Allen and J.D. Romano, Phys. Rev. D 59 (1999) 102001.
• [41] B. Abbott et al. [LIGO Scientific Collaboration], Phys. Rev. D 69, 122004 (2004) [arXiv:gr-qc/0312088].
• [42] B. Abbott et al. [LIGO Collaboration], Phys. Rev. Lett. 95, 221101 (2005) [arXiv:astro-ph/0507254].
• [43] B. Abbott et al. [LIGO Collaboration], Astrophys. J. 659, 918 (2007) [arXiv:astro-ph/0608606].
• [44] B. Abbott et al. [LIGO Scientific Collaboration], Phys. Rev. D 76, 082003 (2007) [arXiv:astro-ph/0703234].
• [45] B. Abbott et al. [ALLEGRO Collaboration and LIGO Scientific Collaboration], Phys. Rev. D 76, 022001 (2007) [arXiv:gr-qc/0703068].
• [46] A. Lommen, J. Romano and G. Woan, in preparation.
• [47] B. Allen, J. D. E. Creighton, E. E. Flanagan and J. D. Romano, Phys. Rev. D 67, 122002 (2003) [arXiv:gr-qc/0205015].
• [48] W.J. Kaufmann, Nature 227 (1970) 157.
• [49] R.W. Hellings and G.S. Downs, Astrophys. J. 265 (1983) L39.
• [50] O. Malaspinas and R. Sturani, Class. Quantum Grav. 23 (2006) 319.
• [51] G.J. Feldman and R.D Cousins, Phys. Rev. D 57, 3873 (1998) 3873.
• [52] T.J. Loredo in Maximum Entropy and Bayesian Methods, ed. P.F. Fougere (Dordrecht: Kluwer), 81.
• [53] P.C. Gregory and T.J. Loredo, Astrophys. J. 398 (1992) 146.
• [54] S. Mitra, S. Dhurandhar, T. Souradeep, A. Lazzarini, V. Mandic, S. Bose and S. Ballmer, Phys. Rev. D 77 (2008) 042002.
• [55] N. R. Lomb, Astrophys. Space Sci. 39, 447 (1976).
• [56] J. D. Scargle, Astrophys. J. 263, 835 (1982).
• [57] J. H. Horne and S. L. Baliunas, Astrophys. J. 302, 757 (1986).
• [58] M. G. House and P. D. Mountcastle, Proceedings of the IEEE Radar Conference, 63-67 (2002).
• [59] R.W. Hellings, Phys. Rev. D 23 (1981) 832.