© The Authors 2023

1 Introduction

Pulsar timing arrays (PTAs; Sazhin 1978; Detweiler 1979) search for gravitational waves (GWs) at nanohertz (nHz) frequencies by observing a suite of stable, rapidly rotating millisecond pulsars (MSPs; Backer 1993) over decadal timescales (Foster & Backer 1990; Desvignes et al. 2016; Verbiest et al. 2016). These GWs may be produced by inspiralling supermassive black holes binaries (SMBHBs; see e.g. Sesana et al. 2004; Sesana 2013), cosmic strings (Damour & Vilenkin 2001; Grishchuk 2005; Boyle & Buonanno 2008), phase transitions in the early Universe (Schwaller 2015), quantum fluctuations in the primordial gravitational field (Grishchuk 2005; Lasky et al. 2016), or from a primordial magnetic field (see Caprini & Figueroa 2018, and references therein).

Several PTA experiments are currently operational: the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; McLaughlin 2013), the Parkes Pulsar Timing Array (PPTA; Manchester et al. 2013), the Indian Pulsar Timing Array (InPTA; Joshi et al. 2022), the MeerKAT Pulsar Timing Array (MPTA) based in South Africa (Miles et al. 2023), and the Chinese Pulsar Timing Array (CPTA; Lee 2016). These regional PTAs are also organised under the International Pulsar Timing Array (IPTA; Verbiest et al. 2016; Perera et al. 2019). Although individual GW sources remain an exciting prospect for detection with PTAs (see e.g. Babak et al. 2016; Falxa et al. 2023, and references therein) it is likely that the first signal detected by PTAs will be a stochastic GW background (GWB) arising from multiple overlapping sources (Rosado et al. 2015).

The recent detection of a common red noise process by EPTA, NANOGrav, PPTA and IPTA (Chen et al. 2021; Alam et al. 2021b; Goncharov et al. 2022; Antoniadis et al. 2022, respectively) suggests that the characteristic spatial quadrupole correlation of the GWB (Hellings & Downs 1983) can become detectable with a modest expansion of PTA datasets. In this article, we present the latest EPTA data release, henceforth referred to as EPTA DR2. The release contains high-precision time-of-arrival (TOA) data for 25 MSPs, as well as the corresponding timing models. We also make available the full suite of software libraries required to reproduce these timing models. This EPTA data release paper is accompanied by two papers (EPTA Collaboration & InPTA Collaboration 2023a,b), reporting the modelling of stochastic noise processes present in individual pulsar datasets and the search for GW signals, respectively.

The article is organised as follows; in Sect. 2 we introduce the MSP selection process and provide details on the EPTA observing systems and observations. In Sect. 3, we describe the data curation and combination process, and in Sect. 4 we present the timing solutions for the 25 pulsars in our dataset. In Sect. 5, we discuss the implications of these results. We conclude with a brief summary in Sect. 6.

2 Observations and data processing

The EPTA uses data from six European radio telescopes: the Effelsberg 100 m radio telescope (EFF) in Germany, the 76 m Lovell Telescope at Jodrell Bank Observatory (JBO) in the United Kingdom, the large radio telescope operated by the Nançay Radio Observatory (NRT) in France, the 64 m Sardinia Radio Telescope (SRT) operated by the Italian National Institute for Astrophysics (INAF) through the Astronomical Observatory of Cagliari (OAC), and the Westerbork Synthesis Radio Telescope (WSRT) operated by ASTRON, the Netherlands Institute for Radio Astronomy. In addition, these telescopes regularly operate as the Large European Array for Pulsars (LEAP), which offers an equivalent diameter of up to 194 m (Bassa et al. 2016b). In this data release, we also incorporated, for the first time, data from the Mark II telescope at JBO.

The first EPTA data release (henceforth DR1; Desvignes et al. 2016) only included data from legacy data recording systems. Most of these made use of the incoherent dedispersion scheme implemented on custom-built hardware, placing a fundamental limit on achievable timing precision. In EPTA DR2, we added data from next-generation coherent dedispersion recording systems that offer a significant increase in bandwidth and sensitivity at each telescope. These new backends use hardware based on field-programmable gate arrays (FPGAs) to carry out the conversion of the electrical signal into a digital data stream and to apply polyphase filterbanks and, in some cases, some pre-filtering of the recorded band to reject known radio frequency interference (RFI). For most of the next-generation EPTA recording backends, this data processing step is implemented on the second generation of the Reconfigurable Open Architecture Computing Hardware (ROACH) platform, developed by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER; Hickish et al. 2016). In the following sub-sections, we describe the pulsar recording systems and data processing schemes for each telescope in more detail (see Table 1 for an overview).

2.1 Effelsberg Radio Telescope

The Effelsberg 100 m single-dish radio telescope located near Bad Munstereifel is maintained and operated by the Max Planck Institut für Radioastronomie. The Effelsberg Berkeley Pulsar Processor (EBPP; Backer et al. 1997) was used to record observations up to 2009. A ROACH-based backend became the main EFF pulsar backend in 2011 (Lazarus et al. 2016). This system is capable of coherently dedispersing and folding dual-polarisation data streams over a total bandwidth of up to ~500 MHz.

The observations in the Effelsberg Radio Telescope were conducted in L, S, and C bands¹. Depending on the central frequency, different receivers were used, as summarised in Lazarus et al. (2016). L-band observations were performed with the P200mm and P217mm receivers of the telescope. These offer bandwidths of 140 MHz and 240 MHz, respectively.

S -band observations around 2600 MHz were made with the S110mm receiver. The receiver frequency coverage increased from 80 MHz to 300 MHz after the second half of 2018. C-band observations were made with the S60mm receiver, which offers 500 MHz of bandwidth.

The typical integration length per source was 30 min, which included a 2 min scan of the noise diode for polarisation calibration. Data reduction of the folded coherently dedispersed pulsar archives, including flux and polarisation calibrations, was carried out via the COASTGUARD pipeline, which is part of the TOASTER software library (see Lazarus et al. 2016, for a detailed description). The removal of RFI was carried out with the pazi command in PSRCHIVE. For each pulsar, we produced a profile template at each observing band by adding the three observations with the highest signal-to-noise ratio (S/N). Data taken with the P200mm and P217mm receivers were analysed with the same template. Each observation was integrated over frequency, time, and polarisation to produce integrated pulse profiles. For L-band data, observations were divided into two subbands, except for PSRs J1738+0333, J1843−1113 and J2322+2057, which were typically weak and required integration over the entire band. TOAs were then produced by cross-correlating the pulse profiles with the respective pulse-profile templates (see Lazarus et al. 2016, for details).

2.2 Lovell Telescope

The Lovell Telescope is located at the JBO, in Cheshire, UK. It is a 76 m (250 ft) parabolic dish on an altitude-azimuth mount. The telescope is operated by the Jodrell Bank Centre for Astrophysics at the University of Manchester. Since January 2009, pulsar data have been processed using a clone of the Digital Filter Bank (DFB) developed by the Australian National Telescope Facility. From April 2011 onwards, data were simultaneously processed using the DFB and a ROACH board.

Although the DFB initially operated over a bandwidth of 128 MHz centred on 1400 MHz, from September 2009 onwards, the observing frequency coverage increased to 512 MHz centred at 1520MHz, of which the central 384MHz were typically used. ROACH observations covered a 512 MHz band centred on 1532 MHz. The edges of the band were masked, leaving a total of 400 MHz of usable bandwidth.

The typical integration time per source was varied depending on the pulsar and epoch, with the median observation times per pulsar ranging between ~ 10 min to 55 min. These observations were time-stamped using an on-site hydrogen maser clock and then corrected to Coordinated Universal Time (UTC) using recorded offsets between local time kept by the hydrogen maser and GPS time. At JBO, the observations were typically not flux or polarisation calibrated.

For the ROACH backend, data streams from two orthogonal polarisations were sampled at the Nyquist rate and digitised as 8-bit numbers. The 512 MHz band was split into 16 MHz sub-bands by a 32-channel polyphase filter. The signal in each 16 MHz sub-band was then dedispersed and folded into 1024 bins at the pulse period in real-time. Each of the 16 MHz sub-bands was further divided into 64 channels, each 0.25 MHz wide and recorded on disk in 10 s long sub-integrations. The lowest and highest 56 MHz of the band were discarded due to persistent and known contamination by RFI sources at these frequencies. Therefore, the total bandwidth recorded on disk was 400 MHz. For the DFB, the two orthogonal polarisation data streams, covering 512 MHz, were incoherently dedispersed and folded into 1024 bins at the pulse period using 0.5 MHz wide channels and 10 s sub-integrations. In this case, the lowest and highest 64 MHz of the total band were discarded due to contamination by RFI.

The mitigation of RFI in both DFB and ROACH data was carried out using a median filtering algorithm, which was followed by manual inspection. Since November 2011, ROACH data have also been cleaned in real-time using a spectral kurtosis method (Morello et al. 2022). For each pulsar observation, the pipeline produced archive files with various frequency and time resolutions. These were fully frequency-averaged archives with full time resolution (10 s long sub-integrations); fully time-averaged profiles with full frequency resolution (0.25 MHz wide and 0.50 MHz wide channels for the ROACH and DFB, respectively); archives that were partially averaged in time and frequency (1 min long sub-integrations and 8 MHz-wide and 12 MHz-wide channels for the ROACH and DFB, respectively) and finally; profiles that were fully averaged in time and frequency.

TOAs were formed for each observation using the pat command in PSRCHIVE (van Straten et al. 2012), which employs the FDM algorithm (Verbiest et al. 2016) using the template profile for the appropriate frequency. In the work presented here, for each epoch, we divided the partially (time-)averaged files to produce two archives, one for each subband centred on 1420 MHz and 1620 MHz, respectively. These files were then averaged in frequency and time individually, to finally produce two TOAs spanning the full observation duration.

2.3 Mark II Telescope

The Mark II Telescope is also located at the JBO, in Cheshire, UK. It is an altitude-azimuth mount instrument, with an elliptical dish with a major axis of 38.1 m (125 ft) and a minor axis of 25.4m (83.3 ft). From August 2017 onwards data have also been recorded for PSR J1713+0747 using this telescope with the ROACH-board-based backend described above. The data were processed in the same way as for the Lovell Telescope.

2.4 Nançay Radio Telescope

Regular EPTA timing observations were conducted with the NRT from late 2004. These observations were made using the L-and S -band receivers of the telescope, with a frequency coverage of 1.1 GHz to 1.8 GHz, and 1.7 GHz to 3.5 GHz, respectively. From late 2004 until early 2014, the legacy Berkeley-Orléans-Nançay (BON) backend (Cognard & Theureau 2006) was used to record the pulsar timing data included in the EPTA DR1 dataset (for a detailed description, see Desvignes et al. 2016). Starting in August 2011, the Nançay Ultimate Pulsar Processing Instrument (NUPPI) became the primary pulsar timing backend (Cognard et al. 2013; Liu et al. 2020). The EPTA DR2 dataset includes data collected with the NUPPI backend between August 2011 and early 2021. NUPPI observations have durations ranging from 20 min to 80 min, and cover a frequency bandwidth of 512 MHz channelised into 128 channels that are coherently dedispersed in real-time. Observations with the L-band receiver were made at a central frequency of 1484 MHz. Those with the S -band receiver were generally centred on 2539 MHz and occasionally on 1854 and 2154 MHz. Data collected before MJD 57924 were time-stamped using a local rubidium clock, and later corrected to GPS standard time stamped using recorded offsets between the clock and the Paris Observatory Universal Time. Data collected after this epoch were directly stamped with the GPS time standard, as the backend was locked to the GPS signal. Data were calibrated for polarisation using a short scan on a reference noise diode, conducted prior to each observation, with the SingleAxis² method of PSRCHIVE, to correct for the differential phase and amplitude between the two polarisations. Since late 2019, this simple calibration scheme was further complemented by regular observations of bright polarised pulsars in a mode where the horn rotated by ~180º across the 1 h observation, enabling a better determination of the polarimetric response of the NRT at the epoch of the observations. These observations and the procedure followed for analysing them are described in Guillemot et al. (2023). Automatic RFI mitigation was performed on polarisation-averaged archives with the full frequency and time resolution available (typically, 4 MHz and 5 s to 60 s time resolution), using a Python script based on the surgical RFI-cleaning algorithm of the COASTGUARD software package (Lazarus et al. 2016). Observations corrupted by strong RFI, calibration issues, incidental backend faults, or those that contained no visible pulsar signal, were discarded. We formed template profiles with the four Stokes parameters and with four frequency sub-bands by integrating the eight highest S/N observations, smoothing the average using a wavelet smoothing method. Finally, we used these polarimetric profiles and the matrix template matching (MTM) technique implemented in pat (van Straten 2006) to extract TOAs from the NUPPI observations. By modelling the transformation between polarised light curves, the MTM method corrects potential polarisation artefacts (caused by, e.g., calibration issues), yielding more accurate TOAs. The procedure developed at Nançay to improve the polarisation calibration of the NUPPI data will be presented in detail in Guillemot et al. (2023).

2.5 Sardinia Radio Telescope

The Sardinia Radio Telescope is a 64 m parabolic dish with an altitude-azimuth mount. The receiver that is used primarily for EPTA observations is the dual-band L/P receiver (original P-band = 305 MHz to 410 MHz and L-band = 1.3–1.8 GHz). Data during the EPTA observations were mainly acquired with a DFB as well as a ROACH backend, the latter of which is a copy of those installed at JBO and Effelsberg. Both backends are capable of performing real-time folding of the incoming data, each with a bandwidth of 500 MHz. Additionally, an 8-node CPU cluster installed with ROACH allows the baseband recording of the full LEAP bandwidth (128 MHz). The data were time-stamped using a local hydrogen maser. EPTA observations at the SRT started in March 2014 and went on until July 2016, and were resumed from May 2018 onwards. The intermediate gap was due to repairs of the L/P receiver, refurbishment of the active surface of the dish, and relocation of the control room along with the digital instruments to a permanent structure. Most of these observations were conducted as part of the LEAP observing programme. For the work in this paper, these data were directly integrated into the LEAP data products, rather than being included as independent pulsar timing data (see Sect. 2.7).

2.6 Westerbork Synthesis Radio Telescope

The WSRT was an equatorial mount linear array consisting of 14 dishes of 25 m diameter which were coherently added to form a 93 m equivalent dish. Pulsar observations formed a key science operation and were carried out using the Pulsar Machine (Pumal; Strom 2002) backend from 1999 until mid-2010. From 2007 to mid-2015, pulsar observations were performed using the second generation Pulsar Machine (PuMaII; Karuppusamy et al. 2008) which was a coherent dedispersion backend capable of simultaneously recording data across a total bandwidth of 64 MHz to 128 MHz. For both PuMa backends, the WSRT was equipped with the turret-mounted Multi-Frequency FrontEnds (MFFEs; Tan 1991), which allowed for rapid changes in observing frequency. From 2013 onwards, the effective size of WSRT for pulsar observations gradually decreased due to dishes being removed from the array as the infrastructure was converted to adapt the telescopes for the new phased array feeds; the APERture Tile In Focus (APERTIF; van Cappellen et al. 2022). To account for the decrease in sensitivity, the observing time per pulsar was steadily increased during this transition period, which ended on 15 June 2015 when the WSRT officially ended its multiband operations.

The PumaI backend was used in conjunction with the MFFEs to observe pulsars at 328 MHz, 382 MHz, 800 MHz, 1420 MHz and 2200 MHz, with the data recorded in a custom format. TOAs were created from the PuMaI backend by first converting the observations to ASCII total intensity profiles and then cross-correlating them with templates as described in Desvignes et al. (2016). The PuMaII backend operated across a total bandwidth of 80 MHz at 350 MHz, and 160 MHz at 1380 MHz and 2273 MHz. The total bandwidth was subdivided into an overlapped polyphase filter bank scheme with individual channels of 10 MHz and 20 MHz each with an overlap of 1 .25 MHz and 0 MHz at 350 MHz and 1380 MHz, respectively. Besides multi-band operation, PuMaII was also operated in LEAP mode (see Sect. 2.7) centred at 1398 MHz, 16.25 MHz of each 20 MHz used to match other telescopes; with a total bandwidth for timing observations set at 130 MHz.

The duration of pulsar observations at the WSRT varied with pulsar and epoch and lasted between 15 min and 45 min. However, at later epochs, the duration per pulsar was extended to account for the loss of sensitivity because fewer dishes were available. These observations were time-stamped to UTC at the observatory using a GPS-referenced local hydrogen maser. Daily averages of the local maser clock were also compared with other observatory masers as part of regular, very large baseline interferometry (VLBI) recordings. The local clock at WSRT was affected by a discrete offset of 10 s between ~56 808 MJD to 56 903 MJD due to a failure of the local Network Time Protocol (Mills 1997) server.

To achieve phase coherence between the dishes, standard interferometric coherent phasing of the array was performed at the start of each observing run. By observing an astronomical calibrator, the signal amplitudes for each polarisation were equalised for all individual antennas. This maximised the vector sum of each polarised data stream, leading to very stable polarisation behaviour, especially at 1380 MHz and 2273 MHz which was stable for at least a couple of days. No further calibration (flux or polarisation) was performed, although the polarisation response was cross-validated against calibrated data from the other telescopes when the WSRT was operated in LEAP mode.

The PuMAII backend recorded raw data in psrdada³ format, which were dedispersed using DM values stored in TEMPO-style ephemeris files. PuMaII produced PSRCHIVE archives for each band with a final sub-integration length of 10 s and 64 channels of 10 MHz and 20 MHz, depending on the observing band. The data were Nyquist sampled for each channel at a resolution of 25 ns leading to a variable number of bins for the pulsars; thus, the fastest pulsars have 256 bins, while the slowest have 8192.

The excision of the RFI in WSRT was performed using a custom tool based on PSRCHIVE, followed by manual inspection. To compensate for decreased amplitude response at the band edges and obtain continuous frequency coverage, files per frequency band were split and re-added, providing effective bandwidths of 70 MHz at 350 MHz, 160 MHz at 1380 MHz, 130 MHz at 1396 MHz (LEAP mode) and 150 MHz at 2200 MHz.

For each of the band-averaged datasets, individual observations were further reprocessed for RFI removal using ITERATIVE_CLEANER⁴, inspired by algorithms from the COASTGUARD (Lazarus et al. 2016) Pulsar Processing Suite. Each observation was then averaged over frequency and time to produce the pulse profiles. Templates for each pulsar were produced by averaging the smallest number of profiles that contribute ~90% of the S/N of the entire dataset, and TOAs were generated following the recommendations in Appendix A of (Verbiest et al. 2016), by using pat with the FDM algorithm.

2.7 The Large European Array for Pulsars (LEAP)

LEAP performs simultaneous monthly observations of more than twenty MSPs at 1.4 GHz with the five EPTA telescopes mentioned above (Bassa et al. 2016b). During each observation, the telescopes switch between the pulsar and a nearby phase calibrator. The typical exposure time per pair of pulsar and phase calibrator is 45–60 and 2–3 min, respectively. All phase calibrators were selected from the VLBA calibrator catalogue⁵. During the observations, baseband data were recorded at the Nyquist rate with 8-bit sampling and 8 × 16 MHz subbands. The data were later assembled at JBO where they were correlated, calibrated for polarisation, and coherently combined with a dedicated software correlator (Smits et al. 2017), using Effelsberg as a geometrical and time reference. Narrowband RFI from a particular telescope was also removed before the data were combined. The combined baseband data of each pulsar were then coherently dedispersed and folded to form 10 s sub-integrations with a 1 MHz frequency resolution, using the DSPSR software package. Data were visually inspected to remove any remaining impulse RFI. Next, for each pulsar, the data from each individual observation were averaged in time and frequency and used to calculate the TOA of the integrated profile with the MCMC implementation of the canonical template matching scheme (Taylor 1992; Verbiest et al. 2016). The template was created using the observation with the highest S/N (which was eventually not included in the timing dataset) using a wavelet smoothing scheme to remove the radiometer noise on top of the pulsar signal. Most of the processing of the dedispersed, folded data was carried out with the PSRCHIVE software package.

3 Data preparation, combination, and timing analysis

The EPTA DR2 includes TEMPO2-compatible pulsar ephemerides and TOA files, with the latter produced as described in the previous section and references therein. Following the customised processing steps at individual telescopes, described in Sect. 2, the final TOA sets were transferred to a central repository, for which a continuous integration scheme was developed, providing quick-look plots of the timing residuals for inspection. The data from the central repository were then combined in parallel, using standard manual steps as presented in Verbiest et al. (2016); Desvignes et al. (2016) as well as a semi-automated combination scheme. In both cases, the results were manually inspected and cross-verified for consistency, with the final data release containing data from the manual process. In the following paragraphs, we describe source selection, combination steps, and the resulting combined dataset in more detail.

3.1 Source selection

Due to the large number of telescope and back-end combinations across EPTA observatories, the heterogeneous recording and processing schemes, and the complex RFI environment in most telescopes, curation and vetting of EPTA data require long lead times. Furthermore, the noise modelling for individual pulsars is computationally expensive, requiring a detailed and iterative analysis of the possible noise models that may be applicable for each pulsar dataset. For these reasons, we adopted a source selection scheme which maximises the detectability of a stochastic GWB through the $\text{S}/{\text{NA}}_2^{}$ statistic of Rosado et al. (2015); Speri et al. (2023),

$$\text{S}/{\text{NA}}_2^{}=2{\displaystyle \sum _{a>b}\frac{{\text{Γ}}_{ab}^2{S}^2\left(f\right)T_{ab}}{P_a\left(f\right)P_b\left(f\right)}\text{d}f,}$$(1)

taking into account the fact that each pulsar contributes differently to the PTA response, due to its inherently distinct noise properties, Pj(ƒ). Here Γ_ab represents the overlap reduction function that translates the mean spectral density of an isotropic stochastic red-noise process, S (f), to cross-correlation power between pulsars a and b that have been observed for a common duration of T_ab. The simple ranking produced through this scheme was then improved by applying the coupling matrix formalism introduced in Roebber (2019) as adapted by Speri et al. (2023) to prioritise pulsars that maximise the response to HD-like correlations while maintaining the ability to distinguish between competing dipolar and monopolar signals. Using this methodology, we found that a subset of 25 pulsars out of the 42 included in the DR1, were sufficient to recover at least 90 % to 98 % of the full array sensitivity to a simulated stochastic GW background with an amplitude of 3 × 10⁻¹⁵ at a frequency of 1 yr⁻¹, and a spectral index of γ = 13/3. The same subset of pulsars would also recover at least 95% of the total sensitivity to possible individually resolvable, monochromatic gravitational wave sources across all frequencies. These 25 pulsars comprise the EPTA DR2. Their distribution on the sky can be seen in Fig. 1.

Fig. 1 Sky projection of the 25 pulsars included in the EPTA DR2 dataset.

3.2 Combination of the dataset

We followed the timing and combination steps described in Verbiest et al. (2016) and Desvignes et al. (2016) to combine the data across telescopes. For each pulsar, data from different telescopes were combined using TEMPO2 to form the joint dataset, starting with parameters from Desvignes et al. (2016) and using a summary TOA file, following Perera et al. (2019), Verbiest et al. (2016) and Desvignes et al. (2016). To align the data from different observing systems, we fitted for an arbitrary phase offset (commonly referred to as jump) for each sub-band and backend combination, using the NUPPI sub-band centred at 1420 MHz as our reference dataset. For a small number of individual backend datasets, discrete time offsets were detected and estimated using multi-pulsar analysis. These were also removed using the TEMPO2 TIME keyword. During certain observing runs, data were collected using both legacy and new backends, or in both single-telescope and LEAP modes. As these observations represented the same signal and noise, we eliminated the older backend and non-LEAP data. However, we kept the data outside the LEAP bands to better constrain dispersion-measure (DM) variations. The overlapping data were removed after the jump values were determined. With this final set of curated TOAs, we then produced initial timing solutions for each pulsar using TEMPO2. For these solutions, the timing parameters were fitted for iteratively using TEMPO2, until the linearised timing solution converged. For each of the pulsars, we then investigated the likelihood of introducing new timing parameters that were not fitted for in DR1, using a 5σ detection threshold, as well as a number of F-statistic and information criteria based tests.

All initial timing models include the spin frequency and its derivative, DM and its first and second derivative, the astrometric parameters (position, proper motions, and in several cases the annual parallax). For binary pulsars, we included fits for five Keplerian parameters and a selection of post-Keplerian (pK) parameters, depending on the pulsar. The full set of parameters included in each timing model is listed in Tables B.1-B.7. We note that we used equatorial coordinates to fit for the positions of most pulsars, except for PSRs J0030+0451, J1022+1001 and J1730−2304, for which we used ecliptic coordinates, as their ecliptic latitude is less than 1 deg. We used the DE440 version of the JPL solar system ephemeris (Park et al. 2021) and TT(BIPM2021) (Petit 2009)⁶ as our reference clock standard. We also applied the default spherical of the Solar Wind electron density model implemented in TEMPO2 to correct for solar wind-related DM variations, fixing the average density in the ecliptic plane at 1 au to 7.9 cm⁻³ (Madison et al. 2019) following Tiburzi et al. (2021), except for PSR J1022+1001, as described in Sect. 3.4.

The combination scheme described above produced the full EPTA DR2 dataset, an overview of which can be found in Fig. 2. This dataset is used in the pulsar timing analysis presented below, as well as in associated work, namely the single-pulsar noise modelling in EPTA Collaboration & InPTA Collaboration (2023a), and the search for GWs in EPTA Collaboration & InPTA Collaboration (2023b). Based on the full DR2 dataset, we also produced additional dataset versions for GW searches. Details for these versions can be found in Appendix A, as well as in EPTA Collaboration & InPTA Collaboration (2023a,b).

3.3 Outlier analysis

The EPTA DR2 dataset was checked for outliers using the following procedures. The first step to eliminate outliers was performed when compiling single telescope data, either by custom automated data flagging or manual inspection. After initial combination, outstanding outliers, such as TOAs with residuals offset by more than 10 times the root mean square (rms) of the timing solutions, were flagged and the observation archives reinspected. We found that these were typically associated with low S/N and were therefore removed. Additionally, we removed observations with known calibration issues as well as those displaying corrupted polarisation profiles or systematic trends in single-epoch timing residuals. A similar analysis was carried out using the semi-automated analysis, including tests from expected correlations such as excess contribution from the Solar Wind or epoch-wise offsets. Finally, we inspected the whitened residuals using the results of the noise modelling analysis (see the next subsection), which revealed that all remaining TOAs were within 5σ of the whitened timing residuals, indicating that no additional outlier removal was needed.

3.4 Bayesian timing analysis

To obtain the final timing solution for each pulsar, we performed a Bayesian timing analysis on the combined dataset using the TEMPONEST toolkit (Lentati et al. 2014). For this step, the timing measurements obtained with TEMPO2 were provided as initial guesses. TEMPONEST is based on the software packages TEMPO2 (Hobbs et al. 2006) and MULTINEST (Feroz et al. 2009). It explores the parameter space of a non-linear pulsar timing models using nested sampling (Skilling 2004) based on Bayesian inference (Lentati et al. 2014), to provide robust estimates for the timing parameters.

For the analysis of each pulsar, in addition to the deterministic pulsar timing parameters, we included a set of stochastic parameters to characterise the noise components in the dataset. Within our adopted framework, the white noise is described by a multiplicative factor E_f (EFAC) that accounts for the possible underestimation of TOA uncertainties, and a factor E_q that is added in quadrature (EQUAD) to model any other possible source of noise, such as pulse phase jitter and systematics (e.g. Liu et al. 2011, 2012). These are related to the uncertainty of the TOA measurement, 0%, as follows:

$$\sigma =\sqrt{E_{\text{q}}^2+E_{\text{f}}^2{\sigma }_{\text{r}}^2.}$$(2)

Such pairs of noise parameters were introduced for each back-end/frontent combination. These combinations correspond to different -group flags in the TOA files (see Table 1). The long-term red-noise processes in the data were described by two types of models that account for both the chromatic and achromatic noise components, respectively. In both cases, the noise components were modelled as a stationary, stochastic process with a power-law spectrum given by:

$$S\left(f\right)=\frac{A^2}{C_0}{\left(\frac{f}{f_r}\right)}^{-\gamma },$$(3)

where S (f), A, γ and f_r correspond to the spectral density of the power as a function of frequency f, the spectral amplitude, the spectral index, and the reference frequency (set at 1 yr⁻¹), respectively. The spectrum in both cases has a low-frequency cutoff, determined by the inverse of the total time span of the dataset. The constant C₀ was set to 1 for the chromatic noise component and to 12π² for the achromatic term. The spectral amplitude of the chromatic noise component is proportional to the observing radio frequency as A ∝ v^−α, where the chromatic index α is either 2 for modelling the features induced by the DM variations, or 4 for capturing scattering variation in the data. For each pulsar, we used a customised noise model obtained by following the procedures outlined in Chalumeau et al. (2022), which employs a Bayesian model-selection framework to determine the components in the noise model and the number of frequency bins used to sample the spectrum of each component. For 21 pulsars, the noise model includes a chromatic component to model DM variation, while for 10 pulsars it includes an achromatic component. Only one pulsar (PSR J1600−3053) favours the inclusion of a chromatic component for modelling scattering variation. J2322+2057 is the only pulsar in our sample that does not show any evidence for the presence of red noise. The noise components for each pulsar are summarised in the tables of Appendix B. More details of our customised noise models can be found in the accompanying paper on noise-modelling in this series (EPTA Collaboration & InPTA Collaboration 2023a).

During the analysis with TEMPONEST, both timing and noise parameters were sampled simultaneously, while the time offsets were analytically marginalised. All model parameters were sampled with uniform priors, except for EQUAD and the amplitude of the red-noise processes, which used a log-uniform prior. For PSR J1022+1001, whose line of sight lies very close to the ecliptic plane, we also fitted the solar wind density at the orbit distance of the Earth, instead of fixing its value to 7.9 cm⁻³ as was done for the other pulsars. This helped to obtain an estimate of the timing parallax consistent with other pulsar timing and VLBI measurements (see Sect. 4.3 for more details).

Fig. 2 Overview of the full EPTA DR2 dataset. Empty circles denote legacy data, while filled squares show new EPTA backends. Vertical lines bound the range of the ‘DR2new’ dataset. Readers can refer to Sect. 3 and Appendix A for details.

4 Results

Table 2 provides a summary of the EPTA DR2 dataset, while the full timing solutions are listed in Tables B.1-B.7. Our timeseries have time spans of 13.6 to 24.5 yr and contain ~800 to ~6000 TOAs. With the exception of PSR J1909−3744, all timing solutions use data from at least two telescopes. In the following, we describe some new timing measurements in DR2 and compare them with the results by Desvignes et al. (2016); Perera et al. (2019); Alam et al. (2021a); Reardon et al. (2021); henceforth referred to as EPTA DR1, IPTA DR2, NG12 and PPTA DR2, respectively. Unless otherwise noted, quoted values and uncertainties represent the median and 68.3% confidence intervals (C.I.) of the one-dimensional marginalised posterior distribution function.

For 12 of the 25 pulsars in DR2, thanks to the extended data and increased sensitivity of the next-generation receivers, we measure several new timing parameters in comparison to DR1. Some of these measurements are also reported here for the first time. These results are discussed in detail in the following subsections. In addition, DR2 includes improved timing solutions for the binary systems J0751 + 1807, J0900−3144, J1024−0719, J1455−3330, J1713+0747, J1857+0943, J1909−3744 and J1910+1256, as well as for the isolated pulsars J0030+0451, J1744−1134, J1843−1113, J2124−3358, and J2322+2057, even though no additional timing parameters were measured compared to EPTA DR1. For PSR J1024–0719, where we only fit the second spin frequency derivative as in previous work (Bassa et al. 2016c), we performed additional investigations to explore the measurability of higher order spin frequency derivatives (see Sect. 5.5). For PSR J1857+0943, we omitted the fit for secular variation of the orbital projected semi-major axis ($\dot{x}$) which was poorly constrained from the timing analysis. Meanwhile, a rough estimate of this parameter may still be inferred by searching for the annual orbital parallax (Sect. 5.4).

4.1 PSR J0613–0200

PSR J0613–0200 is a binary pulsar with a white dwarf (WD) companion (Bassa et al. 2016a). The new timing solution includes three pK parameters, namely the secular change of the orbital period, Ṗ_b, the third harmonic of the Shapiro delay (SD), h₃, and the ratio between the third and fourth SD harmonics, ζ (Freire &Wex 2010). The uncertainty of the measurement on Ṗ_b has improved by a factor of three compared to EPTA DR1, while estimates for the SD parameters with EPTA data are reported here for the first time. The astrometric solution includes a significant measurement of the timing parallax. These measurements agree with the constraints in NG12 and PPTA DR2, but are in slight tension with previous EPTA estimates in DR1. The results are summarised in Table B.1.

4.2 PSR J1012+5307

PSR J1012+5307 is a 1.8M_⊙ pulsar in a 14.4h orbit around a low-mass WD companion (van Kerkwijk et al. 1996; Lazaridis et al. 2009; Antoniadis et al. 2016; Ding et al. 2020; Mata Sánchez et al. 2020). The DR2 dataset for this pulsar now includes more than 5000 TOAs with an average precision of ~5 µs. Previous mass estimates (Antoniadis et al. 2016) for the system indicate that the expected SD amplitude is now larger than the EPTA timing precision for the pulsar (5 µs × 5000^−1/2 ≃ 0.06 µs). This prompted us to include the orthometric SD parameters in our fit, resulting in a significant detection of h₃ and h₄ for the first time (see Table B.2). The updated timing solution also includes a highly significant measurement of Ṗ_b and an improved estimate of the parallax, which is consistent with previous timing estimates (Lazaridis et al. 2009; Desvignes et al. 2016; Perera et al. 2019), but in tension with the estimates from VLBI and Gaia (Ding et al. 2020; Antoniadis 2020, 2021). Possible reasons for this discrepancy may be related to systematics between dynamical and kinematic reference frames (Liu et al. 2023), or covariance between the parallax, DM variability, and solar-wind timing signatures. These will be investigated in a future publication.

4.3 PSR J1022+1001

In the EPTA DR1 dataset, the only measurable pK parameter for this system was the advance of periastron, $\dot{\omega }$. The updated estimate for this parameter is in agreement with the previous value, with an uncertainty improved by five times. The new solution includes three additional pK parameters, namely the two SD parameters and Ṗ_b, all constrained with a significance greater than 3σ (see Table B.2). In particular, Ṗ_b is measured for the first time in this system. Both the SD and solar wind amplitude measurements are in agreement with those reported in Reardon et al. (2021). The updated timing parallax is slightly lower than the VLBI estimate for the system (Deller et al. 2019).

4.4 PSR J1455−3330

This MSP is in a 76 d binary system with a WD companion. The measured $\dot{x}$ is −1.98(6) × 10⁻¹⁴ ls s⁻¹ (see Table B.2), which is in agreement with the NG12 result. In addition, we report a tentative measurement of Ṗ_b for the first time. Further timing analysis may lead to a clear detection of this parameter. The measured parallax agrees with IPTA DR2 and EPTA DR1.

4.5 PSR J1600−3053

PSR J1600–3053 is a binary MSP in orbit with a WD companion of an orbital period of 14.3 day. Using the EPTA DR2 data, we obtain a measurement of the pK parameter Ṗ_b (see Table B.3) for the first time. The pK parameters measured in this system, including Ṗ_b, $\dot{x}$, $\dot{\omega }$ and SD, are all consistent with those reported by NG12 and PPTA DR2.

4.6 PSR J1640+2224

PSR J1640+2224 is an MSP with a WD companion in a 175-day orbit (Deng et al. 2020). The new timing solution was derived using the DDH model in TEMPO2 and includes significant measurements of the timing parallax, SD parameters, and the first derivative of the orbital period. The SD parameters suggest that the pulsar has a low mass, although the measurement uncertainties are still too large to provide a stringent constraint. The parallax estimate is consistent with that derived with VLBI and the timing solutions by other PTAs (see Sect. 5.1). The secular orbital period change can be attributed to kinematic effects, although at face value our estimate is in slight tension with expectations, given the position and proper motion of the pulsar (see Sect. 5.3). The uncertainty for the $\dot{x}$ estimate is improved by three (five) times compared to EPTA DR1 (IPTA DR2). The timing parameters are listed in Table B.3.

4.7 PSR J1730−2304

The proper motion along the ecliptic latitude was not measured in either EPTA DR1, IPTA DR2 or PPTA DR2, due to the proximity of this isolated pulsar to the ecliptic plane. Here, we report the first tentative measurement of this parameter, −4.4(1.8) mas yr⁻¹. The revised timing solution is explained in Table B.3.

4.8 PSR J1738+0333

A pulsar timing analysis of PSR J1738+0333 based on data collected with the Arecibo and EPTA telescopes resulted in a significant detection of the timing parallax and first orbital period derivative (Freire et al. 2012). The latter is consistent with the general relativity prediction given the masses of the pulsar and its companion, which have been measured using optical spectroscopy (Antoniadis et al. 2012). Here, we report for the first time a 4σ measurement of Ṗ_b, = −3.0(7) × 10⁻¹⁴, using EPTA-only data. This estimate is approximately 2σ higher than that of Freire et al. (2012) (−1.7(3) × 10¹⁴), which could be due to systematics introduced by the long integration times typically used in EPTA observations (~0.5–1 h). The latter correspond to a significant fraction of the 7.5 h binary orbit and may result in smearing of the Ṗ_b, timing signature.

4.9 PSR J1751−2857

This 110-day period binary pulsar (Stairs et al. 2005) is only monitored by the EPTA. Compared to EPTA DR1, the new solution (Table B.4) provides significantly improved estimates for all timing parameters. We also report a marginal detection of the timing parallax; ϖ = 1.1(4) mas.

4.10 PSR J1801−1417

PSR J1801−1417 is an isolated pulsar and only monitored by the EPTA. Here, we report a marginal detection of the timing parallax signature at approximately 3σ (see Tables 3 and B.4).

4.11 PSR J1804−2717

PSR J1804−2717 is a binary MSP with a WD companion of an orbital period of 11.1 day. It is currently only monitored by the EPTA. The new timing solution includes a significant measurement of the timing parallax, as well as proper motion parameters that are seven times more precise compared to the EPTA DR1 solution (see Tables 3 and B.5).

4.12 PSRJ1911+1347

PSR J1911 + 1347 is a 4.63-ms isolated MSP. The astrometric solution for this pulsar now includes a parallax measurement (see Table 3). This is also the first parallax measurement reported in this pulsar. The new timing solution is consistent with EPTA DR1 and IPTA DR2, with a significant improvement in the measurement precision of all timing parameters.

4.13 PSR J1918−0642

We report an improved measurement of the SD signature, with h₃ = 8.3(2) × 10⁻⁷ µs, and ζ = 0.908(9), which agrees well with previous estimates. Compared to EPTA DR1, the timing solution now also includes a measured timing parallax (as seen in Table B.6) that is consistent with the value reported in NG12, as well as the VLBI parallax reported by Ding et al. (2023), as can be seen in Table 3.

5 Discussion

5.1 Parallaxes and distances

We report timing parallaxes for 20 of the systems contained in this data release. For PSRs J1640+2224, J1751−2857, J1801−1417, J1804−2717, J1911+1347, and J1918−0642 the parallax was not measured in EPTA DR1, while for PSRs J0900−3144, J1738+0333, J1843−1113, J1910+1256, J2322+2057, we still do not detect a parallax signature. In what follows we briefly compare our new results with EPTA DR1 and the estimates obtained by NANOGrav and VLBI (Vigeland et al. 2018; Deller et al. 2019; Alam et al. 2021b; Ding et al. 2023).

For sources monitored by both EPTA and NANOGrav, we find that the parallax estimates are in agreement, with the EPTA generally being more precise, due to the improved sampling resulting from the larger number of observations. Unlike NG12, we do not find a significant parallax signature for PSR J2322+2057.

To obtain distance estimates we invert the parallaxes using Lutz-Kelker-like volume and luminosity priors (Lutz & Kelker 1973; Binney & Merrifield 1998), as described in Verbiest et al. (2010, 2012); Igoshev et al. (2016)⁷. More specifically we use a volumetric prior given by,

$$P_D\left(\varpi \right)\propto \frac{R^{1.9}}{{\varpi }^4}\exp \left(-\frac{\left|\sin b\right|}{\varpi E}-5\frac{R-R_0}{R_0}\right),$$(4)

where E is the scale height, assumed to be 500 pc for MSPs, b is the Galactic latitude of the source, R₀ is the distance of the sun from the galactic centre (8.5 kpc), and R is the distance between the source and the galactic centre as,

$$R=\sqrt{R_0^2+\frac{\cos b}{\varpi }-2R_0\frac{\cos b\cos l}{\varpi },}$$(5)

where l is the Galactic longitude of the source. For the luminosity prior we adopt,

$$P_L\left(\varpi \right)\propto \frac{1}{\varpi }\exp \left(-0.5{\left[\frac{\log S_{1400}+1.1-2\log \varpi }{0.9}\right]}^2\right),$$(6)

where S₁₄₀₀ is the mean flux density of the pulsar at 1400 MHz. The use of these priors ensures that the posterior distribution function of the distance is well behaved, so that statistical moments can be defined even in the presence of large uncertainties.

A summary of the distance estimates obtained with the method described above is given in Table 3. Four MSPs in our samples, PSRs J1640+2224, J1751−2857, J1801−1417 and J1804−2717, have a poorly constrained parallax and thus the corresponding distance estimates depend sensitively on the prior assumptions and are highly asymmetric. Transverse velocities, calculated using the bias-corrected distances, and the measured proper motions, are also listed in Table 3. The proper motions derived here are in most cases consistent with VLBI estimates, but have approximately three times smaller uncertainties (Deller et al. 2019). For PSR J1022+1001, where numerous authors have reported inconsistent parallax determinations from pulsar timing and VLBI campaigns, we find a strong correlation between the parallax and the adopted model for the solar wind electron density. This will be explored in detail by a future work (Liu et al., in prep.).

5.2 Pulsar mass measurements

We detect the SD signature with a significance greater than 3σ in nine systems. Figure 3 shows the marginalised posterior probability functions for the masses of the corresponding pulsars. These estimates are based solely on SD posteriors and do not take into account all available information, such as additional pK parameters or independent measurements from optical spectroscopy. A more robust analysis of the astrophysical parameters for EPTA DR2 pulsars will be presented in a future publication.

We find that the mass estimates derived from the SD posteriors are generally consistent with the measurements obtained with previous iterations of the EPTA data, as well as those derived with independent datasets (e.g. see Desvignes et al. 2016; Verbiest et al. 2016; Perera et al. 2018, 2019; Liu et al. 2020; Reardon et al. 2021; Alam et al. 2021a). In most cases, the EPTA DR2 constraints are significantly more precise.

An SD detection for PSRs 1012+5307 and J1022+1001 is reported for the first time. For the former system, the SD parameter posteriors are not yet constrained with sufficient precision to provide informative constraints on the pulsar mass. If one considers the mass ratio of the system, q ≡ m_p/m_c = 10.44(11) (Mata Sánchez et al. 2020), and discards the SD of posterior samples that correspond to physically meaningless values for the companion mass (e.g. m_c/M_⊙ ≲ 0.14 or ≳ 0.3), the pulsar mass is constrained to be in the 1.7–2.0 M_⊙ range. As the precision of h₃ and h₄ improves with more data, the SD signature will ultimately provide an independent measurement of the component masses. Combined with the parallax, Ṗ_b, and WD spectroscopic estimates, this will significantly improve the constraints on gravitational dipole radiation, also providing a valuable test for low-mass WD cooling models.

For PSR J1022+1001, the posterior distribution of the SD parameters suggests that the pulsar mass is relatively low. For the corresponding constraints on the component masses, the expected magnitude for the relativistic periastron advance is ${\dot{\omega }}_{\text{GR}}\simeq 0.01°{\text{yr}}^{-1}$, which is ~2.5σ higher than the measured value for this parameter. This can be explained either by an extremely low pulsar mass (~1 M_⊙), or by systematics and non-relativistic contributions to the measured $\dot{\omega }$.

Fig. 3 Posterior probability distributions for mp, for the 9 systems with significant SD measurements. Solid and dashed lines represent the median and 68.3% C.I. of the marginalised posterior distribution function. In each panel, the inset provides a zoomed view around the peak of the distribution.

5.3 Secular variation of the orbital period

Using EPTA DR2, we measured binary orbital period derivatives for a total of ten pulsars. The measurements in PSRs J1022+1001, J1600−3053, J1640+2224, J1713+0747 (all with ≳3σ significance), PSR J1455–3330 (marginal), and J1738+0333 are new compared to those obtained with the EPTA DR1 dataset. Overall, the measured values reported here are consistent with EPTA DR1 (see Table 4), although with improved uncertainties.

GW emission causes a change in the orbital period in binary systems (Peters 1964). The observed Ṗ_b however, may also include extrinsic contributions such as distance- and location-dependent kinematic effects. We determined these extrinsic contributions to the observed Ṗ_b using the bias-corrected distances given in Table 3. For PSRs J1600–3053, J1909–3744 and J1738+0333, we used the reference parallax values from Reardon et al. (2021), Liu et al. (2020) and Freire et al. (2012), respectively (shown in Table 3). We consider two kinematic contributions. The first is the so-called Shklovskii effect (Shklovskii 1970), which scales with the proper motion of a binary

$$\frac{{\dot{P}}_{\text{b,Shk}}}{P_{\text{b}}}=\frac{{\mu }^{2}D}{c},$$(7)

where d is the distance, µ is the proper measured motion, and c is the speed of light, respectively. The second is caused by differences in the acceleration by the galactic gravitational potential at the solar system and binary pulsar locations (Damour & Taylor 1991). To correct for this, we used the Milky Way potential model of McMillan (2017), which considers the gas disc and the halo density profile, and is calibrated against Galactic maser sources.

We then calculated the intrinsic orbital period change (Ṗ_b,int) by subtracting the kinematic contributions from the observed value,

$${\dot{P}}_{\text{b,Int}}={\dot{P}}_{\text{b,Obs}}-{\dot{P}}_{\text{b,Shk}}-{\dot{P}}_{\text{b,Gal}}.$$(8)

From Table 4, it can be seen that for all pulsars except PSRs J0751+1807, J1012+5307, J1640+2224 and J1738+0333, the derived Ṗ_b,int is consistent with zero within 2σ.

We also calculated the change in the orbital period by GW damping (Ṗ_b,GW) using the mass estimates presented in Sect. 5.2, where available. We note that here we used the measurement of the SD parameter h₃ and assumed a 1.4 M_⊙ pulsar mass to estimate Ṗ_b,GW for PSR J0751+1807. The masses of PSR J1738+0333 were directly taken from Antoniadis et al. (2012). For PSRs J0751+1807, J1012+5307 and J1738+0333, the derived Ṗ_b,int is consistent with the predicted Ṗ_b,GW, suggesting that gravitational damping can account for the non-zero intrinsic secular variation of orbital period.

For PSR J1640+2224, the timing parallax and corresponding distance estimate agree well with those inferred using VLBI measurements (Ding et al. 2023). Therefore, the marginal ~2.7σ detection of Ṗ_b,int, which cannot be accounted for by GW damping, is unlikely to be due to an incorrect distance estimate. Although the most probable explanation for the discrepancy is the low-significance measurement of Ṗ_b,Obs (~4.9σ), another intriguing explanation could be that the system experiences an additional acceleration due to a nearby object.

For pulsars with a Ṗ_b,int consistent with either zero or anticipated contribution from GW damping, we calculated the kinematic distances Ṗ_b,obs, by assuming a net zero Ṗ_b (Bell & Bailes 1996). More specifically, we estimate the kinematic distances as

$$D_{\text{kin}}=\frac{c\left({\dot{P}}_{\text{b,Obs}}-{\dot{P}}_{\text{b,Gal}}-{\dot{P}}_{\text{b,GW}}\right)}{{\mu }^2{P}_{\text{b}}}.$$(9)

For PSR J1455–3330, the masses were not measured with our data and the kinetic distance is determined without considering the orbital period derivative GW damping (which should nonetheless be negligible considering its 76-day orbital period). The results are summarised in Table 4. Comparison between the kinematic distance and the bias-corrected parallax distance indicates that they are, in general, consistent with each other.

5.4 Searches for annual orbital parallax signatures

For binary pulsars, the proper motion of the system changes the viewing geometry of the orbit with respect to the Earth. This effect induces an apparent secular variation in the projected semi-major axis of the binary orbit (Kopeikin 1996)⁸:

$$x=x_{\mathrm{int}}\left[1+\cot i\left({\mu }_{\alpha }\cos \text{Ω}-{\mu }_{\delta }\sin \text{Ω}\right)\left(t-t_0\right)\right],$$(10)

which gives

$$\frac{\dot{x}}{x_{\mathrm{int}}}=\mu \cot i\cos \left({\theta }_{\mu }+\text{Ω}\right).$$(11)

Here $\mu \equiv \sqrt{{\mu }_{\alpha }^2+{\mu }_{\delta }^{29}}$, Ω is the longitude of the ascending node of the orbit and θ_µ is the position angle of the proper motion on the sky. It can be seen that if the proper motion and orbital inclination angle are measured, a measurement of $\dot{x}$ can then be used to determine Ω, with an ambiguity of cos(–θ_µ – Ω) = cos(θ_µ + Ω). Typically, the orbital inclination is obtained from the SD parameter sin i, with a π – i ambiguity. Therefore, the determination of Ω can have four possible pairs of solutions for the orbital inclination and ascending node (e.g. Liu et al. 2020). Meanwhile, for binary pulsars, the annual motion of the Earth around the Sun also introduces an apparently periodic variation in the viewing geometry of the orbit, an effect known as the annual orbital parallax. The variation in the projected semi-major axis as a result of this effect can be written as (Kopeikin 1995):

$$x=x_{\mathrm{int}}\left[1-\frac{\cot i}{d}\left({\text{Δ}}_{I_0}\cos \text{Ω}-{\text{Δ}}_{J_0}\sin \text{Ω}\right)\right],$$(12)

where

$${\text{Δ}}_{I_0}=-X\sin \alpha +Y\cos \alpha ,$$(13)

$${\text{Δ}}_{J_0}=-X\sin \delta \cos \alpha -Y\sin \delta \sin \alpha +Z\cos \delta .$$(14)

Here, r = {X, Y, Z} is the position vector of the Earth with respect to the barycentre of the solar system, and (α, δ) are the spherical coordinates of the barycentre of the binary system. For nearby binary pulsars in a wide orbit, the annual orbital parallax can be measurable, giving a unique pair of solutions for the inclination and ascending node (e.g. for PSR J1713+0747). We investigated sources with measured $\dot{x}$ and SD parameters to search for this annual orbital parallax signatures. These sources include PSRs J1600–3053, J1857+0943¹⁰, J1640+2224, J1022+1001, J1012+5307 (Fig. 4). We used the T2 binary model in TEMPO2, where the annual orbital parallax effect is described with the KOM and KIN parameters, corresponding to the longitude of the ascending node (Ω) and inclination angle (i), respectively. The SD parameter s = sin i is then treated as a function of KIN. We mapped the KOM-KIN space with the TEMPONEST toolkit, following a scheme similar to described in Desvignes et al. (2016). In detail, we fixed the set of white noise parameters to their maximum likelihood values and chose to analytically marginalise over the spin and astrometric parameters. We set customised uniform linear priors for KOM, KIN and M_c. For KIN and M_c the prior ranges are set to include all possible values of these parameters that are allowed from the SD measurements (see tables in Appendix B). For KOM the prior range was set from 0° to 360°. The sampling was conducted with the constant efficiency option turned off, in order to map the multimodal parameter space more effectively.

For PSR J1600–3053, four possible solutions were found in the mapping of KIN and KOM using DR1 dataset (Desvignes et al. 2016), which was consistent with a negligible annual orbital parallax parallax. With the EPTA DR2 data, we are now able to single out one solution for these two parameters: $\text{Ω\hspace{0.17em}=\hspace{0.17em}}{82}_{-8}^{+12}\mathrm{deg}$, $i={110.4}_{-0.7}^{+0.9}\mathrm{deg}$. Based on the number of samples (~337k) in the posterior distribution, the logarithmic likelihood ratio of this solution to the others is ≳5.5. This means that the signature of the annual orbital parallax is clearly detected in this binary pulsar system.

For PSR J1640+2224, two possible solutions for KOM and KIN have been found; ${\text{Ω=5}}_{-2}^{+3}\mathrm{deg}$, $i={111}_{-5}^{+4}\mathrm{deg}$, and $\text{Ω}=-{332}_{-5}^{+3}\mathrm{deg}$, $i={66}_{-3}^{+4}\mathrm{deg}$, respectively. The logarithmic likelihood ratio of these two islands is approximately 2.1, which favours the former solution. Thus, there is weak evidence for the detection of the annual orbital parallax signal.

In PSR J1012+5307, from the KOM-KIN mapping, the longitude of the ascending node and the orbital inclination are restricted to two regions. The 2σ range of Ω in these two regions, are 47–142 and 218–318 deg, respectively. Although the most likely value of Ω is around 150 deg, the preference for this area of the solution is unclear.

For PSR J1022+1001, four solutions in KOM and KOM have been found. These are 1). ${\text{Ω=17}}_{-3}^{+2}\mathrm{deg}$, $i={56}_{-4}^{+4}\mathrm{deg}$. ${\text{Ω=50}}_{-2}^{+3}\mathrm{deg}$, $i={56}_{-4}^{+4}\mathrm{deg}$, 3). $\text{Ω}={197}_{-3}^{+2}\mathrm{deg}$, $i={125}_{-4}^{+4}\mathrm{deg}$, 4). $\text{Ω}={229}_{-2}^{+3}\mathrm{deg}$, $i={125}_{-4}^{+4}\mathrm{deg}$. Although two of them are slightly preferred in probability, they cannot be distinguished with high confidence, which means that the signal of annual orbital parallax is not detected.

Desvignes et al. (2016) did not report measurable secular variation of the projected semi-major axis in J1857+0943, and with the new data set we also do not measure significant $\dot{x}$. Still, as in Desvignes et al. (2016), we nonetheless attempted to search for an annual orbital parallax to see if the constraint on the longtidue of the ascending node can be improved. In Desvignes et al. (2016), there was only a tentative constraint on the ascending node. With the new dataset, there is now a tentative preference for Ω in the range of 100 to 200 deg. For the orbital inclination, the posterior space around 88 deg is slightly preferred, compared to that peaked around 92 deg. However, there is no clear evidence for neither of the two solutions.

The analysis of PSR J1713+0747 yields Ω = 91.1 ± 0.5 deg and i = 71.3 ± 0.2 deg as discussed in Sect. 4. These are consistent with previous analyses (Zhu et al. 2019; Alam et al. 2021b; Reardon et al. 2021), as well as with the results based on EPTA DR1 (Desvignes et al. 2016).

Fig. 4 Posterior distributions of KOM and KIN from the search for annual-orbital parallax effect, for PSRs J1600–3053 (top left), J1640+2224 (top right), J1022+1001 (middle left), J1012+5307 (middle right), and J1857+0943 (bottom). The solid, dashed and dotted lines in the side panels stand for the median, 1σ confidence interval and 2σ confident interval of the distribution. The stars in the plots remark the maximum likelihood values of the two parameters each obtained from the posterior distribution in the side panel.

5.5 High order spin frequency derivatives in J1024-0719

There is a significant measurement of the second spin frequency derivative in PSR J1024–0719 (Bassa et al. 2016c; Kaplan et al. 2016). Based on these results, it is anticipated that this pulsar is in orbit with a main-sequence star named 2MASS J10243869–0719190, of an orbital period longer than 200 yr (Bassa et al. 2016c). Here with the new EPTA dataset, in addition to our timing analysis presented in Table B.2, we conducted another round of analysis including modelling of the third, fourth and fifth spin frequency derivatives. These derivatives, if measurable, can be used to derive the properties of the Keplerian orbit of the pulsar (Bassa et al. 2016c). Table 5 shows the results of this analysis, and a comparison with previous work. As can be seen, there is a tentative change in the measured second spin frequency derivative compared to previous results (even when the reference epoch of the spin frequency and the number of modelled spin frequency derivatives are both the same). This is expected given that the EPTA DR2 dataset covers a slightly longer timeline/orbit phase of the binary system.

It is also intriguing to see that with the new and more sensitive data, there is now a tentative detection of the fifth spin frequency derivative. This though needs to be confirmed by future observations.

6 Conclusions

In this paper, we presented the EPTA DR2 dataset and results from a combined frequentist and Bayesian timing analyses. This dataset contains high-precision timing data for 25 MSPs, collected with the five largest radio telescopes in Europe along with LEAP. The DR2 dataset combines data from EPTA DR1 (Desvignes et al. 2016) with those recorded by a new generation of data acquisition systems. The dataset extends to the beginning of 2021 and has baselines ranging from 14 to 25 yr.

We conducted timing analysis of the dataset based on a Bayesian framework to measure the timing parameters of the pulsars. This has yielded a collection of new measurements including annual parallaxes, secular variation of orbital period, Shapiro delay and so forth in some pulsars. We also derived a group of astrophysical parameters of these pulsars, including distances, transverse velocities, binary masses, and annual orbital parallaxes. The DR2 dataset builds the foundation for searching for GW signals. The results of this search are reported in an accompanying publication (EPTA Collaboration & InPTA Collaboration 2023b).

Acknowledgements

The European Pulsar Timing Array (EPTA) is a collaboration between European and partner institutes, namely ASTRON (NL), INAF/Osservatorio di Cagliari (IT), Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), the University of Manchester (UK), the University of Birmingham (UK), the University of East Anglia (UK), the University of Bielefeld (GER), the University of Paris (FRA), the University of Milan-Bicocca (IT), the Foundation for Research and Technology, Hellas (GR), and Peking University (CHN), with the aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves. An Advanced Grant of the European Research Council allowed to implement the Large European Array for Pulsars (LEAP) under Grant Agreement Number 227947 (PI: M. Kramer). The EPTA is part of the International Pulsar Timing Array (IPTA); we thank our IPTA colleagues for their support and help with this paper and the external Detection Committee members for their work on the Detection Checklist. Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a Consolidated Grant (ST/T000414/1) from the UK’s Science and Technology Facilities Council (STFC). I.C.N. is also supported by the STFC doctoral training grant ST/T506291/1. The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from “Programme National de Cosmologie and Galaxies” (PNCG), and “Programme National Hautes Energies” (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). The Sardinia Radio Telescope (SRT) is funded by the Department of University and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as a National Facility by the National Institute for Astrophysics (INAF). The work is supported by the National SKA programme of China (2020SKA0120100), Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. This work is also supported as part of the “LEGACY” MPG-CAS collaboration on low-frequency gravitational wave astronomy. J.A. acknowledges support from the European Commission (Grant Agreement number: 101094354). J.A. and S.Cha. were partially supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of the “Science and Society - Action Always strive for excellence – Theodoras Papazoglou” (Project Number: 01431). A.C. acknowledges support from the Paris Île-de-France Region. A.C., A.F., A.Se., A.Sa., E.B., D.I., G.M.S., M.Bo. acknowledge financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). G.D., K.Li., R.K. and M.K. acknowledge support from European Research Council (ERC) Synergy Grant “BlackHoleCam”, Grant Agreement Number 610058. This work is supported by the ERC Advanced Grant “LEAP”, Grant Agreement Number 227947 (PI M. Kramer). A.V. and P.R.B. are supported by the UK’s Science and Technology Facilities Council (STFC; grant ST/W000946/1). A.V. also acknowledges the support of the Royal Society and Wolfson Foundation. J.P.W.V. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through thew Heisenberg programme (Project No. 433075039) and by the NSF through AccelNet award #2114721. N.K.P. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer PO 2758/1–1, through the Walter-Benjamin programme. A.Sa. thanks the Alexander von Humboldt foundation in Germany for a Humboldt fellowship for postdoctoral researchers. A.Po., D.P. and M.Bu. acknowledge support from the research grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016 (Italy). R.N.C. acknowledges financial support from the Special Account for Research Funds of the Hellenic Open University (ELKE-HOU) under the research programme “GRAVPUL” (K.E.-80383/grant agreement 319/10-10-2022: PI N. A. B. Gizani). E.v.d.W., C.G.B. and G.H.J. acknowledge support from the Dutch National Science Agenda, NWA Startimpuls – 400.17.608. B.G. is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN. L.S. acknowledges the use of the HPC system Cobra at the Max Planck Computing and Data Facility. The work presented here is a culmination of many years of data analysis as well as software and instrument development. In particular, we thank Drs. N. D’Amico, P.C.C. Freire, R. van Haasteren, C. Jordan, K. Lazaridis, P. Lazarus, L. Lentati, O. Löhmer and R. Smits for their past contributions. We also thank Dr. N. Wex for supporting the calculations of the galactic acceleration as well as the related discussions. The EPTA is also grateful to staff at its observatories and telescopes who have made the continued observations possible. Author contributions: The European Pulsar Timing Array: JA, SB, ASBN, CGB, ABe, MBo, EB, PRB, MBu, RNC, AC, DJC, SCha, SChe, IC, GD, MF, RDF, AF, JRG, BG, EG, JMG, LG, YJG, HH, FI, DIV, JJan, JJaw, GHJ, AJ, RK, EFK, MJK, MK, MAK, KLa, KJL, KLi, YL, AGL, JWM, RAM, MBM, ICN, APa, BBPP, DP, APe, NKP, APo, HQL, ASa, SAS, ASe, GS, LS, RS, BWS, SCS, GT, CT, EvdW, AV, VVK, JPWV, JW, LW and ZW. The EPTA is a multi-decade effort and all authors have contributed through conceptualisation, funding acquisition, data-curation, methodology, software and hardware developments as well as (aspects of) the continued running of the observational campaigns, which includes writing and proofreading observing proposals, evaluating observations and observing systems, mentoring students, developing science cases. All authors also helped in (aspects of) verification of the data, analysis and results as well as in finalising the paper draft. Specific contributions from individual EPTA members are listed in the CRediT (https://credit.niso.org/) format below. The Indian Pulsar Timing Array: PA, SA, MB, ABa, SDa, DD, SDe, NDB, CD, YG, SH, BCJ, FK, DK, TK, NK, MAK, YM, KN, AKP, TP, PR, JS, ASr, MS, ASu, KT and PT. InPTA members contributed in uGMRT observations and data reduction to create the InPTA data set which is employed while assembling the DR2full+ and DR2new+ data sets. J.Jan., K.Li., G.M.S. equally share the correspondence of the paper. CRediT statement: Conceptualisation: APa, APo, AV, BWS, CGB, CT, GHJ, GMS, GT, IC, JA, JJan, JPWV, JW, JWM, KJL, KLi, MK. Methodology: APa, AV, DJC, GMS, IC, JA, JJan, JPWV, JWM, KJL, KLi, LG, MK. Software: AC, AJ, APa, CGB, DJC, GMS, IC, JA, JJan, JJaw, JPWV, KJL, KLi, LG, MJK, RK. Validation: AC, APa, CGB, CT, GMS, GT, IC, JA, JJan, JPWV, JWM, KLi, LG. Formal Analysis: APa, CGB, DJC, DP, EvdW, GHJ, GMS, JA, JJan, JPWV, JWM, KLi. Investigation: APa, APo, BWS, CGB, DJC, DP, GMS, GT, IC, JA, JJan, JPWV, JWM, KLi, LG, MBM, MBu, MJK, RK. Resources: APa, APe, APo, BWS, GHJ, GMS, GT, HH, IC, JA, JJan, JPWV, JWM, KJL, KLi, LG, MJK, MK, RK. Data Curation: AC, AJ, APa, BWS, CGB, DJC, DP, EG, EvdW, GHJ, GMS, GT, HH, IC, JA, JJan, JPWV, JWM, KLi, LG, MBM, MBu, MJK, MK, NKP, RK, SChe, YJG. Writing – Original Draft: APa, GMS, JA, JJan, KLi, LG. Writing – Review and Editing: AC, AF, APa, APo, DJC, EB, EFK, GHJ, GMS, GT, JA, JJan, JPWV, JWM, KLi, MK, SChe, VVK. Visualisation: APa, GMS, JA, JJan, KLi. Supervision: APo, ASe, AV, BWS, CGB, DJC, EFK, GHJ, GMS, GT, IC, JA, JPWV, KJL, KLi, LG, MJK, MK, VVK. Project Administration: APo, ASe, AV, BWS, CGB, CT, GHJ, GMS, GT, IC, JJan, JPWV, JWM, KLi, LG, MK. Funding Acquisition: APe, APo, ASe, BWS, GHJ, GT, IC, JA, JJan, LG, MJK, MK. Data Availability: Files containing the TOAs and pulsar timing models can be downloaded from https://epta.pages.in2p3.fr/epta-dr2/. The same site also provides the software environments necessary to recreate the timing analysis presented in this article.

Appendix A Dataset versions

The full EPTA DR2 dataset, hereafter referred to as “DR2full”, is the primary source of data for the pulsar timing analyses presented in this paper, as well as the accompanying articles on the single-pulsar noise analysis (EPTA Collaboration & InPTA Collaboration 2023a) and the search for correlated signals in (EPTA Collaboration & InPTA Collaboration 2023b). Apart from the DR2full, additional datasets were created, guided by iterative analyses carried out during the detailed investigation of single pulsar noise models in EPTA Collaboration & InPTA Collaboration (2023a), as well as multiple analyses carried out to search for artefacts and spurious signals in the common signal searches carried out in EPTA Collaboration & InPTA Collaboration (2023b).

The first of these, the ‘DR2full+’, consists of a combination of the DR2full dataset with the first data release of the InPTA (henceforth InPTA-DR1 Tarafdar et al. 2022) for an overlapping set of 10 pulsars, marked with an asterisk in Table 2. This dataset aims to improve the overall sensitivity of the DR2full to DM-variation linked noise-processes. The combination was performed using the ‘narrowband’ TOAs from the InPTA-DR1 for two observing bands centred around 500 MHz and 1460 MHz, respectively. Unlike the EPTA data, the InPTA-DR1 data requires optimisation based on the flux density of the pulsar being observed, as well as, observing setup specific details leading to changes in the number of sub-bands generated, ranging from a minimum of 4 to a maximum of 16. We refer interested readers to Tarafdar et al. (2022) for the full details on the InPTA-DR1 data processing. To combine the InPTA-DR1 data with the EPTA DR2 data, for each pulsar we used the ephemeris produced from the EPTA DR2 dataset and followed the same steps as described in Section 3. Specifically, we fitted for phase offset for each subband in the InPTA data individually with respect to the reference system in the EPTA data (see Section 3).

The single-pulsar noise analysis of these additional sets of data are presented in EPTA Collaboration & InPTA Collaboration (2023a). For the ‘DR2full+’ dataset, we performed single-pulsar Bayesian timing analysis following the same description as in Section 3.4, using customised noise models determined by EPTA Collaboration & InPTA Collaboration (2023a). The results of this analysis were later used for the gravitational wave searches presented in EPTA Collaboration & InPTA Collaboration (2023b).

A significant difference between the EPTA DR1 and DR2 datasets lies in the use of the technique of coherent dedispersion. This allows for a far more accurate modelling of the frequency-dependent DM delays, leading to the recovery of a sharper profile and thus, improved timing performance. To test the improvements derived from using only coherently dedispersed data, the ‘DR2new’ dataset was created. This dataset spans, approximately, the final 10.3 yr. Similar to the DR2full dataset, we also appended InPTA-DR1 data to the DR2new to produce the ‘DR2new+’ dataset.

Appendix B Pulsar Ephemerides

The measurements of the pulsar timing parameters obtained from the Bayesian timing analysis are presented in the following pages, in Tables B.1–B.7 for all 25 EPTA DR2 pulsars. The measurement values are the medians of the posterior distributions and the errors denote the 1-σ confidence intervals for the respective parameters.

Appendix C Timing residuals

In Figs. C.1–Figs. C.5, we present the timing residuals of the 25 pulsars in the EPTA DR2. Residuals both before and after the application of red-noise subtraction are shown.

Fig. C.1 Timing residuals of PSRs J0030+0030, J0613-0200, J0751 + 1807, J0900−3144, J1012+5307. For each pulsar, the residuals before and after subtraction of DM and monochromatic red noise are shown. The squares, circles and triangles represent P-band, L-band and S/C-band observations, respectively (see Table 2 for frequency coverage of each band). The blue/filled and black/unfilled symbols indicate the new backend data and those from EPTA DR1, respectively.

Fig. C.2 Timing residuals of PSRs J1022+1001, J1024−0719, J1455−3330, J1600−3053, J1640+2224. Figure style is the same as Figure C.1.

Fig. C.3 Timing residuals of PSRs J1713 +0747, J1730−2304, J1738+0333, J1744−1134, J1751−2857. Figure style is the same as Figure C.1.

Fig. C.4 Timing residuals of PSRs J1801−1417, J1804−2717, J1843−1113, J1857+0943, J1909−3744. Figure style is the same as Figure C.1.

Fig. C.5 Timing residuals of PSRs J1910+1256, J1911+1347, J1918−0642, J2124−3358, J2322+2057. Figure style is the same as Figure C.1.

References

All Tables

All Figures

	Fig. 1 Sky projection of the 25 pulsars included in the EPTA DR2 dataset.
In the text

	Fig. 2 Overview of the full EPTA DR2 dataset. Empty circles denote legacy data, while filled squares show new EPTA backends. Vertical lines bound the range of the ‘DR2new’ dataset. Readers can refer to Sect. 3 and Appendix A for details.
In the text

	Fig. 3 Posterior probability distributions for mp, for the 9 systems with significant SD measurements. Solid and dashed lines represent the median and 68.3% C.I. of the marginalised posterior distribution function. In each panel, the inset provides a zoomed view around the peak of the distribution.
In the text

Fig. 4 Posterior distributions of KOM and KIN from the search for annual-orbital parallax effect, for PSRs J1600–3053 (top left), J1640+2224 (top right), J1022+1001 (middle left), J1012+5307 (middle right), and J1857+0943 (bottom). The solid, dashed and dotted lines in the side panels stand for the median, 1σ confidence interval and 2σ confident interval of the distribution. The stars in the plots remark the maximum likelihood values of the two parameters each obtained from the posterior distribution in the side panel.

In the text

Fig. C.1 Timing residuals of PSRs J0030+0030, J0613-0200, J0751 + 1807, J0900−3144, J1012+5307. For each pulsar, the residuals before and after subtraction of DM and monochromatic red noise are shown. The squares, circles and triangles represent P-band, L-band and S/C-band observations, respectively (see Table 2 for frequency coverage of each band). The blue/filled and black/unfilled symbols indicate the new backend data and those from EPTA DR1, respectively.

In the text

	Fig. C.2 Timing residuals of PSRs J1022+1001, J1024−0719, J1455−3330, J1600−3053, J1640+2224. Figure style is the same as Figure C.1.
In the text

	Fig. C.3 Timing residuals of PSRs J1713 +0747, J1730−2304, J1738+0333, J1744−1134, J1751−2857. Figure style is the same as Figure C.1.
In the text

	Fig. C.4 Timing residuals of PSRs J1801−1417, J1804−2717, J1843−1113, J1857+0943, J1909−3744. Figure style is the same as Figure C.1.
In the text

	Fig. C.5 Timing residuals of PSRs J1910+1256, J1911+1347, J1918−0642, J2124−3358, J2322+2057. Figure style is the same as Figure C.1.
In the text