The Australian Chinese European Indian and North American pulsar timing array PTA collaborations recently reported at varying levels evidence for the presence of a nanohertz gravitationalwave background GWB. Given that each PTA made different choices in modeling their data we perform a comparison of the GWB and individual pulsar noise parameters across the results reported from the PTAs that constitute the International Pulsar Timing Array IPTA. We show that despite making different modeling choices there is no significant difference in the GWB parameters that are measured by the different PTAs agreeing within 1. The pulsar noise parameters are also consistent between different PTAs for the majority of the pulsars included in these analyses. We bridge the differences in modeling choices by adopting a standardized noise model for all pulsars and PTAs finding that under this model there is a reduction in the tension in the pulsar noise parameters. As part of this reanalysis we extended each PTAs data set by adding extra pulsars that were not timed by that PTA. Under these extensions we find better constraints on the GWB amplitude and a higher signaltonoise ratio for the HellingsDowns correlations. These extensions serve as a prelude to the benefits offered by a full combination of data across all pulsars in the IPTA i.e. the IPTAs Data Release 3 which will involve not just adding in additional pulsars but also including data from all three PTAs where any given pulsar is timed by more than a single PTA.
['Gravitational waves', 'Pulsars', '678', '1306', 'Astrophysics - High Energy Astrophysical Phenomena', 'General Relativity and Quantum Cosmology']
Comparing Recent Pulsar Timing Array Results on the Nanohertz Stochastic Gravitationalwave Background
{'Comparing recent PTA results on the nanohertz stochastic gravitational wave background': 'The International Pulsar Timing Array Collaboration', 'ABSTRACT': "The Australian, Chinese, European, Indian, and North American pulsar timing array (PTA) collaborations recently reported, at varying levels, evidence for the presence of a nanohertz gravitational wave background (GWB). Given that each PTA made different choices in modeling their data, we perform a comparison of the GWB and individual pulsar noise parameters across the results reported from the PTAs that constitute the International Pulsar Timing Array (IPTA). We show that despite making different modeling choices, there is no significant difference in the GWB parameters that are measured by the different PTAs, agreeing within 1 σ . The pulsar noise parameters are also consistent between different PTAs for the majority of the pulsars included in these analyses. We bridge the differences in modeling choices by adopting a standardized noise model for all pulsars and PTAs, finding that under this model there is a reduction in the tension in the pulsar noise parameters. As part of this reanalysis, we 'extended' each PTA's data set by adding extra pulsars that were not timed by that PTA. Under these extensions, we find better constraints on the GWB amplitude and a higher signal-to-noise ratio for the Hellings and Downs correlations. These extensions serve as a prelude to the benefits offered by a full combination of data across all pulsars in the IPTA, i.e., the IPTA's Data Release 3, which will involve not just adding in additional pulsars, but also including data from all three PTAs where any given pulsar is timed by more than as single PTA.", '1. INTRODUCTION': "Pulsar timing arrays (PTAs) seek to detect low frequency gravitational waves (GW) by monitoring a collection of millisecond radio pulsars (Foster & Backer 1990). When a GW is incident on a PTA it induces shifts in the times of arrival of radio pulses. These shifts are correlated between pairs of pulsars depending on their angular separation, known as the Helling-Downs (HD) correlations (Hellings & Downs 1983).  \nThe most likely source of low frequency GWs are supermassive black hole binaries (SMBHB), although cosmological and other more exotic sources are also possible (Burke-Spolaor et al. 2019, and references therein). It is expected that an ensemble of SMBHBs can generate a stochastic gravitational wave background (GWB) that could be detected first, followed by the detection of individually resolvable SMBHB sources (Rosado et al. 2015). The spectrum of a stochastic GWB of SMBHB origin is  \nCorresponding author: Paul T. Baker, Aur'elien Chalumeau, Nihan S. Pol paul.baker@nanograv.org, aurelien.chalumeau@unimib.it,  \nnihan.pol@nanograv.org  \naffected by the evolution of SMBHBs and their environments (Taylor et al. 2017). In the case of purely GWdriven evolution the pulsar timing residuals induced by the GWB follow a power law with a spectral index of γ = 13 / 3 in the convention used in this work, although more realistic backgrounds can deviate from this expectation (Phinney 2001; Sesana et al. 2008; B'ecsy et al. 2022).  \nCurrently several independent groups organized by region are operating PTAs. The International Pulsar Timing Array (IPTA, Verbiest et al. 2016) is a consortium of four PTAs: The European PTA (EPTA, Desvignes et al. 2016), the Indian PTA (InPTA, Joshi et al. 2018), the North American Nanohertz Observatory for Gravitational waves (NANOGrav, Ransom et al. 2019), and the Australia based Parkes PTA(PPTA, Manchester et al. 2013). Additionally the Chinese PTA (CPTA, Lee 2016) and MeerKAT PTA (Miles et al. 2023) are not yet IPTA members, but have observer status within the IPTA.  \nIn 2020 and 2021, the EPTA, NANOGrav, and PPTA reported the discovery of a common uncorrelated red noise (CURN) process in their data, but showed no evidence for or against HD correlations (Chen et al. 2021;  \nArzoumanian et al. 2020; Goncharov et al. 2021). These results were supported by the analysis of the second data release of the IPTA, a combination of older data from these three PTAs, where a consistent CURN process was also found (Perera et al. 2019; Antoniadis et al. 2022).  \nIn Antoniadis et al. (2023a, hereafter EPTA+InPTA), Agazie et al. (2023a, hereafter NANOGrav), and Reardon et al. (2023a, hereafter PPTA), IPTA members presented analyses of their most recent data sets and report evidence for an HD correlated GWB with varying levels of significance. In this work we compare the results from these three studies: describing key features and results of the published data sets in section 2, outlining the signal and noise models in section 3, comparing the properties of the GWB using published posterior samples and a new factorized analysis in section 4, and comparing individual pulsar noise properties in section 5. Finally, in section 6 we extend each PTA by adding in pulsars that are not timed by that PTA in a 'pseudo-IPTA' combination using the factorized likelihood method. Note that this is different from a full data combination, currently underway as part of IPTA's third data release (IPTA-DR3), where the data from every PTA will be combined for each single pulsar using a unified timing and noise model.  \nThe CPTA also reported on an HD correlated GWB in Xu et al. (2023), that appears broadly consistent with the other PTAs. Unfortunately, we did not have access to those data and were unable to include their results in our study.", '2.1. EPTA+InPTA': 'EPTA+InPTA reported results from several permutations of the second data release of the EPTA (Antoniadis et al. 2023b). In this work we consider the analysis using the data EPTA+InPTA named EPTA DR2new+. This data set contains the most recent 10.3 yr of EPTA observations for 25 pulsars. For 10 of these pulsars, about 3.5 yr of additional observations from InPTA were combined. The resulting combination has a total time span of about 11 yr.  \nThe older EPTA data were excluded, because they were collected using legacy observing systems. These observations had narrow radio bandwidths and were mostly collected at L-band (1400 MHz) only. The lack of radio frequency coverage led to a covariance between dispersion measure variations and achromatic red noise, potentially polluting the GW information.  \nEPTA+InPTA reported a Bayes factor of about 60 in favor of HD correlated GWB over CURN. This corresponds to a false alarm probability of about 10 -3  \n(3 σ ). Assuming a power law with fixed spectral index of 13 / 3, EPTA+InPTA recovered an amplitude (2 . 5 ± 0 . 7) × 10 -15 at a reference frequency of 1 yr -1 (median and 90% C.I.).', '2.2. NANOGrav': 'NANOGrav analyzed the NANOGrav 15 yr data set, which contains observations of 68 pulsars (Agazie et al. 2023b). One pulsar was excluded from the GWB analysis, having less than 3 yr of observations.  \nNANOGrav reported a Bayes factor of about 200 in favor of HD correlated GWB over CURN. This corresponds to a false alarm probability of about 10 -3 or 5 × 10 -5 (3 -4 σ ), depending on the background estimation method. Assuming a power law with fixed spectral index of 13 / 3, NANOGrav recovered an amplitude 2 . 4 +0 . 7 -0 . 6 × 10 -15 at a reference frequency of 1 yr -1 (median and 90% C.I.).', '2.3. PPTA': 'PPTA analyzed the third data release of the PPTA, using observations of 32 pulsars over 18 years (Zic et al. 2023). Two pulsars were excluded from the GWB analysis because of the lack of data for one, and the presence of strong steep red noise for the other. The full data release contains additional legacy data, which was also excluded from the GWB analysis.  \nPPTA reported a Bayes factor of about 1.5 in favor of HDcorrelated GWB over CURN and a false alarm probability of 0 . 02 (2 σ ) based on a likelihood ratio statistic. Assuming a power law with fixed spectral index of 13 / 3, PPTA recovered an amplitude 2 . 0 +0 . 3 -0 . 2 × 10 -15 at a reference frequency of 1 yr -1 (median and 68% C.I.).', '2.4. Analysis methods': 'The comparisons of published results reported in this work used samples collected by each PTA from relevant posterior distributions. EPTA+InPTA, NANOGrav, and PPTA all used the enterprise software package to evaluate the likelihood and priors (Ellis et al. 2020; Johnson et al. 2023), and all three PTAs used PTMCMCSampler to collect Markov chain Monte Carlo (MCMC) samples (Ellis & van Haasteren 2017). All new MCMC runs for this work were conducted using the same enterprise and PTMCMCSampler software stack.', '3. SIGNAL AND NOISE MODELS': 'In three series of papers the PTAs presented a variety of analyses using different signal and noise models. Here we discuss the basic models used and point out where different PTAs do things differently. For a comprehensive look at the noise modeling done by each PTA  \nsee Antoniadis et al. (2023c), Agazie et al. (2023c), and Reardon et al. (2023b).', '3.1. Published analyses': 'Each PTA modeled the GWB as a power law Fourier basis Gaussian process (GP) with HD correlations (Lentati et al. 2013; van Haasteren et al. 2009) . The hyper parameters log 10 A HD and γ HD set the characteristic amplitude and spectral index of the power law, which acts as a prior on the Fourier coefficients, which describe the spectrum of the GWB. Each PTA set the fundamental frequency of their Fourier basis as the inverse total observation time, 1 /T obs , however each PTA had a different T obs resulting in different sampled frequencies. Each PTA performed an analysis to determine the number of frequencies across which there was a significant signal in their data, and based on this used 9, 14 and 28 frequencies in the GWB model for EPTA+InPTA, NANOGrav, and PPTA, respectively. This corresponds to frequency ranges of about 3 -28, 2 -30, and 1 . 7 -47 nHz, for the same.  \nEach PTA used a similar power law Fourier basis GP to model intrinsic pulsar red noise (RN). Like the GWB, this noise is achromatic: it does not depend on the radio frequency of the observation. Each pulsar had two parameters log 10 A RN and γ RN , which characterize the power law. The fundamental Fourier frequency for RN matched that for GWB for NANOGrav, while it varied from pulsar to pulsar for EPTA+InPTA and PPTA, corresponding to the observation time of each pulsar. NANOGrav used 30 frequencies to model RN and included RN for every pulsar in their array. EPTA+InPTA performed Bayesian model selection to determine the noise models and related optimal number of frequencies to use for each pulsar. Pulsars with steeper RN spectra required fewer frequencies, the minimum being 10, whereas pulsars with shallower RN spectra required more, the maximum chosen being nearly 150. For some pulsars there was no intrinsic RN detected, so this model was not used for those pulsars. PPTA set the maximum RN frequency to 1 / (240 day). For some pulsars PPTA added an additional high-frequency achromatic noise, modeled as a power law with, maximum frequency 1 / (30 day). This combined model had a steep low frequency spectrum, and a shallow (but not flat) high frequency spectrum.  \nOne of the main sources of noise in pulsar timing data are dispersion measure (DM) variations that induce a chromatic (i.e., dependent on the radio frequency ν ) delay in the pulse times of arrival proportional to ν -2 . EPTA+InPTA and PPTA used a power law Fourier basis GP to model stochastic DM variations, known  \nas DMGP. Again each pulsar gained two parameters log 10 A DM and γ DM , which characterized the power law spectrum of deviations away from an initial deterministic fit. The deterministic component included a second order Taylor series fit to the DM and a solar wind model, which took into account the relative position of the sun and pulsar when an observation was made.  \nEPTA+InPTA performed model selection to determine which pulsars experienced measurable DM variations and to set the number of Fourier frequencies used. PPTA set a maximum frequency of 1 / (60 day).  \nFor seven and one of their pulsars, respectively, PPTA and EPTA+InPTA included an additional nondispersive chromatic variation model, scattering noise, which modeled time delays with a ν -4 radio-frequency dependence. This model used the same power law GP framework. For each pulsar with scattering noise PPTA also added a band noise model, an additional achromatic RN GP that affects only the lowest radio frequency band, ν < 960 MHz. This band noise accounted for excess noise that arose due to assuming a fixed chromatic scaling, ν -4 , in the scatter component.  \nThe main NANOGrav results we considered in this work did not use a DMGP, instead using the DMX model to account for DM variations. DMX is a piecewise fit that determines a constant DM for each group of observations that occur near each other in time, effectively fitting a new DM for each pulsar every ∼ 6 days. This fit was part of the pulsar timing model and was analytically marginalized as part of the GWB search.  \nPPTA also included additional system noise to model time-varying instrumental noise for five pulsars. This model used the achromatic power law GP framework, but applied only to observations made with a particular observing system.', '3.2. Factorized likelihood analysis': 'As described above, the three data sets were processed using distinct timing pipelines, different modeling choices, and even different Fourier bases when the same models were used. To alleviate these procedural differences, we used the factorized likelihood approach (Taylor et al. 2022) to model all the pulsars from the different data sets using a standardized noise model. Despite the possibility that this choice might be suboptimal for some pulsars, this allowed us to make a more direct comparison between the properties of the noise and the GWB as seen by the difference PTAs.  \nIn the factorized likelihood analyses, we used the NANOGrav 15 yr data set (Agazie et al. 2023b), the PPTA DR3 (Reardon et al. 2023a), and the DR2New+ data set, which is a combination of the newest EPTA  \ndata with the InPTA data (Antoniadis et al. 2023a). For these data sets, we selected a global time span based on the longest observed pulsar in the data set, T span = 18 . 9 yrs. This T span was used to define the fundamental frequency for all models with a Fourier domain GP. The frequencies in all Fourier bases were integer multiples f = n/T span , where n = 1 , 2 , 3 , . . . , n max .  \nThe standardized noise model contained four pieces: intrinsic pulsar RN, interstellar DM variations, a deterministic solar wind model, and a fixed spectral index CURN to account for the GWB. The intrinsic RN was modeled using n max = 30 Fourier components for all pulsars. Each pulsar had two hyperparameters characterizing the spectrum: amplitude, log 10 A RN , and spectral index, γ RN . We modeled the DM variations in each pulsar using a n max = 100 component power law GP (e.g., Antoniadis et al. 2023c). This necessitated the removal of the DMX parameters in the pulsar timing models provided by NANOGrav, which were replaced by Taylor series expansions of DM up to a second order polynomial term, matching the DM modeling done by EPTA+InPTA and PPTA. Each pulsar had two additional hyperparameters characterizing the DM variation spectrum: amplitude, log 10 A DM , and spectral index, γ DM . The deterministic solar wind was modeled as described in Hazboun et al. (2022), with the solar wind free electron density at the earth, n Earth , as the only free parameter. The pulsar J1713+0747 has shown evidence for multiple chromatic events, causing sudden radio frequency dependent delays that relax over time, and two of these events are contained within the data sets analyzed here. We modeled these events as a deterministic exponential delay, where the amplitude, epoch, recovery time scale, and radio frequency dependence (i.e., ν -x ) are parameters of the model (Antoniadis et al. 2022, 2023c). Finally, the CURN was modeled as a power law process with n max = 15 Fourier components and a fixed spectral index of γ CURN = 13 / 3. The amplitude, log 10 A CURN , was left as a free parameter. After the individual pulsar analyses were completed the posterior distributions for the individual log 10 A CURN s were combined to determine a truly joint posterior for the CURN amplitude.', '4.1. Comparing the published GWB measurements': "In previous IPTA work we compared the posterior distributions of GWB parameters as observed by different PTAs using the Mahalanobis distance, which assumes the posteriors to be multivariate Gaussian distributions (Antoniadis et al. 2022). In this work we adopted a nonGaussian tension metric proposed by Raveri & Doux  \n(2021) to assess the tension between the posterior distributions, obtained by different PTAs. To compute the metric, one first constructs a difference distribution by drawing a sample from each posterior distribution and taking the difference between the two, resulting in a single point of the difference distribution. This process is repeated until the difference distribution is sufficiently mapped. Thereafter, one determines the probability that the difference is zero by integration. The tension metric is then the Gaussian-equivalent ' σ ' value corresponding to the zero-difference probability. We performed a detailed comparison of all GWB and noise parameters using a modified version of the tensiometer 1 package, some of which is presented here. The full comparison including visualizations for all GWB and noise parameters is available as supplementary material 2 .  \nIn EPTA+InPTA, NANOGrav, and PPTA each PTA searched for an HD correlated GWB with a power law spectrum described by a spectral slope, γ , and characteristic amplitude, log 10 A . The 2D marginal posterior for these parameters from each PTA is shown in the right panel of Figure 1. To determine the joint posterior, we computed the simple intersection of the three individual PTA posterior distributions, by multiplying them and renormalizing. We did not weight the individual distributions in any way during this process. This is not, strictly speaking, fair as the posteriors are not independent, some of the pulsars were observed by multiple PTAs. However, this simplistic combination shows that combining the output of the three experiments can provide an improvement in parameter estimation.  \nThe left panel of Figure 1 shows the free spectral posteriors from each PTA. The individual violins show the measured HD correlated GWB amplitude at the corresponding frequency. Additionally, the joint power law GWB posterior spectrum is plotted using the median parameters, showing it to be in good agreement with the individual free spectral results.  \nFigure 2 shows the difference distributions for the power law GWB parameters for each pair of PTAs as computed by tensiometer . There is excellent agreement between all three pairs as the tension metric is less than 1 σ in all. Note that EPTA+InPTA had the widest posterior distribution resulting in wider difference distributions for the two combinations involving those data. While the numerical difference between NANOGrav and PPTA is small, the tension is larger owing to the narrower posterior distributions.  \n<!-- image -->  \nFigure 1. Left : Free spectral posteriors for each PTA showing the measured HD correlated GWB power in several frequency bins under no spectral shape assumption. Each PTA used a different Fourier basis set by their own maximum observing time. The dashed line shows a power law spectrum as determined by the joint 2D power law posterior median. Right : 2D posterior for HD correlated power law GWB parameters. Contours show 68, 95, and 99.7% of the posterior mass. The vertical dotted line is at γ = 13 / 3.  \n<!-- image -->  \nFigure 2. Difference distributions for GWB parameters between pairs of PTAs as computed by tensiometer . The contours show 68 and 95% of the distribution mass.  \n<!-- image -->", '4.2. Comparing the GWB sensitivity of PTAs': "A commonly used measure of GW detector performance is a frequency-dependent 'sensitivity curve'. This metric, which estimates the smallest amplitude of a GW induced signal that a detector would detect, is often used in the GW community to assess detector performance (see Moore et al. 2014; Hazboun et al. 2019b; Kaiser & McWilliams 2021, and references therein). The hasasia (Hazboun et al. 2019a) package offers a means to efficiently compute such curves for PTAs. Specifically, the sensitivity curves we compare here are the sensitivity to interpulsar cross-correlations induced in the PTA by a GWB. As input, hasasia uses the original time of arrival data and the median noise parameters for all noise  \nprocesses, including the GWB auto-correlations which act as noise when trying to detect the cross-correlations.  \nIn order to generate sensitivity curves for EPTA+InPTA and PPTA, we made a few modifications to hasasia . This is because hasasia accounts for white noise and achromatic RN only. For analyses like NANOGrav, which modeled DM variations using DMX (which appears in the timing model) this is sufficient, but it is not for analyses that use DMGP, like  \nFigure 3. Estimated sensitivity to the characteristic strain induced by a GWB as a function of GW frequency. The dashed line shows a power law spectrum as determined by the joint 2D power law posterior median from the right panel of Figure 1.  \n<!-- image -->  \nEPTA+InPTA and PPTA. The system and band noise models used by PPTA must also be accounted for 3 .  \nThe three resulting sensitivity curves are presented in Figure 3 and show the relative sensitivity of the PTA data sets. The general behavior of the curves follows the simple expectations based on the intrinsic properties of each data set. The differing low frequency sensitivity supports the reported evidence from each PTA for the presence of an HD correlated GWB signal. The NANOGrav data set shows the best low frequency sensitivity, reaching slightly lower frequencies than EPTA+InPTA due to the longer observing baseline of analyzed data. The PPTA data set which spans the longest time extends to the lowest frequencies but does not achieve the same low frequency sensitivity as the other two. At the higher frequencies, the EPTA+InPTA and PPTA data sets are more sensitive than NANOGrav due to their higher observing cadences with observations occurring every ∼ 3 and ∼ 7 days, respectively. NANOGrav suffers at the high frequency end due to its lower observing cadence of roughly 30 or 14 days, depending on the pulsar.  \nFigure 4. Amplitudes of CURN recovered using the factorized likelihood method. Extended data sets are also shown, where the pulsars of one data set are added to the another, without repeating pulsars. The boxes contain 68% of the distribution mass, and the center line marks the median. The whiskers contain 95% of the distribution mass.  \n<!-- image -->  \n4.3. Comparison using standardized noise models  \n4.3.1. Common uncorrelated red noise amplitude  \nIn order to make a fair comparison of the observed GWB properties we reanalyzed each PTA's data using the standardized noise models described in subsection 3.2. In place of a full Bayesian analysis searching for HD correlated GWB, which is computationally expensive, we used the factorized likelihood approach. Using this method individual pulsars were analyzed independently in parallel and the results combined to arrive at a posterior distribution for the amplitude of CURN, log 10 A CURN , assuming a fixed spectral index of γ CURN = 13 / 3. This method did not include interpulsar cross-correlations, but acted as a good proxy to quickly determine the spectral properties of a common signal like the GWB.  \nWe applied the standardized noise model described in subsection 3.2 to every pulsar with a timing baseline longer than 3 years. Only 4 pulsars were dropped due to this time cutoff: J0614 -3329 from NANOGrav and J0900 -3144, J1741+1351, and J1902 -5105 were dropped from PPTA. Each pulsar was independently analyzed and the posteriors for the pulsars from a given PTA were combined, resulting in the CURN posteriors shown in color in Figure 4. There is broad agreement between the PTAs, and these new results agree well with the fixed spectral index GWB amplitude reported by the PTAs and stated in section 2. The black boxes in Figure 4 are based on extending the individual PTA data sets and will be discussed in subsection 6.2.", '4.3.2. Dropout factors': "The factorized likelihood method can also be used to determine a pulsar's level of support, or lack thereof, for the CURN seen by the rest of the array (Aggarwal et al. 2019; Taylor et al. 2022; Arzoumanian et al. 2020). This is done by using leave-one-out cross validation, comparing the ratio of Bayesian evidence for a CURN seen by the entire array versus that of the array without the pulsar in question. This calculation gives a dropout factor (DF) for each pulsar, where a DF greater than one means that the pulsar is consistent with the common signal seen by the rest of the PTA, and a DF less than one means that the pulsar is in tension with the rest of the PTA. A DF near one means that the pulsar does not have strong support for or against the common signal.  \nThe left panel of Figure 5 shows the distribution of DFs for each PTA. Many pulsars, particularly those with short observing baselines, had a DF near unity. The right panel of Figure 5 shows the DFs for all pulsars observed by more than one PTA, with uncertainties determined by bootstrap resampling (which are very small). There are notable discrepancies for several pulsars examined. In some cases the differences were driven by observation time. If a PTA had a short observing baseline for a pulsar, it should not be sensitive to low frequency effects, and it should have no support for or against CURN. This was seen in J0030+0451 where PPTA have only 4 yr of data, and the other two PTAs have more than 10 yr. In J2124 -3358 the case was reversed with NANOGrav having only 3 . 5 yr of data. Three additional pulsars J1744 -1134, J1600 -3053, and J1857+0943 were observed with long timing baselines by each PTA, but this analysis resulted in disagreement in the DF among the PTAs. Understanding these differences will be an important use for the forthcoming IPTA-DR3. In some cases there was broad agreement, with pulsars like J1909-3744 and J1024-0719 having DFs greater than 1 and less than 1 for all PTAs, respectively. Interestingly, J1713+0747 had a high dropout factor when using the PPTA data in this analysis, while in PPTA it had the lowest support for CURN. The difference is likely attributed to noise modeling. Our standardized noise model differs from the one employed in PPTA, where several additional components were modeled.", '4.4. Comparing interpulsar correlations': 'Included in both the EPTA and NANOGrav analyses were searches for cross-correlation signals using a frequentist approach called the optimal statistic (Anholm et al. 2009; Chamberlin et al. 2015). This approach estimates the power of a cross-correlated GWB signal, ˆ A 2 gw ,  \nby using a variance weighted least squares fit of individual pulsar pair correlated power. This results in a signal-to-noise ratio (SNR) for a cross-correlated GWB. The noise marginalized optimal statistic (NMOS) takes into account the posterior spread of the noise and CURN resulting in a distribution of ˆ A 2 gw and SNRs (Vigeland et al. 2018).  \nNote that Romano et al. (2021) showed that the ˆ A 2 gw recovery is biased in cases where the amplitude has a comparable effect on the residuals as the noise. This can be attributed to the fact that the optimal statistic does not account for pulsar pair covariances from a common GWB. However, the use of SNR as a detection statistic is still valid. As a result, we will focus only on the SNR calculated using the optimal statistic framework.  \nThe process of using a factorized likelihood analysis with the NMOS has been detailed within Taylor et al. (2022). In this case we chose to marginalize over the CURN and the intrinsic pulsar RN only, while holding the DMGP parameters fixed at their maximum likelihood values to reduce the computational complexity of the process. The resulting SNR distributions for each PTA are shown with colors in Figure 6. We measuree a median SNR of 2.2, 4.9, and 1.0 for EPTA+InPTA, NANOGrav, and PPTA, respectively. The black boxes show extended data sets that will be discussed in subsection 6.2.  \nWhen comparing these NMOS results with those published by NANOGrav and EPTA+InPTA, there are some discrepancies. For EPTA+InPTA there was a decrease in SNR from around 4 to just over 2, while for NANOGrav, there was an increase in SNR from roughly 4 to 5. Both of these differences could be attributed to the differences in pulsar noise modeling. Our analysis modeled each pulsar with the standardized noise model, using a Fourier basis set with a time span of just under 19 years. This is in contrast to the EPTA+InPTA published results in which they used data-driven customized noise models for each pulsar. This, combined with the fact that the time span used to define the Fourier basis here was significantly longer than the EPTA+InPTA data set may explain the overall lower SNR. For NANOGrav, our use of DMGP differed from the fiducial DMX model and could account for our observed increase in SNR.', '5. COMPARING PULSAR NOISE PROPERTIES': "While white noise arises during a particular observation and is partially associated with instrumental effects, other noise and signals can be astrophysical in nature and should be seen by all PTAs. As discussed in section 3 each PTA had its own methods to model noise.  \nFigure 5. left : The distribution of dropout factors (DF) for each PTA. DF is a measure of a pulsar's support for ( > 1) or against ( < 1) CURN. The majority of pulsars have a DF of one, neither supporting nor rejecting the CURN seen by the rest of the array. right : DFs for pulsars that are observed by more than one PTA, allowing for data set comparison under the standardized noise model. There are serious discrepancies for many of the pulsars, indicating the importance for fully combining the data set in IPTA-DR3  \n<!-- image -->  \nFigure 6. The noise marginalized optimal statistic distributions of signal-to-noise ratio (SNR) for each single and combination of PTAs represented as box and whiskers. The boxes contain 50% of the distribution mass, and the center line marks the median. The whiskers contain 95% of the distribution mass. Note that the SNR measured for the EPTA+InPTA dataset is lower than that reported in EPTA+InPTA due to our choice of Fourier basis frequencies being sub-optimal for that dataset.  \n<!-- image -->  \nWhatever the method, each PTA accounted for the two main expected noise components: achromatic RN intrinsic to each pulsar and DM variations, respectively referred as 'RN' and 'DM' in what follows. Checking for consistency in these recovered noise model parameter  \nposteriors acts as a further check of agreement between individual PTAs. It should be noted that our standardized noise models differ from the noise modeling done by each PTA. In particular they are simpler than those used by EPTA+InPTA and PPTA and may not be the ideal noise model for many of the pulsars considered.  \nThere are 11 pulsars that were present in all the three data sets, and 10, 1, and 5 pulsars were shared among EPTA+InPTA and NANOGrav, EPTA+InPTA and PPTA, and NANOGrav and PPTA, respectively (cf. Figure 9). Here we compare their noise properties measured with each data set (1) as reported in the published posterior samples from noise analyses and (2) from new noise analyses that used the standardized models described in subsection 3.2, but omitting the CURN component, marginalizing over solar wind electron density error from the timing model, and setting the fundamental Fourier frequency to be 1 /T all , with T all the time duration between the first and the last observing epoch among all data sets.  \nTable 1 shows the tension metric computed for the 27 pulsars which were observed by multiple PTAs. For each pair of PTAs, the tensions for the RN (∆RN) and the DM variations (∆DM) 2D posterior distributions as reported in the published results (left) and obtained with the new analysis (right) are shown. The large majority of cases showed consistent results among the PTAs for both RN and DM variations. For most cases, the consistency was improved with the new analysis using standardized noise models. For some pulsars, significant tension can be explained by differences in observing time  \nspan or differences in observed radio bandwidth. The tensions ≥ 3 σ are highlighted in the table, and the main discrepancies are discussed in more detail below.  \nWhen comparing the published results of EPTA+InPTA and NANOGrav, only the RN of PSR J1012+5307 has a significant tension (here > 3 . 5 σ ), which is greatly reduced when using the standardized noise models. Figure 7 shows the time-domain realisations of the achromatic RN from the new analysis, which shows common features between the two data sets that are consistent with those reported in Antoniadis et al. (2023c). The tension of 2 . 3 σ for RN of PSR J1857+0943 reflects the fact that the achromatic red noise is well constrained for NANOGrav but not for EPTA+InPTA. The tension between DM variations are measured at 1 . 1 σ and ≥ 3 . 5 σ , respectively, for PSRs J2322+2057 and J1730 -2304. The latter is likely due to the difference in observing time span that is less than 4 years for NANOGrav and more than 10 years for EPTA+InPTA. A detailed explanation for the remaining pulsars with tensions ≥ 1 σ : PSRs J0030+0451, J0613 -0200, J1022+1001, J1640+2224, J1713+0747, J1744 -1134 and J2124 -3358, is provided below.  \nComparing the published results of EPTA+InPTA and PPTA, we observed 2 and 6 tensions larger than 1 σ for RN and DM variations, respectively. The use of standardized noise models significantly reduced the tension in RN for both PSRs J0900 -3144 and J1713+0747, and the tension in DM for PSRs J1022+1001, J1713+0747, J1744 -1134, J1909-3744 and J2124 -3358. However, it did not improve DM consistency for PSR J0613 -0200. Interestingly, the tension grew for the RN of PSRs J1022+1001 and J1909 -3744 and for the DM of J1600 -3053 and J1857+0943. Nevertheless, the timedomain realizations for the three latter were highly consistent between EPTA+InPTA and PPTA. The slight discrepancy ( < 2 σ ) in the noise parameter posteriors for these pulsars was likely caused by the significant difference of observing time spans.  \nThe tension between the published RN posteriors in NANOGrav and PPTA were larger than 1 σ for 6 pulsars, and significantly reduced after using standardized noise models for 5 of them. This clearly shows how important the choice of noise modeling is on the observed RN, which for the noise analysis implicitly included any GWB signal. The new noise analysis resulted in more tension for PSRs J0613 -0200 and J1744 -1134, which are discussed below. The shorter observing time span of NANOGrav for PSRs J1730 -2304 ( ∼ 3 yrs vs. > 17 yrs) and J0437 -4715 ( ∼ 5 yrs against > 15 yrs) was very likely the cause of the large tension for these two pulsars. For PSR J1832 -0836, the tension of 2 . 4 σ for  \nDM is likely caused by a lack of cadence and radio frequency coverage for PPTA.  \nLet us now focus on the main remaining inconsistencies:  \n- · PSR J0030+0451 -This pulsar is known for its low ecliptic latitude, yielding significant solar wind effects that contribute to DM variations. The measured RN was consistent among the 3 data sets ( < 1 σ ). However, the posterior distributions for the DM were fully unconstrained for EPTA+InPTA and PPTA and highly constrained to a flat power law for the NANOGrav. The lack of sensitivity was likely caused by the low radio frequency resolution for EPTA+InPTA and due to the short observing time span for PPTA ( ∼ 4 yrs).\n- · PSR J0613 -0200 -The DM variations were highly consistent between NANOGrav and PPTA (0 . 7 σ ), but EPTA+InPTA displays a flatter spectrum power law. This difference was likely caused by a longer time span of low radio frequency ( < 1 GHz) observations for the former two PTAs. The tension metric for the RN showed greater consistency between EPTA+InPTA and PPTA than for NANOGrav. However, the posterior distribution was poorly constrained for EPTA+InPTA and was constrained to differing values in the other two, as seen in Figure 7. Despite this, the time-domain realizations in Figure 7 show a consistent long term trend for NANOGrav and PPTA.\n- · PSR J1022 + 1001 -Observations of this pulsar are more affected by the solar wind due to its low ecliptic latitude (Tiburzi et al. 2021). It also exhibits behavior consistent with profile evolution (Padmanabh et al. 2020, and references therein). Both of these factors made its chromatic and achromatic noise components more difficult to model than other pulsars. When using the standardized noise models, the DM variations between EPTA+InPTA and PPTA were broadly consistent (2 σ ). The short observing time span for NANOGrav ( ∼ 5 yr) compared with the two other ( > 10 yr) was a potential cause for the observed tension, yielding a consistent spectral slope, but a higher amplitude at 1 yr -1 . The RN was poorly constrained for both EPTA+InPTA and NANOGrav. When using standardized noise models, the tension between EPTA+InPTA and PPTA increased from < 0 . 1 to 2 . 5 σ , with the posterior distribution changing from unconstrained to constrained to a flat power law for PPTA.  \n- · PSR J1640+2224 -The tension metric for the RN was only 0 . 2 σ between EPTA+InPTA and NANOGrav, where the former was unconstrained and the latter found a slightly constrained steep power law. However, the DM variations had a tension of 2 . 3 σ , where EPTA+InPTA had a broader constraint for a very flat power law. As expected, we observed a better constrained posterior distribution on DM variations for NANOGrav, which had better radio frequency coverage, with lowfrequency data ( < 500 MHz) for all of the observing time span. This pulsar was not part of PPTA DR3.\n- · PSR J1713+0747 -The DM variation parameters were highly constrained and mainly consistent for the three data sets. The tension of 1 . 9 σ between NANOGrav and PPTA reflects the difference in the amplitude that was slightly larger for the latter. All the tensions for the RN are also less than 2 σ , but slight differences were observed. The RN power law was flatter for EPTA+InPTA compared with the two others. Despite the very low tension ( < 0 . 1 σ ) between NANOGrav and PPTA, the posterior distributions contained visible differences: while the first was constrained to a single peak, the second was bimodal, with one mode consistent with NANOGrav, and the second favoring a steeper power law ( γ > 4 . 5).\n- · PSR J1744 -1134 -The constraints on DM variations were very consistent among the three data sets. Despite the low tension metric for the RN, we observed different behavior for each data set: unconstrained for EPTA+InPTA, steep power law broadly constrained for NANOGrav and flat power law highly constrained for PPTA.\n- · PSR J2124 -3358 -The RN was consistent among the three data sets, being fully unconstrained for EPTA+InPTA and NANOGrav, while poorly constrained to a steep power law with PPTA. This behavior appears to be driven by the data set time spans which are ∼ 4, ∼ 10 and > 17 years for NANOGrav, EPTA+InPTA and PPTA, respectively. However, the DM variations were fully unconstrained for PPTA and poorly constrained for the two others, favoring flat power laws.  \nThe use of standardized noise models allowed us to show high consistency among the three data sets. The remaining discrepancies discussed above will be studied  \nin detail during the preparation for the full data combination of the IPTA-DR3. All plots produced for this analysis are available as supplementary material 4 .", '6.1. Combined free spectral refit': "A realistic GWB spectrum will likely have more complicated features than a simple power law. The γ = 13 / 3 power law is the simplest model for a population of circular SMBHBs that are driven towards inspiral purely due to GW emission (Phinney 2001), although realistic backgrounds will deviate from this primitive spectral model (Sesana et al. 2008; B'ecsy et al. 2022). For example, three-body interactions between the SMBHB and stars in the loss cone and interactions between the SMBHB and its circumbinary gas disk will also drive the inspiral at larger SMBHB separations (Sesana 2013; BurkeSpolaor et al. 2019). This attenuates the GWB spectrum at lower frequencies, creating a 'turnover' spectrum (Sampson et al. 2015). Detecting such a turnover could further resolve the 'last parsec problem' (Milosavljevi'c & Merritt 2003; Khan et al. 2013; Vasiliev et al. 2014; Agazie et al. 2023d; Antoniadis et al. 2023d), though we note that other sources have been proposed as a source of the GWB that can also produce a turnover spectrum (Afzal et al. 2023; Antoniadis et al. 2023d). This turnover spectrum can be modeled using two power laws with different spectral indices and a bend frequency where the transition occurs.  \nTypically, the GWB is assumed to be a Gaussian process. However, if SMBHBs are the source of the signal, it is possible that the occupation fraction of SMBHBs emitting in each frequency bin is too small and will cause a deviation from the Gaussianity assumption. To account for this possibility, we can apply a weight factor to each frequency of the standard power law which is distributed by a Student's t distribution, modeling nonGaussian deviations from the power law spectrum. This model is known as a ' t -process' (Agazie et al. 2023a). A localized peak in the spectrum will result in a large weight measured at that frequency.  \nThe forthcoming IPTA-DR3 combined data set will allow for a precise search for turnovers or deviations away from a pure power law in the GWB spectrum. However, we can make a pseudo-IPTA combination using the fast spectral refit methods of Lamb et al. (2023), where we can estimate the recovered spectral index and amplitude  \nTable 1. Tension metric (in Gaussian equivalent ' σ ' units) between the data sets for noise model parameters, sorted by observing PTA. If one PTA did not observe the pulsar, 'Ø' is shown. If one PTA did not use the relevant noise model ' -' is shown. Values left of the divider are calculated using the published posteriors from individual PTAs, which make different noise modeling choices. Values right of the divider are calculated using the new noise analysis, in which data from all PTAs are analyzed using the same noise models. Instances with tension ≥ 3 σ are in boldface.  \nby simultaneously refitting to the published free spectral posteriors of the three data sets (see Figure 1).  \nWe used the ceffyl (Lamb et al. 2023) software package to represent each posterior with highly optimised kernel density estimators (KDEs). For a given set of spectral parameters, we computed the probability density of the GWB spectral model at each frequency according to the KDE. The effective likelihood is simply the product of these probability densities over frequency bins. We confirmed that the posterior recovered when fitting a power law using this method is consistent with the joint posterior shown in subsection 4.1.  \nThe posteriors for this joint Bayesian fitting are shown in Figure 8. We recovered consistent posteriors across the turnover, t -process and power law models, which suggests that no spectral features beyond a power law  \nare favored by the data. The turnover spectrum had a broadened posterior, showing more support for the γ = 13 / 3 power law. However, this was because the weak constraints on the first frequency bins of each data set allowed some support for a turnover near the low frequency cutoff, and hence a steeper higher-frequency spectral index than for a pure power law. The t -process did not show significant support for any deviations away from the power law, except at high frequencies where the spectrum was white noise dominated. There are some caveats to this result. As explained previously, each PTA used different methods to model their noise sources, such as DM variations. In addition, some pulsars were included in two or more of the data sets and are therefore being 'double counted' in this pseudo-IPTA combination.  \nFigure 7. left: Time-domain realizations of achromatic red noise for PSR J1012+5307 measured with EPTA+InPTA and NANOGrav, using the standardized noise models. right: Time-domain realizations of achromatic red noise for PSR J0613 -0200 measured with EPTA+InPTA, NANOGrav and PPTA, also using the standardized noise models. The plot at the right displays the 2D marginalized posterior of the RN amplitude A RN and spectral slope γ RN (68, 95, 99.7% contours). For both time-domain figures, the colored areas represent 2 σ credible intervals obtained from 300 realizations, and the colored dots show the medians at observation epochs.  \n<!-- image -->  \nFigure 8. Simultaneously fitting GWB spectral models to the free spectra of the three data sets. Left: 68%- and 95%-credible regions for the recovered GWB amplitude and spectral index for a power law, turnover, and t -process model. They are all consistent, however the turnover model is significantly less constrained and has more support for a γ = 13 / 3 process. Right: The free spectral posteriors from each data set overlaid by power laws constructed from random draws from the ceffyl -generated posteriors. Colors correspond to the models on the left panel.  \n<!-- image -->", '6.2. Extending individual PTAs using factorized likelihood': "The flexibility and speed of the factorized likelihood method introduced in subsection 3.2 provides an immediate approach to extending individual PTAs by sequentially adding pulsars observed by one PTA to those from another (being careful to only have one version of each pulsar) to achieve a pseudo-IPTA data set. We did this pseudo-IPTA construction in a piece-wise fashion, where we first chose a 'base' PTA data set, and then added in pulsars from other data sets that are not timed by the base data set.  \nSince we only add the pulsars not already present in the data set, this addition operation is not commutative. For example, adding the new pulsars from EPTA+InPTA to NANOGrav means that the pulsars in common will use NANOGrav data. The number of  \noverlapping pulsars for each PTA are shown in Figure 9 where EPTA+InPTA is represented by the blue circle, NANOGrav by orange, and PPTA by green. Combining all 3 data sets results in 6 different permutations to consider. This is similar to the approach of Mingarelli et al. (in prep) to combine individual pulsars from different PTAs in a 'lite' combination.  \nThe effect of adding pulsars within the factorized likelihood analysis, producing a pseudo-IPTA data set, can be seen in the set of CURN posteriors which are shown in black in Figure 4. In each case, adding pulsars drawn from the other PTAs resulted in a more constrained CURN amplitude posterior than given by the original PTA alone, with a median thinning of the 68% credible interval by 15%. As a direct result of this, the Bayes factor for CURN over intrinsic RN alone also increased, with a median increase of 7 orders of magnitude. This is  \nFigure 9. An unweighted Venn diagram showing the overlapping pulsars between each PTA's data sets. Note that for the factorized likelihood analyses 4 pulsars were dropped from the PPTA DR3 due to baselines for those pulsars being shorter than 3 years, but are included in this diagram for completeness.  \n<!-- image -->  \nconsistent with the findings of the IPTA-DR2 analysis (Antoniadis et al. 2022), where improvements to parameter estimation and detection significance are observed when PTA data sets are extended.  \nStarting from each of the factorized likelihood CURN posterior from the above analysis, we again used the NMOS to determine the SNR distributions for each pseudo-IPTA data set. This method shows how adding in pulsars from other PTAs affects the significance of correlated power.  \nFigure 6 shows each individual PTA's NMOS SNR recovery alongside each possible pseudo-IPTA data set. These SNR distributions show that it did not matter which PTA we started with, adding additional pulsars will always result in a higher SNR. This is consistent with scaling relations for the optimal statistic found in Siemens et al. (2013), and shows the strong promise of a full data combination.  \nComparing the different combinations with each other we see that there is a preference to the ordering of PTA additions. Starting with NANOGrav data first in the process results in the highest SNRs, and adding EPTA+InPTA data before PPTA results in higher SNRs. These preferences are consistent with each PTA's individual SNR.", '7. CONCLUSIONS': "We performed a comparison of the noise models used and the properties of the GWB recovered in the most recent data-sets of EPTA+InPTA, NANOGrav, and PPTA. We found that a majority of the noise parameters were consistent between the three different data sets. Where there were significant ( > 3 σ ) differences, they could be attributed to different time spans and ca-  \nences as well as the lack of frequency coverage in individual PTAs. These tensions also reduced significantly when standardized noise models were adopted. We also calculated and compared the relative sensitivity of the three data sets using sensitivity curves, showing that the different levels of reported evidence for the GWB were consistent with the sensitivities of each of these data sets.  \nDespite the different noise models used by each PTA, the GWB posterior distributions for all three PTA data sets were consistent within 1 σ as calculated by the tension metric. The GWB spectra measured from all three PTA data sets were consistent with a single, 'joint' power law spectrum, with no evidence for deviations from this power law spectral template.  \nFinally, we extended each of the three data sets using the factorized likelihood approach, to achieve a pseudoIPTA analysis. The amplitudes of the CURN process estimated from a global posterior using different permutations of PTA extensions were found to be broadly consistent with each other as well as with those reported by the individual PTAs. The addition of pulsars to any PTA also resulted in an increase in the measured HD SNR.  \nThe members of the IPTA, along with the MeerKAT PTA (Miles et al. 2023), are currently in the process of combining their most recent data sets, which will become IPTA-DR3. The comparisons presented here motivate improvements and best practices to be adopted in a unified analysis for the ongoing data combination by the IPTA. Based on our work we believe that choosing the right noise model for each pulsar through Bayesian model section will be an important step for future IPTA analyses. In particular great care must be taken with the pulsars J0030+0451, J0613 -0200, J1022+1001, J1640+2224, J1713+0747, J1744 -1134 and J2124 -3358 where discrepancies between PTAs persisted when using the standardized noise model. While the results presented here adopted a much simplified approach as opposed to a true combination, these already hint at an enhancement in the significance of GWB detection in the full DR3 over those reported by EPTA+InPTA, NANOGrav, and PPTA individually.", 'ACKNOWLEDGMENTS': "Author contributions. An alphabetical-order author list was used for this paper to recognize the vast extent of the time and effort by many people that went into this project. All authors contributed to the activities of the IPTA collaboration leading to the work presented here, and reviewed the text and figures prior to the paper's submission. This project was organized as part of  \nthe Gravitational Wave Analysis Working Group of the IPTA by PTB, ACh, and NSP.  \nSDan and PTB lead the GW comparisons with contributions from LD, SDe, and AG. ACh lead the noise comparisons with contributions from BBL, SDan, CD, FK, DDeb, and CMFM. GS lead the GW sensitivity analysis with contributions from KGr, NKP, JSH, and KEW. KAG, LSch, and NSP performed the factorized likelihood and optimal statistic analyses. WGL performed the free spectral refitting.  \nPTB, ACh, NSP, SDan, KAG, LSch, WGL, BCJ, and GS wrote the text and generated figures.  \nAcknowledgments. The European Pulsar Timing Array (EPTA) is a collaboration between European and partner institutes, namely ASTRON (NL), INAF/Osservatorio di Cagliari (IT), Max-PlanckInstitut fur Radioastronomie (GER), Nan¸cay/Paris Observatory (FRA), the University of Manchester (UK), the University of Birmingham (UK), the University of East Anglia (UK), the University of Bielefeld (GER), the University of Paris (FRA), the University of MilanBicocca (IT), the Foundation for Research and Technology, Hellas (GR), and Peking University (CHN), with the aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves. The Indian Pulsar Timing Array (InPTA) is an Indo-Japanese collaboration that routinely employs TIFR's upgraded Giant Metrewave Radio Telescope for monitoring a set of IPTA pulsars. The NANOGrav collaboration receives support from National Science Foundation (NSF) Physics Frontiers Center award Nos. 1430284 and 2020265, the Gordon and Betty Moore Foundation, NSF AccelNet award No. 2114721, an NSERC Discovery Grant, and CIFAR. Part of this research was undertaken as part of the Australian Research Council (ARC) Centre of Excellence for Gravitational Wave Discovery (OzGrav) under grant CE170100004.  \nAn Advanced Grant of the European Research Council allowed to implement the Large European Array for Pulsars (LEAP) under Grant Agreement No. 227947 (PI M. Kramer). Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut fur 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 (STFC). The Nan¸cay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We  \nacknowledge 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 a National Facility by the National Institute for Astrophysics (INAF). The Arecibo Observatory is a facility of the NSF operated under cooperative agreement (AST-1744119) by the University of Central Florida (UCF) in alliance with Universidad Ana G. M'endez (UAGM) and Yang Enterprises (YEI), Inc. The Green Bank Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. The National Radio Astronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. Murriyang, the Parkes 64 m radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42), which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Wiradjuri People as the Traditional Owners of the Observatory site.  \nThis work made use of the OzSTAR national facility at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS). This work was conducted in part using the resources of the Advanced Computing Center for Research and Education (ACCRE) at Vanderbilt University, Nashville, TN. This work used resources of the IN2P3 Computing Center (CC-IN2P3 Lyon/Villeurbanne - France) funded by the Centre National de la Recherche Scientifique.  \nThe work is supported by the National SKA programme of China (2020SKA0120100), Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. This work is also supported as part of the 'LEGACY' MPG-CAS collaboration on lowfrequency gravitational wave astronomy. JA acknowledges support from the European Commission (Grant Agreement No. 101094354). JA and SChan were partially supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (HFRI) under the 2nd Call of the 'Science  \nand Society - Action Always strive for excellence Theodoros Papazoglou' (Project No. 01431). ACh acknowledges support from the Paris ˆ Ile-de-France Region. ACh, AF, ASe, ASa, EB, DI, GMS, MBo acknowledge financial support provided under the European Union's H2020 ERC Consolidator Grant 'Binary Massive Black Hole Astrophysics' (B Massive, Grant Agreement: 818691). GD, KLi, RK and MKr acknowledge support from European Research Council (ERC) Synergy Grant 'BlackHoleCam', Grant Agreement No. 610058. ICN is supported by the STFC doctoral training grant ST/T506291/1. AV and PRB are supported by the UK's Science and Technology Facilities Council (STFC; grant ST/W000946/1). AV also acknowledges the support of the Royal Society and Wolfson Foundation. NKP is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer PO 2758/1-1, through the Walter-Benjamin programme. ASa thanks the Alexander von Humboldt foundation in Germany for a Humboldt fellowship for postdoctoral researchers. APo, DP and MBu acknowledge support from the research grant 'iPeska' (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016 (Italy). RNC acknowledges financial support from the Special Account for Research Funds of the Hellenic Open University (ELKE-HOU) under the research programme 'GRAVPUL' (grant agreement 319/10-10-2022). EvdW, CGB and GHJ acknowledge support from the Dutch National Science Agenda, NWA Startimpuls - 400.17.608. BG is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, No. 2017SYRTCN.  \nBCJ, YG, YM, SDan, AG and PR acknowledge the support of the Department of Atomic Energy, Government of India, under Project Identification # RTI 4002. BCJ, YG and YM acknowledge support of the Department of Atomic Energy, Government of India, under project No. 12-R&D-TFR-5.02-0700, while SD, AG and PR acknowledge support of the Department of Atomic Energy, Government of India, under project No. 12R&D-TFR-5.02-0200. KT is partially supported by JSPS KAKENHI Grant Nos. 20H00180, 21H01130, and 21H04467, Bilateral Joint Research Projects of JSPS, and the ISM Cooperative Research Program (2021ISMCRP-2017). AKP is supported by CSIR fellowship Grant No. 09/0079(15784)/2022-EMR-I. SH is supported by JSPS KAKENHI Grant No. 20J20509. KN is supported by the Birla Institute of Technology & Science Institute fellowship. AmS is supported by CSIR fellowship Grant No. 09/1001(12656)/2021-EMR-I and T-641  \n(DST-ICPS). TK is partially supported by the JSPS Overseas Challenge Program for Young Researchers. DDeb acknowledges the support from the Department of Atomic Energy, Government of India through Apex Project - Advance Research and Education in Mathematical Sciences at IMSc.  \nLB acknowledges support from the National Science Foundation under award AST-1909933 and from the Research Corporation for Science Advancement under Cottrell Scholar Award No. 27553. SB gratefully acknowledges the support of a Sloan Fellowship, and the support of NSF under award #1815664. The work of RB, RC, DD, NLa, XS, JPS, and JT is partly supported by the George and Hannah Bolinger Memorial Fund in the College of Science at Oregon State University. MC, PP, and SRT acknowledge support from NSF AST-2007993. MC and NSP were supported by the Vanderbilt Initiative in Data Intensive Astrophysics (VIDA) Fellowship. KCh, ADJ, and MV acknowledge support from the Caltech and Jet Propulsion Laboratory President's and Director's Research and Development Fund. KCh and ADJ acknowledge support from the Sloan Foundation. Support for this work was provided by the NSF through the Grote Reber Fellowship Program administered by Associated Universities, Inc./National Radio Astronomy Observatory. Support for HTC is 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. KCr is supported by a UBC Four Year Fellowship (6456). MED acknowledges support from the Naval Research Laboratory by NASA under contract S-15633Y. TD and MTL are supported by an NSF Astronomy and Astrophysics Grant (AAG) award No. 2009468. ECF is supported by NASA under award No. 80GSFC21M0002. GEF, SCS, and SJV are supported by NSF award PHY2011772. KAG and SRT acknowledge support from an NSF CAREER award #2146016. The Flatiron Institute is supported by the Simons Foundation. SH is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1745301. NLa acknowledges the support from Larry W. Martin and Joyce B. O'Neill Endowed Fellowship in the College of Science at Oregon State University. 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 (80NM0018D0004). DRL and MAMc are supported by NSF #1458952. MAMc is supported by NSF #2009425. CMFM was supported in part by the National Science Foundation  \nunder Grants No. NSF PHY-1748958 and AST-2106552. AMi is supported by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy - EXC 2121 Quantum Universe - 390833306. PN acknowledges support from the BHI, funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation. The Dunlap Institute is funded by an endowment established by the David Dunlap family and the University of Toronto. KDO was supported in part by NSF Grant No. 2207267. TTP acknowledges support from the Extragalactic Astrophysics Research Group at Eotvos Lor'and University, funded by the Eotvos Lor'and Research Network (ELKH), which was used during the development of this research. SMR and IHS are CIFAR Fellows. Portions of this work performed at NRL were supported by ONR 6.1 basic research funding. JDR also acknowledges support from start-up funds from Texas  \nTech University. JS is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388, and acknowledges previous support by the NSF under award 1847938. CU acknowledges support from BGU (Kreitman fellowship), and the Council for Higher Education and Israel Academy of Sciences and Humanities (Excellence fellowship). CAW acknowledges support from CIERA, the Adler Planetarium, and the Brinson Foundation through a CIERA-Adler postdoctoral fellowship. OY is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE2139292.  \nRMS acknowledges support through ARC Future Fellowship FT190100155. SDa is the recipient of an Australian Research Council Discovery Early Career Award (DE210101738) funded by the Australian Government. YL acknowledges support of the Simons Investigator Grant 827103.", 'REFERENCES': "Afzal, A., Agazie, G., Anumarlapudi, A., et al. 2023, The Astrophysical Journal Letters, 951, L11,  \ndoi: 10.3847/2041-8213/acdc91  \nAgazie, G., et al. 2023a, Astrophys. J. Lett., 951, L8,  \ndoi: 10.3847/2041-8213/acdac6  \n- -. 2023b, Astrophys. J. Lett., 951, L9,  \ndoi: 10.3847/2041-8213/acda9a  \n- -. 2023c, Astrophys. J. Lett., 951, L10,  \ndoi: 10.3847/2041-8213/acda88  \n- Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023d, The NANOGrav 15-year Data Set: Constraints on Supermassive Black Hole Binaries from the Gravitational Wave Background. https://arxiv.org/abs/2306.16220  \nAggarwal, K., Arzoumanian, Z., Baker, P. T., et al. 2019, ApJ, 880, 116, doi: 10.3847/1538-4357/ab2236  \nAnholm, M., Ballmer, S., Creighton, J. D., Price, L. R., & Siemens, X. 2009, Physical Review D - Particles, Fields, Gravitation and Cosmology, 79,  \ndoi: 10.1103/PhysRevD.79.084030  \nAntoniadis, J., et al. 2022, Mon. Not. Roy. Astron. Soc., 510, 4873, doi: 10.1093/mnras/stab3418  \n- -. 2023a, arXiv e-prints.  \nhttps://arxiv.org/abs/2306.16214  \n- -. 2023b, arXiv e-prints,  \ndoi: 10.1051/0004-6361/202346841  \n- -. 2023c, arXiv e-prints. https://arxiv.org/abs/2306.16225  \nAntoniadis, J., Arumugam, P., Arumugam, S., et al. 2023d, The second data release from the European Pulsar Timing Array: V. Implications for massive black holes, dark matter and the early Universe.  \nhttps://arxiv.org/abs/2306.16227  \nArzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020,  \nApJL, 905, L34, doi: 10.3847/2041-8213/abd401  \nB'ecsy, B., Cornish, N. J., & Kelley, L. Z. 2022, Astrophys.  \n- J., 941, 119, doi: 10.3847/1538-4357/aca1b2  \nBurke-Spolaor, S., et al. 2019, Astron. Astrophys. Rev., 27,  \n5, doi: 10.1007/s00159-019-0115-7  \n- Chamberlin, S. J., Creighton, J. D., Siemens, X., et al. 2015, Physical Review D, 91, 044048, doi: 10.1103/PhysRevD.91.044048\n- Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970, doi: 10.1093/mnras/stab2833  \nDesvignes, G., et al. 2016, Mon. Not. Roy. Astron. Soc., 458, 3341, doi: 10.1093/mnras/stw483  \nEllis, J., & van Haasteren, R. 2017, jellis18/PTMCMCSampler: Official Release, 1.0.0, Zenodo, doi: 10.5281/zenodo.1037579  \nEllis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2020, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE, Zenodo, doi: 10.5281/zenodo.4059815  \nFoster, R. S., & Backer, D. C. 1990, ApJ, 361, 300,  \ndoi: 10.1086/169195  \nGoncharov, B., Shannon, R. M., Reardon, D. J., et al. 2021, arXiv e-prints, arXiv:2107.12112.  \nhttps://arxiv.org/abs/2107.12112  \n12", 'All Authors and Affiliations': "G. Agazie, 1 J. Antoniadis, 2, 3 A. Anumarlapudi, 1 A. M. Archibald, 4 P. Arumugam, 5 S. Arumugam, 6 Z. Arzoumanian, 7 J. Askew, 8, 9 S. Babak, 10 M. Bagchi, 11, 12 M. Bailes, 8, 9 A.-S. Bak Nielsen, 3, 13 P. T. Baker, 14 C. G. Bassa, 15 A. Bathula, 16 B. B'ecsy, 17 A. Berthereau, 18, 19 N. D. R. Bhat, 20 L. Blecha, 21 M. Bonetti, 22, 23, 24 E. Bortolas, 22, 23, 24 A. Brazier, 25, 26 P. R. Brook, 27 M. Burgay, 28 S. Burke-Spolaor, 29, 30, ∗ R. Burnette, 17 R. N. Caballero, 31 A. Cameron, 8, 9 R. Case, 17 A. Chalumeau, 22 D. J. Champion, 3 S. Chanlaridis, 2 M. Charisi, 32 S. Chatterjee, 25 K. Chatziioannou, 33 B. D. Cheeseboro, 29, 30 S. Chen, 34 Z.-C. Chen, 35, 36 I. Cognard, 18, 19 T. Cohen, 37 W. A. Coles, 38 J. M. Cordes, 25 N. J. Cornish, 39 F. Crawford, 40 H. T. Cromartie, 25, † K. Crowter, 41 M. Curyglyph[suppress]lo, 42 C. J. Cutler, 43, 33 S. Dai, 44 S. Dandapat, 45 D. Deb, 11 M. E. DeCesar, 46 D. DeGan, 17 P. B. Demorest, 47 H. Deng, 17 S. Desai, 48 G. Desvignes, 3 L. Dey, 29, 30 N. Dhanda-Batra, 49 V. Di Marco, 50, 9 T. Dolch, 51, 52 B. Drachler, 53, 54 C. Dwivedi, 55 J. A. Ellis, 56 M. Falxa, 10, 18 Y. Feng, 57 R. D. Ferdman, 58 E. C. Ferrara, 59, 60, 61 W. Fiore, 29, 30 E. Fonseca, 29, 30 A. Franchini, 22, 23 G. E. Freedman, 1 J. R. Gair, 62 N. Garver-Daniels, 29, 30 P. A. Gentile, 29, 30 K. A. Gersbach, 32 J. Glaser, 29, 30 D. C. Good, 63 B. Goncharov, 64, 65 A. Gopakumar, 45 E. Graikou, 3 J.-M. Grießmeier, 18, 19 L. Guillemot, 18, 19 K. Gultekin, 66 Y. J. Guo, 3 Y. Gupta, 67 K. Grunthal, 3 J. S. Hazboun, 17 S. Hisano, 68 G. B. Hobbs, 69 S. Hourihane, 33 H. Hu, 3 F. Iraci, 70, 28 K. Islo, 1 D. Izquierdo-Villalba, 22, 23 J. Jang, 3 J. Jawor, 3 G. H. Janssen, 15, 71 R. J. Jennings, 29, 30 A. Jessner, 3 A. D. Johnson, 1, 33 M. L. Jones, 1 B. C. Joshi, 67, 5 A. R. Kaiser, 29, 30 D. L. Kaplan, 1 A. Kapur, 72, 69 F. Kareem, 73, 74 R. Karuppusamy, 3 E. F. Keane, 75 M. J. Keith, 76 L. Z. Kelley, 77 M. Kerr, 78 J. S. Key, 79 D. Kharbanda, 48 T. Kikunaga, 68 T. C. Klein, 1 N. Kolhe, 80 M. Kramer, 3, 76 M. A. Krishnakumar, 3, 13 A. Kulkarni, 8, 9 N. Laal, 17 K. Lackeos, 3 M. T. Lam, 81, 53, 54 W. G. Lamb, 32 B. B. Larsen, 82 T. J. W. Lazio, 43 K. J. Lee, 83, 84 Y. Levin, 85, 86, 50 N. Lewandowska, 87 T. B. Littenberg, 88 K. Liu, 3 T. Liu, 29, 30 Y. Liu, 84, 13 A. Lommen, 89 D. R. Lorimer, 29, 30 M. E. Lower, 69 J. Luo, 90, ‡ R. Luo, 69, 91 R. S. Lynch, 92 A. G. Lyne, 76 C.-P. Ma, 77, 93 Y. Maan, 67 D. R. Madison, 94 R. A. Main, 3 R. N. Manchester, 69 R. Mandow, 72, 69 M. A. Mattson, 29, 30 A. McEwen, 1 J. W. McKee, 95, 96 M. A. McLaughlin, 29, 30 N. McMann, 32 B. W. Meyers, 41, 20 P. M. Meyers, 33 M. B. Mickaliger, 76 M. Miles, 8, 9 C. M. F. Mingarelli, 82 A. Mitridate, 97 P. Natarajan, 98, 99 R. S. Nathan, 50, 9 C. Ng, 100 D. J. Nice, 101 I. C. Nit¸u, 76 K. Nobleson, 102 S. K. Ocker, 25 K. D. Olum, 103 S. Osglyph[suppress]lowski, 104 A. K. Paladi, 105 A. Parthasarathy, 3 T. T. Pennucci, 106 B. B. P. Perera, 107 D. Perrodin, 28 A. Petiteau, 108, 10 P. Petrov, 32 N. S. Pol, 32 N. K. Porayko, 22, 3 A. Possenti, 28 T. Prabu, 109 H. Quelquejay Leclere, 10 H. A. Radovan, 110 P. Rana, 45 S. M. Ransom, 111 P. S. Ray, 78 D. J. Reardon, 8, 9 A. F. Rogers, 112 J. D. Romano, 113 C. J. Russell, 114 A. Samajdar, 115 S. A. Sanidas, 76 S. C. Sardesai, 1 A. Schmiedekamp, 116 C. Schmiedekamp, 116 K. Schmitz, 117 L. Schult, 32 A. Sesana, 22, 23, 24 G. Shaifullah, 22, 23, 28 R. M. Shannon, 8, 9 B. J. Shapiro-Albert, 29, 30, 118 X. Siemens, 17, 1 J. Simon, 119, § J. Singha, 5 M. S. Siwek, 120 L. Speri, 62 R. Spiewak, 121 A. Srivastava, 48 I. H. Stairs, 41 B. W. Stappers, 76 D. R. Stinebring, 122 K. Stovall, 47 J. P. Sun, 17 M. Surnis, 123 S. C. Susarla, 124 A. Susobhanan, 1 J. K. Swiggum, 101 K. Takahashi, 125, 126 P. Tarafdar, 11 J. Taylor, 17 S. R. Taylor, 32 G. Theureau, 18, 19, 127 E. Thrane, 50, 9 N. Thyagarajan, 128 C. Tiburzi, 28 L. Toomey, 69 J. E. Turner, 29, 30 C. Unal, 129, 130 M. Vallisneri, 43, 33 E. van der Wateren, 15, 71 R. van Haasteren, 131 A. Vecchio, 27 V. Venkatraman Krishnan, 3 J. P. W. Verbiest, 132 S. J. Vigeland, 1 H. M. Wahl, 29, 30 S. Wang, 133 Q. Wang, 32 C. A. Witt, 134, 135 J. Wang, 13, 136, 35 L. Wang, 76 K. E. Wayt, 17 Z. Wu, 84, 13 O. Young, 53, 54 L. Zhang, 137, 69 S. Zhang, 138, 69 X.-J. Zhu 35 And A. Zic 69, 72  \nThe International Pulsar Timing Array Collaboration  \n1 Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee,  \nP.O. Box 413, Milwaukee, WI 53201, USA  \n2 Institute of Astrophysics, FORTH, N. Plastira 100, 70013, Heraklion, Greece  \n3 Max-Planck-Institut fur Radioastronomie, Auf dem Hugel 69, 53121 Bonn  \n4 Newcastle University, NE1 7RU, UK  \n5 Department of Physics, Indian Institute of Technology Roorkee, Roorkee-247667, India  \n6 Department of Electrical Engineering, IIT Hyderabad, Kandi, Telangana 502284, India  \n7 X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA  \n8 Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia  \n9 Australia Research Council Centre for Excellence for Gravitational Wave Discovery (OzGrav)  \n10 Universit'e Paris Cit'e, CNRS, Astroparticule et Cosmologie, 75013 Paris, France  \n11 The Institute of Mathematical Sciences, C.I.T. Campus, Taramani, Chennai 600113, India  \nHomi Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai 400094, India  \n13 Fakultat fur Physik, Universitat Bielefeld, Postfach 100131, 33501 Bielefeld, Germany  \n25  \n55  \n22  \n30  \n15  \n- 14 Department of Physics and Astronomy, Widener University, One University Place, Chester, PA 19013, USA  \nASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands  \n16  \nDepartment of Physical Sciences, Indian Institute of Science Education and Research, Mohali, Punjab 140306, India  \n- 17 Department of Physics, Oregon State University, Corvallis, OR 97331, USA\n- 18 Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Universit'e d'Orl'eans / CNRS, 45071 Orl'eans Cedex 02, France\n- 19 Observatoire Radioastronomique de Nan¸cay, Observatoire de Paris, Universit'e PSL, Universit'e d'Orl'eans, CNRS, 18330 Nan¸cay, France  \n20 International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia  \n21 Physics Department, University of Florida, Gainesville, FL 32611, USA  \nDipartimento di Fisica 'G. Occhialini', Universit'a degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy  \n23 INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy  \n- 24 INAF - Osservatorio Astronomico di Brera, via Brera 20, I-20121 Milano, Italy  \nCornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA  \n26  \nCornell Center for Advanced Computing, Cornell University, Ithaca, NY 14853, USA  \n- 27 Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK  \n28 INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius (CA), Italy  \n29  \nDepartment of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  \nCenter for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV  \n26505, USA  \n- 31 Hellenic Open University, School of Science and Technology, 26335 Patras, Greece\n- 32 Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA  \n33  \nDivision of Physics, Mathematics, and Astronomy, California Institute of Technology, Pasadena, CA 91125, USA  \n34 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing, 100871 China  \n35  \nAdvanced Institute of Natural Sciences, Beijing Normal University, Zhuhai 519087, China  \n- 36 Department of Physics and Synergistic Innovation Center for Quantum Effects and Applications, Hunan Normal University, Changsha, Hunan 410081, China  \n37  \n40  \nDepartment of Physics, New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA  \n38  \nElectrical and Computer Engineering, University of California at San Diego, La Jolla, California, U.S.A.  \n39  \nDepartment of Physics, Montana State University, Bozeman, MT 59717, USA  \nDepartment of Physics and Astronomy, Franklin & Marshall College, P.O. Box 3003, Lancaster, PA 17604, USA  \n- 41 Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada 42 Astronomical Observatory, University of Warsaw, Aleje Ujazdowskie 4, 00-478 Warsaw, Poland  \n43  \nJet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA  \n44  \nSchool of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia  \n- 45 Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, Homi Bhabha Road, Navy Nagar, Colaba, Mumbai 400005, India  \n46 George Mason University, resident at the Naval Research Laboratory, Washington, DC 20375, USA  \n47  \nNational Radio Astronomy Observatory, 1003 Lopezville Rd., Socorro, NM 87801, USA  \n- 48 Department of Physics, IIT Hyderabad, Kandi, Telangana 502284, India\n- 49 Department of Physics and Astrophysics, University of Delhi, Delhi 110007, India  \n50 School of Physics and Astronomy, Monash University, VIC 3800, Australia  \n- 51 Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA  \n52 Eureka Scientific, 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA  \n53 School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA  \n54 Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA  \nDepartment of Earth and Space Sciences, Indian Institute of Space Science and Technology, Valiamala, Thiruvananthapuram, Kerala  \n695547,India  \n56 Bionic Health, 800 Park Offices Drive, Research Triangle Park, NC 27709  \n57 Research Center for Intelligent Computing Platforms, Zhejiang Laboratory, Hangzhou 311100, China  \n58 School of Physics, Faculty of Science, University of East Anglia, Norwich NR4 7TJ, UK  \n59 Department of Astronomy, University of Maryland, College Park, MD 20742  \nCenter for Research and Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD 20771  \n60  \n61 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  \nMax Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Muhlenberg 1, 14476 Potsdam, Germany  \n63  \nDepartment of Physics & Astronomy, University of Montana, 32 Campus Drive, Missoula MT, 59812  \n64 Gran Sasso Science Institute (GSSI), I-67100 L'Aquila, Italy  \n65 INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy  \n62  \n- 66 Department of Astronomy and Astrophysics, University of Michigan, Ann Arbor, MI 48109, USA 67 National Centre for Radio Astrophysics, Pune University Campus, Pune 411007, India 68 Kumamoto University, Graduate School of Science and Technology, Kumamoto, 860-8555, Japan 69 Australia Telescope National Facility, CSIRO, Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia 70 Universit'a di Cagliari, Dipartimento di Fisica, S.P. Monserrato-Sestu Km 0,700 - 09042 Monserrato (CA), Italy 71 Department of Astrophysics/IMAPP, Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands 72 Department of Physics and Astronomy and MQ Research Centre in Astronomy, Astrophysics and Astrophotonics, Macquarie University, NSW 2109, Australia 73 Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur, 741246, India 74 Center of Excellence in Space Sciences India, Indian Institute of Science Education and Research Kolkata, 741246, India 75 School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 PN40, Ireland 76 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK 77 Department of Astronomy, University of California, Berkeley, 501 Campbell Hall #3411, Berkeley, CA 94720, USA 78 Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA 79 University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA 80 Department of Physics, St. Xavier's College (Autonomous), Mumbai 400001, India 81 SETI Institute, 339 N Bernardo Ave Suite 200, Mountain View, CA 94043, USA 82 Department of Physics, Yale University, New Haven, CT 06520, USA 83 Department of Astronomy, School of Physics, Peking University, Beijing 100871, P. R. China 84 National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, P. R. China 85 Physics Department and Columbia Astrophysics Laboratory, Columbia University, 538 West 120th Street, New York, NY 10027, USA 86 Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, NY10011, USA 87 Department of Physics, State University of New York at Oswego, Oswego, NY, 13126, USA 88 NASA Marshall Space Flight Center, Huntsville, AL 35812, USA 89 Department of Physics and Astronomy, Haverford College, Haverford, PA 19041, USA 90 Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada 91 Department of Astronomy, School of Physics and Materials Science, Guangzhou University, Guangzhou 510006, China 92 Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA 93 Department of Physics, University of California, Berkeley, CA 94720, USA 94 Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA 95 E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK 96 Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK 97 Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany 98 Department of Astronomy, Yale University, 52 Hillhouse Ave, New Haven, CT 06511 99 Black Hole Initiative, Harvard University, 20 Garden Street, Cambridge, MA 02138 100 Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George St., Toronto, ON M5S 3H4, Canada 101 Department of Physics, Lafayette College, Easton, PA 18042, USA 102 Department of Physics, BITS Pilani Hyderabad Campus, Hyderabad 500078, Telangana, India 103 Institute of Cosmology, Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA 104 Manly Astrophysics, 15/41-42 East Esplanade, Manly, NSW 2095, Australia 105 Joint Astronomy Programme, Indian Institute of Science, Bengaluru, Karnataka, 560012, India 106 Institute of Physics and Astronomy, Eotvos Lor'and University, P'azm'any P. s. 1/A, 1117 Budapest, Hungary 107 Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA 108 IRFU, CEA, Universit'e Paris-Saclay, F-91191 Gif-sur-Yvette, France 109 Raman Research Institute India, Bengaluru, Karnataka, 560080, India 110 Department of Physics, University of Puerto Rico, Mayaguez, PR 00681, USA 111 National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA 112 Institute for Radio Astronomy & Space Research, Auckland University of Technology, Private Bag 92006, Auckland 1142, New Zealand 113 Department of Physics, Texas Tech University, Box 41051, Lubbock, TX 79409, USA 114 CSIRO Scientific Computing, Australian Technology Park, Locked Bag 9013, Alexandria, NSW 1435, Australia 115 Institut fur Physik und Astronomie, Universitat Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476, Potsdam, Germany 116 Department of Physics, Penn State Abington, Abington, PA 19001, USA 117 Institute for Theoretical Physics, University of Munster, 48149 Munster, Germany 118 Giant Army, 915A 17th Ave, Seattle WA 98122 119 Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA  \n120 Center for Astrophysics, Harvard University, 60 Garden St, Cambridge, MA 02138  \n- 121 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK 122 Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA\n- 123 Department of Physics, IISER Bhopal, Bhopal Bypass Road, Bhauri, Bhopal 462066, Madhya Pradesh, India  \n124 Ollscoil na Gaillimhe - University of Galway, University Road, Galway, H91 TK33, Ireland  \n- 125 Division of Natural Science, Faculty of Advanced Science and Technology, Kumamoto University, 2-39-1 Kurokami, Kumamoto 860-8555, Japan\n- 126 International Research Organization for Advanced Science and Technology, Kumamoto University, 2-39-1 Kurokami, Kumamoto 860-8555, Japan\n- 127 Laboratoire Univers et Th'eories LUTh, Observatoire de Paris, Universit'e PSL, CNRS, Universit'e de Paris, 92190 Meudon, France 128 CSIRO, Space & Astronomy, P. O. Box 1130, Bentley, WA 6102, Australia\n- 129 Department of Physics, Ben-Gurion University of the Negev, Be'er Sheva 84105, Israel  \n130 Feza Gursey Institute, Bogazici University, Kandilli, 34684, Istanbul, Turkey  \n- 131 Max-Planck-Institut fur Gravitationsphysik (Albert-Einstein-Institut), Callinstrasse 38, D-30167, Hannover, Germany\n- 132 Florida Space Institute, University of Central Florida, 12354 Research Parkway, Partnership 1 Building, Suite 214, Orlando, 32826-0650, FL, USA\n- 133\n- Xinjiang Astronomical Observatory, Chinese Academy of Sciences, Urumqi, Xinjiang 830011, China 134 Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208 135 Adler Planetarium, 1300 S. DuSable Lake Shore Dr., Chicago, IL 60605, USA\n- 136 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), 44780 Bochum, Germany 137 National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100101, People's Republic of China  \n138 Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008, China"}