# The MeerKAT Pulsar Timing Array: The first search for gravitational waves with the MeerKAT radio telescope

Matthew T Miles, Ryan M Shannon, Daniel J Reardon, Matthew Bailes, David J Champion, Marisa Geyer, Pratyasha Gitika, Kathrin Grunthal, Michael J Keith, Michael Kramer, et al.

► To cite this version:

Matthew T Miles, Ryan M Shannon, Daniel J Reardon, Matthew Bailes, David J Champion, et al. The MeerKAT Pulsar Timing Array: The first search for gravitational waves with the MeerKAT radio telescope. Monthly Notices of the Royal Astronomical Society, 2024, 536 (2), pp.1489-1500. 10.1093/mnras/stae2571. hal-04821144

HAL Id: hal-04821144

https://hal.science/hal-04821144v1

Submitted on 5 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The MeerKAT Pulsar Timing Array: the first search for gravitational waves with the MeerKAT radio telescope

Matthew T. Miles<sup>1,2,★</sup> Ryan M. Shannon<sup>1,2</sup> Daniel J. Reardon<sup>1,2</sup> Matthew Bailes,<sup>1,2</sup>

David J. Champion<sup>1,3</sup> Marisa Geyer,<sup>4</sup> Pratyasha Gitika,<sup>1,2</sup> Kathrin Grunthal<sup>1,3</sup> Michael J. Keith<sup>1,5</sup>

Michael Kramer<sup>3,5</sup> Atharva D. Kulkarni<sup>1,2</sup> Rowina S. Nathan<sup>2,6</sup> Aditya Parthasarathy<sup>1,3,7,8</sup>

Jaikhomba Singha<sup>1,4</sup> Gilles Theureau,<sup>9,10</sup> Eric Thrane<sup>1,2,6</sup> Federico Abbate<sup>1,3,11</sup> Sarah Buchner<sup>1,12</sup>

Andrew D. Cameron,<sup>1,2</sup> Fernando Camilo,<sup>12</sup> Beatrice E. Moreschi,<sup>11,13</sup> Golam Shaifullah<sup>1,11,13</sup>

Mohsen Shamohammadi<sup>1,2</sup> Andrea Possenti<sup>11</sup> and Vivek Venkatraman Krishnan<sup>1,3</sup>

Accepted 2024 October 25. Received 2024 September 25; in original form 2024 August 5

## ABSTRACT

Pulsar timing arrays search for nanohertz-frequency gravitational waves by regularly observing ensembles of millisecond pulsars over many years to look for correlated timing residuals. Recently the first evidence for a stochastic gravitational wave background has been presented by the major arrays, with varying levels of significance (~2σ–4σ). In this paper, we present the results of background searches with the MeerKAT Pulsar Timing Array. Although of limited duration (4.5 yr), the ~250 000 arrival times with a median error of just 3 µs on 83 pulsars make it very sensitive to spatial correlations. Detection of a gravitational wave background requires careful modelling of noise processes to ensure that any correlations represent a fit to the underlying background and not other misspecified processes. Under different assumptions about noise processes, we can produce either what appear to be compelling Hellings–Downs correlations of high significance (3σ–3.4σ) with a spectrum close to that which is predicted, or surprisingly, under slightly different assumptions, ones that are insignificant. This appears to be related to the fact that many of the highest precision MeerKAT Pulsar Timing Array pulsars are in close proximity and dominate the detection statistics. The sky-averaged characteristic strain amplitude of the correlated signal in our most significant model is $h_{c,yr} = 7.5_{-0.9}^{+0.8} \times 10^{-15}$ measured at a spectral index of $\alpha = -0.26$, decreasing to $h_{c,yr} = 4.8_{-0.9}^{+0.8} \times 10^{-15}$ when assessed at the predicted $\alpha = -2/3$. These data will be valuable as the International Pulsar Timing Array project explores the significance of gravitational wave detections and their dependence on the assumed noise models.

**Key words:** gravitational waves – methods: data analysis – stars: black holes – pulsars: general.

## 1 INTRODUCTION

The first direct detection of gravitational waves (GWs; Abbott et al. 2016) marked a rare occurrence in physics and astronomy: the

beginning of a new field, and a new way to study the Universe. Almost a century passed between the initial prediction of GWs by Einstein (Einstein 1916) and the observation of perturbed space-time from the coalescence of two black holes 30 times more massive than the Sun. During this time the very existence of GWs had been called into question (at one stage by Einstein himself), and early attempts to detect them had been unsuccessful (Cervantes-

<footnote n="1"><sup>1</sup>Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="2"><sup>2</sup>OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="3"><sup>3</sup>MPI für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="4"><sup>4</sup>High Energy Physics, Cosmology & Astrophysics Theory (HEPCAT) Group, Department of Mathematics and Applied Mathematics, University of Cape Town, Cape Town 7700, South Africa</footnote>

<footnote n="5"><sup>5</sup>Jodrell Bank Centre for Astrophysics, University of Manchester, Alan-Turing Building, Oxford Street, Manchester M13 9PL, UK</footnote>

<footnote n="6"><sup>6</sup>School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Australia</footnote>

<footnote n="7"><sup>7</sup>ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensdijk 4, NL-7991 PD Dwingeloo, the Netherlands</footnote>

<footnote n="8"><sup>8</sup>Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, NL-1098 XH Amsterdam, the Netherlands</footnote>

<footnote n="9"><sup>9</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E UMR7328, Université d'Orléans, CNRS, CNES, OSUC, Observatoire de Paris, F-45071 Orléans, France</footnote>

<footnote n="10"><sup>10</sup>Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris Cité, CNRS, F-92190 Meudon, France</footnote>

<footnote n="11"><sup>11</sup>INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius (CA), Italy</footnote>

<footnote n="12"><sup>12</sup>South African Radio Astronomy Observatory, 2 Fir Street, Black River Park, Observatory 7925, South Africa</footnote>

<footnote n="13"><sup>13</sup>Dipartimento di Fisica 'G. Occhialini', Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote>*E-mail: matthewmiles@swin.edu.au</footnote>
Cota, Galindo-Uribarri & Smoot 2016). The Hulse–Taylor binary pulsar (Hulse & Taylor 1975) had provided compelling indirect evidence that GWs were emitted at the rate implied by the general theory of relativity (Taylor & Weisberg 1982). Decades later, there now exists a catalogue of GW detections that continues to grow (Abbott et al. 2023) for observations limited to the audio band, corresponding to GWs that can be observed at frequencies between ~10 and ~1000 Hz. Pulsar timing arrays (PTAs) offer an alternate window into the GW-bright Universe, sensitive to GWs at nanohertz (nHz) frequencies. Potential sources of these include inspiralling supermassive black hole binaries (SMBHBs; Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2004; Kocsis & Sesana 2011; Ravi et al. 2012; Roedig et al. 2012; Taylor, Simon & Sampson 2017b), cosmic strings (Kibble 1976; Ölmez, Mandic & Siemens 2010; Sanidas, Battye & Stappers 2012; Lentati et al. 2015; Arzoumanian et al. 2018), cosmological phase transitions (Starobinsky 1980; Grishchuk 2005), and quantum fluctuations in the early Universe (Maggiore 2000; Lasky et al. 2016). Of these, GWs from SMBHBs possess the strongest theoretical motivation. Thus, the strongest most probable signal is expected to be the incoherent superposition of all gravitationally radiating, inspiralling SMBHBs in the observable Universe, referred to as the stochastic gravitational wave background (SGWB).

PTAs are regularly observed ensembles of millisecond pulsars (MSPs) in which the arrival times of pulses emitted by pulsars are measured over years to decades (Foster & Backer 1990). Of all pulsars, MSPs are the most rotationally stable, allowing the times of arrival (ToAs) of their pulses to be predicted to a sub-microsecond precision in some systems. The predictability of their emission allows for a wide range of insights, from the properties of pulsar local environments (Wolszczan & Frail 1992; Wolszczan 1994), the nature of the nuclear equation of state (Demorest et al. 2010; Antoniadis et al. 2013; Fonseca et al. 2021; Miller et al. 2021; Riley et al. 2021), to stringent tests of gravity (Stairs 2003; Kramer et al. 2021). It also enables an ensemble of pulsars (a PTA) to be searched for common signals such as the influence of an SGWB. This signal is thought to emerge in PTA data sets both as a noise process with consistent spectral properties in individual pulsar arrival times (the autocorrelated signal) and as a temporally and angularly correlated signal between pulsars that is a function of their angular separation (the cross-correlated signal). The autocorrelated signal is expected to emerge first in the data of most PTAs (Siemens et al. 2013), and a signal that is representative of this has been observed in multiple PTAs (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov et al. 2021; Antoniadis et al. 2022; Xu et al. 2023). However, it has also been shown that such a signal can arise spuriously (Goncharov et al. 2022; Zic et al. 2022). As such, it is only a possible indicator that an SGWB signal may be present in the data, rather than evidence for it. Instead, the signature that is sought after is the cross-correlated signal, which shows angular correlations that are the consequence of the quadrupolar nature of GWs, distinct from other correlated signals in PTA data sets (Tiburzi et al. 2016). The detection of an SGWB is the primary focus of PTA collaborations, with searches having been performed by the European PTA (EPTA; Janssen et al. 2008), the Parkes PTA (PPTA; Manchester et al. 2013), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; Jenet et al. 2009), the Chinese PTA (CPTA; Xu et al. 2023), and also in a joint effort through the International PTA (IPTA; Hobbs et al. 2010). The searches have not yet exceeded IPTA-nominated thresholds to claim a detection (Allen et al. 2023). However, strong evidence (3σ–4σ) has recently been shown for an angularly correlated signal (Agazie et al. 2023a; EPTA Collaboration 2023; Reardon et al.

2023a; Xu et al. 2023), suggesting a definitive detection could be imminent. The primary target of GW searches with PTAs is an isotropic background. However, it is possible that the background may emerge as an anisotropic signal from the influence of bright SMBHBs (Taylor & Gair 2013; Mingarelli et al. 2023). In the era of initial detection of an SGWB, whether it is anisotropic may not be immediately obvious, and identifying this is a goal that could benefit from the combination of data from multiple PTAs by the IPTA. Individual inspiralling SMBHBs can also emit GWs in the nHz-frequency band that can potentially be detected by PTAs. These continuous gravitational waves (CWs) are expected to initially appear in PTA data as non-evolving signals, confined to a single frequency characteristic of their emission. Their evolution occurs over a time-scale much longer than what current PTAs can observe, making them detectable as static signals in PTA data. While an SGWB is thought to be the strongest signal in PTA data, bright individual binaries may also be detectable in tandem with, or soon after, a confirmed detection of an SGWB (Sesana, Vecchio & Volonteri 2009; Rosado, Sesana & Gair 2015).

In this paper, we describe the first searches for GWs with the MeerKAT Pulsar Timing Array (MPTA), reporting evidence for nHz-frequency GWs in a 4.5-yr MPTA data set that appears dependent on the assumptions that we make. In Section 2, we describe the methodology, data set, and noise models used for this analysis. The results of the searches are presented in Section 3. The implications of these are discussed in Section 4, and we state our conclusions in Section 5.

## 2 METHODOLOGY

The data set that is used for this analysis is described in Miles et al. (2024) that extends on from Miles et al. (2023), comprising 4.5 yr of observations recorded with the MeerKAT L-band receiver (856–1712 MHz).

### 2.1 Search for the correlated signature of an SGWB

To determine the presence of all signals in the MPTA data, we use the standard PTA likelihood:

$$ \mathcal{L}(d|\theta) = \frac{\exp\left(-\frac{1}{2}\delta t^T \mathbf{C}^{-1} \delta t\right)}{\sqrt{\det(2\pi\mathbf{C})}}, \qquad (1) $$

where $\delta t$ is a vector of timing residuals, $\mathbf{C}$ is the $N_{\text{TOA}} \times N_{\text{TOA}}$ covariance matrix of the data, and $\theta$ describes the parameters of the model (van Haasteren et al. 2009).

We assume that the GW is described by a power law so that its power spectral density (in the residuals) is of the form

$$ P(f; A, \gamma) = \frac{A^2}{12\pi^2} \left(\frac{f}{f_c}\right)^{-\gamma} \text{yr}^3, \qquad (2) $$

where A is the characteristic strain amplitude of the signal referenced to a frequency of $f_c$, f is GW frequency, and $\gamma$ is the spectral index of the process. We set the reference frequency to be 1 yr$^{-1}$, as is standard in PTA analysis.

The angularly correlated component of this signal must also be modelled, which for an isotropic background is defined by the Hellings–Downs correlation function (Hellings & Downs 1983):

$$ \Gamma_{a,b}(\zeta) = \frac{1}{2} - \frac{1}{4} \left( \frac{1 - \cos \zeta}{2} \right) + \frac{3}{2} \left( \frac{1 - \cos \zeta}{2} \right) \ln \left( \frac{1 - \cos \zeta}{2} \right), \qquad (3) $$
where $\zeta$ is the angular separation of a pair of pulsars, a and b, within the PTA. For simplicity, we define the signal expected from an SGWB as the product of the two, such that

$$P_{HD} = P(f; A, \gamma)\Gamma_{a,b}(\zeta). \quad (4)$$

We assess the presence of this signal in the data through Bayesian, frequentist, and pseudo-Bayesian frameworks. In the frequentist framework, we make use of the optimal statistic (OS; Anholm et al. 2009; Chamberlin et al. 2015). While doing this, we employ noise marginalization using the output of a Bayesian search for an uncorrelated common signal using comprehensive noise models. In Miles et al. (2024), a free spectrum analysis of the common signal identified in the data is assessed, finding that the amplitude measurement is constrained only in the three lowest independent frequencies. Due to this, we undertake the frequentist analyses using only these frequencies, the lowest of which corresponds to $1/T$, where T is the observing span of this MPTA data set. Between $1/T$ and $4/T$, where the PTA is no longer sensitive, these harmonics correspond to a range of frequencies between 7.1 and 28.4 nHz.

We assess the significance of the frequentist analysis by employing bootstrap methods that are commonly used in PTA analyses, known as sky and phase scrambling (Cornish & Sampson 2016; Taylor et al. 2017a). In the Bayesian framework, we employ the Savage-Dickey method (Dickey 1971) to compute a Bayes factor (B) to select the model that is preferred by the data and assess its significance. To complement these analysis streams, we also provide an assessment that acts as an intermediate between both, which we refer to as a pseudo-Bayesian pairwise analysis (Reardon et al. 2023a).

## 2.2 Search for a single source

In addition to searching for correlated signals emerging from an SGWB, we also search over the MPTA data set for the presence of a single GW source that is emitting a non-evolving CW. We construct this model assuming an individual SMBHB, emitting a CW at a single frequency. We search for a circular binary both in the presence and absence of an SGWB. As only the three lowest frequency bins are constrained when assessing the common uncorrelated signal via a free spectrum analysis, if a CW is present in the MPTA data, it is possible that spectral leakage could result in the signal being spread between these frequencies. This might result in difficulties detecting the presence of both an SGWB and CW signal if it existed in the most sensitive frequencies of the array; however, for consistency with analyses performed by other PTAs, we search for both simultaneously.

To model the CW signal, we follow the approaches taken in Arzoumanian et al. (2023) and Agazie et al. (2023b), which we summarize here. The signal of a single GW emitting source can be modelled in pulsar timing data to be

$$s(t) = F^{+}(\theta, \phi, \psi)[s_{+}(t_{p}) - s_{+}(t)] + F^{\times}(\theta, \phi, \psi)[s_{\times}(t_{p}) - s_{\times}(t)]. \quad (5)$$

The two terms, $s_{+,x}(t)$ and $s_{+,x}(t_p)$, are the Earth and pulsar terms of the signal, with + and × representing plus and cross polarization modes of GW emission from general relativity. The F factors in equation (5) are the antenna pattern function of the pulsar response, which is dependent on the location of the single source $(\theta, \phi)$ and the polarization angle of the GW $(\psi)$.

For a circular binary, the terms are

$$s_{+}(t) = \frac{M^{5/3}}{d_{L}\omega(t)^{1/3}} \sin 2\Phi(t)(1 + \cos^{2}i) \quad (6)$$

and

$$s_{\times}(t) = \frac{M^{5/3}}{d_{L}\omega(t)^{1/3}} 2\cos 2\Phi(t) \cos i, \quad (7)$$

where $i$ is the inclination angle of the SMBHB, $d_L$ is the luminosity distance, $M$ is the chirp mass of the binary, $\omega(t)$ is the angular frequency, and $\Phi(t)$ is the GW phase. The angular frequency depends on the physical parameters of the binary and evolves with time, although it is not plausible for this evolution to be detectable within the length of this MPTA data set:

$$\omega(t) = 2\pi f_{GW} \left[ 1 - \frac{256}{5} M^{5/3} (2\pi f_{GW})^{8/3} (t - t_0) \right]^{-3/8}, \quad (8)$$

where $f_{GW}$ is the GW frequency at time $t_0$. A similar transformation can be performed for the phase

$$\Phi(t) = \Phi_0 + \frac{1}{32} M^{-5/3} [(2\pi f_{GW})^{-5/3} - \omega(t)^{-5/3}], \quad (9)$$

where $\Phi_0$ is the initial Earth term phase. Finally, the luminosity distance can be replaced to include an expression describing the overall amplitude of the signal in terms of characteristic strain ($h_c$):

$$d_L = \frac{2M^{5/3}(\pi f_{GW})^{2/3}}{h_c}. \quad (10)$$

These substitutions allow a GW signal to be parametrized with only eight parameters, in conjunction with the intrinsic terms connected to the pulsars in the PTA, and the optional inclusion of a term describing an SGWB. Due to generally imprecise PTA constraints on pulsar distances, we forego the inclusion of the pulsar terms in this analysis. A more complete derivation of this signal, as well as a derivation of the more complex case of an eccentric binary, can be found in Taylor et al. (2016).

## 2.3 Pulsar noise models

To assess the contribution of a GW signal in the data, the comprehensive MPTA noise models described in Miles et al. (2024) are used in a Bayesian and frequentist methodology, consistent with other PTA analyses (Agazie et al. 2023a; EPTA Collaboration 2023), as well as a pseudo-Bayesian method as performed by the PPTA (Reardon et al. 2023a). We assess the presence of a common signal not only under the assumption of these noise models but also including additional achromatic red noise in an effort to mitigate potential misspecification of subthreshold processes that were not favoured for inclusion in Miles et al. (2024). Finally, we assess the presence of the signal using a slightly modified version of this model, in which a noise process that is not well understood is removed from the set of models. This is described further below.

We use a set of subscripts and superscripts to describe the model of intrinsic pulsar noise processes that are used in this assessment: DATA to refer to the model described by the data-derived noise processes defined in Miles et al. (2024), ER to refer to the DATA model with additional achromatic red noise processes included, and ALT to refer to the modestly altered ER model. In the frequentist analyses of these models, we employ a solution derived from the one-dimensional probabilistic maximum for each varied parameter. In some cases, the values chosen using this method are not located in the densest regions of the joint multidimensional posterior distribution. Where this is the case, the values of a process are chosen from the densest regions of the multidimensional posterior that are coincident with the one-dimensional probabilistic maximum of the best-constrained parameter (usually the amplitude).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** A comparison between time realizations of a common achromatic process for two MPTA pulsars, PSR J1909−3744 (top) and PSR J2241−5236 (centre), against the timing residuals of PSR J2129−5721 (bottom). The time realizations of the common achromatic process that we display here appear to be highly correlated, as is expected for pulsars with low angular separations under the assumption of an SGWB. In addition, the realizations appear highly correlated with the unaltered timing residuals of PSR J2129−5721, which is at low angular separations with both of these pulsars.</figcaption></figure>

Hereafter, we refer to this as the FP (from *fixed parameter*) solution. If the posterior of the noise process is not well constrained following an ensemble search (i.e. a Bayesian analysis of the entire PTA), the process is not included in the downstream frequentist analysis.<fnref n="1" /> This is done to mitigate the impact of unfavoured processes on the overall assessment of a common signal. Largely, the processes that are found to be unconstrained during ensemble searches are consistent across the models we consider. In a few cases, a process would no longer be favoured in only one of the assessments. As the properties of a common signal are not significantly different between the models, these cases were monitored carefully and assessed for further inclusion based on the relative support for the process between the single pulsar noise analysis and the ensemble analysis.

The ALT model differs from the ER model by excluding the extraordinarily steep (high spectral index) achromatic red noise process identified in PSR J2129−5721. While this process was not supported for inclusion by the analysis method utilized in Miles et al. (2024), its amplitude is constrained in the ensemble Bayesian assessment of the ER models. This is unexpected, as the PPTA does not report any constraints on an achromatic red noise process for this pulsar over almost two decades of observations (Reardon et al. 2023b). Furthermore, the process is unlike any other found in the MPTA: the constraint on its spectral index extends to $\gamma > 10$, the process is only constrained in its lowest frequency bin, and the amplitude of the process in that bin is heavily degenerate with a common signal across the PTA. Due to the unusual and inconsistent

nature of this noise process, we do not believe that it is well modelled using an achromatic red noise process; instead, it is likely a rare and confounding case of noise misspecification. The unusual nature of this noise is further demonstrated in Fig. 1, in which the unaltered timing residuals of the pulsar are juxtaposed against time realizations of a common achromatic process in two other pulsars that are at low angular separations to PSR J2129−5721, providing circumstantial evidence that the noise processes are common between the pulsars. We provide brief descriptions of these models in Table 1.

## 3 RESULTS

We search for two correlated signals: an isotropic SGWB and a single radiating GW source. For the remainder of the paper, we employ the subscript CURN to refer to a common signal without angular correlations, SMBHB for angular correlations representative of a single radiating source, and HD to represent an angularly correlated signal with correlations expected from an isotropic SGWB.

In Miles et al. (2024), we reported the presence of an uncorrelated common process in the data using the ER model, with the characteristic strain amplitude of the power-law signal measured to be $\log_{10} A_{\text{CURN}}^{\text{ER}} = -14.25_{-0.36}^{+0.21}$ and the spectral index measured to be $\gamma_{\text{CURN}}^{\text{ER}} = 3.60_{-0.89}^{+1.31}$ at a Bayes factor of $\ln(B) = 3.17$ in favour of a CURN in the data. For completeness, we also provide equivalent constraints of the data-driven model, measured at an amplitude $\log_{10} A_{\text{CURN}}^{\text{DATA}} = -14.25_{-0.34}^{+0.21}$ with a spectral index of $\gamma_{\text{CURN}}^{\text{DATA}} = 3.52_{-0.90}^{+1.12}$.

### 3.1 Frequentist assessment of isotropic angular correlations

First, we search for this signal under our model assumptions using the frequentist OS. Using the DATA model, we measure an all-sky angularly correlated signal with an apparent signal-to-noise ratio (S/N) of 4.6. We note that this statistic is not suitable to use as a metric of the significance of the signal, and discuss the interpretation below. For this model, this measurement is equivalent to a characteristic strain amplitude of $\hat{A}_{\text{DATA}}^2 = (5.7 \pm 1.2) \times 10^{-29}$ ($\log_{10} A_{\text{DATA}} = -14.12_{-0.05}^{+0.04}$), assuming $\gamma_{\text{CURN}}^{\text{DATA}} = 3.52$. The presence of Hellings-Downs correlations at these values is apparent when individual pulsar pairs are averaged into bins in pulsar pair separation angle, as shown in Fig. 2, displaying the cross-correlated power in 15 equally distributed bins for the 3403 pulsar pairs in the MPTA. We note that both the OS amplitude and S/N are calculated on a per-pulsar-pair basis rather than from the binned averages presented in Fig. 2. The correlations that are displayed have been calculated using an algorithm described in Allen & Romano (2023), which accounts for the covariances between the pulsar pairs stemming from the similarity of their positions in the sky and the corresponding redundancy of their measurements.

In an effort to be conservative in our investigation of any spatially correlated signals, we also present the correlations for the ER and ALT models. These correlations are calculated in the same manner as the DATA model, using only the marginally different assumptions in the pulsar noise models. Assessing the ER model, we measure an OS amplitude of $\hat{A}_{\text{ER}}^2 = (3.7 \pm 1.9) \times 10^{-29}$, the corresponding correlations of which are presented in Fig. 3(a). For the ALT model, we measure an OS amplitude of $\hat{A}_{\text{ALT}}^2 = (7.1 \pm 1.8) \times 10^{-29}$, a statistically marginal increase in the amplitude resulting from minimal differences between the two models that can be seen in Fig. 3(b).

We assess the significance of the OS after noise marginalization as described in Vigeland et al. (2018). This method has been shown

<footnote n="1">¹In practice this is done by setting the amplitude of the noise process to an immeasurably small value of $\log_{10} A = -20$, which is five to six orders of magnitude smaller than a signal that could be measured in the data set.</footnote>
<table><caption>**Table 1.** Models considered in this work. We provide the number of free parameters that are used in each model and a brief description.</caption><thead><tr><th>Model</th><th>Free parameters</th><th>Model summary</th></tr></thead><tbody><tr><td>DATA</td><td>390</td><td>Pulsar noise processes as described in Miles et al. (2024)</td></tr><tr><td>ER</td><td>532</td><td>An extension to the DATA model with additional achromatic red noise processes when not supported in DATA</td></tr><tr><td>ALT</td><td>530</td><td>The ER model, with the exclusion of the achromatic red noise in PSR J2129-5721</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** All-sky interpulsar correlations of the MPTA data set as a function of the angular separation between pairs of pulsars, displayed in 15 bins across the angular separation range. The $\hat{A}_{\text{DATA}}^2 = (5.65 \pm 1.2) \times 10^{-29}$ corresponds to a dimensionless GW amplitude of $7.5_{-0.9}^{+0.8} \times 10^{-15}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** All-sky interpulsar correlations of the MPTA data set for the two alternative models: ER (left) and ALT (right). Using the ER model, we measure an OS amplitude to be $\hat{A}_{\text{ER}}^2 = (3.7 \pm 1.9) \times 10^{-29}$, and for the ALT model we measure an OS amplitude to be $\hat{A}_{\text{ALT}}^2 = (7.1 \pm 1.8) \times 10^{-29}$. While the models are very similar, there is a marked difference in the correlations that stems primarily from the lack of achromatic red noise in PSR J2129-5721.</figcaption></figure>

to obtain a more reliable measurement of the cross-correlated amplitude. To perform this analysis, we use the posterior samples from a Bayesian analysis which models the intrinsic pulsar noise processes in the MPTA, while also searching for a common uncorrelated noise process with a variable spectral index. We perform this assessment on the DATA and ALT models as the ER model did not provide

a sufficiently high S/N to warrant further frequentist analyses. The distribution of the S/N attained through this marginalization is shown in Fig. 4. In this search, we find a mean S/N and a standard deviation of 0.96 and 0.81 for the DATA model, and a corresponding 0.97 and 0.94 for the ALT model. We note that while the mean values of the noise marginalized OS S/N are considerably lower than the FP value
<figure><img src="image_5.png" /><figcaption>**Figure 4.** S/N distribution attained from the noise marginalized OS. The distribution was calculated from 1000 samples taken from the Bayesian search for a CURN. The mean S/N attained through this method is considerably lower than that attained using the FP values.</figcaption></figure>

S/N, both distributions have long tails towards larger S/Ns where the FP solution resides.

We note that the S/N that is reported through this analysis cannot be confidently used to assess detection probability. Instead, alternate methods are required to empirically derive a significance. We calculate a probability (i.e. a *p*-value), using noise realizations from the pulsar timing data to establish a null space to which the observed correlations can be compared. In standard PTA analyses, a combination of *phase shifts* and *sky scrambling* is employed (Cornish & Sampson 2016; Taylor et al. 2017a).

Assessing the DATA model first, we find that 102 out of ~1.3 × 10⁵ randomized pulsar sky positions result in an S/N that is greater than the value we have measured, corresponding to a *p*-value of approximately 8 × 10⁻⁴. Through phase shifting we found a similar significance at a *p*-value of 6.4 × 10⁻⁴ from 1.9 × 10⁶ samples. Assessing the ALT model in a similar manner, we find *p*-values on the order of 1.4 × 10⁻³ and 1.3 × 10⁻³ through sky and phase scrambling, respectively. The distributions of these bootstrapped methods are shown in Fig. 5 and correspond to ~3.2σ–3.4σ, representing evidence for nHz-frequency GWs assuming the aforementioned FP values. While these values are intriguing, they were measured using the FP values derived from Bayesian analyses, and decrease significantly when we instead apply the mean values from the noise marginalization process. The null distributions displayed in Fig. 5 were also calculated using FP values. Where we assume the median values of the Bayesian analyses to establish the null distributions instead, we find comparable estimates for the significance.

We note that variability between pulsars in achieved timing precision and intrinsic pulsar noise processes reduces the number of *effective* pulsar pairs in a PTA. For the MPTA, taking this into account for the DATA model results in only 608 effective pulsar pairs in contrast to the total 3403, and for the ALT model only 482, impacting the suitability of standard scrambling techniques (Di Marco et al. 2023). To account for this, an additional *super-scramble* technique has been developed (Di Marco et al. 2023) which allows for a greater number of independent realizations. We have not employed this as subsequent work by Di Marco et al. (2024) suggests that correlated scrambles probably produce conservative values even in the presence of misspecified noise.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Distribution of the sky (green) and phase (blue) scrambles calculated from the MPTA data, for the DATA (top) and ALT (bottom) assumptions. The black dashed line corresponds to the S/N calculated from the OS using the FP values, the orange line corresponds to the mean of the noise marginalized S/N distribution. The mean noise marginalized S/N is considerably less significant than that attained using the FP values.</figcaption></figure>

## 3.2 Bayesian assessment of isotropic angular correlations

For both the DATA and the ER models, we perform ensemble Bayesian analyses, searching for a common, achromatic signal possessing HD correlations in the MPTA data set. We use the PTMCMC sampler (Haasteren et al. 2009; van Haasteren & Levin 2013; Lentati et al. 2014; van Haasteren & Vallisneri 2014; Ellis & van Haasteren 2017) to determine the posterior distributions of all time-correlated processes while holding the white noise parameters fixed at the values reported in Miles et al. (2024). We consider both models to demonstrate the suitability of our noise model selection technique, and to investigate if the DATA model is more representative than the conservative ER model. For the DATA model, we estimate an amplitude of $log_{10} A_{\text{HD}}^{\text{DATA}} = -14.28_{-0.30}^{+0.23}$ at a spectral index of $\gamma_{\text{HD}}^{\text{DATA}} = 4.50_{-0.93}^{+1.00}$, consistent with the complementary CURN constraints for this model. The ER model is likewise consistent with the CURN-derived measurements of the common signal in the MPTA data, measured at an amplitude of $log_{10} A_{\text{HD}}^{\text{ER}} = -14.31_{-0.42}^{+0.21}$ and spectral index $\gamma_{\text{HD}}^{\text{ER}} = 3.40_{-0.98}^{+1.45}$.

We assess the evidence for HD correlations by comparing the Savage–Dickey Bayes factors that are calculated for Bayesian analyses assessing a CURN signal, and one assessing a common signal possessing HD correlations. If an HD process is supported, it should have a higher Bayes factor, with the ratio of Bayes factors representing the evidence in support for (or against) Hellings–Downs
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Marginalized posterior distributions of the logarithm of the strain amplitude (log₁₀ A<sub>HD</sub>) and spectral index (γ<sub>HD</sub>) of a Hellings–Downs correlated common signal in the MPTA data, assessed both with the DATA (purple) and ER (gold) models. The contours are the 1σ, 2σ, and 3σ confidence regions of the posterior distributions, and the values reported above the one-dimensional posteriors are the median and 1σ values of the DATA model. The spectral indices of both models are consistent with that expected of an SGWB, which we have overlaid for comparison (grey, dashed).</figcaption></figure>

correlations. In neither case we find support for Hellings-Downs correlations in this MPTA data set. For the ER model, we report ln B<sup>HD+CURN</sup><sub>CURN</sub> = -0.21. For the DATA model, we find that there is not sufficient posterior support at low amplitudes to calculate a direct Bayes factor via the Savage-Dickey method. We can, however, provide limits on the Bayes factors. Comparing the results in this way, we find that both models possess comparable lower limits to the PTA intrinsic red noise (IRN) at ln B<sup>HD+/CURN</sup><sub>IRN</sub> < 13.4. As both models share this lower limit, there is a similar ambiguity in assessing which is favoured by the data.

However, the constraints on the common spectral properties do
not appear any more pronounced when considering HD correlations
as opposed to not. The posteriors of these processes are shown in
Fig. 6.

**3.3 Pseudo-Bayesian assessment of isotropic angular correlations**

The MPTA possesses the largest number of pulsar pairs of any
individual PTA, but the shortest observing span over which to
constrain stochastic processes. Due to this, the *pseudo-Bayesian*
pairwise analysis offers us an opportunity to assess the correlations
that may be present in the MPTA in a manner that could be more
suitable than the other analyses we present here: one that places a
greater emphasis on the correlations between the pairs in the array,
rather than on the stochastic processes themselves.

By sampling the shared correlations and amplitudes between each
of the 3403 pulsar pairs, we can factorize the likelihoods of these
constraints and redistribute them into bins corresponding to their
angular separations on the sky. From this, we are able to calculate

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Factorized correlations of each pulsar pair in the MPTA under the assumptions of the ER model. Through this technique, we recover marginal support for Hellings-Downs-like correlations with the ER model.</figcaption></figure>

a likelihood ratio ($L_h^{HD}$) from the factorized posteriors to determine
if they are more suited to HD correlations than not, following the
technique presented in Reardon et al. (2023a).

The number of pulsar pairs in the MPTA may make this analysis
more prone to error propagation than in other PTA experiments.
While the number of pulsar pairs is likely to make the MPTA
sensitive to interpulsar correlations, the process of factorizing hun-
dreds of posteriors together may result in otherwise insubstantial
uncertainties in the posterior constraints compounding, and result
in non-representative factorized correlations. Given this, we only
assess the conservative ER model, which exhibits the smallest margin
of support for HD correlations in the frequentist analysis. We find
marginal support for HD correlations, at $L_0^{HD} \sim 7$. The factorized
angular correlations derived from this technique are found in
Fig. 7.

**3.4 Searching for a single source**

The ambiguities we observe in our results indicate that the signal
might not be well modelled by an isotropic background. Given this,
we also search for a CW as the progenitor of the common signal.
To evaluate this, we have searched this MPTA data set for a single
SMBHB, using the methods described in Section 2.2. In Fig. 8, we
show the marginalized posterior distributions of a Bayesian analysis
searching for both a CURN and an isolated SMBHB. By utilizing the
hypermodel framework of the PTMCMC sampler, we can measure
a Bayes factor directly. We find that the data strongly favour a model
described by both a CW and a CURN as opposed to only that of a
CW, with a Bayes factor of ln B<sup>SMBHB+CURN</sup><sub>SMBHB</sub> = 9.12. By taking the
Savage-Dickey Bayes factor of both signals, we can also calculate
the degree to which the data prefers the CURN model as opposed to
the SMBHB model. We find this to be ln B<sup>CURN</sup><sub>SMBHB</sub> = 7.76.

While the data do not support the presence of an individual source,
it is worth noting that some posterior parameters show interesting
constraints, particularly the frequency of the GW signal, which may
indicate a preferred emission frequency in the MPTA data set.

**4 DISCUSSION**

In this analysis, we have investigated the presence of nHz-frequency
GWs through the assessment of apparent interpulsar correlations well
described by the Hellings-Downs overlap reduction function.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Marginalized posterior distributions of the CW and CURN signal parameters using the ALT model. We utilize the hypermodel framework of the PTMCMC sampler to naively calculate a Bayes factor between the models, represented by the parameter labelled *nmodel*. The median and $1\sigma$ values of the CW and CURN parameters are displayed above their one-dimensional posterior distributions. We note that a smoothing factor has been applied to the one-dimensional distributions.</figcaption></figure>

**4.1 Independent assessment of nHz-frequency gravitational waves**

The common uncorrelated noise process recovered using the DATA model is consistent with that which is recovered using the ER model in Miles et al. (2024). As discussed in Miles et al. (2024), this is inconsistent with the recovered common noise properties recently published by other PTAs by at least ~1.4σ. It remains unclear if the amplitude measurement of the CURN is physical, or if it is the result of the short MPTA data set as compared to other PTAs. If it is physical, the larger amplitude could be an extension of the

amplitude growth observed in the PPTA data (Reardon et al. 2023a), and implied in the analysis performed by the EPTA + InPTA (EPTA Collaboration 2023). If the amplitude growth is of a physical origin, the relative temporal independence of the data used in this analysis would explain why we have measured a larger amplitude.

The analysis in Section 3 presents an assessment of nHz-frequency GWs with a data set that is largely temporally independent from other searches that have been performed with the techniques utilized in this work. Using the DATA model, the interpulsar correlations recovered for the 83 pulsars imply a dimensionless strain amplitude
of $h_c = 7.5 \times 10^{-15}$, with a reported S/N of 4.6 when fixing the noise models at the FP values. Through the use of sky scrambling and phase shifting methods, we calculate the false alarm probabilities of this measurement to be $p = 8 \times 10^{-4}$ and $6.4 \times 10^{-4}$, respectively.

However, through a noise-marginalized approach, where the co-variance of the pulsar intrinsic noise processes is taken into account in the calculation of the OS, we find that the mean S/N for both the DATA and ER models is ~1. This is considerably lower than the S/N reported using the FP values, and implies that the convergence on the spectral properties of some MPTA noise processes may be too broad to be well described by the FP solution. The similarly broad global constraints of the common signal may also play a part in reducing the mean S/N, as this describes a large parameter space that the several hundred free parameters we consider could find support in.

The short time span that the MPTA possesses results in this signal only being significant in three frequency bins, and the analysis of the signal using the frequentist approach is limited to these. The measurement from the OS analysis, where the FP spectral index is employed, is equivalent to an amplitude of $log_{10} A_{DATA} = -14.12^{+0.04}_{-0.05}$, which is consistent within uncertainties of the DATA model measurement of the CURN, $log_{10} A_{CURN}^{DATA} = -14.25^{+0.21}_{-0.34}$. Assuming that the FP value solution is correct, we conclude that the CURN is dominated by the GW signal and not by misspecified pulsar noise (Zic et al. 2022).

Assuming a spectral index of $13/3$, the value associated with an isotropic SGWB formed of inspiralling SMBHBs, the OS amplitude is measured at $log_{10} A_{DATA} = -14.32^{+0.07}_{-0.08}$, slightly smaller than the FP solution due to the modestly steeper spectral index. Further, the Bayesian analysis of HD correlations in the MPTA data set reports a similar amplitude of the common signal, at $log_{10} A_{HD}^{DATA} = -14.28^{+0.23}_{-0.30}$. The consistency between these measurements, although the Bayes factor between them remains inconclusive, is supportive that the signal could be of a GW origin.

During this analysis, the *pseudo-Bayesian* pairwise technique was also employed. Using our most conservative model, we assessed if an analysis stream that is more dependent on the large number of pulsar pairs in the MPTA would be more sensitive to a correlated signal in the data than other techniques. Evidence for the presence of HD correlations in the data was found at $L_{HD}^{HD} \sim 7$, which could be interpreted as marginal support for the correlations where the Bayesian assessment in Section 3.2 had none. However, there is no clear indication that the correlations returned by this technique are more pronounced, or better constrained, than the frequentist FP value solution for this model.

The temporal independence of the MPTA data set potentially provides further support for the presence of detectable gravitational radiation in the nHz-frequency band. Naively it may be expected that MPTA pulsars, and therefore our data set, are also relatively independent as many are at southern declinations as opposed to the largely northern declination pulsars observed by most other PTAs. However, the most precisely timed pulsars will influence these results the most, and many of these (such as PSR J1909−3744) are observed by several PTAs.

For an isotropic background of GWs emitted by SMBHBs, the amplitude that we have measured is thought to be largely a function of the masses of the SMBHBs that make up the ensemble distribution (Sesana 2013; McWilliams, Ostriker & Pretorius 2014). If this is indeed the source of the signal we observe, the amplitude that we have measured is within the 95 per cent confidence limit of models presented in Sesana (2013), when we consider the lowest frequency resolved by the MPTA data set. However, recent work by Somalwar & Ravi (2023) has presented that, in simulations where SMBH seeds

are able to establish their mass early, the predicted amplitude of the background may be larger by as much as a factor of 3. Given this, it may not be unexpected to measure the amplitude of the SGWB that we observe.

## 4.2 Under the assumption of a single source

The ambiguity of the analyses in this data set led us to consider if a single source emitting CWs could explain the apparent angular correlations. A search was performed on the MPTA data for a circularly inspiralling SMBHB emitting CWs; however, there was a clear preference for a model that not only contained a continuous waveform, but one that also included a common stochastic signal in the form of a CURN. The marginalized posterior distribution of this model in Fig. 8 shows a clear preference for the CURN + CW model, and the posterior distribution is unconstrained for most of the CW model parameters. Of note, however, is the marginal support for a CW at higher frequencies shown in the one-dimensional posterior distribution of the GW frequency ($log_{10} f_{GW}$). This constraint also appears to be covariant with the sky-location CW parameter distributions, certainly motivating a more thorough assessment of alternate explanations of a correlated signal other than an isotropic GW background (see Grunthal et al. 2024).

## 4.3 Comparison to other PTAs

In our search for a GW background using the OS, we find the characteristic strain amplitude in the MPTA data is $log_{10} A_{DATA} = -14.12^{+0.04}_{-0.05}$ assuming FP values of the common signal in the data. In Miles et al. (2024), we discussed that fixing the spectral index of the common uncorrelated process at the nominal value of $13/3$, associated with an isotropic SGWB formed from an ensemble of inspiralling SMBHBs, resulted in a minimum deviation of $1.35\sigma$ as compared to other PTA experiments.

Assessing the OS amplitude under these assumptions, we measure a characteristic strain log-amplitude of $log_{10} A_{DATA} = -14.32^{+0.07}_{-0.08}$, consistent with our previous result within $1\sigma$ uncertainties. However, we acknowledge the apparent sensitivity of the OS amplitude to the parameter values and models that are chosen, which can be observed between Figs 2 and 3. The Bayesian assessments of a common signal with HD correlations do not appear to share this dependence between our choice of models (Fig. 6), and so we compare this instead to the results of other PTAs. Isolating the Bayesian search of the DATA model to a small window (0.1) around $\gamma = 4.333$, we measure a log-amplitude of $log_{10} A_{DATA} = -14.39^{+0.11}_{-0.11}$, consistent with our OS measurement but centred at a marginally lower value.

The EPTA + InPTA collaboration measured the amplitude of the signal in their data at a fixed spectral index using the OS for several versions of their data set. Of these measurements, the largest is $log_{10} A = -14.48^{+0.08}_{-0.11}$ (EPTA Collaboration 2023), for a data set including only their most recent 9 yr of data. The PPTA did not use an OS analysis in their most recent search, so we compare to the strain amplitude measured in their pairwise search for an angularly correlated signal, which they report to be $log_{10} A = -14.69^{+0.05}_{-0.05}$ (Reardon et al. 2023a). While NANOGrav used an OS approach, they did not report a strain amplitude from this assessment. Instead, they report a strain amplitude from their Bayesian analysis of $log_{10} A = -14.62^{+0.11}_{-0.12}$ (Agazie et al. 2023a). The MPTA measurement is larger in amplitude than all values reported by other PTAs. In comparison to the PTA that reported the next largest amplitude, the EPTA + InPTA, our measured amplitude
overlaps within reported uncertainties. However, in the case of the PPTA, which reported the smallest amplitude, we are at a far greater $2.5\sigma$ discrepancy.

While the differences between the amplitude we present here and that taken by the PPTA are the largest, the PPTA also reported an apparent growth in the strain amplitude of the common signal found in their data. This was identified by splitting their entire data set into overlapping 9 and 6 yr intervals, which they referred to as slices. The strain amplitude in their most recent 6 yr time slice was measured to be $\log_{10} A = -14.57_{-0.13}^{+0.11}$, just $1.2\sigma$ away from the MPTA measurement. This could suggest that shorter data sets may result in signal measurements with larger amplitudes, however, there is a relative growth in the signal amplitude through the equal slices of data reported in Reardon et al. (2023a), implying a temporal dependence which may be affecting the MPTA measurements.

The PTA data analysis with the most temporal overlap with the MPTA is the Chinese Pulsar Timing Array Data Release 1 (Xu et al. 2023). While they provide a constraint on the CURN in their analysis (which is discussed in Miles et al., 2024), they use an alternate statistic to assess the presence of cross-correlated signals that do not provide an estimate of the characteristic strain amplitude. It is not possible to directly compare our results to theirs.

## 4.4 Interpreting the MPTA results

Using methods that are standard practice in PTA experiments, we have reported evidence of a correlated signal in a 4.5-yr MPTA data set, with interpulsar correlations indicative of a GW background. The signal is found, using data-informed models describing astrophysically motivated noise processes, is visually apparent across the angular separation range, and has a significance equivalent to ~3.4σ. When using a more conservative model, we find a comparable significance of ~3.2σ. An assessment of the signal under both of these model assumptions results in visually apparent HD correlations when using FP values, providing confidence that the signal may be of GW origin. However, assessing the most conservative model that we consider (ER) does not result in apparent correlations.

While this result is likely to be statistically robust, especially as the background-estimation techniques used to derive the significance are thought to be conservative (Di Marco et al. 2024), other aspects of our analysis do not show the same support for the signal in the data. Marginalizing over correlated noise constraints diminishes the significance of the signal, and the results from Bayesian ensemble searches for a common signal possessing HD correlations are similarly insignificant. The difference in significance between our results is notable, even if the constraints on the amplitude measurements are consistent with each other.

The cause of these differences is not well understood. It may be that the short data span of the MPTA places a heavy toll on our GW searches and noise analysis. Each of the methods we employ here is predicated on some form of Bayesian analysis that requires the ability to precisely and accurately constrain the noise processes. While the sensitivity to these processes should increase quickly for the MPTA given the comparatively greater timing precision than other PTA experiments, the processes still require several significant frequencies to be confidently constrained. For many pulsars in the MPTA data set, this is likely not yet achieved, leading to broad posterior constraints. There is also considerable covariance between different families of noise models (e.g. dispersion measure noise and other chromatic noise with some alternate chromaticity), which can lead to further confusion in the analysis.

The broad constraints on the noise processes result in many possible combinations of parameter values that well model both the intrinsic pulsar noise and any common noise process across the array. In many of these combinations, the parameter values will differ from our chosen values, some more so than others. In these cases, it can be expected that the description of the intrinsic pulsar noise does not provide probabilistic support for a Hellings-Downs correlated signal. Naturally, a PTA that has a greater number of noise parameters that are prone to weaker constraints, such as the MPTA, will be more susceptible to this.

The variability of our results with different noise models highlights a promising area for future study. For future work, we envision a single, unified noise model with sufficient flexibility to describe the noise in MPTA pulsars without artificially inflating error bars with overly conservative assumptions. Such a model could employ hierarchical modelling to model the population properties of pulsars (Goncharov et al. 2022; van Haasteren 2024). It could marginalize over multiple submodels, some of which include additional noise processes, e.g. achromatic red noise. In this way, we can use the data to help determine the model complexity. The ambiguity between model selection, and the desire for conservatism when constructing said models, may mean the correct approach is to marginalize over all models we consider. The strategy we have employed in the design of our models is perhaps a crude analogue to this, which could naturally be developed to involve all unmodelled processes across a full PTA analysis. We leave this for future work.

While the specific issues with this analysis are not well understood, it should be possible to resolve the nature of any signal present in the MPTA data set soon. First, the ability to discriminate signals should increase rapidly as MPTA observations continue. The strength of red noise processes (both a GW background and the processes most covariant with it) increases quickly with time, relative to white noise processes associated with ToA measurement uncertainty and pulse jitter noise (Haasteren & Levin 2013). Using the methodology outlined in Siemens et al. (2013), Spiewak et al. (2022) demonstrated that for a background comparable in amplitude to that which we measure here, the MPTA is currently well described by the moderate or weak signal regimes. As such, even modest increases in the length of MPTA data sets are fractionally large compared to the other, more established, PTA experiments.

Secondly, it is possible that the results presented here motivate the need for new pulsar timing analysis methods that are better able to assess the presence of a stochastic background in the presence of a cacophony of noise processes. While power-law Gaussian processes are commonly used to model stochastic processes in pulsar timing analysis, it is unclear if they are appropriate for both achromatic and chromatic noise (e.g. Cordes & Downs 1985; Coles et al. 2015). We see evidence for this through the identification of Gaussian and annual variations in some of the pulsars in our sample. Similarly, population-level inference of noise properties might be able to reduce covariance between intrinsic spin noise, chromatic noise, and any GW signal (Di Marco et al. 2024).

Lastly, the nature of the signals will benefit through comparison and combination with other PTA data sets. Many of the best-timed pulsars in the MPTA are observed by other PTA experiments. The combination of the data analysed here and presented in Miles et al. (2024) in the planned IPTA third data release will greatly assist the efforts of all collaborations in determining the nature of the various stochastic signals.
## 5 CONCLUSIONS

In this analysis, we have presented an assessment of the presence of nanohertz-frequency GWs in the MPTA, using a data set largely independent from other PTA experiments, sensitive at different GW frequencies. Using the data from 83 pulsars in the MPTA, observations of which were collected over a period of 4.5 yr from the MeerKAT radio telescope, we find an angularly correlated signal consistent with the expected Hellings-Downs correlations. Using our most significant model, the amplitude of the cross-correlated signal is estimated to be $A_{GW} = 7.5_{-0.9}^{+0.8} \times 10^{-15}$ measured at the FP spectral index of the process, and $A_{GW} = 4.8_{-0.9}^{+0.8} \times 10^{-15}$ where measured at the predicted spectral index of a SMBHB-dominated SGWB, $\gamma = 13/3$. The S/N measured at FP values for the presence of this signal in the data is 4.6, with false alarm probabilities corresponding to $p = 8 \times 10^{-4}$ and $6.4 \times 10^{-4}$, calculated using the background estimation sky and phase scrambling techniques, respectively. We also undertake a search for a single radiating source in the form of a loud SMBHB, however, we find it is not supported compared to a model including a common uncorrelated signal. The false alarm probabilities that we report under the assumption of an isotropic background correspond to a significance of $\sim 3.4\sigma$, although it is important to note that these results are particularly sensitive to assumptions made in the construction of the intrinsic pulsar noise models. Further, the application of practices such as noise marginalization and Bayesian ensemble analyses used in the assessment of the MPTA results greatly diminishes their significance, implying that the data are ambiguous to the presence of Hellings-Downs correlations. Given these caveats, interpreting the nature of the apparent correlated signal in the MPTA data is non-trivial. Further work using greater amounts of MPTA data, and in combination with IPTA data, is required to build confidence in our results. Nevertheless, our results demonstrate the ability of sensitive, high-cadence PTA experiments to broaden the nanohertz-frequency GW window, and enable access to new GW sources and GW source classes, while also increasing the sensitivity of international efforts.

## ACKNOWLEDGEMENTS

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory (SARAO), which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time-space reference systems department of the Paris Observatory. MeerTime data are stored and processed on the OzStar and Ngarrgu Tindebeek supercomputers, operated by the Swinburne University of Technology. This work was undertaken as part of the Australian Research Council (ARC) Centre of Excellence for Gravitational Wave Discovery (CE170100004 and CE230100016). RMS acknowledges support through ARC Future Fellowship FT190100155. GT acknowledges financial support from 'Programme National de Cosmologie et Galaxies' (PNCG), 'Programme National Hautes Energies' (PNHE), and 'Programme National Gravitation, Références, Astronomie, Métrologie' funded by INSU,CNRS-IN2P3-INP, CEA, and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. MK acknowledges support from the MPG and the CAS-MPG Legacy Programme. KG acknowledges continuing valuable support from the Max-Planck Society. KG also acknowledges the support from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the

Universities of Bonn and Cologne. VVK acknowledges financial support from the European Research Council (ERC) starting grant 'COMPACT' (grant agreement number 101078094). FA acknowledges that part of the research activities described in this paper were carried out with the contribution of the NextGenerationEU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 - Education and Research, Component 2 - From Research to Business (M4C2), Investment Line 3.1 - Strengthening and Creation of Research Infrastructures, Project IR0000034 - 'STILES - Strengthening the Italian Leadership in ELT and SKA'. Pulsar research at Jodrell Bank Centre for Astrophysics was supported by an STFC consolidated grants (ST/T000414/1 and ST/X001229/1). JS acknowledges funding from the South African Research Chairs Initiative of The Department of Science and Technology and the National Research Foundation of South Africa. PG acknowledges support through the Swinburne University of Technology (SUT) stipend, Swinburne University Postgraduate Research Award (SUPRA). AP acknowledges financial support from the European Research Council (ERC) starting grant 'GIGA' (grant agreement no. 101116134) and through the NWO-I Veni fellowship. Funding was provided for the PTUSE machines by INAF - Osservatorio Astronomico di Cagliari under the Grant 'iPeska' (PI: A. Possenti). Funding was provided for the PTUSE machines by the Max-Planck-Institut für Radioastronomie (MPIfR), also supported by the MPG-CAS LEGACY programme, Swinburne University of Technology, and the Australian SKA office. This project utilized the MeerTime data portal, which was supported by Nick Swainston and the ADACS team. We acknowledge and pay respects to the Elders and Traditional Owners of the land on which the Australian institutions stand, the Bunurong and Wurundjeri Peoples of the Kulin Nation. We thank Rutger van Haasteren, Steve Taylor, Siyuan Chen, and Levi Schult for their useful comments and discussions.

## DATA AVAILABILITY

All data used in this work are available courtesy of AAO Data Central (https://datacentral.org.au/) at https://doi.org/10.57891/j0vh-5g31. The data provided includes sub-banded ToAs, the full data archives used to construct this data release, and the ephemerides that have been used to perform timing. Also included are the frequency-resolved portraits used to calculate the ToAs used for this work.

The archives and portraits are in PSRFITS file format. The ephemerides are in a standard *ascii* text file format, and the arrival times are supplied as IFF data.

## REFERENCES

* Abbott B. P. et al., 2016, Phys. Rev. Lett., 116, 061102

* Abbott R. et al., 2023, Phys. Rev. X, 13, 041039

* Agazie G. et al., 2023a, ApJ, 951, L8

* Agazie G. et al., 2023b, ApJ, 951, L50

* Allen B., Romano J. D., 2023, Phys. Rev. D, 108, 043026

* Allen B., Dhurandhar S., Gupta Y., McLaughlin M., Natarajan P., Shannon R. M., Thrane E., Vecchio A., 2023, preprint (arXiv:2304.04767)

* Anholm M., Ballmer S., Creighton J. D. E., Price L. R., Siemens X., 2009, Phys. Rev. D, 79, 084030

* Antoniadis J. et al., 2013, Science, 340, 448

* Antoniadis J. et al., 2022, MNRAS, 510, 4873

* Arzoumanian Z. et al., 2018, ApJ, 859, 47

* Arzoumanian Z. et al., 2020, ApJ, 905, L34

* Arzoumanian Z. et al., 2023, ApJ, 951, L28

* Cervantes-Cota J., Galindo-Uribarri S., Smoot G., 2016, Universe, 2, 22
Chamberlin S. J., Creighton J. D. E., Siemens X., Demorest P., Ellis J., Price
L. R., Romano J. D., 2015, Phys. Rev. D, 91, 044048

Chen S. et al., 2021, MNRAS, 508, 4970

Coles W. A. et al., 2015, ApJ, 808, 113

Cordes J. M., Downs G. S., 1985, ApJS, 59, 343

Cornish N. J., Sampson L., 2016, Phys. Rev. D, 93, 104047

Demorest P. B., Pennucci T., Ransom S. M., Roberts M. S. E., Hessels J. W.
T., 2010, Nature, 467, 1081

Dickey J. M., 1971, Ann. Math. Stat., 42, 204

Di Marco V., Zic A., Miles M. T., Reardon D. J., Thrane E., Shannon R. M., 2023, ApJ, 956, 14

Di Marco V., Zic A., Shannon R. M., Thrane E., 2024, MNRAS, 532, 4026

Einstein A., 1916, Ann. Phys., 354, 769

Ellis J., van Haasteren R., 2017, jellis18/PTMCMCSampler: Official Re-
lease (1.0.0). Zenodo. Available at: https://doi.org/10.5281/zenodo.103
7579

EPTA Collaboration, 2023, A&A, 678, A50

Fonseca E. et al., 2021, ApJ, 915, L12

Foster R. S., Backer D. C., 1990, ApJ, 361, 300

Goncharov B. et al., 2021, ApJ, 917, L19

Goncharov B. et al., 2022, ApJ, 932, L22

Grishchuk L. P., 2005, Phys. Usp., 48, 1235

Grunthal K. et al., 2024, MNRAS, 536, 1501

Hellings R. W., Downs G. S., 1983, ApJ, 265, L39

Hobbs G. et al., 2010, Class. Quantum Grav., 27, 084013

Hulse R. A., Taylor J. H., 1975, ApJ, 195, L51

Jaffe A. H., Backer D. C., 2003, ApJ, 583, 616

Janssen G. H., Stappers B. W., Kramer M., Purver M., Jessner A., Cognard
I., 2008, in Bassa C., Wang Z., Cumming A., Kaspi V. M., eds, AIP Conf.
Proc. Vol. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and
More. Am. Inst. Phys., New York, p. 633

Jenet F. et al., 2009, preprint (arXiv:0909.1058)

Kibble T. W. B., 1976, J. Phys. A: Math. Gen., 9, 1387

Kocsis B., Sesana A., 2011, MNRAS, 411, 1467

Kramer M. et al., 2021, Phys. Rev. X, 11, 041050

Lasky P. D. et al., 2016, Phys. Rev. X, 6, 011035

Lentati L., Alexander P., Hobson M. P., Feroz F., van Haasteren R., Lee K.
J., Shannon R. M., 2014, MNRAS, 437, 3004

Lentati L. et al., 2015, MNRAS, 453, 2576

McWilliams S. T., Ostriker J. P., Pretorius F., 2014, ApJ, 789, 156

Maggiore M., 2000, preprint (arXiv:gr-qc/0008027)

Manchester R. N. et al., 2013, Publ. Astron. Soc. Aust., 30, e017

Miles M. T. et al., 2023, MNRAS, 519, 3976

Miles M. T. et al., 2024, MNRAS, 536, 1467

Miller M. C. et al., 2021, ApJ, 918, L28

Mingarelli C. M. F. et al., 2017, Nat. Astron., 1, 886

Ölmez S., Mandic V., Siemens X., 2010, Phys. Rev. D, 81, 104028

Rajagopal M., Romani R. W., 1995, ApJ, 446, 543

Ravi V., Wyithe J. S. B., Hobbs G., Shannon R. M., Manchester R. N., Yardley
D. R. B., Keith M. J., 2012, ApJ, 761, 84

Reardon D. J. et al., 2023a, ApJ, 951, L6

Reardon D. J. et al., 2023b, ApJ, 951, L7

Riley T. E. et al., 2021, ApJ, 918, L27

Roedig C., Sesana A., Dotti M., Cuadra J., Amaro-Seoane P., Haardt F., 2012, A&A, 545, A127

Rosado P. A., Sesana A., Gair J., 2015, MNRAS, 451, 2417

Sanidas S. A., Battye R. A., Stappers B. W., 2012, Phys. Rev. D, 85, 122003

Sesana A., 2013, MNRAS, 433, L1

Sesana A., Haardt F., Madau P., Volonteri M., 2004, ApJ, 611, 623

Sesana A., Vecchio A., Volonteri M., 2009, MNRAS, 394, 2255

Siemens X., Ellis J., Jenet F., Romano J. D., 2013, Class. Quantum Grav., 30,
224015

Somalwar J. J., Ravi V., 2023, preprint (arXiv:2306.00898)

Spiewak R. et al., 2022, Publ. Astron. Soc. Aust., 39, e027

Stairs I. H., 2003, Living Rev. Relativ., 6, 5

Starobinsky A. A., 1980, Phys. Lett. B, 91, 99

Taylor S. R., Gair J. R., 2013, Phys. Rev. D, 88, 084001

Taylor J. H., Weisberg J. M., 1982, ApJ, 253, 908

Taylor S. R., Huerta E. A., Gair J. R., McWilliams S. T., 2016, ApJ, 817, 70

Taylor S. R., Lentati L., Babak S., Brem P., Gair J. R., Sesana A., Vecchio
A., 2017a, Phys. Rev. D, 95, 042002

Taylor S. R., Simon J., Sampson L., 2017b, Phys. Rev. Lett., 118, 181102

Tiburzi C. et al., 2016, MNRAS, 455, 4339

van Haasteren R., 2024, ApJS, 273, 23

van Haasteren R., Levin Y., 2013, MNRAS, 428, 1147

van Haasteren R., Vallisneri M., 2014, Phys. Rev. D, 90, 104012

van Haasteren R., Levin Y., McDonald P., Lu T., 2009, MNRAS, 395, 1005

Vigeland S. J., Islo K., Taylor S. R., Ellis J. A., 2018, Phys. Rev. D, 98,
044003

Wolszczan A., 1994, Science, 264, 538

Wolszczan A., Frail D. A., 1992, Nature, 355, 145

Wyithe J. S. B., Loeb A., 2003, ApJ, 590, 691

Xu H. et al., 2023, Res. Astron. Astrophys., 23, 075024

Zic A. et al., 2022, MNRAS, 516, 410

This paper has been typeset from a TeX/Latex file prepared by the author.