1. Introduction
After the direct detection of gravitational waves (GWs) from a binary black hole [1] and a binary neutron star [2] mergers by LIGO-Virgo, other types of GW sources are yet to be identified. Of particular interest is the stochastic gravitational-wave background (SGWB) produced by the superposition of a large number of independent GW signals from compact binary coalescences. While the ground-based interferometers are sensitive to GWs from Hz to kHz, a pulsar timing array (PTA) [3–5], which regularly monitors the time of arrivals (TOAs) of radio pulses from an array of stable millisecond pulsars, offers a unique and powerful probe to correlated signals at low frequencies from nHz to μHz. The SGWB sources for PTAs could come from the inspiral of supermassive black hole binaries (SMBHBs) [6–8], the first-order phase transition [9, 10], and the scalar-induced GWs [11–14], etc. It is expected that the SGWB from SMBHBs will be the first GW signal to be detected with PTAs [15].
There are three major PTAs with accumulated pulsar-timing data of more than a decade, namely the European Pulsar Timing Array (EPTA) [16], the North American Nanoherz Observatory for Gravitational Waves (NANOGrav) [17], and the Parkes Pulsar Timing Array (PPTA) [18]. These collaborations support the International Pulsar Timing Array (IPTA) [19, 20]. Over the last decades, PTAs have accumulated increasingly sensitive data sets, and the null-detection of GWs with PTAs has successfully constrained various astrophysical scenarios, such as cosmic strings [21–24], SGWBs from SMBHBs with power-law spectra [21, 22, 25], and primordial black holes [26], etc. It is widely expected that the inkling of an SGWB will first manifest as the emergence of a spatially uncorrelated common-spectrum process (UCP) among all pulsars, and culminate in the appearance of the spatial correlations that unambiguously signify the detection of an SGWB.
In a recent analysis, the NANOGrav collaboration found strong evidence for a stochastic common-spectrum process modeled by a power-law spectrum in their 12.5 year data set [27]. However, there was no statistically significant evidence for the tensor transverse (TT) spatial correlations which are deemed to be necessary to claim an SGWB detection consistent with general relativity. Authors in [28] then reanalyzed the NANOGrav 12.5 year data set and found strong Bayesian evidence that the common-spectrum process reported by NANOGrav collaboration has the scalar transverse (ST) spatial correlations which can originate from a general metric theory of gravity. Later on, the PPTA collaboration analyzed their second data release (DR2) and also found a common-spectrum process but without significant evidence for, or against, the TT spatial correlations [29]. Furthermore, we [30] searched for the non-tensorial polarizations in the PPTA DR2, and found no significant evidence supporting the existence of alternative polarizations, thus constraining the amplitude of each polarization mode. Note that the 95% upper limit on the amplitude of the ST mode from [30] is consistent with the result from [28].
The IPTA collaboration has published their second data release comprising of 65 pulsars in total with a timespan as long as about three decades [31]. The IPTA DR2 is an invaluable complement to the NANOGrav 12.5 year data set and PPTA DR2 in the sense of more pulsars included and a longer timespan. In this letter, we aim to search for the SGWB signal modeled as a power-law spectrum in the IPTA DR2 with a particular interest in exploring if the ST spatial correlations are present in the data set or not.
2. The data set and methodology
The IPTA DR2 [31] was created by combining published data from individual PTA data releases, including EPTA DR1 [32], NANOGrav 9 year data set [33], and PPTA DR1 [18]. It consists of 65 pulsars and provides a better sky coverage compared to the IPTA DR1 [34]. There are two data combination versions in the IPTA DR2, namely VersionA and VersionB. The two versions are different in modeling the dispersion measure (DM) variation and handling the noise properties of pulsars [31]. In particular, the noise parameters of the pulsars in VersionB are re-estimated based on the IPTA data combination, while VersionA uses the previously constrained values from other PTA data sets [31]. It has been shown that the overall time-dependent DM variations modeled by these two methods are largely consistent with each other [31]. In this work, we use the VersionB data release by choosing only pulsars with a timing baseline greater than three years in our analyses and excluding the pulsar J1939+2134 because of its complicated DM variation and timing noise [18, 35, 36]. Therefore, all results in this work are based on 52 pulsars that meet the requirements.
The GWs will manifest as the unexplained residuals in the pulsar TOAs after subtracting a deterministic timing model that accounts for the pulsar spin behavior and the geometric effects due to the motion of the pulsar and the Earth [3, 5]. For two pulsars a and b, the cross-power spectral density of the timing residuals induced by an SGWB at frequency f is [37–39]
$\begin{eqnarray}{S}_{{ab}}^{P}(f)=\displaystyle \frac{{h}_{c,P}^{2}}{12{\pi }^{2}{f}^{3}}{{\rm{\Gamma }}}_{{ab}}^{P}(f),\end{eqnarray}$
where ${h}_{c}^{P}(f)$ is the characteristic strain of the polarization mode P. In this work, we only search for the TT and ST polarization modes. The TT mode is the usual TT mode containing the ‘+' and ‘×' polarizations, and the ST mode is the spin-0 breathing polarization. The overlap functions Γab(f) for the TT and ST polarizations can be approximated by [37, 40] $\begin{eqnarray}{{\rm{\Gamma }}}_{{ab}}^{\mathrm{TT}}(f)=\displaystyle \frac{1}{2}(1+{\delta }_{{ab}})+\displaystyle \frac{3}{2}{k}_{{ab}}\left(\mathrm{ln}{k}_{{ab}}-\tfrac{1}{6}\right),\end{eqnarray}$
$\begin{eqnarray}{{\rm{\Gamma }}}_{{ab}}^{\mathrm{ST}}(f)=\displaystyle \frac{1}{8}\left(3+4{\delta }_{{ab}}+\cos {\xi }_{{ab}}\right),\end{eqnarray}$
where δab is the Kronecker delta symbol, ξab is the angle between pulsars a and b, and ${k}_{{ab}}\equiv (1-\cos {\xi }_{{ab}})/2$. The TT overlap function ${{\rm{\Gamma }}}_{{ab}}^{\mathrm{TT}}$ is also known as the Hellings and Downs [40] or quadrupolar correlations. Note that the TT and ST overlap functions are related to each other by ${{\rm{\Gamma }}}_{{ab}}^{\mathrm{TT}}\,={{\rm{\Gamma }}}_{{ab}}^{\mathrm{ST}}+(3/2){k}_{{ab}}\mathrm{ln}{k}_{{ab}}$. Here, we also consider a parameterized overlap function $\begin{eqnarray}{{\rm{\Gamma }}}_{{ab}}(f)=\displaystyle \frac{1}{8}\left(3+4{\delta }_{{ab}}+\cos {\xi }_{{ab}}\right)+\displaystyle \frac{\alpha }{2}{k}_{{ab}}\mathrm{ln}{k}_{{ab}},\end{eqnarray}$
which reduces to ${{\rm{\Gamma }}}_{{ab}}^{\mathrm{ST}}$ when α = 0, and ${{\rm{\Gamma }}}_{{ab}}^{\mathrm{TT}}$ when α = 3.We take the characteristic strain ${h}_{c}^{P}(f)$ to be a power-law form, which can originate from the SGWB produced by a population of inspiraling SMBHBs [6–8]. Assuming the binaries are in circular orbits and the orbital decay is dominated by the GW emission, the cross-power spectral density can be approximately estimated by [41]
$\begin{eqnarray}{S}_{{ab}}^{P}(f)={{\rm{\Gamma }}}_{{ab}}^{P}\displaystyle \frac{{{ \mathcal A }}_{P}^{2}}{12{\pi }^{2}}{\left(\tfrac{f}{{f}_{\mathrm{yr}}}\right)}^{-{\gamma }_{P}}{f}_{\mathrm{yr}}^{-3},\end{eqnarray}$
where ${{ \mathcal A }}_{P}$ is the GW amplitude of the polarization mode P, and fyr = 1 year−1. The power-law index γP is 13/3 for the TT polarization, and 5 for the ST polarization. The dimensionless GW energy density parameter per logarithm frequency normalized by the critical energy density for the polarization mode P is related to ${{ \mathcal A }}_{P}$ by [42] $\begin{eqnarray}{{\rm{\Omega }}}_{\mathrm{GW}}^{P}(f)=\displaystyle \frac{2{\pi }^{2}}{3{H}_{0}^{2}}{f}^{2}{h}_{c,P}^{2}=\displaystyle \frac{2{\pi }^{2}{f}_{\mathrm{yr}}^{2}}{3{H}_{0}^{2}}{{ \mathcal A }}_{P}^{2}{\left(\tfrac{f}{{f}_{\mathrm{yr}}}\right)}^{5-{\gamma }_{P}},\end{eqnarray}$
where ${H}_{0}=67.4\,{\mathrm{kmsec}}^{-1}{\,\mathrm{Mpc}}^{-1}$ is the Hubble constant taken from Planck 2018 [43].We now briefly describe the noise model used in our analyses. After subtracting the timing model from the TOAs, the timing residuals δt of each single pulsar are contributed from a number of sources by (see e.g. [36])TEMPO2 [45, 46] through libstempo 6(6https://vallis.github.io/libstempo) interface, and ε is a vector denoting small offsets for the parameters of timing model. The second term $\delta {{\boldsymbol{t}}}_{\mathrm{SN}}$ is the stochastic contribution from the red spin noise (SN) intrinsic to each pulsar and is modeled by a power law with 30 frequency components. The third term δtDM denotes the stochastic contribution from the DM noise which is also modeled by a power law with 30 frequency components. Unlike the SN, the DM noise is dependent upon the radio frequency whose information is added to the Fourier basis components. The fourth term δtyrDM is the stochastic contribution caused by annual DM variation described by a deterministic yearly sinusoid (see e.g. [36]). The fifth term δtWN represents the stochastic contribution due to the white noise (WN), including a scale parameter on the TOA uncertainties (EFAC), an added variance (EQUAD), and a per-epoch variance (ECORR) for each backend/receiver system (see e.g. [47]). We include separate EFACs and EQUADs for all the backend/receiver-dependent PTA data sets, and separate ECORRs for the backend/receiver-dependent NANOGrav data sets. The last term δtCP is the stochastic contribution due to the common-spectrum process (such as an SGWB) with the cross-power spectral density given by equation (5 ). In the analyses, we use 10 frequency components roughly starting from 1.09 × 10−9 Hz to 1.09 × 10−8 Hz for the common-spectrum processes. For pulsar J1713+0747, we also include a chromatic exponential dip to model the sudden change in dispersion when the signal passes through the interstellar medium during propagation [36].
$\begin{eqnarray}\delta {\boldsymbol{t}}=M{\boldsymbol{\epsilon }}+\delta {{\boldsymbol{t}}}_{\mathrm{SN}}+\delta {{\boldsymbol{t}}}_{\mathrm{DM}}+\delta {{\boldsymbol{t}}}_{\mathrm{yrDM}}+\delta {{\boldsymbol{t}}}_{\mathrm{WN}}+\delta {{\boldsymbol{t}}}_{\mathrm{CP}}.\end{eqnarray}$
The first term Mε accounts for the inaccuracies in the subtraction of the timing model (see e.g. [44]), in which M is the timing model design matrix obtained from We use the latest JPL solar system ephemeris (SSE) DE438 [48] as the fiducial SSE as opposed to the DE436 [49] that was used to create the IPTA DR2. To extract information from the data, we perform similar Bayesian parameter inferences based on the methodology in [22, 27]. The model parameters and their prior distributions are summarized in table 1. We first perform the parameter estimations for each single pulsar without including the stochastic contribution from the common-spectrum process (i.e. the δtCP term in equation (7 )). To reduce the computational costs, we then fix the white noise parameters to their max likelihood values from single-pulsar analysis. We use enterprise [50] and enterprise _extension [51] software packages to calculate the likelihood and Bayes factors and use PTMCMCSampler [52] package to do the Markov chain Monte Carlo sampling. Similar to [27, 53], we use draws from empirical distributions to sample the parameters from SN, DM noise, and annual DM variation, with the distributions based on the posteriors obtained from the single-pulsar Bayesian analysis, thus reducing the number of samples needed for the chains to burn in.
Table 1. Parameters and their prior distributions used in the analyses. |
| Parameter | Description | Prior | Comments |
|---|---|---|---|
| Red noise | |||
| ${A}_{\mathrm{SN}}$ | Red-noise power-law amplitude | Log-Uniform [−20, −11] | One parameter per pulsar |
| ${\gamma }_{\mathrm{SN}}$ | Red-noise power-law spectral index | Uniform [0, 7] | One parameter per pulsar |
| DM noise | |||
| ADM | Red-noise power-law amplitude | Log-Uniform [−20, −11] | One parameter per pulsar |
| γDM | Red-noise power-law spectral index | Uniform [0, 7] | One parameter per pulsar |
| Annual DM variation | |||
| ${{ \mathcal A }}_{Y}$ | Annual DM variation amplitude | Log-Uniform [−10, −2] | One parameter per annual event |
| φY | Annual DM variation phase | Uniform [0, 2π] | One parameter per annual event |
| DM exponential dip | |||
| ${{ \mathcal A }}_{E}$ | Exponential dip amplitude | Log-Uniform [−10, −2] | One parameter for pulsar J1713+0747 |
| tE[MJD] | Time of the event | Uniform [54500, 55 000] | One parameter for pulsar J1713+0747 |
| τE[MJD] | Relaxation time for the dip | Log-Uniform [0, 2.5] | One parameter for pulsar J1713+0747 |
| White noise | |||
| Ek | EFAC per backend/receiver system | Uniform [0, 10] | Single-pulsar analysis only |
| Qk[s] | EQUAD per backend/receiver system | Log-Uniform [−8.5, −5] | Single-pulsar analysis only |
| Jk[s] | ECORR per backend/receiver system | Log-Uniform [ − 8.5, − 5] | Single-pulsar analysis only |
| Common-spectrum process | |||
| ${{ \mathcal A }}_{\mathrm{UCP}}$ | UCP power-law amplitude | Log-Uniform [−18, −14] | One parameter for PTA |
| Log-Uniform [−18, −11](γUCP varied) | One parameter for PTA | ||
| γUCP | UCP power-law spectral index | Delta function (γUCP = 13/3) | Fixed |
| Uniform [0, 7] (γUCP varied) | One parameter for PTA | ||
| ${{ \mathcal A }}_{\mathrm{TT}}$ | GW amplitude of TT polarization | Log-Uniform [−18, −14] | One parameter for PTA |
| ${{ \mathcal A }}_{\mathrm{ST}}$ | GW amplitude of ST polarization | Log-Uniform [−18, −14] | One parameter for PTA |
| α | Parameter in equation ( | Uniform [−10, 10] | One parameter for PTA |
3. Results and discussion
Our analyses are mainly based on the Bayesian inference in which the Bayes factor is used to quantify the model selection scores. The Bayes factor is defined as
$\begin{eqnarray}{ \mathcal B }{ \mathcal F }\equiv \displaystyle \frac{\Pr ({ \mathcal D }| {{ \mathcal M }}_{2})}{\Pr ({ \mathcal D }| {{ \mathcal M }}_{1})},\end{eqnarray}$
where $\Pr ({ \mathcal D }| { \mathcal M })$ denotes the probability that the data ${ \mathcal D }$ are produced under the assumption of model ${ \mathcal M }$. Model ${{ \mathcal M }}_{2}$ is preferred if the ${ \mathcal B }{ \mathcal F }$ is sufficiently large. An interpretation of the ${ \mathcal B }{ \mathcal F }$ in model comparison given by [54] can be found in table 2.Table 2. An interpretation of the Bayes factor in determining which model is favored, as given by [54]. |
| ${ \mathcal B }{ \mathcal F }$ | $\mathrm{ln}{ \mathcal B }{ \mathcal F }$ | Strength of evidence |
|---|---|---|
| < 1 | < 0 | Negative |
| 1 − 3 | 0 − 1 | Not worth more than a bare mention |
| 3 − 20 | 1 − 3 | Positive |
| 20 − 150 | 3 − 5 | Strong |
| > 150 | > 5 | Very strong |
In table 3, we summarize the Bayes factors of various models. The $\mathrm{ln}{ \mathcal B }{ \mathcal F }$ of the UCP model with γUCP = 13/3 versus the noise only (NO) model without any common-spectrum process is 10.5, indicating significant Bayesian evidence for a common-spectrum process in the IPTA DR2. When allowing the power-index to vary, the γUCP shows a relatively broad distribution as can be seen in figure 1. The median value and 90% equal-tailed credible intervals are ${\mathrm{log}}_{10}{{ \mathcal A }}_{\mathrm{UCP}}=-{14.26}_{-0.47}^{+0.41}$ and ${\gamma }_{\mathrm{UCP}}={3.76}_{-0.95}^{+0.96}$.
Figure 1. One and two-dimensional marginalized posteriors of the amplitude ${{ \mathcal A }}_{\mathrm{UCP}}$ and the power-index γUCP obtained from the UCP model with γUCP allowed to vary. We show both the 1σ and 2σ contours in the two-dimensional plot. |
Table 3. The $\mathrm{ln}{ \mathcal B }{ \mathcal F }$ between pairs of models. The digit in the parentheses represents the uncertainty on the last quoted digit. |
| UCP versus NO | TT versus UCP | ST versus UCP | ST+TT versus UCP |
|---|---|---|---|
| 10.5(9) | 2.53(3) | 3.39(4) | 3.63(4) |
From table 3, the $\mathrm{ln}{ \mathcal B }{ \mathcal F }$ of the ST model versus the UCP (with γUCP = 13/3) model is 3.39, indicating ‘strong' Bayesian evidence of the ST correlations in IPTA DR2 according to table 3. We obtain the amplitude as ${{ \mathcal A }}_{\mathrm{ST}}={1.29}_{-0.44}^{+0.51}\,\times {10}^{-15}$ or equivalently ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{ST}}={2.31}_{-1.30}^{+2.19}\times {10}^{-9}$, at frequency of 1 year−1. This result is consistent with the one reported in [28] where ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{ST}}={1.54}_{-0.71}^{+1.20}\times {10}^{-9}$ from the NANOGrav 12.5 year data set. On the other hand, the $\mathrm{ln}{ \mathcal B }{ \mathcal F }$ of the TT model versus the UCP (with γUCP = 13/3) model is 2.53, indicating ‘positive' Bayesian evidence of the TT correlations in IPTA DR2. As the evidence for the TT correlations is insignificant, we do not report the median value of the ${{ \mathcal A }}_{\mathrm{TT}}$ parameter here. Both the posteriors of the ${{ \mathcal A }}_{\mathrm{TT}}$ and ${{ \mathcal A }}_{\mathrm{ST}}$ parameters obtained respectively from the TT and ST models are shown in figure 2.
Figure 2. Marginalized posteriors of ${{ \mathcal A }}_{\mathrm{TT}}$ and ${{ \mathcal A }}_{\mathrm{ST}}$ parameters obtained from the TT and ST models, respectively. |
Furthermore, we consider an ST+TT model where both the ST and TT spatial correlations are taken into account simultaneously. The $\mathrm{ln}{ \mathcal B }{ \mathcal F }$ of the ST+TT model versus the UCP (with γUCP = 13/3) model is 3.63, indicating no significant evidence for a second common-spectrum process with TT correlations on top of the ST-correlated process. The posterior distributions of the ST and TT amplitudes obtained from the ST+TT model are shown in figure 3, with the amplitude of ST mode from this model is ${{ \mathcal A }}_{\mathrm{ST}}={1.11}_{-1.10}^{+0.63}\times {10}^{-15}$ or equivalently ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{ST}}\,={1.71}_{-1.70}^{+2.48}\times {10}^{-9}$, at frequency of 1/year. Note that this result is consistent with the one obtained from the ST model and that reported in [28] where ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{ST}}={1.54}_{-0.71}^{+1.20}\times {10}^{-9}$ from the NANOGrav 12.5 year data set.
Figure 3. One and two-dimensional marginalized posteriors of ST and TT amplitudes obtained from the ST+TT model. We show both the 1σ and 2σ contours in the two-dimensional plot. |
To evaluate the potential deviation of the ST correlations in IPTA DR2, we also consider a model in which the overlap function is parameterized by the α parameter as equation (4 ). Specifically, α = 0 corresponds to the ST correlations, while α = 3 corresponds to the TT correlations. The marginalized posteriors for the α parameter are shown in figure 4, indicating the IPTA DR2 can be very well described by the ST correlations, but TT correlations are excluded by the 90% credible regions.
Figure 4. Bayesian posteriors for the α parameter in the model with a parameterized overlap function of equation ( |
Based on the above discussions, we conclude that there is strong Bayesian evidence for the ST correlations that could possibly come from some modified gravities with extra scalar polarization modes (see e.g. [55]), but no significant evidence for the TT7(7According to the private communication with Boris Goncharov, we acknowledge that the official IPTA DR2 publication on the search for the TT-mode SGWB is now being submitted. That publication does not report evidence for spatial correlations.) correlations in the IPTA DR2. Therefore, we place 95% upper limits on the amplitude of TT polarization mode as ${{ \mathcal A }}_{\mathrm{TT}}\lesssim 3.95\times {10}^{-15}$ and ${{ \mathcal A }}_{\mathrm{TT}}\lesssim 3.27\times {10}^{-15}$, or equivalently, the 95% upper limits for the energy density parameter per logarithm frequency are ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{TT}}\lesssim 2.16\times {10}^{-8}$ and ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{TT}}\lesssim 1.48\times {10}^{-8}$, from the TT and ST+TT models, respectively.
Note added.—Recently the IPTA collaboration also performed a search for the SGWB in their DR2 data set [56], and confirmed our findings that IPTA DR2 has strong evidence for the UCP, but no significant evidence for the TT mode.