The stochastic background: scaling laws and time to detection for pulsar timing arrays
We derive scaling laws for the signal-to-noise ratio of the optimal cross-correlation statistic, and show that the large power-law increase of the signal-to-noise ratio as a function of the the observation time that is usually assumed holds only at early times. After enough time has elapsed, pulsar timing arrays enter a new regime where the signal to noise only scales as . In addition, in this regime the quality of the pulsar timing data and the cadence become relatively un-important. This occurs because the lowest frequencies of the pulsar timing residuals become gravitational-wave dominated. Pulsar timing arrays enter this regime more quickly than one might naively suspect. For yr observations and typical stochastic background amplitudes, pulsars with residual RMSs of less than about s are already in that regime. The best strategy to increase the detectability of the background in this regime is to increase the number of pulsars in the array. We also perform realistic simulations of the NANOGrav pulsar timing array, which through an aggressive pulsar survey campaign adds new millisecond pulsars regularly to its array, and show that a detection is possible within a decade, and could occur as early as 2016.
Pulsar timing arrays (PTAs) are sensitive to gravitational waves (GWs) in the nano-Hz band. The most promising sources at those frequencies are super-massive binary black holes (SMBBHs) that coalesce when galaxies merge. The superposition of GWs from the SMBBH mergers that have taken place throughout the history of the universe forms a stochastic background of gravitational waves [1, 2, 3, 4, 5, 6, 7, 8]. These systems may also be observable through individual periodic signals [9, 10, 11, 12, 13] and bursts [14, 15]. A number of data analysis techniques have been implemented to search pulsar timing data for the stochastic background [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], periodic signals [31, 32, 33, 34, 35, 36, 37, 38], and bursts .
Although the stochastic background produces random changes in the times-of-arrival (TOAs) of pulses from individual pulsars, its effects on the cross-correlation of timing residuals of two pulsars only depends on the angular separation between them . In this paper we study the scaling properties of the optimal cross-correlation statistic introduced in . This cross-correlation is optimized for the spectrum of the GW background, and accounts for the intrinsic noise of the pulsars. Using this statistic we also calculate the detection probability for the NANOGrav pulsar timing array, finding that a detection could occur as early as 2016, and is likely before the end of the decade.
This paper is organized as follows. In Section 2 we review the cross-correlation statistic and its signal-to-noise ratio (SNR). In Section 3 we derive analytic scaling laws for the SNR of a simplified PTA as a function of the various properties and characteristics of the array. In Section 4 we perform realistic numerical simulations, including red spin noise, of the NANOGrav pulsar timing array. We conclude in Section 5.
2 The optimal cross-correlation statistic
In  we showed that for a timing array consisting of pulsars, a simple linear (though non-invertible) transformation allows us to write the likelihood function for the residuals as a multivariate Gaussian
is a vector of the timing residuals time-series, , for all pulsars, is a vector of noise and signal parameters,
is the covariance matrix of the residuals, is the covariance matrix for the pre-fit underlying Gaussian process which contains the gravitational waves along with other sources of noise, and
is the timing model fitting matrix . The covariance matrix for the timing residuals is the block matrix
are the auto-covariance and cross-covariance matrices for each set of residuals. These matrices can be constructed by taking inverse Fourier transforms of the frequency domain auto- and cross-power spectra, and acting on them with the -matrices,
where are matrix component indices, is the power spectrum of the th pulsar, is the GW power spectrum, and are the Hellings-Downs coefficients. This treatment is somewhat different from previous Bayesian analyses of the likelihood [22, 14, 24]. Our treatment amounts to a conditional pdf whereas [22, 14, 24] used a marginalized pdf. See  for more details.
In , we showed that in the weak-signal limit the likelihood ratio maximized over amplitudes leads to the optimal cross-correlation statistic for a PTA. The optimal statistic is
where is the sum over all pulsar pairs. The normalization factor is chosen so that is also the maximum likelihood estimator for the amplitude of the stochastic background , and the amplitude independent cross-correlation matrix is defined by
In the absence of a cross-correlated signal (or if the signal is weak) the expectation value of vanishes and its standard deviation is 
so if in a particular realization we measure a value of the optimal statistic, the SNR for that realization is
This definition of the SNR measures our degree of surprise (in standard deviations away from 0), of finding cross-correlations in our data. The average SNR is
It is worth pointing out that this SNR only involves the cross-correlation terms of the likelihood, and does not include the auto-correlation terms. It is therefore quite a stringent test of the presence of the stochastic background compared to, say, a likelihood based analysis that contains contributions to the likelihood from both auto- and cross-correlations [42, 28, 29, 30]. There is a very good reason for doing this: red noise in the auto-correlations could arise from intrinsic spin noise, frequency noise, or from propagation in the interstellar medium. The only way to be certain that we have detected a stochastic background of gravitational waves is to ensure the data are cross-correlated as we expect. From a Bayesian perspective we could achieve the same thing by comparing the evidence of a model with GW red noise (with correlated and uncorrelated components), with that of a model with only uncorrelated red noise. We can think of the SNR of the optimal cross-correlation statistic as a proxy for such an analysis.
3 Scaling laws for the cross-correlation statistic
Equation 14 for the SNR can be evaluated in the time domain for any PTA configuration (cadence, intrinsic red and white noises, etc.), along with a GW background amplitude and slope. The auto-covariance and cross-covariance matrices can be constructed by taking inverse Fourier transforms of the power spectra and cross-spectra, respectively, and accounting for fitting of the timing model in the time-domain as shown in Eqs. 8 and 9.
To understand how the SNR scales with the various properties of the PTA, however, it is easier to work in the frequency domain. The timing model, involves, among other things, subtracting out a quadratic (by fitting the period and period derivative out of the TOAs). In the frequency domain this can be crudely approximated by a low frequency cutoff at , where is the time-span of the data. We can express the trace in the sum of Eq. 14 in the frequency domain as follows 
where are the Hellings Downs coefficients and the gravitational wave spectrum is
where is some reference frequency which we will take here to be 1 yr. This means we can write the average SNR for the PTA as
If we assume there is only white noise in our pulsars (in addition to GWs), and that the data are evenly sampled with cadence , namely , the integral in Eq. 17 can be written in terms of hypergeometric functions as follows
where , and is the hypergeometric function. We will investigate the use of these solutions in future work.
For now let’s consider the more simple situation where all the pulsars in our timing array have the same intrinsic white noise with RMS , and that the data are evenly sampled with cadence , i.e. we let . This means the average SNR can be written as the product
It follows from this expression that the average SNR scales like the number of pulsars in the timing array, since . The integral we need to evaluate in Eq. 19 is
which we show below can be expressed in terms of hypergeometric functions.
Before we do this, it is useful to consider two regimes in which Eq. 20 can be easily evaluated. The first is the weak-signal limit where the intrinsic white noise dominates for all frequencies . In this case Eq. 20 becomes,
because , so that the SNR becomes
The opposite regime, the strong-signal limit, occurs when the GW power is larger than the white noise, , for all frequencies . In this case Eq. 20 simply becomes,
where we’ve set the high frequency cutoff to be the Nyquist frequency. So the SNR becomes
There is an intermediate regime, where the low frequencies of the background are above the white noise level, but the highest frequencies are below it. In this regime and . The integral involved in the SNR calculation, Eq. 20, can be written as
and each of these integrals evaluated via
The argument of the hypergeometric function for the second integral of the RHS of Eq. 25 is small, and the hypergeometric function can be approximated to be unity. For the first integral, the situation is somewhat more complicated, the argument is large and we need to make an asymptotic expansion of the hypergeometric function. We give details of this calculation elsewhere , but the result is
Figure 1 shows a plot of the gravitational wave power spectrum for all three regimes as well as the white noise. The transition from the weak signal regime to the intermediate regime occurs when the lowest frequency bin of the background becomes larger than the level of the white noise, i.e. when , the condition on the RMS is then
For typical pulsar timing experiment durations of yr, and cadence of about yr, for a background with amplitude and a spectral index like the one we expect for the SMBBH background (), the pulsar timing array is in the weak signal limit only if the pulsars have white noise RMSs greater than about ns. Note that of the 17 pulsars analyzed in the 5 yr data set analyzed in , 9 have residuals with RMSs smaller than 215 ns.
Pulsar timing arrays currently have at least a few pulsars with RMS residuals considerably smaller than that, and Eq. 22 cannot be used. Table 1 shows the white noise RMS thresholds between the weak and intermediate regimes for low-, mid-, and upper ranges of the amplitude of the background and observation times of 5, 10, 15, and 20 years.
It is worth pointing out that even though the lowest frequencies of the GW background may be larger than the white noise, the red noise induced residuals need not be large. The red noise induced residuals are
for , so
For , , and yr, the red noise induced RMS is ns, which is considerably smaller than the 215 ns RMS of white noise below which we enter the intermediate regime (see Table 1). Note that the pulsar with the smallest residuals analyzed in , J1713+0747, had a timing residual RMS of 30 ns.
|5 yr||10 yr||15 yr||20 yr|
In practice the strong-signal regime is irrelevant. We enter this regime when the highest frequency bin of the background is larger than the white noise, i.e. that . In terms of the RMS this condition is
Using a cadence of yr, for a background with amplitude and a spectral index like the one we expect for the SMBBH background, , the condition on the RMS is s. The constraint can be loosened only by reducing the cadence, which reduces the highest frequencies probed by the experiment.
To summarize, the SNR scales with PTA characteristics in two different ways depending on whether the lowest frequencies of the stochastic background power spectrum are above or below the white noise level. When for all frequencies the SNR is given by
If the PTA begins in this regime, after enough time passes, the condition is no longer satisfied and the lowest frequency bins of the background become larger than the white noise. At this stage the SNR scaling changes to
Figure 2 shows the average SNR versus time for a GW stochastic background calculated in the time domain with Eq. 14, with , , , along with the the two scaling laws Eqs. 33 and 34 at late times. For the timing model we have used quadratic subtraction (fitting for an overall phase, period, and period derivative). In our frequency domain derivation we approximated the fitting of the timing model with a low frequency cutoff of . In this plot we have taken for both Eq. 33 and 34 which reflects the fact that working in the frequency domain with a low frequency cutoff at does not capture the effects of the timing model completely, though it gives the correct dependence on the PTA parameters.
The difference in the dependence on the observation time, cadence, white noise RMS of the pulsars, and amplitude of the background of the two regimes is striking, and critical for PTA optimization. The combination appears in both cases but with two different powers. While increasing the cadence and improving the white noise RMS helps greatly in the weak signal limit, their impact on the SNR in the intermediate regime is not as significant. The reason for this is that the SNR is dominated by the lowest frequency bins. When these bins become gravitational-wave dominated, the dominant contribution to the noise is the uncorrelated pulsar-term GW red noise, rather than the white noise, and this changes the scaling dramatically: the un-correlated part of the gravitational-wave stochastic background is interfering with our ability to detect the correlated piece via cross-correlations. The most effective way to beat down the uncorrelated part of the GW signal is to add more pulsars to the PTA.
PTA data sets are now reaching 20 year time spans with pulsar RMSs at the level of a few s, so the regime where Eq. 34 is valid is probably already relevant or will soon be. This is critical when considering the kinds of decisions the pulsar timing community will have to make in terms of observing strategies to ensure a prompt and confident detection of the stochastic GW background.
In Fig. 3 we illustrate this point. We plot the SNR for a number of different PTA configurations and a SMBBH stochastic background with amplitude and slope . The canonical PTA configuration, with pulsars with ns, and a cadence of yr is shown in gray. The dashed curve shows the same PTA with a cadence of yr, and though at early times the SNR for the gray higher-cadence curve is larger than the dashed curve, at late times the cadence becomes un-important as can be gleaned from Eqs. 33 and 34. The solid curve shows a PTA with pulsars timed with a cadence of 5 points per year and an RMS of ns. Both solid curves involve the same total observation time (). More pulsars and smaller cadence lead to an SNR that is significantly greater and a more confident detection in the long term. The dash-dot curve shows a PTA with pulsars with s residuals and a cadence of 5 point per year. Though initially this PTA has a smaller SNR than the pulsar PTAs, at late times it leads to a larger SNR because the RMS and cadence become unimportant compared to the number of pulsars.
The work we have described in this section began with our attempts to understand (through simulations) the time it would take a realistic PTA to detect the gravitational-wave background. We expected that the weak-signal approximation would be valid for all times, and in particular that we would benefit from the dramatic increase in the SNR with time as predicted by Eq. 33. What we found instead was that after some time the SNR grows only as the square root of the observing time. The reason for that, we now understand, is that the PTA enters the intermediate regime.
4 Time to detection for realistic pulsar timing arrays
In this section we discuss the detection potential of the NANOGrav pulsar timing array . NANOGrav started collecting data on 17 millisecond pulsars in 2005, and upper limits on the GW background using these pulsars timed up to 2010 were published last year . Through various surveys, NANOGrav has since increased the number of pulsars timed to 35, at about a rate of 3 new pulsars added every year. About half of the pulsars are timed at Green Bank, and the other half at Arecibo.
In our simulations we have used the epoch-averaged combined residuals for the 17 pulsars in Table 2 of . New high-bandwidth pulsar timing back-ends have been installed at the Green Bank Telescope (in 2010), and Arecibo (in 2012). The new back ends have a bandwith that is about a factor of 10 larger than the previous ones, and we might expect an increase in TOA precision of about a factor of 3. In fact, on average only about a factor of 2 improvement has been observed, presumably due to jitter noise . In our simulations we have assumed a factor of 2 improvement from the new hardware starting in 2010 for pulsars timed at Green Bank and 2012 for pulsars timed at Arecibo. For new pulsars added to the array we have assumed an epoch combined RMS of 200 ns prior to hardware improvements (this is the median of the 17 pulsars in ). We have taken conservative lower (), middle (), and upper () values for the amplitude of the stochastic background generated by supermassive binary black holes . At each time we perform 1000 injections of simulated signals. The detection probability is the fraction of injections that exceed an SNR of 3. The variance of the SNR is quite large in the presence of a signal , so looking at the SNR threshold alone is not sufficient, and the 90% detection probability is typically reached when the average SNR is significantly larger than 3.
Although most NANOGrav pulsars do not show strong evidence of red noise, in our simulations we have included the effects of red spin noise (with spectral index of ) . We resort to simulations for this because we have not found a means to include the effects of this type of noise into our analytic estimates of the SNR.
Figure 4 shows the results of our simulations. The three panels show red spin noise that induces an RMS of 0 ns, 5 ns, and 10 ns at 5 yrs. The figure shows that for all but the most pessimistic scenario (lowest amplitude and largest intrinsic spin noise), a detection will occur by 2023, and could occur as early as 2016.
This prompt detection is possible because of the addition of new pulsars to the timing array. The importance of adding new pulsars is illustrated in Fig. 5. The top three panels show the detection probability as a function of time for the 17 pulsars in , but adding no additional pulsars. The bottom three panels show a pulsar timing array with just 6 pulsars, with white noise RMSs of 10 ns. Clearly, increasing the number of pulsars is critical to our ability to confidently detect the stochastic background.
We have produced analytic expressions for the SNR of the optimal cross-correlation statistic for stochastic background searches as a function of the various PTA properties. We find that the weak signal limit scalings where do not apply at all times. Once enough time has elapsed the lowest frequencies of pulsar timing residuals become GW-dominated and the scaling changes to . In addition, the dependence on the cadence and the RMS of the residuals weakens significantly, and the best strategy to confidently detect the stochastic background is to increase the number of pulsars in the PTA.
PTAs may enter this new regime more quickly than one might have thought naively. Of the 17 pulsars in the 5 year data set analyzed in , 9 are already in this regime for typical values of the amplitude of the stochastic background (). Some current PTA data sets span close to 20 yrs, and for the same amplitude of the background their lowest frequency bins are GW-dominated if the RMSs of the timing residuals are smaller than a few s.
We have also made realistic simulations of the NANOGrav PTA. NANOGrav is currently timing 35 pulsars and adding more to its array at a rate of about 3 per year. Where possible, we have used the measured RMSs of millisecond pulsars, and for future discoveries we have used the median RMS of the existing set. Most NANOGrav pulsars do not show strong evidence for intrinsic red noise. Nevertheless, to be conservative, we have included the effects of red spin noise in our simulations. We find that a confident detection of the background could occur as early as 2016.
It is worth pointing out that although for stochastic background searches the RMS of timing residuals becomes un-important at late times, there are very good reasons to improve on the RMS of pulsars. For example, for other types of GW signals such as continuous waves and bursts, the sensitivity is dominated by the pulsars with the smallest residuals.
- Lommen A N and Backer D C 2001 ApJ 562 297–302
- Jaffe A H and Backer D C 2003 ApJ 583 616–631
- Wyithe J S B and Loeb A 2003 ApJ 590 691–706
- Volonteri M, Haardt F and Madau P 2003 ApJ 582 559–573
- Enoki M, Inoue K T, Nagashima M and Sugiyama N 2004 ApJ 615 19–28
- Sesana A, Vecchio A and Colacino C N 2008 MNRAS 390 192–209
- Sesana A 2012 arXiv:1211.5375
- McWilliams S T, Ostriker J P and Pretorius F 2012
- Sesana A, Vecchio A and Volonteri M 2009 MNRAS 394 2255–2265
- Sesana A and Vecchio A 2010 Phys. Rev. D 81 104008–+
- Roedig C and Sesana A 2011 arXiv:1111.3742
- Ravi V, Wyithe J S B, Hobbs G, Shannon R M, Manchester R N, Yardley D R B and Keith M J 2012
- Mingarelli C, Grover K, Sidery T, Smith R and Vecchio A 2012 Phys.Rev.Lett. 109 081104
- van Haasteren R and Levin Y 2010 MNRAS 401 2372–2378
- Cordes J and Jenet F 2012 Astrophys.J. 752 54
- Detweiler S 1979 ApJ 234 1100–1104
- Stinebring D R, Ryba M F, Taylor J H and Romani R W 1990 Physical Review Letters 65 285–288
- Lommen A N 2002 Neutron Stars, Pulsars, and Supernova Remnants ed W Becker, H Lesch, & J Trümper pp 114–+ (Preprint arXiv:astro-ph/0208572)
- Jenet F A, Hobbs G B, Lee K J and Manchester R N 2005 ApJL 625 L123–L126
- Jenet F A, Hobbs G B, van Straten W, Manchester R N, Bailes M, Verbiest J P W, Edwards R T, Hotan A W, Sarkissian J M and Ord S M 2006 ApJ 653 1571–1576
- Anholm M, Ballmer S, Creighton J D E, Price L R and Siemens X 2009 Phys. Rev. D 79 084030–+
- van Haasteren R, Levin Y, McDonald P and Lu T 2009 MNRAS 395 1005–1014
- Yardley D R B, Coles W A, Hobbs G B, Verbiest J P W, Manchester R N, van Straten W, Jenet F A, Bailes M, Bhat N D R, Burke-Spolaor S, Champion D J, Hotan A W, Oslowski S, Reynolds J E and Sarkissian J M 2011 MNRAS 414 1777–1787
- van Haasteren R, Levin Y, Janssen G H, Lazaridis K, Kramer M, Stappers B W, Desvignes G, Purver M B, Lyne A G, Ferdman R D, Jessner A, Cognard I, Theureau G, D’Amico N, Possenti A, Burgay M, Corongiu A, Hessels J W T, Smits R and Verbiest J P W 2011 MNRAS 414 3117–3128
- Cordes J and Shannon R 2012 Astrophys.J. 750 89
- Demorest P B, Ferdman R D, Gonzalez M E, Nice D, Ransom S, Stairs I H, Arzoumanian Z, Brazier A, Burke-Spolaor S, Chamberlin S J, Cordes J M, Ellis J, Finn L S, Freire P, Giampanis S, Jenet F, Kaspi V M, Lazio J, Lommen A N, McLaughlin M, Palliyaguru N, Perrodin D, Shannon R M, Siemens X, Stinebring D, Swiggum J and Zhu W W 2012 Astrophys.J. 762 (2013) 94
- van Haasteren R 2013 MNRAS 429 55–62
- Lentati L, Alexander P, Hobson M P, Taylor S and Balan S T 2012
- Taylor S R, Gair J R and Lentati L 2012
- Ellis J, Siemens X and van Haasteren R 2013 ArXiv e-prints
- Jenet F A, Lommen A, Larson S L and Wen L 2004 ApJ 606 799–803
- Yardley D R B, Hobbs G B, Jenet F A, Verbiest J P W, Wen Z L, Manchester R N, Coles W A, van Straten W, Bailes M, Bhat N D R, Burke-Spolaor S, Champion D J, Hotan A W and Sarkissian J M 2010 MNRAS 407 669–680
- Corbin V and Cornish N J 2010 arXiv:1008.1782
- Lee K J, Wex N, Kramer M, Stappers B W, Bassa C G, Janssen G H, Karuppusamy R and Smits R 2011 MNRAS 414 3251–3264
- Ellis J A, Jenet F A and McLaughlin M A 2012 ApJ 753 96
- Babak S and Sesana A 2012 Phys. Rev. D 85 044034
- Ellis J A, Siemens X and Creighton J D E 2012 ApJ 756 175
- Petiteau A, Babak S, Sesana A and de Araújo M 2013 Phys. Rev. D 87(6) 064036 URL http://link.aps.org/doi/10.1103/PhysRevD.87.064036
- Finn L S and Lommen A N 2010 ApJ 718 1400–1415
- Hellings R W and Downs G S 1983 ApJL 265 L39–L42
- Demorest P B 2007 Measuring the gravitational wave background using precision pulsar timing Ph.D. thesis University of California, Berkeley
- van Haasteren R, Levin Y, McDonald P and Lu T 2009 MNRAS 395 1005–1014
- S Chamberlin J Creighton P D J E J R X S 2013 In preparation
- Jenet F, Finn L S, Lazio J, Lommen A, McLaughlin M, Stairs I, Stinebring D, Verbiest J, Archibald A, Arzoumanian Z, Backer D, Cordes J, Demorest P, Ferdman R, Freire P, Gonzalez M, Kaspi V, Kondratiev V, Lorimer D, Lynch R, Nice D, Ransom S, Shannon R and Siemens X 2009 arXiv:0909.1058
- Cordes J M and Shannon R M 2012 ApJ 750 89