# Periodicity search in the timing of the 25 millisecond pulsars from the second data release of the European Pulsar Timing Array

Iuliana Nitu, Michael Keith, David Champion, Ismaël Cognard, Gregory Desvignes, Lucas Guillemot, Yanjun Guo, Huanchen Hu, Jiwoong Jang, Jedrzej Jawor, et al.

► To cite this version:

Iuliana Nitu, Michael Keith, David Champion, Ismaël Cognard, Gregory Desvignes, et al. Periodicity search in the timing of the 25 millisecond pulsars from the second data release of the European Pulsar Timing Array. Monthly Notices of the Royal Astronomical Society, 2024, 534 (3), pp.1753-1762. 10.1093/mnras/stae2162. hal-04685584

HAL Id: hal-04685584

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

Submitted on 13 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Periodicity search in the timing of the 25 millisecond pulsars from the second data release of the European Pulsar Timing Array

Iuliana C. Nițu<sup>1,1*</sup> Michael J. Keith<sup>1,1*</sup> David J. Champion<sup>1,2</sup> Ismaël Cognard<sup>1,3,4</sup>

Gregory Desvignes<sup>1,2</sup> Lucas Guillemot<sup>1,3,4</sup> Yanjun Guo<sup>2</sup> Huanchen Hu<sup>1,2</sup> Jiwoong Jang<sup>1,2</sup>

Jedrzej Jawor<sup>1,2</sup> Ramesh Karuppusamy<sup>1,2</sup> Evan F. Keane<sup>1,5</sup> Michael Kramer<sup>2</sup> Kristen Lackeos<sup>1,2</sup>

Kuo Liu<sup>1,2,6</sup> Robert A. Main<sup>2</sup> Delphine Perrodin<sup>1,7</sup> Nataliya K. Porayko<sup>2</sup> Golam M. Shaifullah<sup>1,7,8,9</sup>
and Gilles Theureau<sup>1,3,4,10</sup>

<sup>1</sup>Department of Physics and Astronomy, Jodrell Bank Centre for Astrophysics, The University of Manchester, Manchester M13 9PL, UK

<sup>2</sup>Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>3</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E UMR7328, Université d'Orléans, CNRS, F-45071 Orléans, France

<sup>4</sup>Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, CNRS, Université d'Orléans, F-18330 Nançay, France

<sup>5</sup>School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 PN40, Ireland

<sup>6</sup>Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China

<sup>7</sup>INAF – Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy

<sup>8</sup>Dipartimento di Fisica ‘G. Occhialini’, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

<sup>9</sup>INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

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

Accepted 2024 September 12. Received 2024 August 18; in original form 2024 March 27

## ABSTRACT

In this work, we investigated the presence of strictly periodic, as well as quasi-periodic signals, in the timing of the 25 millisecond pulsars (MSPs) from the European Pulsar Timing Array Second Data Release data set. This is especially interesting in the context of the recent hints of a gravitational wave background in these data, and the necessary further study of red noise timing processes, which are known to behave quasi-periodically in some normal pulsars. We used Bayesian timing models developed through the RUN ENTERPRISE pipeline: a strict periodicity was modelled as the influence of a planetary companion on the pulsar, while a quasi-periodicity was represented as a Fourier-domain Gaussian process. We found that neither model would clearly improve the timing models of the 25 MSPs in this data set. This implies that noise and parameter estimates are unlikely to be biased by the presence of a (quasi-)periodicity in the timing data. Nevertheless, the results for PSRs J1744–1134 and J1012+5307 suggest that the standard noise models for these pulsars may not be sufficient. We also measure upper limits for the projected masses of planetary companions around each of the 25 pulsars. The data of PSR J1909–3744 yielded the best mass limits, such that we constrained the 95 percentile to ~2 × 10⁻⁴ M<sub>⊕</sub> (roughly the mass of the dwarf planet Ceres) for orbital periods between 5 d and 17 yr. These are the best pulsar planet-mass limits to date.

**Key words:** methods: data analysis – planets and satellites: detection – pulsars: general.

## 1 INTRODUCTION

Pulsar timing relies on unambiguously counting every rotation of a pulsar. The observed periodic pulsar emission is considered to be representative of its rotation, which is also assumed to be stable within the time-scale of each observing epoch (e.g. Liu et al. 2012). In practice, pulsar timing involves using a physical model of the pulsar and of the propagating medium to predict the rotational phase at each measured time of arrival (ToA hereafter), and comparing this to the observed phase at some carefully defined fiducial point in the pulsar’s rotation. The difference between the observed and the modelled phase

gives the *timing residual*, often referred to simply as the *residual*. If the model is complete, the residuals should be characterized by white noise with zero mean. However, excess noise is observed in many pulsar residuals in the form of long-term, time-correlated variations (Cordes & Downs 1985; Hobbs, Lyne & Kramer 2010; Parthasarathy et al. 2019), indicating that there may be additional stochastic effects that need to be included in the timing model. These can be achromatic (independent of observing frequency) spin red noise, which has been found to sometimes look quasi-periodic (QP) in slow pulsars, and shows correlations with the observed pulse shapes in the radio emission (Lyne et al. 2010). Moreover, stochastic variations in the interstellar medium may cause chromatic red noise in the timing residuals.

Nevertheless, the high stability of pulses, particularly from the recycled millisecond pulsars (MSPs; e.g. Backer et al. 1982; Verbiest

<footnote>* E-mail: iuliana.nitu29@gmail.com (ICN); michael.keith@manchester.ac.uk (MJK)</footnote>
et al. 2009), makes pulsar timing an invaluable tool for studying a range of astrophysical and cosmological phenomena, such as the structure of the interstellar medium or the solar wind, nuclear matter in exotic conditions, tests of General Relativity (e.g. Manchester 2017), and searches for dark matter and gravitational waves (GWs) from the early Universe (e.g. EPTA Collaboration 2024). Pulsar Timing Arrays (PTAs) represent a network of precisely and frequently timed pulsars, distributed across the Galaxy, with the purpose of detecting GWs in correlated pulsar timing signals. PTAs are expected to be sensitive to space-time distortions caused by GWs of nanohertz frequencies (Sazhin 1978; Detweiler 1979). In this frequency range, the main contribution is believed to come from a stochastic GW background (GWB) created by the incoherent superposition of signals from in-spiralling supermassive black hole binaries (SMBHBs; Rajagopal & Romani 1995). In addition, continuous GWs from strong individual SMBHB sources may also be seen in this regime (Estabrook & Wahlquist 1975), as well as other, more exotic theoretically predicted GW sources such as cosmic strings or a relic GWB (Kibble 1976; Grishchuk 2005). Several collaborations have been developed with the specific purpose of detecting these GWs, including the European Pulsar Timing Array (EPTA; Janssen et al. 2008; Kramer & Champion 2013), the Parkes Pulsar Timing Array (PPTA; Hobbs et al. 2010; Manchester et al. 2013), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; Demorest et al. 2013; NANOGrav Collaboration 2015), the Indian Pulsar Timing Array (InPTA; Joshi et al. 2018), the Chinese Pulsar Timing Array (CPTA; Lee 2016), and the MeerTime Pulsar Timing Array (Spiewak et al. 2022). These groups also work together as part of the International Pulsar Timing Array (IPTA; Hobbs et al. 2010; Manchester et al. 2013; Verbiest et al. 2016).

Recently, results on the latest search for a GWB signature were published concurrently by the EPTA + InPTA (EPTA Collaboration 2023c), PPTA (Reardon et al. 2023), NANOGrav (Agazie et al. 2023), and CPTA (Xu et al. 2023). These reported an emerging evidence of a stochastic GWB in their data sets, ranged between $2\sigma$ and $5\sigma$ in significance, depending on the PTA. A comparison and initial effort to combine these analyses can be found in The International Pulsar Timing Array Collaboration (2024). However, these analyses do not yet meet the requirements for a clear GWB detection. Furthermore, there are still effects in the data of all PTAs that are not yet fully understood, such as how to best model red noise variations in the timing residuals (e.g. EPTA Collaboration 2023b).

In this work, we study the individual-pulsar timing models of the EPTA Second Data Release (DR2), as used in the GWB search (EPTA Collaboration 2023a,b,c). To perform the analyses, we use the Bayesian toolkit `RUN fatENTERPRISE` (Keith, Nițu & Liu 2022) based on the pulsar timing frameworks `ENTERPRISE` (Ellis et al. 2019) and `TEMPO2` (Edwards, Hobbs & Manchester 2006). Specifically, we investigate the consequences on the timing of each of the 25 MSPs of adding a strictly periodic, planet-like component, as well as a Fourier-domain QP component to the Bayesian model fitting. The latter is particularly interesting since, as shown in Keith & Nițu (2023) for slow pulsars, unmodelled QP behaviours present in pulsar timing data can affect the robustness of some parameter and noise estimates.

Furthermore, the sensitivity of this data set provides excellent constraints on the masses of any putative pulsar planetary companions. Although there have been several systematic searches in the timing of pulsars (Thorsett & Phillips 1992; Kerr et al. 2015; Behrens et al. 2020; Nițu et al. 2022), only six pulsars have been confirmed to host planetary-mass companions. The most famous, PSR B1257+12, has three low-mass companions, of 0.020(2), 4.3(2), and 3.9(2) $M_⊕$,

with orbital periods of 25.262(3), 66.5419(1), and 98.2114(2) d, respectively (Wolszczan & Frail 1992; Wolszczan 1994; Konacki & Wolszczan 2003). PSR B1620-26 is in a triple system located in the globular cluster M4, also containing a white dwarf and a ~2.5 Jupiter-mass planet (~800 $M_⊕$) of orbital period ~36 500 d (Thorsett et al. 1999). PSRs J1719-1438 (Bailes et al. 2011), J0636+5128 (Stovall et al. 2014), J1311-3430 (Pletsch et al. 2012; Romani et al. 2012), and J2322-2650 (Spiewak et al. 2018) each have one ‘diamond planet’<fnref n="1" /> companion of mass between ~1–10 Jupiter-mass (~300–3200 $M_⊕$) and orbital period < 1 d. The rarity of pulsar planetary companions is likely a consequence of the extreme conditions in which pulsars form. There is currently no clear mechanism(s) to creating these systems – for an overview of proposed scenarios and formation paths, see for example Podsiadlowski (1993), Phillips & Thorsett (1994) and the Introduction of Nițu et al. (2022). Recently, Nițu et al. (2022) performed the largest scale search for pulsar planetary companions using the timing data sets of 800 pulsars observed at the Jodrell Bank Observatory (JBO). They estimated that at most 0.5 per cent of all pulsars are expected to host Earth-mass planets, concluding that planets around pulsars must be extremely rare. Separately, Behrens et al. (2020) conducted a search on the NANOGrav 11-yr data set, and found no evidence of planetary companions around any of the 45 MSPs. They also estimated the planet-mass sensitivity of the NANOGrav data set, as a function of orbital period. In this analysis, we follow the analysis method in Nițu et al. (2022) to search for and assess the sensitivity of the EPTA DR2 data set to the influence of planetary companions. We then also directly compare our results with those of Behrens et al. (2020).

This paper is structured as follows. In Section 2, the properties of the EPTA DR2 data set are summarized. Section 3 describes the main properties of the timing model used throughout this work. In Section 4, we outline the planet-fitting method, and present and discuss the corresponding results. Section 5 presents the set-up and discusses the results of the QP fitting. In Section 6, we summarize our conclusions.

## 2 DATA SET

The EPTA DR2 is one of the current state-of-the-art pulsar timing data sets, having been recently used in the search for a GWB signature (EPTA Collaboration 2023c). It contains high-precision pulsar timing data from 25 MSPs, collected over 25 yr, with five large telescopes in Europe: the 100-m Effelsberg Telescope (in Germany), Jodrell Bank Observatory’s 76-m Lovell Telescope (in the United Kingdom), Nançay Radio Observatory’s large Radio Telescope (NRT; in France), the Astronomical Observatory of Cagliari’s 64-m Sardinia Radio Telescope (SRT; in Italy), and the Westerbork Synthesis Radio Telescope (WSRT; in the Netherlands). Once a month, these telescopes also functioned collectively, as the Large European Array for Pulsars (LEAP), which is equivalent to a 194-m sixth interferometric telescope in the EPTA (Bassa et al. 2016). Most of the observations part of the EPTA DR2 are at ‘L-band’ frequencies (1–2 GHz) and above, with bandwidths of up to 512 MHz. A limited number of observations are centred at lower frequencies of 350 MHz. For a detailed description of the properties of the EPTA DR2 data set, see Chen et al. (2021) and EPTA Collaboration (2023a).

<footnote n="1">¹The so-called diamond planets are ultra-low mass carbon white dwarfs, believed to be the remains of a disrupted stellar companion (e.g. Bailes, Lyne & Shemar 1991).</footnote>
## 3 CORE TIMING MODEL

Throughout this work, we use the RUNNVENTERPRISE pipeline (Keith et al. 2022) on the timing data of each pulsar. In this Bayesian framework, all the components of the desired model are fit for simultaneously. In the following sections, we discuss, separately, two additional model components which have not been considered in the standard GWB analyses: a simple periodicity (planet) component (Section 4); and a QP Gaussian process component (Section 5). In both cases, however, the analysis includes the same typical ‘core’ (base) timing model, which is summarized in this section.

The core timing model used here follows the same principles as in the individual pulsar analysis of the EPTA DR2 data set (EPTA Collaboration 2023b). The deterministic components of the pulsar model (i.e. the effects of spin frequency and its derivatives, pulsar position, known binary companions, etc.) are fit for simultaneously with the noise components. The parameters describing the former are generally marginalized over for efficiency. The noise model consists of a white noise component as well as between one and three Fourier-domain red noise terms. The white noise is modelled using the EFAC and EQUAD parameters which scale and add in quadrature to the estimated measurement error (see e.g. Verbiest et al. 2016 for more details); these account for unmodelled instrumental errors and intrinsic pulse jitter in the arrival times (Liu et al. 2012; Parthasarathy et al. 2021). The EPTA DR2 data set is composed of primarily narrow-band observations and does not use the ECORR parameter (van Haasteren & Vallisneri 2014). The choice of red noise model makes use of the model selection process in EPTA Collaboration (2023b). The red noise components are chosen from a combination of: (i) an achromatic red noise term, (ii) a chromatic Dispersion Measure (DM) term inducing a timing delay with an $f_{\text{obs}}^{-2}$ dependence on the observing frequency $f_{\text{obs}}$, and (iii) a chromatic scattering variation (SV) term, with an $f_{\text{obs}}^{-4}$ dependence. These are all implemented following the method in Lentati et al. (2014) and using power-law priors for the corresponding power spectral density (PSD), i.e.

$$P(f) = \frac{A^2}{K} \left( \frac{f}{1 \text{ yr}^{-1}} \right)^{-\gamma}, \qquad (1)$$

where f is the Fourier frequency, A and γ are the power-law amplitude and slope (index), respectively, and K is a scale factor for each process. In total there are six hyperparameters: the slopes for each process, γ<sub>red</sub>, γ<sub>DM</sub>, γ<sub>SV</sub>; and the corresponding log-amplitudes, log<sub>10</sub>A<sub>red</sub>, log<sub>10</sub>A<sub>DM</sub>, log<sub>10</sub>A<sub>SV</sub>. For the achromatic and SV processes, K = 12π<sup>2</sup>, and for the DM process K = k<sub>DM</sub><sup>2</sup>, with the DM constant k<sub>DM</sub> = 2.41 × 10<sup>-4</sup>cm<sup>-3</sup>pc MHz<sup>2</sup>s<sup>-1</sup>. For more details see EPTA Collaboration (2023b). The Fourier basis has equally spaced frequencies $f_n = n/T_{\text{span}}$, with n ∈ {1, 2, ..., N<sub>c</sub>}, T<sub>span</sub> the total time span of the data, and N<sub>c</sub> the number of Fourier components. In this study, we set N<sub>c</sub> for each of the achromatic, DM, and scattering red noise to the ‘optimal’ values for this data set as found by EPTA Collaboration (2023b).

## 4 FITTING FOR A PLANET INFLUENCE

### 4.1 Set-up

We use the Bayesian planet fitting model implemented in RUNNVENTERPRISE, which is introduced and described in Nițu et al. (2022). In short, this is parametrized by the (fitted) orbital parameters:

projected mass<fnref n="2" /> of the planet, m sin i (where m is the planetary mass and i is the inclination of the orbit), orbital period, P<sub>b</sub>, eccentricity, e, argument of periapsis, ω, and the phase φ of the planet on the orbit with respect to the periastron crossing, defined at a reference time t<sub>ref</sub> = 55000 MJD. These parameters all contribute to the Rømer delay, which quantifies the planetary influence on the pulsar signal for a ToA t. This can be expressed as (e.g. Blandford & Teukolsky 1976)

$$\Delta_R(t) = C P_b^{2/3} m \sin i \left[ (\cos E(t) - e) \sin \omega + \sin E(t) \sqrt{1 - e^2 \cos \omega} \right], \qquad (2)$$

where C is a constant of proportionality. We set the mass of each pulsar to be fixed at a characteristic value of M<sub>PSR</sub> = 1.4 M<sub>☉</sub> (e.g. Lattimer 2012), and to be much larger than the mass of the planet, such that C ≃ 23.4 d<sup>-2/3</sup>M<sub>☉</sub><sup>-1</sup> ms with the Rømer delay generally expressed in microseconds, P<sub>b</sub> in days, and m in Earth masses (M<sub>☉</sub>). The eccentricity anomaly E(t) in equation (2) is related to the true (A<sub>T</sub>(t)) and mean (M(t)) anomalies by

$$\cos E(t) = \frac{e + \cos A_{\text{T}}(t)}{1 + e \cos A_{\text{T}}(t)} \qquad (3)$$

and

$$M(t) \equiv \frac{2\pi}{P_b} (t - t_0) = E(t) - e \sin E(t), \qquad (4)$$

where t<sub>0</sub> is the time of closest periastron approach. To find an equation for E(t) that is only a function of t and the fitted orbital parameters (in this case e and φ), we start from the definition of φ, i.e. A<sub>T</sub>(t<sub>ref</sub>) = 2πφ. It follows that E(t<sub>ref</sub>) can be obtained from φ, e, and equation (3) at t ≡ t<sub>ref</sub>. Further, from equation (4) at t ≡ t<sub>ref</sub>, we can compute

$$t_0 = t_{\text{ref}} - \frac{P_b}{2\pi} [E(t_{\text{ref}}) - e \sin E(t_{\text{ref}})]. \qquad (5)$$

Finally, E(t) can consequently be estimated for any t by solving the non-linear expression in equation (4).

The Bayesian priors of the five orbital parameters determining the planet influence (m sin i, P<sub>b</sub>, e, ω, and φ) are set-up as follows. Uniform priors are used for ω ∈ [0, 2π) and φ ∈ [0, 1). Log-uniform priors are used to explore the parameter space for e ∈ [0, 0.9], and the projected mass m sin i ∈ [10<sup>-5</sup>, 10<sup>-1</sup>] M<sub>☉</sub>. Note that we search for lower planetary masses than in the Nițu et al. (2022) analysis since we expect the EPTA DR2 data set to be more sensitive to such influences than the Jodrell Bank Observatory data set by itself. We split the parameter space of P<sub>b</sub> into 11 period bins, to thoroughly explore the large prior range; the bounds of these are (5.0, 10.1, 21.3, 42.5, 85, 170, 340, 390, 780, 1560, 3120, 6240) d. 10 of these bins are log-uniformly spaced, while one narrower bin (340–390 d) is considered around P<sub>b</sub> = 1 yr to account for the loss of sensitivity due to also fitting for the pulsar position and parallax (which have a 1-yr periodicity) in the same analysis.

The EPTA DR2 timing data for each of the 25 MSPs are initially processed through the pipeline RUNNVENTERPRISE, including the ‘core’ and planet model, using the PYTHON Markov chain Monte Carlo (MCMC) sampler EMCEE (Foreman-Mackey et al. 2013, 2019). From the posteriors, planetary-mass limits are then estimated at the 95 per cent threshold for each period bin in each pulsar. Potential planet candidates are selected based on a 3σ ‘flag’ in the linear mass

<footnote n="2">²Note that whenever we talk about a planetary-companion ‘mass’ in this work, we mean the projected mass, as it is not possible to disentangle the inclination dependence in this type of analysis without additional information.</footnote>
<table><caption>**Table 1.** A summary of the main properties of the periodicity analysis on the four MSPs showing $3\sigma$ detections in the planet-mass posterior. The second column shows the natural-logarithm Bayes factor (evidence) in favour of including the planet model to the 'core' timing model. The mean and standard deviation values are given for the planetary projected mass, $m \sin i$, and orbital period, $P_b$. The uncertainties are given in standard parenthetical notation, i.e. representing one standard deviation in the last digit. The inverse of the maximum-likelihood fundamental frequency of the QP analysis ($f_{qp}$) is given for comparison ('maxL'); the standard deviation ('stdev') of the $1/f_{qp}$ posterior is also included, but should only be used to give an idea of the spread of values; see also Fig. 7 for the shape of these posteriors. The total time of the observations, $T_{span}$, is included for a comparison with the quoted periodicities. Finally, the choice of noise models (achromatic 'Red' noise, 'DM' Noise, and scattering variation 'SV') used in the analysis of each pulsar is shown, as determined by EPTA Collaboration (2023b).</caption><thead><tr><th rowspan="2">PSR</th><th rowspan="2">lnB</th><th colspan="2">m sin i (10<sup>-4</sup> M<sub>☉</sub>)</th><th colspan="2">P<sub>b</sub>(d)</th><th colspan="2">1/f<sub>qp</sub>(d)</th><th rowspan="2">T<sub>span</sub> [d]</th><th colspan="3">Noise models</th></tr><tr><th>EMCEE | DYNESTY</th><th>DYNESTY</th><th>EMCEE | DYNESTY</th><th>DYNESTY</th><th>maxL</th><th>stdev</th><th>Red</th><th>DM</th><th>SV</th></tr></thead><tbody><tr><td>J0751+1807</td><td>3.1(8)</td><td>1.5(6) | 1.6(6)</td><td></td><td>2490(230)</td><td>2460(210)</td><td>1940</td><td>630</td><td>8835</td><td>✗</td><td>✓</td><td>✗</td></tr><tr><td>J1012+5307</td><td>0.5(9)</td><td>1.6(7) | 1.0(9)</td><td></td><td>760(30)</td><td>730(130)</td><td>770</td><td>560</td><td>8648</td><td>✓</td><td>✓</td><td>✗</td></tr><tr><td>J1744-1134</td><td>0.4(8)</td><td>0.6(2) | 0.4(3)</td><td></td><td>1550(70)</td><td>1480(200)</td><td>1460</td><td>500</td><td>8770</td><td>✓</td><td>✓</td><td>✗</td></tr><tr><td>J1918-0642</td><td>1.7(7)</td><td>1.5(1.4) | 2.4(1.3)</td><td></td><td>3000(570)</td><td>2940(330)</td><td>720</td><td>480</td><td>7199</td><td>✗</td><td>✓</td><td>✗</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The timing residuals, without removing the possible planet influence, of the four MSPs that showed a $3\sigma$ flag in the planet fitting (black data points), and the corresponding maximum-likelihood Rømer delay (green continuous curve), for comparison. Note that for PSR J1744-1134 only, the green curve also includes the lowest two Fourier-frequency components of the achromatic red noise to aid the direct comparison. See text for more details on how these residuals have been computed.</figcaption></figure>

posterior, i.e. when the mean of the distribution is more than three times the standard deviation. For more details on these processes, see Nițu et al. (2022). Moreover, for each of these ‘flagged’ pulsars, a further nested-sampling run is performed, using the PYTHON sampler DYNESTY (Speagle 2020). This nested sampling method yields the Bayesian evidence in favour of including the planet model, as well as providing a consistency check for the parameter values obtained

from the MCMC run. Further investigations are then made into the credibility of the potential planet detection, by inspecting the residuals, and the PSD plots in this context.

## 4.2 Results and discussion

### 4.2.1 The $3\sigma$ flagged MSPs

For 4 of the 25 MSPs in this study, the planet-mass posterior distribution in one period bin has the property that the mean is more than three standard deviations away from zero. We call these the ‘flagged’ MSPs, or those showing ‘$3\sigma$ detections’. Three of these MSPs are known to be in stellar binary systems, with orbital periods: 0.26 d for PSR J0751+1807, 0.60 d for PSR J1012+5307, and 10.91 d for PSR J1918-0642. The former two are below the minimum orbital period that we search for in this analysis (5 d) and are therefore unlikely to affect this substantially. The latter may reduce the sensitivity of the lowest period bins we use. The fourth ‘flagged’ MSP, PSR J1744-1134, is a solitary pulsar.

Table 1 shows the projected mass and orbital period of the potential planet companion for each pulsar, as estimated from the posteriors of the separate MCMC (EMCEE) and nested sampling (DYNESTY) runs, respectively. The eccentricity is not quoted as the posterior distributions generally recover the prior for all pulsars. The estimated log-Bayes evidence in favour of including a planet influence in the timing model for each pulsar is also shown in Table 1. As all log-Bayes factors are within four standard deviations of 0, these values suggest that there is no support for the planet model from the Bayes factors.

We note that the masses shown in Table 1 are roughly two orders of magnitude below even the smallest known pulsar planet, i.e. the 0.02-M<sub>☉</sub> companion of PSR B1257+12 (Wolszczan & Frail 1992; Wolszczan 1994; Konacki & Wolszczan 2003). On the other hand, the orbital periods are two orders of magnitude above those of the planets of PSR B1257+12 (<100 d), such that the inferred pulsar-planet distances would be similarly larger in these systems. This is somewhat expected, as our data set is highly sensitive to these kinds of behaviours. Interestingly, the 800-M<sub>☉</sub> planet in the triple-system B1620-26 has an orbital period roughly an order of magnitude above those of our ‘flagged’ pulsars, three of which are also known to be in binary systems (Thorsett et al. 1999).

To further understand these potential planet detections and their cause, we look at how the fitted planet influence compares to the residuals of each of the four pulsars (Fig. 1), as well as the shape of the PSD of the achromatic red noise modelled in these residuals in the absence of the planet fitting (Fig. 2).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The estimated PSD as a function of Fourier frequency (black continuous line); the best-fitting power-law red noise (diagonal red dotted line), or white noise level (horizontal blue dot-dashed line) is shown, as well as the orbital period detected (vertical green line), all for each of the four ‘flagged’ MSPs. The grey dashed vertical line represents the 1/yr frequency.</figcaption></figure>

The choice of noise models used for each of the four MSPs is summarized in Table 1, as per the modelling selected in EPTA Collaboration (2023b). Correspondingly, the residuals are computed by subtracting from the observed ToAs the modelled DM time series and the best deterministic timing model, except for the planet influence. This is enough to show the left-over planet influence. Note that the modelled noise that is removed is the best-fitting model to the data in the simultaneous Bayesian fit which included the planet influence. Thus the resulting residuals should only contain white noise, any achromatic red noise, and the planet influence, according to the best timing model; these residuals can therefore be compared to the analytical Rømer delay as computed from the best-fitting orbital parameters. Note that for PSR J1744−1134 – the only one of the four to have significant achromatic red noise – we compare the residuals obtained as above with the analytical Rømer delay to which we also add the lowest two Fourier-frequency achromatic red noise components. These comparisons are shown in Fig. 1.

In these plots, perhaps the most noticeable property is that the
possible planet influence is generally of small amplitude, similar in
value to the white noise level, and therefore largely indistinguishable
from it. This is not unexpected following the mass and orbital period
values quoted in Table 1, as a companion of mass ~10⁻⁴ M⊕ orbiting
with a period of order 1 yr induces a maximum variation of order
0.1 µs in the pulsar ToAs (estimated from the Rømer delay as
in equation 2). For the two pulsars of highest Bayesian evidence

(J0751+1807 and J1918−0642) the detected periodicity is large
compared to the total time span of the data, such that less than four
full periodicities are seen in the analysed data set. It is therefore not
clear whether this observed variation will continue to behave as a
simple periodicity, or will take the form of a red noise like process
in a longer data set.

It is also interesting to investigate the shape of the fitted achromatic
red noise in the absence of the planet fitting model, particularly in the
Fourier domain. We expect any noticeable periodicity in the residuals
to show as a peak in the corresponding PSD of these residuals,
which may also bias a simple power-law fit. To estimate these PSD
shapes from the residuals, we employ the widely used ‘Cholesky
method’ as described in Coles et al. (2011), and also adopted in
CHOLSPECTRA. In short, this is based on estimating the covariance
matrix of the residuals assuming the power-law form of the PSD
(as in equation 1), then using the Cholesky decomposition on this
covariance matrix to determine the transformation that ‘whitens’ the
residuals. Thus the fitting problem becomes a simple ordinary least-
squares on uncorrelated data, and the PSD can be straightforwardly
estimated from the best-fitting parameters. The appropriate residuals
are obtained by subtracting the best-fitting deterministic timing
model, as well as any fitted DM and SV noise, from the pulsar
signal; note that this is now in the case of a model *not* including a
planet influence. These PSD estimates are seen in Fig. 2, together
with the best power-law model (or lack there-of), and the maximum-
likelihood $P_b$ from the planet-fitting runs.

Similarly to the findings from the residual plots of these pulsars,
any tentative peak in the PSD at the 'detection' orbital period cannot
be unambiguously distinguished from a stochastic variation around
the power-law prior. While for PSRs J0751+1807 and J1918-0642
the best noise model as per EPTA Collaboration (2023b) did not
include an achromatic term, the estimated PSD appears to have some
excess power at very low Fourier frequencies compared to a flat
('white') power. For both of these pulsars, comparing the DM-only
model to a combination of DM and achromatic noise favours the
combined model, but with an insignificant Bayes factor of $\ln B \sim 1$.
The planet-like periodicities found for these two pulsars are similar
to the time span of their data, and the evidence for the planetary
model is similar to that for an achromatic power-law noise model.
Therefore, we argue that the simplest conclusion is that these pulsars
have a low level of achromatic noise that was not included in the base
model, which can equally well be modelled by a single low-frequency
sinusoid or a power-law process. If we repeat the analysis with the
addition of an achromatic power-law model, there is no evidence in
favour of a planetary companion. Given the prevalence of power-law
red noise in pulsars, and the scarcity of planetary companions, it
seems most plausible that this is simply unmodelled red noise.

To summarize, given the discussed properties and mass limits
obtained, there is no evidence of planetary companions in this data
set, and the planet-like periodicities found in this analysis are likely
an artefact of the choice of noise models. However, the relatively large
achromatic red noise power on time-scales close to the dataspan is
worth investigating further in future work. Interestingly, the PSRs
J1012+5307 and J1744-1134 highlighted in this analysis are also
characterized as having some level of complex behaviours in the
recent noise analysis of the EPTA DR2 data set (EPTA Collaboration
2023b). Indeed, we notice in the PSD shapes of these pulsars that
there is additional complexity with respect to a single power-law
model. EPTA Collaboration (2023b) propose that these effects may
be due to, for example, a non-stationarity of the stochastic red noise,
or perhaps some unknown instrumental effects. They also suggest
that the observed noise properties in these pulsars should be further
<figure><img src="image_4.png" /></figure>

**Figure 3.** The 95 per cent planet-mass limits for each orbital period bin, as estimated from the planet fitting results for the four ‘flagged’ MSPs (black step line). The properties in the top right are the spin-period of the pulsar and its derivative, as well as the orbital period of the known binary companion. The black squares represent the maximum-likelihood values in each period bin, while the stars represent $3\sigma$ detections. The additional runs at particular period bins were informed by the initial results, and are shown in horizontal grey lines for PSRs J1012+5307 and J1744–1134. The vertical green dashed lines show the fundamental periodicity and harmonics of the known stellar companions of the respective pulsars. The pink line shows the sensitivity from the NANOGrav 11-yr data set, as estimated by Behrens et al. (2020). In PSR J1744–1134, the blue line represents the mass limits as estimated only from the JBO data, as given in Nițu et al. (2022).

studied with a different data set such as the upcoming IPTA data combination, in an attempt to fully understand their behaviours.

### 4.2.2 Companion mass limits

It is also interesting to consider the planet-mass upper limits that can be inferred from the planet search, as we expect this data set to be highly sensitive to such influences. For each of the 11 period bins that we perform our analysis on, we can estimate the 95 per cent upper limit on the projected mass, simply from the posterior distribution of the mass. Note that the sampling gives the log-mass posterior, which is then converted to the linear-mass posterior for this estimate.

Fig. 3 shows the estimated 95 per cent mass limits for the four ‘flagged’ MSPs, while Fig. 4 shows the mass limits of the most sensitive pulsar in this data set, PSR J1909–3744. Fig. 5 shows the results for the other 20 MSPs. Overall, the mass limits shown here

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

**Figure 4.** The 95 per cent planet-mass limits for each orbital period bin, as estimated from the planet fitting results for the most-sensitive MSPs in this data set, PSR J1909–3744. The pink continuous line shows the sensitivity from the NANOGrav 11-yr data set, as estimated by Behrens et al. (2020).

are remarkably low. For example, any possible planet companion around any of the 25 MSPs is highly unlikely to have a projected mass higher than $10^{-3} M_\oplus$, for orbital periods between roughly 20 d and 17 yr. This is excluding the small area of parameter space where $P_b \approx 1$ yr, where sensitivity is lost due to fitting for the proper motion and distance to the pulsar; we further ignore this in the following discussion.

PSR J1744–1134 – which is a solitary MSP – was also part of the planet search on JBO data in Nițu et al. (2022). We can therefore directly compare the sensitivity of the JBO data by itself and the EPTA DR2, and this is shown in Fig. 3 for this pulsar; note that the JBO data are included in the EPTA DR2. As can be seen, for PSR J1744–1134 the EPTA timing data are more than an order of magnitude more sensitive to a planet influence than the JBO data by itself, while showing similar trends with different orbital periods.

The best mass limits found in our analysis are those estimated from the data of PSR J1909–3744 (shown in Fig. 4), which constrain the 95 percentile of a planet mass to be lower than $2 \times 10^{-4} M_\oplus$ for all orbital periods investigated (5 d and 17 yr). For context, this value is approximately equal to the mass of the dwarf planet Ceres ($1.6 \times 10^{-4} M_\oplus$; Park et al. 2019), and roughly a tenth of the mass of Pluto ($2.2 \times 10^{-3} M_\oplus$; Stern et al. 2015). Furthermore, this is roughly an order of magnitude better than previously published state-of-the-art sensitivity limits, also plotted in Fig. 3 for comparison. This sensitivity curve was derived by Behrens et al. (2020) using simulations, for a ‘typical’ MSP in the NANOGrav 11-yr data set. Roughly half of the 25 MSPs show mass limits better than those computed by Behrens et al. (2020), while the rest are still within one order of magnitude of these.

We note that, although we have not specifically searched for a population of ‘diamond’ planets – since they would be at orbital periods of <1 d (Spiewak et al. 2018), and this is not robust with our observational cadence – this can still be ruled out in our data set. This is because, as the residuals have noise of about 10 µs, any planet-like influence corresponding to an orbital period <1 d would be obvious above the noise for a projected mass > 0.4 $M_\oplus$. A diamond-planet is expected to be above a Jupiter-mass (i.e. roughly 300 $M_\oplus$), which would therefore make it immediately obvious in these residuals.

## 5 FITTING FOR A QUASI-PERIODIC GAUSSIAN PROCESS

### 5.1 Set-up

It is well established that QP timing noise is prevalent in the ‘normal’ pulsar population (Hobbs et al. 2010), thought to be due to multimodal switching of magnetospheric processes (Lyne et al.
<figure><img src="image_6.png" /><figcaption>Figure 5. 95 per cent mass limits for the 20 not-‘flagged’ MSPs. See the caption of Fig. 3 for more details.</figcaption></figure>
2010). While power-law models can perform well in such conditions, strong QP noise can lead to errors in parameter estimation and overestimation of the power-law noise model (Keith & Nițu 2023). If timing noise in MSPs is governed by similar QP processes then there is potential to better characterize the noise and consequently increase sensitivity to GW signals at the lowest frequencies. Therefore, we search for a QP process in the timing data of each of the 25 MSPs in the EPTA DR2 data set. The same ‘core’ timing model is used, supplemented with the Fourier-domain Gaussian-process QP model described in Keith & Nițu (2023), and implemented in RUNNVENTERPRISE, which is outlined in the following.

The full timing model is fit simultaneously using a Bayesian method and the MCMC sampler EMCEE, as before. The shape of the QP model used was motivated by the observed frequency-domain behaviours of slow pulsars (see e.g. Lyne et al. 2010). The QP effect is represented in Fourier-domain using the same infrastructure as the power-law red noise models. The characteristic PSD of the spin-frequency derivative of the pulsar, $\dot{v}$, is described by a sum of $N_{\text{harm}}$ Gaussian-function terms. Each Gaussian function is centred at harmonically related frequencies $f_k = k f_{\text{qp}}$, with $k \in \{1, \dots, N_{\text{harm}}\}$, where $f_{\text{qp}} \equiv f_1$ is the fundamental frequency of the quasi-periodicity. The Gaussian functions have increasing widths, characterized by a standard deviation $\sigma_k = f_k \sigma = k f_{\text{qp}} \sigma$, where $\sigma$ is a dimensionless ‘fractional’ standard deviation, and is the same for all harmonically related terms. Furthermore, the amplitude of each Gaussian function is also decreasing exponentially, by a factor of $\exp[-(k-1)/\lambda]$, where in practice the value of $\lambda$ quantifies the number of significant harmonics.

Consequently, the PSD of the residuals for the QP process is described as

$$P_{\text{qp}}(f) = \begin{cases} R_{\text{qp}} P_{\text{pl}}(f_{\text{qp}}) q(f) \left(\frac{f}{f_{\text{qp}}}\right)^{-4} & , f \ge f_{\text{cut}} \\ 0 & , f < f_{\text{cut}}. \end{cases} \quad (6)$$

for a Fourier frequency $f$ and a threshold

$$f_{\text{cut}} = 0.5 f_{\text{qp}} \left(1 - \sqrt{1 - 16\sigma^2}\right) \quad (7)$$

defined by the local minima of $P_{\text{qp}}(f)$ to avoid an unwanted increase at very low frequencies. Further,

$$P_{\text{pl}}(f) = \frac{A_{\text{red}}^2}{12\pi^2} \left(\frac{f}{1\text{yr}^{-1}}\right)^{-\gamma_{\text{red}}} \text{yr}^3 \quad (8)$$

is the power-law PSD for achromatic red noise (equivalent to equation 1), and the quantity

$$q(f) = \sum_{k=1}^{N_{\text{harm}}} \frac{1}{k} \exp\left[-\frac{(k-1)}{\lambda}\right] \exp\left[-\frac{(f-k f_{\text{qp}})^2}{2k^2 f_{\text{qp}}^2 \sigma^2}\right] \quad (9)$$

encapsulates the QP-type variability. In practice, we set $N_{\text{harm}} = 10$ for simplicity, since we expect fewer than 10 harmonics to contribute significantly to the total signal in all cases. The $f^{-4}$ dependence in equation (6) is due to the transformation between the $\dot{v}$ and the residuals, as $\dot{v} \propto -\ddot{v} \propto f^2 r$. For more details on this choice of QP model, see Keith & Nițu (2023). The hyperparameters describing this process are $R_{\text{qp}}$, $f_{\text{qp}}$, $\lambda$, and $\sigma$. We set uniform priors on the log of the ratio, $\log_{10} R_{\text{qp}}$, the central periodicity, $1/f_{\text{qp}}$, and on $\lambda$ and $\sigma$ as defined in Table 2. The upper bound on $1/f_{\text{qp}}$ is chosen to prevent it becoming degenerate with the power-law noise process, and the lower bound is chosen to keep $f_{\text{qp}}$ within the bounds of the Fourier basis.

The QP model is designed to be used alongside the power-law red noise model, as described in Section 3, such that the full achromatic

<table><caption>**Table 2.** Uniform prior bounds for the QP process model. $N_c$ is the number of Fourier coefficients in the achromatic power-law noise model.</caption><thead><tr><th>Parameter</th><th>Lower</th><th>Upper</th></tr></thead><tbody><tr><td>log<sub>10</sub> R<sub>qp</sub></td><td>-2</td><td>3.5</td></tr><tr><td>1/f<sub>qp</sub></td><td>T<sub>span</sub>/4</td><td>1.2T<sub>span</sub>/N<sub>c</sub></td></tr><tr><td>λ</td><td>0.01</td><td>10</td></tr><tr><td>σ</td><td>10<sup>-3</sup></td><td>0.2</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The functional form of the total QP and power-law PSD of the residuals, as well as its two separate components, as given by equations (8) and (6). The parameters used are $A_{\text{pl}} = 3.9 \times 10^{-10}$, $\gamma = 4.3$, $R_{\text{qp}} = 501.2$, $f_{\text{qp}} = 0.39 \text{ yr}^{-1}$, $\sigma = 0.047$, and $\lambda = 0.7$.</figcaption></figure>

red noise is modelled according to a PSD of $P_{\text{qp}} + P_{\text{pl}}$. Fig. 6 illustrates the shape of the functional form of $P_{\text{qp}}$ and $P_{\text{pl}}$ for some example parameters.

## 5.2 Results and discussion

In short, our analysis concluded that none of the MSPs are well-characterized by our QP model. In all cases, the posterior distributions of the parameters $\lambda$ and $\sigma$ are almost identical to their respective uniform prior. A similar situation was seen for parameter $R_{\text{qp}}$, which quantifies the strength of the QP power compared to the power-law red noise power. While the prior for $R_{\text{qp}}$ was chosen to be log-uniform, its posterior distribution generally recovers this prior, but with a decrease at large values; the only exception to this is perhaps the case of PSRs J0751+1807 and J1744–1134. This is not unexpected, as it illustrates a (not very constraining) upper limit for any existing QP-type power; the posterior distributions suggest that it is unlikely that $R_{\text{qp}} > 10^3$ for any of the studied MSPs.

The fourth hyperparameter, $f_{\text{qp}}$, describes the fundamental frequency of the QP process. For most MSPs in this analysis, the posterior of $f_{\text{qp}}$ is also unconstraining. However, for PSRs J0751+1807, J1012+5307, and J1744–1134 there was some level of preference for particular frequencies. These pulsars are included in the subset of four MSPs that were ‘flagged’ in the planet search in Section 4.2.1. We therefore concentrate on the results of the QP fitting for the previously considered four MSPs. The most interesting, and most constrained properties of the QP model are the fundamental period of the QP process (given by $1/f_{\text{qp}}$), and the power amplitude at this Fourier frequency. According to equations (6) and (8), we can write this amplitude at $f_{\text{qp}}$ as

$$A_{\text{qp}} = R_{\text{qp}} \frac{A_{\text{red}}^2}{12\pi^2} \left( \frac{f_{\text{qp}}}{1 \text{yr}^{-1}} \right)^{-\gamma_{\text{red}}} \quad (10)$$

Thus, we straightforwardly derive the posterior of $\log_{10} A_{\text{qp}}$ from the posteriors of $\log_{10} R_{\text{qp}}$, $\log_{10} A_{\text{red}}$, $f_{\text{qp}}$, and $\gamma_{\text{red}}$. We note that although
<figure><img src="image_8.png" /><figcaption>**Figure 7.** A subset of the corner-plots for the four pulsars highlighted in the periodicity search of Section 4. The properties of the QP model fitted are summarized through the fundamental periodicity $1/f_{qp}$ and the power amplitude at this corresponding Fourier frequency. The contour levels are at the (0.5, 1, 1.5, 2)σ equivalent, as used in the library CORNER.PY. The purple vertical lines correspond to the maximum-likelihood values of each parameter.</figcaption></figure>

the prior on $\log_{10} R_{qp}$ is uniform, the effective prior on the derived parameter $\log_{10} A_{qp}$ is more complex. Fig. 7 shows the posteriors of $1/f_{qp}$ and $\log_{10} A_{qp}$ and the relationship between them for the four mentioned pulsars. Of the four MSPs, PSR J1012+5307 shows the clearest peak in the $1/f_{qp}$ posterior. The corresponding posterior of PSR J0751+1807 shows a preference for the maximum of the prior, which is set to be a quarter of the total observing span. Further, the QP-period posterior of J1744–1134 is wide, but nearly Gaussian-shaped. On the other hand, the posterior of PSR J1918–0642 does not significantly differ from the prior. None of the amplitude posteriors of the four MSPs are strictly Gaussian shaped, although the results for PSRs J1012+5307 and J1918–0642 are reasonably close. For PSRs J0751+1807 and J1744–1134, the posterior distributions have a long tail at low amplitudes, with an abrupt cut-off near $A_{qp} \sim 10^{-26} \text{ yr}^3$.

Furthermore, Table 1 shows a direct comparison between the period found in the planet fitting, and the inverse of the maximum-likelihood fundamental frequency of this QP fitting. The standard deviation estimated from the posterior of $1/f_{qp}$ is also included to give an idea of the width of the posterior distribution, but this should be considered in conjunction with the shape of the posterior, which is generally not exactly Gaussian (as shown in Fig. 7). From the results in Table 1, we can see that in all but PSR J1918–0642 the maximum-likelihood QP fundamental periodicity is similar to the planetary orbital period. In the case of PSRs J0751+1807 and J1744–1134, the large periodicity as a fraction of the total data span suggests that both the planet and the QP fitting are trying to account for the large power observed at high time-scales. For PSRs J1012+5307, the periodicity of approximately 2 yr is found in the QP fitting as well, although not well-described by this QP model.

# 6 CONCLUSIONS

Overall, we conclude that none of the timing models of the 25 MSPs in the EPTA DR2 would clearly benefit from the addition of a QP/periodic process. In particular, the lack of a detectable Fourier-domain QP Gaussian process in these data means we do not expect there to be a bias in parameters or noise estimates due to this kind of behaviour, as was found by Keith & Nițu (2023).

Four pulsars are highlighted as showing potential periodic or QP processes in our analysis. Two (PSRs J0751+1807 and J1918–0642) are pulsars for which the original EPTA analysis selected against achromatic red noise, but we find it likely that the periodic signal we measure is caused by unmodelled red noise. The other two (PSRs J1744–1134 and J1012+5307) are among the pulsars suggested to have ‘complex’ behaviours in the EPTA noise analysis (EPTA Collaboration 2023b). Our analysis also hints that a pure power-law process may not be sufficient, however, neither the periodic nor the QP models we trialled seem to meaningfully improve the results. These pulsars will certainly continue to be studied in the near future with EPTA and IPTA analyses, as we attempt to better specify the noise models for PTAs.

Finally, the planet-fitting analysis on these 25 MSPs allowed us to put highly constraining limits on the masses of any planetary companions orbiting these pulsars. The timing data of PSR J1909–3744, which yielded the best mass limits, allowed us to constrain the 95 percentile to approximately the mass of the dwarf planet Ceres (~2 × $10^{-4} M_\oplus$) for orbital periods between 5 d and 17 yr. These limits are more than an order of magnitude improved compared to the previous sensitivity curve estimated by Behrens et al. (2020) for the NANOGrav 11-yr data, and are the best planet-mass limits from pulsar timing to date.

# ACKNOWLEDGEMENTS

ICN was supported by the STFC doctoral training grant ST/T506291/1.

Part of the EPTA data used in 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 (STFC). The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from ‘Programme National de Cosmologie and Galaxies’ (PNCG), and ‘Programme National Hautes Energies’ (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA, and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). The Sardinia Radio Telescope (SRT) is funded by the Department of University and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as a National Facility by the National Institute for Astrophysics (INAF).

# DATA AVAILABILITY

The EPTA DR2 data underlying the work in this paper are available at https://doi.org/10.5281/zenodo.8164424, as per EPTA Collaboration (2023a).
## REFERENCES

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

Backer D. C., Kulkarni S. R., Heiles C., Davis M. M., Goss W. M., 1982, *Nature*, **300**, 615

Bailes M., Lyne A. G., Shemar S. L., 1991, *Nature*, **352**, 311

Bailes M. et al., 2011, *Science*, **333**, 1717

Bassa C. G. et al., 2016, *MNRAS*, **456**, 2196

Behrens E. A. et al., 2020, *ApJ*, **893**, L8

Blandford R., Teukolsky S. A., 1976, *ApJ*, **205**, 580

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

Coles W., Hobbs G., Champion D. J., Manchester R. N., Verbiest J. P. W., 2011, *MNRAS*, **418**, 561

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

Demorest P. B. et al., 2013, *ApJ*, **762**, 94

Detweiler S., 1979, *ApJ*, **234**, 1100

EPTA Collaboration, 2023a, *A&A*, **678**, A48

EPTA Collaboration, 2023b, *A&A*, **678**, A49

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

EPTA Collaboration, 2024, *A&A*, **685**, A94

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

Ellis J. A., Vallisneri M., Taylor S. R., Baker P. T., 2019, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust Pulsar Inference SuitE, zenodo

Estabrook F. B., Wahlquist H. D., 1975, *Gen. Relativ. Gravit.*, **6**, 439

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, *PASP*, **125**, 306

Foreman-Mackey D. et al., 2019, *J. Open Source Softw.*, **4**, 1864

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

Hobbs G., Lyne A. G., Kramer M., 2010, *MNRAS*, **402**, 1027

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

Joshi B. C. et al., 2018, *JA&A*, **39**, 51

Keith M. J., Nițu I. C., 2023, *MNRAS*, **523**, 4603

Keith M. J., Nițu I. C., Liu Y., 2022, *run_enterprise*, zenodo

Kerr M., Johnston S., Hobbs G., Shannon R. M., 2015, *ApJ*, **809**, L11

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

Konacki M., Wolszczan A., 2003, *ApJ*, **591**, L147

Kramer M., Champion D. J., 2013, *Class. Quantum Gravity*, **30**, 224009

Lattimer J. M., 2012, *Annu. Rev. Nucl. Part. Sci.*, **62**, 485

Lee K. J., 2016, in *Qain L., Li D., eds, ASP Conf. Ser. Vol. 502, Frontiers in Radio Astronomy and FAST Early Sciences Symposium 2015. Astron. Soc. Pac., San Francisco*, p. 19

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

Liu K., Keane E. F., Lee K. J., Kramer M., Cordes J. M., Purver M. B., 2012, *MNRAS*, **420**, 361

Lyne A., Hobbs G., Kramer M., Stairs I., Stappers B., 2010, *Science*, **329**, 408

Manchester R. N., 2017, *J. Phys.: Conf. Ser.*, **932**, 012002

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

NANOGlav Collaboration, 2015, *ApJ*, **813**, 65

Nițu I. C., Keith M. J., Stappers B. W., Lyne A. G., Mickaliger M. B., 2022, *MNRAS*, **512**, 2446

Park R. S. et al., 2019, *Icarus*, **319**, 812

Parthasarathy A. et al., 2019, *MNRAS*, **489**, 3810

Parthasarathy A. et al., 2021, *MNRAS*, **502**, 407

Phillips J. A., Thorsett S. E., 1994, *Ap&SS*, **212**, 91

Pletsch H. J. et al., 2012, *Science*, **338**, 1314

Podsiadlowski P., 1993, in *Phillips J. A., Thorsett S. E., Kulkarni S. R., eds, ASP Conf. Ser. Vol. 36, Planets Around Pulsars. Astron. Soc. Pac., San Francisco*, p. 149

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

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

Romani R. W., Filippenko A. V., Silverman J. M., Cenko S. B., Greiner J., Rau A., Elliott J., Pletsch H. J., 2012, *ApJ*, **760**, L36

Sazhin M. V., 1978, *SvA*, **22**, 36

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

Spiewak R. et al., 2018, *MNRAS*, **475**, 469

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

Stern S. A. et al., 2015, *Science*, **350**, aad1815

Stovall K. et al., 2014, *ApJ*, **791**, 67

The International Pulsar Timing Array Collaboration, 2024, *ApJ*, **966**, 105

Thorsett S. E., Phillips J. A., 1992, *ApJ*, **387**, L69

Thorsett S. E., Arzoumanian Z., Camilo F., Lyne A. G., 1999, *ApJ*, **523**, 763

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

Verbiest J. P. W. et al., 2009, *MNRAS*, **400**, 951

Verbiest J. P. W. et al., 2016, *MNRAS*, **458**, 1267

Wolszczan A., 1994, *Science*, **264**, 538

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

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

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