The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background

ABSTRACT

We searched for an isotropic stochastic gravitational wave background in the second data release of the International Pulsar Timing Array, a global collaboration synthesizing decadal-length pulsar-timing campaigns in North America, Europe, and Australia. In our reference search for a power-law strain spectrum of the form |$h_c = A(f/1\, \mathrm{yr}^{-1})^{\alpha }$|⁠, we found strong evidence for a spectrally similar low-frequency stochastic process of amplitude |$A = 3.8^{+6.3}_{-2.5}\times 10^{-15}$| and spectral index α = −0.5 ± 0.5, where the uncertainties represent 95 per cent credible regions, using information from the auto- and cross-correlation terms between the pulsars in the array. For a spectral index of α = −2/3, as expected from a population of inspiralling supermassive black hole binaries, the recovered amplitude is |$A = 2.8^{+1.2}_{-0.8}\times 10^{-15}$|⁠. None the less, no significant evidence of the Hellings–Downs correlations that would indicate a gravitational-wave origin was found. We also analysed the constituent data from the individual pulsar timing arrays in a consistent way, and clearly demonstrate that the combined international data set is more sensitive. Furthermore, we demonstrate that this combined data set produces comparable constraints to recent single-array data sets which have more data than the constituent parts of the combination. Future international data releases will deliver increased sensitivity to gravitational wave radiation, and significantly increase the detection probability.

gravitational waves, methods: data analysis, pulsars: general

1 INTRODUCTION

Inspiralling supermassive black hole binaries (SMBHBs) with masses larger than 10⁷M_⊙ are expected to generate the strongest gravitational-wave (GW) signals in the Universe. The incoherent superposition of all of these inspiralling SMBHBs should generate a stochastic GW background (GWB) that is the strongest in the nanohertz frequency band (e.g. Rajagopal & Romani 1995; Jaffe & Backer 2003; Sesana, Vecchio & Colacino 2008; Burke-Spolaor et al. 2019). Other sources that could also produce a stochastic background in the nanohertz band are cosmic strings (e.g. Ölmez, Mandic & Siemens 2010), cosmological phase transitions, and a primordial background produced by quantum fluctuations in the gravitational field in the early universe (e.g. Grishchuk 2005; Lasky et al. 2016). For comparison, the Laser Interferometer Gravitational-Wave Observatory (LIGO) and the Virgo Collaboration, which are terrestrial GW detectors and have detected GWs from merging stellar mass black holes and neutron stars (e.g. Abbott et al. 2019, 2021), are only sensitive to GW signals that are ten orders of magnitude higher in frequency than PTAs.

A nanohertz GWB can be detected using a precisely timed ensemble of millisecond pulsars (Sazhin 1978; Detweiler 1979), called a pulsar timing array (PTA; Foster & Backer 1990). The GWs distort the space-time between the Earth and pulsars, changing their proper distance, thereby leading to a measurable deviation of the pulsar pulse arrival times. Since such effects cannot be detected with confidence using only one pulsar, PTAs leverage the imprint of spatially correlated timing deviations between pulsars which are separated by kiloparsec distances across the galaxy, yet are subject to the common influence of the GWB.

An isotropic GWB manifests itself as a long time-scale, low-frequency (or red) common signal across the pulsars in a PTA. This common signal is characterized by the common spectrum and the interpulsar spatial correlations. For an isotropic GWB, these spatial correlations are unique, referred to as the Hellings & Downs (1983) (HD) correlations, and thus are considered to be the ‘smoking gun’ signature for the presence of a GWB (Tiburzi et al. 2016) in any PTA data set. The spectral amplitude of this common signal is determined by the characteristic strain, h_c(f), of the GWB, which itself is a function of the physics sourcing the GWB (e.g. SMBHB masses, merger time-scale, and number density) (e.g. Sesana 2013a; Kelley et al. 2017; Chen, Sesana & Conselice 2019). Thus, precise spectral characterization of the GWB will allow us to extract the underlying astrophysics of the background, as well as distinguish between different sources of the GWB (e.g. Pol et al. 2021).

The ability to detect GWs relies on, among other things, the number of pulsars available to cross-correlate in GWB searches, and on the length of each pulsar data set (Siemens et al. 2013). Improvements in both of these parameters increases the detection significance of the GWB signal which in turn allows for better constraints on the parameters of the GWB spectrum (Pol et al. 2021). Hence, international efforts spanning decades from the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), North American Nanohertz Observatory for Gravitational Waves (NANOGrav; Arzoumanian et al. 2016), and the Parkes Pulsar Timing array (PPTA; Manchester et al. 2013), as well as newer PTAs such as the Indian Pulsar Timing Array (InPTA; Joshi et al. 2018), Chinese Pulsar Timing Array (CPTA; Lee 2016), and with the MeerKAT Interferometer in South Africa (Bailes et al. 2020) share and combine their data to form the International PTA (IPTA; Hobbs et al. 2010).

In this spirit of international collaboration, the IPTA has produced two data sets to date. The first IPTA data release (DR1; Verbiest et al. 2016) consisted of 44 millisecond pulsars and yielded no conclusive detection of a GWB. The second IPTA data release (DR2; Perera et al. 2019) consists of 65 pulsars and is the focus of this analysis. The pulsars in DR2 have data sets spanning 0.5−30 yr. For the first time, we process the data subsets from each individual PTA and search for a GWB in a self-consistent way, thus enabling us to make a fair comparison of respective PTA constraints.

Recently, a spatially uncorrelated (pulsar-weighted-average) spectrally similar common process or common-spectrum process (CP) was detected in the NANOGrav 12.5-yr data set (Arzoumanian et al. 2020), the second data release of the Parkes Pulsar Timing Array (Goncharov et al. 2021), and the EPTA six-pulsar data set of the second data release (Chen et al. 2021). The process is modelled as an additional time-correlated term with the same power spectrum in all of the pulsars. However, there is little evidence to support the existence of spatial HD correlations in any of these data sets. We compare the IPTA DR2 constraints on the GWB with those obtained from these analyses.

The paper is organized as follows: in Section 2 we give an overview of the second IPTA data release, hereafter referred to as DR2. We describe our data analysis methods in Section 3, and give our results in Section 4. Caveats and implications of our analysis and results are discussed in Section 5, including the astrophysical interpretation of a potential GWB. The conclusion is given in Section 6.

2 IPTA DATA RELEASE 2

IPTA DR2 includes a combination of timing data from the following individual PTA data releases: the EPTA data release 1.0 (Desvignes et al. 2016), the NANOGrav 9-yr data set (Arzoumanian et al. 2015), and the PPTA first data release (Manchester et al. 2013) and its extended version (Reardon et al. 2016). The EPTA data set includes high-precision timing data from 42 MSPs obtained with the largest radio telescopes in Europe – Effelsberg telescope, Lovell telescope, Nançay telescope, and Westerbork Synthesis telescope – covering data from 1996 to 2015 with a time baseline between 7 and 18 yr. In addition to these data, archival timing data of PSR J1939+2134 since 1994 was included. The NANOGrav 9-yr data set includes high-precision timing observations from 37 MSPs obtained with the Robert C. Byrd Green Bank Telescope and the Arecibo telescope, spanning a time baseline between 0.6 and 9.2 yr, covering the data from 2004 to 2013. In addition, the long-term timing data of PSR J1713+0747 from Zhu et al. (2015) and the data of PSRs J1857+0943 and J1939+2134 from 1984 through 1992 (Kaspi, Taylor & Ryba 1994) were included. The PPTA data set includes high-precision timing observations from 20 MSPs obtained with the Parkes radio telescope (also known as Murriyang) from 2004 to 2011. IPTA DR2 also included single frequency band (1.4 GHz/L-band) Parkes telescope legacy data obtained since 1994. The additional 3.0 GHz timing data reported in Shannon et al. (2015) for PSRs J0437−4715, J1744−1134, J1713+0747, and J1909−3744 were also included in the data set. In total, the timing data from 65 MSPs were included in IPTA DR2, which has 21 more source than the IPTA DR1 (Verbiest et al. 2016). There are 27 and 7 MSPs in IPTA DR2 with a timing baseline >10 yr and >20 yr, respectively. All pulsars were observed at multiple frequencies. All EPTA and PPTA observations were averaged in time and frequency to obtain a single time-of-arrival (TOA) for each receiver and observation. The NANOGrav observations were averaged in time and included sub-band information, i.e. averaged in frequency to maintain a frequency resolution ranging from 1.5 to 12.5 MHz depending on the receiver and backend instrument combination, resulted in a single TOA for each frequency channel. More details about the constituent PTA data sets can be found in Perera et al. (2019).

The different data sets for a given pulsar in IPTA DR2 were combined by fitting for time offsets, referred to as JUMPs, in the timing model to account for any systematic delays between data sets. The highest weighted data set with the lowest sum of TOA uncertainties was used as the reference data set in this process. The timing models of pulsars included astrometric parameters, rotational frequency information, dispersion measure information, and Keplerian and Post-Keplerian parameters if the pulsar is in a binary system. For NANOGrav observed pulsars, ‘FD’ parameters were included to minimize the effect of frequency-dependent profile variations of pulsars (see Arzoumanian et al. 2015). IPTA DR2 produced two data set versions depending on different methods of handling the dispersion measures (DM) variations of pulsars over time (VersionA and VersionB – see Perera et al. 2019, for details). In VersionA, the DM variations of pulsars were determined using DMMODEL described in Keith et al. (2013) and the noise parameters for different data sets were directly taken from their original data releases. In VersionB, the DM variations were modelled using the first two time derivatives of the DM and a time-correlated stochastic DM process in the timing model. The noise parameters were also re-estimated based on the new IPTA data combination in this version. We use VersionB for this work.

3 DATA ANALYSIS METHODS

In this work, we follow the conventions established by other pulsar timing array data analyses (i.e. Lentati et al. 2015; Arzoumanian et al. 2016, 2020; Goncharov et al. 2021). The multivariate Gaussian likelihood |$\mathcal {L}(\boldsymbol {\delta t} | \boldsymbol {\theta })$| is employed to model noise and signal contributions, parametrized by |$\boldsymbol {\theta }$|⁠, to the observed timing residuals. Our likelihood was of the same form as other PTA analyses, (e.g. Arzoumanian et al. 2015). We used enterprise (Ellis et al. 2020) to evaluate the likelihood and priors, and PTMCMCSampler (Ellis & van Haasteren 2017) to perform a Markov chain Monte Carlo (MCMC) simulation, drawing samples from the posterior probability distribution. Model selection was performed via the product-space sampling method (Carlin & Chib 1995; Hee et al. 2016). Additionally, we used the Savage-Dickey approximation to the Bayesian evidence ratio when appropriate.

3.1 Noise models

For each pulsar, we modelled the TOAs with a combination of four processes: the timing model, white noise, intrinsic red noise, and DM variations. Deterministic contributions from the timing model, described in Section 2, were analytically marginalized (van Haasteren et al. 2009). The time-uncorrelated white noise was modelled with EFAC and EQUAD, and ECORR parameters for NANOGrav pulsars with their sub-banded TOAs (definitions of EFAC, EQUAD, and ECORR can be found in, e.g. Verbiest et al. 2016; Perera et al. 2019). Every observing receiver and backend system combination is given its own set of white noise parameters. The time-correlated red noise process (e.g. pulsar spin noise; Shannon & Cordes 2010) and stochastic DM variations (Keith et al. 2013) were modelled as Fourier basis Gaussian processes. In each case Fourier spectrum coefficients were modelled as power laws,

$$\begin{eqnarray*} P_{\rm RN}(f) &= \frac{{A_{\rm RN}}^2}{12 \pi ^2} f_\mathrm{yr}^{-3} \left(\frac{f}{f_\mathrm{yr}}\right)^{-\gamma _{\rm RN}}, \end{eqnarray*}$$

(1)

$$\begin{eqnarray*} P_{\rm DM}(f,\nu) &=\frac{{A_{\rm DM}}^2}{12 \pi ^2} f_\mathrm{yr}^{-3} \left(\frac{f}{f_\mathrm{yr}}\right)^{-\gamma _{\rm DM}} \left(\frac{1400\ {\rm MHz}}{\nu }\right)^2 , \end{eqnarray*}$$

(2)

where A is the power-law amplitude, γ its spectral index, and f_yr = 1 yr⁻¹ ≈ 3.17 × 10⁻⁸ Hz. The difference between these two processes lies in the radio frequency ν dependence. Intrinsic red noise is achromatic, i.e. frequency independent, while DM variations follow a ν⁻² dependence (e.g. Lentati et al. 2016).

Despite all MSPs in IPTA DR2 exhibiting high rotational stability, such that the marginalized timing model, red, and DM noise terms are in general sufficient, certain pulsars have been found to experience timing events that need to be included in their data model. Of interest to this analysis is PSR J1713+0747, which was observed to experience multiple sudden drops in apparent DM with an exponential recovery (Demorest et al. 2013; Lam et al. 2018; Goncharov et al. 2021). Only the first such event lies within the time-span of the IPTA DR2 and was included as an additional deterministic term to the full noise model of PSR J1713+0747. The amplitude, epoch, and recovery time-scale of the DM exponential dip are sampled simultaneously with the pulsar red and DM noise terms.

Lentati et al. (2016) found additional sources of red noise in IPTA DR1. These include radio frequency band-dependent and observing system-dependent terms, which may affect measurements of P_RN and P_DM, if not modelled. It is possible that mismodelling these effects can bias recovery of the CP. More prescriptive models for the CP should be less affected by this bias. Recent PTA analyses have included more complex red noise and DM variation models, where different pulsars in the array use different models (e.g. Aggarwal et al. 2019; Goncharov et al. 2021). In the name of computational efficiency, we opted to use the same power-law models for all pulsars except when absolutely necessary, as is the case for PSR J1713+0747.

IPTA DR2, being the combination of data from multiple telescopes and many observing systems, has larger model parameter space than its constituent data sets. The large number of model parameters and TOAs increases the computational complexity of the analysis. As we searched for long-term processes, such as the GWB, we limited our analysis to pulsars whose observation time exceeded 3 yr. This reduced the number of pulsars from the full 65 in DR2 to 53. Additionally, we fixed the white noise parameters (EFAC, EQUAD, and ECORR) to median a posteriori values from single pulsar analyses. Both of these choices reduced the analysis parameter space to a more manageable size.

3.2 Common-spectrum process models

In addition to modelling noise intrinsic to the individual pulsars, we also include a red CP that is present in all of the pulsars. The source of this process could be the GWB, or any other common noise that manifests itself in all pulsars, such as clock errors (Caballero et al. 2016; Hobbs et al. 2020) or errors in the Roemer delays from Solar-system ephemeris (SSE) systematics (Tiburzi et al. 2016; Vallisneri et al. 2020). The choice of red noise priors also affects the recovery of a CP due to covariance between pulsar intrinsic red noise and the CP (Hazboun et al. 2020b; Goncharov et al. 2021). Each of these effects can be distinguished by a unique pattern of spatial cross-correlations between pulsars. The cross-power spectral density is defined as,

$$\begin{eqnarray*} S_{ab}(f) = \Gamma _{ab} P_{\rm CP}(f), \end{eqnarray*}$$

(3)

where P_CP is the common-spectrum process and Γ_ab is the overlap reduction function (ORF) describing the interpulsar correlations.

For some analyses, we did not account for any interpulsar correlations, taking Γ_ab = δ_ab to be the identity matrix. For others, we also included different choices of non-diagonal ORFs, such as the quadrupolar Hellings & Downs (1983) correlations that describe a GWB, dipolar correlations associated with SSE errors, or monopolar correlations, Γ_ab = 1, associated with clock errors. In some cases we split the diagonal autocorrelation part of the ORF from the off-diagonal cross-correlation part, treating them as independent processes as a consistency check. When modelling the CP using only the autocorrelations, it is possible to analyse the data from each pulsar independently, then recombine the results to achieve a joint posterior on the CP. We refer to this as the factorized likelihood approach.

We modelled the CP using a Fourier basis Gaussian process, using basis frequencies f = 1/T, 2/T, …, where T is the time-span between the earliest and latest observation in the data set. We model the power spectrum of the CP as a power law using equation (1), replacing the pulsar noise amplitude and spectral index with those from the common process A_CP and γ_CP. In this parameterization of the power spectrum, the characteristic strain spectrum for the GWB is

$$\begin{eqnarray*} h_c(f) = A_\mathrm{CP} \left(\frac{f}{f_\mathrm{yr}} \right)^{(3-\gamma _\mathrm{CP})/2}. \end{eqnarray*}$$

(4)

In some cases we fixed γ_CP = 13/3, equivalent to α = (3−γ_CP)/2 = −2/3, the expected spectrum for a GWB composed of circular supermassive binary black holes (Phinney 2001), and in others we left γ_CP as a free parameter. To determine the number of Fourier frequencies used in the power-law CP model, we fit the power spectrum with a broken power-law model. The broken power law is the sum of the standard, red power law and a white spectrum. This is implemented as a single spectrum with a fixed spectral index at low frequencies that smoothly transitions into a flat, white noise dominated spectrum at high frequencies:

$$\begin{eqnarray*} P(f) = \frac{{A_\mathrm{CP}}^2}{12\pi ^2} f_\mathrm{yr}^{-3} \left(\frac{f}{f_\mathrm{yr}}\right)^{-\gamma _\mathrm{CP}} \left[ 1 + \left(\frac{f}{f_\mathrm{bend}}\right)^{1/\kappa }\right]^{\kappa \, \gamma _\mathrm{CP}}, \end{eqnarray*}$$

(5)

where f_bend is the frequency where the spectral index of the power spectrum changes and κ controls the smoothness of the transition. In this model |$P(f) \sim f^{-\gamma _\mathrm{CP}}$| for f ≪ f_bend, and P(f) constant for f ≫ f_bend. As a verification of our power-law models, we performed a free spectral analysis, where the power at each frequency is fit independently rather than being constrained by a particular spectral shape, P(f).

3.3 Frequentist analyses

As a comparison for our primary Bayesian data analysis pipeline, we performed a frequentist analysis using the noise-marginalized optimal statistic. The optimal statistic is an estimator for the amplitude of the GWB based on the interpulsar correlations (Anholm et al. 2009; Demorest et al. 2013; Chamberlin et al. 2015). Its original derivation assumed the pulsars have no intrinsic red noise. The noise-marginalized optimal statistic uses posterior samples from the Bayesian data analysis to marginalize over the pulsars’ red noise. It has been shown to more accurately estimate the amplitude of the GWB when the pulsars have intrinsic red noise (Vigeland et al. 2018), as is the case in IPTA DR2.

4 RESULTS

The IPTA DR2 data set with its large number of pulsars (53 for this work), long time-span, and various independent observing systems offers a wealth of different analysis opportunities. Here we present a selection of analyses to give a complete picture of what IPTA DR2 teaches us and how it compares to other PTA data sets.

4.1 IPTA DR2 data set

We first show results from the full combined IPTA DR2 data set using methods which differ in their spectral modelling and choice of ORF. The full standard GWB search uses all the available information from the auto- as well as the cross-terms (which are assumed to follow the HD correlation). As the cross-terms, which come from the interpulsar correlations, are the defining feature of the GWB, we insist on their presence in order to confidently claim a detection. However, the autocorrelations are initially the dominant source of information, especially for spectral parameter estimation, and the detection of power in them is considered to be the first hint of a GWB (Romano et al. 2021). It is important to emphasize that detecting the autocorrelations alone is insufficient to claim a detection of a GWB.

4.1.1 Common-spectrum process

To begin, we apply a free spectral model for the CP, measuring the amount of power at each sampling frequency independently, up to 30 frequency bins, the result of which is shown in the top panel of Fig. 1. The common red noise power can be converted into GW strain using equation (4). This alternative representation of the IPTA DR2 analysis can be found in the lower panel of Fig. 1. For reference, we include the predicted sensitivity curve, made using hasasia (Hazboun, Romano & Smith 2019a,b) and the measured white noise parameters of the DR2 data set. Note that the noise power spectral density used in this curve only contains TOA errors, EFAC, EQUAD, and ECORR, and does not contain any estimates for the red noise, as for many pulsars it is difficult, and in fact the point of this analysis, to disentangle intrinsic red noise from a GWB. Hence the low-frequency end of the sensitivity represents a ‘best case’ scenario for comparison. The lowest frequency bin corresponds to the longest time-span, and only two pulsars, J1939+2134 (29 yr) and J1857+0943 (28 yr), have observation baselines sufficient to probe this frequency. However, both have significant RN with spectral indices ∼3.3 (Perera et al. 2019). Therefore, it is not surprising that we do not confidently detect any power there, evidenced by the wide tail extending to low power and median of ∼10⁻⁷ s. However, the second, third, fifth, and eighth frequency bins show power well above the expected sensitivity of IPTA DR2. This could either be the emergence of a GWB or some other unmodelled noise process.

Figure 1.

IPTA DR2 free spectrum and characteristic strain: The top panel shows the power in terms of time of arrival residuals at frequencies in the nanohertz band for the full PTA. The maximum likelihood power law is shown overlaid on posteriors for free spectral parameters, a generic model that measures power at various frequencies without imposing any empirical model. The bottom panel shows the power law and free spectral information from above converted into units of characteristic strain, i.e. the noise power measured in the same units as GW amplitude. The additional line shows the characteristic strain for the detector using the noise parameters for the pulsars. The lower limit of all violins is a result of the lower bound of the prior range for each frequency component.

The CP power spectrum can be modelled with a simple power law using equation (1). Arzoumanian et al. (2020) have shown that the choice of the number of modelled frequencies can affect the constraints on the power-law amplitude and spectral index. Thus, we apply the broken power-law model from equation (5) to find the optimal number of frequency bins for the analysis. Fig. 2 shows the marginalized posterior on the bend frequency. We can identify a clear peak at the 13th frequency in N/T, corresponding to a frequency of 1.4 × 10⁻⁸ Hz, indicated by the orange dashed line. For the remainder of this work, we will limit the search to use only the lowest 13 frequencies with the simple power-law model. This produces constraints equivalent to an analysis using the broken power law, but in a simpler and computationally efficient way.

Figure 2.

Bend frequency: The posterior for the bend frequency parameter in a broken power law search is shown. The peak of the posterior is at the 13th frequency, 13/T for the data set, denoted by the dashed vertical line.

We have also verified that the addition of the BAYESEPHEM SSE model (Vallisneri et al. 2020) to the analysis does not change our results significantly from an analysis that fixed the SSE to DE438. The nearly 30 yr of time-span allows for the separation of the SSE effects from other correlated signals (Vallisneri et al. 2020). For simplicity, we will only show results with DE438, unless stated otherwise.

Fig. 3 compares the results when using two different ORFs. The model that uses only the autocorrelation terms, which we denote CP in Table 1, is very strongly favoured over a model with only intrinsic pulsar noise and no common-spectrum process with log₁₀ Bayes factor of 8.2. Despite the large Bayes factor in favour of the CP, this does not suffice to claim a GWB detection, as we have only used the autocorrelations. This strong evidence only indicates that a number of pulsars have red noise with similar spectral characteristics. We must turn to the cross-correlations to determine if this CP is HD correlated as a GWB should be. Using the full HD ORF containing both auto- and cross-correlations, we find only middling evidence in favour of the auto + cross HD model. The log₁₀ Bayes factor for the full HD model compared to the autocorrelated only CP is 0.3, as shown in Table 1.

Figure 3.

Comparison of common-spectrum process parameters when using autocorrelations only and the full auto + cross-correlated HD model. Left-hand panel: 2D posterior for common-spectrum process power-law parameters. The green lines mark γ = 13/3 and A_CP = 2.8 × 10⁻¹⁵, while the contours represent the 1–, 2–, and 3–σ confidence intervals. Right-hand panel: 1D posterior for common-spectrum process power-law amplitude, using fixed spectral index γ = 13/3.

Table 1.

Bayes factors model comparison: The table shows the logarithmic Bayes factors for a number of model comparisons from the hypermodel and factorized likelihood (marked with an asterisk*) analyses. The preferred model is on the left-hand side of the two models. The brackets indicate the uncertainty in the last digit of the Bayes factors.

Model comparison	|$\log _{10}\rm {BF}$|
HD versus CP	0.3111(6)
CP versus pulsar noise	8.2*
CP versus monopole	4.67(2)
CP versus dipole	2.28(3)

Fig. 3(a) shows the 2D posterior contours of these two models are in relatively good agreement. A small shift towards lower amplitudes and higher spectral index can be seen when using the full HD ORF with both auto- and cross-correlations. Using the autocorrelation terms only, we find |$A_{\rm CP} = 5.1^{+6.7}_{-3.1} \times 10^{-15}$| and γ_CP = 3.9 ± 0.9, where the errors represent 95 per cent credible regions. Using the full HD ORF we can constrain the CP power law to |$A_{\rm CP} = 3.8^{+6.3}_{-2.5} \times 10^{-15}$| and γ_CP = 4.0 ± 0.9.

When we fix the power spectrum index to γ = 13/3 as shown in Fig. 3(b), it is clear that full HD model finds a systematically lower amplitude. In this case we find an amplitude of A_CP = 3.2 ± 1.0 × 10⁻¹⁵ for the autocorrelation only analysis and an amplitude of |$A_{\rm CP} = 2.8^{+1.2}_{-0.9} \times 10^{-15}$| using the full HD ORF, where the uncertainties represent the 95 per cent credible regions. These results are in broad agreement with published constraints on the CP. A more detailed comparison can be found in Section 4.3.

4.1.2 Split ORF analysis

Similar to how we may consider the autocorrelation parts of the ORF alone, the full ORF can be split into two independent processes. In this case the autocorrelation and the cross-correlation parts each have their own independent amplitude, as was done in Arzoumanian et al. (2020). In the HD ORF the autocorrelation part is Γ_aa = 1 and the cross-correlation parts are suppressed by at least a factor of 2, Γ_ab < 0.5. This makes the cross-correlations harder to constrain. Fig. 4 shows the posteriors for the two amplitudes of a split ORF analysis for fixed γ = 13/3, compared to the full auto + cross-correlation model. The cross-correlations do not have sufficient precision to place constraints on the amplitude of the GWB. However, they do place a 95 per cent upper limit of 3.6 × 10⁻¹⁵ on the GWB when the prior choice is taken into account. This is consistent with the amplitude derived using the full auto- and cross-correlation model in Fig. 4. The autocorrelation terms are much more informative. Combining the information from both shifts the amplitude towards lower values. This shows that the cross-terms can contribute to the full GWB search, even if they provide less information. The autocorrelations are more likely to be affected by intrinsic pulsar noise. Using a more sophisticated noise model for each pulsar can help produce a more robust estimate on the amplitude of any CP.

Figure 4.

The constraints from the split ORF analysis on the common process amplitude with spectral index fixed to γ = 13/3. The posterior from the autocorrelations is well constrained, while the posterior from the cross-correlations is unconstrained, but prefers amplitudes slightly smaller than the auto-only analysis. The effect of this is seen in the amplitude posterior from the full auto + cross-correlation model, where the posterior peaks at lower amplitude than the auto-only analysis.

4.1.3 Optimal statistic

Fig. 5 shows the amplitudes and signal-to-noise ratio (S/N) that are recovered by the pulsar noise marginalized optimal statistic (OS) method, which uses cross-correlations only. We find no evidence for a dipolar correlated process, as the amplitude and S/N for this model are centered on 0. SSE systematics are expected to manifest at specific frequencies related to the celestial bodies. The IPTA DR2 data set is long enough to probe lower frequencies that should be less affected by SSE errors (Vallisneri et al. 2020). The S/N |$= 0.6^{+1.2}_{-0.8}$| for the Hellings-Downs correlation is insufficient to claim a detection. This is consistent with the Bayesian model selection. The HD amplitude from the OS seems to be in tension with the Bayesian results for the autocorrelated CP, but consistent with the Bayesian results for the full HD model. This strengthens the case that the cross-terms have a significant role to play in parameter estimation as well as detection confidence. Finally, the OS has the largest S/N |$= 2.0^{+1.8}_{-1.4}$| for a monopole with a small amplitude. This can be due to the complexity of IPTA DR2 and some amount of unmodelled noise.

Figure 5.

Results from the noise marginalized optimal statistic. Top: Optimal statistic, A², distribution for monopole, dipole, and HD ORFs. The relevant posteriors from the Bayesian split ORF analysis are also shown for comparison. Bottom: The signal-to-noise (S/N) distribution for monopole, dipole, and HD ORFs.

As the spatial correlations are not well constrained, see Fig. 6, both the HD and monopolar correlation can fit the data. We have binned pulsar pairs according to their angular separation. Increasing the number of pulsars in the array as well as better timing of pulsars can help us to tighten the constraints.

Figure 6.

Cross-correlation ORF curve from the optimal statistic. The black points indicate the amount of cross-correlation for a given angular separation. Due to the large number of pulsar pairs, we have binned multiple pairs with similar angular separation. The blue and orange dashed lines show the best-fitting values for the HD and monopole correlations.

We test the significance of the OS S/N by performing two analyses which estimate the false alarm rate for a given S/N. The so-called phase-shifts and sky scrambles (Cornish & Sampson 2016; Taylor et al. 2017) break the correlations between the pulsars leaving the red noise power in the pulsars, but removing evidence for spatial correlations. By analysing the phase-shifted and sky scrambled data we can determine the rate of observing a particular S/N for a type of correlations in data that has none. These two false alarm studies result in p values, Table 2, too high to conclude that there is evidence for any correlations. It is possible that the measured HD S/N can therefore arise by chance from a common process with no spatial correlation.

Table 2.

The p-values calculated from various false alarm analyses of the data set. The measured values of the S/N are compared to the distribution of 10k analyses where the correlations are broken with phase shifts and sky scrambles. Since a monopolar spatial correlation is uniform across the sky, sky scrambles are unable to break the correlations. Hence only the phase shift p value is quoted.

	p, phase shift	p, sky scramble
HD ORF	0.25	0.24
Monopole	0.03	−

4.2 IPTA DR2 data subsets

With the large volume of the IPTA DR2 data set, we can look at different subsets to investigate hints of the origin and evolution of the CP signal across different slicings of the full data set.

4.2.1 Pulsar-based selection

Since a PTA is made from a number of single pulsars, we can look at how each pulsar contributes to the CP by itself. The dropout factor gives a measure on how consistent a given pulsar’s intrinsic red noise is with the CP by comparing a model with and without the CP for the pulsar (see e.g. Arzoumanian et al. 2020). The dropout factors for each pulsar computed using both the traditional hypermodel and factorized likelihood approaches are shown in Fig. 7. About 20 pulsars have factors >1, while only three slightly disfavour the CP, with the remaining pulsars displaying indifference.

Figure 7.

Individual pulsar consistency with common-spectrum process, error bars represent 95 per cent credible intervals. Pulsars with dropout factors >1 contribute to the detection of the CP. Dropout factors of ∼1 correspond to no evidence for or against the CP, usually due to higher white noise levels and/or shorter observation time-spans. Pulsars with dropout factors <1 are in tension with the CP.

Monte Carlo sampling uncertainties on the dropout factors (computed either way) can be estimated through statistical bootstrapping (Efron & Tibshirani 1994). In the hypermodel dropout analysis, the MCMC chain is re-sampled with replacement, generating a new statistical realization of the sampled chain that is exceptionally unlikely to be identical to the original chain. This process was repeated 10³ times, generating many realizations of the MCMC chain from which the dropout factors were computed. Hence a distribution of dropout factors over bootstrap realizations was generated for each pulsar, allowing us to compute median values and 95 per cent confidence intervals. A similar procedure was performed for the factorized likelihood approach. For a given bootstrap realization, each individual pulsar’s MCMC chain was re-sampled with replacement. With the re-sampled CP posteriors for each pulsar, the factorized-likelihood approach pieces together the dropout factor by iteratively removing pulsars from the array, and all by using bootstrapped pulsar CP posterior chains. This process was repeated for 10³ bootstrap realizations across 25 different combinations of meta-parameters used in the factorized-likelihood dropout factor calculation. The end result is that the median dropout factor and |$95{{\ \rm per\ cent}}$| confidence intervals were computed from a total of 2.5 × 10⁴ factorized-likelihood dropout values for each pulsar. As is seen in Fig. 7, all dropout factors are consistent between both techniques. The vast majority of pulsars have dropout factors with overlapping error bars from both methods, and those that do not are within a few sigma of each other. Those pulsars that show the largest disparity are ones for which there were MCMC sampling inefficiencies that manifested in different stages of the dropout factor calculation, e.g. in the Savage-Dickey density ratio, or in the integral of the (N−1) array’s CP likelihood over the posterior of a given pulsar (Taylor et al., in preparation).

The modularity and speed allowed by the factorized likelihood method can be used to approximate different combinations of pulsars within the array. These sub-arrays are a useful way to verify and understand the results we see in the full array. We created four sub-arrays, consisting of pulsars with the highest/lowest dropout factors, longest/shortest time-spans. These pulsars were selected by sorting all pulsars in the array by their dropout and time-span characteristics, and then taking the top half of 27 pulsars and the bottom half of 26 pulsars. The Savage-Dickey density ratio was calculated for these sub-arrays to compare them to that of the full array. The sub-array made up of the top half of pulsars according to dropout factor had a Savage-Dickey density ratio of 5.6 × 10⁹, an order of magnitude larger than that of the unaltered array, 1.6 × 10⁸. The corresponding sub-array of the bottom half of pulsars based on dropout had a density ratio of 1.8. When comparing the pulsars with the longest and shortest time-spans, the sub-array with the longer time-spans had a Savage-Dickey density ratio of 1.8 × 10⁷ and the shorter time-span sub-array had a density ratio of 3.6. These results were not surprising, as the dropout factor is a method of measuring the evidence for the array’s common process in a particular pulsar, so by removing those that have low (high) dropout factors, the evidence for the CP will increase (decrease). Similarly, pulsars with short time-spans are not sensitive to the lowest frequencies explored by the array when in combination with longer time-span pulsars.

4.2.2 Splitting IPTA DR2 by time

To test the evolution of the common-spectrum process, the DR2 is split into two data sets that have equally long time-spans, i.e. cutting the DR2 in two time slices (in a similar manner as Hazboun et al. 2020a). The two data sets are not fully equivalent though: the early part contains only 19 pulsars and is mostly dominated by single-radio frequency observations, while the second part has data from all 53 pulsars as well as multiradio frequency coverage and higher quality timing measurements. Each data set is then analysed separately. We find that the first half gives little information to the CP, with a broad power law 2D posterior contour that still encompasses the contour from the full data set. The second half contains the majority of information and produces almost identical constraints as the full data set. This is the expected evolution as the quality of the data set gradually improves over time (Hazboun et al. 2020a).

4.2.3 Constituent data sets

We can also select the data that were provided by the constituent PTA collaborations to get three data subsets: EPTA, NANOGrav, and PPTA. As each data subset has a different time-span, we set a frequency cutoff at 1.4 × 10⁻⁸Hz to limit the number of frequencies for the analyses. Fig. 8 shows that IPTA DR2 produces the tightest constraints on the CP power law compared to the constituent data sets. While the PPTA data is still consistent with a upper limit, some support for a common red noise can be found with the EPTA and NANOGrav data. The free spectra also show consistency with a power-law model that spans across all three constituent data sets and IPTA DR2.

Figure 8.

Comparison of IPTA DR2 to constituent data sets. Left-hand panel: Free spectral common-spectrum process model. The different sampling frequencies are a result of the constituent subsets covering different time-spans. In all cases the lowest frequency is the inverse observation time. Right-hand panel: 2D posterior for CP parameters log-amplitude and spectral index, where the contours represent 1–, 2–, and 3–σ confidence intervals. The combined IPTA DR2 data set is more constraining that its parts.

4.3 Comparison with other recent data sets

Since the data from the regional PTAs were combined to form the IPTA DR2, the regional PTAs have continued to collect data and improve their data analysis methodology. We can compare the results using the older IPTA DR2 data set and the most recent data sets from PPTA, NANOGrav, and EPTA. Compared to the constituent PTA data sets, the PPTA DR2 expands by ∼ 3 yr and 7 pulsars, the recent NANOGrav data set includes ∼ 4 more years and 10 new pulsars, the EPTA DR2 adds ∼ 7 yr for 6 pulsars (Kerr et al. 2020; Alam et al. 2021; Chen et al. 2021).

The published free spectral and power-law model recoveries can be found in Fig. 9. For simplicity, we also show the recovered amplitudes at the reference frequency of 1/(1yr) and fixed γ_CP = 13/3 in Fig. 10. The Mahalanobis distance D_M acts as a generalization to compute the n-dimensional sigma deviation between two distributions (Mahalanobis 1936),

$$\begin{eqnarray*} D_M = \sqrt{(\vec{\mu _1} - \vec{\mu _2}) \Sigma ^{-1} (\vec{\mu _1} - \vec{\mu _2})} \, , \end{eqnarray*}$$

(6)

where |$\vec{\mu _{1}}$| and |$\vec{\mu _{2}}$| are the mean vectors of the multivariate distributions to be compared and Σ = Σ₁+Σ₂ is the joint covariance. To quantify the overlap and consistency of the power-law parameters as determined using each data set, the Mahalanobis distance between the 2D posterior distributions are computed in Table 3. Despite some differences the posteriors overlap better than 3-sigma for all pairs of distributions.

Figure 9.

Comparison of IPTA DR2 to other recent data sets. Left-hand panel: Free spectral common-spectrum process model. The inclusion of legacy data not used in recent PTA analyses allows IPTA DR2 to reach lower frequencies despite missing the most recently collected data. Right-hand panel: 2D posterior for CP parameters log-amplitude and spectral index, where the contours represent the 1–, 2–, and 3–σ confidence intervals. All recent data sets are in broad agreement on the characteristics of a common-spectrum process.

Figure 10.

CP amplitude posteriors for fixed spectral index, γ = 13/3. IPTA DR2 and EPTA DR2 find a systematically higher amplitude for the common-spectrum process than NANOGrav 12.5 yr and PPTA DR2, although the disagreement is not substantial.

Table 3.

Mahalanobis distance between CP parameters (log-amplitude and spectral index) for each pair of PTAs. For all cases, there is less than 3-sigma separation.

	EPTA	PPTA	NANOGrav
IPTA	0.6	2.6	2.6
EPTA	–	2.3	2.4
PPTA	–	–	1.4

IPTA DR2, using older observations, still shows similar features as the NANOGrav 12.5, 6-pulsar EPTA DR2 and PPTA DR2 analyses, which have added a significant amount of new data to the regional PTA data sets. A future combination of these data sets will boost the total PTA sensitivity in the same way IPTA DR2 is more sensitive than its constituent data sets. Future combined IPTA data sets will be important for investigating the origin of this common-spectrum process.

5 DISCUSSION AND OUTLOOK

5.1 Source of the common-spectrum process

The first IPTA data release did not show signs of a common-spectrum temporally correlated process, but set an upper limit of 1.7 × 10⁻¹⁵ instead. This appears to be in tension with our results from analysis of the second data release with a CP amplitude of 2.8 × 10⁻¹⁵. However, there are two major differences to point out: (1) the different choice of priors for the pulsar red, DM, and common noise (Hazboun et al. 2020b) and (2) the DR1 upper limit was computed without the use of an SSE uncertainty model (Vallisneri et al. 2020). Both of which have been shown to lead to an increase in the upper limit, alleviating tensions between the DR1 and DR2 CP amplitudes.

As in other recent PTA analyses, we find strong evidence in favour of the CP over the noise only hypothesis. It is important to note that (1) the lack of support for GW-like spatial correlations prohibits any claims of GW detection, however (2) this type of evidence for a similar red noise is expected to precede a detection of spatial correlations (Siemens et al. 2013; Pol et al. 2021; Romano et al. 2021).

Goncharov et al. (2021) recently demonstrated that the common-spectrum process model is favoured over the noise-only hypothesis when the noise spectra cluster in a similar range, and it is not favoured anymore when the noise spectra are drawn from the prior distribution. Because we know that the employed prior distribution for red noise parameters is not representative, it is possible that the evidence we find for a common-spectrum process is caused by a rejection of a null hypothesis rather than by all pulsars exhibiting the spatially uncorrelated component of a GWB.

Thus, it is important to examine the single pulsar red noise in detail. We have looked at constraints on the simple power-law models for the pulsars used in the CP search. In general, pulsars with detectable intrinsic noise have comparable or larger noise than the CP; pulsars without red noise typically have large amount of white noise, such that the CP is ‘hidden‘. One noticeable exception is PSR J2317+1439, whose noise spectrum falls clearly below the CP, see also its low dropout factor in Fig. 7.

As the search for the common spectrum can be influenced by pulsar intrinsic noise, especially in an inhomogeneous data set, the crucial analysis has to consider information from the cross-correlations. It should be noted that the median amplitudes are slightly different in the analyses with and without spatial correlations, 2.8 × 10⁻¹⁵ versus 3.2 × 10⁻¹⁵. One can also note the stark difference in the posterior for the split ORF analysis, Fig. 4 and the optimal statistic analyses versus the Bayesian uncorrelated analysis, Fig. 5. In other analyses (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov et al. 2021) the amplitudes between the two analyses are more in line with one another. The difference here could be in part due to the very long baselines of only a handful of pulsars. This legacy data allows only scant opportunity for correlations amongst those few pulsars, while the long baselines allow the detection of autocorrelated power, even in noisier data. Another possibility is that there is unaccounted for noise in individual pulsars that is contaminating the signal. Advanced models to take these pulsar noises into account for the GWB search have been shown to affect the individual pulsar red noise (e.g. Arzoumanian et al. 2020; Goncharov et al. 2021; Chalumeau et al., 2022).

5.2 Other correlated signals

Spatial correlations in pulsar timing data have been studied in depth in the literature (Tiburzi et al. 2016; Roebber 2019), and their consideration is an important part of any GW detection procedure. While GWs induce a quadrupole-dominated set of correlations there are other types of spatial correlations between pulsar data sets (Tiburzi et al. 2016; Roebber & Holder 2017; Roebber 2019). Monopolar spatial correlations, i.e. all pulsars seeing the same shifts in residuals irrespective of sky position, can manifest from clock errors, either in the BIPM clock standards, or the various observatory clocks used across the world (Hobbs et al. 2020). Dipolar spatial correlations can manifest from the error in measurement of processes where the motion of the Earth in the Solar system is important (Caballero et al. 2018). This is most direct in the modelling of the Solar system barycenter frame of reference, into which all pulsar TOAs are transformed. Errors in solar wind modelling can also add dipolar correlations (Tiburzi et al. 2016). While monopole and dipole correlations are theoretically orthogonal to HD correlations in the space of overlap reduction functions, real data with noise can result in some of these modes mixing (Roebber 2019). This mixing could erroneously be detected as a GWB.

The polarization content of the GWB could also deviate from the two tensor transverse (TT) modes predicted by general relativity which lead to the HD spatial correlations (Arzoumanian et al. 2021b). Deviations from general relativity would result in a correlation pattern that differs from HD. We would like to emphasize that the current data set does not allow us to draw any conclusions on the presence of spatial correlations. As can be seen in Fig. 6 the uncertainties on the spatial correlation coefficients Γ_ab(ξ) determined by the optimal statistic analysis are large. For most angular bins the correlation is indistinguishable from zero, corresponding to the uncorrelated CP. Close to submission of this paper, we noticed a preprint by Chen, Wu & Huang (2021). They have analysed the IPTA DR2 searching for alternative GW polarizations and claim evidence for spatial correlations induced by a scalar transverse (ST) polarization mode. Chen, Wu & Huang (2021) also report a Bayes factor in favour of HD correlations (TT modes) compared to an uncorrelated CP about six times larger than we find: 12 against our 2. Even though we have not searched for the ST mode, we would like to highlight that the reported high evidence of spatial correlations in Chen, Wu & Huang (2021) is contrary to what we conclude using the same data set. The scalar transverse ORF is positive definite and should be accompanied by positive evidence for a monopolar correlation in analyses using only the cross-terms (Arzoumanian et al. 2021b). This is the case in our optimal statistic analysis where the monopolar correlations have the largest S/N and smallest false alarm p value. In this sense finding some evidence in favour of ST correlations is not too surprising, however, we find no conclusive evidence in favour of any correlation pattern. Using more information in our analysis with both auto- and cross-terms disfavours monopolar correlations compared to an uncorrelated CP. Therefore, the conclusions of Chen, Wu & Huang (2021) need to be taken with caution. Finally, Chen, Wu & Huang (2021) find a Bayes factor in favour of the common process to be several orders of magnitude smaller (⁠|$\log _{10}\rm {BF}=4.6$| compared to 8.2) than we do. They use a different pulsar noise model, including an additional sinusoidal annual DM variation in all pulsars, which could account for some, but likely not all, of the differences.

5.3 Astrophysical interpretation

The nHz GWB is generally thought to be dominated by GW emission from SMBHBs (Burke-Spolaor et al. 2019), with the most massive local SMBHBs expected to be individually observable in the next 5-10 yr (Mingarelli et al. 2017). Given that these systems are just a local subset of the cosmological SMBHB population producing the GWB, their local number density, Φ_{BHB, 0}, should correlate with the amplitude of the GWB. The GWB amplitude induced by a cosmological population of circularized SMBHBs can be expressed in geometric units (where G = c = 1) as (e.g. Phinney 2001; Sesana et al. 2008; Sesana 2013b)

$$\begin{eqnarray*} A_{{\rm GWB}}^{2} = \frac{4}{3 \pi ^{1 / 3}} \iiint \, {\rm d}M_{1} {\rm d}z {\rm d}q \frac{\mathcal {M}^{5 / 3}}{\left(1 + z\right)^{1 / 3}} \frac{d^{3}{\Phi _{\rm {BHB}}}}{{\rm d}M_{1} {\rm d}z {\rm d}q}, \end{eqnarray*}$$

(7)

where M₁ is the mass of the primary SMBH, M₂ the mass of the secondary, q = M₂/M₁ ≤ 1 is the mass ratio, |$\mathcal {M}^{5 / 3} = [q (1 + q)^{-2}] M^{5 / 3}$| is the SMBHB chirp mass with total binary mass M = M₁+M₂, and d³Φ_BHB/(dM₁dzdq) is the differential comoving number density of SMBHBs per unit M₁, z, and q.

To determine the local number density of SMBHBs implied by a A_GWB ≈ 2.8 × 10⁻¹⁵ background we use the quasar-based SMBHB model from Casey-Clyde et al. (2021), which assumes proportionality between SMBHB and quasar populations (which may be triggered by galaxy major mergers; Stemo et al. 2021) over mass and redshift. This has the effect of setting |$A_{{\rm GWB}} \propto \sqrt{\Phi _{{\rm BHB}, 0}}$|⁠, so that A_GWB directly implies Φ_{BHB, 0}. To check coverage of the entire signal from SMBHBs over mass and redshift, we parametrize A_GWB = A_GWB(Φ_{BHB, 0}, M_{1, min}, z_max), where M_{1, min} and z_max are the minimum primary SMBH mass and maximum redshift, respectively, in equation (7) (Casey-Clyde et al. 2021). We plot this parametrization of the GWB compared to various strain measurements in Fig. 11, including an A_GWB ≈ 2.8 × 10⁻¹⁵ signal (bold contour, grey isosurface). The 2D panels of Fig. 11 show three representative slices from this parameter space (one along each axis), with contours denoting their intersection with isosurfaces of constant GWB signal amplitude. The 3D plot shows the bottom right 2D panel in this 3D parameter space, along with its intersection with an A_GWB ≈ 2.8 × 10⁻¹⁵ isosurface. We find that recovery of a background amplitude like A_CP requires |$\Phi _{\rm {BHB}, 0} \approx 1.5 \times 10^{-5}\,\,\rm {Mpc}^{-3}$| (corresponding to the bottom right 2D panel in Fig. 11), roughly an order of magnitude larger than the |$\sim 1.6 \times 10^{-6}\,\,\rm {Mpc}^{-3}$| number density implied by Mingarelli et al. (2017).

Figure 11.

The GWB characteristic strain as a function of the local SMBHB number density, Φ_{BHB, 0}, and the minimum primary BH mass, |$M_{\rm {BH}, 1 \min }$|⁠, and maximum redshift, z_max, of the population contributing |$\gtrsim 95 {{\ \rm per\ cent}}$| of the GWB signal. Left-hand panel: Three representative slices of the strain in this parameter space (one along each axis), with solid contours showing their intersection with isosurfaces of constant strain value (A_CP shown in bold). Left-hand panel: 3D visualization of the |$z_{\rm {max}} - M_{\rm {BH}, 1, \min }$| panel from the left and its intersection with an A_GWB = 2.8 × 10⁻¹⁵ isosurface (grey).

Besides this new quasar-based method the standard approach to determining the local number density Φ_{BHB, 0} is to model d³Φ_BHB/(dM₁dzdq) using major mergers and empirically observed galaxy and black hole relations (e.g. Simon & Burke-Spolaor 2016; Chen et al. 2019; Middleton et al. 2021). Following the methods from Middleton et al. (2021) we analyse the IPTA DR2 CP amplitude. Fig. 12 compares the spread of amplitude and spectral index from the IPTA DR2 CP against values recovered from realizations of SMBHB population simulations. The original population constraints using the NANOGrav (Arzoumanian et al. 2020) frequency bins are shown by the grey shaded area. Repeating the analysis with the frequency coverage of the IPTA DR2 gives the purple shaded contours. As we reach into lower frequencies the simulations become more constrained towards the expected spectral index γ = 13/3. Limiting the SMBHB chirp mass |$\mathcal {M} \gt 10^{8.5} M_\odot$| in the integral of equation (7) we get |$\Phi _{\rm {BHB}, 0} \approx 3.0 \times 10^{-5}\,\,\rm {Mpc}^{-3}$|⁠, which is about a factor of 20 times larger than the number density from Mingarelli et al. (2017).

Figure 12.

Comparison of power-law constraints versus theoretical SMBHB populations. The 2D amplitude, spectral index constraints of the CPs from Fig. 9 are compared to the region of parameters recovered from a large number of realizations of SMBHB population simulations using astrophysical relations from (M21 Pop, Middleton et al. 2021) shown in grey contours and this work (IPTA DR2 Pop) shown in purples contours.

5.4 Outlook for other GWBs

Although we have evidence for a single common process, whose amplitude and spectral index are consistent with predicted values from a population of SMBHBs with a number density Φ₀ ≈ 10⁻⁵Mpc⁻³, other sources could also be plausible, for example, primordial black holes (e.g. De Luca, Franciolini & Riotto 2021; Kohri & Terada 2021; Vaskonen & Veermäe 2021), cosmic strings (e.g. Blanco-Pillado, Olum & Wachter 2021; Blasi, Brdar & Schmitz 2021; Ellis & Lewicki 2021) or phase transitions (e.g. Arzoumanian et al. 2021a, and references therein). These sources could produce a GWB consistent with the CP contours from PTA data.

Pol et al. (2021) have shown that the initial confident detection of a GWB including HD correlation will place very stringent constraints on the properties of the possible source of GWBs. It is also possible to have several backgrounds affecting the data, splitting the total common power into several components. A detailed study into how well we can separate multicomponent GWBs is underway.

5.5 Modern noise mitigation

The PTA community continues to develop new data analysis strategies towards the detection of gravitational waves in pulsar timing data. Aggarwal et al. (2019, 2020) showed that unmodelled noise features in a single pulsar could leak into the gravitational wave channel for deterministic continuous GW and GW memory signals, respectively. More recently immense effort was put into the development of individualized noise models for PPTA pulsars, demonstrating that sensitivity can be gained from better modelling (Goncharov et al. 2021). Similar advanced noise modelling efforts are currently underway in both NANOGrav (Arzoumanian et al., in preparation) and EPTA (Chalumeau et al., 2022). More sophisticated noise modelling is important, because many types of noise can add steep spectral index, low frequency power to pulsar data sets, complicating GWB recovery. For example, noise from fluctuations due to the interstellar medium or time-correlated noise from long-term instrumental effects in telescope systems, which could arise from a combination of popularization miscalibration and secular changes in receiver gain.

6 CONCLUSION

This work shows the immense potential of combining the global efforts of PTA collaborations into one data set. Figs 1 and 8 show that IPTA DR2 is significantly more sensitive than any of the constituent data sets from which it is constructed. While the data in DR2 have now been superseded by more up-to-date efforts, (Kerr et al. 2020; Alam et al. 2021; Chen et al. 2021), Fig. 9 shows that the sensitivity from combining these older data sets is comparable with these newer single PTA data releases.

The conclusions of this analysis are broadly similar to the various GWB analyses carried out by the NANOGrav, the EPTA and PPTA (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov et al. 2021). All of these data sets favour the CP model over one with only intrinsic red noise in the individual pulsars. None of these data sets shows clear support for the spatial correlations indicative of a GWB. Therefore a detection of a GWB cannot be claimed. The strong detection of red noise that broadly matches the spectral characteristic of a GWB from SMBHBs before there is support for spatial correlations is expected from our understanding of the change in sensitivity of PTAs (Siemens et al. 2013; Romano et al. 2021). As was shown in Pol et al. (2021), if the power in the autocorrelations of these pulsars is the first sign of the GWB then evidence for the spatial correlations should follow in upcoming data sets. If the next individual PTA data sets show increased support, but are short of detection thresholds, then combining them into an IPTA data set could immediately result in a data set with a significant detection. Such a combination will have the longest time-span, largest number of pulsars and independent observing systems, and will thus enable a robust GWB search.

ACKNOWLEDGEMENTS

The International Pulsar Timing Array (IPTA) is a consortium of existing Pulsar Timing Array collaborations, namely, the European Pulsar Timing Array (EPTA), North American Nanohertz Observatory for Gravitational Waves (NANOGrav), Parkes Pulsar Timing Array (PPTA), and the recent addition of the Indian Pulsar Timing Array (InPTA). Observing collaborations from China and South Africa are also part of the IPTA.

The EPTA is a collaboration between European and partner institutes with the aim to provide high precision pulsar timing to work towards the direct detection of low-frequency gravitational waves. An Advanced Grant of the European Research Council to implement the Large European Array for Pulsars (LEAP) also provides funding. Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a Consolidated Grant (ST/T000414/1) from the UK’s Science and Technology Facilities Council. The Nançay radio Observatory is operated by the Paris Observatory, associated to the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from ‘Programme National de Cosmologie and Galaxies’ (PNCG), and ‘Programme National Hautes Energies’ (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). The Sardinia Radio Telescope (SRT) is funded by the Department of University and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as National Facility by the National Institute for Astrophysics (INAF).

The NANOGrav Physics Frontiers Center is supported by NSF award number 1430284. The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Arecibo Observatory is a facility of the National Science Foundation operated under cooperative agreement by the University of Central Florida in alliance with Yang Enterprises, Inc. and Universidad Metropolitana.

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.

JA acknowleges support by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of ‘Science and Society’ Action Always strive for excellence – ‘Theodoros Papazoglou’ (Project Number: 01431). SBS acknowledges generous support by the NSF through grant AST-1815664. The work is supported by National SKA program of China 2020SKA0120100, Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. JACC was supported in part by NASA CT Space Grant PTE Federal Award Number 80NSSC20M0129. CMFM and JACC are also supported by the National Science Foundation’s NANOGrav Physics Frontier Center, Award Number 2020265. AC acknowledges support from the Paris Île-de-France Region. Support for HTC was provided by NASA through the NASA Hubble Fellowship Program grant HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. SD is the recipient of an Australian Research Council Discovery Early Career Award (DE210101738) funded by the Australian Government. GD, RK and MKr acknowledge support from European Research Council (ERC) Synergy Grant ‘BlackHoleCam’ Grant Agreement Number 610058 and ERC Advanced Grant ‘LEAP’ Grant Agreement Number 337062. TD is supported by the NSF AAG award number 2009468. ECF is supported by NASA under award number 80GSFC17M0002.002. BG is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN. Portions of this work performed at the Naval Research Laboratory is supported by NASA and ONR 6.1 basic research funding. MTL acknowledges support received from NSF AAG award number 200968. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. JWMK is a CITA Postdoctoral Fellow: This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), (funding reference CITA 490888-16). ASa, ASe, and GS acknowledge financial support provided under the European Union’s H2020 ERC Consolidator Grant ‘Binary Massive Black Hole Astrophysic’ (B Massive, Grant Agreement: 818691). RMS acknowledges support through Australian Research Council Future Fellowship FT190100155. JS acknowledges support from the JPL R&TD program. This research was funded partially by the Australian Government through the Australian Research Council (ARC), grants CE170100004 (OzGrav) and FL150100148. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. SRT acknowledges support from NSF grants AST-2007993 and PHY-2020265. SRT also acknowledges support from the Vanderbilt University College of Arts & Science Dean’s Faculty Fellowship program. AV acknowledges the support of the Royal Society and Wolfson Foundation. JPWV acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039).

Author Contributions: An alphabetical-order author list was used for this paper to recognize the large amount of work and efforts contributed by many people within the IPTA consortium. All authors have contributed to the work via data creation, combination or analysis or contributions to the paper. The data set was created by a team led by BBPP, MED, PBD, MKe, LL, DJN, SO, SMR, and MJK, by combining data sets from the constituent PTAs (EPTA, NANOGrav, and PPTA). The data analysis has been the task of the Gravitational Wave Analysis Working Group (WG) of the IPTA. Many analyses have been performed over more than 2 yr by PTB, AC, SChe, JAE, JMG, JSH, KI, ADJ, WGL, CMFM, NSP, NKP, DJR, LSc, JS, LSp, SRT, and SJV under the leadership of several WG chairs over this time period: PTB, SChe, JSH, PDL, CMFM, DJR, and SRT.

Initial tests, checks of the data quality and exploratory analyses were made by PTB, SChe, JAE, JMG, JSH, KI, CMFM, DJR, JS, SRT, and SJV. Although not much of the work is presented in this paper, it helped us tremendously in understanding the subtle details of the data set and built confidence in our analysis. The final results shown come from runs performed by AC, SChe, JSH, ADJ, WGL, NSP, NKP, ASa, LSc, GS, LSp, and SRT. The single pulsar noise analysis was done by SChe. The Bayesian parameter estimation were done by NSP, JSH, and WGL. The optimal statistic analysis was run by JSH and NSP, with help from AC, NKP, and LSp for the phase-shifts and sky-scramble null-distributions. ADJ and JSH have contributed to the Bayes factor computations. The factorized likelihood method has been applied to the data by LSc and SRT. NSP has computed the dropout factors. SChe performed the analyses of the constituent data sets, with cross-checks from ASa and GS. PTB lead the comparison between the IPTA DR2 and recent PTA data sets, which were graciously provided by JS for the NANOGrav 12.5 yr, BG for the PPTA DR2 and SChe for the EPTA DR2. The astrophysical interpretation was lead by JACC and CMFM with contributions from SChe and ASe.

SChe lead the coordination of the paper writing with the help of JSH and PTB. PTB, JACC, SChe, MED, BG, JSH, CMFM, BBPP, NSP, LSc, RMS, SRT, and SJV contributed to the writing of the paper. The figures were created by NSP, JSH, JACC, and SChe.

DATA AVAILABILITY

The timing data used can be found on https://gitlab.com/IPTA/DR2. Selected chain files for the analyses presented in this paper can be found on https://ipta4gw.org and 10.5281/zenodo.5787557.

Author notes

NASA Hubble Fellowship Program Einstein Postdoctoral Fellow.