# Combining the second data release of the European Pulsar Timing Array with low-frequency pulsar data

F Iraci, A. Chalumeau, C Tiburzi, J.P.W Verbiest, A Possenti, S.C Susarla, M.A Krishnakumar, G.M Shaifullah, J Antoniadis, M Bagchi, et al.

► To cite this version:

F Iraci, A. Chalumeau, C Tiburzi, J.P.W Verbiest, A Possenti, et al. Combining the second data release of the European Pulsar Timing Array with low-frequency pulsar data. Astron.Astrophys., 2025, 704, pp.A109. 10.1051/0004-6361/202555516. hal-05447007

HAL Id: hal-05447007

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

Submitted on 8 Jan 2026

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 CC BY 4.0 - Attribution - International License
# Combining the second data release of the European Pulsar Timing Array with low-frequency pulsar data

F. Iraci<sup>1,2,*</sup>, A. Chalumeau<sup>3,**</sup>, C. Tiburzi<sup>2</sup>, J. P. W. Verbiest<sup>4</sup>, A. Possenti<sup>2</sup>, S. C. Susarla<sup>5</sup>, M. A. Krishnakumar<sup>6</sup>, G. M. Shaifullah<sup>78</sup>, J. Antoniadis<sup>910</sup>, M. Bagchi<sup>1311</sup>, C. Bassa<sup>3</sup>, R. N. Caballero<sup>15</sup>, B. Cecconi<sup>1617</sup>, S. Chen<sup>1819</sup>, S. Chowdhury<sup>11</sup>, B. Ciardi<sup>14</sup>, I. Cognard<sup>20</sup>, S. Corbel<sup>21</sup>, S. Desai<sup>23</sup>, D. Deb<sup>11</sup>, J. Girard<sup>16</sup>, A. Golden<sup>5</sup>, J.-M. Grießmeier<sup>2017</sup>, L. Guillemot<sup>2017</sup>, M. Hoeft<sup>22</sup>, H. Hu<sup>10</sup>, F. Jankowski<sup>20</sup>, G. Janssen<sup>324</sup>, B. C. Joshi<sup>625</sup>, S. Kala<sup>11</sup>, E. Keane<sup>12</sup>, K. Nobelson<sup>2634</sup>, A. Konovalenko<sup>27</sup>, I. Kravtsov<sup>2720</sup>, M. Kramer<sup>10</sup>, K. Liu<sup>1819</sup>, A. Parthasarathy<sup>33110</sup>, P. Rana<sup>2930</sup>, D. Schwarz<sup>32</sup>, J. Singha<sup>28</sup>, A. Srivastava<sup>3323</sup>, K. Takahashi<sup>2634</sup>, P. Tarafdar<sup>2</sup>, G. Theureau<sup>2017</sup>, O. Ulyanov<sup>27</sup>, C. Vocks<sup>35</sup>, J. Wang<sup>36</sup>, V. Zakharenko<sup>27</sup>, and P. Zarka<sup>16</sup>

*(Affiliations can be found after the references)*

Received XX; accepted XX

## ABSTRACT

**Context.** Low radio frequency data are highly valuable for enhancing the sensitivity of pulsar timing arrays (PTAs) to propagation effects, such as dispersion measure (DM) variations. These low-frequency observations are particularly sensitive to DM fluctuations and can therefore significantly improve noise characterization in PTA datasets, which is essential for detecting the stochastic gravitational wave background (GWB).

**Aims.** For this work we incorporated for the first time low-frequency observations from LOFAR (100–200 MHz) and NenuFAR (30–90 MHz) into a PTA context by combining them with the most recent data release from the European and Indian PTAs (in particular, with the subsample labeled DR2new+, which includes only data from the new backends). This new combined dataset, labeled DR2low, consists of 12 pulsars observed over a time span of ~11 years, with radio frequencies spanning the range 30–2500 MHz. The expanded frequency coverage of DR2low enables us to update and refine the noise models of DR2new+, and this is crucial in order to increase the PTA sensitivity when searching for the stochastic gravitational wave background, which is the primary goal of PTA observations. This work is a milestone in the integration of low-frequency data into the upcoming third data release of the International PTA, which is posed to achieve the 5σ detection of the GWB.

**Methods.** We used the pulsar timing software packages LIBSTEMPO and ENTERPRISE to perform a noise analysis of DR2low. At first, we applied a standard noise model including red noise (RN) and time-variable dispersion measure (DMv) as power laws, with Fourier components up to 30 and 100 frequencies, respectively. Next, we performed a fully Bayesian model selection to identify the favored noise model for each pulsar and compute the Bayes factors across all combinations of RN, DMv, and a noise term with a chromatic index of 4 (CN<sub>4</sub>). Finally, we carried out a detailed analysis on the choice of the chromatic index for CN<sub>4</sub> and the contribution of the solar wind.

**Results.** The comparison between DR2low and DR2new+ using the standard noise model highlights the benefits of including low-frequency data. In particular, the additional frequency coverage improves the constraints on the DM variations and helps disentangle the DM and RN noise components in most pulsars. Through a Bayesian model selection, we found that the RN is required in the final model for 10 out of 12 pulsars, compared to only 5 in the DR2new+ dataset. The improved sensitivity to plasma effects provided by DR2low also favors the identification of significant CN<sub>4</sub> in eight pulsars, while none showed such evidence in DR2new+. The chromatic index of this process is consistent with four of the five pulsars, while two (PSRs J0030+0451 and J1022+1001) show significant deviations from such a value. We attribute this discrepancy to unmodeled contributions from the solar wind, especially because of the high DM sensitivity of LOFAR and NenuFAR and the high observing cadence provided by these datasets near solar conjunction. A dedicated analysis confirms that the current solar wind model fails to fully capture the observed delay, and residual power is absorbed into the DM component of the model.

**Key words.** pulsars – gravitational waves – propagation effects

## 1. Introduction

Pulsar timing array (PTA) experiments exploit the clock-like precision of the periodic broadband radio emission received from pulsars to detect low-frequency (on the order of the nHz order) gravitational waves (GWs). Unpredicted fluctuations in such a precise series of radiation pulses reflect any kind of perturbation that the pulsar signal experiences, ranging from irregularities in the pulsar's rotation to the influence of Galactic plasma and GWs. A stochastic and isotropic gravitational wave background (GWB), generated by a population of supermassive black hole

binaries (SMBHBs), introduces a temporally correlated noise in the time-of-arrival (ToA) fluctuations—also known as timing residuals—relative to the deterministic timing model (Sesana 2013). In the most simple case of circular SMBHBs, this noise has a red power spectrum, characterized by a steep spectral index of -13/3 (Phinney 2001). Moreover, these GW-induced fluctuations follow a spatial correlation among pulsars, referred to as the Hellings & Downs curve (Hellings & Downs 1983). To date, eight PTA collaborations exist in the world: the European PTA (EPTA Desvignes et al. 2016), the North-American Nanohertz Observatory for Gravitational waves (NANOGrav Ransom et al. 2019), the Parkes PTA (Manchester et al. 2013), the Indian PTA (InPTA Joshi et al. 2018), the Chinese PTA (CPTA Lee 2016),

<footnote>* e-mail: francesco.iraci@inaf.it</footnote>

<footnote>** e-mail: chalumeau@astron.nl</footnote>
the MeerKAT PTA (MeerTIME Miles et al. 2023), the African Pulsar Timing (APT<fnref n="1" />) consortium, and the International PTA (IPTA Verbiest et al. 2016). In 2023, five PTA collaborations presented results supporting the very first evidence for the presence of a GWB in PTA data (EPTA Collaboration et al. 2023b; Agazie et al. 2023; Reardon et al. 2023; Xu et al. 2023), followed a year later by the sixth collaboration (Miles et al. 2025). Although none of these studies proved the GWB presence with a significance higher than the required 5σ threshold (Allen et al. 2023), the results will be improved by the IPTA (formed by EPTA, NANOGrav, PPTA, InPTA, MeerTIME, and APT), which will combine the data from all of the individual PTAs. The EPTA will not only contribute with its second data release (described in EPTA Collaboration et al. (2023a) and EPTA Collaboration et al. (2023c)), but also by enhancing the IPTA dataset. The combination of PTA datasets within the IPTA framework improves the overall electromagnetic frequency-domain coverage, increases the sky coverage by including more pulsars, and benefits from a higher observing cadence. However, certain noise contributions remain challenging to model, particularly due to the limited availability of ultra-low-frequency observations spanning long time baselines.

The delay experienced by a radio wave as it propagates through the ionized interstellar medium (IISM) is inversely proportional to its frequency. Consequently, low-frequency pulsar data are especially sensitive to plasma-induced effects within the IISM through which the signal travels. The most relevant propagation effects for PTAs are dispersion and, albeit less impacting, scattering. Dispersion indicates a phenomenon for which the group velocity $v_g$ of the propagating light changes, depending on its radio frequency, due to the variation in the refractive index $\mu$. In particular,

$$ \mu \approx \sqrt{1 - \left(\frac{f_p}{v}\right)^2}, \qquad (1) $$

where $f_p$ is the plasma frequency and $v$ the observing (radio) frequency. The frequency-dependent delay $\delta t$ can be expressed as

$$ \delta t \approx \mathcal{D} \frac{\mathrm{DM}}{v^2}, \qquad (2) $$

where $\mathcal{D}$ is the dispersion constant, and DM the dispersion measure (DM), expressed as

$$ DM = \int_{\mathrm{LoS}} n_e dl, \qquad (3) $$

where LoS is the line of sight and $n_e$ the electron density.

Scattering implies that the inhomogeneity distribution in the IISM generates refraction of the incoming wavefronts, so that their propagation path changes. Hence, a certain fraction of rays reaches the observer with a delay relative to those that are not deviated and follow the LoS, causing the pulse profile to broaden and develop an exponential tail. The frequency-dependence of this tail is parameterized by the scattering timescale $\tau$, which scales as

$$ \tau \propto v^{-\gamma_s}, \qquad (4) $$

where $v$ is the observing frequency and $\gamma_s$ is the chromatic index. The parameter $\gamma_s$ varies with the kind of structure that induces the scattering and with the turbulence of the IISM. It is predicted

to be 4.4 for an extended Kolmogorov medium, although measurements deviate somewhat from this value (Bhat et al. 2004).

The parameters DM and $\tau$ (or even $\gamma_s$) can be time-dependent, as most pulsars are fast-moving objects due to the kick they receive during the supernova explosion (Hansen & Phinney 1997; Hobbs et al. 2005; Kapil et al. 2023). This motion causes the pulsar signal to cross different regions of the IISM, leading primarily to variable dispersion and, to a lesser degree, scattering, both of which can contribute significantly to red noise in PTAs. Since both of these phenomena are inversely dependent on the radio frequency of the propagating radiation, low-frequency data are fundamental to disentangling them from the achromatic contribution to the noise budget, including the GWB signal (Larsen et al. 2024; Iraci et al. 2024).

In this article we present the combination of the second data releases of the EPTA and InPTA (hereafter DR2) with the corresponding datasets from the LOw Frequency ARray (LOFAR, van Haarlem, M. P. et al. 2013) and the New extension in Nançay upgrading LOFAR (NenuFAR, Bondonneau et al. 2021), along with the associated noise analysis. In Section 2 we review the characteristics of DR2 and of the datasets collected using LOFAR and NenuFAR, and present their combination and basic timing properties in Section 3. In Section 4 we outline the approach and tools used to characterize the noise in the combined dataset. In Section 5 we present the results of the basic noise runs and compare them with those from EPTA DR2. In Section 6 we describe the more advanced model selection analysis, and highlight the differences with respect to the EPTA DR2 results. In Section 7 we examine the impact of the solar wind on the noise analysis. Finally, we draw our conclusions in Section 8.

## 2. Data description

### 2.1. EPTA-InPTA DR2

The data combination and noise analysis for the EPTA/InPTA are presented in EPTA Collaboration et al. (2023a) and EPTA Collaboration et al. (2023c) respectively. In this data release, EPTA and InPTA combined data collected from seven main European and Indian radio telescopes, i.e., the Nançay Radio Telescope (NRT), the 100-m Telescope at Effelsberg, the Lovell Telescope, the upgraded Giant Meterwave Radio Telescope (GMRT), the Westerbork Synthesis Radio Telescope, the Sardinia Radio Telescope and “Large European Array for Pulsars” (LEAP), from mid-1990s to 2021. Instead of analyzing the 42 pulsars presented in the first European data release (Desvignes et al. 2016), the authors selected 25 pulsars only, following an approach to maximize the sensitivity to the GWB and continuous GWs with a minimum number of pulsars (Speri et al. 2022). These pulsars and their characteristics are reported in Tables B1-B7 of EPTA Collaboration et al. (2023a). The data were combined using time offsets known as JUMPs with respect to the reference dataset, chosen to be the L-band NUPPI backend from the NRT. After that, initial timing solutions were obtained using the TEMPO2 timing software (Hobbs et al. 2006), and refined using the TEMPONEST Bayesian toolkit (Lentati et al. 2014) to introduce white noise (WN) parameters, red noise (RN; e.g., originating from rotational instabilities), and chromatic noise models (i.e., long-period noise that has some dependence on the observing frequency).

By considering only the new EPTA backends—more precise and more uniformly distributed across the frequency bands—combined with the InPTA dataset, the so-called DR2new+ dataset was constructed. This combination showed ev-

<footnote n="1">¹ https://africanpulsartiming.github.io</footnote>
idence for a GWB with a statistical significance exceeding 3.5σ (EPTA Collaboration et al. 2023b). The importance of radio-frequency coverage in the EPTA data was further investigated in (Ferranti, I. et al. 2025), which demonstrated that poor frequency coverage can significantly reduce the sensitivity of PTAs to a GWB signal. As discussed in Section 3, this work focuses exclusively on the DR2new+ dataset.

## 2.2. LOFAR

LOFAR (van Haarlem, M. P. et al. 2013) is a European low-frequency interferometer operating from approximately 30 to 240 MHz. Its largest concentration of modules (called stations) is located in the Netherlands, which also hosts the LOFAR core (i.e., a subsection of the interferometer that can form its own tied-array beam), while other international stations are distributed in Germany (six), Poland (three), France, the UK, Ireland, Sweden and Latvia (one each). Each station is an interferometer in itself, and is provided with a field of high-band antennae (HBA) operating from 100 to 240 MHz, and low-band antennae (LBA) operating below 100 MHz. Each international station can be detached from the largest interferometer and used as a stand-alone telescope.

Since 2013 a large pulsar monitoring campaign has been conducted with a main sub-band of the HBA fields (from about 100 to 200 MHz band) of the LOFAR core (LC) and the six German international stations (named DE601, DE602, DE603, DE604, DE605 and DE609) – some details are given by Konradtiev et al. (2016); Bilous et al. (2016); Donner et al. (2020); Donner (2022). More than 100 pulsars were observed with a monthly to bimonthly cadence with the core, and with a weekly cadence with the German stations, and similar observing programs started in many other countries hosting the international stations. Among the observed sources, 12 millisecond pulsars (MSPs) are in common with the twenty-five presented in EPTA DR2.

To calculate the ToAs from the described LOFAR dataset, we followed the methodology outlined in Donner et al. (2020) and Tiburzi et al. (2021). This involved constructing a frequency-resolved, high signal-to-noise ratio (S/N) profile template using either a single bright observation or an average of several high S/N profiles. We then cross-correlated each channel from every observation with the corresponding template profile at the that frequency, using the method described by Taylor (1992) but with the uncertainties calculated through a Monte-Carlo routine as recommended by Verbiest et al. (2016). Based on the work by Donner et al. (2020), we generally used eight frequency channels (and we therefore calculated 8 initial ToAs per observation per pulsar), but this can vary depending on the pulsar's S/N. The ToAs are then checked for outliers using a Huber regression-based algorithm described in (Tiburzi et al. 2019).

## 2.3. NenuFAR

NenuFAR (Bondonneau et al. 2021) is a low-frequency French interferometer that covers the 10–85 MHz range, featuring enhanced sensitivity and a significantly improved passband compared to the LBA field of any LOFAR station, including those in the core. As part of the NenuFAR Pulsar Key Project, a pulsar monitoring program has been ongoing since 2019, following a bi-weekly to monthly cadence and targeting approximately 40 pulsars — including two MSPs also present in the EPTA DR2 list: PSRs J0030+0451 and J1022+1001.

To calculate the ToAs for the NenuFAR pulsars, we first constructed a frequency-resolved, high S/N template profile based on the DM-corrected average of all observations taken when the pulsar was located more than 45° from the Sun. This selection helps avoid contamination from solar dispersion effects, which are non-trivial to correct. As with the LOFAR dataset, ToAs are then obtained by a) partial averaging in frequency of both the observations and the template in order to achieve eight frequency channels with reasonable S/N; b) channel-wise cross correlating the frequency resolved template with the individual observations; and c) cleaning of outliers using the aforementioned software based on the Huber-regression scheme.

## 3. Data combination and basic timing

In the context of PTA studies, combining data from different observatories is essential to increase the overall data volume (number of observations), the frequency coverage and the time span, thereby improving sensitivity to the stochastic GWB (Siemens et al. 2013). This section outlines the procedures used to combine the low radio frequency datasets of LOFAR and NenuFAR with DR2new+ and describes the basic timing analysis performed to prepare the data for the noise analysis. In Figure 1 we display the sky positions of the 25 MSPs from the DR2new+, emphasizing the names of the 12 sources studied in this work. For the sake of clarity, the dataset composed of the combination of DR2new+, LOFAR, and NenuFAR is referred to as DR2low in the rest of the document.

To address systematic time delays between observing systems, we fit for time offsets (JUMPs) for each observing station. Specifically, we included one JUMP for each international LOFAR station, one for the LOFAR core, and one for NenuFAR when present. Another important aspect to address when combining frequency-resolved data (one ToA per frequency channel, timed against the relevant frequency channels of a frequency-resolved template) is the difference between the DMs that were used to de-disperse the template profiles. More precisely put, when timing frequency-resolved data against a frequency-resolved template, the ToAs of each frequency channel are measured relative to a phase reference that is defined by the relevant frequency channel of the template. Through dedispersion of the template profile, the relative phases of the different channels are changed with respect to one another, introducing a frequency dependence of the ToA's phase reference. This frequency-dependent phase reference is different for each template used and is defined by the combination of the DM used to dedisperse the template and the actual interstellar DM affecting the data that were used to construct the template. Consequently, this effect is fully covariant (and indistinguishable from) a different DM affecting the datasets derived from a given template profile. In order to account for this effect, we therefore applied a constant DM offset to all but one of the frequency-resolved datasets. Specifically, in our analysis, frequency-resolved ToAs were available for the LOFAR, NenuFAR, and InPTA data, the latter of which was used as a reference for the other two to compute the DM offset.

Once the datasets have been prepared, we started from the EPTA/InPTA timing model and included the additional JUMPs. After fitting the combined dataset, we observed that only a few key parameters exhibit significant changes—namely, the DM and its derivative (if present), the rotational period and its derivative, and the electron density at 1 AU to model the solar wind contribution. All other timing model parameters remain consis-
<figure><img src="image_2.png" /><figcaption>Fig. 1. Mollweide sky projection of the 25 DR2new+ pulsars (blue stars) in equatorial coordinates. The red circles and green squares display sources observed by respectively LOFAR and NenuFAR. The names of the pulsars studied in this work are written next to the position. The orange line (and area) are representing the Sun’s position (± 20 degrees angle) along the year, denoting the influence zone of the solar wind.</figcaption></figure>

tent with their initial values, likely due to the larger ToA uncertainties in the low-frequency data.

In Figures D.1 and D.2, we show the frequency coverage and timing residuals for each of the 12 combined pulsars. The ToAs, timing models, and low-frequency templates for all the pulsars are available on Zenodo<fnref n="2" />.

# 4. Noise analysis: Methodology

## 4.1. Framework and tools

Having the combined dataset, we performed various iterations of noise analysis (EPTA Collaboration et al. 2023c), compared the results after including low radio frequency (< 300 MHz) data, and assessed its improvement on constraining the models. The noise analyses are performed in a Bayesian framework, using the software LIBSTEMPO<fnref n="3" /> to extract data products from TEMPO2. The software ENTERPRISE (Ellis et al. 2020) and its complementary module ENTERPRISE_EXTENSIONS (Taylor et al. 2021) are used to model the noise components and build the likelihood and prior functions, applied to the MCMC sampler PTMCMCSAMPLER (Ellis & van Haasteren 2017).

The remaining signals from the timing analysis described in Section 3, referred to as timing residuals $\delta t$, contain any effects impacting ToAs that are not included in the deterministic timing model. The analysis is performed by including the noise components in the multivariate Gaussian time-domain likelihood function, which thus assumes a Gaussian distribution of the noise-reduced timing residuals, as

$$ \log L(\delta t | \theta, \theta') = -\frac{1}{2} ([\delta t - y(\theta')]^T C^{-1}(\theta) [\delta t - y(\theta')]) + \text{cst.}, $$

(5)

where the timing residuals $\delta t$ are reduced by any deterministic time-domain waveform $y$ with parameters $\theta'$, and where the covariance matrix is composed of the sum of all the covariance matrices from all the stochastic component $i$ with parameters $\theta_i$ such that $C(\theta) = \sum_{\lambda} C^{\lambda}(\theta_{\lambda})$. The following subsection describes the noise components used in this work and shows how they are included in the likelihood function.

## 4.2. Model components

The noise components used in this work are similar to those from EPTA Collaboration et al. (2023c), and allow a one-to-one comparison. However, given the high sensitivity of low-frequency data to the solar wind dispersion, we add a refined model (Hazboun et al. 2022; Susarla et al. 2024) to include its constant and time-varying components that are fitted simultaneously with other noise components during the noise analysis. Table (A.1) describes all the components, their main properties and their parameter priors.

### 4.2.1. White noise

In the case of a perfect description of the ToAs from the timing model, the timing residuals divided by their uncertainties would follow a standard normal distribution. However, the measurement of the ToAs and their uncertainties (that estimate the so-called radiometer noise) are in practice subject to errors which have been found to be system-dependent. These errors are most commonly stochastic and time-uncorrelated, thus being "white noise" since they correspond to a flat power spectral density (PSD).

In the case of the observing systems of DR2new+, we use the same approach as in EPTA Collaboration et al. (2023c), with an

<footnote n="2">² https://doi.org/10.5281/zenodo.16794085</footnote>

<footnote n="3">³ https://github.com/vallis/libstempo</footnote>
Error FACtor (EFAC) applied to the provided ToA uncertainty $\sigma_{init,i}$ and an Error QUADRature (EQUAD) term added to them in quadrature, for each observing system $\alpha$. For the ToAs provided by LOFAR and NenuFAR, an Error CORRelated term (ECORR) is also added for each observing system to account for a white noise that is fully correlated among ToAs measured at the same epoch but at different radio frequencies, this type of error is typically associated with pulse phase jitter, but also any frequency-correlated systematics (Cordes & Shannon 2010; Shannon & Cordes 2012; Liu et al. 2012) may affect ECORR value. The white noise terms are added to the likelihood covariance matrix C

$$C_{\alpha}^{\text{WN}} = \begin{cases} (\text{EF}_{\alpha}^2 \times \sigma_{\text{init},i}^2 + \text{EQ}_{\alpha}^2) \delta_{ij} \delta_{\alpha\beta}, & \text{if } \alpha \notin \{\text{LOFAR, NenuFAR}\}, \\ (\text{EF}_{\alpha}^2 \times \sigma_{\text{init},i}^2 + \text{EQ}_{\alpha}^2) \delta_{ij} + \text{EC}_{\alpha}^2 \delta_{\text{E(i)E(j)}} \delta_{\alpha\beta}, & \text{else,} \end{cases} \quad (6)$$

where EF, EQ, and EC stand for EFAC, EQUAD, and ECORR, respectively, the indices i, j and $\alpha, \beta$ range over the ToAs and the observing systems respectively and E(i) indexes the epoch of the ToA i. The symbol $\delta$ refers to the Kronecker delta, as in the following parts of this document.

### 4.2.2. Time-correlated stochastic signals

The time-correlated stochastic signals are modeled as Gaussian processes (GPs) following the approach described in van Haasteren & Vallisneri (2014), using the covariance-function expansion to express the likelihood component of the noise $\lambda$ as

$$C^{\lambda,GP}(t_i, t_j) = \sum_{\mu, \nu} F_{\mu}^{\lambda}(t) \Phi_{\mu\nu}^{\lambda} F_{\nu}^{\lambda}(t'), \quad (7)$$

where $F$ is a set of basis functions, $\Phi$ is the covariance matrix of the Gaussian process in the weight-space view (Rasmussen & Williams 2006), and $\mu, \nu$ both range over the number of basis functions.

In this work, this approach is used to model both the achromatic red noise and time-varying chromatic signals arising from the heterogeneous nature of the IISM. These include DM variations not captured by the DM1 and DM2 parameters, as well as an additional chromatic delay with a $v^{-4}$ frequency dependence, consistent with the scattering expected under the thin screen approximation with same-size inhomogeneities (Cordes & Shannon 2010).<fnref n="4" /> These three noise components (respectively labeled RN, DMv and CN₄) are modeled in the PSD domain from an incomplete Fourier basis approach where the basis functions are built as a set of sine and cosine functions for each PSD frequency $f$, and the Gaussian process covariance matrices are defined as

$$\Phi_{\mu\nu}^{\lambda} = S^{\lambda}(f_{\mu}) \Delta f \delta_{\mu\nu}, \quad (8)$$

where $S(f_\mu)$ is the PSDs amplitude at frequency $f_\mu$ and $\Delta f$ is the spectral resolution, equal to the inverse of the total time span of the time series.

Unless it is specified, the PSD of these three components are modeled as power laws:

$$S^{\lambda}(f) = \frac{A_{\lambda}^{2}}{12\pi^{2}} \left(\frac{f}{f_{\text{yr}}}\right)^{-\gamma_{\lambda}} f_{\text{yr}}^{-3}. \quad (9)$$

Here $f_{\text{yr}}$ refers to the frequency of $1\text{yr}^{-1}$ and $A_\lambda$, $\gamma_\lambda$ are respectively the amplitude referred at $f_{\text{yr}}$ and the spectral index of the noise component $\lambda$.

In our case, the different noise components are $\lambda = \{\text{RN, DMv, CN}_4\}$, which differ from their chromaticity (i.e., dependence on the observed radio-frequency $v$) that is added in the basis functions as $F_\mu^\lambda(t) \propto (v/v_{\text{ref}})^{-\chi}$, where $\chi = 0, 2, 4$ for RN, DMv and CN₄ respectively, and $v_{\text{ref}}$ being a reference frequency set at 1.4 GHz in this work. While the chromatic index $\chi$ is fixed for these components, it is also possible to fit for this parameter simultaneously with other parameters and hyper-parameters of the model.

The timing model parameter uncertainties are also modeled as GPs also by using the covariance-function expansion (cf. Eq. 7), but they are analytically marginalized over. The basis functions are the timing model parameter design matrix and the Gaussian process covariance matrix is a diagonal matrix with very large values (Chalumeau et al. 2022).

### 4.2.3. Solar wind

The solar wind (SW) is an outflow of charged particles from the Sun which can significantly affect the DM especially for observations of pulsars near the ecliptic plane. Properly modeling the impact of SW is crucial for precise timing analysis, especially in presence of very low-frequency data. To account for this, we adopted a detailed SW model following Susarla et al. (2024) which contains a deterministic and a stochastic component.

We include a deterministic signal which accounts for the electron content $n_e$ as a function of the angular separation between the pulsar and the Sun sampling it from

$$DM_{\text{SW}} = n_e \frac{\rho}{r_e \sin \rho} [1 \text{ AU}]^2, \quad (10)$$

where $\rho$ is the pulsar-Sun-observer angle, $n_e$ is the electron density at 1AU (denoted as NE_SW in TEMPO2), and $r_e$ is the observatory-Sun distance (Tiburzi et al. 2019; Edwards et al. 2006). Additionally, we account for secular variations in the solar wind by incorporating a stochastic SW signature modeled with GPs. This approach, implemented as SWGP, is available in the chromatic noise module of enterprise_extensions. It uses a power spectral density described by the following power-law formula:

$$S_{\text{SW}} = A_{\text{SW}}^2 \left(\frac{f}{\text{yr}^{-1}}\right)^{-\gamma_{\text{SW}}} \text{yr}^{-3}. \quad (11)$$

Here $S_{\text{SW}}$ is the power spectral density for SWGP, $A_{\text{SW}}$ is the SW amplitude at a frequency of $1/1\text{yr}$, $\gamma_{\text{SW}}$ is the spectral index, and $f$ is the Fourier frequency. We use the default value for the number of Fourier bins that goes up to $1.5\text{yr}^{-1}$.

In this work, the SW model is specific to each pulsar and primarily depends on factors such as the ecliptic latitude, the ToA precision, and the observing cadence. In Table A.2 we summarize how we modeled the SW for each pulsar, and provide a comparison with the model used in the EPTA DR2.

<footnote n="4">⁴ We note that the chromaticity of IISM-induced delays strongly depends on the specific properties of the electron distribution and may vary between pulsars. For instance, a $v^{-4.4}$ dependence is expected for a turbulent extended medium described by Kolmogorov theory, while refractive propagation effects could produce a steeper scaling, up to $v^{-6.4}$ (Shannon & Cordes 2016).</footnote>
In the EPTA analysis, a fixed value of NE_SW = 7.9 was used for all pulsars, except PSR J1022+1001, and its error was marginalized over during the noise analysis, as for the other timing model parameters (the limitations of this approach have already been addressed by Liu et al. (2025)). For the combined dataset, the treatment of the SW contribution varies across pulsars depending on their ecliptic latitude (ELAT) and the availability of high-cadence observations. Pulsars with very low ELAT, such as PSRs J1022+1001 and J0030+0451, clearly exhibit solar wind signatures in their timing residuals. Although Susarla et al. (2024) recommend including only the deterministic component of the SW model for pulsars with ELAT below 3°, we also included the stochastic term, SWGP, in these two cases, due to the availability of NenuFAR ToAs, which were not part of the previous analysis. We applied the same extended model (including SWGP) to PSR J1744-1134, whose ELAT is above 3°, and to PSR J0751+1807, which was not included in Susarla et al. (2024). The only exception is PSR J1730-2304, for which we excluded the stochastic SW contribution due to the lack of high-cadence observations around solar conjunction epochs.

### 4.2.4. Exponential dip

It is known a priori for PSR J1713+0747 that timing residuals contain multiple exponential dips (Coles et al. 2015; Lam et al. 2018; Xu et al. 2021): a brief radio-frequency dependent drop (where ToAs are measured before the time predicted by the timing model) followed by an exponential recovery for several weeks/months. Depending on the corresponding chromatic index $\chi^E$, the origin of these events are either considered to be related with a drop in the IISM content in the line-of-sight (Lam et al. 2018, if $\chi^E = 2$), or related to magnetospheric effects if $\chi < 2$ and if profile changes are observed at the time of the event. The latter scenario is favored for the only event present in the dataset used in this work that happened in 2016 (~ MJD 57510), where Goncharov et al. (2020) found a significant profile change at this epoch and a chromatic index measured at $1.15^{+0.18}_{-0.19}$ with data from the Murriyang (Parkes) radio telescope, also being consistent with the value of $1.00^{+0.56}_{-0.49}$ measured with EPTA data (Chalumeau et al. 2022). As in EPTA Collaboration et al. (2023c), we model this effect as a deterministic signal with a waveform defined as

$$y^E(A_E, \tau_E, t_{0,E}, \chi_E) = \begin{cases} 0, & \text{if } t < t_{0,E} \\ A_E (\nu/1.4 \text{ GHz})^{-\chi_E} \exp\left(-\frac{t-t_{0,E}}{\tau_E}\right), & \text{if } t \ge t_{0,E} \end{cases}, \quad (12)$$

where $A^E$, $\chi^E$, $t^E$, and $\tau^E$ are respectively the amplitude (in residual units), the chromatic index, the reference epoch (in MJD) and the relaxation time (in days) of the exponential dip modeled for a time $t$ and a radio-frequency $\nu$.

### 4.3. Parameter estimation and model selection

In this work, we performed several analyses per pulsar, including parameter estimation and model selection methods, in order to (1) get information on the properties of the chosen noise components described below, (2) analyze the correlations between them, and (3) obtain the most favored combination among the

noise components, which could be compared against the results published in EPTA Collaboration et al. (2023c). This enabled us to understand the implications of including LOFAR and NenuFAR data with DR2new+.

For the sake of clarity, since we consider multiple noise models per pulsar, we describe the noise models from their Gaussian process components (apart from the timing model uncertainties), and the number of PSD frequencies used. As an example, the noise model commonly adopted in PTA analyses includes marginalized timing model parameter uncertainties, white noise, an achromatic RN component, and DM variations. The RN and DM components are both modeled as a power-law processes, with 30 and 100 frequency bins respectively. This noise model will be referred to as RN30+DMv100 throughout this document.

For the model selection, we applied a similar approach as the one described in EPTA Collaboration et al. (2023c), where we defined a set of candidate noise models, fit for the number of frequency bins (i.e., $\nu$ in Eq. 7) used to model the power-law PSD GPs, and compute a Bayes factor among the different chosen models with their most favored number of frequency bins. The selected models are WN (white noise only), WN+RN, WN+DMv, WN+RN+DMv, WN+DMv+CN$_4$, WN+RN+DMv+CN$_4$, identical to the models in EPTA Collaboration et al. (2023c).

In all runs, we included the SW model described in Table A.2, which consists of a deterministic component, n_earth, and a stochastic term, SWGP for a specific set of pulsars. Previous analyses based on LOFAR-only data indicated strong evidence in favor of this choice, with significantly improved model fit when the SWGP component is included. While this model may not fully capture all aspects of the solar wind signal-particularly around annual harmonics-it generally leads to a more accurate description of the data, with notably reduced white noise levels in several systems.

## 5. Noise analysis - Results with standard models

We performed a parameter estimation analysis following the approach described in Section 4 on two datasets, DR2new+ and DR2low, by using standard noise models, composed of TM, WN, RN30, DMv100 and the SW following the description above.

In Figure 2 we show the credible intervals (at confidence levels equivalent to 1, 2 and $3-\sigma$ in a Gaussian distribution) of the noise parameter posterior distributions for DR2new+ (orange) and DR2low (blue). For a better clarity, the color is lighter for the $3-\sigma$ contour levels. The left-hand column presents the DM variations posterior distributions, with $\log_{10} A_{DM}$ on the y-axis and $\gamma_{DM}$ on the x-axis, while the right-hand column displays the corresponding parameters for the RN. For reference, we display the two principal directions<fnref n="5" /> of the posteriors distributions for the same parameters obtained in Agazie et al. (2024) on the NANOGrav 15-year dataset (Agazie et al. 2023) and PPTA DR3 (Zic et al. 2023) with solid and dashed lines respectively. The green vertical lines indicate the spectral index values expected from Kolmogorov turbulence (for the DM variations; left column), and the one predicted from a population of GW-driven supermassive black hole binaries with circular orbits (for the RN; right column).

<footnote n="5">⁵ the principal directions are evaluated by projecting the distributions onto the eigenvectors, with asymmetric scaling factors corresponding to $3-\sigma$ in a Gaussian distribution.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Spectral index vs. amplitude posterior distributions of the power-law PSD for DMv (left) and RN (right) obtained for 12 pulsars from DR2low (blue) and DR2new+ (orange) using the standard noise models. The green vertical lines show the predicted spectral index from a Kolmogorov turbulence in the IISM for DMv (left), and from circular and GW-driven SMBHBs (right). The black crosses display the 99.7% credible intervals of posterior distributions obtained with NANOGrav-15yr and PPTA DR3, displaying the principal directions.</figcaption></figure>

## 5.1. Observations for the DM variations

### 5.1.1. Comparing precision between DR2new+ and DR2low

As expected, adding very low-frequency data from LOFAR and NenuFAR provides stronger constraints of the DMv parameters (i.e., DR2low having narrower posterior distributions than DR2new+) for all the 12 pulsars, apart for PSRs J1022+1001, J1713+0747 and J1918-0642 for which they remain similar. In particular, we observed a very significant improvement for PSRs J0030+0451, J1024-0719, J1730-2304 and J1738+0333, which remained largely unconstrained in the DR2new+ analysis.

### 5.1.2. Consistency between DR2new+ and DR2low

The new constraints remain fully consistent with DR2new+ for 7 of the 12 pulsars: PSRs J0030+0451, J0751+1807, J1024-0719, J1730-2304, J1738+0333, J1857+0943 and J1918-0642, while a slight shift toward lower amplitudes is observed for PSR J1713+0747, with a tension of only ~ 0.77σ according to the description of a tension metric described in EPTA Collaboration et al. (2023c) based on Raveri & Doux (2021). For PSR J1918-0642, we actually observed a slight increase of the wide posterior distribution at the high spectral index-low amplitudes in the parameter space. However, the constraints are significantly different for PSRs J1012+5307, J1022+1001, J1640+2224, and J1744-1134, where the posterior distributions are shifted toward higher spectral indices (from about 1 to 2) and lower amplitudes (logarithmic values ~0.3 to 0.5). These discrepancies may arise from several factors. First, the DR2new+ dataset mainly includes observations at high frequencies, which may limit its ability to disentangle chromatic (DMv) and achromatic (RN) noise components. For PSR J1640+2224, DR2new+ contains only data at L and S bands, making it inherently unreliable for DM characterization. Furthermore, PSRs J1022+1001, J1744-1134, and even J1012+5307—despite its relatively high ELAT of 38.76°—are all affected by the SW signal, which, as discussed in Section 7, biases the DM recovery. Another, although less-likely, explanation could be found in DM chromaticity (Cordes et al. 2016).

### 5.1.3. Comparison with predicted spectral index

The spectral indices inferred from both datasets are consistent with the expected Kolmogorov value of 8/3 (indicated by green vertical lines) for six pulsars: PSRs J0751+1807, J1024-0719, J1730-2304, J1738+0333, J1857+0943, and J1918-0642. Interestingly, while this was not the case with DR2new+, the spectral index becomes consistent with 8/3 for PSRs J1012+5307 and J1640+2224, and shows a mild agreement for J1713+0747. Although the constraints for PSRs J0030+0451, J1022+1001, and J1744-1134 remain below the expected value, the distributions for the latter two exhibit an increasing trend. While the increase in spectral index from DR2new+ to DR2low could be attributed to frequency-dependent ray path averaging, we highlight the significant impact of solar wind effects on these three pulsars. Such effects may bias the estimation of DM variation parameters due to limitations in current solar wind models (cf. Section 7), potentially resulting in spectral index constraints below the expected 8/3 value. In the context of pulsar timing data, a spectral index lower than the true value typically indicates excess power at higher frequencies in the PSD, which may arise from unmodeled contributions.
to dispersion variations, such as those induced by the solar wind.

### 5.1.4. Comparison with NANOGrav 15-yr

We noticed a clear consistency of the posterior distributions from DR2low with the constraints obtained with NANOGrav 15-yr dataset (black solid lines) for seven of the eleven considered pulsars: PSRs J0030+0451, J1024-0719, J1640+2224, J1713+0747, J1738+0333, J1857+0943 and J1918-0642, providing an improvement in the robustness of our results. The constraints for PSR J1744-1144 are consistent for the spectral index, but DR2low has a lower amplitude. We also observed a strong difference for PSR J1012+5307, where NANOGrav 15-yr appears to be more consistent with DR2new+, with higher amplitude and lower spectral index. Such property could be explained by the fact that NANOGrav 15-yr shares the same observing radio frequencies with DR2new+ for this pulsar. Although, the NANOGrav dataset does not fully overlap in time with DR2low since it contains ToAs from 2004 for this pulsar (starting almost 7 years before DR2low). Thus, any effects implying non-stationary DM variations would yield to different spectral properties as well. We attributed the large differences for PSRs J1022+1001 and J1730-2304 to the short time span of NANOGrav for these two pulsars, with about 5 and 4 years respectively.

### 5.1.5. Comparison with PPTA DR3

The posterior distributions from PPTA DR3 (dashed lines) are fully consistent with DR2low for four of the seven overlapping pulsars: PSRs J1024-0719, J1713+0747, J1730-2304 and J1857+0943, while mildly consistent with PSR J1744-1134, such as for NANOGrav. However, the posterior distribution for PSR J1022+1001 appears to be more consistent with DR2new+, with a lower spectral index than DR2low. The very poor constraints for PSR J0030+0451 might be due to the short time span for PPTA data, with about 3 years.

## 5.2. Observations for the RN

### 5.2.1. Comparing precision between DR2new+ and DR2low

We noted a large impact on the RN constraints after adding data from LOFAR and NenuFAR. The posterior distributions become significantly narrower for PSRs J1022+1001, J1713+0747, J1730-2304, J1744-1134 and J1918-0642, where even the $3-\sigma$ contour levels appear now constrained, while it was only the case for PSR J1713+0747 with DR2new+. We also observed an improvement of precision for PSRs J0030+0451, J0751+1807, J1640+2224 and J1857+0943, where a peak appears more prominent in contour levels $\le 2-\sigma$. Furthermore, we noted an improvement in the upper limits for PSRs J1024-0719 and J1738+0333, for which including LOFAR data helps reducing the covariances between RN and DMv (see the following paragraph for more details). We noted a similar precision between both datasets only for PSR J1012+5307.

### 5.2.2. Consistency between DR2new+ and DR2low

The RN distributions from DR2low are fully consistent with DR2new+ for PSRs J0030+0451 and J1713+0747, where DR2low results in a more prominent peak. However, we obtained a constrained RN for PSRs J0751+1807, J1640+2224, J1730-2304, J1744-1134, J1857+0943 and J1918-0642, while it appears fully unconstrained for DR2new+. This improved constraint on the RN component is a new result, but it remains highly consistent with DR2new+ for PSRs J0751+1807, J1640+2224, J1730-2304, and J1857+0943, with tensions of only $0.35\sigma$, $0.06\sigma$, $0.04\sigma$, and $0.18\sigma$, respectively. It is in mild tension for PSR J1918-0642 ($1.02\sigma$) and PSR J1744-1134 ($1.32\sigma$), where DR2low yields a RN highly constrained at low spectral index ($\gamma_{RN} \sim 1.1$) and high amplitudes ($\log_{10}A_{RN} \sim -13.3$). Conversely, while RN was slightly constrained for PSRs J1024-0719 and J1738+0333 with DR2new+, it appears significantly reduced with DR2low. This behavior could be observed in Figure 3 (top and bottom plots for PSRs J1024-0719 and J1738+0333 respectively), where we observed strong "L-shaped" correlations between RN and DMv, as described in Ferranti, I. et al. (2025) in the case of poor radio frequency coverage. Adding lower frequencies and thus improving the precision on DM variations has removed these correlations, avoiding a false positive for the RN in these cases. For PSR J1012+5307, the distribution is slightly shifted compared to DR2new+, mainly toward a higher spectral index ($\sim$ from 1 to 2), with a mild tension of about $1.46\sigma$. Interestingly, we observed for the first time a bimodal posterior distribution for the RN of PSR J1022+1001, which lies inside the poorly constrained distribution for DR2new+ (tension at only $0.48\sigma$). We observed a peak constrained with a low spectral index ($\gamma_{RN} \sim 0$), while the second one favors a steeper power law ($\gamma_{RN} \sim 5$). This feature could reveal a disentanglement between short-term signals with low spectral index, that could be related with pulse shape variations (Liu et al. 2015; Padmanabh et al. 2020), and longer term processes that could be from different origin (e.g., nHz gravitational waves or spin noise). It also means that a careful treatment may need to be applied to properly model the RN of this pulsar, as it was already pointed out in EPTA Collaboration et al. (2023c).

### 5.2.3. Comparison with a predicted value

As described in Section 4, the achromatic red noise might have multiple origins (e.g., pulsar spin noise, nanoHertz gravitational waves), which implies some caveats about comparing it against predicted values from single pulsar noise analyses. Nevertheless, it is in practice interesting to apply such comparisons with predicted value for the GWB (see the recent summary by Verbiest et al. 2024), in particular regarding the slight tensions from the most recent results published by PTAs. In this case, we compared the constraints on the spectral index of the RN against the predicted value of 13/3 from the simplistic case of a GWB caused by a population of GW-driven and circular SMBHBs (green vertical line). In the case of DR2low, the spectral index value of 13/3 is within the $1\sigma$ contour levels for PSRs J1024-0719, J1640+2224, J1738+0333, J1857+0943 and J1918-0642, and withing the $2\sigma$ contour levels for PSRs J0030+0451, J0751+1807, J1022+1001 (second peak), J1713+0747 and J1730-2304. It is however strongly inconsistent with the lower spectral indices of PSRs J1012+5307 and J1744-1134 (with 1D marginalized medians resp. equal to 1.75
and 1.10), for which the major contribution to the RN might be
from another origin (e.g., pulsar-related effects).

5.2.4. Comparison with other PTAs

As mentioned in Subsection 5.2, the short time spans for PSRs J1730-2304 and J1022+1001 with NANOGrav 15-yr (about 4 and 5 years respectively), and for PSR J0030+0451 (~3 years) with PPTA DR3 are the cause of their large uncertainties. We found high consistency in the RN properties between DR2low and other PTAs for PSRs J0030+0451, J1012+5307, J1022+1001 (where only the peak at low spectral index shows up for PPTA), J1024-0719, J1640+2224, J1730-2304, J1738+0333, J1857+0943 and J1918-0642, with tensions lower than 0.6σ. We observed mild tensions for PSR J1713+0747 (1.85σ and 1.50σ for NANOGrav and PPTA, respectively), where the constraints of DR2low imply larger amplitudes than the other PTAs. We emphasize that the shorter time spans for DR2low (half of the other's datasets) could be the cause of such inconsistencies. Given the differences in the timing and noise models employed, we recommend a thorough investigation of the combined datasets to evaluate their consistency for these pulsars in future studies. For PSR J1744-1134, we observed consistent results with PPTA (~0.43σ) but a moderate tension with NANOGrav (~2.18σ) which favors a higher spectral index and a lower amplitude as obtained in DR2new+. Inconsistencies between PTAs were already found in Agazie et al. (2024) for this pulsar, and hence a combination of the datasets in an IPTA context might help understand these features better.

**6. Noise analysis - Results from model selection**

We performed an analysis aimed at finding the noise model that is most favored by the DR2low data, following a similar approach as EPTA Collaboration et al. (2023c). The posterior distributions of ECORR parameters for every pulsars/systems are shown in Figure B.1, where only five pulsars have unconstrained distributions (PSRs J0571+1807, J1024-0719, J1738+0333, J1857+0943 and J1918-0642). The SW parameters of PSRs J0030+0451, J0751+1807, J1022+1001 and J1744-1134 are presented in Figure C.1. In this section, we first intend to provide a comparison with the results that were obtained with DR2new+, then interpret the results for DR2low.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Posterior distributions of spectral indices (y) and amplitudes (A) of DMv and RN for PSRs J1024-0719 (top) and J1738+0333 (bottom) obtained with DR2low (blue) and DR2new+ (red), using the standard noise models.</figcaption></figure>

6.1. Comparing the favored models of DR2new+ and DR2low

Table 1 presents the final noise models obtained with DR2low (left) from this work, and DR2new+ (right) from the work in EPTA Collaboration et al. (2023c). It explicitly describes the noise components (RN, DMv and CN₄) modeled as GPs with a power-law power spectral density, along with the corresponding favored number of frequency bins. It is important to note that model selection was not applied to PSR J1022+1001 for the EPTA DR2 datasets (including DR2new+), as it was found that a single power-law model was insufficient to adequately describe the data across both low and high frequencies. Consequently, we excluded this pulsar from the comparison with DR2new+, as the RN30+DMv100 model was retained for this dataset.

We observed a significant change in the final noise models
between the two datasets, where the favored models are the same
only for PSR J1024-0719, with just DMv. Firstly, we observed
that the RN is favored for 10 of the 12 pulsars with DR2low,

while it was only the case for five of them with DR2new+. This
behavior supports the notion that incorporating low-frequency
data can enhance PTA sensitivity. As demonstrated in EPTA Col-
laboration et al. (2023b), the inclusion of InPTA DR1—featuring
observations in the 300–500 MHz range—alongside EPTA DR2
led to a modest increase in the evidence for a GWB signal. In
the case of LOFAR and NenuFAR, we anticipated a stronger im-
provement since the frequencies are lower (< 200 MHz), and
also because it contains a longer time span and higher cadence,
allowing us to better sample signals in the power spectral do-
main. Additionally, we observed the presence of higher order
chromatic noise for eight pulsars, while none were observed with
DR2new+. Lastly, while the DM variations are favored for all pul-
sars apart from PSR J0030+0451 with DR2new+, they are now
present for all of them with DR2low. This property is not sur-
prising, however, thanks to the high sensitivity of low-frequency
data toward DM variations.
<table><caption>**Table 1.** Comparison of favored noise models between DR2low and DR2new+.</caption><thead><tr><th>Pulsar</th><th>DR2low</th><th>DR2new+</th></tr></thead><tbody><tr><td>J0030+0451</td><td>RN13 DMv100 CN<sub>4</sub>77</td><td>RN10 X X</td></tr><tr><td>J0751+1807</td><td>X DMv92 CN<sub>4</sub>38</td><td>X DMv50 X</td></tr><tr><td>J1012+5307</td><td>RN87 DMv20 CN<sub>4</sub>100</td><td>RN100 DMv47 X</td></tr><tr><td>J1022+1001</td><td>RN98 DMv99 CN<sub>4</sub>99</td><td>RN30 DMv100 X</td></tr><tr><td>J1024-0719</td><td>X DMv18 X</td><td>X DMv15 X</td></tr><tr><td>J1640+2224</td><td>RN54 DMv40 CN<sub>4</sub>40</td><td>X DMv47 X</td></tr><tr><td>J1713+0747</td><td>RN73 DMv74 CN<sub>4</sub>12</td><td>RN73 DMv72 X</td></tr><tr><td>J1730-2304</td><td>RN18 DMv52 X</td><td>X DMv11 X</td></tr><tr><td>J1738+0333</td><td>RN29 DMv28 X</td><td>X DMv24 X</td></tr><tr><td>J1744-1134</td><td>RN22 DMv100 CN<sub>4</sub>31</td><td>RN19 DMv98 X</td></tr><tr><td>J1857+0943</td><td>RN87 DMv99 X</td><td>X DMv92 X</td></tr><tr><td>J1918-0642</td><td>RN36 DMv24 CN<sub>4</sub>42</td><td>X DMv75 X</td></tr></tbody></table>

For the selected number of frequency bins, we observed comparable results for the RN component, with relatively low values (~10–22) for PSRs J0030+0451 and J1744–1134, and significantly higher values (~87–100) for PSRs J1012+5307 and J1713+0747—indicating the presence of constrained power at higher Fourier frequencies. The outcomes are also similar for the DM variations with PSRs J1024–0719, J1640+2224, J1713+0747, J1738+0333, J1744–1134 and J1857+0943, while significant differences appear for PSRs J0751+1807, J1012+5307, J1730–2304 and J1918+0642.

## 6.2. Discussing the results from the noise model selection

Table 2 provides the estimated amplitude (set to a reference frequency of 1yr⁻¹, and 1.4 GHz for DMv and CN₄) and slope (medians with 95% confidence intervals) for the selected components. The last column displays the Bayes factors in favor of the selected model, against a model composed of RN+DMv, with favored number of frequency bins specifically obtained for this model.

As described in the previous subsection, the selected models include all the noise components (RN, DMv and CN₄) for seven pulsars: PSRs J0030+0451, J1012+5307, J1022+1001, J1640+2224, J1713+5307, J1744–1134 and J1918–0642. For these, we observed very high Bayes factors (> 100) in favor of including CN₄ over RN+DMv, apart for PSRs J1713+0747 and J1918–0642, with only $B_{RN+DMv+CN_4}^{RN+DMv} = 16.0$ and 1.7 respectively. The additional chromatic noise is also included in the favored model for PSR J0751+1807 (DMv + CN₄), where the apparent RN in Figure 2 becomes highly unconstrained once CN₄ is included. We still noticed low Bayes factors $B_{RN+DMv}^{DMv+CN_4} = 9.4$ and $B_{RN+DMv+CN_4}^{DMv+CN_4} = 1.2$, indicating only mild support for the final model. It is similar for PSR J1024–0719, where the favored model DMv is poorly favored against RN+DMv: $B_{RN+DMv+CN_4}^{DMv} = 1.4$. This shows a lack of support for the presence of RN in the data, also being consistent for this pulsar with a similar analysis applied on the PPTA DR2 (Goncharov et al. 2020). For the three remaining pulsars, PSRs J1730–2304, J1738+0333 and J1857+0943, the selected models remain composed of RN+DMv.

### 6.2.1. Comparison with standard noise models

Here we compare our results with those obtained using the standard noise models – namely TM, WN, RN30, and DM100, as defined in Section 5. The estimation of RN (if present) and DMv parameters are highly consistent with the ones obtained with the standard noise models (cf. Figure 2) for PSRs J0030+0451, J0751+1807, J1024–0719, J1713+0747, J1730–2304, J1738+0333, J1744–1134 and J1857+0943. Mild differences are observed for PSR J1012+5307, with DMv being slightly shifted toward a lower spectral index, and PSR J1918–0642, where the RN and DMv are respectively less and more constrained. For the case of PSR J1022+1001, while DMv remains fully consistent, the RN becomes highly constrained to the peak at the low spectral index (cf. Figure 2), due to the inclusion of many frequency bins (RN98). However, we emphasize that such noise model is expected to be insufficient to properly include the power of the RN at low PSD frequencies. A dedicated study for this pulsar is still required to ensure a proper characterization of the noise. Lastly, we observed significant changes on both RN and DMv parameter posterior distributions for PSR J1640+2224, where the RN gets highly constrained at a low spectral index and high amplitude, while the DMv displays a bimodal shape, with a main peak consistent with the standard noise model and a secondary peak favoring high spectral indices ($\gamma > 7$) and low amplitudes ($\log_{10}A < -15$).

### 6.2.2. Examining the value of the chromatic index

To validate the choice of the chromatic index of CN₄, we performed an additional parameter estimation using the same model as in Table 1, but instead of fixing its chromatic index, we treat it as a free parameter $\chi$. Figure 4 displays the posterior distributions of the chromatic indices obtained for all of the pulsars that include CN₄ in the noise model and the error bars display the 99.7% credible interval of the distribution. We observed a strong consistency with the predicted value of 4 from the thin screen approximation (shown as a black horizontal line in the figure) for PSR J1012+5307, with $\chi = 4.00^{+0.8}_{-0.6}$ (median with 68% credible intervals). The results are also consistent with this predicted value for PSRs J0751+1807, J1713+0747 and J1744–1134, that all display constraints favoring higher values, and also for PSR J1640+2224, where the posterior distribution is favoring slightly lower values. Instead, the chromatic
<table><caption>**Table 2.** Amplitude and spectral index of RN, DMv, and CN₄ from DR2low.</caption><thead><tr><th rowspan="2">Pulsar</th><th colspan="2">Red noise</th><th colspan="2">DM variations</th><th colspan="2">Chromatic noise</th><th rowspan="2">B<sub>Fav</sub><sub>RN+DMv</sub></th></tr><tr><th>log<sub>10</sub> A</th><th>γ</th><th>log<sub>10</sub> A</th><th>γ</th><th>log<sub>10</sub> A</th><th>γ</th></tr></thead><tbody><tr><td>J0030+0451</td><td>-14.79<sup>+1.24</sup><sub>-2.02</sub></td><td>4.61<sup>+2.26</sup><sub>-3.56</sub></td><td>-13.49<sup>+0.10</sup><sub>-0.09</sub></td><td>0.96<sup>+0.30</sup><sub>-0.27</sub></td><td>-17.15<sup>+0.48</sup><sub>-0.77</sub></td><td>3.12<sup>+1.78</sup><sub>-1.16</sub></td><td>≥ 10<sup>3</sup></td></tr><tr><td>J0751+1807</td><td>X</td><td>X</td><td>-13.37<sup>+0.24</sup><sub>-0.64</sub></td><td>2.69<sup>+1.72</sup><sub>-0.77</sub></td><td>-15.15<sup>+0.17</sup><sub>-0.47</sub></td><td>0.74<sup>+1.57</sup><sub>-0.68</sub></td><td>9.4</td></tr><tr><td>J1012+5307</td><td>-12.92<sup>+0.14</sup><sub>-0.13</sub></td><td>1.41<sup>+0.50</sup><sub>-0.48</sub></td><td>-13.73<sup>+0.14</sup><sub>-0.12</sub></td><td>1.72<sup>+0.72</sup><sub>-0.65</sub></td><td>-16.30<sup>+0.16</sup><sub>-0.22</sub></td><td>1.21<sup>+0.53</sup><sub>-0.72</sub></td><td>≥ 10<sup>3</sup></td></tr><tr><td>J1022+1001</td><td>-13.04<sup>+0.13</sup><sub>-0.12</sub></td><td>0.15<sup>+0.39</sup><sub>-0.15</sub></td><td>-13.27<sup>+0.10</sup><sub>-0.09</sub></td><td>1.59<sup>+0.34</sup><sub>-0.32</sub></td><td>-16.52<sup>+0.20</sup><sub>-0.54</sub></td><td>1.30<sup>+1.47</sup><sub>-0.67</sub></td><td>≥ 10<sup>3</sup></td></tr><tr><td>J1024-0719</td><td>X</td><td>X</td><td>-13.74<sup>+0.33</sup><sub>-0.41</sub></td><td>3.38<sup>+1.51</sup><sub>-1.10</sub></td><td>X</td><td>X</td><td>1.4</td></tr><tr><td>J1640+2224</td><td>-13.26<sup>+0.18</sup><sub>-2.27</sub></td><td>0.28<sup>+1.49</sup><sub>-0.27</sub></td><td>-14.26<sup>+0.55</sup><sub>-1.77</sub></td><td>3.19<sup>+3.71</sup><sub>-1.83</sub></td><td>-15.93<sup>+0.19</sup><sub>-1.78</sub></td><td>1.79<sup>+4.13</sup><sub>-0.86</sub></td><td>≥ 10<sup>3</sup></td></tr><tr><td>J1713+0747</td><td>-13.87<sup>+0.16</sup><sub>-0.43</sub></td><td>2.07<sup>+1.83</sup><sub>-0.63</sub></td><td>-13.80<sup>+0.11</sup><sub>-0.12</sub></td><td>1.72<sup>+0.62</sup><sub>-0.50</sub></td><td>-15.88<sup>+0.39</sup><sub>-1.23</sub></td><td>2.30<sup>+2.53</sup><sub>-1.89</sub></td><td>16.0</td></tr><tr><td>J1730-2304</td><td>-13.83<sup>+0.63</sup><sub>-1.63</sub></td><td>3.20<sup>+3.46</sup><sub>-2.43</sub></td><td>-13.64<sup>+0.30</sup><sub>-0.78</sub></td><td>2.19<sup>+2.77</sup><sub>-1.21</sub></td><td>X</td><td>X</td><td>X</td></tr><tr><td>J1738+0333</td><td>-15.10<sup>+2.21</sup><sub>-4.60</sub></td><td>3.72<sup>+3.12</sup><sub>-3.52</sub></td><td>-13.05<sup>+0.28</sup><sub>-0.51</sub></td><td>2.45<sup>+1.52</sup><sub>-0.89</sub></td><td>X</td><td>X</td><td>X</td></tr><tr><td>J1744-1134</td><td>-13.34<sup>+0.22</sup><sub>-0.30</sub></td><td>1.05<sup>+1.20</sup><sub>-0.86</sub></td><td>-13.68<sup>+0.10</sup><sub>-0.10</sub></td><td>1.69<sup>+0.38</sup><sub>-0.32</sub></td><td>-15.90<sup>+0.17</sup><sub>-0.25</sub></td><td>0.53<sup>+0.92</sup><sub>-0.49</sub></td><td>124.0</td></tr><tr><td>J1857+0943</td><td>-15.62<sup>+1.93</sup><sub>-4.15</sub></td><td>3.72<sup>+3.09</sup><sub>-3.48</sub></td><td>-13.31<sup>+0.14</sup><sub>-0.26</sub></td><td>1.88<sup>+1.11</sup><sub>-0.53</sub></td><td>X</td><td>X</td><td>X</td></tr><tr><td>J1918-0642</td><td>-15.26<sup>+1.90</sup><sub>-4.49</sub></td><td>3.45<sup>+3.30</sup><sub>-3.21</sub></td><td>-14.48<sup>+1.21</sup><sub>-1.15</sub></td><td>5.03<sup>+1.86</sup><sub>-3.15</sub></td><td>-16.33<sup>+1.28</sup><sub>-3.29</sub></td><td>3.34<sup>+3.43</sup><sub>-2.96</sub></td><td>1.7</td></tr></tbody></table>

Notes. Median values with 95% confidence intervals are shown. The estimates are obtained using the favored noise models listed in Table 1.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Posterior distribution of the chromatic index (blue violins) and their 99.7% credible intervals (errorbars) obtained for the pulsars that favor the presence of CN₄ in their noise model. Here we use the favored noise model and replace the CN₄ noise component by a Gaussian process power-law PSD signal with the same number of frequency bins, but a free chromatic index, with a prior defined as $U(0, 10)$. The black horizontal line displays the fixed value for the case of CN₄.</figcaption></figure>

index of PSR J1918-0642 appears constrained at lower values, with a main peak centered at $\chi = 1.6^{+0.4}_{-0.4}$. However, we also observed a secondary peak centered at 4, leading to a 99.7% credible interval consistent with the predicted value. Thus, we suggest conducting a dedicated analysis of this pulsar, particularly when incorporating low radio frequency data. The two remaining pulsars, PSRs J0030+0451 and J1022+1001 provide chromatic indices that are highly inconsistent with thin screen turbulence theory, being respectively $\chi = 3.0^{+0.1}_{-0.1}$ and $0.92^{+0.09}_{-0.08}$. We emphasize the presence of strong dispersion variations caused by the solar wind for both pulsars, which would induce a reduction of the measured value in case of imperfect modeling (cf. next section) since DMv has a chromaticity of 2. Furthermore, we could also expect a lower value for PSR J1022+1001 as we pointed a bimodal distribution of RN parameters that is poorly modeled with a single power law. Thus, an imperfect treatment of RN, for which the chromaticity is equal to zero, would lead to reduction of the measured chromatic index.

## 7. Impact of solar wind

In this section, we study in more detail the impact of the SW model on PSRs J0030+0451 and J1022+1001. Specifically, we aim to investigate the differences in the noise properties when the SW signature is absent by inspecting the free spectra of DM variations and RN, where we fit for the power spectral density amplitude for each bin independently. For this reason, we analyzed two truncated datasets per pulsar by removing the DR2low data taken closer than 45° and 75° to the Sun, respectively. The 45° threshold follows Tiburzi et al. (2019), where it is shown that the SW effects are well modeled beyond this angular distance from the Sun (see their Figure 12). The more conservative 75° cut is introduced because NenuFAR— whose data were absent in Tiburzi et al. (2019)—offers an even higher sensitivity to DM effects than LOFAR. This stricter cut ensures that any residual and unmodeled SW contribution is removed, even in the presence of these highly sensitive datasets.

The adopted noise model is the same as in Table 1; however, in the previous analysis the SW was mitigated following Table A.2, while in this case it is no longer part of the model. In Figure 5 we show the resulting free-spectra. By inspecting the DM spectra of the DR2low we observed that there is a noticeable increase in power at the harmonics of $1\text{yr}^{-1}$ compared to the two cases with solar angle cuts. This suggests that the adopted SW model does not fully capture the entirety of the signal dynamics leading to a periodic leakage into the DM spectrum. Hence, the simple power-law model compensates by flattening its spectral index and increasing the number of Fourier components as shown in Section 5 and 6 for the solar pulsars such as PSRs J0030+0451, J1022+1001 and J1744-1134. This compensatory behavior highlights the limitations of the current SW modeling in accurately characterizing its induced noise. In future works one could explore some improvement of the current SW model by changing the number of Fourier bins and the power-law assumption of the PSD. We also noticed that the 45°- and 75°-cut datasets yield consistent results across almost all bins in the spectrum. In contrast, the RN does not appear to be affected by the SW since no differences in the free-spectra are detected among the three datasets.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Free spectra for DMv (referred to 1.4GHz) and RN for PSR J0030+0451 (top) and PSR J1022+1001 (bottom). Each violin represents the posterior distribution of the RMS amplitude in a given Fourier bin. The blue violins correspond to the full DR2low dataset; the green and orange violins show the results obtained after applying solar angle cuts of 45° and 75°, respectively. The vertical dashed lines indicate the harmonics of the 1yr⁻¹ frequency, which is representative of the SW effect. The number of Fourier bins displayed matches that of the preferred noise model listed in Table 1.</figcaption></figure>

## 8. Conclusions

In this work we combined high-precision timing data from the EPTA DR2new+ with low-frequency observations from LOFAR and NenuFAR. Among the 25 MSPs included in the DR2new+, we identified respectively 12 and 2 sources in common with the LOFAR and NenuFAR datasets. The time span of the pulsars in common almost perfectly overlaps with that of DR2new+. The addition of low-frequency radio data significantly extends the frequency coverage of the dataset, enhancing our sensitivity to propagation effects such as DM variations and higher order chromatic delays. This new dataset is referred to as DR2low. We performed multiple iterations of the noise analysis, following the approach outlined in EPTA Collaboration et al. (2023c), and we compared the results obtained with DR2low against those from DR2new+. As a first step, we adopted a standard noise model consisting of the TM, WN, RN30, and DMv100 components. In addition, we applied SW corrections to those pulsars with low ecliptic latitudes and sufficient observing cadence around solar conjunctions (see Table A.2). The corner plots for the standard noise model are presented in Figure 2.

Overall, the results show a significant improvement in constraining the DM variations, as expected with the increased low-frequency coverage of the dataset. Moreover, this enhanced sensitivity allows us to more clearly distinguish between DM variations and RN, or at least to reduce the correlation between the

two processes. This improvement is evident in Figure 3, where the typical L-shaped correlation between RN and DMv has been mitigated by the inclusion of low-frequency data. This result is consistent with the simulations presented by Ferranti, I. et al. (2025). Furthermore, we showed an improved consistency with the results from NANOGrav 15y and PPTA DR3.

We then performed a more in-depth analysis to find the favored noise model for DR2low. We observed a significant change in the final noise model with respect to DR2new+ (see Table 1). The extended time span and high cadence of the low-frequency data (approximately 10 years) now allow us to better sample signals across the power spectral domain. This improvement enables us to accurately redistribute power into the correct noise components. As a result, we found that ten pulsars have significant RN in the final model in contrast to the five in DR2new+. Moreover, DMv are present in all the 12 pulsars, while for DR2new+ this component was missing for PSR J0030+0451. We also found an additional noise component with a fixed chromatic index of 4 in the final models of 8 out of 12 pulsars, whereas none were detected in DR2new+. To validate the chosen chromatic index, we fit for this parameter ($\chi$) in an additional analysis. Five out of eight pulsars show consistency with the expected value from a thin screen with same-size inhomogeneities, $\chi = 4$. A more detailed analysis is required for PSR J1918−0642, which displays a main peak at $\chi < 2$, but also a secondary peak that is
consistent with $\chi = 4$. For PSRs J0030+0451 and J1022+1001 the chromatic index is inconsistent with 4, possibly due to the strong DM variations caused by the solar wind.

Finally, we conducted a detailed study on the impact of the SW model for PSR J0030+0451 and PSR J1022+1001. We applied two thresholds to remove epochs from DR2low with solar angles smaller than 45° and 75°, respectively. We analyzed the datasets using the same final noise model shown in Table 1, while excluding the SW model in the truncated datasets. Figure 5 shows the free spectra of DMv and RN. We observed clear excess power in the DMv around the harmonics of the $1\text{yr}^{-1}$ frequency. This suggests that the current SW model does not fully capture all the features of this process, causing the excess power to be absorbed into the DMv modeling.

The inclusion of low-frequency data has proven to significantly improve the noise modeling in PTAs. These data help disentangle noise components and offer enhanced sensitivity to DM variations and other chromatic signals. Moreover, the addition of low-frequency observations has been shown to improve the accuracy of astrometric parameters, as demonstrated in Liu et al. (2025). This initial integration of LOFAR and NenuFAR data is a crucial step, as both datasets will be included in the upcoming IPTA DR3 (Good & International Pulsar Timing Array Team 2023), which is set to become the most extensive PTA dataset ever created, with the widest frequency coverage. Furthermore, LOFAR and NenuFAR will also be part of the next EPTA Data Release 3, which will feature a larger sample of pulsars in common between the datasets. In this context, the present analysis lays important groundwork for future studies. In particular, it will be especially interesting to assess the effect of these very low-frequency data on the detection and characterization of a GWB in upcoming data releases. The extended frequency coverage provided by LOFAR and NenuFAR has the potential to improve sensitivity to the GWB and will be an important aspect to explore in future work.

Finally, it is important to highlight the significance of future low-frequency observations, particularly with SKA-low. The expected improvements in sensitivity and time coverage provided by SKA-low will be indispensable for advancing pulsar timing array studies and for ensuring the continued success of gravitational wave detection efforts.

**Acknowledgements.** FI is supported by the University of Cagliari (IT). FI, CT, AP are supported by the Istituto Nazionale di Astrofisica. JPWV acknowledges support from NSF AccelNet award No. 2114721. SCS acknowledges the support of a College of Science and Engineering University of Galway Postgraduate Scholarship. GMS acknowledge financial support provided under the European Union's H2020 ERC Consolidator Grant "Binary Massive Black Hole Astrophysics" (B Massive, Grant Agreement: 818691). J.A. acknowledges support from the European Commission under project ARGOS-CDS (Grant Agreement number: 101094354) MB acknowledges the support from the Department of Atomic Energy, Government of India through Apex-I Project - Advance Research and Education in Mathematical Sciences at IMSc. This paper is partially based on data obtained using the NenuFAR radio-telescope. The development of NenuFAR has been supported by personnel and funding from: Observatoire Radioastronomique de Nançay, CNRS-INSU, Observatoire de Paris-PSL, Université d'Orléans, Observatoire des Sciences de l'Univers en Région Centre, Région Centre-Val de Loire, DIM-ACAV and DIM-ACAV + of Région Ile-de-France, Agence Nationale de la Recherche. We acknowledge the use of the Nançay Data Centre computing facility (CDN - Centre de Données de Nançay). The CDN is hosted by the Observatoire Radioastronomique de Nançay (ORN) in partnership with Observatoire de Paris, Université d'Orléans, OSUC, and the CNRS. The CDN is supported by the Région Centre-Val de Loire, département du Cher. The Nançay Radio Observatory (ORN) is operated by Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS) and Université d'Orléans. LOFAR, the Low-Frequency Array designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a

joint scientific policy. BCJ acknowledges the support from the Raja Ramanna Chair fellowship of the Department of Atomic Energy, Government of India (grant 3/3401 Atomic Energy Research 00 004 Research and Development 27 02 31 1002//2/2023/RRC/R&D-II/13886) and the support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI 4002 and under project No. 12-R&D-TFR-5.02-0700. IK acknowledges the support of Collège de France by means of "PAUSE -Solidarité Ukraine" program and he would like to thank the Paris Observatory and the CNRS (LPC2E) for being great host organizations for the PAUSE program. This work is supported by 100101 Key Laboratory of Radio Astronomy and Technology (Chinese Academy of Science). AC and APa acknowledges financial support from the European Research Council (ERC) starting grant 'GIGA' (grant agreement number: 101116134) and through the NWO-I Veni fellowship. DJS acknowledges support by the project "NRW-Cluster for data intensive radio astronomy: Big Bang to Big Data (B3D)" funded through the programme "Profilbildung 2020", an initiative of the Ministry of Culture and Science of the State of North Rhine-Westphalia. KT is partially supported by JSPS KAKENHI Grant Numbers 20H00180, 23K20868, and 21H04467, and Bilateral Joint Research Projects of JSPS. JW acknowledges support by the BMBF Verbundforschung under the grant 05A23PC2. JS acknowledges the support from the University of Cape Town Vice Chancellor's Future Leaders 2030 Awards programme and the South African Research Chairs Initiative of the Department of Science and Technology and the National Research Foundation. DD acknowledges the Department of Atomic Energy, Government of India's support through 'Apex Project-Advance Research and Education in Mathematical Sciences' at The Institute of Mathematical Sciences. PR acknowledges the financial assistance of the South African Radio Astronomy Observatory (SARAO) towards this research (www.sarao.ac.za).

## References

* Agazie, G., Alam, M. F., Anumarlapudi, A., et al. 2023, The Astrophysical Journal Letters, 951, L9

* Agazie, G., Antoniadis, J., Anumarlapudi, A., et al. 2024, The Astrophysical Journal, 966, 105

* Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023, ApJ, 951, L8

* Allen, B., Dhurandhar, S., Gupta, Y., et al. 2023, arXiv e-prints, arXiv:2304.04767

* Bhat, N. D. R., Cordes, J. M., Camilo, F., Nice, D. J., & Lorimer, D. R. 2004, ApJ, 605, 759

* Bilous, A. V., Kondratiev, V. I., Kramer, M., et al. 2016, A&A, 591, A134

* Bondonneau, L., Grießmeier, J. M., Theureau, G., et al. 2021, A&A, 652, A34

* Chalumeau, A., Babak, S., Petiteau, A., et al. 2022, MNRAS, 509, 5538

* Coles, W. A., Kerr, M., Shannon, R. M., et al. 2015, ApJ, 808, 113

* Cordes, J. M. & Shannon, R. M. 2010, arXiv e-prints, arXiv:1010.3785

* Cordes, J. M., Shannon, R. M., & Stinebring, D. R. 2016, The Astrophysical Journal, 817, 16

* Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3341

* Donner, J. Y. 2022, PhD thesis, Universität Bielefeld, Bielefeld (Germany)

* Donner, J. Y., Verbiest, J. P. W., Tiburzi, C., et al. 2020, A&A, 644, A153

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

* Ellis, J. & van Haasteren, R. 2017, jellis18/PTMCMCSampler: Official Release

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

* EPTA Collaboration, Antoniadis, J., Babak, S., et al. 2023a, A&A, 678, A48

* EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023b, A&A, 678, A50

* EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023c, A&A, 678, A49

* Ferranti, I., Falxa, M., Sesana, A., et al. 2025, A&A, 694, A38

* Goncharov, B., Reardon, D. J., Shannon, R. M., et al. 2020, Monthly Notices of the Royal Astronomical Society, 502, 478

* Good, D. & International Pulsar Timing Array Team. 2023, in American Astronomical Society Meeting Abstracts, Vol. 241, American Astronomical Society Meeting Abstracts, 438.02

* Hansen, B. M. S. & Phinney, E. S. 1997, MNRAS, 291, 569

* Hazboun, J. S., Simon, J., Madison, D. R., et al. 2022, ApJ, 929, 39

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

* Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974

* Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655

* Iraci, F., Chalumeau, A., Tiburzi, C., et al. 2024, A&A, 692, A170

* Joshi, B. C., Arumugasamy, P., Bagchi, M., et al. 2018, Journal of Astrophysics and Astronomy, 39, 51

* Kapil, V., Mandel, I., Berti, E., & Müller, B. 2023, MNRAS, 519, 5893

* Kondratiev, V. I., Verbiest, J. P. W., Hessels, J. W. T., et al. 2016, A&A, 585, A128
Lam, M. T., Ellis, J. A., Grillo, G., et al. 2018, ApJ, 861, 132

Larsen, B., Mingarelli, C. M. F., Hazboun, J. S., et al. 2024, ApJ, 972, 49

Lee, K. J. 2016, in Astronomical Society of the Pacific Conference Series, Vol. 502, Frontiers in Radio Astronomy and FAST Early Sciences Symposium 2015, ed. L. Qain & D. Li, 19

Lentati, L., Alexander, P., Hobson, M. P., et al. 2014, MNRAS, 437, 3004

Liu, K., Karuppusamy, R., Lee, K. J., et al. 2015, Monthly Notices of the Royal Astronomical Society, 449, 1158

Liu, K., Keane, E. F., Lee, K. J., et al. 2012, MNRAS, 420, 361

Liu, K., Parthasarathy, A., Keith, M., et al. 2025, MNRAS, 536, 2603

Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017

Miles, M. T., Shannon, R. M., Bailes, M., et al. 2023, MNRAS, 519, 3976

Miles, M. T., Shannon, R. M., Reardon, D. J., et al. 2025, MNRAS, 536, 1489

Padmanabh, P. V., Barr, E. D., Champion, D. J., et al. 2020, Monthly Notices of the Royal Astronomical Society, 500, 1178

Phinney, E. S. 2001, arXiv e-prints, astro

Ransom, S., Brazier, A., Chatterjee, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 195

Rasmussen, C. & Williams, C. 2006, Gaussian Processes for Machine Learning, Adaptative computation and machine learning series (University Press Group Limited)

Raveri, M. & Doux, C. 2021, Physical Review D, 104, 043504

Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6

Sesana, A. 2013, MNRAS, 433, L1

Shannon, R. M. & Cordes, J. M. 2012, ApJ, 761, 64

Shannon, R. M. & Cordes, J. M. 2016, Monthly Notices of the Royal Astronomical Society, 464, 2075

Siemens, X., Ellis, J., Jenet, F., & Romano, J. D. 2013, Classical and Quantum Gravity, 30, 224015

Speri, L., Porayko, N. K., Falxa, M., et al. 2022, Monthly Notices of the Royal Astronomical Society, 518, 1802

Susarla, S. C., Chalumeau, A., Tiburzi, C., et al. 2024, A&A, 692, A18

Taylor, J. H. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 117

Taylor, S. R., Baker, P. T., Hazboun, J. S., Simon, J., & Vigeland, S. J. 2021, enterprise_extensions, v2.4.3

Tiburzi, C., Shaifullah, G. M., Bassa, C. G., et al. 2021, A&A, 647, A84

Tiburzi, C., Verbiest, J. P. W., Shaifullah, G. M., et al. 2019, MNRAS, 487, 394

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

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

Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267

Verbiest, J. P. W., Vigeland, S. J., Porayko, N. K., Chen, S., & Reardon, D. J. 2024, Results in Physics, 61, 107719

Xu, H., Chen, S., Guo, Y., et al. 2023, Research in Astronomy and Astrophysics, 23, 075024

Xu, H., Huang, Y. X., Burgay, M., et al. 2021, The Astronomer's Telegram, 14642, 1

Zic, A., Reardon, D. J., Kapur, A., et al. 2023, PASA, 40, e049

14 Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str 1, D-85741 Garching, Germany

15 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China

16 LIRA, Observatoire de Paris, CNRS, PSL, Sorbonne U., U. Paris Cité, Meudon, France

17 ORN, Observatoire de Paris, CNRS, PSL, U. Orléans, Nançay, France

18 Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, P. R. China

19 Key Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, Beijing 100101, P. R. China

20 LPC2E, OSUC, Univ Orleans, CNRS, CNES, Observatoire de Paris, F-45071 Orleans, France

21 Université Paris Cité and Université Paris Saclay, CEA, CNRS, AIM, 91190 Gif-sur-Yvette, France

22 Thüringer Landessternwarte, Sternwarte 5, 07778 Tautenburg, Germany

23 Department of Physics, IIT Hyderabad, Kandi Telangana 502284, India

24 Department of Astrophysics/IMAPP, Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

25 Indian Institute of Technology, Roorkee 247667 India

26 Faculty of Advanced Science and Technology, Kumamoto University, 2-39-1 Kurokami, Kumamoto 860-8555, Japan

27 Institute of Radio Astronomy of the National Academy of Sciences of Ukraine, Mystetstv St. 4, Kharkiv, Ukraine 61002

28 High Energy Physics, Cosmology & Astrophysics Theory (HEP-CAT) Group, Department of Mathematics and Applied Mathematics, University of Cape Town, Cape Town 7700, South Africa

29 Department of Astronomy, University of Cape Town, Cape Town 7700, South Africa

30 Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, Colaba, Mumbai 400005

31 Anton Pannekoeck Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

32 Fakultät für Physik, Universität Bielefeld, Postfach 100131, 33501 Bielefeld, Germany

33 Department of Physics, GLA University, Mathura 281406, India

34 International Research Organization for Advanced Science and Technology, Kumamoto University, 2-39-1 Kurokami, Kumamoto 860-8555, Japan

35 Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany

36 Ruhr-Universität Bochum, Fakultät für Physik und Astronomie, Astronomisches Institut (AIRUB), 44801 Bochum, Germany

<footnote n="1">¹ Dipartimento di Fisica, Università di Cagliari, Cittadella Universitaria, I-09042 Monserrato (CA), Italy</footnote>

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

<footnote n="3">³ ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands</footnote>

<footnote n="4">⁴ Florida Space Institute, University of Central Florida, 12354 Research Parkway, Orlando, FL 32826, USA</footnote>

<footnote n="5">⁵ Physics, School of Natural Sciences, Ollscoil na Gaillimhe — University of Galway, University Road, Galway, H91 TK33, Ireland</footnote>

<footnote n="6">⁶ National Centre for Radio Astrophysics, Pune University Campus, Pune 411007, India</footnote>

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

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

<footnote n="9">⁹ Institute of Astrophysics, Foundation for Research & Technology – Hellas (FORTH), GR-70013 Heraklion, Greece</footnote>

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

<footnote n="11">¹¹ The Institute of Mathematical Sciences, C. I. T. Campus, Taramani, Chennai 600113, India</footnote>

<footnote n="12">¹² School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 PN40, Ireland</footnote>

<footnote n="13">¹³ Homi Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai 400094, India</footnote>
# Appendix A: Parameters of the noise models used in this work

<table><caption>**Table A.1.** Noise components and corresponding abbreviations.</caption><thead><tr><th>Noise component (abrev.)</th><th>Parameters</th><th>Priors</th></tr></thead><tbody><tr><td rowspan="3">White noise (WN)</td><td>EFAC</td><td>U(0.01, 10)</td></tr><tr><td>EQUAD</td><td>log<sub>10</sub> U(10<sup>-9</sup>, 10<sup>-5</sup>)</td></tr><tr><td>ECORR</td><td>log<sub>10</sub> U(10<sup>-9</sup>, 10<sup>-2</sup>)</td></tr><tr><td rowspan="2">Achromatic red noise (RN)</td><td>A<sub>RN</sub></td><td>log<sub>10</sub> U(10<sup>-20</sup>, 10<sup>-11</sup>)</td></tr><tr><td>γ<sub>RN</sub></td><td>U(0, 7)</td></tr><tr><td rowspan="2">DM variations (DMv)</td><td>A<sub>DMv</sub></td><td>log<sub>10</sub> U(10<sup>-20</sup>, 10<sup>-11</sup>)</td></tr><tr><td>γ<sub>DMv</sub></td><td>U(0, 7)</td></tr><tr><td rowspan="2">Chromatic noise (CN<sub>4</sub>)</td><td>A<sub>CN<sub>4</sub></sub></td><td>log<sub>10</sub> U(10<sup>-20</sup>, 10<sup>-11</sup>)</td></tr><tr><td>γ<sub>CN<sub>4</sub></sub></td><td>U(0, 7)</td></tr><tr><td rowspan="4">Exponential dip (E)</td><td>A<sub>E</sub></td><td>log<sub>10</sub> U(10<sup>-10</sup>, 10<sup>-2</sup>)</td></tr><tr><td>τ<sub>E</sub></td><td>log<sub>10</sub> U(10<sup>0</sup>, 10<sup>2.5</sup>)</td></tr><tr><td>t<sub>0,E</sub></td><td>U(57490.0, 57530.0)</td></tr><tr><td>χ<sub>E</sub></td><td>U(0, 10)</td></tr><tr><td rowspan="3">Solar wind dispersion (SW)</td><td>n<sub>earth</sub></td><td>U(0, 30)</td></tr><tr><td>A<sub>SW</sub></td><td>log<sub>10</sub> U(10<sup>-10</sup>, 10<sup>1</sup>)</td></tr><tr><td>γ<sub>SW</sub></td><td>U(-2, -1)</td></tr></tbody></table>

**Notes.** The table lists the noise components used in this work with their abbreviations (first column), parameter names (second column), and prior definitions (third column). $U(a, b)$ denotes a uniform prior between $a$ and $b$.

<table><caption>**Table A.2.** Modeling approach for dispersion caused by solar winds.</caption><thead><tr><th rowspan="2">Pulsar</th><th rowspan="2">ELAT (deg.)</th><th colspan="2">SW modeling – DR2low</th><th colspan="2">SW modeling – DR2new+</th></tr><tr><th>Timing - NE_SW</th><th>Noise</th><th>Timing - NE_SW</th><th>Noise</th></tr></thead><tbody><tr><td>J1022+1001</td><td>-0.06</td><td>9.51 ± 0.03 (fitted)</td><td>n_earth + SWv</td><td>7.2 ± 0.5 (fitted)</td><td>marg.</td></tr><tr><td>J1730-2304</td><td>0.19</td><td>9.9 ± 0.4 (fitted)</td><td>n_earth</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J0030+0451</td><td>1.45</td><td>8.20 ± 0.03 (fitted)</td><td>n_earth + SWv</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J0751+1807</td><td>-2.81</td><td>7.9 (fixed)</td><td>n_earth + SWv</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1744-1134</td><td>11.81</td><td>3.62 ± 0.08 (fitted)</td><td>n_earth + SWv</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1918-0642</td><td>15.35</td><td>4 (fixed)</td><td>marg.</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1024-0719</td><td>-16.04</td><td>4 (fixed)</td><td>marg.</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1738+0333</td><td>26.88</td><td>4 (fixed)</td><td>marg.</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1713+0747</td><td>30.70</td><td>4 (fixed)</td><td>marg.</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1857+0943</td><td>32.32</td><td>4 (fixed)</td><td>marg.</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1012+5307</td><td>38.76</td><td>4 (fixed)</td><td>marg.</td><td>7.9 (fixed)</td><td>marg.</td></tr><tr><td>J1640+2224</td><td>44.06</td><td>4 (fixed)</td><td>marg.</td><td>7.9 (fixed)</td><td>marg.</td></tr></tbody></table>

**Notes.** The modeling strategy follows Susarla et al. (2024). The pulsar names and ecliptic latitudes (in degrees) are listed in the first and second columns. For both DR2low (third major column) and DR2new+ (last major column), the table shows the adopted methods for timing and noise analyses. In the timing analysis, the constant term is either fixed or fitted with TEMPO2. In the noise analysis, the constant term is included by either marginalizing over its error (marg.; see Section 4) or fitting for it ($n_{\text{earth}}$). The time-varying component (SWv) represents fluctuations in the electron density at 1 AU, modeled as a power-law PSD parameterized by its amplitude $A_{\text{SW}}$, referred to $f = 1 \text{ yr}^{-1}$, and spectral index $\gamma_{\text{SW}}$.
Appendix B: Estimation of the ECORR parameters

<figure><img src="image_7.png" /><figcaption>**Fig. B.1.** Posterior distribution of the ECORR parameters (logarithmic scale) and their 99.7% credible intervals (errorbars) obtained for the LOFAR systems for each pulsar.</figcaption></figure>
Appendix C: Solar wind parameters

<figure><img src="image_8.png" /><figcaption>Fig. C.1. Posterior distribution of the SW parameters for PSRs J0030+0451, J0751+1807, J1022+1001, and J1744-1134.</figcaption></figure>
Appendix D: Timing residuals and observing radio frequencies of DR2low

<figure><img src="image_9.png" /><figcaption>**Fig. D.1.** Time in year and MJD (resp. for the top and bottom x-axes) vs. observing frequency (in MHz) shown in logarithmic scale for the 12 pulsars considered in this work displayed in green, orange, and blue respectively for the DR2new+, LOFAR, and NenuFAR data. The gray vertical areas emphasize the periods when the pulsar-sun angles are equal to or lower than 20 degrees.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. D.2.** Same as Fig. D.1, but displaying the post-fit timing residuals (in µs) for the y-axis.</figcaption></figure>