Does a 'stochastic' background of gravitational waves exist in the pulsar timing band?
Abstract
We investigate the effects of gravitational waves (GWs) from a simulated population of binary supermassive black holes (SMBHs) on pulsar timing array datasets. We construct a distribution describing the binary SMBH population from an existing semianalytic galaxy formation model. Using realizations of the binary SMBH population generated from this distribution, we simulate pulsar timing datasets with GWinduced variations. We find that the statistics of these variations do not correspond to an isotropic, stochastic GW background. The “Hellings & Downs” correlations between simulated datasets for different pulsars are recovered on average, though the scatter of the correlation estimates is greater than expected for an isotropic, stochastic GW background. These results are attributable to the fact that just a few GW sources dominate the GWinduced variations in every Fourier frequency bin of a 5year dataset. Current constraints on the amplitude of the GW signal from binary SMBHs will be biased. Individual binary systems are likely to be detectable in 5year pulsar timing array datasets where the noise is dominated by GWinduced variations. Searches for GWs in pulsar timing array data therefore need to account for the effects of individual sources of GWs.
Subject headings:
black hole physics — galaxies: evolution — gravitational waves — methods: data analysis1. Introduction
Detecting and performing science with gravitational waves (GWs) is currently a major goal of experimental astrophysics. Pulsar timing arrays (PTAs, Hellings & Downs, 1983; Foster & Backer, 1990) employ contemporaneous timing observations of millisecond radio pulsars in order to search for the effects of GWs perturbing the spacetime metric along each pulsarEarth line of sight. PTAs are complementary to other GW detection experiments, such as groundbased and spacebased interferometers, in that they are sensitive to GWs in a different frequency band (currently nHz, i.e., periods of five years to a few months). In this frequency band, the most promising astrophysical sources of GWs are binary supermassive black holes (SMBHs).
SMBHs, with masses in the range to ,
(1) 
where is the characteristic strain as a function of GW frequency, , at the Earth, and is the characteristic strain at a frequency of . The characteristic strain is the GW strain per logarithmic frequency interval, and is given by
(2) 
where is the onesided strain power spectral density (PSD).
The combined GW signal from binary SMBHs is widely assumed to form an isotropic, stochastic GW background. The value of is used to specify the amplitude of the background. The most recent predictions for the value of were made by Sesana et al. (2008), who found a likely range of to . This range of predictions was derived by considering a variety of SMBH seeding, accretion and feedback scenarios, as well as uncertainties in the galaxy merger rate and in the SMBH mass function.
PTA projects are based on observations of periodic pulses of radio emission from pulsars using large radio telescopes. These observations lead to measurements of pulse times of arrival (ToAs) at the observatories. ToA measurements are typically conducted every few weeks over many years. Pulsar timing involves parameters of a physical model for the ToAs being fitted to the measured ToAs, with the differences between the observed ToAs and the model predictions being the timing residuals. ToAs can be modeled using combinations of deterministic and stochastic processes.
GWs incident on the Earth (and on the pulsars) cause shifts in the measured pulse frequencies of the pulsars (Sazhin, 1978; Detweiler, 1979). For a pulsar, indexed by , with intrinsic rotation frequency , consider a GWinduced shift, , to this frequency. This shift is a function both of time, , and the Earthpulsar direction vector, . The resulting discrete timeseries of GWinduced variations to the ToAs, (the subscript indicates that is sampled at times ), is given by (Detweiler, 1979)
(3) 
For any GW signal, the expected values of zerolag crosscorrelations between the timeseries for different pulsars can be specified. For any stochastic GW signal, the expected value of the normalized correlation for each pulsar pair is expressed in terms of the angular separations between the pulsars as (Hellings & Downs, 1983; Jenet et al., 2005)
(4) 
where , is the angular separation between pulsars and , and the Kroenecker delta, , is unity if , and zero otherwise. This function is known as the Hellings & Downs curve. An isolated source of GWs will give rise to correlations that are different from the the Hellings & Downs curve. Measurements of correlations between pulsar timing datasets that are attributable to the effects of GWs are necessary for the detection of GWs with PTAs (Jenet et al., 2005; Yardley et al., 2011; Demorest et al., 2012).
The prospect of detecting or constraining the amplitude of a background of GWs from binary SMBHs has been the primary rationale for the development of the PTA concept. Various works have placed upper bounds on the value of for a background with the characteristic strain spectral form of Equation 1 (Jenet et al., 2006; van Haasteren et al., 2011; Demorest et al., 2012). The best published upper bound (van Haasteren et al., 2011) finds that with 95% confidence. All analysis methods developed to study the combined GW signal from binary SMBHs with PTAs assume that the timeseries for multiple pulsars can be described as a specific stochastic process. We describe the exact nature of this assumption in §2.
In this paper, we elucidate the statistical nature of the ToA variations induced by GWs from binary SMBHs. We accomplish this by modeling the GW signal from the predicted population of binary SMBHs, and by simulating realizations of corresponding to realizations of the GW signal. This study is critical to the validity of interpreting published upper limits on as representative of limits on the mean characteristic strain spectrum of GWs predicted to arise from binary SMBHs. Our results are also important for the optimization of GW detection techniques with PTAs. In §3, we outline our method of simulating pulsar timing datasets including GWs from the predicted population of binary SMBHs. Our analysis and results are presented in §4 and §5, and we discuss the interpretation and implications of our results in §6. We present our conclusions in §7.
Throughout this work, we assume a CDM concordance cosmology based on a combined analysis of the firstyear WMAP data release (Spergel et al., 2003) and the 2dF Galaxy Redshift Survey (Colless et al., 2001), with , , , and km s Mpc. Although these parameter values have since been superseded by more recent observations, we adopt them in order to remain consistent with the model we use for the binary SMBH population (Guo et al., 2011). A list of important symbols in this paper is shown in Table 1, along with the sections of the text in which they are introduced.
Symbol  Section  Description 

1  GWinduced ToA variation for pulsar at time .  
1  Expected zerolag normalized crosscorrelation between and timeseries.  
2  Expected PSD of timeseries.  
2  Periodogram estimator of at frequency .  
3.1  GW strain amplitude divided by frequency dependence.  
3.1  Binned distribution of binary SMBHs derived from realizations of the Millennium and MillenniumII coalescence lists.  
3.1  Average of 1000 realizations of .  
3.1  Analytic fit to .  
3.2  derived in terms of .  
3.2  Expected GW characteristic strain spectrum derived in terms of .  
4  1 ns rms ToA variation at time .  
4  Sum of and .  
4  Expected PSD of timeseries.  
4  Periodogram estimator of at frequency .  
5  Estimator of . 
2. The current model for ToA variations induced by gravitational waves from binary SMBHs
ToA variations induced by GWs from binary SMBHs () are commonly modeled among the PTA community as a widesense stationary random Gaussian process. This is based on the hypothesis that many GW sources forming a GW background contribute to the ToA variations, resulting in a statistical process governed by the central limit theorem. While the nature of the random Gaussian model for has been extensively described elsewhere (e.g., van Haasteren et al., 2009), we summarize it here for completeness.
The key property of a random Gaussian process is that a linear combination of samples will have a joint normal distribution function. Different samples need not be statistically independent. The distribution of samples from a (zeromean) random Gaussian process is characterized by the covariance matrix of the samples. Consider a vector, , containing samples of . That is,
(5) 
Let be another vector defined similarly to , corresponding to a pulsar , containing simultaneouslyobtained samples of . Under the random Gaussian assumption, the joint probability distribution of the samples in and , which we denote , is given by
(6) 
Here, is an matrix containing the covariances between the samples of and ; that is, element of is given by the covariance between and . As the Gaussian process is widesense stationary, each element of depends only on the time difference between samples and for pulsar and respectively. Elements of are sampled from a covariance function, , between the GWinduced ToA variations for pulsars and . This covariance function is defined by the inverse Fourier transform of the onesided PSD, , of GWinduced ToA variations for a given pulsar:
(7) 
Here, denotes a complex Fourier transform, and is a timelag. The PSDs of the GWinduced ToA variations for all pulsars are equivalent, and given by (Jenet et al., 2006)
(8) 
for a GW signal with the expected characteristic strain spectrum .
The above discussion applies equivalently if pulsar and pulsar are the same pulsar, or if they are different pulsars. The Hellings & Downs factor , defined in Equation 4, accounts for the correlation between GWinduced TOA variations for different pulsars. If has the form in Equation 1, we have . The GWinduced ToA variations for each pulsar will therefore be a “red” process. In this work, we only consider timeseries with finite lengths .
We are interested in comparing a new model for with the random Gaussian model described in this section. To this end, we need to be able to simulate realizations of as a random Gaussian process. Multiple PTA groups test their data analysis algorithms by simulating realizations of using the GWbkgrd plugin (Hobbs et al., 2009) to the tempo2 pulsar timing package (Hobbs et al., 2006). While this plugin does not explicitly generate random Gaussian realizations of by construction, we and others (van Haasteren et al., 2011; Demorest et al., 2012) have checked that it approximates a random Gaussian process well.
In the plugin GWbkgrd, a number of GW oscillators, , are simulated between GW frequencies and , with the normally distributed and GW polarization amplitudes set to be purely real with zero mean, variance
(9) 
and frequency probability distribution, , given by
(10) 
ToA variations calculated for a given pulsar at different times corresponding to GWs from each of these oscillators are summed to produce a realization of the timeseries. The frequency limits and are generally chosen respectively to be much less than the and much greater than the Nyquist frequency corresponding to the minimum sampling interval.
We make a distinction between the expected PSD of ToA variations induced by GWs from binary SMBHs, as defined in Equation 8, and estimates of this PSD based on realizations of the ToA variations. A commonlyused nonparametric, unbiased estimator of the PSD of a timeseries is the periodogram (Schuster, 1898). The periodogram, , of is defined as
(11) 
where DFT denotes a discrete Fourier transform. We adopt the following standard definition for the DFT of samples of :
(12) 
where in this case. The DFT is evaluated for frequencies
(13) 
where is the interval (assumed to be constant) between samples of . Throughout this work, we estimate the PSD, , of realizations of by evaluating .
3. Simulating pulsar ToA variations accounting for binary SMBH population characteristics
In this section, we describe a new method of simulating ToA variations caused by GWs from the predicted population
of binary SMBHs. Various works have presented
models for the cosmic demographics of binary SMBHs (Dotti et al., 2012). More recently, such efforts have been based on analytic
prescriptions for baryon physics applied to dark matter halo merger tree catalogues from Nbody simulations (Guo et al., 2011, hereafter G11).
In particular, the G11 prescriptions were applied to merger trees from both the Millennium (Springel et al., 2005) and the
MillenniumII (BoylanKolchin et al., 2009) simulations. The Millennium and MillenniumII simulations follow the evolution of dark matter
structures, using the same physical prescriptions and number of particles. The MillenniumII simulation was however carried out in a
comoving cubic volume with onefifth the side length as that of the Millennium simulation, with the aim of resolving smallerscale dark matter structures than the
Millennium simulation.
The G11 model is the latest in a series (Springel et al., 2005; Croton et al., 2006; De Lucia & Blaizot, 2007) of semianalytic prescriptions applied to the Millennium simulations. A host of observables of galaxies at low redshifts are reproduced, along with the redshiftevolution of the quasar population and star formation. Of most relevance here is that the model also traces the SMBH population, reproducing the SMBHgalaxy scaling relations in their slopes, normalization and scatters, as well as the inferred SMBH mass function (Marulli et al., 2008).
We base our description of the binary SMBH population emitting GWs on the prediction for the SMBHSMBH coalescence rate from the G11 model. We fit an analytic function to the distribution of binary SMBHs, and randomly draw GW sources from this distribution to produce realizations of the GW sky corresponding to binary SMBHs. We then add the effect of each GW source to simulated pulsar ToA datasets in order to analyze the GWinduced ToA variations.
This work is different from previous attempts to model the GW signal from binary SMBHs. Initial attempts (e.g., Jaffe & Backer, 2003) to predict the mean GW characteristic strain spectrum from binary SMBHs used empirical determinations of the galaxy merger rate and the SMBH mass function. Wyithe & Loeb (2003) predicted the GW spectrum by analytically following the dark matter halo merger hierarchy in an extended PressSchechter framework, and by deriving the SMBH coalescence rate by relating the SMBH masses to the halo circular velocities. Sesana et al. (2008) considered the possible range of predictions of the characteristic strain spectrum, using Monte Carlo realizations of dark matter halo merger trees and various prescriptions for SMBH growth.
The key difference between the present work and previous calculations of the GW signal from binary SMBHs is that we are chiefly concerned with the statistics of . Our approach to modeling the binary SMBH population is similar to Sesana et al. (2009) in our use of mock galaxy catalogues derived from analytic prescriptions applied to the Millennium simulations. However, whereas Sesana et al. (2009) modeled the SMBH population by using empirical SMBHgalaxy scaling relations combined with (earlier) mock catalogues, we utilize SMBHs modeled by G11 in a selfconsistent framework which reproduces the relevant observables.
3.1. Modeling the distribution of binary SMBHs
As in the previous works discussed above, we consider all binary SMBHs to be in circular orbits, and use expressions for the resulting GW emission presented by Thorne (1987). We briefly discuss the assumption of circular orbits in §6.3. The strain amplitude, , at frequency of GWs from a circular binary, averaged over all orientations and polarizations, is given by:
(14) 
where is the universal gravitational constant, is the chirp mass of the binary, is the vacuum speed of light, is the redshift, and is the comoving distance to the binary. The evolution of the received GW frequency with observed time, , is determined by
(15) 
The restframe binary orbital frequency is given by . We assume, after Hughes (2002), that the maximum received frequency, , is attained at a binary separation of three Schwarzschild radii, corresponding to the last stable orbit:
(16) 
assuming a mass ratio of unity.
The mock galaxy catalogues resulting from the G11 model are available online
(17) 
where
(18) 
and is the skyintegrated comoving volume shell between redshifts and . Also, , and the derivative was obtained from Equation 15. is the predicted distribution of binary SMBHs in (which corresponds to the frequencyindependent GW “power”, or squared strain amplitude) and the observed GW frequency.
For chirp masses below , the limited capability of the Millennium simulation to resolve lowmass halos caused an underprediction of the chirp mass function as compared to the MillenniumII simulation. In order to ensure a complete chirp mass function, we included binary SMBHs with from the Millennium list of coalescence events, and binaries with from the MillenniumII list.
Some degrees of randomization in the coalescence lists were possible. First, in cases where more than two SMBHs coalesced to form a single SMBH between redshift snapshots, the merger order was not specified. In these instances we randomized over merger order. Second, a spherical comoving volume shell between any pair of redshifts less than could be contained within the simulation volume. Some Millennium redshift snapshots exist at , and the comoving volume shells between redshifts corresponding to these snapshots enclose some SMBHSMBH coalescence events in the G11 model. An observer located at the center of the Millennium simulation volume would observe only a fraction of the total list of events in the G11 model at , and an observer located elsewhere in the volume would observe a different selection of events. This is not the case for the MillenniumII simulation, where the volume was too small to enclose any comoving volume shells between redshift snapshots. For each realization of the Millennium (but not the MillenniumII) coalescence list, we therefore specified randomlyplaced spherical shells within the simulation box to select binary SMBHs at these redshifts. For coalescence events at , the corresponding comoving volume shells between redshift snapshots were smaller than the Millennium simulation volume, though not enclosed by it. For these coalescence events, we randomly included binaries in the Millennium coalescence list according to probabilities given by the ratios between the volumes of the comoving shells and the Millennium simulation volume. We used 1000 realizations of the Millennium and MillenniumII coalescence lists to form realizations of the binary SMBH distribution .
In generating realizations of the distribution , we assumed that every SMBHSMBH coalescence event in the G11 model
catalogues was the result of a binary SMBH system that had decayed through GW emission. The G11 model included the assumption
that, upon the merger of two galaxies with central SMBHs, the SMBHs coalesced in every case, before accretion onto the newlyformed
SMBH.
The 1000 realizations of were averaged to form a distribution . We fitted with an analytic function which could be used to generate random realizations of the observable binary SMBH population. We did not use realizations of as realizations of the binary SMBH population because the distributions were binned for computational purposes. A fourparameter function,
(19) 
with free parameters , , and , was found to fit well. We performed the fit on the logarithm of the data to approximate linearity in the fitting procedure. The bestfit parameters are given in Table 2. The frequencyexponent was held fixed at , as predicted by Equations 15 and 17.
3.2. Realizations of pulsar ToAs with GWinduced variations
For a set of binary SMBHs drawn from the distribution , we simulated a corresponding timeseries by summing the contributions from each individual binary. Details of the method used to calculate these contributions are presented in Hobbs et al. (2009). For each binary, we randomized over the right ascension and declination, the orbital inclination angle, the orientation of the line of nodes, and the orbital phase angle at the line of nodes. A new publiclyavailable tempo2 plugin, “addAllSMBHBs”, was written to perform this simulation.
For most of the present work, we did not use tempo2 to fit timing model parameters. Instead, we made use of the tools available for spectral analysis of timing residuals. In our simulations, the “timing residuals” corresponded exactly to given the absence of timing model fitting.
We consider it important to emphasize the distinction between the timeseries and the timing residuals resulting from analyses of observed ToA datasets. Consider a set of observed ToAs that exactly match a particular timing model, except for the addition of GWinduced variations (a timeseries). Given that an observer does not actually possess any prior knowledge of the timing model parameters, the observer will fit the model parameters to the ToAs. The resulting timing residuals will not be equivalent to . This is because the variations in the ToAs can alter the apparent pulsar timing parameters. For example, the presence of a timeseries consisting of a sinusoidal signal with a period of one year will alter the apparent pulsar position.
In order to investigate the statistics of given our model for the binary SMBH population, we first used tempo2 to generate 500 ToAs spanning 5 years exactly corresponding to the PSR J04374715 timing model (Manchester et al., 2012). We then added realizations of the timeseries evaluated at the observed ToAs to these datasets, that is, with yr and yr. The pulsar distance was set to 1 kpc, and the position was held fixed for all simulations. As the binary SMBHs used to produce realizations of had randomized positions and orientations, allowing the pulsar position to vary between realizations would not alter our results. The results presented in this paper are not dependent on the timing model used or on the pulsar distance from the Earth. We also added Gaussian white noise variations with 1 ns rms to the ToAs. This white noise component is much smaller than is usually observed in ToA datasets, but was necessary to smooth over machine precision errors.
Parameter  Value 

4.878  
1.72249  
0.3473 
We included GW sources between Hz and Hz in our simulations of . The lower frequency cutoff was chosen to be less than one fifth of . The upper frequency cutoff was chosen to be greater than . We assumed, after previous works, that all GW sources between these frequency cutoffs are nonevolving over a 5year timespan, i.e., they do not evolve in frequency by more than (5 years). The upper bound on the values of sources, , was set by the last stable orbit of binary SMBHs. We identified this bound by fitting a power law to the high edge of the distribution. The lower bound on , , was set by the lowest nonzero value in . This value corresponds to a binary SMBH containing two components at . The distribution included more than GW sources within this domain; the vast computational cost involved makes it impossible to simulate with this many sources. Fortunately, the shape of the distribution was such that, at a given frequency, the highest sources contributed most to (we return to this point below), defined in terms of (using Equation 8) as:
(20) 
The average characteristic strain spectrum derived from is
(21) 
We found a function, , such that
(22) 
Thus, the GW sources in the domain contribute, on average, 90% of at every frequency. Between Hz and Hz, this amounted to sources. We refer to this domain as the “90% domain”. The 90% domain, along with , and , is shown in Figure 1.
We approximated the total number of sources () in the domain between and and between Hz and Hz as constant. For a given realization of , the actual number of sources in the 90% domain is governed by binomial statistics. We therefore drew a (binomial)random number of sources from the 90% domain, and added contributions from each of them to each realization of . We assumed that the sources remaining in the distribution with , contributing on average 10% to at every frequency, resulted in a stochastic contribution to governed by the central limit theorem. We therefore simulated them as described in §2 using the method of simulating ToA variations corresponding to a GW background implemented in the tempo2 plugin GWbkgrd. We simulated sources between Hz and Hz using the tempo2 method, with the characteristic strain spectrum given by . For each realization of , we added contributions from the GW sources drawn from the 90% domain of the distribution, and from the GW sources corresponding to the remaining (on average) 10% of drawn using the tempo2 method. Shifts in the measured pulse frequencies caused by metric perturbations at both the Earth and the pulsar (i.e., the “Earth term” and the “pulsar term”) were included in our simulations.
The top panel of Figure 2 shows in the spectral bins. We also show a characteristic strain spectrum of the form in Equation 1 with . This value of can be taken to be the prediction from the G11 model, as it corresponds to the characteristic strain in the lowest () spectral bin of a 5year dataset. This particular curvature in the curve, also predicted by Wyithe & Loeb (2003), is caused by the frequencydependence of , which represents the bound beyond which binary SMBHs have crossed the last stable orbit. This curvature is caused by different processes from the curvature reported by Sesana et al. (2008). We show the mean characteristic strain spectrum for GWs from all binary SMBHs in the predicted distribution . Sesana et al. (2008) fitted a broken powerlaw to realizations of the characteristic strain spectrum, accounting for various randomizations over the source population. In particular, Sesana et al. (2008) randomized over the existence of “fractional” sources in every frequency bin of a fiducial dataset, and also excluded the strongest single source in every frequency bin in an attempt to isolate the background signal. The smaller number of sources per unit frequency at higher GW frequencies, combined with the greater contributions to the signal from the strongest single sources in frequency bins at higher frequencies, both resulted in the curved characteristic strain spectra presented by Sesana et al. (2008).
The bottom panel of Figure 2 shows the mean numbers of highest sources that contribute 90% and 50% of in these frequency bins. A small number of sources contribute a large fraction of at every frequency. In the frequency bin, the highest sources contribute on average 90% of , and only 30 sources on average contribute 50% of . At frequencies Hz, the strongest source, on average, contributes more than 90% of in each frequency bin. This is a consequence of the shallow powerlaw nature of the distribution of the GW sources in the distribution.
In this work, we compare the Millenniumbased simulations of with simulations of created using the tempo2 method described in §2. To this end, we simulated ToAs as before, but added realizations of corresponding to oscillators simulated using the tempo2 plugin GWbkgrd. These oscillators were simulated as described in §2 with a mean characteristic strain spectrum given by . We refer to this latter method of simulating as Case H09, after Hobbs et al. (2009). Simulations of using GW sources drawn from will be referred to as Case R12 after the present work.
4. Fourierspectral analysis, and results
In this section, we consider the differences between the cases in the distributions of the periodograms, , evaluated for realizations of for a single pulsar. This is motivated by the results in Figure 2, in particular that the number of GW sources per spectral bin that contribute 90% of varies significantly with frequency. The Case R12 simulations are intended to more accurately represent the effects of GWs from binary SMBHs on ToA datasets than the Case H09 simulations. Example realizations of 5year timeseries in both cases are shown in Figure 3. The timeseries appear quite similar: realizations in both cases are dominated by lowfrequency components. Values of up to 1 s are also present in one realization.
Instead of directly measuring , the added white noise component in the simulated ToAs required us to analyze the periodograms of a timeseries, , given by
(23) 
where is a timeseries of Gaussian white 1 ns rms ToA variations as discussed above. The PSD of , , is significantly red, with a spectral index of (see Equations 1, 20), and is expected to dominate the PSD of at low frequencies. We used the generalized leastsquares algorithm described in Coles et al. (2011) to measure the periodograms, of realizations of . This method requires knowledge of the autocovariance function of , which we obtained as in Equation 7 using the inverse DFT of the known PSD of , , given by
(24) 
In the following, we consider the distributions of in the lower spectral bins, where , to be approximately equivalent to the distributions of .
We produced 1000 realizations of in Case R12 and in Case H09, and measured for each realization. In the top panel of Figure 4, we show the averages of measured from each of the Case R12 and Case H09 realizations, along with the expected PSDs of and . Arbitrarily chosen single realizations of in each case are also shown. The means of the periodograms in both cases are clearly the same, and equivalent to the predicted PSD, , given in Equation 24. Though this is as expected, it is both a check of the simulations of , and a demonstration of the ability of our PSD estimation method to measure steep red spectra without bias.
In contrast, the distributions of in the frequency bins where are different between the cases. The single realizations of in each case shown in the top panel of Figure 4 begin to hint at these differences. In most spectral bins, the Case R12 periodogram is below the Case H09 periodogram. That this is a genuine trend is confirmed in the bottom panel of Figure 4. Here, we depict various “percentile” periodograms of the distributions of in each of Case R12 and Case H09. The percentile periodograms may be interpreted as contours of equivalent percentiles of the periodogram distributions in different spectral bins. For example, the “50%” percentile periodogram links the 50th percentile points of the distributions of periodogram values in each spectral bin. Below the 95th percentile, all Case R12 percentile periodograms lie below Case H09 percentile periodograms. This implies that in most spectral bins, most measurements of a periodogram in Case R12 will, like the individual ones shown in the top panel of Figure 4, be below most Case H09 periodograms. However, the 95th percentile periodograms are essentially equivalent, and the maximum value Case R12 periodogram is well above the maximum value Case H09 periodogram. These implied “tails” at high values in the Case R12 periodogram distributions in each spectral bin are highlighted in Figure 5, which depicts the distributions of the in the spectral bins indicated by the vertical lines in the bottom panel of Figure 4. The distributions are shown as the fractions of Case R12 and Case H09 periodograms at or above a given value.
Figure 5 also shows that the Case R12 periodogram distribution in the spectral bin has a longer tail relative to the Case H09 distribution, as compared to the spectral bin. This effect is also evident in the bottom panel of Figure 4, in that the fractional differences between the percentile periodograms are greater at the upper end of the GWdominated frequency regime. This is consistent with the number of GW sources per spectral bin included in the Case R12 simulations going down with increasing frequency, as shown in Figure 2.
In summary, approximating with as discussed above, we find that:

In most spectral bins, most realizations of in Case R12 will be below most realizations of in Case H09.

The maximum possible values of in Case R12 will be higher than the maximum possible values of in Case H09.
5. Estimates of correlations between GWinduced ToA variation timeseries
Hellings & Downs (1983) showed that the average values of correlations between for different pulsars, for a stochastic GW signal, will always be given by the Hellings & Downs curve (Equation 4). The distributions of individual estimates of these correlations, much like the distributions of the PSD estimator considered above, will however depend on the nature of the GW signal. Here, we characterize the distributions of estimates of these correlations for multiple pulsar pairs in each case discussed above. We simulated 100 realizations of each of Case H09 and Case R12 ToAs as described in §3.2 for pulsars at the positions of each of the 20 pulsars timed by the Parkes PTA, using the timing models specific to each pulsar (Manchester et al., 2012). For each realization, the same set of GW sources was used to simulate GWinduced ToA variations for each pulsar. The pulsar distances were set arbitrarily between 1 kpc and 20 kpc.
We estimated the correlations between timeseries and , , for each pulsar pair in each realization of Case R12 and Case H09 ToAs. No autocorrelations were estimated. A frequencydomain estimation technique, based on the method of Yardley et al. (2011), was used. This technique will be detailed elsewhere (Hobbs et al., in preparation). We refer to our estimates of as . These estimates were performed using timeseries, rather than timeseries (see Equation 23), and the estimation technique was optimized using the expected PSD of given in Equation 24. The technique removes bestfit linear and quadratic terms from each timeseries using the standard tempo2 leastsquares fitting algorithm. This mimics the effect of fitting pulse frequency and frequencyderivative terms to the simulated ToAs and then analyzing the timing residuals.
Following the removal of linear and quadratic terms from one of the 100 Case R12 realizations of , a single GW source was found to dominate the residual timeseries. We show the corresponding timeseries for two of the 20 simulated pulsars for this realization in Figure 6. The left panel of this Figure shows the large sinusoidal oscillations induced by the source, and the right panel shows example Case R12 realizations of that are not dominated by an individual source. It is possible that an individual GW source with a period greater than the 5year dataspan could dominate the realizations of in the righthand panel of Figure 6. The ToA variations induced by such a source would however be absorbed in the removal of the linear and quadratic terms from the timeseries.
We averaged all measurements of for each pulsar pair from the Case R12 realizations, besides the one clearly dominated by an individual source. The realization dominated by an individual source added a large amount of scatter to the average Case R12 correlations, and was left out of the average to enable a better comparison between the cases. We also averaged the Case H09 measurements of for each pair in 99 arbitrarilychosen realizations. The average measurements of are shown for both cases in Figure 7. The functional form of the Hellings & Downs curve is recovered in both Case R12 and Case H09. However, the Case R12 estimates are significantly more scattered about the expected values of the correlations than the Case H09 estimates.
The increased scatter in the Case R12 correlations with respect to the Case H09 correlations in Figure 7 is caused by outlying estimates in only a few realizations of ToAs. This is shown in Figure 8, where we display the histograms of the measurements between simulated ToA datasets for PSR J04374715 and PSR J06130200 in each case. Correlation estimates were possible because we normalized the estimated covariances between timeseries using the expected crossPSD between the timeseries. While most measurements in both cases are concentrated about the expected value of , a few Case R12 measurements are significantly displaced. This is consistent with the results of §4. We also stress that the large scatter of the estimator common to both cases is expected, and intrinsic to the GW signal.
6. Discussion
We have shown that the ToA variations induced by GWs from the predicted binary SMBH population are not consistent with the model described in §2. The random Gaussian model for described in §2 and approximated in Case H09 is reasonable given that a large number of GW sources are expected to contribute to the GWinduced ToA variations. That is, the values of at all times are the sums of many random variables. An argument based on the classical central limit theorem would suggest that would then be Gaussian random at every time . It is apparent, however, that such a central limit theorembased argument does not apply to the Case R12 realizations of . This is because of the nature of the GW sources contributing to in Case R12.
In Case R12, a few sources contribute most of the PSD of at every frequency, as shown in Figure 2. These sources are rare, because they are found at the high tail of the source distribution. The estimators we consider in this work, and , are dominated in Case R12 by a few GW sources that need not occur in every realization of the timeseries. This is why the distributions of these estimators are different between the cases. The quantities that we estimate, and , are used to define the covariance matrix of the GWinduced ToA variations (see Equation 7). We have therefore shown that the ToA variations induced by GWs from binary SMBHs are dominated by the effects of a few strong, rare sources and cannot be accurately modeled using the random Gaussian process discussed in §2.
6.1. Implications of our results for experiments focused on a GW background
Current PTA data analysis techniques use assumptions about the statistics of to attempt to estimate or constrain the amplitude of the
characteristic strain spectrum of GWs from binary SMBHs. We consider the implications of our results for a selection of techniques in turn. We assume, in this
discussion, that our results for the statistics of GWinduced ToA variations would apply even if the normalization of the GW characteristic strain spectrum
, which we refer to as the GW amplitude,
We summarize a few key techniques here:

Jenet et al. (2005) describe a statistic which measures the degree of correlation between estimates of from ToA data, and the expected functional form of . The expected detection significance, which is estimated under the assumption that the GWinduced ToA variations are Gaussian random, saturates at high values of the GW amplitude once the variance of the statistic is dominated by the stochasticity of the GW signal (see our Figure 8 and related discussion).

Jenet et al. (2006) constrain the amplitude of the GW characteristic strain spectrum from binary SMBHs by estimating the maximum possible GW signal present in measured data, under the assumption that the data could be modeled using a white noise process and GWinduced ToA variations. A statistic that estimates the GW background amplitude from individual pulsars was measured, and compared to the simulated distributions of the statistic for different GW amplitudes. The simulated statistic distributions were created from simulated ToA datasets with GWinduced ToA variations included using the tempo2 plugin GWbkgrd.

van Haasteren et al. (2009) present a Bayesian parameter estimation method for the GW characteristic strain spectral amplitude.
^{8} van Haasteren et al. (2011) used this method to constrain the GW amplitude. This method requires an evaluation of the likelihood of the parameters used to model the ToA datasets, which include the GW amplitude. The likelihood is the probability distribution of the data given the model parameters. The GW amplitude is used to calculate the covariance matrix, , between the GWinduced ToA variations for pulsars and (see Equations 7, 8). This covariance matrix in turn is used to define the PTA likelihood, assuming that the GWinduced ToA variations can be modeled as a random Gaussian process. 
Demorest et al. (2012) use a PTA likelihood similar to the work of van Haasteren et al. (2009) to constrain the GW amplitude by evaluating the distribution of a maximum likelihood estimator for the amplitude. They also use a method similar in concept to Jenet et al. (2005) to attempt to detect the GW signal from binary SMBHs.
First, the nonGaussianity of the GWinduced ToA variations means that the estimate of the intrinsic GWinduced variance of the Jenet et al. (2005) statistic will be incorrect. This will affect estimates of the detection significance, particularly in the “strong signal regime”, where the effects of GWs in the ToAs are large compared to all other noise processes. Second, the limit on the GW amplitude placed by Jenet et al. (2006) will be biased. The Jenet et al. (2006) limit was placed by finding the GW amplitude for which 95% of simulated statistic values were above the measured value. The distribution of their statistic derived using our simulations would be different. Ruling out a GW amplitude using the Jenet et al. (2006) technique does not necessarily rule out a GW signal corresponding to our simulations with the same confidence. Finally, our results indicate that the likelihoods evaluated by van Haasteren et al. (2009) and Demorest et al. (2012) will also be biased, leading to a similar effect on GW amplitude constraints made using their methods. A definitive statement on the magnitude of the consequences of our simulations for current constraints on the GW amplitude cannot be made, however, without fully considering the various PTA data analysis methods.
6.2. Single GW source detection prospects
We have established that in every frequency bin of the 5year datasets we consider, a few strong GW sources dominate the
PSD of . This means that we cannot consider the expected GW signal from binary SMBHs to form a
background.
The exact number of resolvable GW sources given a GW background level for PTAs depends on the particular search method. For example, Boyle & Pen (2010) suggest that a PTA composed of pulsars could resolve up to sources per frequency bin. In Figure 9, we present a simple indication of the expected amplitudes of strong individual sources in the spectral bins of our fiducial 5year dataset. Using 300 Case R12 realizations of the GW source population, we found the mean strain amplitudes in each spectral bin of the strongest three GW sources. We express these amplitudes as multiples of the mean summed strain amplitude, , of the remaining sources. The errors in the values were not included in the error bars as they were very small.
If we consider the sources besides the strongest three in a spectral bin to form a “background”,
Blind searches for single GW sources with PTAs are therefore important. PTA data analysis methods that attempt to detect an isotropic component will not optimally recover the entirety of the GW signal from binary SMBHs, and could perhaps miss a large component of the signal for some realizations of the GW source population. A careful consideration of the efficacy of GW background detection methods as compared to search methods for single sources, given the predicted source characteristics, is required.
6.3. Limitations of our approach to modelling the GW signal from binary SMBHs
A shortcoming of our approach towards modeling , and indeed of all predictions for the GW signal from binary SMBHs to date, is the assumption of circular orbits for all binaries. Recent work (e.g., Sesana, 2010; Preto et al., 2011; Khan et al., 2011) suggests that binary SMBHs emitting GWs in the PTA frequency regime will have highly eccentric orbits. The candidate binary SMBH OJ 287 (e.g., Valtonen et al., 2011) is in fact modeled with an orbital eccentricity of 0.7. The GW waveform of an eccentric binary radiating in the PTA band spans many frequencies, and does not follow the frequencytime relation of Equation 15, or the dependence of the GW strain amplitude on frequency of Equation 14. Therefore, if most binary SMBHs radiating in the PTA band are eccentric, the predicted mean spectral slope of the characteristic strain spectrum (Equation 1) will change. We will also need to account for the binary eccentricity distribution in our predictions of the statistics of GWinduced ToA variations. Other effects that require further investigation are the effects of gas and stars on binaries. We do not include these effects in our model because the understanding of their combined contributions towards specifying the binary SMBH population is not sufficiently advanced.
7. Conclusions
We have used a sophisticated model for galaxy evolution (Guo et al., 2011) to predict the distribution of binary SMBHs radiating GWs in the PTA frequency band. By drawing lists of GW sources from this distribution, we simulate the effects of GWs from binary SMBHs on 5year pulsar ToA datasets. We compare these simulations (Case R12) with simulated pulsar datasets containing the effects of an equivalentamplitude GW signal modeled as a random Gaussian process (Case H09). We estimate the PSDs of the simulated GWinduced ToA variation timeseries, and the correlations between these timeseries for different pulsars. We find that the distributions of the PSD estimators of the realizations of the GWinduced ToA variations are different between the cases in every frequency bin, although the mean estimated PSDs are the same in each case. While in Case R12 the estimated PSDs are concentrated at lower values than in Case H09, the Case R12 estimations extend to higher PSD values than the Case H09 estimations. We also find that the functional form of the Hellings & Downs curve is recovered on average in both cases. The correlations between the GWinduced ToA variation timeseries for different pulsars in Case R12 are, however, significantly more scattered about the expected values than in Case H09. We interpret our results in terms of the influence of strong individual GW sources on the ToAs in Case R12.
We conclude the following:

The effects of GWs from binary SMBHs on pulsar ToAs cannot be accurately modeled using existing methods, i.e, as a random Gaussian process. This is because a few GW sources dominate the PSD of the GWinduced ToA variations at all frequencies, with reducing numbers of sources contributing equivalent PSD fractions in higher frequency bins.

Our results directly affect existing PTA data analysis methods aimed at detecting or estimating the parameters of the GW signal from binary SMBHs. The projected detection significance will be biased.

The prospects for single GW source detection are strong. Individual sources could potentially be resolved in all GWdominated frequency bins of a 5year dataset.
We emphasize that future searches for GW signals from binary SMBHs in pulsar datasets need to be sensitive to both individual sources as well as a GW background.
Footnotes
 affiliation: Also at CSIRO Astronomy and Space Science, Australia Telescope National Facility, PO Box 76, Epping, NSW 1710, Australia
 affiliation: Email address: v.vikram.ravi@gmail.com
 This mass range approximately corresponds to the current sample of SMBHs with dynamical mass measurements (cf. McConnell et al., 2011).
 The MillenniumII simulation however does not reproduce largerscale structures as well as the Millennium simulation.
 http://www.mpagarching.mpg.de/millennium/
 There are various mechanisms by which extreme mass ratio binary SMBH systems and triple or higherorder systems can avoid coalescence (e.g., Volonteri et al., 2003).
 We make a distinction between this amplitude and the parameter introduced in Equation 1, because does not have exactly the same form as given in Equation 1.
 Though their method also estimates the spectral index of the GW characteristic strain spectrum, we assume marginalization over this parameter in our discussion here.
 This result is analogous to the case of the extragalactic background light (Domínguez et al., 2011), where the summed electromagnetic radiation from AGN and starforming galaxies is dominated by strong individual sources, behind which myriad further objects combine to form an apparently isotropic background too uniform to be resolved by current telescopes.
 This is by no means a rigorous definition of a background relative to the number of sources. The exact definition is dependent on the single source search method and the characteristics of the PTA.
References
 Babak, S., & Sesana, A. 2012, Phys. Rev. D, 85, 044034
 BoylanKolchin, M., Springel, V., White, S. D. M., Jenkins, A., & Lemson, G. 2009, MNRAS, 398, 1150
 Boyle, L., & Pen, U.L. 2010, arXiv:1010.4337
 Coles, W., Hobbs, G., Champion, D. J., Manchester, R. N., & Verbiest, J. P. W. 2011, MNRAS, 418, 561
 Colless, M., Dalton, G., Maddox, S., Sutherland, W., Norberg, P., Cole, S., et al. 2001, MNRAS, 328, 1039
 Corbin, V., & Cornish, N. J. 2010, arXiv:1008.1782
 Croton, D. J., Springel, V., White, S. D. M., De Lucia, G., Frenk, C. S., Gao, L, et al. 2006, MNRAS, 365, 11
 De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2
 Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., Nice, D., Ransom, S., Stairs, I. H., et al. 2012, arXiv:1201.6641
 Detweiler, S. 1979, ApJ, 234, 1100
 Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33
 Domínguez, A., Primack, J. R., Rosario, D. J., Prada, F., Gilmore, R. C., Faber, S. M., et al. 2011, MNRAS, 410, 2556
 Dotti, M., Sesana, A., & Decarli, R. 2012, Advances in Astronomy, 2012
 Ellis, J., Siemens, X., & Creighton, J. 2012, arXiv:1204.4218
 Foster, R. S., & Backer, D. C. 1990, ApJ, 361, 300
 Gitti, M., Brighenti, F., & McNamara, B. R. 2012, Advances in Astronomy, 2012
 Guo, Q., White, S., BoylanKolchin, M., De Lucia, G., Kauffmann, G., Lemson, G., et al. 2011, MNRAS, 413, 101
 Hellings, R. W., & Downs, G. S. 1983, ApJ, 265, L39
 Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655
 Hobbs, G. B., Jenet, F., Lee, K. J., Verbiest, J. P. W., Yardley, D. R. B., Manchester, R. N., et al. 2009, MNRAS, 394, 1945
 Hughes, S. A. 2002, MNRAS, 331, 805
 Jaffe, A. H., & Backer, D. C. 2003, ApJ, 583, 616
 Jenet, F. A., Hobbs, G. B., Lee, K. J., & Manchester, R. N. 2005, ApJ, 625, L123
 Jenet, F. A., Hobbs, G. B., van Straten, W., Manchester, R. N., Bailes, M., Verbiest, J. P. W., et al. 2006, ApJ, 653, 1571
 Khan, F. M., Just, A., & Merritt, D. 2011, ApJ, 732, 89
 Lee, K. J., Wex, N., Kramer, M., Stappers, B. W., Bassa, C. G., Janssen, G. H., et al. 2011, MNRAS, 414, 3251
 Lemson, G., & Virgo Consortium 2006, arXiv:astroph/0608019
 Manchester, R. N., et al. 2012, Proc. Astr. Soc. Aust., submitted
 Marulli, F., Bonoli, S., Branchini, E., Moscardini, L., & Springel, V. 2008, MNRAS, 385, 1846
 McConnell, N. J., Ma, C.P., Gebhardt, K., Wright, S. A., Murphy, J. D., Lauer, T. R., et al. 2011, Nature, 480, 215
 Phinney, E. S. 2001, arXiv:astroph/0108028
 Preto, M., Berentzen, I., Berczik, P., & Spurzem, R. 2011, ApJ, 732, L26
 Sazhin, M. V. 1978, Soviet Ast., 22, 36
 Schuster, A. 1898, Terrestrial Magnetism and Atmospheric Electricity, 3, 13
 Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192
 Sesana, A., Vecchio, A., & Volonteri, M. 2009, MNRAS, 394, 2255
 Sesana, A. 2010, ApJ, 719, 851
 Spergel, D. N., Verde, L., Peiris, H. V., Komatsu, E., Nolta, M. R., Bennett, C. L., et al. 2003, ApJS, 148, 175
 Springel, V., White, S. D. M., Jenkins, A., Frenk, C. S., Yoshida, N., Gao, L.,et al. 2005, Nature, 435, 629
 Thorne, K. S. 1987, Three Hundred Years of Gravitation, 330
 Valtonen, M. J., Lehto, H. J., Takalo, L. O., & Sillanpää, A. 2011, ApJ, 729, 33
 van Haasteren, R., Levin, Y., McDonald, P., & Lu, T. 2009, MNRAS, 395, 1005
 van Haasteren, R., Levin, Y., Janssen, G. H., Lazaridis, K., Kramer, M., Stappers, B. W.,et al. 2011, MNRAS, 414, 3117
 Verbiest, J. P. W., Bailes, M., Coles, W. A., Hobbs, G. B., van Straten, W., Champion, D. J., et al. 2009, MNRAS, 400, 951
 Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559
 Wyithe, J. S. B., & Loeb, A. 2003, ApJ, 590, 691
 Yardley, D. R. B., Hobbs, G. B., Jenet, F. A., Verbiest, J. P. W., Wen, Z. L., Manchester, R. N.,et al. 2010, MNRAS, 407, 669
 Yardley, D. R. B., Coles, W. A., Hobbs, G. B., Verbiest, J. P. W., Manchester, R. N., van Straten, W., et al. 2011, MNRAS, 414, 1777