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

HAL
open science

# Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array

M Falxa, S Babak, P.T Baker, B Bécsy, A Chalumeau, S Chen, Z Chen, N.J Cornish, L Guillemot, J.S Hazboun, et al.

► To cite this version:

M Falxa, S Babak, P.T Baker, B Bécsy, A Chalumeau, et al. Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array. Mon.Not.Roy.Astron.Soc., 2023, 521 (4), pp.5077-5086. 10.1093/mnras/stad812. hal-04055845

HAL Id: hal-04055845

https://hal.science/hal-04055845

Submitted on 22 Apr 2024

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.
# Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array

M Falxa, S Babak, P.T Baker, B Bécsy, A Chalumeau, S Chen, Z Chen, N.J Cornish, Lucas Guillemot, J.S Hazboun, et al.

► To cite this version:

M Falxa, S Babak, P.T Baker, B Bécsy, A Chalumeau, et al. Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array. Monthly Notices of the Royal Astronomical Society, 2023, 521 (4), pp.5077-5086. 10.1093/mnras/stad812. hal-04055845

HAL Id: hal-04055845

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

Submitted on 22 Apr 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array

M. Falxa,¹★ S. Babak,¹ P. T. Baker,² B. Bécsy,³ A. Chalumeau,⁴ S. Chen,⁵ Z. Chen,⁶ N. J. Cornish,³ L. Guillemot,⁷,⁸ J. S. Hazboun,⁹ C. M. F. Mingarelli,¹⁰,¹¹ A. Parthasarathy,¹² A. Petiteau,¹,¹³ N. S. Pol,¹⁴ A. Sesana,⁴,¹⁵ S. B. Spolaor,¹⁶ S. R. Taylor,¹⁴ G. Theureau,⁷,⁸,¹⁷ M. Vallisneri,¹⁸ S. J. Vigeland,¹⁹ C. A. Witt,²⁰,²¹ X. Zhu,²² J. Antoniadis,¹²,²³ Z. Arzoumanian,²⁴ M. Bailes,²⁵,²⁶ N. D. R. Bhat,²⁷ L. Blecha,²⁸ A. Brazier,²⁹ P. R. Brook,¹⁶ N. Caballero,²⁵,⁵ A. D. Cameron,³⁰ J. A. Casey-Clyde,¹¹ D. Champion,¹² M. Charisi,¹⁴ S. Chatterjee,²⁹ I. Cognard,⁷,⁸ J. M. Cordes,²⁹ F. Crawford,³¹ H. T. Cromartie,²⁹ K. Crowter,³² S. Dai,³³ M. E. DeCesar,³⁴ P. B. Demorest,³⁵ G. Desvignes,¹² T. Dolch,³⁶,³⁷ B. Drachler,³⁸ Y. Feng,³⁹ E. C. Ferrara,⁴⁰ W. Fiore,¹⁶ E. Fonseca,⁴¹ N. Garver-Daniels,¹⁶ J. Glaser,⁴² B. Goncharov,⁴³ D. C. Good,¹¹,¹⁰ J. Griessmeier,⁷ Y. J. Guo,¹² K. Gültekin,⁴⁴ G. Hobbs,⁴⁵ H. Hu,¹² K. Islo,¹⁹ J. Jang,¹² R. J. Jennings,¹⁶ A. D. Johnson,⁴⁶ M. L. Jones,¹⁹ J. Kaczmarek,⁴⁷ A. R. Kaiser,¹⁶ D. L. Kaplan,¹⁹ M. Keith,⁴⁸ L. Z. Kelley,²⁰ M. Kerr,⁴⁹ J. S. Key,⁵⁰ N. Laal,⁵¹ M. T. Lam,³⁸ W. G. Lamb,¹⁴ T. J. W. Lazio,¹⁸ K. Liu,¹² T. Liu,¹⁶ J. Luo,⁵² R. S. Lynch,⁵³ D. R. Madison,⁵⁴ R. Main,¹² R. Manchester,⁴⁵ A. McEwen,¹⁹ J. McKee,⁵⁵,⁵⁶ M. A. McLaughlin,¹⁶ C. Ng,⁵⁷ D. J. Nice,⁵⁸ S. Ocker,²⁹ K. D. Olum,⁵⁹ S. Osłowski,⁶⁰ T. T. Pennucci,⁶¹ B. B. P. Perera,⁶² D. Perrodin,⁶³ N. Porayko,¹² A. Possenti,⁶³ H. Quelquejay-Leclere,¹ S. M. Ransom,⁶⁴ P. S. Ray,⁶⁵ D. J. Reardon,³⁰ C. J. Russell,⁶⁶ A. Samajdar,⁴ J. Sarkissian,⁶⁷ L. Schult,¹⁴ G. Shaifullah,⁴,¹⁵,⁶³ R. M. Shannon,³⁰ B. J. Shapiro-Albert,⁶⁸ X. Siemens,⁹ J. J. Simon,⁶⁹ M. Siwek,⁷⁰ T. L. Smith,⁷¹ L. Speri,⁷² R. Spiewak,⁴⁸ I. H. Stairs,³² B. Stappers,⁴⁸ D. R. Stinebring,⁷³ J. K. Swiggum,⁵⁸ C. Tiburzi,⁶³ J. Turner,¹⁶ A. Vecchio,⁷⁴ J. P. W. Verbiest,⁷⁵,¹² H. Wahl,¹⁶ S. Q. Wang,⁷⁶ J. Wang,⁷⁷ J. Wang,⁷⁵ Z. Wu,⁷⁵ L. Zhang,⁷⁸ and S. Zhang⁷⁹ (IPTA Collaboration)

*Affiliations are listed at the end of the paper*

Accepted 2023 February 24. Received 2023 February 17; in original form 2022 November 16

## ABSTRACT

The International Pulsar Timing Array 2nd data release is the combination of data sets from worldwide collaborations. In this study, we search for continuous waves: gravitational wave signals produced by individual supermassive black hole binaries in the local universe. We consider binaries on circular orbits and neglect the evolution of orbital frequency over the observational span. We find no evidence for such signals and set sky averaged 95 per cent upper limits on their amplitude $h_{95}$. The most sensitive frequency is 10 nHz with $h_{95} = 9.1 \times 10^{-15}$. We achieved the best upper limit to date at low and high frequencies of the PTA band thanks to improved effective cadence of observations. In our analysis, we have taken into account the recently discovered common red noise process, which has an impact at low frequencies. We also find that the peculiar noise features present in some pulsars data must be taken into account to reduce the false alarm. We show that using custom noise models is essential in searching for continuous gravitational wave signals and setting the upper limit.

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

## 1 INTRODUCTION

The goal of the Pulsar Timing Array (PTA) collaborations is to detect gravitational wave (GW) signals in the nanohertz band,

where we expect to see a gravitational wave background (GWB) produced by the superposition of GW signals from the population of supermassive black hole binaries (SMBHBs) (Jaffe & Backer 2003; Sesana, Vecchio & Colacino 2008; Maiorano, Paolis & Nucita 2021). Some individual SMBHBs might be brighter than most and stand above the stochastic signal; those are individually resolvable sources (Sesana, Vecchio & Volonteri 2009; Rosado, Sesana & Gair

<footnote>* E-mail: falxa@apc.in2p3.fr</footnote>
# Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array

M. Falxa<sup>1,★</sup> S. Babak,<sup>1</sup> P. T. Baker,<sup>2</sup> B. Bécsy,<sup>3</sup> A. Chalumeau,<sup>4</sup> S. Chen<sup>5</sup>,<sup>6</sup> Z. Chen<sup>6</sup>,<sup>7</sup> N. J. Cornish,<sup>3</sup>
L. Guillemot,<sup>7,8</sup> J. S. Hazboun<sup>9</sup>,<sup>10,11</sup> C. M. F. Mingarelli<sup>10</sup>,<sup>12</sup> A. Parthasarathy<sup>12</sup>,<sup>13</sup> A. Petiteau,<sup>1,13</sup>
N. S. Pol,<sup>14</sup> A. Sesana,<sup>4,15</sup> S. B. Spolaor,<sup>16</sup> S. R. Taylor,<sup>14</sup> G. Theureau,<sup>7,8,17</sup> M. Vallisneri<sup>16</sup>,<sup>18</sup>
S. J. Vigeland,<sup>19</sup> C. A. Witt,<sup>20,21</sup> X. Zhu<sup>16</sup>,<sup>22</sup> J. Antoniadis,<sup>12,23</sup> Z. Arzoumanian,<sup>24</sup> M. Bailes,<sup>25,26</sup>
N. D. R. Bhat<sup>16</sup>,<sup>27</sup> L. Blecha,<sup>28</sup> A. Brazier,<sup>29</sup> P. R. Brook,<sup>16</sup> N. Caballero<sup>5</sup>,<sup>30</sup> A. D. Cameron<sup>16</sup>,<sup>30</sup>
J. A. Casey-Clyde,<sup>11</sup> D. Champion<sup>12</sup>,<sup>13</sup> M. Charisi,<sup>14</sup> S. Chatterjee,<sup>29</sup> I. Cognard,<sup>7,8</sup> J. M. Cordes,<sup>29</sup>
F. Crawford<sup>31</sup>,<sup>32</sup> H. T. Cromartie,<sup>29</sup> K. Crowter,<sup>32</sup> S. Dai<sup>33</sup>,<sup>34</sup> M. E. DeCesar,<sup>34</sup> P. B. Demorest,<sup>35</sup>
G. Desvignes<sup>12</sup>,<sup>36,37</sup> B. Drachler,<sup>38</sup> Y. Feng,<sup>39</sup> E. C. Ferrara,<sup>40</sup> W. Fiore,<sup>16</sup> E. Fonseca,<sup>41</sup>
N. Garver-Daniels,<sup>16</sup> J. Glaser,<sup>42</sup> B. Goncharov<sup>43</sup>,<sup>44</sup> D. C. Good,<sup>11,10</sup> J. Griessmeier,<sup>7</sup> Y. J. Guo,<sup>12</sup>
K. Gültekin<sup>44</sup>,<sup>45</sup> G. Hobbs<sup>45</sup>,<sup>46</sup> H. Hu<sup>12</sup>,<sup>47</sup> K. Islo,<sup>19</sup> J. Jang,<sup>12</sup> R. J. Jennings,<sup>16</sup> A. D. Johnson,<sup>46</sup>
M. L. Jones,<sup>19</sup> J. Kaczmarek,<sup>47</sup> A. R. Kaiser,<sup>16</sup> D. L. Kaplan<sup>19</sup>,<sup>48</sup> M. Keith<sup>48</sup>,<sup>49</sup> L. Z. Kelley<sup>16</sup>,<sup>50</sup>
M. Kerr<sup>49</sup>,<sup>50</sup> J. S. Key,<sup>50</sup> N. Laal,<sup>51</sup> M. T. Lam,<sup>38</sup> W. G. Lamb,<sup>14</sup> T. J. W. Lazio,<sup>18</sup> K. Liu,<sup>12</sup> T. Liu,<sup>16</sup>
J. Luo,<sup>52</sup> R. S. Lynch,<sup>53</sup> D. R. Madison,<sup>54</sup> R. Main<sup>12</sup>,<sup>55</sup> R. Manchester,<sup>45</sup> A. McEwen,<sup>19</sup> J. McKee,<sup>55,56</sup>
M. A. McLaughlin<sup>16</sup>,<sup>57</sup> C. Ng,<sup>58</sup> D. J. Nice,<sup>58</sup> S. Ocker<sup>29</sup>,<sup>59</sup> K. D. Olum<sup>59</sup>,<sup>60</sup> T.
T. Pennucci,<sup>61</sup> B. B. P. Perera<sup>62</sup>,<sup>63</sup> D. Perrodin,<sup>63</sup> N. Porayko,<sup>12</sup> A. Possenti,<sup>63</sup> H. Quelquejay-Leclere,<sup>1</sup>
S. M. Ransom<sup>64</sup>,<sup>65</sup> P. S. Ray<sup>65</sup>,<sup>66</sup> D. J. Reardon<sup>30</sup>,<sup>67</sup> C. J. Russell<sup>66</sup>,<sup>68</sup> A. Samajdar,<sup>4</sup> J. Sarkissian,<sup>67</sup>
L. Schult,<sup>14</sup> G. Shaifullah<sup>4,15,63</sup>,<sup>64</sup> R. M. Shannon<sup>30</sup>,<sup>69</sup> B. J. Shapiro-Albert,<sup>68</sup> X. Siemens,<sup>9</sup> J. J. Simon,<sup>69</sup>
M. Siwek<sup>70</sup>,<sup>71</sup> T. L. Smith,<sup>71</sup> L. Speri<sup>72</sup>,<sup>73</sup> R. Spiewak<sup>48</sup>,<sup>74</sup> I. H. Stairs,<sup>32</sup> B. Stappers<sup>48</sup>,<sup>75</sup>
D. R. Stinebring,<sup>73</sup> J. K. Swiggum,<sup>58</sup> C. Tiburzi,<sup>63</sup> J. Turner,<sup>16</sup> A. Vecchio,<sup>74</sup> J. P. W. Verbiest<sup>75,12</sup>,<sup>76</sup>
H. Wahl,<sup>16</sup> S. Q. Wang,<sup>76</sup> J. Wang,<sup>77</sup> J. Wang,<sup>75</sup> Z. Wu,<sup>75</sup> L. Zhang,<sup>78</sup> and S. Zhang<sup>79</sup> (IPTA Collaboration)

*Affiliations are listed at the end of the paper*

Accepted 2023 February 24. Received 2023 February 17; in original form 2022 November 16

## ABSTRACT

The International Pulsar Timing Array 2nd data release is the combination of data sets from worldwide collaborations. In this study, we search for continuous waves: gravitational wave signals produced by individual supermassive black hole binaries in the local universe. We consider binaries on circular orbits and neglect the evolution of orbital frequency over the observational span. We find no evidence for such signals and set sky averaged 95 per cent upper limits on their amplitude $h_{95}$. The most sensitive frequency is 10 nHz with $h_{95} = 9.1 \times 10^{-15}$. We achieved the best upper limit to date at low and high frequencies of the PTA band thanks to improved effective cadence of observations. In our analysis, we have taken into account the recently discovered common red noise process, which has an impact at low frequencies. We also find that the peculiar noise features present in some pulsars data must be taken into account to reduce the false alarm. We show that using custom noise models is essential in searching for continuous gravitational wave signals and setting the upper limit.

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

## 1 INTRODUCTION

The goal of the Pulsar Timing Array (PTA) collaborations is to detect gravitational wave (GW) signals in the nanohertz band,

where we expect to see a gravitational wave background (GWB) produced by the superposition of GW signals from the population of supermassive black hole binaries (SMBHBs) (Jaffe & Backer 2003; Sesana, Vecchio & Colacino 2008; Maiorano, Paolis & Nucita 2021). Some individual SMBHBs might be brighter than most and stand above the stochastic signal; those are individually resolvable sources (Sesana, Vecchio & Volonteri 2009; Rosado, Sesana & Gair

<footnote>* E-mail: falxa@apc.in2p3.fr</footnote>
2015). The binaries detectable in the PTA band are in the orbits with the period from a few months to a few years and emit almost monochromatic GWs continuously during decades; we refer to those signals as continuous GWs (CGWs; Corbin & Cornish 2010; Ellis, Siemens & Creighton 2012; Babak et al. 2016; Aggarwal et al. 2019).

The GWs affect propagation of the radio emission from millisec-
ond pulsars leaving an imprint in the time of arrival (TOA) of pulses
observed with the radio telescopes. CGWs impact TOAs from all
observed millisecond pulsars in a deterministic manner characterized
by parameters of the SMBHBs. In this work, we consider the data
combined by the International Pulsar Timing Array (IPTA). IPTA is
a consortium of NANOGRAV (NANOGrav Collaboration 2015),
European Pulsar Timing Array (EPTA; Desvignes et al. 2016),
Australian (PPTA; Manchester et al. 2013), and Indian Pulsar Timing
Array (InPTA; Tarafdar et al. 2022) collaborations. In particular, we
analyse the second data released by IPTA (IPTA DR2) described in
details in Perera et al. (2019).

Recently, PTA collaborations have reported on the discovery of the common red noise signal, that is the stochastic signal with the spectral shape common to all pulsars in the array. Its high statistical significance was demonstrated independently by three collaborations (Arzoumanian et al. 2020; Goncharov et al. 2021; Chen et al. 2021) and, with even higher statistical confidence, was assessed using the IPTA DR2 (Antoniadis et al. 2022). We do not yet know the nature of this process, and its interpretation as GW background is inconclusive: the data are not informative enough to resolve the Helllings-Downs spatial correlations (Hellings & Downs 1983), which should be present in the case of the GW signal.

In this work, we search for continuous GWs which could be
present in the data in addition to the stochastic GWB. Following
the steps of previous studies (Zhu et al. 2014; Babak et al. 2016;
Aggarwal et al. 2019; Becsy & Cornish 2020; Becsy, Cornish &
Digman 2022; Arzoumanian et al. 2023), we search for a single
GW signal produced by an SMBHB binary in a circular orbit.
In our study, we neglect the pulsar terms during the search and
setting an upper limit on GW amplitude. However, we do an in-
depth analysis of the (weak) candidate events identified as plausible
GW signals. In the follow-up analysis on the restricted parameter
space (frequency and sky position), we extend our model to include
(i) pulsar term, (ii) eccentricity, (iii) extend the model beyond the
assumption of a single source. For the first time, we have included in
the analysis the common red component as a part of the total noise
model.

The main results of the paper can be summarized as follows. We did not detect any CGW signal and set an upper limit on GW amplitude. We have found that the noise model plays a crucial role in interpreting PTA observations. The detailed analysis performed on the most promising candidate event revealed that it could be explained by a time-correlated high-frequency noise in one of the pulsars.

The paper is organized as follows. In the next Section, we will
briefly describe the IPTA DR2 data set and the data model used
in the analysis. Most of the material needed for this Section is
available in the literature, and we heavily refer to it, keeping only
parts which are necessary for further presentation. In Section 3, we
describe the methodology which we have followed to get our results
presented in Section 4. In Section 5, we give a detailed follow-up
study of a most promising GW candidate event and demonstrate
the importance of noise modelling at high frequencies. We con-
clude with Section 6. Throughout the paper, we adopt geometrical
units G = c = 1.

2 IPTA DATA RELEASE 2 AND THE DATA
MODEL

## 2.1 IPTA DR2 data set

The IPTA DR2 consists of 65 stable MSPs with the duration of
observations up to 30 yr (Perera et al. 2019; Antoniadis et al.
2022).<fnref n="1" /> It combines the pulsar timing data acquired by IPTA
(Desvignes et al. 2016), North American Nanohertz Observatory
for Gravitational Waves (NANOGrav; Arzoumanian et al. 2016),
and the Parkes Pulsar Timing Array (Manchester et al. 2013). The
combined data is superior to the data sets of each collaboration: (i) it
has better sky coverage providing better localization of GW signals,
(ii) allows better decoupling and identification of noise components
due to increased number of observing backends, and (iii) reduces the
number of gaps in the data due to absence of observations. We have
already observed the improvement in the detection of the common
red noise process in Antoniadis et al. (2022) by using IPTA data.

The combined data set was analyzed to extract the properties
of individual pulsars (pulse frequency, spin-down, parallax, etc.)
by fitting a timing model (TM) that predicts the TOAs (Edwards,
Hobbs & Manchester 2006). Differences in predicted TOAs and
measured TOAs in the data set form the timing residuals. The
residuals are the result of various noise processes as well as the
interaction of the radio emission with GWs, which is the main subject
of this work.

## 2.2 Noise model

The noise of each pulsar data is modelled using the Gaussian process and split into several components [see Antoniadis et al. (2022) for details]. The white noise (WN), that estimates the TOA measurement errors, quantifies the radiometer noise in the receiver backend system and models the jitter noise which is intrinsic to the pulsar (statistics of pulse-to-pulse variation). The TM corrects deterministic TOA perturbations of physical origin (Edwards et al. 2006). Even though we fit the TOAs for the TM before we start the analysis, the fit might not be perfect and leave behind some residuals which we assume to be small and use a linear model<fnref n="2" /> to describe TM-generated errors. The low-frequency part of the data is strongly affected by red noise, which is a time-correlated process which power spectral density (P(f)) we describe as a simple power law. We distinguish achromatic red noise (RN, P<sub>RN</sub><sup>a</sup>) intrinsic to each pulsar (denoted by a subscript a) due to stochastic variations in the rotation of a pulsar and chromatic (i.e. dependent on the radio frequency at which pulses are observed) dispersion measurements variations (DM, P<sub>DM</sub><sup>a</sup>) noise caused by time-varying interstellar medium properties along the line of sight. Those two processes are described as

$$
P_{RN}(f) = \frac{A_{RN}^2}{12\pi^2} f_{yr}^{-3} \left(\frac{f}{f_{yr}}\right)^{-\gamma_{RN}}, \quad (1)
$$

$$
P_{\mathrm{DM}}(f) = \frac{A_{\mathrm{DM}}^2}{12\pi^2} f_{\mathrm{yr}}^{-3} \left(\frac{f}{f_{\mathrm{yr}}}\right)^{-\gamma_{\mathrm{DM}}} \left(\frac{1400 \, \mathrm{MHz}}{v}\right)^2, \quad (2)
$$

where v is radio observation frequency and $f_{yr} = \text{yr}^{-1}$ [see Anto-
niadis et al. (2022) for details]. Note that {$A_{RN/DM}$, $\gamma_{RN/DM}$} are
individual for each pulsar and we omit the pulsar index $\alpha$ to ease the
notations. In addition, for the first time, we will also add the common
red noise which is firmly established recently (Arzoumanian et al.

<footnote n="1">¹We used only the 53 pulsars of Antoniadis et al. (2022) for our analyses.</footnote>

<footnote n="2">²Linear in deviations from already determined parameters</footnote>
2015). The binaries detectable in the PTA band are in the orbits with the period from a few months to a few years and emit almost monochromatic GWs continuously during decades; we refer to those signals as continuous GWs (CGWs; Corbin & Cornish 2010; Ellis, Siemens & Creighton 2012; Babak et al. 2016; Aggarwal et al. 2019).

The GWs affect propagation of the radio emission from millisec-
ond pulsars leaving an imprint in the time of arrival (TOA) of pulses
observed with the radio telescopes. CGWs impact TOAs from all
observed millisecond pulsars in a deterministic manner characterized
by parameters of the SMBHBs. In this work, we consider the data
combined by the International Pulsar Timing Array (IPTA). IPTA is
a consortium of NANOGRAV (NANOGrav Collaboration 2015),
European Pulsar Timing Array (EPTA; Desvignes et al. 2016),
Australian (PPTA; Manchester et al. 2013), and Indian Pulsar Timing
Array (InPTA; Tarafdar et al. 2022) collaborations. In particular, we
analyse the second data released by IPTA (IPTA DR2) described in
details in Perera et al. (2019).

Recently, PTA collaborations have reported on the discovery of the common red noise signal, that is the stochastic signal with the spectral shape common to all pulsars in the array. Its high statistical significance was demonstrated independently by three collaborations (Arzoumanian et al. 2020; Goncharov et al. 2021; Chen et al. 2021) and, with even higher statistical confidence, was assessed using the IPTA DR2 (Antoniadis et al. 2022). We do not yet know the nature of this process, and its interpretation as GW background is inconclusive: the data are not informative enough to resolve the Helllings-Downs spatial correlations (Hellings & Downs 1983), which should be present in the case of the GW signal.

In this work, we search for continuous GWs which could be
present in the data in addition to the stochastic GWB. Following
the steps of previous studies (Zhu et al. 2014; Babak et al. 2016;
Aggarwal et al. 2019; Becsy & Cornish 2020; Becsy, Cornish &
Digman 2022; Arzoumanian et al. 2023), we search for a single
GW signal produced by an SMBHB binary in a circular orbit.
In our study, we neglect the pulsar terms during the search and
setting an upper limit on GW amplitude. However, we do an in-
depth analysis of the (weak) candidate events identified as plausible
GW signals. In the follow-up analysis on the restricted parameter
space (frequency and sky position), we extend our model to include
(i) pulsar term, (ii) eccentricity, (iii) extend the model beyond the
assumption of a single source. For the first time, we have included in
the analysis the common red component as a part of the total noise
model.

The main results of the paper can be summarized as follows. We did not detect any CGW signal and set an upper limit on GW amplitude. We have found that the noise model plays a crucial role in interpreting PTA observations. The detailed analysis performed on the most promising candidate event revealed that it could be explained by a time-correlated high-frequency noise in one of the pulsars.

The paper is organized as follows. In the next Section, we will
briefly describe the IPTA DR2 data set and the data model used
in the analysis. Most of the material needed for this Section is
available in the literature, and we heavily refer to it, keeping only
parts which are necessary for further presentation. In Section 3, we
describe the methodology which we have followed to get our results
presented in Section 4. In Section 5, we give a detailed follow-up
study of a most promising GW candidate event and demonstrate
the importance of noise modelling at high frequencies. We con-
clude with Section 6. Throughout the paper, we adopt geometrical
units G = c = 1.

2 IPTA DATA RELEASE 2 AND THE DATA
MODEL

2.1 IPTA DR2 data set

The IPTA DR2 consists of 65 stable MSPs with the duration of
observations up to 30 yr (Perera et al. 2019; Antoniadis et al.
2022).<fnref n="1" /> It combines the pulsar timing data acquired by IPTA
(Desvignes et al. 2016), North American Nanohertz Observatory
for Gravitational Waves (NANOGrav; Arzoumanian et al. 2016),
and the Parkes Pulsar Timing Array (Manchester et al. 2013). The
combined data is superior to the data sets of each collaboration: (i) it
has better sky coverage providing better localization of GW signals,
(ii) allows better decoupling and identification of noise components
due to increased number of observing backends, and (iii) reduces the
number of gaps in the data due to absence of observations. We have
already observed the improvement in the detection of the common
red noise process in Antoniadis et al. (2022) by using IPTA data.

The combined data set was analyzed to extract the properties
of individual pulsars (pulse frequency, spin-down, parallax, etc.)
by fitting a timing model (TM) that predicts the TOAs (Edwards,
Hobbs & Manchester 2006). Differences in predicted TOAs and
measured TOAs in the data set form the timing residuals. The
residuals are the result of various noise processes as well as the
interaction of the radio emission with GWs, which is the main subject
of this work.

## 2.2 Noise model

The noise of each pulsar data is modelled using the Gaussian process and split into several components [see Antoniadis et al. (2022) for details]. The white noise (WN), that estimates the TOA measurement errors, quantifies the radiometer noise in the receiver backend system and models the jitter noise which is intrinsic to the pulsar (statistics of pulse-to-pulse variation). The TM corrects deterministic TOA perturbations of physical origin (Edwards et al. 2006). Even though we fit the TOAs for the TM before we start the analysis, the fit might not be perfect and leave behind some residuals which we assume to be small and use a linear model<fnref n="2" /> to describe TM-generated errors. The low-frequency part of the data is strongly affected by red noise, which is a time-correlated process which power spectral density (P(f)) we describe as a simple power law. We distinguish achromatic red noise (RN, P<sub>RN</sub><sup>a</sup>) intrinsic to each pulsar (denoted by a subscript a) due to stochastic variations in the rotation of a pulsar and chromatic (i.e. dependent on the radio frequency at which pulses are observed) dispersion measurements variations (DM, P<sub>DM</sub><sup>a</sup>) noise caused by time-varying interstellar medium properties along the line of sight. Those two processes are described as

$$
P_{RN}(f) = \frac{A_{RN}^2}{12\pi^2} f_{yr}^{-3} \left(\frac{f}{f_{yr}}\right)^{-\gamma_{RN}}, \qquad (1)
$$

$$
P_{DM}(f) = \frac{A_{DM}^2}{12\pi^2} f_{yr}^{-3} \left(\frac{f}{f_{yr}}\right)^{-\gamma_{DM}} \left(\frac{1400 \text{ MHz}}{v}\right)^2, \qquad (2)
$$

where v is radio observation frequency and $f_{yr} = \text{yr}^{-1}$ [see Antoniadis et al. (2022) for details]. Note that {$A_{RN/DM}$, $\gamma_{RN/DM}$} are individual for each pulsar and we omit the pulsar index $\alpha$ to ease the notations. In addition, for the first time, we will also add the common red noise which is firmly established recently (Arzoumanian et al.

<footnote n="1">¹We used only the 53 pulsars of Antoniadis et al. (2022) for our analyses.</footnote>

<footnote n="2">²Linear in deviations from already determined parameters</footnote>
2020; Chen et al. 2021; Goncharov et al. 2021; Antoniadis et al. 2022), we denote it as {$A_{cm}, \gamma_{cm}$}. This is a red noise component with the spectral properties shared across all pulsars in the array. The nature of this signal is still unclear, there is not enough evidence to support its GW origin, and it is a subject of current active investigations; for now, we call it 'crn'.

All observations are translated to the Solar system barycenter (SSB) frame. The transformation from the Earth's based frame to SSB relies on a precise knowledge of the Solar system ephemeris (SSE): in this work, we use DE438 ephemeris (Folkner & Park 2018). It was noted that there could be unaccounted systematic errors in the SSE, which could be mistaken for a stochastic GW signal. We have included BAYESEPHEM (Vallisneri et al. 2020) in the data model to mitigate those potential errors. BAYESEPHEM is a phenomenological model that varies the orbital elements of major external planets and takes into account possible systematics in SSE; note that it might also absorb part of the stochastic GW signal.

The GW background would require including the Hellings-Downs correlations between pulsars in the data model. We do neglect cross-correlation terms in our analysis, reducing in practice the GW background to the detected CRN. This is justified because: (i) any cross-correlation present in the data is rather weak; otherwise, it would have been detected (ii) the autocorrelation part is captured already in the CRN process that we include in the model.

All in all, the model of the timing residuals is a superposition of the noise components described above. In addition, we will assume (and test this hypothesis) that the data contain a deterministic continuous GW signal $s_{CGW}(t)$:

$$ \vec{\delta}t = \mathbf{M}\vec{\epsilon} + n_{\mathrm{WN}} + n_{\mathrm{RN}} + n_{\mathrm{CRN}} + n_{\mathrm{DM}} + s_{\mathrm{SSE}}(t) + s_{\mathrm{CGW}}(t). \quad (3) $$

The TM contribution is described by $\mathbf{M}\vec{\epsilon}$ where $\mathbf{M}$ is a design matrix and $\vec{\epsilon}$ are the linear corrections to the TM parameters; $n_{WN}$, $n_{RN}$, $n_{CRN}$, $n_{DM}$ are the components described above and correspond to the WN, the individual red noise, the common red noise, and the dispersion measurement variations noise; the signal $s_{SSE}(t)$ denotes the BAYESEPHEM model for SSE systematics.

We base our analysis on the noise model derived for each pulsar independently. This approximation assumes that the GW signal contribution is sub-dominant and neglected in modelling each pulsar data. In fact, the GW signal will be absorbed into the RN during this step and should be decoupled when we analyse the full array allowing RN parameters to vary. For the main analysis, we fix parameters of the WN component to their maximum-likelihood values obtained from the single pulsar noise analysis. It was shown (Chalumeau et al. 2021; Antoniadis et al. 2022) that this does not affect the result for the stochastic GW signal search, and we assume the same for continuous GW search. This assumption tremendously reduces the parameter space, which otherwise would be computationally intractable.

We model each noise component as a Gaussian process (van Haasteren & Vallisneri 2014) using sin and cos as basis functions evaluated at $f_k = k/T_{obs}$ Fourier frequencies, where $T_{obs}$ is time-span of observations:

$$ n(t) = \sum_{k}^{N} X_k \cos(2\pi t f_k) + Y_k \sin(2\pi t f_k), \quad (4) $$

where the $X_k$ and $Y_k$ are the Gaussian-distributed weights with the covariance matrix defined by the power spectral density of the noise. In our approach, we marginalize over the weights. In previous PTA analysis (Zhu et al. 2014; Babak et al. 2016; Aggarwal et al. 2019), the number of Fourier components ($N$) was fixed to 30 for both RN and DM.

This choice was motivated by computational savings and that those noise components mainly contribute at the low-frequency end of PTA sensitivity. However, the recent study based on the Bayesian model selection has shown that the noise models with specific values of $N$ for each noise component and each pulsar are better supported by the data (Chalumeau et al. 2021).

## 2.3 Continuous gravitational waves

The concept of detecting GWs with PTA was formulated in Sazhin (1978); Detweiler (1979). The response to a deterministic GW signal can be written as

$$ s_a(t, \omega) = \int_0^t \frac{1}{2} \frac{\hat{p}_a^j \hat{p}_a^j}{1 + \hat{p}_a \cdot \hat{\Omega}} \Delta h_{ij}(t') dt', \quad (5) $$

where $\hat{p}_a$ is the unit vector pointing to the pulsar $a$ in the sky, $\hat{\Omega}$ is the direction of GW propagation and $h_{ij}$ is the GW strain in the transverse-traceless gauge ($i$ and $j$ are the spatial indices). The response depends on the GW strain at two instances of time: the time of emission of electromagnetic signal and the time of its reception:

$$ \Delta h_{ij}(t) = h_{ij}(t - \tau_a) - h_{ij}(t), \quad (6) $$

where

$$ \tau_a = L_a(1 + \hat{\Omega} \cdot \hat{p}_a) \quad (7) $$

and $L_a$ is the distance to the pulsar $a$. The time difference in the strain corresponds to the light travel time between the Earth and the pulsar with a geometrical factor. The corresponding two terms in the expression of the timing residuals are usually referred to as Earth $s_e$ and pulsar $s_p$ terms:

$$ s_a(t) = s_{p,a}(t - \tau_a) - s_{e,a}(t). \quad (8) $$

The strain amplitude of a GW produced by a circular binary system is given by

$$ h_{ij}(t, \Omega) = \sum_{A=+,-}^{A^+} e_{ij}^A(\hat{\Omega}) h_A(t), \quad (9) $$

where $e_{ij}^{+,-}$ are two GW polarization tensor defined as

$$ e_{ij}^{+}(\hat{\Omega}) = \hat{m}_i \hat{m}_j - \hat{n}_i \hat{n}_j, \quad (10) $$

$$ e_{ij}^{+}(\hat{\Omega}) = \hat{m}_i \hat{m}_j - \hat{n}_i \hat{n}_j, \quad (11) $$

$$ \hat{\Omega} = -\sin\theta \cos\phi\hat{x} - \sin\theta \sin\phi\hat{y} - \cos\theta\hat{z} \quad (12) $$

and the unit vectors are

$$ \hat{m} = -\sin\phi\hat{x} + \cos\phi\hat{y}, \quad (13) $$

$$ \hat{n} = -\cos\theta \cos\phi\hat{x} - \cos\theta \sin\phi\hat{y} + \sin\theta\hat{z}, \quad (14) $$

where $(\theta, \phi)$ are the polar coordinates of the GW source sky location.

Plugging these expressions in equation (5), we obtain the timing residuals expected in the PTA data for a CGW signal coming from a circular SMBHB:

$$ s_a(t, \hat{\Omega}) = \sum_{A} F^A(\hat{\Omega})[s_A(t) - s_A(t - \tau_a)] \quad (15) $$

with:

<footnote n="3">³Note, NG have used a different model for DM, namely DMX, which is not decomposed in Fourier basis functions (NANOGrav Collaboration 2015; Lam et al. 2016).</footnote>
2020; Chen et al. 2021; Goncharov et al. 2021; Antoniadis et al. 2022), we denote it as {$A_{cm}, \gamma_{cm}$}. This is a red noise component with the spectral properties shared across all pulsars in the array. The nature of this signal is still unclear, there is not enough evidence to support its GW origin, and it is a subject of current active investigations; for now, we call it 'crn'.

All observations are translated to the Solar system barycenter (SSB) frame. The transformation from the Earth's based frame to SSB relies on a precise knowledge of the Solar system ephemeris (SSE): in this work, we use DE438 ephemeris (Folkner & Park 2018). It was noted that there could be unaccounted systematic errors in the SSE, which could be mistaken for a stochastic GW signal. We have included BAYESEPHEM (Vallisneri et al. 2020) in the data model to mitigate those potential errors. BAYESEPHEM is a phenomenological model that varies the orbital elements of major external planets and takes into account possible systematics in SSE; note that it might also absorb part of the stochastic GW signal.

The GW background would require including the Hellings-Downs correlations between pulsars in the data model. We do neglect cross-correlation terms in our analysis, reducing in practice the GW background to the detected CRN. This is justified because: (i) any cross-correlation present in the data is rather weak; otherwise, it would have been detected (ii) the autocorrelation part is captured already in the CRN process that we include in the model.

All in all, the model of the timing residuals is a superposition of the noise components described above. In addition, we will assume (and test this hypothesis) that the data contain a deterministic continuous GW signal $s_{CGW}(t)$:

$$ \vec{\delta}t = \mathbf{M}\vec{\epsilon} + n_{\text{WN}} + n_{\text{RN}} + n_{\text{CRN}} + n_{\text{DM}} + s_{\text{SSE}}(t) + s_{\text{CGW}}(t). \quad (3) $$

The TM contribution is described by $\mathbf{M}\vec{\epsilon}$ where $\mathbf{M}$ is a design matrix and $\vec{\epsilon}$ are the linear corrections to the TM parameters; $n_{WN}$, $n_{RN}$, $n_{CRN}$, $n_{DM}$ are the components described above and correspond to the WN, the individual red noise, the common red noise, and the dispersion measurement variations noise; the signal $s_{SSE}(t)$ denotes the BAYESEPHEM model for SSE systematics.

We base our analysis on the noise model derived for each pulsar independently. This approximation assumes that the GW signal contribution is sub-dominant and neglected in modelling each pulsar data. In fact, the GW signal will be absorbed into the RN during this step and should be decoupled when we analyse the full array allowing RN parameters to vary. For the main analysis, we fix parameters of the WN component to their maximum-likelihood values obtained from the single pulsar noise analysis. It was shown (Chalumeau et al. 2021; Antoniadis et al. 2022) that this does not affect the result for the stochastic GW signal search, and we assume the same for continuous GW search. This assumption tremendously reduces the parameter space, which otherwise would be computationally intractable.

We model each noise component as a Gaussian process (van Haasteren & Vallisneri 2014) using sin and cos as basis functions evaluated at $f_k = k/T_{\text{obs}}$ Fourier frequencies, where $T_{\text{obs}}$ is time-span of observations:

$$ n(t) = \sum_{k}^{N} X_k \cos(2\pi t f_k) + Y_k \sin(2\pi t f_k), \quad (4) $$

where the $X_k$ and $Y_k$ are the Gaussian-distributed weights with the covariance matrix defined by the power spectral density of the noise. In our approach, we marginalize over the weights. In previous PTA analysis (Zhu et al. 2014; Babak et al. 2016; Aggarwal et al. 2019), the number of Fourier components ($N$) was fixed to 30 for both RN and DM.<fnref n="3" />

This choice was motivated by computational savings and that those noise components mainly contribute at the low-frequency end of PTA sensitivity. However, the recent study based on the Bayesian model selection has shown that the noise models with specific values of $N$ for each noise component and each pulsar are better supported by the data (Chalumeau et al. 2021).

## 2.3 Continuous gravitational waves

The concept of detecting GWs with PTA was formulated in Sazhin (1978); Detweiler (1979). The response to a deterministic GW signal can be written as

$$ s_a(t, \omega) = \int_0^t \frac{1}{2} \frac{\hat{p}_a^j \hat{p}_a^j}{1 + \hat{p}_a \cdot \hat{\Omega}} \Delta h_{ij}(t') dt', \quad (5) $$

where $\hat{p}_a$ is the unit vector pointing to the pulsar $a$ in the sky, $\hat{\Omega}$ is the direction of GW propagation and $h_{ij}$ is the GW strain in the transverse-traceless gauge ($i$ and $j$ are the spatial indices). The response depends on the GW strain at two instances of time: the time of emission of electromagnetic signal and the time of its reception:

$$ \Delta h_{ij}(t) = h_{ij}(t - \tau_a) - h_{ij}(t), \quad (6) $$

where

$$ \tau_a = L_a(1 + \hat{\Omega} \cdot \hat{p}_a) \quad (7) $$

and $L_a$ is the distance to the pulsar $a$. The time difference in the strain corresponds to the light travel time between the Earth and the pulsar with a geometrical factor. The corresponding two terms in the expression of the timing residuals are usually referred to as Earth $s_e$ and pulsar $s_p$ terms:

$$ s_a(t) = s_{p,a}(t - \tau_a) - s_{e,a}(t). \quad (8) $$

The strain amplitude of a GW produced by a circular binary system is given by

$$ h_{ij}(t, \Omega) = \sum_{A=+,-} e_{ij}^{A}(\hat{\Omega}) h_{A}(t), \quad (9) $$

where $e_{ij}^{+,-}$ are two GW polarization tensor defined as

$$ e_{ij}^{+}(\hat{\Omega}) = \hat{m}_i \hat{m}_j - \hat{n}_i \hat{n}_j, \quad (10) $$

$$ e_{ij}^{+}(\hat{\Omega}) = \hat{m}_i \hat{m}_j - \hat{n}_i \hat{n}_j, \quad (11) $$

$$ \hat{\Omega} = -\sin\theta \cos\phi\hat{x} - \sin\theta \sin\phi\hat{y} - \cos\theta\hat{z} \quad (12) $$

and the unit vectors are

$$ \hat{m} = -\sin\phi\hat{x} + \cos\phi\hat{y}, \quad (13) $$

$$ \hat{n} = -\cos\theta \cos\phi\hat{x} - \cos\theta \sin\phi\hat{y} + \sin\theta\hat{z}, \quad (14) $$

where $(\theta, \phi)$ are the polar coordinates of the GW source sky location.

Plugging these expressions in equation (5), we obtain the timing residuals expected in the PTA data for a CGW signal coming from a circular SMBHB:

$$ s_a(t, \hat{\Omega}) = \sum_{A} F^A(\hat{\Omega})[s_A(t) - s_A(t - \tau_a)] \quad (15) $$

with:

<footnote n="3">³Note, NG have used a different model for DM, namely DMX, which is not decomposed in Fourier basis functions (NANOGrav Collaboration 2015; Lam et al. 2016).</footnote>
$$s_+(t) = \frac{M^{5/3}}{d_L \omega(t)^{1/3}} [-\sin[2\Phi(t)](1 + \cos^2 t) \cos 2\psi - 2\cos[2\Phi(t)] \cos t \sin 2\psi], \quad (16)$$

$$s_{\times}(t) = \frac{M^{5/3}}{d_L \omega(t)^{1/3}} [-\sin[2\Phi(t)](1 + \cos^2 t) \cos 2\psi + 2\cos[2\Phi(t)] \cos t \sin 2\psi], \quad (17)$$

where $M$ is the chirp mass, $d_L$ the luminosity distance, $\omega(t)$ the CGW orbital angular frequency, $t$ is the orbital inclination to the line of sight, $\psi$ is a polarization angle, and $\Phi(t)$ is the phase of CGW.

The $F^A$ are the antenna pattern functions (Sesana & Vecchio 2010; Babak & Sesana 2012; Ellis et al. 2012; Taylor et al. 2016) given as

$$F^+(\hat{\Omega}) = \frac{1}{2} \frac{(\hat{m} \cdot \hat{p})^2 - (\hat{n} \cdot \hat{p})^2}{1 + \hat{\Omega} \cdot \hat{p}}, \quad (18)$$

$$F^{\times}(\hat{\Omega}) = \frac{(\hat{m} \cdot \hat{p})(\hat{n} \cdot \hat{p})}{1 + \hat{\Omega} \cdot \hat{p}}. \quad (19)$$

In this work, we neglect the pulsar term considering it as a part of the noise, assuming that the source has evolved sufficiently over $\tau_\alpha$ to move the pulsar term off the earth-term frequency. Including pulsar term should improve the parameter estimation but comes with a huge price of the increase in the complexity of the likelihood surface and the dimensionality of parameter space [2 additional parameters per pulsar for phase and frequency of pulsar term, e.g. see Corbin & Cornish (2010)]. We foresee the possibility of following up the candidate events (identified using the earth term only) with the extended signal model (pulsar term, eccentric orbit) on the reduced parameter space. We also neglect the evolution of the GW frequency ($\omega_o = 2\pi f_{gw}$) over the observation time. The frequency evolution becomes potentially measurable for the heavy sources emitting at frequency $\ge 10^{-7}$ Hz, neglecting the frequency evolution does not prevent us from detecting the sources but introduces a bias in the measured GW frequency (overestimating it), for more details see conclusion in Petiteau et al. (2013). So the CGW phase takes a very simple form:

$$\Phi(t) = \omega_0 t + \phi_0/2, \quad (20)$$

where $\phi_0$ is initial orbital phase. Finally, the CGW amplitude $h$ is a function of $M$, $d_L$ and $f_{gw}$ given by

$$h = \frac{2M^{5/3}(\pi f_{gw})^{-2/3}}{d_L}. \quad (21)$$

We consider the model containing only one CGW signal. This model still detects multiple CGW if they are present in the data at the non-overlapping Fourier frequencies (see Babak & Sesana 2012 for discussion). If we find more than one candidate with sufficient statistical significance as potential GW sources, we will conduct additional investigations extending our model to include several CGWs. We start the analysis with 1 CGW source characterized by 7 parameters summarized in Table 1 together with their prior range (we always assume a uniform prior)$^4$.

## 3 METHOD

We work within the Bayesian framework and start with running the search for the CGW signal. As mentioned above, we sample

<table><caption>**Table 1.** List of the CGW parameters as defined in our model with their respective ranges.</caption><thead><tr><th>CGW parameter</th><th>Range</th></tr></thead><tbody><tr><td>log<sub>10</sub>h</td><td>[-18, -11]</td></tr><tr><td>f<sub>gw</sub> (Hz)</td><td>[10<sup>-9</sup>, 10<sup>-7</sup>]</td></tr><tr><td>φ<sub>0</sub></td><td>[0, 2π]</td></tr><tr><td>cos t</td><td>[-1, 1]</td></tr><tr><td>ψ</td><td>[0, π]</td></tr><tr><td>θ</td><td>[0, π]</td></tr><tr><td>φ</td><td>[0, 2π]</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** 2D sky sensitivity of IPTA for CGW signal around most sensitive frequency. The figure was obtained computing the 95 per cent upper limit for $h$ on 12 patches across the sky. Gaussian interpolation was used to smoothen it across boundaries.</figcaption></figure>

parameters of CGW together with the noise parameters for RN and DM. We keep the white WN parameters fixed and marginalize over the TM and BAYESEPHEM parameters. We made two runs: with and without CRN, to check how much it affects our result. We always use Markov chain Monte Carlo (MCMC) sampler (Ellis & van Haasteren 2017), and we use ENTERPRISE (Ellis et al. 2020) software to construct the models and compute the likelihood and prior probability.

We extensively use the single pulsar noise explorations runs performed before the main analysis. We have converted posteriors for the RN and DM into 2D histograms ({A<sub>RN, DM</sub>, γ<sub>RN, DM</sub>} and use them as one of a proposal for those parameters. This empirical proposal improves the efficiency of MCMC and reduces the autocorrelation length of the chain (Aggarwal et al. 2019).

During the search, we compute the Bayes factor (BF), comparing the null model (noise only) against the model where we have a CGW signal on the logarithmically spaced frequencies. In the absence of the detection, we proceed to setting an upper limit, building a 95 per cent sensitivity curve. During the upper limit analysis, we used a uniform prior on the amplitude of the GW signal.

IPTA data contains 53-millisecond pulsars; however, not all of them are equally sensitive to the CGW. We have selected 21 pulsars which, on average, recover 95 per cent of the array's total signal-to-noise ratio (SNR) to CGW. The selected pulsars are depicted on the projected sky map in Fig. 1 as red stars; we have used large green stars and annotation for the four best timers. The ranking method is briefly outlined in Babak et al. (2016) and in greater detail in Speri et al. (2022). This significantly reduces the computational cost without much affecting the final results.$^5$

<footnote>$^4$For setting an upper limit we use uniform prior on the amplitude of GW strain</footnote>

<footnote>$^5$We have analyzed 21 pulsar, the noise model for each pulsar is characterized by 4 parameters (an amplitude and a spectral index for RN and DM); Bayesphem adds 11 parameters. Finally, we have 7 parameters for the continuous wave and 2 for the common red noise, giving a total of 104 parameters to sample.</footnote>
$$s_{+}(t) = \frac{M^{5/3}}{d_{L}\omega(t)^{1/3}} [-\sin[2\Phi(t)](1 + \cos^{2}\iota)\cos 2\psi - 2\cos[2\Phi(t)]\cos \iota \sin 2\psi], \quad (16)$$

$$s_{\times}(t) = \frac{M^{5/3}}{d_{L}\omega(t)^{1/3}} [-\sin[2\Phi(t)](1 + \cos^{2}\iota)\cos 2\psi + 2\cos[2\Phi(t)]\cos \iota \sin 2\psi], \quad (17)$$

where $M$ is the chirp mass, $d_L$ the luminosity distance, $\omega(t)$ the CGW orbital angular frequency, $\iota$ is the orbital inclination to the line of sight, $\psi$ is a polarization angle, and $\Phi(t)$ is the phase of CGW.

The $F^A$ are the antenna pattern functions (Sesana & Vecchio 2010; Babak & Sesana 2012; Ellis et al. 2012; Taylor et al. 2016) given as

$$F^{+}(\hat{\Omega}) = \frac{1}{2} \frac{(\hat{m} \cdot \hat{p})^2 - (\hat{m} \cdot \hat{p})^2}{1 + \hat{\Omega} \cdot \hat{p}}, \quad (18)$$

$$F^{\times}(\hat{\Omega}) = \frac{(\hat{m} \cdot \hat{p})(\hat{n} \cdot \hat{p})}{1 + \hat{\Omega} \cdot \hat{p}}. \quad (19)$$

In this work, we neglect the pulsar term considering it as a part of the noise, assuming that the source has evolved sufficiently over $\tau_\alpha$ to move the pulsar term off the earth-term frequency. Including pulsar term should improve the parameter estimation but comes with a huge price of the increase in the complexity of the likelihood surface and the dimensionality of parameter space [2 additional parameters per pulsar for phase and frequency of pulsar term, e.g. see Corbin & Cornish (2010)]. We foresee the possibility of following up the candidate events (identified using the earth term only) with the extended signal model (pulsar term, eccentric orbit) on the reduced parameter space. We also neglect the evolution of the GW frequency ($\omega_o = 2\pi f_{gw}$) over the observation time. The frequency evolution becomes potentially measurable for the heavy sources emitting at frequency $\ge 10^{-7}$ Hz, neglecting the frequency evolution does not prevent us from detecting the sources but introduces a bias in the measured GW frequency (overestimating it), for more details see conclusion in Petiteau et al. (2013). So the CGW phase takes a very simple form:

$$\Phi(t) = \omega_0 t + \phi_0/2, \quad (20)$$

where $\phi_0$ is initial orbital phase. Finally, the CGW amplitude $h$ is a function of $M$, $d_L$ and $f_{gw}$ given by

$$h = \frac{2M^{5/3}(\pi f_{gw})^{-2/3}}{d_L}. \quad (21)$$

We consider the model containing only one CGW signal. This model still detects multiple CGW if they are present in the data at the non-overlapping Fourier frequencies (see Babak & Sesana 2012 for discussion). If we find more than one candidate with sufficient statistical significance as potential GW sources, we will conduct additional investigations extending our model to include several CGWs. We start the analysis with 1 CGW source characterized by 7 parameters summarized in Table 1 together with their prior range (we always assume a uniform prior)$^4$.

## 3 METHOD

We work within the Bayesian framework and start with running the search for the CGW signal. As mentioned above, we sample

<table><caption>**Table 1.** List of the CGW parameters as defined in our model with their respective ranges.</caption><thead><tr><th>CGW parameter</th><th>Range</th></tr></thead><tbody><tr><td>log<sub>10</sub>h</td><td>[-18, -11]</td></tr><tr><td>f<sub>gw</sub> (Hz)</td><td>[10<sup>-9</sup>, 10<sup>-7</sup>]</td></tr><tr><td>φ<sub>0</sub></td><td>[0, 2π]</td></tr><tr><td>cos ι</td><td>[-1, 1]</td></tr><tr><td>ψ</td><td>[0, π]</td></tr><tr><td>θ</td><td>[0, π]</td></tr><tr><td>φ</td><td>[0, 2π]</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** 2D sky sensitivity of IPTA for CGW signal around most sensitive frequency. The figure was obtained computing the 95 per cent upper limit for $h$ on 12 patches across the sky. Gaussian interpolation was used to smoothen it across boundaries.</figcaption></figure>

parameters of CGW together with the noise parameters for RN and DM. We keep the white WN parameters fixed and marginalize over the TM and BAYESEPHEM parameters. We made two runs: with and without CRN, to check how much it affects our result. We always use Markov chain Monte Carlo (MCMC) sampler (Ellis & van Haasteren 2017), and we use ENTERPRISE (Ellis et al. 2020) software to construct the models and compute the likelihood and prior probability.

We extensively use the single pulsar noise explorations runs performed before the main analysis. We have converted posteriors for the RN and DM into 2D histograms ($\{A_{RN, DM}, \gamma_{RN, DM}\}$) and use them as one of a proposal for those parameters. This empirical proposal improves the efficiency of MCMC and reduces the autocorrelation length of the chain (Aggarwal et al. 2019).

During the search, we compute the Bayes factor (BF), comparing the null model (noise only) against the model where we have a CGW signal on the logarithmically spaced frequencies. In the absence of the detection, we proceed to setting an upper limit, building a 95 per cent sensitivity curve. During the upper limit analysis, we used a uniform prior on the amplitude of the GW signal.

IPTA data contains 53-millisecond pulsars; however, not all of them are equally sensitive to the CGW. We have selected 21 pulsars which, on average, recover 95 per cent of the array's total signal-to-noise ratio (SNR) to CGW. The selected pulsars are depicted on the projected sky map in Fig. 1 as red stars; we have used large green stars and annotation for the four best timers. The ranking method is briefly outlined in Babak et al. (2016) and in greater detail in Speri et al. (2022). This significantly reduces the computational cost without much affecting the final results.$^5$

<footnote>$^4$For setting an upper limit we use uniform prior on the amplitude of GW strain</footnote>

<footnote>$^5$We have analyzed 21 pulsar, the noise model for each pulsar is characterized by 4 parameters (an amplitude and a spectral index for RN and DM); Bayesphem adds 11 parameters. Finally, we have 7 parameters for the continuous wave and 2 for the common red noise, giving a total of 104 parameters to sample.</footnote>
Historically, we have performed the search with the noise model with the uniform settings across all pulsars in the array; namely, we have used 30 frequency bins for modelling both RN and DM processes [as it was done in preceding work Arzoumanian et al. (2014); Aggarwal et al. (2019)]. Using this uniform setting we have obtained quite erroneous results, and after a long investigation, we realized that this noise model does not adequately describe the observational data (see Section 5). We have switched to another noise model where we have used the custom-made noise model for the six best EPTA pulsars (see Chalumeau et al. 2021 for details), and for other pulsars, we changed the number of used frequency bins: RN30DM100 – 30 bins for the RN and 100 bins for DM modelling. We will continue using this short-hand notation for the noise model, showing explicitly the number of Fourier frequencies used by the base functions in the Gaussian process describing the corresponding noise. The results presented in the next section were obtained using the ‘custom’ made noise model; we postpone the detailed discussion on the noise model selection and influence until Section 5.

## 3.1 Model selection

During the search, we consider two competing models: noise only and noise with 1 CGW signal. We compute the BF to measure which model is preferred by the IPTA DR2 data set. In particular we employ hyper-model jumps to compute BF following the methods outlined in Hee et al. (2015) and implemented in the enterprise extensions (Taylor et al. 2021). In this method, we introduce a hyperparameter which indexes the models; then, we perform sampling inside each model and in this parameter. The BF then is given by the ratio of the number of iterations that the chain spends in each model. For example, if we consider two models $M_A$ and $M_B$ with the hyperparameter $n$, then the BF is the ratio

$$B_B^A = \frac{n_A}{n_B}, \qquad (22)$$

where $n_A$, $n_B$ counts the number of samples in the chain corresponding to the models A, B. In our previous investigations, we have found that this method gives quite a reliable result, e.g. see Chalumeau et al. (2021) where comparison is done against the evidence computation obtained with Dynasty nested sampling (Skilling 2006; Speagle 2020).

Later in the paper we will compute BF between the models where we also vary the noise.

## 3.2 Upper limit

For obtaining the upper limit, we again use MCMC assuming a uniform prior distribution for log $h$ within the bounds [−18, −11] on the fixed set of $f_{gw}$. We use marginalized posterior probability distribution for the CGW amplitude $p(h)$ to set 95 per cent upper limit $h_{95}$ defined as

$$0.95 = \int_{0}^{h_{95}} p(h)dh. \qquad (23)$$

We use a grid of 100 logarithmically spaced frequencies between $10^{-9}$ and $10^{-7}$ Hz. The lowest bound determined by the IPTA observational span $1/T_{obs}$ that is $\approx 10^{-9}$ Hz while the upper bound, $10^{-7}$ Hz, is constrained by reduction in the sensitivity due to response ($\propto f^{-1/3}$) and by our assumption that GW frequency does not evolve. We should take into account the frequency evolution of the GW signal above $10^{-7}$ Hz during the analysis, as discussed in Petiteau et al. (2013). As mentioned above, we compute the upper limit for

<figure><img src="image_3.png" /><figcaption>Figure 2. On the top panel, the upper bound of the 95 per cent central credible region on log $h$ obtained with a log-uniform prior (with CRN). Bottom panel is the associated BF at given CGW frequency. The black dashed line shows where BF is equal to 1. The two peaks around 35 and 70 nHz with BF > 1 are the 1 and 2 yr and should not be taken into account.</figcaption></figure>

two models of noise; with and without CRN. The inclusion of CRN adds two more parameters to the model (an amplitude parameter $A_{CRN}$ and a spectral index $\gamma_{CRN}$), which we sample together with parameters of CGW and other noise components.

## 4 RESULTS

We compute BF between two models $B_{custom}^{custom+CGW}$ for a noise only model $M_{custom}$ and noise + CGW signal model $M_{custom+CGW}$. The subscript ‘custom’ corresponds to the noise model, we have used for these main results and distinguishes it from other noise models considered in the next section. Similarly to the upper limit run, we have used a log-uniform prior on $h$ bounded by [−18, −11], and the BF was computed on the grid of 100 CGW frequencies ($f_{gw}$) between 1 and 100 nHz.

The main result of this paper is summarized in Fig. 2. In the lower panel we plot the BF, and, as one can see, the noise model is usually preferred. There are few spikes where BF reaches 1, this is definitely not a detection; however, PTA efforts probably should monitor carefully those frequencies in the future extended IPTA data sets. The excess in BF was also used as identification of the CGW candidates to follow. These results show that there is nothing to follow:$^6 B_{custom}^{custom+CGW} \le 1$.

The top panel of Fig. 2 shows the upper bound of the 95 per cent central credible region for CGW strain computed at the same set of frequencies. Note the ‘spiky’ features at several frequencies (8.1, 13, 16) nHz corresponding to the outlier indicating potential candidates of CGW. However, the corresponding BFs are (0.95, 0.33, 0.34) indicating no statistical significance in the analysed data.

We have computed the upper limit using uniform in GW strain prior with and without CRN on the fine frequency grid; the results are present in Fig. 3. The upper limit slightly worsens at low frequencies when we add the CRN; this is understandable as we need a higher amplitude of CGW to get the same SNR when raising the noise floor. The most sensitive frequency of IPTA DR2 is 10.2 nHz with $h_{95} = 9.1 \times 10^{-15}$ in both cases with and without CRN component.

In the model that includes CRN, we are allowing slope and amplitude to vary during the sampling. The recovered posterior

<footnote>$^6$Note that the highest spike corresponds to frequencies close to $f_{yr}$ and should be discarded.</footnote>
Historically, we have performed the search with the noise model with the uniform settings across all pulsars in the array; namely, we have used 30 frequency bins for modelling both RN and DM processes [as it was done in preceding work Arzoumanian et al. (2014); Aggarwal et al. (2019)]. Using this uniform setting we have obtained quite erroneous results, and after a long investigation, we realized that this noise model does not adequately describe the observational data (see Section 5). We have switched to another noise model where we have used the custom-made noise model for the six best EPTA pulsars (see Chalumeau et al. 2021 for details), and for other pulsars, we changed the number of used frequency bins: RN30DM100 – 30 bins for the RN and 100 bins for DM modelling. We will continue using this short-hand notation for the noise model, showing explicitly the number of Fourier frequencies used by the base functions in the Gaussian process describing the corresponding noise. The results presented in the next section were obtained using the ‘custom’ made noise model; we postpone the detailed discussion on the noise model selection and influence until Section 5.

## 3.1 Model selection

During the search, we consider two competing models: noise only and noise with 1 CGW signal. We compute the BF to measure which model is preferred by the IPTA DR2 data set. In particular we employ hyper-model jumps to compute BF following the methods outlined in Hee et al. (2015) and implemented in the enterprise extensions (Taylor et al. 2021). In this method, we introduce a hyperparameter which indexes the models; then, we perform sampling inside each model and in this parameter. The BF then is given by the ratio of the number of iterations that the chain spends in each model. For example, if we consider two models $M_A$ and $M_B$ with the hyperparameter $n$, then the BF is the ratio

$$B_B^A = \frac{n_A}{n_B}, \qquad (22)$$

where $n_A$, $n_B$ counts the number of samples in the chain corresponding to the models A, B. In our previous investigations, we have found that this method gives quite a reliable result, e.g. see Chalumeau et al. (2021) where comparison is done against the evidence computation obtained with Dynasty nested sampling (Skilling 2006; Speagle 2020).

Later in the paper we will compute BF between the models where we also vary the noise.

## 3.2 Upper limit

For obtaining the upper limit, we again use MCMC assuming a uniform prior distribution for $log h$ within the bounds [−18, −11] on the fixed set of $f_{gw}$. We use marginalized posterior probability distribution for the CGW amplitude $p(h)$ to set 95 per cent upper limit $h_{95}$ defined as

$$0.95 = \int_{0}^{h_{95}} p(h)dh. \qquad (23)$$

We use a grid of 100 logarithmically spaced frequencies between $10^{-9}$ and $10^{-7}$ Hz. The lowest bound determined by the IPTA observational span $1/T_{obs}$ that is $\approx 10^{-9}$ Hz while the upper bound, $10^{-7}$ Hz, is constrained by reduction in the sensitivity due to response ($\propto f^{-1/3}$) and by our assumption that GW frequency does not evolve. We should take into account the frequency evolution of the GW signal above $10^{-7}$ Hz during the analysis, as discussed in Petiteau et al. (2013). As mentioned above, we compute the upper limit for

<figure><img src="image_3.png" /><figcaption>Figure 2. On the top panel, the upper bound of the 95 per cent central credible region on $log h$ obtained with a log-uniform prior (with CRN). Bottom panel is the associated BF at given CGW frequency. The black dashed line shows where BF is equal to 1. The two peaks around 35 and 70 nHz with BF > 1 are the 1 and 2 yr and should not be taken into account.</figcaption></figure>

two models of noise; with and without CRN. The inclusion of CRN adds two more parameters to the model (an amplitude parameter $A_{CRN}$ and a spectral index $\gamma_{CRN}$), which we sample together with parameters of CGW and other noise components.

## 4 RESULTS

We compute BF between two models $B_{custom}^{custom+CGW}$ for a noise only model $M_{custom}$ and noise + CGW signal model $M_{custom+CGW}$. The subscript ‘custom’ corresponds to the noise model, we have used for these main results and distinguishes it from other noise models considered in the next section. Similarly to the upper limit run, we have used a log-uniform prior on $h$ bounded by [−18, −11], and the BF was computed on the grid of 100 CGW frequencies ($f_{gw}$) between 1 and 100 nHz.

The main result of this paper is summarized in Fig. 2. In the lower panel we plot the BF, and, as one can see, the noise model is usually preferred. There are few spikes where BF reaches 1, this is definitely not a detection; however, PTA efforts probably should monitor carefully those frequencies in the future extended IPTA data sets. The excess in BF was also used as identification of the CGW candidates to follow. These results show that there is nothing to follow:$^6 B_{custom}^{custom+CGW} \le 1$.

The top panel of Fig. 2 shows the upper bound of the 95 per cent central credible region for CGW strain computed at the same set of frequencies. Note the ‘spiky’ features at several frequencies (8.1, 13, 16) nHz corresponding to the outlier indicating potential candidates of CGW. However, the corresponding BFs are (0.95, 0.33, 0.34) indicating no statistical significance in the analysed data.

We have computed the upper limit using uniform in GW strain prior with and without CRN on the fine frequency grid; the results are present in Fig. 3. The upper limit slightly worsens at low frequencies when we add the CRN; this is understandable as we need a higher amplitude of CGW to get the same SNR when raising the noise floor. The most sensitive frequency of IPTA DR2 is 10.2 nHz with $h_{95} = 9.1 \times 10^{-15}$ in both cases with and without CRN component.

In the model that includes CRN, we are allowing slope and amplitude to vary during the sampling. The recovered posterior

<footnote>$^6$Note that the highest spike corresponds to frequencies close to $f_{yr}$ and should be discarded.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** 95 per cent sky averaged upper limit $h_{95}$ on CGW amplitude for models with and without CRN.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** CRN parameters for IPTA DR2 without CGW signal obtained in (orange) and with CGW signal at 1 nHz (blue).</figcaption></figure>

for CRN with CGW at 1 nHz (blue) is compared to the posterior obtained in Antoniadis et al. (2022) (orange) in Fig. 4. We observe that the amplitude of CRN in the model which does not include CGW (Antoniadis et al. 2022) is slightly higher while the slope is almost the same. This could be explained by a partial absorption of the CRN into CGW. At the same time adding the CRN to the model increases the overall noise level (therefore decreasing the SNR of the signal). The interplay between CRN and CGW appears to mitigate the effect of the CRN on the upper limits in Fig. 4.

We have also overplotted the best CGW upper limit available to date based on the analysis of the NANOGRAV 11 yr data set (Aggarwal et al. 2019) as a dashed (green) line. Note that only the nine-year NANOGRAV data set was included in IPTA DR2. As expected, our current results are better at very low frequencies thanks to the longer observation time. Extended sky coverage, improved effective cadence of observation thanks to overlapping timing data (gaps coverage) and the addition of pulsars like J0437–4715 [only

<table><caption>**Table 2.** Number of frequency bins used for individual RN and DM noise.</caption><thead><tr><td>Pulsar name</td><td>RN bins</td><td>DM bins</td></tr></thead><tbody><tr><td>J0613–0200</td><td>10</td><td>30</td></tr><tr><td>J1012+5307</td><td>150</td><td>30</td></tr><tr><td>J1713+0747</td><td>15</td><td>150</td></tr><tr><td>J1744–1134</td><td>10</td><td>100</td></tr><tr><td>J1909–3744</td><td>10</td><td>100</td></tr><tr><td>J0437–4715</td><td>30</td><td>100</td></tr><tr><td>Other pulsars</td><td>30</td><td>100</td></tr></tbody></table>

present in PPTA data Manchester et al. (2013)], which is an excellent timer, is reflected in a much improved upper limit at 100 nHz, where we might expect the first detection of CGWs.

We want to point at the double-peak feature just above 10 nHz (see Fig. 2). This peak remains the same under the prior change (from the uniform in the amplitude to the uniform in the log-amplitude), which often corresponds to a signal present in the data. The Bayes factors for those peaks are low (0.33, 0.34). None the less, we should keep an eye on those frequencies in the next data releases.

## 5 EFFECT OF NOISE MODELLING ON THE CGW SEARCH

In this section, we consider several noise models and compute BF between those models of noise with and without CGW. The main results presented in the previous section were obtained with the custom-made noise model for the best EPTA pulsars (Chalumeau et al. 2021) and with RN30DM100 choice of Fourier frequencies for other pulsars. The custom noise model modifies the number of Fourier frequencies as given in Table 2 and includes additional noise components (like system noise). However, the Fourier basis for CRN is always fixed at 30 frequency bins.

Here, we present results with what was considered ‘standard’ noise settings before this work, namely RN30DM30 model. We have started analysis using this model, and the search quickly converged to a particular sky position at 51 nHz. The first peculiarity of these results is that 51 nHz is very close to Venus orbital frequency, and the second is that the sky position had a bi-modal structure and was located very close to J1012+5307 see Fig. 5. The Bayes factor for this event with the RN30DM30 noise model was $B_{RN30DM30}^{RN30DM30+CGW} = 18$, which is not very high; however, it seriously triggered our attention by being relatively well constrained in the parameter space.

We have launched a set of investigations trying to understand this event. Using samples taken from the candidate’s posterior, we have checked the contribution to the SNR from each pulsar. It came out that the main part comes from J1012+5307, but a few other pulsars (J1713+0747 and J0437–4715) also contributed not negligibly. We have checked that the zero contribution from a very good timer J1909–3744 is expected given the presumed sky position of the event.

We also conducted a set of injections of CGW signals with parameters taken from the candidate’s posteriors using a simulated IPTA data (same TOAs as the real) with WN only (RN and DM are supposed to be sub-dominant at the candidate frequency). We could not reproduce the observed results with injections even when we increased the amplitude of the simulated signal. More specifically, the sky location posteriors that were recovered for the injected CGW signals were not matching the double blob structure around J1012+5307 observed using the real data.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** 95 per cent sky averaged upper limit $h_{95}$ on CGW amplitude for models with and without CRN.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** CRN parameters for IPTA DR2 without CGW signal obtained in (orange) and with CGW signal at 1 nHz (blue).</figcaption></figure>

for CRN with CGW at 1 nHz (blue) is compared to the posterior obtained in Antoniadis et al. (2022) (orange) in Fig. 4. We observe that the amplitude of CRN in the model which does not include CGW (Antoniadis et al. 2022) is slightly higher while the slope is almost the same. This could be explained by a partial absorption of the CRN into CGW. At the same time adding the CRN to the model increases the overall noise level (therefore decreasing the SNR of the signal). The interplay between CRN and CGW appears to mitigate the effect of the CRN on the upper limits in Fig. 4.

We have also overplotted the best CGW upper limit available to date based on the analysis of the NANOGRAV 11 yr data set (Aggarwal et al. 2019) as a dashed (green) line. Note that only the nine-year NANOGRAV data set was included in IPTA DR2. As expected, our current results are better at very low frequencies thanks to the longer observation time. Extended sky coverage, improved effective cadence of observation thanks to overlapping timing data (gaps coverage) and the addition of pulsars like J0437–4715 [only

<table><caption>**Table 2.** Number of frequency bins used for individual RN and DM noise.</caption><thead><tr><td>Pulsar name</td><td>RN bins</td><td>DM bins</td></tr></thead><tbody><tr><td>J0613–0200</td><td>10</td><td>30</td></tr><tr><td>J1012+5307</td><td>150</td><td>30</td></tr><tr><td>J1713+0747</td><td>15</td><td>150</td></tr><tr><td>J1744–1134</td><td>10</td><td>100</td></tr><tr><td>J1909–3744</td><td>10</td><td>100</td></tr><tr><td>J0437–4715</td><td>30</td><td>100</td></tr><tr><td>Other pulsars</td><td>30</td><td>100</td></tr></tbody></table>

present in PPTA data Manchester et al. (2013)], which is an excellent timer, is reflected in a much improved upper limit at 100 nHz, where we might expect the first detection of CGWs.

We want to point at the double-peak feature just above 10 nHz (see Fig. 2). This peak remains the same under the prior change (from the uniform in the amplitude to the uniform in the log-amplitude), which often corresponds to a signal present in the data. The Bayes factors for those peaks are low (0.33, 0.34). None the less, we should keep an eye on those frequencies in the next data releases.

## 5 EFFECT OF NOISE MODELLING ON THE CGW SEARCH

In this section, we consider several noise models and compute BF between those models of noise with and without CGW. The main results presented in the previous section were obtained with the custom-made noise model for the best EPTA pulsars (Chalumeau et al. 2021) and with RN30DM100 choice of Fourier frequencies for other pulsars. The custom noise model modifies the number of Fourier frequencies as given in Table 2 and includes additional noise components (like system noise). However, the Fourier basis for CRN is always fixed at 30 frequency bins.

Here, we present results with what was considered ‘standard’ noise settings before this work, namely RN30DM30 model. We have started analysis using this model, and the search quickly converged to a particular sky position at 51 nHz. The first peculiarity of these results is that 51 nHz is very close to Venus orbital frequency, and the second is that the sky position had a bi-modal structure and was located very close to J1012+5307 see Fig. 5. The Bayes factor for this event with the RN30DM30 noise model was $B_{RN30DM30}^{RN30DM30+CGW} = 18$, which is not very high; however, it seriously triggered our attention by being relatively well constrained in the parameter space.

We have launched a set of investigations trying to understand this event. Using samples taken from the candidate’s posterior, we have checked the contribution to the SNR from each pulsar. It came out that the main part comes from J1012+5307, but a few other pulsars (J1713+0747 and J0437–4715) also contributed not negligibly. We have checked that the zero contribution from a very good timer J1909–3744 is expected given the presumed sky position of the event.

We also conducted a set of injections of CGW signals with parameters taken from the candidate’s posteriors using a simulated IPTA data (same TOAs as the real) with WN only (RN and DM are supposed to be sub-dominant at the candidate frequency). We could not reproduce the observed results with injections even when we increased the amplitude of the simulated signal. More specifically, the sky location posteriors that were recovered for the injected CGW signals were not matching the double blob structure around J1012+5307 observed using the real data.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Posterior distribution for $\log h$, $\cos\theta$, and $\phi$ using noise model $M_{\text{RN30DM}30}$. The black star indicates the sky position of J1012+5307. The dashed lines represent the median values of the parameters. The quantile values correspond to [0.16, 0.5, 0.84].</figcaption></figure>

BAYESEPHEM does not include a contribution from the inner planets to the phenomenological model, so we consulted a group from Observatoire de la Côte d'Azur (INPOP group), inquiring if there could be an error in the Venus orbit picked up in the CGW analysis. We were reassured that the Venus parameters are known with very high precision: this is a simple planet without any moons. However, looking into the future, we probably should extend BAYESEPHEM (or alternative mitigation models) by including the perturbation of orbital elements for the inner planets. We have performed the analysis with a narrow prior around this event with the model including the pulsar term. The Bayes factor has slightly increased $BF = 21$; however, the parameter estimation did not change appreciably. We have tried several runs with an extended model that includes the orbital eccentricity (again with a narrow prior). The eccentric runs show a very poor convergence, but all of them suggest a relatively high eccentricity. Results of both models (eccentric and with pulsar term) imply that the power is not localized at one particular frequency but spread over some finite-size frequency band.

We have turned to several noise models. We have started with varying number of Fourier frequencies used in the Gaussian process for RN and DM and tried $M_{\text{RN30DM}100}$ and $M_{\text{RN100DM}100}$. The BF for CGW with those noise models has increased tremendously (by a factor 100–1000). Finally, we have tried the custom model for the six best EPTA pulsars. Most notable is the peculiar noise model for J1012+5307 (see table 2) which, in combination with the sky position of the candidate event (being next to it), suggests that the explanation might be in the time-correlated high-frequency noise present in that pulsar (see appendix of Chalumeau et al. 2021). This would also be consistent with the results of eccentric runs suggesting that this could be an extended frequency feature. We do not know the origin of that noise; it was found empirically. Including the custom noise model reduced the BF for this event to $B_{\text{custom+CGW}}^{\text{custom}} = 0.95$ and the posterior samples are not anymore constrained (see Fig. 6 and,

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Posterior distribution for $\log h$, $\cos\theta$, and $\phi$ using noise model $M_{\text{custom}}$. The black star indicates the sky position of J1012+5307. The dashed lines represent the median values of the parameters. The quantile values correspond to [0.16, 0.5, 0.84].</figcaption></figure>

please, note that we have used different range for the parameters). We have also computed the Bayes factor for the noise-only models (assuming that GW signal is weak) and found that $B_{\text{custom}}^{\text{RN30DM}30} = < 10^{-2}$ and similar result for RN30DM100 model suggesting that the data prefers by far the custom model.

Custom made noise models for PPTA pulsars were studied in Goncharov et al. (2020). However, the peculiar behaviour that we have found in IPTA data was mainly originating from J1012+5307 which is not observed by PPTA. For that reason, we have chosen to focus on noise models from Chalumeau et al. (2021). On the other hand, the noise model for J0437–4715 was based on Goncharov et al. (2020) where the spectral index for RN is of 3 and the optimal number of frequency bins for DM is (at least) 91. For this pulsar, we chose 30 frequency bins for RN as it was the recommended value for spectral index > 1.5 (see Goncharov et al. 2021).

This was a useful exercise that triggered a set of investigations we would have to do in case of any CGW candidate. In addition, this section shows the importance of custom modelling noise for the best pulsars in the array, especially the noise at high frequencies, which is often partially neglected, assuming that it is dominated by the white (measurements) noise and it does not affect the search for the stochastic GW signal (which is most pronounced at low frequencies). The considered event shows how unmodelled high-frequency noise could conspire for CGW signal.

# 6 CONCLUSION

We have searched for a continuous GW signal in the IPTA DR2 data set. We have used the Bayesian approach and based detection criteria on the Bayes factor. We have shown that using a custom noise model for the six best EPTA pulsars is essential for the correct interpretation of the data. This is especially true for J1012+5307, which exhibits time-correlated noise at high frequencies. We found
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Posterior distribution for $\log h$, $\cos\theta$, and $\phi$ using noise model $M_{\text{RN30DM}30}$. The black star indicates the sky position of J1012+5307. The dashed lines represent the median values of the parameters. The quantile values correspond to [0.16, 0.5, 0.84].</figcaption></figure>

BAYESEPHEM does not include a contribution from the inner planets to the phenomenological model, so we consulted a group from Observatoire de la Côte d'Azur (INPOP group), inquiring if there could be an error in the Venus orbit picked up in the CGW analysis. We were reassured that the Venus parameters are known with very high precision: this is a simple planet without any moons. However, looking into the future, we probably should extend BAYESEPHEM (or alternative mitigation models) by including the perturbation of orbital elements for the inner planets. We have performed the analysis with a narrow prior around this event with the model including the pulsar term. The Bayes factor has slightly increased $BF = 21$; however, the parameter estimation did not change appreciably. We have tried several runs with an extended model that includes the orbital eccentricity (again with a narrow prior). The eccentric runs show a very poor convergence, but all of them suggest a relatively high eccentricity. Results of both models (eccentric and with pulsar term) imply that the power is not localized at one particular frequency but spread over some finite-size frequency band.

We have turned to several noise models. We have started with varying number of Fourier frequencies used in the Gaussian process for RN and DM and tried $M_{\text{RN30DM}100}$ and $M_{\text{RN100DM}100}$. The BF for CGW with those noise models has increased tremendously (by a factor 100–1000). Finally, we have tried the custom model for the six best EPTA pulsars. Most notable is the peculiar noise model for J1012+5307 (see table 2) which, in combination with the sky position of the candidate event (being next to it), suggests that the explanation might be in the time-correlated high-frequency noise present in that pulsar (see appendix of Chalumeau et al. 2021). This would also be consistent with the results of eccentric runs suggesting that this could be an extended frequency feature. We do not know the origin of that noise; it was found empirically. Including the custom noise model reduced the BF for this event to $B_{\text{custom+CGW}}^{\text{custom}} = 0.95$ and the posterior samples are not anymore constrained (see Fig. 6 and,

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Posterior distribution for $\log h$, $\cos\theta$, and $\phi$ using noise model $M_{\text{custom}}$. The black star indicates the sky position of J1012+5307. The dashed lines represent the median values of the parameters. The quantile values correspond to [0.16, 0.5, 0.84].</figcaption></figure>

please, note that we have used different range for the parameters). We have also computed the Bayes factor for the noise-only models (assuming that GW signal is weak) and found that $B_{\text{custom}}^{\text{RN30DM}30} = < 10^{-2}$ and similar result for RN30DM100 model suggesting that the data prefers by far the custom model.

Custom made noise models for PPTA pulsars were studied in Goncharov et al. (2020). However, the peculiar behaviour that we have found in IPTA data was mainly originating from J1012+5307 which is not observed by PPTA. For that reason, we have chosen to focus on noise models from Chalumeau et al. (2021). On the other hand, the noise model for J0437–4715 was based on Goncharov et al. (2020) where the spectral index for RN is of 3 and the optimal number of frequency bins for DM is (at least) 91. For this pulsar, we chose 30 frequency bins for RN as it was the recommended value for spectral index > 1.5 (see Goncharov et al. 2021).

This was a useful exercise that triggered a set of investigations we would have to do in case of any CGW candidate. In addition, this section shows the importance of custom modelling noise for the best pulsars in the array, especially the noise at high frequencies, which is often partially neglected, assuming that it is dominated by the white (measurements) noise and it does not affect the search for the stochastic GW signal (which is most pronounced at low frequencies). The considered event shows how unmodelled high-frequency noise could conspire for CGW signal.

# 6 CONCLUSION

We have searched for a continuous GW signal in the IPTA DR2 data set. We have used the Bayesian approach and based detection criteria on the Bayes factor. We have shown that using a custom noise model for the six best EPTA pulsars is essential for the correct interpretation of the data. This is especially true for J1012+5307, which exhibits time-correlated noise at high frequencies. We found
no CGW candidates using this noise model and proceeded to set the upper limit on GW strain. The addition of CRN in the noise model slightly affects the upper limit by lowering the sensitivity of the array at low frequencies. The most sensitive frequency appears to be around 10 nHz with a 95 per cent sky averaged upper limit for CGW amplitude $h_{95} = 9.1 \times 10^{-15}$ $\text{m}$. The IPTA DR2 shows a much better upper limit than previously set at higher frequencies, making it a promising data set to detect CGW.

During the analysis, we demonstrated the CGW candidate follow-up investigations program, which was an important exercise that should be used in subsequent PTA CGW analysis. The expected CGW signal has low SNR, and its SNR will be slowly accumulated as we get more pulsars and a longer observational span. Modelling noise in pulsar data is essential, especially at high frequencies.

This analysis was limited to circular SMBHBs using only the Earth term. The use of eccentric CGW signal and including the pulsar term might potentially improve the search; however, it brings signal complexity which might make harder the interpretation of the results and increases the parameter space. We are entering the era of very high quality and high cadence radio observations with new instruments like FAST (Hobbs et al. 2019) or SKA (Stappers et al. 2018) with sophisticated data analysis techniques. Additional investigations of the best approach to detecting CGW have to be re-investigated, probably using simulated data and/or an extended CGW signal injection campaign.

## ACKNOWLEDGEMENTS

The IPTA is a consortium of existing PTA collaborations, namely, the EPTA, NANOGrav, PPTA, and the recent addition of the InPTA. Observing collaborations from China and South Africa are also part of the IPTA.

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

The NANOGrav Physics Frontiers Center is supported by the National Science Foundation (NSF) Physics Frontier Center award numbers 1430284 and 2020265. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The

Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Arecibo Observatory is a facility of the National Science Foundation operated under cooperative agreement by the University of Central Florida in alliance with Yang Enterprises, Inc. and Universidad Metropolitana.

The Parkes radio telescope (Murriyang) is part of the Australia Telescope which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO.

JA acknowledges support by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of 'Science and Society' Action Always strive for excellence – 'Theodoros Papazoglou' (Project Number: 01431). SBS acknowledges generous support by the NSF through grant AST-1815664. The work is supported by National SKA program of China 2020SKA0120100, Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. JACC was supported in part by NASA CT Space Grant PTE Federal Award Number 80NSSC20M0129. CMFM and JACC are also supported by the National Science Foundations NANOGrav Physics Frontier Center, Award Number 2020265. The Center for Computational Astrophysics is a division of the Flatiron Institute in New York City, which is supported by the Simons Foundation. This research was supported in part by the National Science Foundation grant AST-2106552. AC acknowledges support from the Paris Île-de-France Region. Support for HTC was provided by NASA through the NASA Hubble Fellowship Program grant HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. SD is the recipient of an Australian Research Council Discovery Early Career Award (DE210101738) funded by the Australian Government. GD, RK, and MKr acknowledge support from European Research Council (ERC) Synergy Grant 'BlackHoleCam' Grant Agreement Number 610058 and ERC Advanced Grant 'LEAP' Grant Agreement Number 337062. TD and MTL acknowledge support received from NSF AAG award number 200968. ECF is supported by NASA under award number 80GSFC17M0002.002. BG is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN. Portions of this work performed at the Naval Research Laboratory is supported by NASA and ONR 6.1 basic research funding. 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. JWM gratefully acknowledges support by the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference CITA 490888-16]. KDO was supported in part by NSF Grant No. 2207267. ASa, ASe, and GS acknowledge financial support provided under the European Union's H2020 ERC Consolidator Grant 'Binary Massive Black Hole Astrophysic' (B Massive, Grant Agreement: 818691). RMS acknowledges support through Australian Research Council Future Fellowship FT190100155. JJS is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388, and this research was supported in part by NSF AST-1847938. This research was funded partially by the Australian Government through the Australian Research Council (ARC), grants CE170100004 (OzGrav) and FL150100148. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. SRT acknowledges support from NSF grants AST-2007993 and PHY-2020265. SRT also acknowledges support from the Vanderbilt University
no CGW candidates using this noise model and proceeded to set the upper limit on GW strain. The addition of CRN in the noise model slightly affects the upper limit by lowering the sensitivity of the array at low frequencies. The most sensitive frequency appears to be around 10 nHz with a 95 per cent sky averaged upper limit for CGW amplitude $h_{95} = 9.1 \times 10^{-15}$ $\text{m}$. The IPTA DR2 shows a much better upper limit than previously set at higher frequencies, making it a promising data set to detect CGW.

During the analysis, we demonstrated the CGW candidate follow-up investigations program, which was an important exercise that should be used in subsequent PTA CGW analysis. The expected CGW signal has low SNR, and its SNR will be slowly accumulated as we get more pulsars and a longer observational span. Modelling noise in pulsar data is essential, especially at high frequencies.

This analysis was limited to circular SMBHBs using only the Earth term. The use of eccentric CGW signal and including the pulsar term might potentially improve the search; however, it brings signal complexity which might make harder the interpretation of the results and increases the parameter space. We are entering the era of very high quality and high cadence radio observations with new instruments like FAST (Hobbs et al. 2019) or SKA (Stappers et al. 2018) with sophisticated data analysis techniques. Additional investigations of the best approach to detecting CGW have to be re-investigated, probably using simulated data and/or an extended CGW signal injection campaign.

## ACKNOWLEDGEMENTS

The IPTA is a consortium of existing PTA collaborations, namely, the EPTA, NANOGrav, PPTA, and the recent addition of the InPTA. Observing collaborations from China and South Africa are also part of the IPTA.

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

The NANOGrav Physics Frontiers Center is supported by the National Science Foundation (NSF) Physics Frontier Center award numbers 1430284 and 2020265. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The

Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Arecibo Observatory is a facility of the National Science Foundation operated under cooperative agreement by the University of Central Florida in alliance with Yang Enterprises, Inc. and Universidad Metropolitana.

The Parkes radio telescope (Murriyang) is part of the Australia Telescope which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO.

JA acknowledges support by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of 'Science and Society' Action Always strive for excellence – 'Theodoros Papazoglou' (Project Number: 01431). SBS acknowledges generous support by the NSF through grant AST-1815664. The work is supported by National SKA program of China 2020SKA0120100, Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. JACC was supported in part by NASA CT Space Grant PTE Federal Award Number 80NSSC20M0129. CMFM and JACC are also supported by the National Science Foundations NANOGrav Physics Frontier Center, Award Number 2020265. The Center for Computational Astrophysics is a division of the Flatiron Institute in New York City, which is supported by the Simons Foundation. This research was supported in part by the National Science Foundation grant AST-2106552. AC acknowledges support from the Paris Île-de-France Region. Support for HTC was provided by NASA through the NASA Hubble Fellowship Program grant HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. SD is the recipient of an Australian Research Council Discovery Early Career Award (DE210101738) funded by the Australian Government. GD, RK, and MKr acknowledge support from European Research Council (ERC) Synergy Grant 'BlackHoleCam' Grant Agreement Number 610058 and ERC Advanced Grant 'LEAP' Grant Agreement Number 337062. TD and MTL acknowledge support received from NSF AAG award number 200968. ECF is supported by NASA under award number 80GSFC17M0002.002. BG is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN. Portions of this work performed at the Naval Research Laboratory is supported by NASA and ONR 6.1 basic research funding. 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. JWM gratefully acknowledges support by the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference CITA 490888-16]. KDO was supported in part by NSF Grant No. 2207267. ASa, ASe, and GS acknowledge financial support provided under the European Union's H2020 ERC Consolidator Grant 'Binary Massive Black Hole Astrophysic' (B Massive, Grant Agreement: 818691). RMS acknowledges support through Australian Research Council Future Fellowship FT190100155. JJS is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388, and this research was supported in part by NSF AST-1847938. This research was funded partially by the Australian Government through the Australian Research Council (ARC), grants CE170100004 (OzGrav) and FL150100148. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. SRT acknowledges support from NSF grants AST-2007993 and PHY-2020265. SRT also acknowledges support from the Vanderbilt University
College of Arts & Science Dean's Faculty Fellowship program. Multiple NANOGrav members acknowledge support of NSF Physics Frontiers Center awards 1430284 and 2020265. SMR is a CIFAR Fellow. AV acknowledges the support of the Royal Society and Wolfson Foundation. JPWV acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039). MED acknowledges support from the National Science Foundation (NSF) Physics Frontier Center award 1430284, and from the Naval Research Laboratory by NASA under contract S-15633Y. ZCC is supported by the National Natural Science Foundation of China (Grant No. 12247176) and the China Postdoctoral Science Foundation Fellowship No. 2022M710429. CAW acknowledges support from CIERA, the Adler Planetarium, and the Brinson Foundation through a CIERA-Adler postdoctoral fellowship.

## DATA AVAILABILITY

The timing data used in this article is available on the IPTA website https://ipta4gw.org/data-release/ (second data release: https://gitlab.com/IPTA/DR2).

## REFERENCES

Aggarwal K. et al., 2019, ApJ, 880, 116

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

Arzoumanian Z. et al., 2014, ApJ, 794, 141

Arzoumanian Z. et al., 2016, ApJ, 821, L13

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

Arzoumanian Z. et al., 2023, preprint (arXiv:2301.03608)

Babak S., Sesana A., 2012, Phys. Rev. D, 85, 044034

Babak S. et al., 2016, MNRAS, 455, 1665

Becsy B., Cornish N. J., 2020, Class. Quantum Gravity, 37, 135011

Becsy B., Cornish N. J., Digman M. C., 2022, Phys. Rev. D, 105, 122003

Chalumeau A., et al., 2021, MNRAS, 509, 5538

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

Corbin V., Cornish N. J., 2010, preprint (arXiv:1008.1782) https://arxiv.org/abs/1008.1782

Desvignes G. et al., 2016, MNRAS, 458, 3341

Detweiler S., 1979, ApJ, 234, 1100

Edwards R. T., Hobbs G. B., Manchester R. N., 2006, MNRAS, 372, 1549

Ellis J., van Haasteren R., 2017, jellis18/PTMCMCSampler: Official Release. https://doi.org/10.5281/zenodo.1037579

Ellis J. A., Siemens X., Creighton J. D. E., 2012, ApJ, 756, 175

Ellis J. A., Vallisneri M., Taylor S. R., Baker P. T., 2020, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust Pulsar Inference SuitE. Zenodo, https://doi.org/10.5281/zenodo.4059815

Folkner W. M., Park R. S., 2018, JPL Inter-Office Memorandum

Goncharov B. et al., 2020, MNRAS, 502, 478

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

Hee S., Handley W. J., Hobson M. P., Lasenby A. N., 2015, MNRAS, 455, 2461

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

Hobbs G., Dai S., Manchester R. N., Shannon R. M., Kerr M., Lee K.-J., Xu R.-X., 2019, Res. Astron. Astrophys., 19, 020

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

Lam M. T. et al., 2016, ApJ, 834, 35

Maiorano M., Paolis F. D., Nucita A. A., 2021, Symmetry, 13, 2418

Manchester R. N. et al., 2013, PASA, 30, e017

NANOGrav Collaboration, 2015, ApJ, 813, 65

Perera B. B. P., et al., 2019, MNRAS, 490, 4666

Petiteau A., Babak S., Sesana A., de Araújo M., 2013, Phys. Rev. D, 87, 064036

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

Sazhin M., 1978, Sov. Astron., 22, 36

Sesana A., Vecchio A., 2010, 104008 Phys. Rev. D, 81

Sesana A., Vecchio A., Colacino C. N., 2008, MNRAS, 390, 192

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

Skilling J., 2006, Bayesian Anal., 1, 833

Speagle J. S., 2020, MNRAS, 493, 3132

Stappers B., Keane E., Kramer M., Possenti A., Stairs I., 2018, Phil. Trans. R. Soc. A Math. Phys. Eng. Sci., 376, 20170293

Speri L. et al., 2022, MNRAS, 518, 1802

Tarafdar P. et al., 2022, The Indian Pulsar Timing Array: First data release., https://arxiv.org/abs/2206.09289

Taylor S. R., Vallisneri M., Ellis J. A., Mingarelli C. M. F., Lazio T. J. W., van Haasteren R., 2016, ApJ, 819, L6

Taylor S. R., Baker P. T., Hazboun J. S., Simon J., Vigeland S. J., 2021, Enterprise extensions. https://github.com/nanograv/enterprise_extensions

Vallisneri M., et al., 2020, ApJ, 893, 112

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

Zhu X. J. et al., 2014, MNRAS, 444, 3709

<footnote n="1">¹CNRS, Astroparticule et Cosmologie, Université de Paris, F-75013 Paris, France</footnote>

<footnote n="2">²Department of Physics and Astronomy, Widener University, One University Place, Chester, PA 19013, USA</footnote>

<footnote n="3">³Department of Physics, Montana State University, Bozeman, MT 59717, USA</footnote>

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

<footnote n="5">⁵Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, P. R. China</footnote>

<footnote n="6">⁶Department of Astronomy, Beijing Normal University, Beijing 100875, China</footnote>

<footnote n="7">⁷Station de Radioastronomie de Nançay, Observatoire de Paris, PSL University, CNRS, Université d'Orléans, F-18330 Nançay, France</footnote>

<footnote n="8">⁸Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E UMR7328, Université d'Orléans, CNRS, F-45071 Orléans, France</footnote>

<footnote n="9">⁹Department of Physics, Oregon State University, Corvallis, OR 97331, USA</footnote>

<footnote n="10">¹⁰Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>

<footnote n="11">¹¹Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA</footnote>

<footnote n="12">¹²Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="13">¹³IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="14">¹⁴Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA</footnote>

<footnote n="15">¹⁵INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="16">¹⁶Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA</footnote>

<footnote n="17">¹⁷Laboratoire Univers et Théories LUTh, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France</footnote>

<footnote n="18">¹⁸Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="19">¹⁹Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, PO Box 413, Milwaukee, WI 53201, USA</footnote>

<footnote n="20">²⁰Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208, USA</footnote>

<footnote n="21">²¹Adler Planetarium, 1300 S. DuSable Lake Shore Dr, Chicago, IL 60605, USA</footnote>

<footnote n="22">²²Advanced Institute of Natural Sciences, Beijing Normal University, Zhuhai 519087, China</footnote>

<footnote n="23">²³Institute of Astrophysics, FORTH, N. Plastira 100, 70013 Heraklion, Greece</footnote>

<footnote n="24">²⁴X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA</footnote>

<footnote n="25">²⁵Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, Victoria 3122, Australia, ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Australia</footnote>
College of Arts & Science Dean's Faculty Fellowship program. Multiple NANOGrav members acknowledge support of NSF Physics Frontiers Center awards 1430284 and 2020265. SMR is a CIFAR Fellow. AV acknowledges the support of the Royal Society and Wolfson Foundation. JPWV acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039). MED acknowledges support from the National Science Foundation (NSF) Physics Frontier Center award 1430284, and from the Naval Research Laboratory by NASA under contract S-15633Y. ZCC is supported by the National Natural Science Foundation of China (Grant No. 12247176) and the China Postdoctoral Science Foundation Fellowship No. 2022M710429. CAW acknowledges support from CIERA, the Adler Planetarium, and the Brinson Foundation through a CIERA-Adler postdoctoral fellowship.

## DATA AVAILABILITY

The timing data used in this article is available on the IPTA website https://ipta4gw.org/data-release/ (second data release: https://gitlab.com/IPTA/DR2).

## REFERENCES

Aggarwal K. et al., 2019, ApJ, 880, 116

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

Arzoumanian Z. et al., 2014, ApJ, 794, 141

Arzoumanian Z. et al., 2016, ApJ, 821, L13

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

Arzoumanian Z. et al., 2023, preprint (arXiv:2301.03608)

Babak S., Sesana A., 2012, Phys. Rev. D, 85, 044034

Babak S. et al., 2016, MNRAS, 455, 1665

Becsy B., Cornish N. J., 2020, Class. Quantum Gravity, 37, 135011

Becsy B., Cornish N. J., Digman M. C., 2022, Phys. Rev. D, 105, 122003

Chalumeau A., et al., 2021, MNRAS, 509, 5538

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

Corbin V., Cornish N. J., 2010, preprint (arXiv:1008.1782) https://arxiv.org/abs/1008.1782

Desvignes G. et al., 2016, MNRAS, 458, 3341

Detweiler S., 1979, ApJ, 234, 1100

Edwards R. T., Hobbs G. B., Manchester R. N., 2006, MNRAS, 372, 1549

Ellis J., van Haasteren R., 2017, jellis18/PTMCMCSampler: Official Release. https://doi.org/10.5281/zenodo.1037579

Ellis J. A., Siemens X., Creighton J. D. E., 2012, ApJ, 756, 175

Ellis J. A., Vallisneri M., Taylor S. R., Baker P. T., 2020, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust Pulsar Inference SuitE. Zenodo, https://doi.org/10.5281/zenodo.4059815

Folkner W. M., Park R. S., 2018, JPL Inter-Office Memorandum

Goncharov B. et al., 2020, MNRAS, 502, 478

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

Hee S., Handley W. J., Hobson M. P., Lasenby A. N., 2015, MNRAS, 455, 2461

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

Hobbs G., Dai S., Manchester R. N., Shannon R. M., Kerr M., Lee K.-J., Xu R.-X., 2019, Res. Astron. Astrophys., 19, 020

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

Lam M. T. et al., 2016, ApJ, 834, 35

Maiorano M., Paolis F. D., Nucita A. A., 2021, Symmetry, 13, 2418

Manchester R. N. et al., 2013, PASA, 30, e017

NANOGrav Collaboration, 2015, ApJ, 813, 65

Perera B. B. P., et al., 2019, MNRAS, 490, 4666

Petiteau A., Babak S., Sesana A., de Araújo M., 2013, Phys. Rev. D, 87, 064036

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

Sazhin M., 1978, Sov. Astron., 22, 36

Sesana A., Vecchio A., 2010, 104008 Phys. Rev. D, 81

Sesana A., Vecchio A., Colacino C. N., 2008, MNRAS, 390, 192

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

Skilling J., 2006, Bayesian Anal., 1, 833

Speagle J. S., 2020, MNRAS, 493, 3132

Stappers B., Keane E., Kramer M., Possenti A., Stairs I., 2018, Phil. Trans. R. Soc. A Math. Phys. Eng. Sci., 376, 20170293

Speri L. et al., 2022, MNRAS, 518, 1802

Tarafdar P. et al., 2022, The Indian Pulsar Timing Array: First data release., https://arxiv.org/abs/2206.09289

Taylor S. R., Vallisneri M., Ellis J. A., Mingarelli C. M. F., Lazio T. J. W., van Haasteren R., 2016, ApJ, 819, L6

Taylor S. R., Baker P. T., Hazboun J. S., Simon J., Vigeland S. J., 2021, Enterprise extensions. https://github.com/nanograv/enterprise_extensions

Vallisneri M., et al., 2020, ApJ, 893, 112

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

Zhu X. J. et al., 2014, MNRAS, 444, 3709

<footnote n="1">¹CNRS, Astroparticule et Cosmologie, Université de Paris, F-75013 Paris, France</footnote>

<footnote n="2">²Department of Physics and Astronomy, Widener University, One University Place, Chester, PA 19013, USA</footnote>

<footnote n="3">³Department of Physics, Montana State University, Bozeman, MT 59717, USA</footnote>

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

<footnote n="5">⁵Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, P. R. China</footnote>

<footnote n="6">⁶Department of Astronomy, Beijing Normal University, Beijing 100875, China</footnote>

<footnote n="7">⁷Station de Radioastronomie de Nançay, Observatoire de Paris, PSL University, CNRS, Université d'Orléans, F-18330 Nançay, France</footnote>

<footnote n="8">⁸Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E UMR7328, Université d'Orléans, CNRS, F-45071 Orléans, France</footnote>

<footnote n="9">⁹Department of Physics, Oregon State University, Corvallis, OR 97331, USA</footnote>

<footnote n="10">¹⁰Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>

<footnote n="11">¹¹Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA</footnote>

<footnote n="12">¹²Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="13">¹³IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="14">¹⁴Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA</footnote>

<footnote n="15">¹⁵INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="16">¹⁶Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA</footnote>

<footnote n="17">¹⁷Laboratoire Univers et Théories LUTh, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France</footnote>

<footnote n="18">¹⁸Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="19">¹⁹Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, PO Box 413, Milwaukee, WI 53201, USA</footnote>

<footnote n="20">²⁰Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208, USA</footnote>

<footnote n="21">²¹Adler Planetarium, 1300 S. DuSable Lake Shore Dr, Chicago, IL 60605, USA</footnote>

<footnote n="22">²²Advanced Institute of Natural Sciences, Beijing Normal University, Zhuhai 519087, China</footnote>

<footnote n="23">²³Institute of Astrophysics, FORTH, N. Plastira 100, 70013 Heraklion, Greece</footnote>

<footnote n="24">²⁴X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA</footnote>

<footnote n="25">²⁵Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, Victoria 3122, Australia, ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Australia</footnote>
<footnote n="26"><sup>26</sup>ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Australia</footnote>

<footnote n="27"><sup>27</sup>International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia</footnote>

<footnote n="28"><sup>28</sup>Department of Physics, University of Florida, 2001 Museum Road, Gainesville, FL 32611-8440, USA</footnote>

<footnote n="29"><sup>29</sup>Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA</footnote>

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

<footnote n="31"><sup>31</sup>Department of Physics and Astronomy, Franklin & Marshall College, PO Box 3003, Lancaster, PA 17604, USA</footnote>

<footnote n="32"><sup>32</sup>Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="33"><sup>33</sup>School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia</footnote>

<footnote n="34"><sup>34</sup>George Mason University, Resident at the Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="35"><sup>35</sup>National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA</footnote>

<footnote n="36"><sup>36</sup>Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA</footnote>

<footnote n="37"><sup>37</sup>Eureka Scientific, 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA</footnote>

<footnote n="38"><sup>38</sup>Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA</footnote>

<footnote n="39"><sup>39</sup>Research Center for Intelligent Computing Platforms, Zhejiang Laboratory, Hangzhou 311100, China</footnote>

<footnote n="40"><sup>40</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="41"><sup>41</sup>Department of Physics and Astronomy, West Virginia University, PO Box 6315, Morgantown, WV 26506, USA</footnote>

<footnote n="42"><sup>42</sup>VWU Center for Gravitational Waves and Cosmology, White Hall, Morgantown, WV 26506, USA</footnote>

<footnote n="43"><sup>43</sup>Gran Sasso Science Institute, Viale Francesco Crispi, 7, I-67100 L'Aquila AQ, Italy</footnote>

<footnote n="44"><sup>44</sup>Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI 48109, USA</footnote>

<footnote n="45"><sup>45</sup>CSIRO Space and Astronomy, Australia Telescope National Facility, PO Box 76, Epping NSW 1710, Australia</footnote>

<footnote n="46"><sup>46</sup>Theoretical AstroPhysics Including Relativity (TAPIR), California Institute of Technology, MC 350-17, Pasadena, CA 91125, USA</footnote>

<footnote n="47"><sup>47</sup>Dominion Radio Astrophysical Observatory, Herzberg Research Centre for Astronomy and Astrophysics, National 27 Research Council Canada, PO Box 248, Penticton, BC V2A 6J9, Canada</footnote>

<footnote n="48"><sup>48</sup>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="49"><sup>49</sup>Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>

<footnote n="50"><sup>50</sup>Department of Physics, University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA</footnote>

<footnote n="51"><sup>51</sup>Oregon State University, Department of Physics, 301 Weniger Hall, Corvallis, OR 97331, USA</footnote>

<footnote n="52"><sup>52</sup>Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="53"><sup>53</sup>Green Bank Observatory, PO Box 2, Green Bank, WV 24944, USA</footnote>

<footnote n="54"><sup>54</sup>Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA</footnote>

<footnote n="55"><sup>55</sup>E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull HU6 7RX, UK</footnote>

<footnote n="56"><sup>56</sup>Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull HU6 7RX, UK</footnote>

<footnote n="57"><sup>57</sup>Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="58"><sup>58</sup>Department of Physics, Lafayette College, Easton, PA 18042, USA</footnote>

<footnote n="59"><sup>59</sup>Institute of Cosmology, Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA</footnote>

<footnote n="60"><sup>60</sup>Manly Astrophysics, 15/41-42 East Esplanade, Manly, NSW 2095, Australia</footnote>

<footnote n="61"><sup>61</sup>Institute of Physics, Eötvös Loránd University, Pázmány P.s. 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="62"><sup>62</sup>Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA</footnote>

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

<footnote n="64"><sup>64</sup>National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA</footnote>

<footnote n="65"><sup>65</sup>U.S. Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="66"><sup>66</sup>CSIRO Scientific Computing, Australian Technology Park, Locked Bag 9013, Alexandria, NSW 1435, Australia</footnote>

<footnote n="67"><sup>67</sup>CSIRO Parkes Observatory, PO Box 276, Parkes, NSW 2870, Australia</footnote>

<footnote n="68"><sup>68</sup>Giant Army, 915A 17th Avenue, Seattle, WA 98122, USA</footnote>

<footnote n="69"><sup>69</sup>Department of Astrophysical and Planetary Sciences, University of Coroado, Boulder, CO 80309, USA</footnote>

<footnote n="70"><sup>70</sup>Center for Astrophysics | Harvard & Smithsonian, USA</footnote>

<footnote n="71"><sup>71</sup>Department of Physics and Astronomy, Swarthmore College, Swarthmore, PA 19081, USA</footnote>

<footnote n="72"><sup>72</sup>Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, D-14476 Potsdam, Germany</footnote>

<footnote n="73"><sup>73</sup>Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA</footnote>

<footnote n="74"><sup>74</sup>Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="75"><sup>75</sup>Fakultät für Physik, Universität Bielefeld, Postfach 100131, D-33501 Bielefeld, Germany</footnote>

<footnote n="76"><sup>76</sup>Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150 Kexueyi St, Xinshi District, Ürümqi, Ürümqi, Xinjiang 830011, China</footnote>

<footnote n="77"><sup>77</sup>Institute of Optoelectronic Technology, Lishui University, Lishui 323000, China</footnote>

<footnote n="78"><sup>78</sup>National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China</footnote>

<footnote n="79"><sup>79</sup>Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023, China</footnote>

<footnote>This paper has been typeset from a TEX/LATEX file prepared by the author.</footnote>
<footnote n="26"><sup>26</sup>ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Australia</footnote>

<footnote n="27"><sup>27</sup>International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia</footnote>

<footnote n="28"><sup>28</sup>Department of Physics, University of Florida, 2001 Museum Road, Gainesville, FL 32611-8440, USA</footnote>

<footnote n="29"><sup>29</sup>Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA</footnote>

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

<footnote n="31"><sup>31</sup>Department of Physics and Astronomy, Franklin & Marshall College, PO Box 3003, Lancaster, PA 17604, USA</footnote>

<footnote n="32"><sup>32</sup>Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="33"><sup>33</sup>School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia</footnote>

<footnote n="34"><sup>34</sup>George Mason University, Resident at the Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="35"><sup>35</sup>National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA</footnote>

<footnote n="36"><sup>36</sup>Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA</footnote>

<footnote n="37"><sup>37</sup>Eureka Scientific, 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA</footnote>

<footnote n="38"><sup>38</sup>Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA</footnote>

<footnote n="39"><sup>39</sup>Research Center for Intelligent Computing Platforms, Zhejiang Laboratory, Hangzhou 311100, China</footnote>

<footnote n="40"><sup>40</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="41"><sup>41</sup>Department of Physics and Astronomy, West Virginia University, PO Box 6315, Morgantown, WV 26506, USA</footnote>

<footnote n="42"><sup>42</sup>VWU Center for Gravitational Waves and Cosmology, White Hall, Morgantown, WV 26506, USA</footnote>

<footnote n="43"><sup>43</sup>Gran Sasso Science Institute, Viale Francesco Crispi, 7, I-67100 L'Aquila AQ, Italy</footnote>

<footnote n="44"><sup>44</sup>Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI 48109, USA</footnote>

<footnote n="45"><sup>45</sup>CSIRO Space and Astronomy, Australia Telescope National Facility, PO Box 76, Epping NSW 1710, Australia</footnote>

<footnote n="46"><sup>46</sup>Theoretical AstroPhysics Including Relativity (TAPIR), California Institute of Technology, MC 350-17, Pasadena, CA 91125, USA</footnote>

<footnote n="47"><sup>47</sup>Dominion Radio Astrophysical Observatory, Herzberg Research Centre for Astronomy and Astrophysics, National 27 Research Council Canada, PO Box 248, Penticton, BC V2A 6J9, Canada</footnote>

<footnote n="48"><sup>48</sup>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="49"><sup>49</sup>Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>

<footnote n="50"><sup>50</sup>Department of Physics, University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA</footnote>

<footnote n="51"><sup>51</sup>Oregon State University, Department of Physics, 301 Weniger Hall, Corvallis, OR 97331, USA</footnote>

<footnote n="52"><sup>52</sup>Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="53"><sup>53</sup>Green Bank Observatory, PO Box 2, Green Bank, WV 24944, USA</footnote>

<footnote n="54"><sup>54</sup>Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA</footnote>

<footnote n="55"><sup>55</sup>E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull HU6 7RX, UK</footnote>

<footnote n="56"><sup>56</sup>Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull HU6 7RX, UK</footnote>

<footnote n="57"><sup>57</sup>Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="58"><sup>58</sup>Department of Physics, Lafayette College, Easton, PA 18042, USA</footnote>

<footnote n="59"><sup>59</sup>Institute of Cosmology, Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA</footnote>

<footnote n="60"><sup>60</sup>Manly Astrophysics, 15/41-42 East Esplanade, Manly, NSW 2095, Australia</footnote>

<footnote n="61"><sup>61</sup>Institute of Physics, Eötvös Loránd University, Pázmány P.s. 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="62"><sup>62</sup>Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612, USA</footnote>

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

<footnote n="64"><sup>64</sup>National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA</footnote>

<footnote n="65"><sup>65</sup>U.S. Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote n="66"><sup>66</sup>CSIRO Scientific Computing, Australian Technology Park, Locked Bag 9013, Alexandria, NSW 1435, Australia</footnote>

<footnote n="67"><sup>67</sup>CSIRO Parkes Observatory, PO Box 276, Parkes, NSW 2870, Australia</footnote>

<footnote n="68"><sup>68</sup>Giant Army, 915A 17th Avenue, Seattle, WA 98122, USA</footnote>

<footnote n="69"><sup>69</sup>Department of Astrophysical and Planetary Sciences, University of Coroado, Boulder, CO 80309, USA</footnote>

<footnote n="70"><sup>70</sup>Center for Astrophysics | Harvard & Smithsonian, USA</footnote>

<footnote n="71"><sup>71</sup>Department of Physics and Astronomy, Swarthmore College, Swarthmore, PA 19081, USA</footnote>

<footnote n="72"><sup>72</sup>Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, D-14476 Potsdam, Germany</footnote>

<footnote n="73"><sup>73</sup>Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA</footnote>

<footnote n="74"><sup>74</sup>Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="75"><sup>75</sup>Fakultät für Physik, Universität Bielefeld, Postfach 100131, D-33501 Bielefeld, Germany</footnote>

<footnote n="76"><sup>76</sup>Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150 Kexueyi St, Xinshi District, Ürümqi, Ürümqi, Xinjiang 830011, China</footnote>

<footnote n="77"><sup>77</sup>Institute of Optoelectronic Technology, Lishui University, Lishui 323000, China</footnote>

<footnote n="78"><sup>78</sup>National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China</footnote>

<footnote n="79"><sup>79</sup>Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023, China</footnote>

<footnote>This paper has been typeset from a TEX/LATEX file prepared by the author.</footnote>