ARE WE THERE YET? TIME TO DETECTION OF NANOHERTZ GRAVITATIONAL WAVES BASED ON PULSAR-TIMING ARRAY LIMITS

With the recent first discovery of gravitational waves (GWs) in the kilohertz frequency band, we are now firmly in the era of GW astronomy (Abbott et al. 2016). For the last 10 years, three international collaborations have been collecting precise timing observations of the most stable millisecond pulsars, with the aim of identifying the imprint of gravitational waves (GWs) in the nanohertz frequency band (Hobbs 2013; Kramer & Champion 2013; McLaughlin 2013). The best-motivated GW source for such pulsar-timing arrays (PTAs) is the stochastic background signal from the cosmological population of gravitationally bound supermassive black hole binaries (SMBHBs; Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003) that are expected to form at the centers of galaxies after these merge (Sesana et al. 2004, 2008). If the binary inspirals are driven purely by GW radiation reaction at orbital separations corresponding to the PTA frequency band, the resulting timing-residual signal has a power-law spectrum with a characteristic exponent,

$$\begin{eqnarray}&&S(f)=\displaystyle \frac{{h}_{c}{(f)}^{2}}{12{\pi }^{2}{f}^{3}}=\displaystyle \frac{{A}_{{\rm{GW}}}^{2}}{12{\pi }^{2}}{\left(\displaystyle \frac{f}{{\mathrm{yr}}^{-1}}\right)}^{-13/3}\;{\mathrm{yr}}^{3},\end{eqnarray} \tag{ 1 }$$

where h_c(f) is the characteristic strain spectrum of the GW background (GWB). Theoretical expectations for the GW amplitude A_GW center around a few 10⁻¹⁵ (Sesana 2013; McWilliams et al. 2014; Ravi et al. 2014). Recent observational upper limits from the three international pulsar-timing consortia (Arzoumanian et al. 2015a; Lentati et al. 2015; Shannon et al. 2015) imply A_GW ≲ 10⁻¹⁵ at 95% confidence. These limits appear to be in contrast with theoretical expectations, so much so that it seems plausible that environmental effects in galactic centers either stall the formation of gravitationally bound SMBHBs or accelerate their inspiral; in both cases, the effective A_GW is reduced at the frequencies where PTAs are most sensitive (Arzoumanian et al. 2015a; Shannon et al. 2015).

Indeed, the PTA detection of the stochastic GWB from SMBHBs, considered imminent only a few years ago (Siemens et al. 2013; McWilliams et al. 2014), now seems to be receding toward the future. Is this really the case? In this Letter, we answer this question quantitatively. Namely, given the upper limit A_ul = 10⁻¹⁵ obtained by the Parkes PTA (PPTA; Shannon et al. 2015), we ask when we can expect to make a positive detection using different PTAs: the PPTA (Hobbs 2013), limited to the four low-timing-noise pulsars used for the upper limit; the North American Nanohertz Observatory for Gravitational-waves (NANOGrav; McLaughlin 2013); the European PTA (EPTA; Kramer & Champion 2013); the International PTA (IPTA; see, e.g., Hobbs et al. 2010); and the pulsar-timing project that will be supported by SKA1 (Janssen et al. 2015).

We adopt the frequentist formalism developed by Rosado et al. (2015), and we characterize the PTAs simply by listing, for each pulsar, the duration of the observation and the levels of measurement and timing noise. While GW searches and upper limits with PTAs have recently been given a Bayesian treatment (van Haasteren et al. 2009; Lentati et al. 2013; van Haasteren & Levin 2013), the frequentist approach based on optimal statistics is both convenient and appropriate for the purpose of this paper, because it dispenses with the simulation of actual data sets, and because the quantification of detection probability (DP) is intrinsically a frequentist statement (see, e.g., Vallisneri 2012). Furthermore, experience shows that frequentist upper limits and detection prospects are rather close to Bayesian results (see, e.g., Arzoumanian et al. 2015a).

We proceed as follows: for a specified PTA and true GWB amplitude (A_true), we compute the DP as a function of observation time, using the detection statistic described below. Detection is a probabilistic endeavor since the actual realization of measurement and timing noise may obscure or expose the underlying GW signal. Indeed, we do not currently have access to A_true, but only to the observed upper limit A_ul. Setting upper limits is also a probabilistic endeavor because different realizations of noise lead to different A_ul given the same A_true, so we use the upper limit optimal statistic (also described below) to compute $p({A}_{\mathrm{ul}}| {A}_{\mathrm{true}})$. By introducing an astrophysically motivated prior p(A_true), we can then use Bayes's theorem to obtain $p({A}_{\mathrm{true}}| {A}_{\mathrm{ul}})\propto p({A}_{\mathrm{ul}}| {A}_{\mathrm{true}})p({A}_{\mathrm{true}})$. Finally, we obtain the expected DP(A_ul) as $\int \mathrm{DP}({A}_{\mathrm{true}})p({A}_{\mathrm{true}}| {A}_{\mathrm{ul}})\;{{dA}}_{\mathrm{true}}$.⁵

We conclude that, while small arrays of a few low-timing-noise pulsars (represented here by the Shannon et al. 2015 configuration) provide the most constraining upper limits on the nanohertz GWB, they are suboptimal for detection in comparison to larger arrays such as NANOGrav, EPTA, full PPTA, and the IPTA, which are expanded regularly with newly discovered pulsars. Indeed, detection and upper limits have rather different demands, and the PPTA itself may employ a ∼20-pulsar array for detection in the near future (Reardon et al. 2016). Encouragingly, we find that larger PTAs have ∼80% probability of making a detection in the next 10 years, even allowing for a reduced GWB signal due to significant binary stalling or environmental influences.

The details of the frequency-domain statistical framework used here can be found in Rosado et al. (2015). In the frequentist context, a statistic X is a single number that summarizes the (noisy) data with respect to the presence of a (GW) signal, so that, on average, a higher A_true yields a higher X, with fluctuations due to the range of possible noise realizations. Given a set of data, there is one observed $\hat{X}$, and setting an upper limit ${A}_{\mathrm{ul}}^{95\%}$ amounts to stating that if A_true were equal to ${A}_{\mathrm{ul}}^{95\%}$ for 95% of noise realizations, the observed X would have been higher than $\hat{X}$.

Detection schemes based on statistics, on the other hand, proceed as follows: one considers separately the case where a signal is present in the data at a certain level, and the case where it is not, and computes the probability distribution of the statistic (over the ensemble of noise realizations) in both cases. The zero-signal distribution is used to set a detection threshold as a function of a false-alarm probability (FAP); the signal distribution is used to compute the probability that for a certain A_true the statistic will exceed the threshold—that is, the DP (or detection efficiency).

In the following, we use cross-correlation statistics for both upper limit and detection considerations, making use of the correlated influence of a GWB signal on pulsars that are widely separated on the sky.

Upper limits—For the discussion of upper limits in this work, we will use the so-called A-statistic of Rosado et al. (2015) modified so that the standard deviation is measured in units of squared GWB amplitude and that the value of the statistic is equal to the injected GWB amplitude on average, as is done in Chamberlin et al. (2015). Assuming that the distribution of this statistic is Gaussian (a reasonable approximation at low signal-to-noise ratios; Chamberlin et al. 2015), we obtain

$$\begin{eqnarray}&&p({A}^{2}| {A}_{{\rm{true}}})=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{0}^{2}}}\mathrm{exp}\left(-\displaystyle \frac{{({A}^{2}-{A}_{{\rm{true}}}^{2})}^{2}}{2{\sigma }_{0}^{2}}\right)\end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray}&&{\sigma }_{0}={A}_{{\rm{true}}}^{-2}{\left(2\displaystyle \sum _{k}\displaystyle \sum _{{ij}}{{\rm{\Gamma }}}_{{ij}}^{2}\displaystyle \frac{{S}_{h0}^{2}}{{P}_{i}{P}_{j}}\right)}^{1/2},\end{eqnarray} \tag{ 3 }$$

where ${\sum }_{k}$ denotes a sum over frequencies f_k, ${\sum }_{{ij}}$ denotes a sum over pulsar pairs, Γ_ij is the overlap reduction function (i.e., this is the Hellings & Downs 1983 curve for an isotropic GWB) for pulsars i and j, ${S}_{h0}\equiv {S}_{h0}({f}_{k})$ is the modeled cross-power spectral density, and ${P}_{i}\equiv {P}_{i}({f}_{k})$ is the intrinsic noise power spectral density. From this Gaussian probability distribution, it can be shown that

$$\begin{eqnarray}&&{A}_{{\rm{ul}}}^{2}={\hat{A}}^{2}+\sqrt{2}{\sigma }_{0}{{\rm{erfc}}}^{-1}(2(1-C)),\end{eqnarray} \tag{ 4 }$$

where $\hat{A}$ is the measured value of the GWB amplitude and C is our upper limit confidence (e.g., 95%). Thus, the distribution $p({A}_{{\rm{ul}}}^{2}| {A}_{{\rm{true}}})$ is a Gaussian distribution with standard deviation σ₀ and mean ${A}_{{\rm{true}}}^{2}+\sqrt{2}{\sigma }_{0}{{\rm{erfc}}}^{-1}(2(1-C))$. Note that to compute $p({A}_{{\rm{ul}}}| {A}_{{\rm{true}}})$ we perform a coordinate transformation to obtain $p({A}_{{\rm{ul}}}| {A}_{{\rm{true}}})$ = $2{A}_{{\rm{ul}}}p({A}_{{\rm{ul}}}^{2}| {A}_{{\rm{true}}})$.

Detection—In the previous section, we used the A-statistic for upper limits in order to compare our results with upper limits computed using the statistic presented in Chamberlin et al. (2015). For assessing our DP, however, we will make use of the B-statistic of Rosado et al. (2015). In this case, we wish to determine the DP

$$\begin{eqnarray}&&{\rm{DP}}=\displaystyle \frac{1}{2}{\rm{erfc}}\left[\displaystyle \frac{\sqrt{2}{\sigma }_{0}{{\rm{erfc}}}^{-1}(2{\alpha }_{0})-{\mu }_{1}}{\sqrt{2}{\sigma }_{1}}\right],\end{eqnarray} \tag{ 5 }$$

where α₀ is the FAP; μ₁, σ₁, and σ₀ (distinct from that used previously) are the mean in the presence of a signal and the standard deviation in the presence (subscript 1) and in the absence (subscript 0) of a signal (see Equations (A16)–(A18) of Rosado et al. 2015 for more details). Both μ₁ and σ₁ grow with increasing signal amplitude. Here, we choose an FAP of 0.13% to match the 3σ detection threshold of Siemens et al. (2013).⁶

Implementation—In the work of Rosado et al. (2015), it is argued that timing model subtraction (namely, subtraction of the quadratic term of the timing model; due to pulsar spin down) can be emulated by not including the lowest frequency in the sum over k frequencies. However, we have found that these analytic results agree much better with simulations including timing model subtraction if all frequencies are included in the sum. Furthermore, to emulate data sets with different time spans, we only include frequencies in the sum that are greater than 1/T_min, where T_min is the shortest time span for the given pulsar pair.

In this work, we always include the GWB power spectral density in the intrinsic pulsar noise P_i. This is meant to mimic a real optimal statistic analysis where the intrinsic noise would be based on single-pulsar analyses in which there is no way to distinguish the intrinsic noise from the GWB power. This will make our analysis conservative in the sense that we are overestimating the intrinsic noise.

In the following, we characterize detection prospects in terms of small and large PTAs, and specifically consider the following five configurations:

• 1.  
  PPTA4, consisting of the four pulsars described in Shannon et al. (2015) with associated measurement (white) and timing (red) noise (Reardon et al. 2016; Shannon et al. 2015).
• 2.  
  NANOGrav+, consisting of the 37 pulsars described in Table 3 of Arzoumanian et al. (2015b) with their associated noise properties, to be expanded in future observations by adding four new pulsars each year with 250 ns TOA measurement noise, in line with current expectations.
• 3.  
  EPTA+, consisting of the 42 pulsars described in Table 3 of Caballero et al. (2015) with their associated noise properties, and regularly adding pulsars as in (2). Some of these additional pulsars are monitored with the Large European Array for Pulsars (LEAP; Kramer & Champion 2013; Bassa et al. 2016), which synthesizes a 194 m steerable dish from the five EPTA telescopes.
• 4.  
  A conservative IPTA+, consisting of the 49 pulsars of the first IPTA data release, described in Verbiest et al. (2016) with measurement and timing noise⁷ properties as in Lentati et al. (2016), and expanding by six new pulsars each year, again with 250 ns TOA measurement noise.
• 5.  
  A theorist's PTA (TPTA), consisting of the toy configuration of Rosado et al. (2015), which may be supported by an advanced radio telescope such as SKA1: 50 pulsars with 100 ns TOA measurement noise and no intrinsic timing noise.

These configurations are summarized in Table 1.⁸

We consider four combined estimates for the amplitude and spectral shape of the stochastic GWB from SMBHBs, ranging from more optimistic (detection wise) to more conservative. For the amplitude, we adopt the observationally motivated lognormal distribution of Sesana (2013) with mean ${\mathrm{log}}_{10}\;{A}_{\mathrm{true}}=-15$ and standard deviation of 0.22; however, we also consider the case where 90% of binaries require longer than ∼9–13.7 Gyr before reaching the PTA band (i.e., the binaries stall; McWilliams et al. 2014; S. Burke-Spolaor & J. Simon 2015, private communication), which corresponds to a mean ${\mathrm{log}}_{10}\;{A}_{\mathrm{true}}=-15.5$, with the same standard deviation. For each of these amplitude priors, we examine the purely GW-driven power-law spectrum of Equation (1), as well as the case where the GW strain has a turnover at f = 1/(11 year) caused by SMBHB interactions with stars in galactic nuclei, which accelerate binary inspiral and remove low-frequency power with respect to the pure power law (see, e.g., Arzoumanian et al. 2015a and Sampson et al. 2015). The 1/(11 year) turnover frequency was chosen to lie just beyond the low-frequency reach of the current PPTA data set since the Shannon et al. (2015) analysis still found ∼9% consistency with the purely GW-driven model of Sesana (2013) for f ≥ 1/(11 year). Similarly, the recent NANOGrav (Arzoumanian et al. 2015a) analysis finds 20% consistency using a nine-year data set (Arzoumanian et al. 2015b). Furthermore, our choice of a turnover at 1/(11 year) corresponds to reasonable assumptions about the density of stellar populations interacting with SMBHBs in galactic nuclei (see Figures 10 and 11 of Arzoumanian et al. 2015a). At very low frequencies, the frequency scaling of the timing-residual spectrum becomes S(f) ∝ 1/f due to the environmental influences.

In the top panel of Figure 1, we show the probability distribution $p({A}_{\mathrm{true}}| {A}_{\mathrm{ul}})$ obtained by setting the 95% A-statistic upper limit to the Shannon et al. (2015) value, 1 × 10⁻¹⁵, and by assuming no-stalling (solid blue line) and 90% stalling (faded solid red line) amplitude priors, with pure power-law true spectra. The dashed lines correspond to turnover spectra.⁹ For comparison, in the bottom panel of Figure 1, we show also $p({A}_{\mathrm{true}}| {A}_{\mathrm{ul}})$ as computed for NANOGrav's optimal-statistic upper limit 1.3 × 10⁻¹⁵ (Arzoumanian et al. 2015a). The underlying astrophysically motivated priors for no-stalling and 90% stalling are shown in each panel as blue and red dash-dotted lines, respectively.

Zoom In Zoom Out Reset image size

Next, we compute the B-statistic GW DPs, Equation (5), for each PTA, as a function of A_true and of the observation time T beyond current data sets. The resulting DPs are shown in Figures 2: the left and right columns correspond to the pure power-law and turnover spectra, respectively.¹⁰

Zoom In Zoom Out Reset image size

As mentioned above, for NANOGrav+/EPTA+ and IPTA+ we add new pulsars at rates of four and six per year, respectively. This annual expansion has a positive impact on the DP since it provides more independent pulsar pairs at differing angular separations, which help discriminate GWs with quadrupolar spatial correlations (the Hellings & Downs 1983 curve, or the more general anisotropic signatures discussed by Mingarelli et al. 2013 and Gair et al. 2014) from non-spatially correlated red-noise processes (Siemens et al. 2013). The unique GW correlation signature is the smoking gun for detection, and we must coordinate future efforts to maximize its measurement.

We obtain the expected DPs ${\langle \mathrm{DP}\rangle }_{{A}_{\mathrm{true}}}$ as a function of time by integrating DP (A_true, T) against the $p({A}_{\mathrm{true}}| {A}_{\mathrm{ul}})$ curves of Figure 1 (top panel). The resulting ${\langle \mathrm{DP}\rangle }_{{A}_{\mathrm{true}}}$ are shown in Figures 3. As in Figure 1, the darker blue and lighter red curves correspond to no-stalling and 90% stalling amplitude priors. For the PPTA4 array, expected DPs remain below 10% throughout the next 20 years of observation, even in the most optimistic GW-signal scenario. This is unsurprising—a PTA consisting of a few exquisitely timed pulsars may provide very tight upper limits, but it will not usually yield convincing detection statistics, which require an array of many pulsars to map out the expected spatial correlation signature.

Zoom In Zoom Out Reset image size

By contrast, large pulsar arrays such as NANOGrav+, EPTA+, and IPTA+ provide high DPs even with strong binary environmental couplings since they allow the quadrupolar spatial correlation signature of a GWB to be mapped by many different pulsar pairings. We expect these results to be qualitatively the same for a full PPTA that regularly adds pulsars to the array. In NANOGrav+, EPTA+, and IPTA+, DPs begin to grow rapidly after only five years of observation beyond current data sets. Binary stalling (the red curves) stunts this growth by three years at most. While the low-frequency turnover reduces DP, its effect is mitigated by the large number of pulsars and the annual catalog expansion.¹¹

Finally, we see that the large number of well-timed pulsars in the TPTA builds highly convincing DPs after only a few years of operation. By the time the pulsars have been observed long enough that the influence of the turnover at f = 1/(11 year) may be noticeable, the DP is already close to unity. The same is true were we to add timing noise at currently known levels to these TPTA pulsars—the DP will already be close to unity when low-frequency timing noise begins to have a significant influence.

We stress that there are several important caveats to our analysis. We have focused on detecting the stochastic background rather than deterministic signals. Even for these, using data from more pulsars is desirable in order to build evidence for a coherent signal. We have assumed for each pulsar a single value of measurement noise that does not improve over time (as occurs in reality when receivers and backends are upgraded), nor do we consider interstellar medium effects such as dispersion (see, e.g., Stinebring 2013). The influence of a GWB spectral turnover depends on its frequency, which is a function of the typical environmental properties of galactic nuclei (either directly or in how these properties drive SMBHB orbital eccentricity)—our choice corresponds to reasonable assumptions about the stellar mass density of SMBHB environments. The timelines for the growth of detection in each PTA are approximate, intended to emphasize the differences between PTAs suited to upper limits versus detection and to demonstrate the influence of various binary stalling and environmental scenarios on DPs.

We conclude by emphasizing the different demands of placing stringent upper limits on the stochastic background versus actually detecting it.

• 1.  
  Highly constraining, astrophysically significant upper limits are achievable with only a few exquisitely timed pulsars, but such a PTA is suboptimal for the detection of a stochastic GWB.
• 2.  
  Timing many pulsars allows for the quadrupolar spatial correlation signature of the SMBHB background to be sampled at many different angular separations, enhancing prospects for detection.
• 3.  
  Adding more pulsars regularly to PTAs will continually improve DP, in addition to the gains already made by timing existing pulsars for longer, and will help to mitigate the deleterious influences of binary stalling and environmental couplings.

It is our pleasure to thank Pablo Rosado, Alberto Sesana, Jonathan Gair, Lindley Lentati, Sarah Burke-Spolaor, Xavier Siemens, Maura McLaughlin, Joseph Romano, and Michael Kramer for very useful suggestions. We also thank the full NANOGrav collaboration for their comments and remarks. S.R.T. was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities through a contract with NASA. M.V. acknowledges support from the JPL RTD program. J.A.E. and R.v.H. acknowledge support by NASA through Einstein Fellowship grants PF4-150120 and PF3-140116, respectively. C.M.F.M. was supported by a Marie Curie International Outgoing Fellowship within the European Union Seventh Framework Programme. This work was supported in part by National Science Foundation Physics Frontier Center award No. 1430284 and by grant PHYS-1066293 and the hospitality of the Aspen Center for Physics. This research was performed at the Jet Propulsion Laboratory, under contract with the National Aeronautics and Space Administration.