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). For this work we incorporated for the first time low-frequency observations from LOFAR (\(100-200\,\mathrm{MHz}\)) and NenuFAR (\(30-90\,\mathrm{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 \(\sim\)​11 years, with radio frequencies spanning the range \(30-2500\,\mathrm{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\(\sigma\) detection of the GWB. 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\(_4\)). Finally, we carried out a detailed analysis on the choice of the chromatic index for CN\(_4\) and the contribution of the solar wind. 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\(_4\) 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.

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 [1]. In the most simple case of circular SMBHBs, this noise has a red power spectrum, characterized by a steep spectral index of -13/3 [2]. Moreover, these GW-induced fluctuations follow a spatial correlation among pulsars, referred to as the Hellings & Downs curve [3]. To date, eight PTA collaborations exist in the world: the European PTA [4], the North-American Nanohertz Observatory for Gravitational waves [5], the Parkes PTA [6], the Indian PTA [7], the Chinese PTA [8], the MeerKAT PTA [9], the African Pulsar Timing (APT³) consortium, and the International PTA [10]. In 2023, five PTA collaborations presented results supporting the very first evidence for the presence of a GWB in PTA data [11]–[14], followed a year later by the sixth collaboration [15]. Although none of these studies proved the GWB presence with a significance higher than the required 5\(\sigma\) threshold [16], 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 [17] and [18]), 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}{\nu} \right )^2},\] where \(f_p\) is the plasma frequency and \(\nu\) the observing (radio) frequency. The frequency-dependent delay \(\delta t\) can be expressed as \[\delta t \approx \mathcal{D} \frac{\mathrm{DM}}{\nu^2},\] where \(\mathcal{D}\) is the dispersion constant, and DM the dispersion measure (DM), expressed as \[DM = \int_{\mathrm{LoS}} n_e dl ,\] 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 \nu^{-\gamma_{s}},\] where \(\nu\) 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 [19].

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 [20]–[22]. 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 [23], [24].

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 [25] and the New extension in Nançay upgrading LOFAR [26], 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 [17] and [18] 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 [4], 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 [27]. These pulsars and their characteristics are reported in Tables B1-B7 of [17]. 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 [28], and refined using the TEMPONEST Bayesian toolkit [29] 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 evidence for a GWB with a statistical significance exceeding \(3.5\sigma\) [11]. The importance of radio-frequency coverage in the EPTA data was further investigated in [30], 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 [25] 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 [31]–[34]. 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 [33] and [35]. 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 [36] but with the uncertainties calculated through a Monte-Carlo routine as recommended by [10]. Based on the work by [33], 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 [37].

2.3 NenuFAR↩︎

NenuFAR [26] 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\(^\circ\) 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 [38]. 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 consistent with their initial values, likely due to the larger ToA uncertainties in the low-frequency data.

In Figures 8 and 9, 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⁴.

4 Noise analysis: Methodology↩︎

4.1 Framework and tools↩︎

Having the combined dataset, we performed various iterations of noise analysis [18], 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⁵ to extract data products from TEMPO2. The software Enterprise [39] and its complementary module Enterprise_Extensions [40] are used to model the noise components and build the likelihood and prior functions, applied to the MCMC sampler PTMCMCSampler [41].

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

\[\rm{log} \;L(\vec{\delta t}|\vec{\theta},\vec{\theta'}) = -\frac{1}{2} \left( \left[ \vec{\delta t} - \vec{y}(\vec{\theta'}) \right]^T \;C^{-1}(\vec{\theta}) \;\left[ \vec{\delta t} - y(\vec{\theta'})\right] \right) + \rm{cst.},\] where the timing residuals \(\vec{\delta t}\) are reduced by any deterministic time-domain waveform \(\vec{y}\) with parameters \(\vec{\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(\vec{\theta}) = \sum_\lambda C^\lambda(\vec{\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 [18], 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 [42], [43] to include its constant and time-varying components that are fitted simultaneously with other noise components during the noise analysis. Table ([tab:parprior]) 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 [18], with an Error FACtor (EFAC) applied to the provided ToA uncertainty \(\sigma_{\rm{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 [44]–[46] may affect ECORR value. The white noise terms are added to the likelihood covariance matrix \(C\)

\[\label{eq:wn} \begin{align} &C^{\rm{WN}}_{\alpha} = \\ &\begin{cases} \left( \texttt{EF}_{\alpha}^2 \times \sigma_{\rm{init, i}}^2 + \texttt{EQ}_{\alpha}^2 \right) \delta_{ij} \;\delta_{\alpha\beta} \text{, if \alpha \notin \{LOFAR, NenuFAR\},}\\ \left[ \left( \texttt{EF}_{\alpha}^2 \times \sigma_{\rm{init, i}}^2 + \texttt{EQ}_{\alpha}^2 \right) \;\delta_{ij} + \texttt{EC}_\alpha^2 \;\delta_{\rm{E}(i)\rm{E}(j)} \right] \;\delta_{\alpha\beta} \text{, else,} \end{cases} \end{align}\tag{1}\] 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 \(\rm{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 [47], using the covariance-function expansion to express the likelihood component of the noise \(\lambda\) as \[\label{eq:GP} C^{\lambda,\rm{GP}}(t_i,t_j) = \sum_{\mu,\nu} \;F^\lambda_\mu(t) \;\Phi^\lambda_{\mu\nu} F^\lambda_\nu(t'),\tag{2}\] where \(F\) is a set of basis functions, \(\Phi\) is the covariance matrix of the Gaussian process in the weight-space view [48], 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 \(\nu^{-4}\) frequency dependence, consistent with the scattering expected under the thin screen approximation with same-size inhomogeneities [44].⁶ These three noise components (respectively labeled RN, DMv and CN\(_4\)) 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

\[\label{eq:GPPSD} \Phi^\lambda_{\mu\nu} = S^\lambda(f_\mu) \; \Delta f \;\delta_{\mu\nu},\tag{3}\] 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: \[\label{mgdfwylk} S^\lambda(f) = \frac{A_\lambda^2}{12\pi^2} \left( \frac{f}{f_{\rm{yr}}}\right)^{-\gamma_\lambda} f_{\rm{yr}}^{-3}.\tag{4}\] Here \(f_{\rm{yr}}\) refers to the frequency of \(1 \rm{yr}^{-1}\) and \(A_\lambda\), \(\gamma_\lambda\) are respectively the amplitude referred at \(f_{\rm{yr}}\) and the spectral index of the noise component \(\lambda\).

In our case, the different noise components are \(\lambda = \{\texttt{RN}, \texttt{DMv}, \texttt{CN_4}\}\), which differ from their chromaticity (i.e., dependence on the observed radio-frequency \(\nu\)) that is added in the basis functions as \(F^\lambda_\mu(t) \propto (\nu / \nu_{\rm{ref}})^{-\chi}\), where \(\chi = 0, 2, 4\) for RN, DMv and CN\(_4\) respectively, and \(\nu_{\rm{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. 2 ), 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 [50].

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 [43] which contains a deterministic and a stochastic component.

We include a deterministic signal which account for the electron content \(n_e\) as a function of the angular separation between the pulsar and the Sun sampling it from \[DM_{\mathrm{SW}} = n_e\frac{\rho}{r_e\sin{\rho}}\left[1~\mathrm{AU}\right]^2,\] where \(\rho\) is the pulsar-Sun-observer angle, \(n_e\) is the electron density at 1\(\mathrm{AU}\) (denoted as NE_SW in TEMPO2), and \(r_e\) is the observatory-Sun distance [37], [51]. 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_{\mathrm{SW}} = A_{\mathrm{SW}}^2\left( \frac{f}{\mathrm{yr}^{-1}} \right)^{-\gamma_{\mathrm{SW}}} \mathrm{yr}^3.\] Here \(S_{\mathrm{SW}}\) is the power spectral density for SWGP, \(A_{\mathrm{SW}}\) is the SW amplitude at a frequency of \(1/1\mathrm{yr}\), \(\gamma_{\mathrm{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.5yr\(^{-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 [tab:SWmodelling] we summarize how we modeled the SW for each pulsar, and provide a comparison with the model used in the EPTA DR2.

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 [52]). 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 [43] recommend including only the deterministic component of the SW model for pulsars with ELAT below \(3^\circ\), 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^\circ\), and to PSR J0751+1807, which was not included in [43]. 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 [53]–[55]: 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 [54], 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 (\(\sim\) MJD \(57510\)), where [56] 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 [50]. As in [18], we model this effect as a deterministic signal with a waveform defined as

\[\label{kjpxmgre} \begin{align} &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 \;\rm{GHz})^{-\chi_E} \;\rm{exp}\left( -\frac{t - t_{0,E}}{\tau_E}\right) \text{, if \;t \geq t_{0,E}} \end{cases}, \end{align}\tag{5}\] where \(A^E\), \(\chi^E\), \(t_0^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 [18]. 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 [18], where we defined a set of candidate noise models, fit for the number of frequency bins (i.e., \(\nu\) in Eq. 2 ) 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 [18].

In all runs, we included the SW model described in Table [tab:SWmodelling], 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⁷ of the posteriors distributions for the same parameters obtained in [57] on the NANOGrav 15-year dataset [58] and PPTA DR3 [59] 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).

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 \(\sim 0.77\sigma\) according to the description of a tension metric described in [18] based on [60]. 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 \(\sim\)​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\(^\circ\)—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 [61].

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 \(\leq 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_{\rm{RN}} \sim 1.1\)) and high amplitudes (\(\rm{log}_{10} A_{\rm{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 [30] 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_{\rm{RN}} \sim 0\)), while the second one favors a steeper power law (\(\gamma_{\rm{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 [62], [63], 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 [18].

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 [64], 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 (\(\sim3\) 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\sigma\). We observed mild tensions for PSR J1713+0747 (\(1.85\sigma\) and \(1.50\sigma\) 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 (\(\sim 0.43\sigma\)) but a moderate tension with NANOGrav (\(\sim 2.18\sigma\)) which favors a higher spectral index and a lower amplitude as obtained in DR2new+. Inconsistencies between PTAs were already found in [57] 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 [18]. The posterior distributions of ECORR parameters for every pulsars/systems are shown in Figure 6, 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 7. In this section, we first intend to provide a comparison with the results that were obtained with DR2new+, then interpret the results for DR2low.

6.1 Comparing the favored models of DR2new+ and DR2low↩︎

Table [tab:favmod] presents the final noise models obtained with DR2low (left) from this work, and DR2new+ (right) from the work in [18]. It explicitly describes the noise components (RN, DMv and CN\(_4\)) 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 [11], 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 improvement 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 domain. 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 pulsars apart from PSR J0030+0451 with DR2new+, they are now present for all of them with DR2low. This property is not surprising, however, thanks to the high sensitivity of low-frequency data toward DM variations.

For the selected number of frequency bins, we observed comparable results for the RN component, with relatively low values (\(\sim\)​10–22) for PSRs J0030+0451 and J1744\(-\)​1134, and significantly higher values (\(\sim\)​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 [tab:favparvals] provides the estimated amplitude (set to a reference frequency of 1\(\rm{yr}^{-1}\), and \(1.4\) GHz for DMv and CN\(_4\)) 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\(_4\)) 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\(_4\) over RN+DMv, apart for PSRs J1713+0747 and J1918\(-\)​0642, with only \(\mathcal{B}^{\rm{RN+DMv+CN_4}}_{\rm{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\(_4\)), where the apparent RN in Figure 2 becomes highly unconstrained once CN\(_4\) is included. We still noticed low Bayes factors \(\mathcal{B}^{\rm{DMv+CN_4}}_{\rm{RN+DMv}} = 9.4\) and \(\mathcal{B}^{\rm{DMv+CN_4}}_{\rm{RN+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: \(\mathcal{B}^{\rm{DMv}}_{\rm{RN+DMv+CN_4}} = 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 [56]. 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 (\(\rm{log}_{10} \rm A < -15\)).

6.2.2 Examining the value of the chromatic index↩︎

To validate the choice of the chromatic index of CN\(_4\), we performed an additional parameter estimation using the same model as in Table [tab:favmod], 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\(_4\) 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.6}^{+0.8}\) (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 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.08}^{+0.09}\). 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\(^\circ\) and 75\(^\circ\) to the Sun, respectively. The 45\(^\circ\) threshold follows [37], 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\(^\circ\) cut is introduced because NenuFAR— whose data were absent in [37]—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 [tab:favmod]; however, in the previous analysis the SW was mitigated following Table [tab:SWmodelling], 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\mathrm{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^\circ\)- and \(75^\circ\)-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.

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 [18], 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 [tab:SWmodelling]). 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 [30]. 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 [tab:favmod]). 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\(^\circ\) and 75\(^\circ\), respectively. We analyzed the datasets using the same final noise model shown in Table [tab:favmod], 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\mathrm{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 [52]. This initial integration of LOFAR and NenuFAR data is a crucial step, as both datasets will be included in the upcoming IPTA DR3 [65], 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.

9 Parameters of the noise models used in this work↩︎

column1-Z = c, cell21 = r=3c, cell51 = r=2c, cell71 = r=2c, cell91 = r=2c, cell111 = r=4c, cell151 = r=3c, vline1,2,3,4,5,6,7 = -, hline1,2,5,7,9,11,15,18 = -, Noise component (abrev.) & Parameters & Priors
White noise (WN) & EFAC & \(\mathcal{U}(0.01,10)\)
& EQUAD & \(\rm{log}_{10} \;\mathcal{U}(10^{-9},10^{-5})\)
& ECORR & \(\rm{log}_{10} \;\mathcal{U}(10^{-9},10^{-2})\)
Achromatic red noise (RN) & \(A_{\rm{RN}}\) & \(\rm{log}_{10} \;\mathcal{U}(10^{-20},10^{-11})\)
& \(\gamma_{\rm{RN}}\) & \(\mathcal{U}(0,7)\)
DM variations (DMv) & \(A_{\rm{DMv}}\) & \(\rm{log}_{10} \;\mathcal{U}(10^{-20},10^{-11})\)
& \(\gamma_{\rm{DMv}}\) & \(\mathcal{U}(0,7)\)
Chromatic noise (CN\(_4\)) & \(A_{\rm{CN_4}}\) & \(\rm{log}_{10} \;\mathcal{U}(10^{-20},10^{-11})\)
& \(\gamma_{\rm{CN_4}}\) & \(\mathcal{U}(0,7)\)
Exponential dip (E) & \(A_E\) & \(\rm{log}_{10} \;\mathcal{U}(10^{-10},10^{-2})\)
& \(\tau_E\) & \(\rm{log}_{10} \;\mathcal{U}(10^0,10^{2.5})\)
& \(t_{0,E}\) & \(\mathcal{U}(57490.0,57530.0)\)
& \(\chi_E\) & \(\mathcal{U}(0,10)\)
Solar wind dispersion (SW) & \(\rm{n_{earth}}\) & \(\mathcal{U}(0,30)\)
& \(A_{\rm{SW}}\) & \(\rm{log}_{10} \;\mathcal{U}(10^{-10}, 10^{1})\)
& \(\gamma_{\rm{SW}}\) & \(\mathcal{U}(-2,-1)\)

column1-Z = c, cell11 = r=2c, cell12 = r=2c, cell13 = c=2c, cell15 = c=2c, vline1,2,3,5,7 = 1-14, vline4,6 = 2-14, hline1,3,15 = -, hline2 = 3-5, hline2 = 5-7, rowsep = 1pt, Pulsar & ELAT (deg.) & SW modeling \(-\) DR2low & & SW modeling \(-\) DR2new+ &
& & Timing - NE_SW & Noise & Timing - NE_SW & Noise
J1022+1001 & \(-0.06\) & \(9.51 \pm 0.03\) (fitted) & n_earth + SWv & \(7.2 \pm 0.5\) (fitted) & marg.
J1730\(-\)​2304 & \(0.19\) & \(9.9 \pm 0.4\) (fitted) & n_earth & 7.9 (fixed) & marg.
J0030+0451 & \(1.45\) & \(8.20 \pm 0.03\) (fitted) & n_earth + SWv & 7.9 (fixed) & marg.
J0751+1807 & \(-2.81\) & 7.9 (fixed) & n_earth + SWv & 7.9 (fixed) & marg.
J1744\(-\)​1134 & \(11.81\) & \(3.62 \pm 0.08\) (fitted) & n_earth + SWv & 7.9 (fixed) & marg.
J1918\(-\)​0642 & \(15.35\) & 4 (fixed) & marg. & 7.9 (fixed) & marg.
J1024\(-\)​0719 & \(-16.04\) & 4 (fixed) & marg. & 7.9 (fixed) & marg.
J1738+0333 & \(26.88\) & 4 (fixed) & marg. & 7.9 (fixed) & marg.
J1713+0747 & \(30.70\) & 4 (fixed) & marg. & 7.9 (fixed) & marg.
J1857+0943 & \(32.32\) & 4 (fixed) & marg. & 7.9 (fixed) & marg.
J1012+5307 & \(38.76\) & 4 (fixed) & marg. & 7.9 (fixed) & marg.
J1640+2224 & \(44.06\) & 4 (fixed) & marg. & 7.9 (fixed) & marg.

10 Estimation of the ECORR parameters↩︎

11 Solar wind parameters↩︎

12 Timing residuals and observing radio frequencies of DR2low↩︎

References↩︎