GRAVITATIONAL WAVES FROM INDIVIDUAL SUPERMASSIVE BLACK HOLE BINARIES IN CIRCULAR ORBITS: LIMITS FROM THE NORTH AMERICAN NANOHERTZ OBSERVATORY FOR GRAVITATIONAL WAVES

It was shown in van Haasteren & Levin (2013) (and modified in Ellis 2013 to include deterministic sources) that the likelihood function for the residuals, marginalized over the timing model parameters, can be written as

$$\begin{equation} p(\delta t|boldsymbol {\phi },boldsymbol {\lambda })=\frac{\exp \left[- \frac{1}{2} \left(\delta t -s \right)^T G(G^T CG)^{-1} G^T\left(\delta t -s\right)\right]}{\sqrt{(2\pi)^{(N_{\rm TOA}-m)}\det (G^TCG)}}, \end{equation} \tag{ 18 }$$

where δt is a vector of the pulsar residuals, $boldsymbol {\phi }$ are parameters that describe the noise in the pulsar residuals, $boldsymbol {\lambda }$ are the parameters that characterize the continuous GW signal,²⁹ and C is the covariance matrix of the noise. In practice, the noise in pulsar timing residuals is non-Gaussian due to ISM scintillation effects which are manifested through a time varying pulse intensity. Nonetheless, the noise in each residual is modeled very well by a Gaussian with zero mean and standard deviation equal to the uncertainty on the TOA. In other words, the noise in the weighted (by the individual TOA errors) residuals is very well approximated as a Gaussian. As will be detailed in the next section, we include these error bar weights in our noise covariance matrix. Thus, assuming the noise is Gaussian³⁰ is justified. Here G is an N_TOA × (N_TOA − m) matrix with N_TOA the number of TOAs and m the number of fitted parameters in the timing model. The derivation of the G-matrix approach can be found in van Haasteren & Levin (2013) and will not be explored here. The matrix G^T is a projection operator that projects our data δt onto the null space of the timing model design matrix M, that is, it projects the data into a subspace orthogonal to the linearized timing model. In the timing analysis used here, DM variation and profile frequency evolution effects are part of the timing model, and these terms are included when constructing the G matrix. In this way we have fully taken into account the timing model fitting procedure.

For this work we will assume that the noise in the residuals between pulsars is uncorrelated. In other words, we are assuming that the stochastic GW background will be negligible compared to the intrinsic noise in each pulsar. In general this is not likely to be a good assumption when we would expect a detection of a single GW source. Furthermore, terrestrial clock errors (Hobbs et al. 2012) and errors in solar system ephemerides (Champion et al. 2010)³¹ can also cause correlations between the noise in the residuals from different pulsars, with however different angular correlation properties on the sky than are expected from GWs. The effects of omitting the correlations in the likelihood function are unknown and will be the subject of future work. Under these assumptions, the likelihood function for the full PTA can be written as

$$\begin{equation} p(\delta t|boldsymbol {\lambda })=\prod _{\alpha =1}^{N_{\rm psr}}p(\delta t_{\alpha }|boldsymbol {\lambda }_{\alpha }), \end{equation} \tag{ 19 }$$

where δt_α and $boldsymbol {\lambda }_{\alpha }$ and the residuals and model parameters for the αth pulsar, respectively, and $boldsymbol {\lambda }$ is the full continuous wave (CW) parameter vector including pulsar distances for all pulsars. In cases where we fix the noise values, we can write the log-likelihood ratio of a model with a single continuous GW to a model with just noise as

$$\begin{equation} \ln \,\Lambda =\sum _{\alpha }^{N_{\rm psr}}\left[ \left(\delta t_{\alpha }|s(boldsymbol {\lambda }_{\alpha })\right)-\frac{1}{2}\left(s(boldsymbol {\lambda }_{\alpha })|s(boldsymbol {\lambda }_{\alpha })\right)\right], \end{equation} \tag{ 20 }$$

where the inner product between two time series x and y is

$$\begin{equation} (x|y)=x^TG(G^T CG)^{-1} G^Ty. \end{equation} \tag{ 21 }$$

In the remainder of the paper we will refer to the signal-to-noise ratio in the following form

$$\begin{equation} \rho =\sqrt{2\langle \ln \,\Lambda \rangle }= \left(\sum _{\alpha =1}^{N_{\rm psr}} (s(boldsymbol {\lambda }_{\alpha })| s(boldsymbol {\lambda }_{\alpha }) )\right)^{1/2}, \end{equation} \tag{ 22 }$$

where the brackets denote the expectation value over many noise realizations.

In the above section we have derived the likelihood function used for our analysis; however, we have not specified the form of the noise covariance matrix C. Previous Bayesian analysis schemes (van Haasteren et al. 2009, 2011; van Haasteren & Levin 2010, 2013; Ellis et al. 2012; van Haasteren 2013; Ellis et al. 2013; Ellis 2013) have used a power-law red noise model and EFAC (constant multiplier on the TOA uncertainties) and EQUAD (additional Gaussian white noise added in quadrature to EFAC noise) parameters to describe the white noise, with a covariance matrix of the form

$$\begin{equation} C = E^2W + \mathcal {Q}^2\mathbb {I} + C^{\rm red}(A_{\rm red}, \gamma _{\rm red}), \end{equation} \tag{ 23 }$$

where E is the EFAC parameter, $W={\rm diag}\lbrace \sigma _i^2\rbrace$, with σ_i the uncertainty on the ith TOA, $\mathcal {Q}$ is the EQUAD parameter, and C_red is an analytic expression of the red noise amplitude A_red and spectral index γ_red. It is worth noting that we use no EFAC or EQUAD parameters in our pulsar timing model fit but instead include them directly in our noise model. The EFAC is simply a parameter that quantifies any additional uncertainty in the TOA uncertainties and the EQUAD parameter quantifies any additional white noise that is not related to the formal TOA uncertainties. In principle, a different EFAC value should be used for each pulsar timing backend as this parameter is related to intrinsic receiver noise; however, in this five-year NANOGrav data set, only one backend per telescope was used.³² Therefore, we are justified in only using one EFAC parameter per pulsar. This noise model is quite general and works well for many pulsars; however, the size of the matrices is quite large (on the order of 10³ × 10³) and inversion is a large bottleneck in the analysis pipelines. Furthermore, current NANOGrav observing schemes produce large sets of multifrequency observations that are essentially simultaneous. One may be tempted to simply perform a weighted average of the TOAs and work with the new reduced data sets but in the Bayesian scheme we must marginalize over the timing model parameters analytically and it is unclear how to carry out this process for epoch-averaged TOAs. Because of this, we have developed a framework to essentially work backward from the marginal likelihood to derive a nearly exact averaging scheme. First we re-write our noise covariance matrix

$$\begin{equation} C = N + U\tilde{\mathcal {C}}U^T, \end{equation} \tag{ 24 }$$

where $\tilde{\mathcal {C}}$ is a q × q reduced covariance matrix with q the number of epochs³³ in our data set, N is a white noise covariance matrix of the EFAC and EQUAD terms, and U is the "exploder" matrix that maps epochs (columns) to the full set of TOAs (rows). If we now make use of this new formalism, the likelihood function is then

$$\begin{equation} p(\delta t|boldsymbol {\phi },boldsymbol {\lambda }) = \frac{\exp \left[ -\frac{1}{2} ((\delta t-s)^T \tilde{N}^{-1} (\delta t-s) - d^T \Sigma ^{-1}d)\right]}{\sqrt{(2\pi)^{n-m}\det (\tilde{\mathcal {C}}) \det (G^TNG) \det (\Sigma)}}, \end{equation} \tag{ 25 }$$

where we have used the Woodbury Lemma³⁴ to compute the inverse and determinant of C, $\tilde{N}^{-1} = G(G^T N G )^{-1} G^T$, $d = U^T\tilde{N}^{-1}(\delta t-s)$, and $\Sigma = (\tilde{\mathcal {C}}^{-1} + U^T\tilde{N}^{-1}U )$. Note that d here are essentially daily averaged residuals. For NANOGrav data sets the number of epochs per pulsar is on the order of 30–100, while the total number of TOAs per pulsar is on the order of 10³, thus the inversions (here N is diagonal and $\tilde{N}^{-1}$ can be pre-computed, thus the only dense matrix inversion is Σ⁻¹) required in this likelihood function scale as q³ as opposed to n³, resulting in computational speed ups of several orders of magnitude. Furthermore, the epoch-averaged covariance matrix $\tilde{\mathcal {C}}$ can take on several forms depending on the red noise model used; however, as long as it is a slowly varying function of the TOAs (i.e., a truly red process) then this formalism is completely valid.

In order to attain further computational speed ups and to gain more control over the low frequency component of our noise model we make use of the methods described in Lentati et al. (2013), but now applied to a single pulsar instead of the full PTA. This method relies on explicitly splitting up the red and white components of the residuals, so that the residuals are now written as

$$\begin{equation} \delta t = M\delta \boldsymbol{\xi }+ n_{\rm white}+ n_{\rm red} + s, \end{equation} \tag{ 26 }$$

where M is the design matrix (an N_TOA × m matrix composed of the functional form of the linearized timing model (column) evaluated at each TOA (row)), n_white and n_red are the white and red components of the residuals, respectively, and $\delta \boldsymbol{\xi }$ are the timing model parameter offsets. It is possible to expand the red noise piece in a Fourier series

$$\begin{equation} n_{\rm red} = \sum _{j=1}^{N_{\rm mode}}\left[ a_j \sin \left(\frac{2\pi j t}{T}\right)+ b_j \cos \left(\frac{2\pi j t}{T}\right)\right] = F\mathbf {a}, \end{equation} \tag{ 27 }$$

where a is a vector of the concatenated sine and cosine amplitudes, T is the total time span of the data, and F is an N_TOA × 2N_mode matrix with alternating sine and cosine terms with N_mode the number of frequencies used. Now, we assume that the underlying ensemble average red noise process is wide-sense stationary and can be completely described by a power spectrum. Then, by orthogonality, the Fourier coefficients a will be diagonal with components

$$\begin{equation} \varphi _{ij} = \langle \mathbf {a} \mathbf {a}^T\rangle _{ij} = \delta _{ij}\varphi _{i}, \end{equation} \tag{ 28 }$$

where there is no sum over i and the elements of φ, denoted {φ_i} are the coefficients of the model power spectrum of the red noise process in the residuals. If the red noise process is wide-sense stationary, then this relation is always true irrespective of the sampling as all information about the uneven sampling here comes from the Fourier design matrix F. Thus, we can write the covariance and epoch-averaged covariance matrices, respectively, as

$$\begin{equation} C = N + F\varphi F^T \end{equation} \tag{ 29 }$$

$$\begin{equation} \tilde{\mathcal {C}} = \tilde{F}\varphi \tilde{F}^T, \end{equation} \tag{ 30 }$$

where $\tilde{F}$ is a q × N_mode matrix and is constructed in the same manner as F but the epoch-averaged TOAs are used as opposed to the full set of TOAs. As is done in Lentati et al. (2013), it is possible to treat each diagonal element of φ as a free parameter; however, for this work we chose to parameterize it by a power law

$$\begin{equation} \varphi _i = \frac{1}{T}\frac{A_{\rm red}^2}{12 \pi ^2}\left(\frac{f_i}{f_{\rm yr}} \right)^{3-\gamma _{\rm red}} f_i^{-3}, \end{equation} \tag{ 31 }$$

where f_i is the ith Fourier frequency assuming Nyquist sampling. In general, any Fourier-based method with finite length data sets and especially with irregular sampling will suffer from spectral leakage whereby power from the lowest frequencies will leak into the higher frequencies. In effect, this makes Fourier-based methods sensitive to the low-frequency cutoff. However, it was shown in van Haasteren & Levin (2013) that including the effects of the timing model (specifically the quadratic spin-down in this case) in our analysis acts as a window function that fully removes any sensitivity to the low-frequency cutoff, thereby also removing nearly all spectral leakage. We have done extensive simulations to test this notion and have found no evidence for spectral leakage and no bias in red noise parameter estimation and waveform reconstruction.

In the course of our single pulsar noise analysis (J. A. Ellis et al. 2014, in preparation) we found that the addition of an extra white noise parameter was needed to accurately describe the data. This new white noise term incorporates a correlation among frequency channels (within a given epoch) while still remaining independent of other epochs. In other words, this white noise term accounts for epoch-to-epoch fluctuations as opposed to fluctuations within an epoch. We defer to another paper the inclusion of pulse-jitter noise from pulsar magnetospheric activity but point out that our inferred extra term may be the same as jitter noise known to be present in all well-studied pulsars (Cordes & Shannon 2010). This parameter is quite easy to incorporate as it is simply an EQUAD-like parameter in the epoch-averaged sense, that is,

$$\begin{equation} J = U\tilde{J}U^T = \mathscr{J}^2U\mathbb {I}_{q}U^T, \end{equation} \tag{ 32 }$$

where $\mathscr{J}$ is our frequency correlated EQUAD parameter and $\mathbb {I}_q$ is the identity matrix in the epoch-averaged space. With this, we have our final noise model with a total covariance matrix of

$$\begin{equation} C = N + U\left(\tilde{F}\varphi \tilde{F}^T + \mathscr{J}^2\mathbb {I}_q \right)U^T, \end{equation} \tag{ 33 }$$

and noise parameter vector

$$\begin{equation} boldsymbol {\phi } = \lbrace E, \mathcal {Q}, \mathscr{J}, A_{\rm red}, \gamma _{\rm red}\rbrace. \end{equation} \tag{ 34 }$$

Throughout the remainder of the paper, this noise model is always used for all pulsars.

The $\mathcal {F}_p$ statistic was first derived from the likelihood function of Equation (18) in Ellis et al. (2012, hereafter ESC12) as a "total-power" frequentist detection statistic. We will not derive the full expression here (see Appendix B for an alternate derivation to ESC12), but rather we will explain its functional form and discuss its statistics. First we define the following harmonic basis functions:

$$\begin{equation} B_{\alpha }^{1}(t)=\frac{1}{\omega _{0}^{1/3}}\sin (2\omega _{0}t) \end{equation} \tag{ 35 }$$

$$\begin{equation} B_{\alpha }^{2}(t)=\frac{1}{\omega _{0}^{1/3}}\cos (2\omega _{0}t), \end{equation} \tag{ 36 }$$

where, again, ω₀ is the orbital angular frequency of the SMBHB. Following ESC12, the $\mathcal {F}_p$ statistic is written as

$$\begin{equation} 2\mathcal {F}_p = \sum _{\alpha =1}^{N_{\rm psr}}P^{i}_{\alpha }Q_{ij}^{\alpha }P_{\alpha }^{j}, \end{equation} \tag{ 37 }$$

where we have assumed Einstein summation notation over latin indices, $P^{i}_{\alpha } = (\delta t|B^i_\alpha (t))$, $Q_{ij}^{\alpha }=(B^{\alpha }_{i}|B^{\alpha }_j)$ and the formal sum is over all pulsars in the array. An intuitive way to think of this statistic is a weighted (by the noise power spectral density) sum of the power spectrum of the residual data done in the time domain by making use of a harmonic time domain basis. It was shown in ESC12 that $2\mathcal {F}_p$ follows a χ² distribution with 2N_psr degrees of freedom and non-centrality parameter ρ² such that

$$\begin{equation} \langle 2\mathcal {F}_p \rangle = 2N_{\rm psr}+ \rho ^2. \end{equation} \tag{ 38 }$$

Figure 1 shows the distribution of the $\mathcal {F}_p$ statistic (top panel) for both real and simulated data as well as a p-value test (bottom panel) for each pulsar, where we compare the single-pulsar $\mathcal {F}_p$ distribution to the expected χ² distribution. To compute the $\mathcal {F}_p$-statistic, we have used the maximum a posteriori noise values obtained in a previous single-pulsar noise analysis to construct the noise covariance matrix. Since we do not have independent realizations of our data, we compute the $\mathcal {F}_p$-statistic for each independent³⁵ frequency bin and then construct a histogram of the results. If our noise model is a good description of the true noise in our data and there is no GW present in the data then this distribution should follow the correct χ² distribution. We see from Figure 1 that the $\mathcal {F}_p$-statistic values do indeed follow a χ² distribution with 2N_psr degrees of freedom. The black (blue) curve in the top panel of Figure 1 shows the aforementioned histogram along with the χ² distribution in the dashed gray (red) line. The p-value that results from a Kolmogorov–Smirnov (KS) test comparing the $2\mathcal {F}_p$ and χ² (with 34 degrees of freedom) distributions is 0.33, showing good agreement between our data and the expected χ² distribution. As a cross-check, we have also simulated 100,000 data sets with the measured noise parameters and have evaluated the $\mathcal {F}_p$ statistic for each. This distribution is plotted as a gray (green) histogram in the figure and it is obvious that this distribution follows a χ² distribution with 34 degrees of freedom nearly perfectly. We have also performed a similar test but for each pulsar separately. In the lower panel of Figure 1 we carry out the same KS test mentioned above but now compute the $\mathcal {F}_p$ statistic values for each pulsar individually and then compare to a χ² with two degrees of freedom. The solid line corresponds to the p-value at which we should reject the null hypothesis that the two distributions are the same with 99.7% confidence. We see that with the exception of one pulsar, J1640+2224, all others lie above this threshold value. This indicates that our noise model for all pulsars, except for J1640+2224, provides a good description of the true noise in the data set. Better noise models for this pulsar are currently being explored (J. A. Ellis et al. 2014, in preparation). Since our full 17 pulsar $\mathcal {F}$-statistic distribution is totally consistent with the expected χ² distribution and PSR J1640+2224 does not bias the overall PTA distribution, we use the standard noise model described in Section 4.2.

Zoom In Zoom Out Reset image size

For the detection problem, we are interested in the false-alarm-probability (FAP), that is, the probability that a measured value $\mathcal {F}_p$ exceeds a given threshold $\mathcal {F}_{p,0}$ when no signal is present. From ESC12, the probability distribution of $\mathcal {F}_p$ when the signal is absent is

$$\begin{equation} p_{0}(\mathcal {F}_p)=\frac{\mathcal {F}_p^{n/2-1}}{(n/2-1)!}\exp (-\mathcal {F}_p), \end{equation} \tag{ 39 }$$

where n is the number of degrees of freedom (2N_psr in this case). The corresponding FAP is then written as

$$\begin{equation} P_{F}(\mathcal {F}_{p,0})=\int _{\mathcal {F}_{p,0}}^{\infty }p_{0}(\mathcal {F}_p)d\mathcal {F}_p =\exp (-\mathcal {F}_{p,0})\sum _{k=0}^{n/2-1}\frac{\mathcal {F}_{p,0}^{k}}{k!}. \end{equation} \tag{ 40 }$$

In a search over GW frequencies (the only free parameter in the $\mathcal {F}_p$-statistic) we will incur a trials factor such that the resulting FAP for the search is

$$\begin{equation} P_{F}^{T}(\mathcal {F}_{p,0})=1-[ 1-P_{F}(\mathcal {F}_{p,0})]^{N_{f}}, \end{equation} \tag{ 41 }$$

where N_f is the number of independent frequencies. For this work we place our detection threshold on $\mathcal {F}_p$ such that the corresponding FAP is less than 10⁻⁴. The results of performing this search on the five-year NANOGrav data set will be presented in the next section.

The Bayesian search pipeline in this work is very similar to that of Ellis (2013, hereafter E13). Here we use an MPI-enabled Parallel Tempered Markov Chain Monte Carlo (PTMCMC) sampler³⁶ (see E13 and Appendix C.1 for details on the implementation). In this work we use two "modes" of operation for the Bayesian search. The first is the most general in which we evaluate the full likelihood function of Equation (18) and allow both the GW parameters, $boldsymbol {\lambda }$, and the noise model parameters, $boldsymbol {\phi }$ to vary simultaneously. In principle, this is the more desirable setup as it allows the uncertainty in our noise model to propagate into the measured GW parameters and also accounts for any correlations between the noise and GW parameters. This mode does require significantly more computational power as the number of search parameters in the MCMC is quite large. The total parameter space consists of eight GW parameters, N_psr pulsar distances, and 5 × N_psr noise parameters; this comes to 110 parameters for the full 17 pulsar array.

The second mode is when we fix the noise parameters to their maximum a posteriori values obtained from a previous single pulsar analysis. All previous GW searches for single sources have been performed in this manner (Jenet et al. 2004; Yardley et al. 2010; Babak & Sesana 2012; Ellis et al. 2012; Petiteau et al. 2013; Ellis 2013), which is justified if the noise model only contains white noise and the GW signal present in any single data set is weak. If the noise model contains only white noise, there will be little to no correlation between the GW parameters and the noise parameters, and if the signal is weak then it will not affect the single pulsar noise analysis. However, there is some evidence of red noise in our pulsars and because of the highly varying noise levels among pulsars, it is likely that a detectable source would be seen in the best timed pulsars individually. Therefore, this type of Bayesian analysis is not robust and could possibly lead to biased results; nonetheless, we will carry out this mode for comparison purposes in this study. Note that we will have the same problem with the $\mathcal {F}_p$ statistic. Possible methods to ameliorate this problem in fixed-noise searches are being explored and will be the subject of a future paper.

In a Bayesian sense, the detection problem is an application of model selection and the upper limit problem is an application of parameter estimation. The upper limit problem will be discussed in the next section. The model selection technique is very similar to that discussed in E13 and we refer the reader there for more details (see Appendix C.2 for more details on the implementation). If we want to perform model selection to claim a detection or compare different waveforms then we can make use of the Bayesian odds ratio between models "A" and "B"

$$\begin{equation} \mathcal {O}=\frac{p(d|\mathcal {H}_A)}{p(d|\mathcal {H}_B)}\frac{p(\mathcal {H}_A)}{p(\mathcal {H}_B)}, \end{equation} \tag{ 42 }$$

where d represents the measured data and $\mathcal {H}_{A}$ and $\mathcal {H}_{B}$ are two alternative hypotheses. The first factor is known as the Bayes Factor, which is our confidence in one model over the other based on the data (henceforth we will denote the Bayes factor as $\mathcal {B}$) and the second term is the prior odds ratio for models A and B, which describes our prior belief in both models. In this paper we consider only the Bayes factor, and assume the prior odds are even. (The choice of the prior odds will determine the false-alarm rate of a detection scheme based on the odds ratio (Vallisneri 2012).

In this work we wish to compare the signal model $\mathcal {H}_1$ parameterized by $\lbrace boldsymbol {\lambda }, boldsymbol {\phi }\rbrace$ to the null-hypothesis model $\mathcal {H}_0$ parameterized by $boldsymbol {\phi }$. When we allow the noise and GW parameters to vary simultaneously we will need to compute the evidence for both models, $\mathcal {H}_1$ and $\mathcal {H}_0$, separately. On the other hand, when we fix the noise parameters we can simply evaluate the likelihood-ratio in (20) and use this to compute the evidence. In this case, since the null-hypothesis model has no free parameters, the Bayes factor is simply given by the evidence computed using the likelihood ratio. In practice, to compute the evidence we make use of thermodynamic integration as detailed in E13 and Appendix C.2. The results of our Bayesian search and verification on injections will be discussed in the next section.

4.4.1. Priors

In a Bayesian analysis, especially when using parallel tempering and thermodynamic integration, it is very important to choose reasonable priors so that we are not exploring areas of parameter space that have been ruled out by previous experiments. We choose isotropic priors on all angular parameters and uniform priors in the log of the chirp mass with $\mathcal {M}\in [10^8,10^{10}]\, {M}_{\odot }$, log of the luminosity distance with d_L ∈ [1, 10⁴] Mpc, and log of the frequency of the GW with f_gw ∈ [6 × 10⁻⁹, 4 × 10⁻⁷] Hz. We impose an additional condition on the strain amplitude such that $h_0(\mathcal {M},d_L,f_{\rm gw}) \le h_{\rm ref}(f_{\rm gw}/f_0)^{2/3}$, where h_ref = 1 × 10⁻¹³ and f₀ = 10⁻⁸ Hz. This value is chosen so that the maximum strain is well above the level of detection. Essentially this is a cheap way to impose a correlated prior on chirp mass, luminosity distance, and GW frequency. The normalization is computed through Monte Carlo integration. For the pulsar distance prior we use the current electromagnetic (EM) measurements either from timing parallax or very long baseline interferometry corresponding to the best measured values taken from Verbiest et al. (2012; 10 pulsars) if available, otherwise, we use the values from the Australia National Telescope Facility (ATNF) pulsar catalog (Manchester et al. 2005)³⁷ which have distances based on dispersion measure and the NE2001 Galactic electron-density model (Cordes & Lazio 2002, 2003). For pulsars without parallax distances we assume a 20% uncertainty on the distance. Using this information, we write the distance prior as follows

$$\begin{equation} p(boldsymbol {L})=\prod _{\alpha =1}^{N_{\rm psr}}\frac{1}{\sqrt{2\pi \sigma _\alpha ^2}}\exp \left(-\frac{\left(L_{\alpha }-L^{\rm EM}_{\alpha }\right)^2}{2\sigma _{\alpha }^2} \right), \end{equation} \tag{ 43 }$$

where $L^{\rm EM}_{\alpha }$ is the best measured distance for the αth pulsar and σ_α is the 1σ uncertainty on that distance measurement. In principle it would be more correct to use a Gaussian prior for the parallax, which is proportional to L⁻¹. If the variance on the parallax is quite large then the corresponding prior on distance will differ significantly, namely it will have a long tail toward higher distances. However, for the pulsars used in this analysis, the distance uncertainty is small enough that the two prior distributions are effectively the same and we are safe in using a Gaussian prior on the pulsar distance itself; however, for future analyses we will move to Gaussian priors in L⁻¹.

For our noise parameters, we use priors that are uniform in the EFAC in the range [0.5, 5], uniform in the log of the EQUAD with EQUAD ∈ [10⁻⁹, 10⁻⁵] s, uniform in the log of the jitter value with the same range as the EQUAD, uniform in the log of the red noise amplitude with A_red ∈ [10⁻¹⁸, 10⁻¹¹], where the amplitude is in GW units, and uniform in the red noise spectral index with γ_red ∈ [1, 7]. We impose a further prior on the red noise such that the variance $\sigma _{\rm red}^2$ is less than the unweighted standard deviation of the pulsar timing residuals, where

$$\begin{eqnarray} \sigma _{\rm red}&=&\int _{1/T}^{\infty }{dfP}(f)\nonumber\\ &=& 2.05\frac{1}{\sqrt{\gamma _{\rm red}-1}}\left(\frac{A_{\rm red}}{10^{-15}}\right)\left(\frac{T}{1\, {\rm yr}}\right)^{\frac{\gamma _{\rm red}-1}{2}}\,\rm {ns}, \end{eqnarray} \tag{ 44 }$$

with T the total observation time and P(f) the power spectrum of the red noise. This prior essentially restricts the model from considering red-noise-dominated residuals, which is a very good approximation (Perrodin et al. 2013; J. A. Ellis et al. 2014, in preparation). This prior is chosen because it leads to much more computationally efficient runs by allowing us to run fewer high temperature chains in the Thermodynamic Integration scheme (see Appendix C.2 for more details). In principle this red noise prior is illegal in the sense that it uses the data (i.e., the variance of the residuals); however, this prior restricts access to an area of parameter space that is not supported by the likelihood function. Therefore, by omitting this area of parameter space the evidence calculation for each model, $\mathcal {H}_1$ and $\mathcal {H}_0$, will be biased slightly low. However, the likelihood function evaluated at this area of parameter space is essentially zero, and this slight bias will be negligible.

First, it is interesting to determine how much each pulsar in the 17 pulsar array will contribute to the overall S/N (signal-to-noise ratio) when a GW is present. In Figure 2 we plot the fraction ρ_α/ρ_total, where ρ_α is the single pulsar S/N, for each pulsar in the array. To compute this fraction we simulate 5000 S/N = 10 GW realizations (with parameters drawn from isotropic distributions in all angles and distributions uniform in the log of chirp mass and frequency) and calculate the single pulsar and total PTA S/N from Equation (22). The black (blue) points in the plot show the mean and standard deviation of the aforementioned ratio for each pulsar and the gray (green) curve is a simple naive scaling of $1/\sigma _{\alpha }^2$, where σ_α is the weighted rms of the αth pulsar's TOA uncertainties. It is obvious that J1713+0747 contributes more than 55% of the S/N on average, and PSRs J1909−3744, J0030+0451, and J0613−0200 contribute ∼10% on average. As we see from the gray (green) curve, this is very consistent with the overall scaling with the inverse of the variance of the noise; however, PSRs J0030+0451 and J0613−0200 carry a higher percentage because they are located opposite the bulk of other pulsars on the sky, and therefore will contribute more to the S/N for GWs coming from that side of the sky due to the antenna pattern response. This calculation does not mean that we advocate only timing the pulsars with the highest timing precision. Although many of the lower timing precision pulsars do not help with continuous GW detection or parameter estimation, they are essential for detection and parameter estimation of a stochastic GW background (see, e.g., Siemens et al. 2013).

Zoom In Zoom Out Reset image size

The fact that one pulsar dominates the total S/N means that it will be harder to make a confident GW detection as we require the same GW signal (with quadrupolar correlations) to be present in all pulsars. In other words, if the GW is only "seen" in one or two pulsars then it is hard to distinguish it from some other effect due to the pulsar timing model, ISM effects or some other systematic effect. This also implies the need to run a Bayesian analysis where both the noise and GW parameters are allowed to vary simultaneously. This does not mean that a continuous GW would not be a valid interpretation of a loud sinusoidal signal in one pulsar, only that statistically we do not have enough information to confidently claim a detection. Furthermore, if we did have a loud detectable signal, parameter estimation would be quite poor with the current NANOGrav PTA as there would be large degeneracies in the sky location (due to the small effective number of detector baselines), making sky localization and binary orientation estimates very poor. However, NANOGrav is currently timing 43 pulsars with microsecond or better precision. Also, new ultra-wideband receivers (DuPlain et al. 2008) have increased timing precision by a factor of ∼1.7 for many of the pulsars in this five-year data set. More pulsars and better timing precision could help ameliorate some of the limitations we have with the five-year data set.

Despite the potential limitations discussed above, we verify the efficacy of our pipeline by running both the frequentist and Bayesian pipeline on a synthetic data set with an injected GW source. To create the synthetic data set we first compute the residuals of our 17 NANOGrav pulsars using the tempo2 (Hobbs et al. 2006) package. Next we subtract these residuals from the site arrival times, thereby producing a new set of arrival times that match our timing model perfectly. To each set of idealized TOAs we then add a Gaussian noise process with the same characteristics as those measured in the real data, and a GW signal using the fully evolving signal model. We then use these new TOAs to produce a set of synthetic residuals. For this simulation we have chosen to inject a signal with S/N 10 and parameters $boldsymbol {\lambda } =\lbrace \theta =2.07, \varphi =5.4, f_{\rm gw}=4\times 10^{-8} \rm {\,Hz}, \mathcal {M}=5\times 10^{8}\, {M}_{\odot }, \mathnormal {d_L}=1.0\rm {\,Mpc}, \psi =0.78, \iota =0.26, \Phi _0=0.53\rbrace$.

The $\mathcal {F}_p$ statistic pipeline was run on this synthetic data set. Since we are treating this injection as if it were a true blind search, we must first run a single-pulsar noise analysis to determine the maximum a posteriori noise parameters; however, since a strong continuous GW and red noise will be covariant we have included a single frequency sinusoid as part of our noise analysis for each pulsar. This is implemented by simply adding a free amplitude and frequency parameter to the noise model discussed above. While this may appear to be special treatment for the injected signal, we have run the same noise model on the real data and find no evidence for any sinusoidal features. After we have obtained the maximum a posteriori noise parameters (not including the sinusoid parameters), we use these values to construct the noise covariance matrix for use in the $\mathcal {F}_p$ statistic as well as the fixed-noise Bayesian search. By performing the noise search with an included sinusoid but not including it in our noise covariance matrix in the subsequent GW analysis, we are sidestepping the problem of the GW being absorbed into red noise parameters.

We have carried out this analysis and the results are shown in Figure 3 where we plot $\mathcal {F}_p$ versus GW frequency when using the measured noise values (black) and the true injected noise values (gray). The vertical dashed line indicates the injected frequency and the horizontal dashed line represents our detection threshold corresponding to an FAP of 10⁻⁴. To compute the total FAP over all frequencies we make use of Equation (41) and choose N_f = 324, resulting in a total FAP of 1.6 × 10⁻⁸ which is a decisive detection. The number of independent frequencies is difficult to calculate when we are using many data sets with very irregular sampling. In this work we have chosen N_f = 324 as this corresponds to f_gw ∈ [1/T, 3.3 × 10⁻⁷] and Δf_gw = 10⁻⁹ Hz. The upper limit on frequency was chosen because our approximate observing cadence is 2.5 weeks⁻¹ and the frequency spacing was chosen by imposing the condition that the autocorrelation function of $\mathcal {F}_p$ when no signal is present drops to half of its maximum value at that frequency lag. This analysis shows us that we can indeed detect a continuous GW if it is present in our data by conducting a fully blind search; however, we also see that our results will not be conclusive as there are several frequencies at which the FAP is above our threshold value. From comparison with the true-noise case, we see that the uncertainty (and residual correlations between GW and noise parameters) in the noise parameters can lead to confusing results. This again is mostly due to the fact that our sensitivity is dominated by a small number of pulsars. Because of this, we caution against using a fixed-noise method to make final detection statements but instead advocate these methods as a first round in a suite of analyses.

Zoom In Zoom Out Reset image size

Both Bayesian pipelines (with and without varying noise parameters) were run on this synthetic data set. For both runs we have used PTMCMC and thermodynamic integration as discussed in Appendix C.2. Due to the large parameter spaces when using the full GW and noise model, we have chosen to use only the pulsars that contribute more than 1% to the injected S/N, resulting in six pulsars, J1713+0747, J1909−3744, B1855+09, J0030+0451, J0613−0200, and J1012+5307. Here we use the same noise parameters as mentioned above for our fixed-noise search. Even though these estimates are different from the true noise parameters, we nonetheless achieve a log-Bayes factor of 27.4 for the fixed-noise search (a log-Bayes factor greater than five is considered decisive evidence). However, as we mentioned earlier, we should not totally trust this level of evidence as it does not fully incorporate our uncertainty in the noise model. When we run an analysis where we allow the noise and GW parameters to vary simultaneously, we only achieve a log-Bayes factor of 5.35. While still decisive, this search is much less sensitive to the GW; nonetheless, this search is the most robust and will be the real test as to whether or not one can trust a real GW detection candidate. Of course these results could change depending on the noise realization or GW parameter combinations. A more detailed study of this is warranted but beyond the scope of this paper. Nonetheless, the large spread of overall noise levels in modern PTAs will most likely make the confident detection of a continuous GW very difficult.

First we will discuss the results of the $\mathcal {F}_p$ statistic search on the real 17 pulsar NANOGrav data. To carry out the analysis we have computed $\mathcal {F}_p$ for many frequencies with f_gw ∈ [1/T, 3.3 × 10⁻⁷] Hz. These frequencies were chosen based on the fact that the approximate cadence is 2.5 week⁻¹. The results of this search are shown in Figure 4, where the solid black line is the value $\mathcal {F}_p$ at each frequency; and the dotted, dash-dotted, and dashed lines are the value of $\mathcal {F}_p$ corresponding to a 1.0%, 0.5%, and 0.1% FAP, respectively, where these values are calculated from Equation (40). Furthermore if we maximize $\mathcal {F}_p$ over frequencies then the total FAP, accounting for the trials factor N_f, is very nearly 1, indicating that we should accept the null hypothesis (no visible GW signal) with very high confidence.

Zoom In Zoom Out Reset image size

We will now briefly discuss the results of both Bayesian searches. To carry out this analysis we have run our PTMCMC and computed the Bayes factors for each case. In the first case we allow the noise parameters and GW parameters to vary and explicitly compute the Bayesian evidence via thermodynamic integration for a model with a GW and noise and a model with just noise. In the second case, we fix the noise parameters to the maximum a posteriori obtained from single pulsar analyses and only allow the GW parameters to vary. As mentioned above, the second case is not reliable since there are likely to be correlations between the GW and noise parameters; however, we give the results of both searches for completeness. As above, in the case of a true continuous GW signal we can get very different results from a fixed-noise versus a varying noise search. However, in our case the log-Bayes factor for searches with and without varying noise parameters is −0.55 and −0.1, respectively, both indicating that there is no evidence for a continuous GW and a model consisting of noise is preferred. We further note that this is completely consistent with our frequentist analysis.

In this section we will outline the procedures used to compute both the frequentist and Bayesian upper limits on the strain amplitude, h₀. First we wish to state that an x% upper limit on the strain amplitude does not mean that we would have detected a signal with that amplitude with 95% confidence, it simply means that the true value of the amplitude is less than the upper limit with 95% probability (this probability measures the frequency of measuring that value of the amplitude in the frequentist case and the degree of belief in the true value of the amplitude in the Bayesian case). In the following sections we will discuss the mathematics of upper limit computation in the frequentist and Bayesian frameworks, and then we will lay out our computational procedure.

5.3.1. Frequentist Approach

From a frequentist viewpoint, the data are random while the signal parameters are fixed but unknown (i.e., we construct probability distributions for the data, or rather some function of the data, given a set of signal parameters), whereas in the Bayesian framework the data are fixed and the signal parameters are uncertain (i.e., we construct probability distributions of the signal parameters given a data set). From the above statement it then follows that frequentist upper limits are derived from integrating the probability distribution of some statistic of the data (the $\mathcal {F}_p$ statistic in this case) at a fixed value of the parameter of interest, and Bayesian upper limits are derived from integrating the probability distribution of the parameter of interest for the given data set.

More formally, the probability distribution of the $\mathcal {F}_p$ statistic given a value of the strain amplitude h₀ is

$$\begin{equation} p(\mathcal {F}_p|h_0)=\int p(\mathcal {F}_p|h_0,\tilde{\lambda },\mathbf {n})p(\tilde{\lambda })p(\mathbf {n})\,d\tilde{\lambda }\,d\mathbf {n}, \end{equation} \tag{ 45 }$$

where $\tilde{\lambda }=\lbrace \theta,\varphi, f_{\rm gw}, \mathcal {M}, \iota, \psi, \Phi _0\rbrace$ is a reduced parameter space vector, $p(\tilde{\lambda })$ is the sampling distribution of $\tilde{\lambda }$ (these sampling distributions are identical to the prior probability distributions in the Bayesian case), n is a noise time series drawn from the distribution

$$\begin{equation} p(\mathbf {n})=\frac{1}{\sqrt{\det {2\pi C}}}\exp {\left(-\frac{1}{2}\mathbf {n}^TC^{-1}\mathbf {n} \right)}, \end{equation} \tag{ 46 }$$

with C the covariance matrix of the noise in the pulsar timing residual time series, and $p(\mathcal {F}_p|h_0,\tilde{\lambda },\mathbf {n})$ is the probability distribution function for the $\mathcal {F}_p$ statistic for given values of h₀ and $\tilde{\lambda }$ and a given noise realization n. An upper limit on h₀ at confidence level α is then computed as

$$\begin{eqnarray} \alpha &=& \int _{\mathcal {F}_{p,0}}^{\infty } p(\mathcal {F}_p|h_0)\,d\mathcal {F}_p\nonumber\\ &=& \left \langle \frac{1}{N} \sum _{i=1}^{N} \left\lbrace \begin{array}{@{}l@{\quad }l@{}} 1 & {\rm if\,\,} \mathcal {F}_{p,i} \ge \mathcal {F}_{p,0}\\ 0 & {\rm otherwise} \end{array}\right \rbrace \right \rangle, \end{eqnarray} \tag{ 47 }$$

where the N observables $\mathcal {F}_{p,i}$ are drawn from the "signal distribution," $p(\mathcal {F}_p|h_0)$, and the average, 〈 · 〉, is over that distribution. In other words, we integrate the probability distribution of the $\mathcal {F}_p$ statistic over the so called "signal space" (i.e., from the measured value $\mathcal {F}_{p,0}$ to infinity), that is, we count the number of signal realizations that gives an $\mathcal {F}_p$ statistic value larger than the one measured in the actual data set. This integral can take on any value α ∈ [0, 1] for a given h₀; therefore, the integral must be repeated with different values of h₀ until α = 0.95 for a 95% upper limit.

In practice, we carry out the following computational procedure.

• 1.  
  Measure the value $\mathcal {F}_{p,0}$ from the real 17 pulsar NANOGrav data set as described in Section 4.3.
• 2.  
  Simulate a synthetic noise vector n = Lw for each pulsar, where L is the Cholesky decomposition of the noise covariance matrix C, and w is a unit variance, zero-mean, vector.
• 3.  
  Choose strain amplitude, h₀ and construct a GW waveform $\mathbf {s}(t,h_0,\tilde{\lambda })$ for each pulsar where the parameters $\tilde{\lambda }$ are drawn from the distribution $p(\tilde{\lambda })$.
• 4.  
  Construct a new set of residuals for each pulsar $\delta t_{\rm sim} = R(\mathbf {n} + \mathbf {s}(t, h_0,\tilde{\lambda }))$, where R is the so-called fitting projection matrix introduced in Demorest et al. (2013) and Ellis et al. (2013).³⁸
• 5.  
  Now measure the value $\mathcal {F}_{p,i}$ for the simulated data set.
• 6.  
  Repeat steps 2–5 10,000 times and measure the number of realizations that result in $\mathcal {F}_{p,i} > \mathcal {F}_{p,0}$.
• 7.  
  Repeat steps 2–6 with different values of h₀ until 95% of simulations result in $\mathcal {F}_{p,i} > \mathcal {F}_{p,0}$.

In the remainder of the paper we will choose to compute upper limits on the strain amplitude as a function of GW frequency or GW sky location at a fixed GW frequency. To facilitate such upper limits we simply fix the parameters (either GW frequency or sky location) when simulating waveforms in step 3.

5.3.2. Bayesian Approach

As mentioned above, in the Bayesian framework we do not rely on simulations as we treat the data as fixed and integrate the posterior PDF of the parameter of interest to compute the upper limit. In principle, a Bayesian upper limit is much more simple and intuitive than a frequentist upper limit. To compute a Bayesian upper limit we compute an integral that is analogous to Equation (47)

$$\begin{eqnarray} \alpha &=& \int _0^{h_{\rm up}}dh_0 \int d\tilde{\lambda }\, dboldsymbol {\phi }\, p(\delta t|h_0,\tilde{\lambda },boldsymbol {\phi })p(h_0)p(\tilde{\lambda })p(boldsymbol {\phi })\nonumber\\ &=& \int _0^{h_{\rm up}}dh_0\,p(\delta t|h_0)p(h_0), \end{eqnarray} \tag{ 48 }$$

where $p(\delta t|h_0,\tilde{\lambda },boldsymbol {\phi })$ is the likelihood function; and p(h₀), $p(\tilde{\lambda })$, $p(boldsymbol {\phi })$ are the prior probability distributions on h₀, $\tilde{\lambda }$, and $boldsymbol {\phi }$, respectively, where $boldsymbol {\phi }$ denotes the noise model parameters. In words, we simply integrate the marginalized posterior distribution of h₀ until the desired credible region corresponding to a probability of α is reached at h = h_up. As in the frequentist case, we want upper limits on the strain amplitude as a function of GW frequency or sky location. In this case we simply fix the parameters and then marginalize over the others. In practice, to compute the Bayesian upper limits we carry out a separate MCMC run for fixed values of frequency and/or fixed sky locations and then compute the 95% upper limit for each. The choice of prior on h₀ is very important and can lead to very different upper limits. Such a detailed analysis of priors is beyond the scope of this work but will be addressed in a future paper. In principle, our prior distribution should come from population synthesis models (Sesana 2013); however, since we wish our upper limits to be informed by our data and not dominated by our prior distribution we use a very conservative³⁹ prior that is uniform in h₀ with h₀ ∈ [0, 10⁻¹¹].

5.3.3. Sky-averaged Strain Upper Limits

In Figure 5 we report the 95% upper limits on the strain amplitude, h₀, as a function of GW frequency computed using the methods described above for the frequentist and Bayesian pipelines. The gray (red), thick black (blue), and thin black (purple) curves are the 95% upper limits on strain amplitude computed using the $\mathcal {F}_p$ statistic, Bayesian method with fixed-noise values, and Bayesian method with varying noise values, respectively. There are several features in the plot that require explanation. First, the decrease in sensitivity at f_gw = 1 yr⁻¹ and f_gw = 2 yr⁻¹ is due to the sky position and parallax fitting in the timing model, respectively. The upward trend at lower frequencies is due to the quadratic spin-down model fit. The noisiness of the frequentist upper limit is due to the fact that the $\mathcal {F}_p$ statistic distribution at higher frequencies is indeed quite noisy when computed using the real data, and since our upper limits compare the value measured in real data to values measured in simulated data, this noisiness is to be expected.

Zoom In Zoom Out Reset image size

If we compare our results to those of Yardley et al. (2010), we see that the upper limits using the five-year NANOGrav data sets are a factor of two to three times more constraining. The main reason for this improvement is the higher timing precision of the NANOGrav data set as compared to the older PPTA data sets (Verbiest et al. 2009). Although the procedures for setting frequentist and Bayesian upper limits are quite different, our results are very similar. In part, this is due to the fact that we have used a uniform prior on the strain amplitude, h₀, making our Bayesian analysis very similar to a pure likelihood analysis. Since the $\mathcal {F}_p$ statistic is just the likelihood (ratio) maximized over amplitudes, we would expect a likelihood analysis to give similar results. Note that the Bayesian upper limits when varying the noise parameters are somewhat less constraining than the fixed-noise case. This is to be expected since at lower frequencies the GW amplitude and red noise amplitude are somewhat correlated and at higher frequencies the GW amplitude and jitter parameter are somewhat correlated. Both correlations will result in slightly worse upper limits on the GW amplitude when allowing the noise parameters to vary.

In Figure 5, the dashed curves indicate lines of constant chirp mass for a source with a distance to the Virgo cluster of 16.5 Mpc and chirp mass of 10⁹ and 10¹⁰, respectively, and the gray (green) squares are the strain amplitude of the loudest GW events in 1000 Monte Carlo realizations using an optimistic phenomenological model of Sesana (2013). The model used here produces a stochastic GW background with dimensionless strain amplitude of ∼2 × 10⁻¹⁵, just below the current upper limits presented in Shannon et al. (2013). Astrophysically, these upper limits tell us that we can essentially rule out any source with $\mathcal {M}\ge 10^{10}\, {M}_{\odot }$ at the distance to the Virgo cluster (16.5 Mpc); however, our horizon distance falls just short of the Virgo cluster for sources with $\mathcal {M}\le 10^9\, {M}_{\odot }$. This limit is not particularly constraining for the Virgo cluster, which generally does not contain galaxies with such large implied masses; however, it is interesting to assess this limit against other nearby massive galaxies (see, Lommen & Backer 2001). Moreover, at our most sensitive frequency of ∼10 nHz, we can convert our upper limit on the GW strain to a lower limit on the luminosity distance of d_L ⩽ 425 Mpc for sources with $\mathcal {M}=10^{10}\, {M}_{\odot }$.

Finally, we see that all sources from Monte Carlo realizations (gray (green) squares) have strain amplitudes below our upper limits indicating that it is very unlikely that we will see a resolvable source at the current sensitivity (consistent with our search results). It is important to note, however, that these strain amplitude upper limits are averaged over sky location and inclination angle (either through marginalization in the Bayesian case, or from Monte Carlo sampling in the frequentist case), both of which play a large part in the overall amplitude of the signal. Therefore, these results have the caveat that they make statements about the average sensitivity to such GW sources; however, it is still unlikely (i.e., probability of detection ≲ 50%) that we could detect even the loudest optimally oriented source shown in Figure 5. For face-on systems (i.e., ι = π/2) and sky location near the best timed pulsars, the overall amplitude of the GW can be ∼five times larger than the averaged strain amplitudes reported here.

5.3.4. Angular Upper Limits

In Figures 6 and 7 we report the 95% lower limit on the luminosity distance as a function of sky location computed using the frequentist and Bayesian techniques, respectively. We have chosen to present our results in terms of the luminosity distance instead of the strain amplitude as it is a true physical parameter and it gives a more intuitive feel as to what the data can constrain. To compute this lower limit we carry out the same procedure as above but we fix the frequency to f_gw = 10⁻⁸ Hz and compute an upper limit on the strain amplitude as a function of sky location; we can then use Equation (17) to convert an upper limit on strain amplitude into a lower limit on luminosity distance. The values in the color bar are calculated assuming a chirp mass of $\mathcal {M}=10^{9}\, {M}_{\odot }$ and a frequency of f_gw = 10⁻⁸ Hz but this can be scaled to determine the minimum luminosity distance for any chirp mass value and GW frequency. In Figures 6 and 7 the white diamonds represent the locations of the 17 NANOGrav pulsars used in the analysis and the black (white) stars are the sky locations of potential GW hotspots (Simon et al. 2014) and possible GW source candidates (Valtonen et al. 2008; Iguchi et al. 2010; Ju et al. 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We will now discuss the features of this sky-dependent upper limit computed using the frequentist $\mathcal {F}_p$ statistic. Firstly, we notice that the overall distribution is quite similar to the antenna pattern response (i.e., 1 + cos μ) as is to be expected in the case of no detection. Due to this, we are most sensitive (larger lower limit on luminosity distance) at sky locations near the best timed pulsars (i.e., J1713+0747, B1855+09, J1909−3744) and least sensitive in the opposite direction. More quantitatively, we note that in the most sensitive areas of the sky we can constrain the luminosity distance d_L ≳ 47 Mpc for $\mathcal {M}=10^{9} \, {M}_{\odot }$. Furthermore, it is possible to constrain the luminosity distance d_L ≳ ∼ 2 Gpc in the most sensitive sky locations if we consider 10¹⁰ M_☉ chirp mass sources. It should be noted that the Bayesian fixed-noise search gives nearly identical results to the fixed-noise frequentist search.

We now move to the sky-dependent upper limit computed using the full Bayesian technique where the GW and noise parameters are varied simultaneously. The first observation that we make is that the overall scale is about a factor of two lower than the fixed-noise frequentist or Bayesian upper limit. At first this may be surprising given the general agreement of the sky-averaged upper limits of Figure 5; however, full Bayesian sky-dependent upper limits exacerbate the problem of relatively few pulsars contributing to the overall PTA sensitivity as shown in Figure 2. Another difference in this upper limit, as opposed to the frequentist upper limit, is that it does not quite match the expected antenna pattern response function. These differences are due to the fact that we are simultaneously varying the GW and noise parameters, and when only one or a few pulsars contribute to the PTA sensitivity, there is a degeneracy between intrinsic red-noise processes in the pulsar and a common GW among all pulsars. In other words, it is very difficult to distinguish between a low-frequency continuous GW and a red noise process if only a small number of pulsars have sufficiently low noise levels to resolve the GW.

Because Bayesian upper limits marginalize or integrate over all parameters except the amplitude, the correlations between the GW and the red noise amplitude will broaden the 1D PDF of the amplitude and thus will result in larger upper limits as opposed to the fixed-noise case. As is clear from Figure 7, the aforementioned effect is very strong for GW sky locations near our best timed pulsars. For example, we are not most sensitive to GWs around the sky location of PSR J1713+0747 because this pulsar contributes a very large percentage of the overall S/N of the GW in this case and thus results in a very large correlation between the GW and red noise amplitudes.

Since, at the moment, we have no way of measuring the noise properties of the pulsars independently of any GWs that may be present in the data, to perform a completely robust upper limit or search we must allow both to vary simultaneously. Given this reality, we must view any fixed-noise results with the caveat that they assume that the noise parameters are measured perfectly and are independent of any GWs in the data.

Unfortunately, many of the GW hotspots and potential SMBHB sources are located at insensitive sky locations, for both frequentist and Bayesian analyses, where our lower limit on distance only allows us to constrain 10¹⁰ M_☉ sources. This fact is a great argument for aggressive pulsar search campaigns and the addition of new pulsars to the PTA at sky locations that are currently insensitive (Burt et al. 2011).

5.3.5. Constraints on the SMBHB Coalescence Rate

A non-detection of continuous GW, as we have presented here, allows us to compute an upper limit on the rate of SMBH coalescences using methods presented in Wen et al. (2011). Since we have made no detections, we assume Poisson statistics for the probability of an event (i.e., a detectable signal) occurring, that is, the probability of no events is e^−〈N〉, where 〈N〉 is the expected number of events. We use this probability distribution function to place a 95% upper limit on the expected number of events such that exp (− N₉₅) = 0.05, telling us that 〈N〉 ⩽ 〈N₉₅〉 = 3. Therefore, if the expected number of events were greater than 3, at least one source would have been detected with 95% probability. Now, following Wen et al. (2011), the expected number of events is

$$\begin{eqnarray} \langle N \rangle &=& \int \frac{d^2R}{d\log _{10}(1+z)d\log _{10}(\mathcal {M}_r)} \left(\frac{df}{dt} \right)^{-1}\nonumber\\ &&\times P_d(\mathcal {M}_r,z,f)\, d\log _{10}(1+z)\, d\log _{10}(\mathcal {M}_r)\,df,\nonumber\\ && \end{eqnarray} \tag{ 49 }$$

where $P_d(\mathcal {M}_r,z,f)$ is the probability of detecting an SMBHB with chirp mass $\mathcal {M} = \mathcal {M}_r(1+z)$, redshift z, and observed GW frequency f. Following the derivation in Wen et al. (2011) and making the assumption that the differential coalescence rate does not vary significantly over the range $\Delta \log _{10}{\mathcal {M}_r}=1$ and Δlog₁₀(1 + z) = 0.2, it is possible to show that

$$\begin{equation} \frac{d^2R}{d\log _{10}(1+z)d\log _{10}(\mathcal {M}_r)} < \frac{15}{\int \left(\frac{df}{dt} \right)^{-1}P_d(\mathcal {M}_r,z,f)\, df}. \end{equation} \tag{ 50 }$$

In order to compute the detection probability P_d, we make use of the $\mathcal {F}_p$ statistic. We use the same method that we have used for the upper limits, except now we compare the value of $\mathcal {F}_p$ computed using simulated data with injections to a specified threshold based on an FAP of 10⁻⁴. We use 10,000 realizations at each value of z and f. After the probability of detection is computed, we numerically integrate the above expression to obtain a limit on the differential coalescence rate. It should be noted that we will be able to place more meaningful constraints on the coalescence rate using upper limits on the amplitude of a stochastic background of SMBHBs; however, this is beyond the scope of this work and will be addressed in a future paper. In Figure 8 we plot our constraints on the differential coalescence rate as a function of redshift. Since we have made the assumption that this differential coalescence rate does not vary significantly over an order of magnitude in chirp mass, the results presented here are for the $\mathcal {M}_r=10^{10}\, {M}_{\odot }$ case. We are unable to place meaningful constraints on less massive systems. The light gray (red) shaded area is constructed using the model presented in Jaffe & Backer (2003) along with measurements from the Sloan Digital Sky Survey (Wen et al. 2009). The medium gray (blue) shaded area is constructed by considering the different galaxy merger rates based on observations (Sesana 2013) along with the most recent MBH–sigma relation from McConnell & Ma (2013). The dashed line comes from an a posteriori implementation of the McConnell & Ma (2013) MBH–sigma relation into the semi-analytic model of Guo et al. (2011) assuming accretion onto both SMBHs before merger. The black (green) shaded region is constructed by using the observed evolution of the galaxy mass function combined with the MBH–M-stars relation from McConnell & Ma (2013) to calibrate an analytical model for evolving the mass function via mergers (McWilliams et al. 2012). The figure shows that the coalescence rate for MBHs of ∼10¹⁰ M_☉ is poorly constrained. This is mostly because of the steepness of the galaxy mass function at such high masses: a small change in the slope results in a large variation in the sparse population of 10¹⁰ M_☉ black holes. The intrinsic scattering (e.g., Gültekin et al. 2009) and poor knowledge of the behavior of the MBH–galaxy relations at the high mass end (e.g., Lauer et al. 2007) add further uncertainties, making the coalescence rate estimate problematic. As is clear from the figure, we are unable to place any constraints on the physical models mentioned above; however, as our GW sensitivity improves with time, we will begin to place meaningful constraints on physical models pertaining to the coalescence rate of SMBHBs.

Zoom In Zoom Out Reset image size

Predicting the future sensitivity of PTAs to continuous GWs is quite difficult and depends on a number of poorly constrained factors that make up the entire noise budget for each pulsar (Cordes & Shannon 2010). A detailed study of these effects on the future performance of a PTAs sensitivity to continuous GWs is beyond the scope of this work. Here we simply derive a rough scaling law for the S/N of a continuous GW measured by a PTA and make a few statements about expected future performance. The square of the S/N is defined to be⁴⁰

$$\begin{equation} \rho ^2 = \sum _{\alpha }(s_{\alpha }|s_{\alpha }) = \sum _{\alpha }2{\rm Re\,}\int _{-\infty } ^{\infty }df \frac{|\tilde{s}_{\alpha }(f)|^2}{S^{n}_{\alpha }(f)}, \end{equation} \tag{ 51 }$$

where $\tilde{s}_{\alpha }(f)$ is the Fourier transform of the GW induced timing residuals and $S_{\alpha }^n(f)$ is the power spectral density of the noise for the αth pulsar. As we mentioned above, the frequency of the GW will not vary over the observation time and the waveform is approximated by a sine wave at a single frequency. We will further assume that the pulsar term is at the same frequency for this scaling law computation. The S/N then becomes

$$\begin{eqnarray} \rho ^2 &\approx& \sum _{\alpha }2{\rm Re\,}\int _{-\infty }^{\infty } df \frac{A(f)^{2}a_{\alpha }^{2}\delta (f-f^{\prime })^{2}}{S^{\rm red}_{\alpha }(f)+\sigma _{\alpha }^{2}/c_{\alpha }}\nonumber\\ &\approx& A(f^{\prime })^{2}\sum _{\alpha } \frac{T_{\alpha }}{2}\frac{a_{\alpha }^{2}}{S^{\rm red}_{\alpha }(f^{\prime })+\sigma _{\alpha }^{2}/c_{\alpha }}, \end{eqnarray} \tag{ 52 }$$

where δ(f − f') is the Dirac delta function, A(f') is the pulsar independent amplitude of the GW, a_α is a geometric factor that depends on the antenna pattern functions for each pulsar, c_α is the observing cadence for each pulsar, and $S^{\rm red}_{\alpha }(f^{\prime })$ and σ_α are the red noise power spectral density at f = f' and white noise rms, respectively. Our sensitivity to continuous GWs is proportional to the S/N, thus this expression for the S/N can serve as a proxy for how our upper limits and sensitivity will improve with various quantities. It is interesting to examine this scaling law in the white noise and red noise dominated regime

$$\begin{equation} \rho _{\rm white}\propto \left(\sum _{\alpha } \frac{T_{\alpha }a_{\alpha }^{2}c_{\alpha }}{\sigma _{\alpha }^{2}}\right)^{1/2} \end{equation} \tag{ 53 }$$

$$\begin{equation} \rho _{\rm red} \propto \left(\sum _{\alpha } \frac{T_{\alpha }a_{\alpha }^{2}}{S^{\rm red}_{\alpha }(f^{\prime })}\right)^{1/2}. \end{equation} \tag{ 54 }$$

The above scaling laws tell us that many pulsars distributed across the sky with high timing precision, high observing cadences, and low red noise levels observed over a long baseline will result in the best possible sensitivity to continuous GWs. New pulsar timing backends at Arecibo and the GBT (DuPlain et al. 2008) give roughly a factor of two higher timing precision for many pulsars, which will translate into an expected upper limit on the amplitude of continuous GWs that is a factor of two more constraining. As we acquire more data on our currently timed and newly discovered pulsars we will gain more sensitivity and will be able to probe to lower frequencies. Access to IPTA data would essentially serve to increase the observing cadence, and thus our sensitivity, since we would have complementary data from many different observatories measured at different times. Furthermore, current pulsar search campaigns (Lynch & Green Bank North Celestial Cap Survey Collaborations 2013) are discovering new MSPs in our least sensitive sky locations (see Figures 6 and 7), which will dramatically increase our overall sky coverage and will allow for better distinction between GW and noise models. Finally, advanced detectors, such as the Square Kilometer Array (SKA; Lazio 2013) are expected to time tens of pulsars at or below the 100 nanosecond level, which will likely solve many of the current problems that we face with poor angular sensitivity and inability to distinguish between single GW source and intrinsic pulsar noise.

In the red noise dominated regime, the cadence of observations and the overall white noise level is negligible and we essentially only gain sensitivity through the addition of new pulsars and continued timing. Here we note again that there is very little evidence for red noise in the five-year NANOGrav data set but this may change in the future as we become more sensitive to the stochastic GW background which would induce a common red noise signal in all pulsars. Of course, red noise from a stochastic background of GWs or from intrinsic pulsar spin noise (Shannon & Cordes 2010) is likely to have a steep spectrum and this red-noise dominated regime would only apply at the lowest frequencies.

For many pulsars in the five-year data set we are in the white noise dominated regime and since the S/N consists of the sum of the inverses of the white noise rms values, we see that only the best timed pulsars will contribute (as we have seen throughout this paper) and one may argue that we should focus all observing time on the best pulsars. However, as was shown in Siemens et al. (2013), our sensitivity to the stochastic background has a significantly weaker dependence on observing cadence and white noise rms but has a linear dependence on the number of pulsars in the array. Thus, it is difficult to realistically optimize a PTA for both continuous and stochastic GW sources.

In this paper we have performed various searches for continuous GWs from non-spinning SMBHBs in circular orbits using both frequentist and Bayesian techniques. Specifically, we have run a fixed-noise frequentist and Bayesian pipeline, as well as a varying noise Bayesian pipeline. In the absence of any detections we have placed upper limits on the strain amplitude of continuous GWs as a function of GW frequency. We have also computed a lower limit on the distance to such SMBHBs as a function of sky location, as well as placing constraints on the differential coalescence rate of such SMBHBs. Our sky-averaged upper limits on strain amplitude as a function of frequency are a factor of ∼three times more constraining than the previously published upper limits (Yardley et al. 2010) and we see good agreement between all three data analysis methods. Although improving, our limits still lie well above the amplitudes of individual sources produced from several realizations of an optimistic SMBHB population. We have shown that with good estimates of the intrinsic noise we can rule out any sources with luminosity distance <425 Mpc and a chirp mass of ∼10¹⁰ M_☉. Furthermore, this limit becomes about a factor of four more constraining for GW sources that are near our most precisely timed pulsars in the sky. Unfortunately we are not yet able to place any constraints on predictions for the coalescence rate of SMBHBs obtained from both theory and observations.

Throughout the paper we have made several statements about what is needed for completely robust data analysis techniques and what will be required from future PTAs in order to secure a confident detection of a continuous GW. These statements can be summarized as follows.

• 1.  
  Currently we have no way to confidently separate intrinsic noise in the residuals from any GW that may be present. Therefore, it is necessary to include both noise and GW parameters in any data analysis pipeline that aims to be truly robust. This is not to say that fixed-noise methods should not be used; instead we advocate a hierarchical approach where the faster fixed-noise methods are used as a first pass and then followed up with a full GW plus noise search. Lastly, a signal with more information, such as that from an eccentric system, could help break this degeneracy between signal and noise models and will be the subject of a future paper.
• 2.  
  Even with simultaneous noise and GW characterization, unless we have several well-timed pulsars (with very similar timing precision on all) with decent sky coverage, a confident detection of a continuous GW is unlikely even if the signal is loud.

In summary, we note that while not as likely as a detection of a stochastic GW background, with continually improving timing precision, the addition of new pulsars to PTAs, and improved data analysis techniques, prospects are good for obtaining astrophysically constraining GW limits, or possibly even a detection of a continuous GW, over the next decade.

In order to facilitate good mixing of the MCMC chains, especially in large parameter spaces, it is very important to have good jump proposals. In our implementation of the PTMCMC algorithm we use a jump proposal that is composed of a randomized cycle of sub-proposals. Here we will briefly outline the different jump proposals used in the cycle.

C.1.1. Correlated Jumps

As described in E13, we use an Adaptive Metropolis (AM) scheme to make correlated jump proposals. In essence, this jump uses the previous points in the chain's history to construct a sample covariance matrix, C_n (Haario et al. 2001), which is then updated throughout the run. In practice we update this covariance matrix every 1000 iterations. Primarily the sample covariance matrix is multiplied by a scale parameter s_d = 2.4²/n_dim, which gives an optimal 25% acceptance rate in the case of Gaussian posterior distributions; however, we will occasionally make small jumps (scale by 0.01) or large jumps (scale by 10). This jump is used in ∼20% of our total jump cycle.

In large parameter spaces, such as those we encounter when modeling the GW and noise simultaneously, the above method can result in very low acceptance and thus, slow convergence. Haario et al. (2005) introduce the Single Componentwise Adaptive Metropolis (SCAM) algorithm in which only one uncorrelated variable is updated in the jump proposal. If the variables are completely uncorrelated, then this method is identical to using the AM algorithm but only updating one parameter. However, if the parameters are correlated, as they are in our problem, we can define a set of uncorrelated parameters

$$\begin{equation} y = U^Tx, \end{equation} \tag{ C1 }$$

where x is our original vector of parameters and U is defined by the eigenvalue decomposition C_n = USU^T, where U is a unitary matrix and S is a diagonal matrix. It is then easy to see that the covariance matrix of y, averaged over many steps in the chain, is

$$\begin{equation} \langle y y^T\rangle = U^T\langle x x^T\rangle U = U^T U S U^TU = S. \end{equation} \tag{ C2 }$$

Since $S={\rm diag}\lbrace \sigma _s^2\rbrace$ is a diagonal matrix, each y represents an uncorrelated parameter. Therefore, we choose an uncorrelated parameter at random and propose the jump

$$\begin{equation} y_{i+1}^j = y_i^j + 2.4\mathcal {N}\left(0,\sigma _s^j\right), \end{equation} \tag{ C3 }$$

where N(0, σ_s) is a zero mean Gaussian deviate with variance $\sigma _s^2$, i is the iteration number, and j is the parameter number. We can then relate the jump in the uncorrelated parameters back to a jump in the correlated parameters

$$\begin{equation} x_{i+1} = Uy_{i+1}. \end{equation} \tag{ C4 }$$

If U is not the identity matrix (i.e., the parameters, x, are correlated) then this means that we will jump in combinations of correlated physical parameters even though we only jump in one uncorrelated component at a time. We have found that jumps of this kind greatly improve mixing when running with a large number of search parameters (e.g., >100). This jump is used in ∼40% of our total jump cycle.

We also employ a third type of correlated jump proposal known as differential evolution (DE; Braak 2006). Differential evolution is a simple genetic algorithm that also makes use of the previous history of samples in the chain. A differential evolution jump can be constructed as follows. First choose, at random, two previous iterations of the chain. Denote the parameter vector at those two new points as x_m and x_n. The DE jump is then

$$\begin{equation} x_{i+1} = x_i + s_{DE}(x_m-x_n), \end{equation} \tag{ C5 }$$

where s_DE is a scale factor which we choose to be s_DE = 2.4²/n_dim and s_DE = 1, each with 50% probability. The first scale factor here is identical to that used in the AM jumps and the second is known as a "mode jump," that is, if x_m and x_n are located at two different modes of the posterior distribution, then the mode jump will result in a jump that stays on the same mode as x_i or jumps to the other mode. For this reason, DE jumps are usually employed if there are strong multimodal structures in the posterior PDF, which we may expect in the case of a weak continuous GW. Also, since we are drawing points from the posterior, then these jumps will also "learn" about any correlations among parameters and will be taken into account in the jump proposal. This jump is used in ∼20% of our total jump cycle.

Finally, we use a special jump proposal for the chirp mass and luminosity distance. This jump makes use of the fact that there is likely to be a large correlation between chirp mass and luminosity distance for nearly non-evolving sources. If the frequency evolution is negligible then the waveform only has information about the combination $\mathcal {M}^{5/3}/d_L$, with a very weak mass dependence in the frequency derivative term. For this jump, we draw a random luminosity distance value from the prior and then resolve for what value of chirp mass will keep the combination $\mathcal {M}^{5/3}/d_L$ constant. This jump is used in ∼5% of our total jump cycle.

C.1.2. Uncorrelated Jumps

C.1.3. Auxiliary Pulsar Mode Jump

Here we will describe our parallel tempering and thermodynamic integration techniques used to calculate the Bayesian evidence. We want our algorithm to then quickly locate the global maxima in the parameter space. To accomplish this in a way that satisfies detailed balance we make use of parallel tempering. This technique involves different chains exploring the parameter space simultaneously, each with a different target distribution

$$\begin{equation} p(\Theta |d,\beta)=p(\Theta)p(d|\Theta)^{\beta }, \end{equation} \tag{ C9 }$$

where β ⩽ 1 is the inverse "temperature" and Θ parameterizes the model. This will essentially flatten out the likelihood surface allowing the chains to more freely explore the entire prior volume. The "hot" chains will inform the "colder" chains and vice versa by proposing parameter swaps between different temperatures. A parameter swap between the ith and jth temperature is accepted with probability α = min (1, H), where the multi-temperature Hastings ratio is

$$\begin{equation} H_{i\rightarrow j}=\frac{p(d|\Theta _i,\beta _j)p(d|\Theta _j,\beta _i)}{p(d|\Theta _i,\beta _i)p(d|\Theta _j,\beta _j)}. \end{equation} \tag{ C10 }$$

By swapping parameter states between different temperatures this ensures rapid location of the global maxima. In practice, we perform swaps only between adjacent temperature chains every 1000 iterations. The true posterior samples will come from the β = 1 chain but the higher temperature chains can be used to evaluate the evidence via thermodynamic integration (see, e.g., Littenberg & Cornish 2009 and references therein). Consider the evidence for a chain with temperature 1/β as part of a partition function

$$\begin{eqnarray} Z(\beta)&=&\int d\Theta \, p(d|\Theta,\mathcal {H},\beta)p(\Theta |\mathcal {H})\nonumber\\ &=&\int d\Theta \, p(d|\Theta,\mathcal {H})^{\beta }p(\Theta |\mathcal {H}). \end{eqnarray} \tag{ C11 }$$

Since the prior is independent of β, we can take the log and integrate over β to obtain

$$\begin{equation} \ln \,p(d|\mathcal {H})=\int _0^1d\beta \, \langle \ln \, p(d|\Theta,\mathcal {H}) \rangle _{\beta }, \end{equation} \tag{ C12 }$$

where $\langle \ln \, p(d|\Theta,\mathcal {H}) \rangle _{\beta }$ is the expectation value of the likelihood for the chain with temperature 1/β. The expectation values are calculated over the post burn-in chains.

In practice, it is important to choose a temperature ladder such that we explore the entire likelihood surface and recover the full integrand of Equation (C12). Here we will closely follow Littenberg & Cornish (2010) in the construction of our temperature ladder and diagnostic techniques. In constructing a temperature ladder to be used with thermodynamic integration it is important to understand that there are two regimes that we are interested in (at least in the GW detection problem). The first regime is the range of temperatures in which the (tempered) likelihood is still in "contact" with the GW, that is, the data still inform on the GW parameters. Since this is where the bulk of the integrand is concentrated when a signal is present, it is very important that we choose a fine temperature spacing here to resolve the point at which the likelihood loses contact with the GW. To do this we choose a geometrically spaced temperature ladder with temperature spacing

$$\begin{equation} \Delta T = 1 + \sqrt{\frac{2}{n_{\rm dim}}}, \end{equation} \tag{ C13 }$$

where n_dim is the number of dimensions in our search. Now that a temperature spacing is defined we must choose a maximum temperature T_max for this regime. This choice is based on the expected maximum S/N of a GW signal in the data. Since $\rho \propto \sqrt{\ln \,p(\delta t | boldsymbol {\lambda }, boldsymbol {\phi })}$, the effective S/N for a chain at temperature T is $\rho _{\rm eff}\propto \rho /\sqrt{T}$, therefore, for a chosen maximum S/N we have

$$\begin{equation} T_{\rm max} = \left(\frac{\rho _{\rm max}}{\rho _{\rm eff,max}}\right)^2, \end{equation} \tag{ C14 }$$

where ρ_{eff, max} is the S/N at which we lose contact with the GW signal. For this work we have chosen ρ_max = 10 and ρ_{eff, max} = 3, corresponding to T_max = 11.1. These values were chosen based on the fact that if we indeed have a signal with ρ ≈ 10, it would have likely been detected in previous PTA data sets, and ρ_{eff, max} is chosen based on trial and error and simulations.

If we were only interested in parameter estimation then we would cut off the temperature ladder here; however, for evidence evaluation we must explore the full parameter space. This is the second temperature regime of evidence evaluation via thermodynamic integration. This is most important when varying the noise parameters as well as the GW parameters. Here, we must choose an overall maximum temperature such that we are effectively sampling from the prior. In other words, the temperature must be sufficiently high such that the average log-likelihood has become constant with respect to increasing temperature. In this regime we choose a more coarse temperature spacing with ΔT = 1.5 and geometric spacing. As noted in Littenberg & Cornish (2010) a good diagnostic to ensure that we are using a high enough temperature is to plot the mean log-likelihood for each temperature chain versus β.

In Figure 9 we plot the mean log-likelihood versus β for GW plus noise (gray (blue)) and noise (black (green)) models. First we notice that at low temperatures (high β) the GW plus noise model fits the data better based on the higher likelihood values (the data has an S/N 10 GW injection) but that it has slightly lower values at high temperature (low β) because of the expanded prior volume due to the GW parameter space. Since the Bayesian evidence is the area under these curves the question that is being answered by computing a Bayes factor is "does the fact that the GW plus noise model fit the data better (low temperature regime) overcome the fact that that model has a larger prior volume (high temperature regime)?" Because of this, it is crucial that we include temperatures high enough so that the average log likelihood becomes constant, indicating that we are sampling the prior distribution. For this work, we find that a maximum temperature of ∼10⁵ is sufficient.

Zoom In Zoom Out Reset image size