Zu-Cheng Chen,Yu-Mei Wu and Qing-Guo Huang,5

1 CAS Key Laboratory of Theoretical Physics,Institute of Theoretical Physics,Chinese Academy of Sciences,Beijing 100190,China

2 School of Physical Sciences,University of Chinese Academy of Sciences,No.19A Yuquan Road,Beijing 100049,China

3 Department of Astronomy,Beijing Normal University,Beijing 100875,China

4 Advanced Institute of Natural Sciences,Beijing Normal University,Zhuhai 519087,China

5 School of Fundamental Physics and Mathematical Sciences,Hangzhou Institute for Advanced Study,UCAS,Hangzhou 310024,China

Abstract We search for isotropic stochastic gravitational-wave background(SGWB)in the International Pulsar Timing Array second data release.By modeling the SGWB as a power-law,we find very strong Bayesian evidence for a common-spectrum process,and further this process has scalar transverse(ST)correlations allowed in general metric theory of gravity as the Bayes factor in favor of the ST-correlated process versus the spatially uncorrelated common-spectrum process is 30±2.The median and the 90% equal-tail amplitudes of ST mode are or equivalently the energy density parameter per logarithm frequency is at frequency of 1 year?1.However,we do not find any statistically significant evidence for the tensor transverse(TT)mode and then place the 95% upper limits as or equivalently at frequency of 1 year?1.

Keywords:stochastic gravitational-wave background,pulsar timing array,beyond general relativity

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 gravitationalwave 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 commonspectrum 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 reestimated 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]

Table 1.Parameters and their prior distributions used in the analyses.

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 TEMPO2[45,46]through libstempo6https://vallis.github.io/libstempointerface,and∈is a vector denoting small offsets for the parameters of timing model.The second termδtSNis 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 δtDMdenotes 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 δtyrDMis the stochastic contribution caused by annual DM variation described by a deterministic yearly sinusoid(see e.g.[36]).The fifth term δtWNrepresents the stochastic contribution due to the white noise(WN),including a scale parameter on the TOA uncertainties(EFAC),an added variance(EQUAD),and a perepoch 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 δtCPis 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?9Hz to 1.09×10?8Hz 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].

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 δtCPterm 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 chainMonte 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 2.An interpretation of the Bayes factor in determining which model is favored,as given by[54].

Table 3.The lnBFbetween pairs of models.The digit in the parentheses represents the uncertainty on the last quoted digit.

Figure 1.One and two-dimensional marginalized posteriors of the amplitude AUCP 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.

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

Figure 2.Marginalized posteriors of ATT and AST parameters obtained from the TT and ST models,respectively.

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.

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 TT7According to the private communication with Boris Goncharov,we acknowledge that the official IPTA DR2 publication on the search for the TTmode 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 asor equivalently,the 95% upper limits for the energy density parameter per logarithm frequency are

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.

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.

Figure 4.Bayesian posteriors for the α parameter in the model with a parameterized overlap function of equation(4).We consider the power-law SGWB spectrum with both γ=13/3 and γ=5.The two vertical dashed lines correspond to the TT(α=3)and ST(α=0)overlap functions,respectively.

Acknowledgments

We acknowledge the use of HPC Cluster of ITP-CAS and HPC Cluster of Tianhe II in National Supercomputing Center in Guangzhou.This work is supported by the National Key Research and Development Program of China Grant No.2020YFC2201502,grants from NSFC(Grant No.11 975 019,11 991 052,12 047 503),Key Research Program of Frontier Sciences,CAS,Grant NO.ZDBS-LY-7009,CAS Project for Young Scientists in Basic Research YSBR-006,the Key Research Program of the Chinese Academy of Sciences(Grant NO.XDPB15).