On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array

While detections of gravitational waves (e.g., Abbott et al. 2016) have been made with ground-based interferometers that are sensitive to hertz–kilohertz gravitational waves, experiments that operate at lower frequencies have yet to identify a signal. Pulsar timing array (PTA) experiments, which monitor and measure arrival times from millisecond pulsars (MSPs), have been established to search for signals in the nanohertz band. This is the domain of the stochastic background from supermassive black hole binaries, which is expected to be the first gravitational-wave signal detected with PTAs (Rosado et al. 2015). The background manifests as a temporally and spatially correlated process in the MSP arrival times. The strain spectrum of such a background is predicted to have the power-law form $h{(f)=A(f/1\,{\mathrm{yr}}^{-1})}^{-2/3}$, where f is the gravitational-wave frequency and A is the strain amplitude at f = 1 yr⁻¹ (Phinney 2001). The amplitude, A, will depend on the demographics of the supermassive black hole population. Astrophysical effects relating to supermassive binary black hole evolution may cause deviations from a power law (e.g., Ravi et al. 2014; Sampson et al. 2015; Taylor et al. 2017b).

A definitive detection of the gravitational-wave background is the presence of specific spatial correlations in the arrival times (Hellings & Downs 1983). Other processes can produce signals with similar temporal properties, with either no (Shannon & Cordes 2010) or different spatial correlation. Tiburzi et al. (2016) and Taylor et al. (2017a) showed the challenges in distinguishing between sources that produce different spatial correlations.

Previous searches for the gravitational-wave background from supermassive black hole binaries have reported limits on the strain amplitude, A, ranging between 1.1 × 10⁻¹⁴ and 1.0 × 10⁻¹⁵ at 95% of either confidence or credibility, where appropriate (Jenet et al. 2006; van Haasteren et al. 2011; Demorest et al. 2013; Shannon et al. 2013; Lentati et al. 2015; NANOGrav Collaboration et al. 2015; Shannon et al. 2015; Arzoumanian et al. 2018). The limits are now known to be affected by systematic uncertainties in the ephemeris of the solar system, which impacts pulsar timing because the arrival times are necessarily referenced to an inertial frame located at solar system barycenter (e.g., Arzoumanian et al. 2018).

In a recent analysis of the NANOGrav 12.5 yr data set, a common-noise process was reported having a Bayes factor greater than 10⁴ (this corresponds to 9 on the commonly used natural logarithmic scale), with the signal persisting even when accounting for uncertainties in the solar system ephemeris (Arzoumanian et al. 2020). We discuss the meaning of the term "common process" in detail later. Here we define the symbol CP to represent the common process as obtained by the hypotheses used in the NANOGrav analysis. Evidence for Hellings–Downs correlations was insignificant.

We have carried out a similar analysis using the Parkes Pulsar Timing Array (PPTA; Manchester et al. 2013) second data release (Kerr et al. 2020). The observations and methodology are described in Section 2. In Section 3 we discuss the results of the searches. In Section 4, we discuss limitations in the methodology. In particular, we demonstrate through simulation that the search methods can spuriously detect a common red process in timing array data sets in which it is absent.

The PPTA project monitors an ensemble of MSPs with the 64 m Parkes radio telescope (also named Murriyang) in New South Wales, Australia. The data used, namely pulse arrival times from the observations, were acquired between 2003 and 2018 and were published as part of the second data release of the project (PPTA-DR2; Kerr et al. 2020). Observations were taken at a cadence of approximately three weeks. At each epoch, data were usually recorded in bands centered at three different radio frequencies in order to correct variations in pulsar dispersion measures (Keith et al. 2013). Data were recorded with a series of digital processing systems, with quality having improved over the course of the project.

We analyzed the data set using methodology that was based on that applied to the NANOGrav 12.5 yr data set (Arzoumanian et al. 2020), which itself was based on NANOGrav Collaboration et al. (2015) and Taylor et al. (2017a). Stochastic signals were modeled as being correlated (red) or uncorrelated (white) in time. We had previously characterized the noise processes for individual pulsars in the PPTA sample (Goncharov et al. 2021). That analysis showed that the PPTA data sets contain a wide variety of noise processes, including instrument-dependent or band-dependent processes. In this work we included red-noise processes in all pulsars, even for those pulsars that showed no evidence for such noise in previous analyses.

As in Goncharov et al. (2021), we assume that the power spectral density of all red processes follows a power law, parameterized such that the amplitude, A, is in units of gravitational-wave strain at 1 yr⁻¹:

$$\begin{eqnarray}&&P(f| A,\gamma )={\rm{\Gamma }}({\zeta }_{{ab}})\displaystyle \frac{{A}^{2}}{12{\pi }^{2}}{\left(\displaystyle \frac{f}{{\mathrm{yr}}^{-1}}\right)}^{-\gamma }{\mathrm{yr}}^{3}.\end{eqnarray} \tag{ 1 }$$

The fluctuation frequency of the pulse arrival time power spectrum is denoted f and the spectral index is − γ. The noise terms are modeled using a Fourier series, starting with a fundamental frequency that is the inverse of the observation span corresponding to the entire pulsar data set, T_obs. We use n_c = 30 harmonics if γ > 1.5 (Goncharov et al. 2020), otherwise we use n_c from the single-pulsar analysis (Goncharov et al. 2021). The overlap reduction function (Γ(ζ_ab ); Finn et al. 2009) characterizes the spatial correlation of the signal, and depends on the angular separation, ζ_ab , of two pulsars a and b with respect to the observer. For an isotropic stochastic background from the gravitational waves of general relativity (Hellings & Downs 1983; Jenet & Romano 2015),

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Gamma }}}_{\mathrm{GWB}}({\zeta }_{{ab}})=\displaystyle \frac{1}{2}-\displaystyle \frac{1}{4}\left(\displaystyle \frac{1-\cos {\zeta }_{{ab}}}{2}\right)\\ \quad +\displaystyle \frac{3}{2}\left(\displaystyle \frac{1-\cos {\zeta }_{{ab}}}{2}\right)\mathrm{ln}\left(\displaystyle \frac{1-\cos {\zeta }_{{ab}}}{2}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

when a ≠ b.

Our analysis proceeded through these steps:

• 1.  
  We first searched for a common power-law, red-spectrum stochastic process (CP1), with an identical power spectrum and unrelated temporal evolution or spatial correlation across pulsars. We emphasize that the timing spectrum of the process is assumed to be a statistically identical ensemble-average power spectrum among pulsars, which would be the case for a gravitational-wave background. ²¹ Throughout our analysis we marginalized over deterministic terms in the timing model (Reardon et al. 2021). This included instrument-dependent offsets ("jumps") of unknown value, as identified by Kerr et al. (2020). Offsets with a priori measured values were held fixed. We trialed both marginalizing over the white-noise parameters and also holding them fixed at their maximum a posteriori values. We obtained consistent results between these approaches.
• 2.  
  Following the NANOGrav analysis we also assumed a power-law model with γ = 13/3. This value has astrophysical interest as it is the expected value for a gravitational-wave background caused by supermassive binary black holes (Phinney 2001). The resulting common power-law, red-noise process is here labeled as CP2. Based on a factorized-likelihood approach, we performed a dropout analysis to evaluate the consistency of individual PPTA DR2 pulsars with the signal identified by CP2 (see Arzoumanian et al. 2020).
• 3.  
  We measured the amplitude of individual Fourier components ${P}_{i}({f}_{i}| {\rho }_{i})={\rho }_{i}^{2}\,{T}_{\mathrm{obs}}$ of a common process, at each frequency f_i separately in order to determine whether the power-law assumption of CP1,2 is valid.
• 4.  
  We searched for evidence that CP2 exhibits spatial correlations from a gravitational-wave background, a monopolar signal, MP, or a dipolar signal, DP. In this analysis we held the white-noise stochastic components fixed at their maximum a posteriori value to reduce the number of parameters in the search and reduce computation time.
• 5.  
  To assess the shape of any spatial correlations, we measured Γ(ζ_ab ) at seven equally separated "node" angles between 0° and 180° inclusive, using linear interpolation to determine Γ(ζ_ab ) from pulsar pairs between the nodes. The interpolant modeling of the PTA correlation curve was first done in Taylor et al. (2013).

2.1. Comparison with the NANOGrav Data Set and Processing Methods

Following the process described above, we obtained strong evidence for CP1, with $\mathrm{ln}{{ \mathcal B }}_{\varnothing }^{\mathrm{CP}1}\gt 15.0$. This implies an odds ratio in favor of CP1 of >3 × 10⁶: 1, if both models had even odds a priori. The results from the parameter estimation are shown in the left panel in Figure 1. We obtain ²² ${\mathrm{log}}_{10}{A}_{\mathrm{CP}1}=-{14.55}_{-0.23}^{+0.10}$ and $\gamma ={4.11}_{-0.41}^{+0.52}$. The results from the NANOGrav analysis are overlaid and have significant overlap, although the NANOGrav analysis preferred a steeper spectral exponent, when the latter analysis was conducted with five Fourier components. Unlike in Arzoumanian et al. (2020), our measurements are consistent when we change the numbers of fluctuation frequencies used to model the common process.

Zoom In Zoom Out Reset image size

The right panel in Figure 1 represents the parameter estimation for the free-spectral model. It is challenging to obtain a complete noise model for the brightest MSP, PSR J0437−4715 (Goncharov et al. 2021), and so we show this spectrum with, and without, the inclusion of this pulsar. We overlay the astrophysically interesting spectrum corresponding to γ = 13/3, along with the frequencies corresponding to the orbital periods of the planets.

In the left-hand panel of Figure 2, we show the posterior distribution for ${\mathrm{log}}_{10}A$, assuming a power-law model with γ = 13/3. The measured CP2 amplitude is ${\mathrm{log}}_{10}A=-14.66\pm 0.07$ corresponding to $A={2.2}_{-0.3}^{+0.4}\times {10}^{-15}$, which is consistent with that measured in the NANOGrav 12.5 yr data set.

Zoom In Zoom Out Reset image size

The NANOGrav analysis attempted to determine which pulsars contributed to this signal by calculating dropout factors, which for our sample are displayed in the right-hand panel of Figure 2. The pulsars with the smallest dropout factors (PSRs J1824−2452A and J1939+2134) are known to have strong timing noise inconsistent in strength with that of other pulsars. However, the pulsars with the highest dropout factors include pulsars with high (PSR J0437−4715) and low timing precision (PSR J1022+1001), and pulsars with shorter timing baselines (PSR J2241−5236). The meaning and use of dropout factors is further discussed in Section 4.3.

The results from the model-independent parameter estimation of the overlap reduction function (obtained assuming γ = 13/3) are provided in Figure 3. The left-hand panel shows the inferred spatial correlations. They were sampled at seven node angles, whereas spatial correlations for other angles were obtained with linear interpolation. The Hellings–Downs relation is overplotted. The right-hand panel shows inferred amplitude of the Hellings–Downs spatially correlated noise and that for CP2. With the entire PPTA sample of pulsars the Bayes factor is $\mathrm{ln}{{ \mathcal B }}_{\mathrm{CP}2}^{\mathrm{HD}}=0.3$, which provides no significant evidence for, or against, Hellings-and-Downs correlations. We note that if PSR J0437−4715 is removed from the sample then the Bayes factor increases to $\mathrm{ln}{{ \mathcal B }}_{\mathrm{CP}2}^{\mathrm{HD}}=1.0$. The data strongly disfavor CP2 having monopole or dipole spatial correlations ($\mathrm{ln}{{ \mathcal B }}_{\mathrm{CP}2}^{\mathrm{MP}}$ and $\mathrm{ln}{{ \mathcal B }}_{\mathrm{CP}2}^{\mathrm{DP}}$ both <−10).

Zoom In Zoom Out Reset image size

Under the same assumptions the Bayesian analyses of both the PPTA and NANOGrav data sets show a preference for models that include a common-noise process in addition to individual noise terms. The CP1 model has consistent spectral index and amplitude between the data sets and therefore we can exclude this signal from being telescope dependent. Given the different strategies employed in mitigating interstellar propagation effects by the two projects (both in terms of choice of observing band and methods for correcting for dispersion-measure variations), it is also unlikely that the noise is associated with the interstellar medium.

However, we are attempting to detect a common-noise process from a single realization of the process in each pulsar. The noise process is strongest at the lowest fluctuation frequency, so the process is being characterized on a timescale comparable to the typical data span. This greatly complicates tests of the noise modeling. Consequently there are a number of caveats for interpreting the CP1 and CP2 results as discussed in the following subsections. We conclude by discussing our search for spatial correlations in the data.

4.1. Are the Models of the Intrinsic Noise Correct and Complete?

4.2. What Assumptions Lead to the Evidence of a Common Process?

The assessment of whether a common red-noise process is present is based on the null hypothesis that all pulsars are affected by independent red timing noise, modeled by a power-law spectrum described by amplitudes and exponents drawn from a uniform prior. The detection hypothesis is that all pulsars also exhibit a red process with the same power-law spectrum, where amplitudes and exponents are drawn from the delta-function prior. There is a possibility that neither of these models nor their priors are sufficient descriptions of the data, which is often referred to as model misspecification (e.g., Davidson & MacKinnon 1981; Müller 2013). In simulations we demonstrate that noise without a statistically identical spectrum between pulsars can be misinterpreted as the common red process.

We simulated 10 realizations of timing residuals for the 26 PPTA-DR2 pulsars, with a range of realistic white noise levels, and injected power-law timing noise models with amplitudes and spectral indices drawn uniformly across approximately several orders of magnitude (${\mathrm{log}}_{10}A$ spanning approximately −16 to −13 and γ spanning 3 to 5). The power spectral densities of the simulated residuals are shown in the left-hand panel of Figure 4. We performed model selection for a model with a common-spectrum process along with intrinsic pulsar timing noise (CP2), against a model with intrinsic timing noise only (∅), and obtained $\mathrm{ln}{{ \mathcal B }}_{\varnothing }^{\mathrm{CP}2}\gt 13.5$ in all realizations, implying that our methodology can detect a "common" process if the properties of the noise are broadly similar but far from identical, with the amplitude of the noise varies by three orders of magnitude. Figure 5 shows the recovered common-noise spectrum and the injected timing noise models. We continued to increase the spread in injected noise amplitudes upward from ≈3 orders of magnitude, and found that common noise is disfavored when the spread in amplitude exceeds 5 orders of magnitude. Thus, the simulated data only favors the correct hypothesis when the range of simulated noise parameters starts to resemble the uniform red noise prior assumed in recent analyses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is physically likely that intrinsic pulsar timing noise has similar, but not identical properties between different pulsars (Shannon & Cordes 2010). As such, a second null hypothesis (not tested by the current analysis) would be that each pulsar has independent red noise, but the properties of that noise cluster in a similar range. This could be examined by assuming that the amplitude and spectral index of the noise terms for the pulsars are not identical, but are drawn from a distribution. The nonuniform noise hypothesis is distinct both from the signal hypothesis with the delta-function prior and from the noise hypothesis with the uniform prior. For example, the new noise prior could be modeled as a Gaussian distribution of width μ_A and μ_γ , and variance ${\sigma }_{A}^{2}$ and ${\sigma }_{\gamma }^{2}$. If the variance is inconsistent with zero it would suggest that the noise processes were not common but just similar. If it is consistent with zero then values could be used to constrain the properties of a common process.

As a gravitational-wave background signal will affect all pulsars, it will be apparent not only as a common spectral process, but also as a "noise floor". The noise level of a given pulsar should not be below this floor apart from statistical fluctuations, including instances where two noise processes cancel each other out. We updated the simulations shown in Figure 4 by including timing noise with identical power-law spectral densities (${\mathrm{log}}_{10}A=-15.51$, γ = 5.5) for 25 pulsars, and a lower, shallower-spectrum timing noise (${\mathrm{log}}_{10}A=-14.38$, γ = 1.5) in the pulsar with the lowest white noise levels. The power spectral densities corresponding to these simulations are shown in the right-hand panel of Figure 4. We performed model selection for CP2 over ∅, and again found significant support for the common-spectrum noise process ($\mathrm{ln}{{ \mathcal B }}_{\varnothing }^{\mathrm{CP}2}\,\gt$ 13.1 across all realizations), which is comparable to the support found in both the NANOGrav and PPTA analyses. As the power spectral density in the lowest frequency channel for the low-noise pulsar is 4 orders of magnitude lower than the common signal, the model selection is not explicitly identifying a noise floor. However, a factorized-likelihood dropout analysis showed that the simulated pulsar with a low noise level was not consistent with the retrieved common-spectrum process (with a dropout factor <1; see the following section).

4.3. Do We Understand which Pulsars Are Contributing to the Evidence?

4.4. Are We Affected by the Choice of Solar System Ephemeris?

4.5. Searching for Spatial Correlations

Under the assumptions of an analysis of the NANOGrav 12.5 yr data set (Arzoumanian et al. 2020), we have detected with high confidence a common-spectrum time-correlated signal in the timing of the 26 PPTA-DR2 millisecond pulsars. However, as noted above, there are some important caveats that need to be addressed before the signature can be confidently attributed to a physical process common to all the pulsars in the array. We do not confirm or rule out that the common-spectrum process is a spatially correlated stochastic gravitational-wave signal. However, the identified process does not possess monopole or dipole correlations and is not caused by errors in masses and trajectories of Mars, Jupiter, and Saturn, that would have resulted in deterministic timing residuals according to bayesephem models (Vallisneri et al. 2020). Proposed follow-up work relates to (1) improving our understanding of the properties of the intrinsic timing noise in millisecond pulsars and (2) identifying the optimal model comparisons and methodologies that can determine whether the noise detected corresponds to a "noise floor" and is identical in all pulsars.

The PPTA project now has nearly three additional years of data obtained with a higher sensitivity wide-band system (Hobbs et al. 2020) that can be added with the data presented to increase timing baselines. We will also combine the PPTA data sets with observations from other observatories as part of the International Pulsar Timing Array project (Perera et al. 2019), with the latter being ideal to continue this work. Such lengthened and more sensitive data sets will allow us to probe timescales significantly longer than the orbital period of Jupiter and closer to that of Saturn. We will be able to compare noise models obtained from a wide range of telescopes and maximize our chance of determining the nature of the red-noise signals that are present in our data.

These data sets will also enable more sensitive and robust searches for the Hellings–Downs spatial correlations necessary to make a definitive detection of the gravitational-wave background.

This work has been carried out by the Parkes Pulsar Timing Array, which is part of the International Pulsar Timing Array. The Parkes radio telescope (Murriyang) is part of the Australia Telescope, which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO. This paper includes archived data obtained through the CSIRO Data Access Portal (data.csiro.au). We acknowledge the use of chainconsumer (Hinton 2016). Parts of this research were conducted by the Australian Research Council (ARC) Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. R.M.S. acknowledges support through ARC future fellowship FT190100155. R.S. acknowledges support through the ARC Laureate fellowship FL150100148. The author list is based on three tiers, which correspond to primary contributors, to members of the collaboration who provided feedback, and to members of the collaboration with significant observing records. The data-processing code that was used in this work is available at github.com/bvgoncharov/correlated_noise_pta_2020.