© The Authors 2025

1 Introduction

Pulsars are rotation-powered neutron stars that emit beams of electromagnetic radiation. These beams are swept across the sky with the rotation of the stars, so that distant observers whose lines of sight are crossed by the beams might detect periodic pulses. The pulsar timing technique consists in modeling the rotation of pulsars by measuring the times of arrival (TOAs) of the pulses detected at a telescope and comparing them to the predictions from a timing model. Differences between measured and predicted TOAs, called timing residuals, contain information about the pulsars’ spin properties, orbital properties (for pulsars in multiple systems), astrometric properties, or information about the interstellar medium or their immediate environment (see, e.g., Chapter 2 of Lorimer & Kramer 2004, for a complete description). Applications of pulsar timing are numerous and diverse: for example, tests of general relativity (GR) in the strong-field regime (Kramer et al. 2021); searches for exoplanets orbiting pulsars (Nitu et al. 2024; Voisin et al. 2025); measurements of pulsar masses and constraints on the equations of state of neutron star masses using these measurements (e.g., Özel & Freire 2016); or searches for low-frequency gravitational waves using pulsar timing arrays (PTAs; see FERMI-LAT Collaboration 2022; Agazie et al. 2023; Reardon et al. 2023; EPTA Collaboration 2023b, for recent results from the different PTA collaborations).

In standard radio pulsar timing analyses, TOAs for a given pulsar are determined by comparing radio observations to a standard profile, assumed to be a noise-free representation of the pulsar’s intrinsic pulse profile shape. The standard profiles can be used to determine TOAs for arbitrarily large fractions of the recorded frequency bandwidths, assuming that the pulse profile does not evolve with frequency within these frequency intervals. Alternatively, the standard profiles can encode frequency evolution of their component(s), leading to more accurate TOAs in the case of large bandwidth observations (the latter wide-band template matching technique is for instance implemented in the PulsePortraiture software library, see Pennucci et al. 2014). Independently of the method used for constructing the template profiles, standard pulsar timing analyses crucially rely on the assumption that the template profiles are faithful representations of the pulsar pulse profiles. Reciprocally, individual observations of a given pulsar must resemble the standard profile and not vary, otherwise the cross-correlation between the individual observations and the standard profile will lead to inaccurate TOAs and introduce systematic timing errors in the TOAs.

If not properly accounted for, pulse profile variations can cause individual observations of a pulsar to differ from the chosen standard profile. In some cases, these variations originate from the pulsar itself, e.g., from pulse jitter, moding, intermittency, or relativistic spin precession (see Padmanabh et al. 2021, and references therein). Profile variations can also result from instrumental effects. In particular, inaccurate polarization calibration can distort recorded pulse profiles and in turn introduce noise in TOA datasets (see for instance van Straten 2013; Foster et al. 2015; Guillemot et al. 2023; Rogers et al. 2024). Such timing noise resulting from instrumental artifacts can degrade pulsar timing analyses, bias parameter measurements, or mimic signals expected from stochastic gravitational wave background emission (Lentati et al. 2016).

Accurate polarization calibration is therefore crucial for any pulsar timing studies, and for studying pulsar polarization properties. Recently, Dey et al. (2024) investigated different methods for calibrating pulsar observations done with the Green Bank Telescope, and found that the simplest method, the ideal feed assumption (hereafter IFA; Manchester et al. 2013) which assumes perfectly linear and orthogonal polarization feeds, led to better pulsar timing results than the more sophisticated measurement equation modeling (MEM, van Straten 2004) and measurement equation template matching (METM, van Straten 2013) methods in which polarimetric responses are respectively determined from the analysis of strongly linearly polarized pulsars observed over wide ranges of parallactic angles, or by comparing observations of a given reference pulsar with a well-calibrated template profile for the same pulsar. Dey et al. (2024) attributed the poorer performance of the other tested methods to potential instabilities of the local polarization reference source, a noise diode coupled to the receiver, and time instabilities in the pulse profile of the pulsar used as reference for the METM method. The IFA, MEM, and METM calibration methods are described in more detail later in this paper. Rogers et al. (2024), on the other hand, compared polarization calibration methods using pulsar data from the Murriyang 64m CSIRO Parkes radio telescope, and found that the METM method significantly improved pulsar timing quality over the IFA calibration scheme. Furthermore, Rogers et al. (2024) found that combining the METM data calibration method with the matrix template matching (MTM) technique (van Straten 2006), which uses all four Stokes components of the pulse profile and template profile when forming TOAs to compensate for residual inconsistencies in the polarimetric profiles, leads to significantly improved pulsar timing residuals. In their conclusions, Rogers et al. (2024) recommended that MTM be used for determining TOAs in PTA experiments. Finally, in previous work (see Guillemot et al. 2023) we significantly improved the polarization calibration of pulsar observations made with the Nançay decimetric Radio Telescope (NRT) in France by determining full polarimetric responses of the instrument from the analysis of observations of a highly linearly polarized pulsar in a special observing mode where the feed horn rotates across the observation, and by employing the MEM calibration model (details of this analysis are given in Sect. 2). The conclusion from this work is that the MEM-based calibration scheme outperforms IFA calibration for the NRT, with much improved polarimetric measurements, contrasting with the conclusions of Dey et al. (2024). In addition, similarly to Rogers et al. (2024) we also found that determining TOAs using the MTM technique significantly improved the NRT pulsar timing quality, which also suggested that this TOA extraction technique should be used in high-precision pulsar timing studies.

In this article, we present new work done following Guillemot et al. (2023) to further improve the polarization calibration of pulsar data taken with the NRT. The results from the observations of a bright highly linearly polarized pulsar with feed horn rotation presented in Guillemot et al. (2023) hinted at an hour angle-dependence of the response of the NRT. We therefore conducted an analysis aimed at determining whether the instrumental response of the NRT depends on the observed direction, and more specifically on the hour angle and declination of the observed pulsar. In Sect. 3, we examine the possibility of properly calibrating NRT data taken before the first observations, in November 2019, of bright highly polarized pulsars with horn rotation. This work, that used the METM calibration technique, led to the construction of a time- and frequency-dependent calibration model that we applied to observations of millisecond pulsars (MSPs) carried out with the NRT. In Sect. 4, we present the new polarimetric profiles for these MSPs as well as an assessment on pulsar timing quality with the calibration scheme. We finally present a summary of our work in Sect. 5.

2 Polarimetric response of the NRT: Possible dependence on the observed direction

The NRT is a meridian transit-type telescope of the Kraus/Ohio State design (see, for instance, Kraus 1960, 1966) equivalent to a parabolic dish with a diameter of approximately 94m, and capable of observing objects with declinations above ~39°. Since August 2011, NRT pulsar observations are conducted with the NUPPI backend which uses Graphics Processing Units (GPUs) to coherently dedisperse and fold the dual linear polarization signals from the receiver in real time, over a total bandwidth of 512 MHz (see, e.g., Desvignes et al. 2011; Cognard et al. 2013, for additional details on NUPPI). Pulsar observations with the NRT are primarily conducted at the central frequency of 1484 MHz. Further details on the NRT and on pulsar observations with the NUPPI backend can be found in Guillemot et al. (2023). In this new article, we consider pulsar observations made with the NUPPI data at 1.4 GHz, leaving analyses of observations at other frequencies or observations conducted with the previous backend, the Berkeley-Orléans-Nançay (BON) pulsar backend, for future work. For all NUPPI data analyses presented in this new article we made use of the PSRCHIVE software library (Hotan et al. 2004), and cleaned the data of radiofrequency interference (RFI) using the Surgical method of the CoastGuard pulsar data analysis library (Lazarus et al. 2016).

In the first article of this series (Guillemot et al. 2023), we presented an improved procedure for calibrating pulsar observations conducted with the NRT, developed in late 2019. Owing to its design, the NRT can track objects for, typically, one to two hours around transit (the maximum allowed tracking time depends on the observed declination). As a consequence, the NRT cannot sample wide ranges of parallactic angles when observing a given object. This limitation, in principle, prevents measurements of the full polarimetric response of the telescope with long observations of bright, highly polarized pulsars over wide ranges of parallactic angles (see, e.g., Stinebring et al. 1984; Johnston 2002; van Straten 2004, for examples of such experiments). Nevertheless, in November 2019 we circumvented this limitation by beginning to carry out regular observations of the bright, highly linearly polarized pulsar J0742-2822 (also known as B0740-28, see Bonsignori-Facondi et al. 1973), in a non-standard observation mode in which the feed horn rotates by ~180° over the course of the 1 h observation. During these observations, the feed horn rotation effectively mimics wide parallactic angle variations, enabling much wider ranges of pulsar polarization angles to be sampled over the hour of observation. In turn, this enabled us to determine the full polarimetric response of the telescope. The calibration solutions obtained from these rotating horn observations enabled us to obtain vastly improved polarimetric results for the selection of millisecond pulsars (MSPs) that was used in our tests, with increased signal-to-noise ratios (S/N) in most cases and polarimetric profiles that were consistent with published results, unlike those obtained with the previously used, simplistic calibration method.

The framework and main equations used for analyzing these special mode observations are described in detail in Guillemot et al. (2023) and references therein. In a nutshell, at the beginning of every pulsar observation conducted with the NRT, a 3.33 Hz noise diode is used to inject a polarized reference signal into the receiver feed horn for ten seconds. Following the noise diode observation, PSR J0742-2822 is observed for approximately 1 h, with the feed horn made to rotate by -180° between the beginning and the end of the observation. The pcm tool of PSRCHIVE is then used to carry out an MEM analysis (referred to as “Reception model” in Guillemot et al. 2023) of the noise diode observation and the observation of the pulsar over wide horn orientation angles to solve for the actual Stokes parameters of the reference noise diode signal, as well as for the values of the absolute gain G, the differential gain between the feeds, γ, the differential phase φ, and the ellipticities e and orientations θ of the two feeds. Since the solutions determined by pcm are degenerate under commutation (see Appendix B of van Straten 2004), we made two assumptions to constrain the mixing between I and V, and between Q and U: we assumed that the NRT feeds have equal ellipticities, implying that there is no mixing between I and V. We also assumed that the misorientation of the first feed is equal to zero. We made the same assumptions for the analyses presented in this new article. By determining the above-mentioned calibration parameters as well as the actual Stokes parameters of the reference noise diode signal, the MEM method makes it possible to go beyond the IFA, which assumes that the receiver feeds are perfectly orthogonally polarized and that the reference source is 100% linearly polarized, and illuminates the feeds equally and in phase. As illustrated in Guillemot et al. (2023), the analysis of observations of PSR J0742-2822 with horn rotation revealed that none of the IFA assumptions are valid in the case of NRT observations.

Although the new calibration scheme led to much improved polarimetric results, slight discrepancies between observed and modeled Stokes parameters (in particular, in the Q and U parameters) in some J0742-2822 observations with horn rotation hinted at time- (or equivalently, hour angle-) dependence of the response of the NRT. Fig. 9 of Guillemot et al. (2023) illustrates this: during the first half of the observation, observed and modeled values of the Stokes Q parameter are slightly inconsistent. The apparent time-dependence of the response of the NRT needed to be investigated as part of subsequent work. In this section we present the work done since then to understand and correct the above-mentioned inconsistencies, and test for potential dependence of the polarimetric response on certain observational parameters, such as hour angle and declination.

We began by investigating our analysis codes, and found an issue with the calculation of the time-varying parallactic for the NRT, in the analysis that led to the results presented in Fig. 9 of Guillemot et al. (2023). In the case of the NRT, the parallactic angle is given by $$\sin \left(\mathrm{\Phi }\right)=\sin \left(\delta \right)\times \sin \left(\mathrm{H}\mathrm{A}\right),$$(1)

where δ denotes the declination of the observed source, HA is the hour angle of the source at a given time, and Φ is the parallactic angle. The above expression is the same as the one given in Guillemot et al. (2023), and was the expression implemented in PSRCHIVE. However, we found that, due to the conventions used in PSRCHIVE for the signs and orientations of angles, the parallactic angle calculated using Eq. (1) needed to be used with the opposite sign in the software. We found that correcting for the sign of Φ for NRT observations in PSRCHIVE significantly improved the results of the analysis of rotating horn observations. In Fig. 1 we show measured and modeled Stokes parameters for an observation of J0742-2822 with horn rotation conducted on MJD 59368 (3 June 2021), before and after correction. As can be seen from the right panel of Fig. 1, measured and modeled values of the Stokes Q parameter in the first half of the observation are much more consistent. The fit results shown in the left panel of Fig. 1 (i.e., before correction) give a reduced χ² statistic of -2.25, while those shown in the right panel give -1.73, indicating a better agreement between the data and the model. The analysis of all other rotating horn observations of J0742-2822 also found a better agreement between the data and the model, after the sign correction. In addition to this correction we made other improvements to the analysis pipeline. First, variations in the ionospheric Rotation Measure (RM) are now taken into account, using RMextract (Mevius 2018) to determine ionospheric RM values for each sub-integration. Second, we now make use the “-Q” option of pcm to indicate that the reference noise diode is coupled after the frontend, as is indeed the case for our observations. We note, however, that the latter two modifications of the analysis pipeline compared to Guillemot et al. (2023) resulted in minor improvements to the analysis results.

Beyond the correction made to the calculation of NRT parallactic angles in PSRCHIVE that significantly improved the agreement between measured and modeled Stokes parameters for rotating horn observations, the initial goal of the study was to test for potential dependencies of the response of the NRT on observation parameters. As was noted in Guillemot et al. (2023), the absolute gain of the NRT is known to vary with the declination of the observed object (see, e.g., Theureau et al. 2005). To test for a dependence of calibration parameters on declination, and also test for a dependence on hour angle (e.g., resulting from the unusual design of the telescope), we therefore conducted observations of pulsars other than J0742-2822 in the special observation mode. We conducted a first series of rotating horn observations of pulsars at different declinations between MJDs 59368 and 59373 (3 and 8 June 2021). During this first series of observations, we sampled pulsar declinations between —28.38° (J0742-2822) and -55.73° (J0454+5543, a.k.a. B0450+55, Damashek et al. 1978). The other pulsars observed during this session were PSRs J0826+2637 (B0823+26, Craft et al. 1968), J0953+0755 (B0950+08, Pilkington et al. 1968), J1239+2453 (B1237+25, Lang 1969), J1705-1906 (B1702-19, Manchester et al. 1978), and J2113+4644 (B2111+46, Davies & Large 1970). The observations were conducted close in time (in this case, within five days), so that we could assume the instrumental response to be similar at all epochs. Table 1 lists the start epochs, declinations, durations, and start and end hour angles of these observations. Figs. 2 and 3 show the results of the MEM analysis of the individual observations and preceding observations of the reference noise diode. As in Fig. 1, left (resp., right) panels show the results obtained before (resp., after) correction of the sign of the NRT parallactic angle. Interestingly, inconsistencies between pulsars for some parameters (particularly, the differential gain γ, the feed orientation parameter θ₁ and the Stokes Q parameter of the reference signal, i.e., parameters that constrain the polarization position angle, which strongly depends on the parallactic angle correction), appear to be significantly reduced by the correction. After correction, calibration and reference signal parameters seem consistent between pulsar observations, with the exception of the absolute gain G, known to vary with declination.

At first glance, the analysis of individual pulsar observations and associated noise diode observations indicated consistent results between pulsars. In order to further test the apparent consistency of calibration parameters, we used pcm to carry out a joint analysis of all seven pulsar observations and associated noise diode observations, conducted prior to each of them. Determining the instrumental response of the telescope using multiple pulsar observations has several advantages:

• If pulsars at different declinations have been observed, potential dependencies of the instrumental response on declination can be searched for and determined.

• Strong RFI conditions can corrupt some of the observations, e.g., leading to entire frequency channels being blanked by RFI excision in these observations. If other observations are less affected then the analysis will determine the polarimetric response in more frequency channels than that of the individual, RFI-corrupted observations.

• In the case of the NRT, every pulsar observation is preceded with a short observation of a noise diode, injecting a reference signal for calibration purposes. If the pulsar observations were made closely in time, then the actual content of the reference noise diode signal can be assumed to be the same for all observations made during the session. The Stokes parameters of the noise diode can thus be determined with more precision and over more frequency channels than using individual pulsar observations.

Prior to this work, the pcm tool of PSRCHIVE did not allow multiple pulsar observations to be analyzed jointly. Therefore, in order to be able to conduct analyses of this type, we added new functionalities to pcm¹.

We conducted an MEM analysis of the dataset comprised of the seven pulsar observations with horn rotation conducted between MJDs 59368 and 59373, and of the seven noise diode observations conducted prior to each of them. Parameters fitted for were the same as in analyses of individual pulsar observations: the absolute gain G, differential gain γ, differential phase φ, feed orientation parameter θ₁, feed ellipticities ε₀ and ε₁ (with ε₀ = ε₁ by assumption), and the Stokes parameters Q, U and V of the reference signal. However, for this analysis the differential gain and phase parameters were modeled as two-dimensional polynomials, of the form: $$p(x,y)=c_0+c_1\left(y\right)x+c_2\left(y\right)x^2+\dots$$(2) $$=c_0+\sum _{i=1}^{N_x}{c}_i\left(y\right)x^i,$$(3)

with, for i ≥ 1: $$c_i\left(y\right)=c_{i,0}+c_{i,1}y+c_{i,2}{y}^2+\dots$$(4) $$=\sum _{j=0}^{N_y}{c}_{i,j}{y}^j.$$(5)

In the above expressions, p (x, y) represents the 2D polynomial, N_x is the degree of the polynomial in x and N_y is the degree of the polynomial in y, and the c terms represent coefficients to be fitted for. In our analysis, x and y respectively corresponded to the pulsar declination δ, and to the hour angle HA at a given time. We tested polynomial orders ranging from 0 (i.e., constant values) to 3 in HA, and ranging from 0 to 6 in declination. We used the same polynomial orders for both the differential gain and phase functions, to limit the number of polynomial order combinations to be tested. In the rest of the article we denote the different analyses we performed using the following nomenclature: “A/B”, meaning that the 2D polynomials representing the differential gain and phase were of order A in declination, and of order B in hour angle.

Initial results from the analysis of the rotating horn observations made between MJDs 59368 and 59373 revealed that the 0/0 model, i.e., the simplest model in which the differential gain and phase do not depend on declination and hour angle, gave the best results in terms of the test statistic we chose for comparing results for the different configurations: a “quasi-AIC” (hereafter qAIC), consisting in an average value of Akaike Information Criterion (AIC; Akaike 1973) scores over the frequency channels in which the fit was successful. A number of metrics can be used for comparing results for the various configurations, in these cases where the number of degrees of freedom varies from one configuration to another: for example, the AIC, the Bayesian Information Criterion (BIC; Schwarz 1978), the Geometric Information Criterion (GIC; Kanatani 1998), or the Stochastic Information Complexity (SIC; Rissanen 1989). However, in our pulsar observations declinations are not uniformly sampled and only N different declination values are tested, where N is the number of pulsars observed, whereas hour angles are sampled uniformly and the number of hour angles tested is much larger than N. As a consequence, strong asymmetries exist between the two dimensions, in terms of the information available to constrain model parameters. The AIC, BIC and GIC metrics, that are in principle designed to penalize model complexity by accounting for the varying numbers of degrees of freedom in the compared models, cannot account for such strong asymmetries in the sampled dimensions. The SIC metric, which we also tested, uses the logarithm of the determinant of the Hessian metrix to penalize model complexity. However, we found that the Hessian tends to be ill-conditioned in the cases where the polynomial order in declination approaches the number of pulsars available to constrain any model of variation along this dimension. Finally, additional complexity emerged from the fact that pcm analyses tended to successfully fit the available data in fewer frequency channels when polynomial orders increased. With all the above-mentioned caveats in mind, we chose to adopt the simple AIC metric in our tests, but dividing AIC scores by the number of frequency channels with successful fits, to form the “quasi-AIC” metric. In future work we will further investigate the issue of the choice of metric for comparing the results from these analyses. We also note that in addition to qAIC scores, we used pulsar observation S/N values as another metric for comparing the different configurations, as is presented later in this section.

In Fig. 4, we show the best-fit calibration parameters and Stokes parameters of the noise diode from the 0/0 analysis of the rotating horn observations made between MJDs 59368 and 59373, and those from the analysis of the observation of J0742-2822 on MJD 56368. In general the results from the two analyses are very similar. In some cases, such as those of the differential gain γ or the ellipticity parameters ε₀ and ε₁ (assumed to be equal), some small differences are found. Nevertheless, bestfit parameters from the 0/0 analysis tend to present smoother variations across the bandwidth than those from the analysis of J0742-2822 and preceding noise diode observation, with less pronounced variability and fewer outliers (see the best-fit values for the orientation parameter θ₁ for a particularly clear example). This suggests that the instrumental response inferred from the 0/0 analysis better represents the actual response than that inferred from the analysis of J0742-2822, and highlights the usefulness of jointly analyzing multiple pulsar observations at once.

Most other analyses resulted in higher qAIC scores than the simplest 0/0 model, with the exceptions of the 0/1, 1/0 and 1/1 configurations which gave very similar qAIC scores. In order to discriminate between these models, we conducted two other series of rotating horn observations of pulsars at different declinations, between MJDs 60240 (23 October 2023) and 60242 (25 October 2023), and on MJD 60461 (31 May 2024). These second and third observation sessions are described in Table 1. The pulsars observed during these two new sessions were the same as in the first, with the exception of PSR J0358+5413 (B0355+54, Manchester et al. 1972) which was added to the target list for the second session. Beyond the fact of simply taking more data, we took the opportunity of these new observing sessions to observe over longer durations than what is typically done with the NRT, and thereby observe over wider HA ranges, close to the limit that can explored with this telescope. With these longer observations, we wanted to increase our sensitivity to potential dependencies of calibration parameters on HA.

We carried out pcm analyses for each of these three observation sessions, testing the same combinations of 2D polynomial orders as described above. The results from these combined analyses are summarized in Table 2. In this table, the “Total AIC” column gives the sum of the AIC scores from the three observation sessions, the “Total n_chan” columns gives the total number of frequency channels in which the fit was successful (since the NUPPI backend splits the 512 MHz of bandwidth into 128 channels of 4 MHz each, the maximum possible value is 3 X 128 = 384), and the “qAIC” column gives the quasi-AIC scores, i.e., the total AIC scores divided by the total number of frequency channels with a successful fit. As can be seen from the table, for configurations with polynomial orders in declination larger than three, the total numbers of frequency channels with a successful fit are smaller than under the simplest configurations, and even much smaller in some cases (see, e.g., the results for the 6/3 analysis), indicating unstable fits. Quasi-AIC scores are also higher, with the exception of the 4/3 configuration which gave the lowest qAIC score among all analyses. In this analysis, however, successful fits were found in much fewer frequency channels than for the simplest analyses, so that we did not consider the 4/3 configuration further. The simplest configurations with polynomial orders in declination of up to two are found to produce the best results, both in terms of qAIC score and number of frequency channels with successful fits.

Taken at face value, Table 2 suggests that 0/1, 0/2, 0/3 and 1/1 configurations produce better results than the simplest model, the 0/0 configuration. To test whether these seemingly better configurations also constitute better calibration solutions for normal mode observations (i.e., regular NRT observations of pulsars with no horn rotation), we analyzed selections of NUPPI observations of MSPs, calibrating the observations with the parameters found by pcm for each of the three sessions and configurations. The selections consisted in MSP observations conducted within 30 days of MJDs 59370, 60240, and 60461, i.e., close in time to the three special mode observation sessions presented in Table 1. We selected normal mode MSP observations with S/N values higher than 100, to limit the sample of observations to analyze². In this work, S/N values were calculated using Eq. (7.1) of Lorimer & Kramer (2004): $$\mathrm{S}/\mathrm{N}=\frac{1}{{\sigma }_{\mathrm{O}\mathrm{f}\mathrm{f}}\sqrt{W_{\mathrm{e}\mathrm{q}}}}\sum _{i=1}^{N_{\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{s}}}(p_i-p_{\mathrm{O}\mathrm{f}\mathrm{f}}).$$(6)

Here p_i denotes the amplitude of the i-th profile bin, σ_Off and p_Off represent the off-pulse standard deviation and mean amplitude, and W_eq is the equivalent width, expressed in number of profile bins, of a top-hat pulse with the same area and peak height as the profile. After the above-mentioned selection cuts, we obtained datasets comprising 52, 55 and 57 MSP observations, respectively.

In Fig. 5, we compare the S/N values of the total intensity pulse profiles for the three MSP datasets under different calibration methods. In the left panel we compare S/N values as obtained when calibrating MSP observations with the 0/0 model, and those obtained using the IFA calibration model. The improved calibration produces increased S/N values in more than 75% of the considered observations, in line with the results of Guillemot et al. (2023). The middle and right panels show comparisons of the 0/0 model S/N values with those obtained with the 0/1 and 1/1 calibration models. S/N values obtained with the 0/0 and 0/1 models are very similar, with differences very close to 0 in most cases. Yet, more than 60% of the S/N values are higher with the 0/0 model than with the 0/1 model. As can be seen from the right panel, S/N values obtained when calibrating the data with the 1/1 model (i.e., with differential gain and phase values that depend on declination) are degraded in almost 80% of cases, with, in particular, large differences at high declination values. We also calibrated the same datasets using 0/2 and 0/3 models, finding that these models degraded S/N values in 79% and 78% of the cases, respectively. We therefore find that the 0/0 model, in which the differential gain and phase parameters do not vary with declination and hour angle, produce better results than the other models on normal mode MSP observations, in terms of S/N values. This suggests that the better qAIC values obtained for the 0/1, 0/2, 0/3 and 1/1 configurations (see Table 2) resulted from overfitting of the data, and not from actual improvements of the models. Finally, we repeated the same analyses as those presented earlier in this section, this time allowing the feed orientation parameter θ₁ and the feed ellipticities ε = ε₀ = ε₁ to vary with hour angle and declination using the same parameterization as presented in Eq. (3), and fitting constant values for the differential gain and phase parameters. The conclusions from this analysis were similar: better qAIC scores were found for a few of the tested configurations (the same configurations as in Table 2) than for the 0/0 one. However, the application of these calibration solutions to normal-mode observations of MSPs made within 30 days of MJDs 59370, 60240, and 60461 again led to degraded S/N values. We therefore take our results to mean that the polarimetric response of the NRT does not seem to depend on hour angle and declination.

Fig. 1 Left: measured and modeled Stokes parameters as a function of time and position angle Φ′ = Φ + α (where Φ denotes the parallactic angle and α is the orientation of the horn) for PSR J0742-2822, observed on MJD 59368. The pulsar was observed for ~1 h with the NUPPI backend at the central frequency of 1484 MHz. Across the observation, the feed horn was made to rotate by 180° The data points represent measured Stokes parameters; the dashed lines represent the modeled Stokes parameters as determined by analyzing the data using the MEM model (see Sect. 2 for the description of the analysis). The data were normalized by the total invariant interval (Britton 2000) to limit apparent variations in the gain caused by scintillation during the observation. The results shown in this panel were obtained by using an expression for the parallactic angle that included a sign error. Right: same results as in the left panel, after correcting for the sign error in the expression for the parallactic angle in the case of NRT observations.

Fig. 2 Left: best-fit calibration parameters as a function of channel index, as determined from the analysis of observations of a selection of pulsars (conducted between MJDs 59368 and 59373; see Table 1 for details) and of observations of the noise diode conducted prior to the pulsar observations. During the pulsar observations, the feed horn was made to rotate by 180° between the beginning and the end of the observation. The absolute gain G is expressed in units of the square root of the noise diode flux density in the considered frequency channel, $c_0=\sqrt{C_0}$. The two top panels show the feed ellipticities, ε0 and ε1, and the feed orientation parameter θ1. Right: best-fit parameters from the analysis of the same observations, after correction of the expression for the parallatic angle to be used in the case of NRT observations.

Fig. 3 Same as Fig. 2, for the Stokes parameters Q, U, and V of the noise diode. Stokes parameters are expressed in units of C0, the noise diode flux density in the considered channel.

Fig. 4 Best-fit calibration parameters as a function of channel index, from two analyses. The blue dots show the results from a combined analysis of the seven pulsar observations with horn rotation conducted between MJDs 59368 and 59373 listed in Table 1 under the 0/0 analysis configuration (see Sect. 2 for details). For comparison, the results from the analysis of one of the seven pulsar observations, in this case that of J0742-2822 on MJD 59368, are shown as orange crosses.

Fig. 5 Comparison of S/N values of normal mode 1.4 GHz observations of millisecond pulsars with NUPPI, calibrated using the ideal feed assumption (IFA) calibration model, and using improved calibration models obtained from 0/0, 0/1, and 1/1 analyses. The middle panel includes a zoomed-in image of the values near zero. MSP observations within 30 days of MJD 59370 are shown as blue circles, those conducted within 30 days of MJD 60240 are shown as orange squares, and those within 30 days of MJD 60461 are displayed as green triangles.

3 Improving polarization calibration over the entire NUPPI dataset

The analyses presented in the previous section led us to the conclusion that the polarimetric response of the NRT does not appear to depend on hour angle and declination, and that calibration solutions derived from joint observations of multiple pulsars with feed horn rotation are similar to those obtained from observations of J0742-2822 (and the noise diode) only, albeit with best-fit parameters determined over more frequency channels and with smoother variations across the available bandwidth, suggesting that they better represent the response of the NRT.

With the approximately monthly observations of J0742-2822 with horn rotation done with the NRT since November 2019, and the three observation sessions described in Table 1 conducted in 2021, 2023 and 2024, it is possible to properly calibrate the NUPPI data taken since late 2019. Nevertheless, as is described in Sect. 3.1 of Guillemot et al. (2023), between the first NUPPI observations in August 2011 and late 2019 a number of interruptions occurred over the timespan of the NUPPI dataset, some of them due to instrumental changes. These instrumental changes, as well as some of the other events, caused the polarimetric response of the NRT to evolve with time, preventing us from applying the improved calibration solutions to NUPPI data taken between 2011 and 2019.

In order to improve the polarization calibration of the earlier NUPPI observations, we employed the METM method introduced by van Straten (2013). Using a polarized reference profile, consisting in a single well-calibrated observation of a pulsar with high S/N, the METM method can be used to derive calibration solutions from observations of the same pulsar at different epochs, by comparison with the reference profile. The calibration solutions derived in this manner can then be used to calibrate other observations, made around the same epochs.

The first analysis step was thus to construct a high S/N well-calibrated polarimetric profile, of a particular pulsar. Ideally, the chosen pulsar should be bright enough for METM to be able to determine accurate calibration solutions for the individual 4 MHz frequency channels. Tests we conducted showed that this condition requires individual observations to have S/N values of several thousands at least. In addition, there need to be numerous observations of the desired pulsar, between August 2011 and the end of the considered NUPPI dataset. Two pulsars stand out, in the case of the NUPPI data archive: J0953+0755 (B0950+08), and J1136+1551 (B1133+16; Pilkington et al. 1968). Among these two objects, J0953+0755 has the highest flux density at 1.4 GHz: 100 ± 40 against 20 ± 10 mJy (Jankowski et al. 2018), and is therefore typically detected with higher S/N values with the NRT. J1136+1551, however, has been observed much more often than J0953+0755 with NUPPI: almost four times more often, so that the density of J1136+1551 observations is larger. Given that both pulsars appeared to be suitable choices for the METM analysis, we conducted analysis tests using both of them.

We selected individual 1.4 GHz NUPPI observations of J0953+0755 and J1136+1551, conducted close in time to the first session described in Table 1 and with high S/N values. The selection led us to choose the MJD 59372 observation of J0953+0755 listed in Table 1, and a normal mode observation of J1136+1551 conducted on MJD 59376 (11 June 2021). The two observations were calibrated using a modified version of the best-fit calibration solution obtained from the 0/0 analysis of the rotating horn observations made between 59368 and 59373. This modified version consisted in a version of the best-fit calibration solution in which parameters for the frequency channels with no data (e.g., due to RFI mitigation) were linearly interpolated using data from the neighboring channels. The selected and calibrated NUPPI observations were integrated in time, and the full frequency resolution of 128 channels was kept. Both observations had high S/N values: ~14 000 for J0953+0755 and ~4480 for J1136+1551.

All other 1.4 GHz NUPPI observations of J0953+0755 and J1136+1551 (and preceding observations of the noise diode) were then compared to the above-described reference polarimetric profiles, using the METM method of pcm. In addition to determining the gain G, differential gain γ, differential phase φ and the Stokes parameters of the reference noise diode, pcm fits for four additional parameters: the difference in the receiver ellipticities, ${\delta }_{\chi }={ϵ}_0-{ϵ}_1$, the sum of the ellipticities, ${\sigma }_{\chi }={ϵ}_0+{ϵ}_1$, the difference in the receiver orientations, ${\delta }_{\theta }={\theta }_0-{\theta }_1$ and the sum of the orientations, σ_θ = θ₀ + θ₁. Inspection of the analysis results for the individual observations of J0953+0755 and J1136+1551 enabled us to observe that the best-fit calibration parameters as determined from analyses of J1136+1551 observations displayed very strong variability across the bandwidth, strong variability from epoch to epoch, and large numbers of outlier values. On the other hand, best-fit results from the analysis of J0953+0755 observations displayed much smoother variations, and with small numbers of outlier values. The much better fits obtained in the case of J0953+0755 led us to choose this pulsar for the subsequent analyses, and discard J1136+1551.

In order to calibrate a normal mode NUPPI observation of a pulsar at a given epoch, one possibility is to use the calibration parameters derived using the METM method from the J0953+0755 observation made closest to that epoch. However, in further investigating the results for individual J0953+0755 observations, we noticed that the frequency variations in the best-fit parameters were consistent over time segments. More specifically, in each of these time segments, the individual parameters appeared to be scaled and/or offset versions of common functions, to be determined. Instead of calibrating pulsar observations using the calibration parameters derived from the J0953+0755 observation conducted closest in time, we thus chose to explore a different strategy: determine the boundaries of the above-mentioned time segments, determine the “archetypal” functions describing the frequency variations of the different calibration parameters over those time segments, and eventually use this information to find the calibration parameters at the considered epoch.

As mentioned above, best-fit results from the analysis of J0953+0755 were more suitable for the analysis presented here than those obtained from the analysis of J1136+1551 observations. Nevertheless, the larger density of J1136+1551 observations made the data for the latter pulsar more appropriate for finding the boundaries of the different time segments. In Fig. 6 we show three of the calibration parameters determined with the METM method as a function of time, for PSR J1136+1551, and for the 64th 4 MHz frequency channel. As can be noted from this figure, NUPPI observations of J1136+1551 have been conducted regularly over the entire dataset, with a small number of gaps that generally coincide with time intervals over which the NRT was inactive (see Sect. 3.1 of Guillemot et al. 2023, for a list of the main interruptions between 2011 and 2023). Vertical black lines indicate J0953+0755 observations, also conducted regularly since August 2011, but with lower cadence. We visually inspected METM results for J1136+1551, searching for time discontinuities in the frequency variations of best-fit parameters. In Fig. 7 we illustrate, with the example of the absolute gain parameter G, the variations seen from one time interval to another in the evolution of calibration parameters with frequency. This step enabled us to determine the boundaries of the above-mentioned time segments, indicated as dashed vertical lines in the figure.

In the following step, for each time segment we selected a single observation and used the calibration parameters for this observation as initial guesses of the “archetype” functions for the corresponding time segment. The single observations we chose as starting points had good fits and small numbers of excised frequency channels. The frequency evolution of the parameters for the chosen observations were spline smoothed to form the initial guesses of the archetypes. Then, the parameters for the other observations in the time segments were fitted to their corresponding archetypes. The gain parameter G was fitted using the following function: $$p^{\mathrm{\prime }}\left(f_i\right)=a\times p\left(f_i\right).$$(7)

Here p(f_i) represents the value of the archetype function at frequency f_i, p′ is the value of the parameter to be fitted and a is a scale parameter. The σ_θ and σ_χ parameters, found to be close to zero in most channels and observations, were fitted to the function $$p^{\mathrm{\prime }}\left(f_i\right)=p\left(f_i\right)+b,$$(8)

where b is an offset parameter. Finally, the other parameters (γ, φ, δ_χ, δ_θ, and Stokes parameters Q, U, and V for the reference noise diode) were fitted as $$p^{\mathrm{\prime }}\left(f_i\right)=a\times p\left(f_i\right)+b,$$(9)

that is, using both a scale and an offset parameter.

In the following analysis step, we refined the archetypes for each parameter and time segment, by forming weighted averages of the calibration parameters within each segment, correcting for the scale and offset parameters a and b for each observation determined from the previous step. This analysis step enabled us to obtain improved parameter archetypes, that in some cases included values for more frequency channels, i.e., the frequency channels that were absent in the chosen starting points but present in other observations from the same time segment. The improved parameter archetypes were then used to determine new estimates of the scale and offset parameters a and b using Eq. (7) to (9), for each calibration parameter and time segment. In Fig. 8, we plot the values of the calibration parameters as a function of frequency for two of the time segments. Parameter values were multiplied and/or offset using the a and b factors determined from this iteration of the analysis. The scaled and offset calibration parameters match the corresponding archetypes well, validating the construction of the archetypes and the measurement of the a and b parameters for each observation. Fig. 8 illustrates the fact that calibration parameters have common frequency evolution properties in each time segment, and that the frequency dependence of the parameters change significantly in other time segments.

As a final step in this procedure, we modeled the time variations of the scale and offset parameters determined in the previous step as Gaussian processes. We used the GPflow library (Matthews et al. 2017), and modeled the time variations of the different parameters using squared exponential kernels. In Fig. 9, we display the best-fit Gaussian processes we obtained for the scale and offset parameters for the Stokes U parameter of the reference noise diode. Absolute values of the scale and offset parameters across the NUPPI dataset strongly depend on the shapes of the archetypes, and are therefore difficult to compare between time segments. However, it can be noted from Fig. 9 that the time variations of the offset and scale parameters are generally well captured by the best-fit Gaussian processes. The modeling of time variations of the offset and scale parameters of the other calibration parameters also gave satisfactory results. In the next section, we present results from the application of calibration solutions predicted by the Gaussian processes obtained from this analysis, to NUPPI data on a sample of MSPs.

Fig. 6 Gain, differential gain and differential phase parameters as a function of time, as determined from the METM analyses of 1.4 GHz NUPPI data on PSR J1136+1551. The best-fit parameters displayed as blue dots are those for the 64th 4 MHz channel recorded by NUPPI, which has a central frequency of 1482 MHz. The solid vertical lines indicate PSR J0953+0755 observation epochs. The dashed vertical lines delimit intervals in which the instrumental response of the NRT is consistent across the bandwidth, between the epochs (see Sect. 3 for details).

Fig. 7 Top: gain parameter G as determined from the METM analysis of individual observations selected within each of the time intervals delimited in Fig. 6. Bottom: gain parameters divided by the value found for the 64th frequency channel. While discontinuities between intervals can be readily seen from the bottom panel, we note that in certain time intervals the discontinuities were more visible in other calibration parameters.

Fig. 8 Calibration parameters as a function of frequency as determined from METM analyses of PSR J0953+0755 observations, in two different time intervals. The best-fit calibration parameters were corrected for the scale and offset factors determined from the analysis described in Sect. 3, to match the archetypes for the corresponding parameter and time segment, displayed as dashed black lines. A different color was used for each of the observations made within the MJD range specified in the title.

4 Application to a sample of MSPs

4.1 Polarimetric profiles

We tested the calibration procedure presented in the previous section by analyzing NUPPI data on a sample of MSPs, whose names and main properties are listed in Table 3. The MSPs in this list are the same as those also used as test cases in Guillemot et al. (2023), with two exceptions: we introduced PSR J0125-2327 (Morello et al. 2019), a highly stable MSP observed with the NRT since mid-2017 with high cadence, and discarded PSR J1713+0747 (Foster et al. 1993), which underwent very significant and abrupt changes in its pulse profile on MJD 59321 (see, e.g., Xu et al. 2021; Singha et al. 2021; Jennings et al. 2022), rendering profile and timing analyses unnecessarily complex for the work presented here. Most of the MSPs listed in Table 3 were included in the EPTA and IPTA second data releases (Perera et al. 2019; EPTA Collaboration 2023a) for gravitational wave searches.

We selected NUPPI observations of the MSPs at 1.4 GHz, and prepared two datasets for each MSP. The first dataset consisted in observations calibrated with the IFA method (i.e., with the SingleAxis method of pcm), that is, using calibration parameters obtained from the analysis of the preceding observations of the noise diode, and assuming ideal polarization feeds and noise diode signal. In the second dataset, observations were calibrated using parameters derived from the analysis presented in Sect. 3. In practice, the Gaussian processes and archetypes for all calibration parameters were used to predict calibration solutions at the epochs of the MSP observations. In the rest of this section we refer to this second dataset as the “METM” dataset, given that the polarization calibration was carried out using the METM method of PSRCHIVE.

In Fig. 10, we show frequency-integrated polarimetric pulse profiles for PSR J1730-2304, from the two datasets. All 1.4 GHz NUPPI observations of the pulsar are shown in the figure. Similar figures for the MSPs listed in Table 3 can be found in Appendix A (see Figs. A.1 to A.4). It is clear, from the figures, that the polarimetric profiles for the METM dataset (displayed in the right panels) are significantly more homogeneous than in the other dataset. The calibration procedure described in Sect. 3, based on the METM method, thus appears to have successfully propagated the improved calibration to the entire NUPPI dataset.

In addition to comparing polarimetric profiles we investigated the impact of the new calibration scheme on the timing of these MSPs. We followed the same procedure as described in Guillemot et al. (2023) for the extraction of TOAs from the two datasets on each pulsar: observations were integrated in time and frequency, to form eight frequency sub-bands of 64 MHz each. Then, template profiles for each of the MSPs were constructed by summing the total intensity profiles for the ten highest S/N observations, and smoothing the results. TOAs were finally extracted by comparing observations with the template profiles, using the Fourier-domain Monte Carlo method implemented in the “FDM” algorithm of the pat tool of PSRCHIVE. In Fig. 11 we show ratios of S/N and TOA uncertainty values for PSR J1730-2304. Figs. B.1 to B.3 in Appendix B show the same figures for the other MSPs. Finally, in Table 3 we quote the median S/N and TOA uncertainty ratios, R₁ and R₂. In almost all tested MSPs, the new polarization calibration scheme increased the S/N values of individual NUPPI observations, leading to R₁ values generally higher than 1. The highest values of R₁ are found for PSRs J1730-2304 and J1744-1134, which are strongly polarized MSPs, thus particularly susceptible to changes in the polarization calibration. The TOA uncertainty ratios, on the other hand, have median values R₂ smaller than 1 in most cases, as a result from the generally higher S/N values in the METM dataset. It can be noted from Figs. 11 and B.1 to B.3 that the S/N and TOA uncertainty ratios all display sudden changes around MJD 58600, with, in general, higher S/N ratios and lower TOA uncertainty ratios after the event. This step coincides with the replacement of a pre-amplifier that occurred in 2019, which significantly modified the response of the NRT. The significantly different polarimetric response of the NRT is also visible in the left panels of Fig. 10 and Figs. A.1 to A.4: most of these figures show two distinct modes for the polarimetric profiles (this is particularly visible in the linear polarization profiles). The two modes correspond to the data taken before and after the pre-amplifier change.

Fig. 9 Scale and offset parameters (blue points; a and b terms from Eq. (9)) as a function of time, as determined from fits of the variations of the Stokes U of the reference noise diode with frequency, for PSR J0953+0755. The dashed red lines represent fits of the scale and offset parameters with Gaussian processes. The orange shaded areas represent 95% confidence intervals on the parameter values predicted by the Gaussian processes. Additional details on this analysis can be found in Sect. 3.

Fig. 10 NUPPI polarimetric pulse profiles for PSR J1730-2304 at 1.4 GHz. The left panel shows the profiles obtained when calibrating the NUPPI data with the ideal feed assumption (IFA) method, and the right panel shows polarimetric profiles obtained with the calibration method presented in Sect. 3, based on the METM method of PSRCHIVE. The black, red, and blue lines respectively represent the total intensity (Stokes parameter I), the linear polarization (L), and the circular polarization (V). All NUPPI observations of J1730-2304 are shown, and were normalized to the maximum value of the total intensity (see Appendix A for additional details on the construction of this figure, and for the polarimetric profiles for the other analyzed MSPs).

Fig. 11 Ratios of S/N and TOA uncertainty values for PSR J1730-2304 as a function of time. Signal-to-noise ratio (S/N) and TOA uncertainty values (σ) were derived from NUPPI observations at 1.4 GHz, calibrated using the IFA method and the calibration method presented in Sect. 3. Individual NUPPI observations were fully integrated in time and integrated in frequency so as to form eight sub-bands of 64 MHz each. TOAs whose uncertainties are shown in this figure were extracted from individual observations using the FDM algorithm (see Sect. 4.1 for additional details on the analysis, and Appendix B for equivalent figures for the other analyzed MSPs).

4.2 Timing

In addition to comparing the NUPPI polarimetric profiles, S/N values and TOA uncertainties for the MSPs in Table 3 with the improved calibration method and with the IFA method, we also investigated the influence of the new calibration method of pulsar timing quality.

As mentioned in Sect. 4.1, TOA datasets were extracted using the FDM algorithm of pat, from NUPPI data calibrated with the IFA method on the one hand, and with the METM method on the other hand. The earlier datasets are referred to as “IFA + FDM” datasets in the rest of the article, while the latter datasets are referred to as “METM + FDM”. Additionally, we extracted TOAs from the same datasets using the MTM technique also implemented in pat, and using versions of the same template profiles that included all four Stokes parameters. The Stokes parameters are used by the MTM technique to model a transformation between the template profile and the considered observations, determining TOAs that are less likely to be biased by distorsions of the pulse profile caused by instrumental artifacts. Doing this, we obtained two other TOA datasets for each MSP. We refer to these two other TOA datasets as “IFA + MTM” and “METM + MTM”.

Timing analyses were conducted using the TEMPO2 (Hobbs et al. 2006) and ENTERPRISE (Ellis et al. 2020) pulsar timing packages. ENTERPRISE was run using the PTMCMC Sampler (Ellis & van Haasteren 2017), applied under a Bayesian approach to estimate posterior distributions for the sampled parameters. Topocentric TOAs were converted to Barycentric Coordinate Time (TCB) using the DE440 solar system ephemeris (Park et al. 2021). Best-fit timing models for the MSPs in Table 3 from the EPTA Data Release 2 (EPTA Collaboration 2023a) were used as starting points for our timing analyses. For each pulsar and TOA dataset, we determined two white noise parameters: an error scale factor (EFAC, Ef ) and an error added in quadrature (EQUAD, Eq), such that $${\sigma }_i^{\mathrm{\prime }}=\sqrt{E_f^2\times {\sigma }_i^2+E_q^2},$$(10)

where σ_i denotes the uncertainty on the i-th TOA, and σ′_i is the EFAC- and EQUAD-corrected TOA uncertainty. We also modeled an achromatic red noise component, in the form of a power law in the power spectral density domain: $$S_{\mathrm{R}\mathrm{N}}\left(f\right)=\frac{A_{\mathrm{R}\mathrm{N}}^2}{12{\pi }^2}{\left(\frac{f}{1{\mathrm{y}\mathrm{r}}^{-1}}\right)}^{-{\gamma }_{\mathrm{R}\mathrm{N}}}{\mathrm{y}\mathrm{r}}^3.$$(11)

Here f represents the frequency of the signal, A_RN is the amplitude of the signal, and γ_RN is its spectral index. We sampled this spectrum using 30 Fourier coefficients.

Finally, we modeled Dispersion Measure (DM) variations using a chromatic red noise component, also defined as a power law in the Power Spectral Density (PSD) domain. The amplitude and spectral index are denoted by A_DM and γ_DM, respectively. Unlike achromatic red noise, DM variations depend on the observing radio frequency. Specifically, the delay induced in the arrival time is given by $\mathrm{\Delta }{t}^{\mathrm{D}\mathrm{M}}\propto {\nu }^{-2}$, where ν represents the radio frequency associated with each TOA. The latter spectrum was sampled using 100 Fourier coefficients. In addition to modeling the above-mentioned white and red noise parameters, we marginalized timing model parameters, following the approach described in Chalumeau et al. (2022). The E_f, γ_RN and γ_DM parameters were sampled using uniform prior ranges, whereas log-uniform prior ranges were used for E_q, A_RN and A_DM parameters.

In Figs. 12–15, we show the values of the weighted rms residuals, W_rms, and of the reduced χ² obtained after the fits of the different datasets, as well as the best-fit E_f, log₁₀ E_q, log₁₀ A_RN, γ_RN, log₁₀ A_DM, and γ_DM parameters for each pulsar and dataset, where the best-fit values and associated uncertainties here correspond to the medians and 68% credible intervals of posterior distributions. Table 4 lists the median values of each of the parameters, for each TOA dataset. Finally, Table C.1 gives the full list of best-fit parameters found for each pulsar and TOA dataset. In the latter two tables as well as in Fig. 12, the values of W_rms and ${\chi }_r^2$ were determined using TOA uncertainties that were not corrected for the best-fit EFAC and EQUAD terms, since, for example, EFAC- and EQUAD-corrected ${\chi }_r^2$ values are equal to 1 by construction.

We note, first of all, that for all tested MSPs but one (PSR J1744-1134), the lowest W_rms values are found for the METM+MTM dataset. Interestingly, the lowest W_rms found in the case of PSR J1744-1134 is for METM+FDM, i.e., with the improved calibration and without using the MTM algorithm for extracting TOAs. The fact that, as was noted in Rogers et al. (2024), MTM is theoretically expected to perform worse than standard template matching for this highly polarized pulsar due to the multiple correlation between the phase shift and the parameters of the Jones transformation matrix may explain the lower W_rms value for this TOA dataset than for the others. We also note that in almost all cases, lower ${\chi }_r^2$ values are obtained for METM+FDM than for IFA+FDM datasets, and lower ${\chi }_r^2$ values are found for METM+MTM than for IFA+MTM datasets.

Best-fit values of E_f are, in general, close to unity for all pulsars and TOA datasets. The lowest median E_f value (see Table 4) is found for the METM+MTM dataset, albeit with a small range of values for the four dataset types. Regarding the EQUAD parameter, on the other hand, we find that the lowest values of log₁₀ E_q are found for METM+MTM in most cases, and the values of log₁₀ E_q for METM+MTM are compatible within error bars with the lowest value when it is not the case. This suggests that the METM+MTM datasets contain the lowest levels of additional white noise.

Detailed comparisons of the A_RN, γ_RN, A_DM, and γ_DM parameters (displayed in Figs. 14 and 15), are made difficult by the fact that the best-fit values found for these parameters often have large uncertainties. Nevertheless, we find that in some cases (e.g., those of PSRs J0125-2327, J1022+1001, J1730-2304 or J1857+0943) significantly lower levels of achromatic noise were found in the analysis of the METM+FDM, IFA+MTM and METM+MTM datasets than in that of the IFA+FDM datasets. The latter TOA datasets, extracted from NUPPI data calibrated with the IFA method and using standard template matching, therefore contain higher levels of red noise, likely caused by imperfect calibration, that the improved calibration scheme presented in Sect. 3 or the usage of the MTM method for determining TOAs attenuated at least partially.

From this analysis we conclude that the enhanced polarization calibration presented in Sect. 3 improves the quality of NRT pulsar timing, by reducing red and white noise caused by imperfect calibration. As already demonstrated in e.g. Guillemot et al. (2023) and Rogers et al. (2024), using the MTM technique for extracting TOAs from polarimetric profiles also significantly improves TOA precision. However, it is interesting to note that the best timing results in our analysis are achieved for the METM+MTM datasets, i.e., combining the improved calibration method and the MEM technique for extracting TOAs. This highlights the usefulness of accurately calibrating pulsar observations, even when the MTM algorithm is used for determining TOAs.

Fig. 12 Weighted rms residual values and reduced χ2 values from the analysis of timing data for different MSPs (see Sect. 4.2 for details on the different datasets and TOA extraction methods used).

Fig. 13 Best-fit EFAC (Ef) and EQUAD (Eq) values for different pulsars and data calibration methods (see Section 4.2 for descriptions of the EFAC and EQUAD parameters).

Fig. 14 Best-fit achromatic red noise power-law amplitudes and spectral indices for different pulsars and data calibration methods.

Fig. 15 Best-fit DM power-law amplitudes and spectral indices for different pulsars and data calibration methods.

5 Summary and conclusions

In this paper, we have presented the continuation of the work initiated in 2019, on the polarization calibration of 1.4 GHz pulsar observations recorded with the NUPPI backend of the NRT. Until November 2019, pulsar observations were calibrated using the IFA method. In Guillemot et al. (2023) we presented an improved polarization calibration method for NRT pulsar observations, based on the MEM analyses of observations of the bright, highly linearly polarized pulsar J0742-2822 in a special mode in which the feed horn rotates by ~180° across the 1 h observation. Following the work presented in Guillemot et al. (2023), we conducted three series of observations of several bright pulsars covering wide declination ranges, to test whether the polarimetric response of the NRT depends on hour angle and declination. The analysis of the data taken as part of these three series of observations was conducted with a modified version of pcm, enabling the joint analysis of multiple pulsar observations at once, and allowing the differential gain, differential phase, feed orientation, or feed ellipticity parameters to vary with hour angle and declination. By applying the best-fit calibration solutions as obtained from these analyses to standard mode observations of MSPs conducted at similar dates, we found that the simplest model, i.e., the model in which the above-mentioned calibration parameters do not vary with hour angle and declination, gave the best results. With the exception of the absolute gain, the instrumental response of the NRT therefore does not appear to depend on the observed direction.

Since the first rotating horn observations of J0742-2822 were made in late 2019, and because the instrumental response of the NRT evolved with time, a different polarization calibration procedure needed to be developed for the calibration of pulsar observations conducted between the first NUPPI observations in August 2011 and late 2019. The solution found involved METM analyses of NUPPI observations of the bright pulsars J0953+0755 and J1136+1551, observed regularly between mid-2011 and 2024. The analysis of J1136+1551 observations enabled us to determine time segments during which the polarimetric response of the NRT appeared to be consistent. Then, based on the calibration parameters determined from METM analyses of J0953+0755 observations we determined “archetypes”, representing the frequency variations of the different calibration parameters over the individual time segments. We finally modeled the time variations of the scale and offset parameters needed to be applied to the parameter archetypes to fit individual METM results on J0953+0755, using Gaussian processes. These Gaussian processes and the parameter archetypes provide a way to account for the abrupt changes in the instrumental response of the NRT. In the future we will investigate the possibility of building calibration solutions based on METM results for more pulsars than J0953+0755 only. By providing denser measurements of the calibration parameters in time, such an analysis should further improve the efficiency of this calibration method.

Finally, we applied the newly devised calibration procedure to NUPPI data on a selection of MSPs. We found that the polarimetric profiles for the selected MSPs were significantly more homogeneous than the profiles obtained with the IFA method, indicating that the new calibration procedure provides a better representation of the polarimetric response of the NRT across the NUPPI dataset. We additionally found that, in most cases, the new calibration procedure increases the S/N values of individual observations while decreasing TOA uncertainties. The timing analysis of four sets of TOAs for each of the MSPs revealed that the new calibration procedure also led to TOAs containing lower levels of red and white noise, and that the best results in terms of pulsar timing quality were achieved when combining the new calibration procedure and the MTM technique for extracting TOAs.

Overall, unlike Dey et al. (2024) but similarly to Rogers et al. (2024) we find that the MEM and METM polarization calibration methods produce significantly better results than the IFA method for NUPPI data. It should be noted that the NUPPI TOAs included in the EPTA Data Release 2 (EPTA Collaboration 2023a), which will also be included in the future IPTA Data Release 3, were extracted using the MTM algorithm, as recommended by Rogers et al. (2024), from NRT observations calibrated with the IFA method for observations conducted before late 2019, and using the MEM method for observations made since. In future pulsar timing studies with the NRT as well as in future EPTA Data Releases we will adopt the new calibration method presented in this article, in addition to extracting TOAs with the MTM method.

The improved polarization calibration and more homogeneous profiles across the NUPPI dataset should, in principle, simplify polarization studies (e.g., RM measurements, analysis of polarimetric properties of pulsars) and significantly improve the results of wide-band template matching on NUPPI observations. In future work we will investigate the latter aspects, and will also look into applying a similar technique as presented in Sect. 3 for the polarization calibration of BON data. Between the first BON observations in 2004 and 2008, no calibration of the polarization data was performed. Then, starting from 2008, observations of the noise diode were conducted prior to each BON observation, so that BON observations could then get calibrated using the IFA approach. Past PTA analyses have shown that BON data display significant noise (see, e.g., Lentati et al. 2016), that may well result from inadequate polarization calibration. We expect that an improved calibration model, similar to the one presented in Sect. 3, will diminish the noise present in the BON data, making pulsar timing analyses using BON data more sensitive and reliable.

Acknowledgements

The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from the “Programme National de Cosmologie et Galaxies” (PNCG) and “Programme National Hautes Energies” (PNHE) of CNRS/INSU, France. A.C. acknowledges financial support provided under the European Union’s Horizon Europe ERC Starting Grant “A Gamma-ray Infrastructure to Advance Gravitational Wave Astrophysics” (GIGA; 101116134). Throughout this work, we made extensive use of the ATNF Pulsar Catalogue database (Manchester et al. 2005), available at: https://www.atnf.csiro.au/research/pulsar/psrcat/. We also made extensive use of the Astropy (Astropy Collaboration 2022), GPflow (Matthews et al. 2017), Mat-plotlib (Hunter 2007), Numpy (Harris et al. 2020), and Scipy (Virtanen et al. 2020) software libraries.

Appendix A Polarimetric profiles

In this section we present 1.4 GHz NUPPI polarimetric pulse profiles for the selection of MSPs considered in our calibration tests (see Sect.3). In these profiles, shown in Figs. A.1 to A.4, the total intensity I, linear polarization L and circular polarization V are respectively shown as black, red and blue lines. Left panels show the profiles as obtained after calibrating the NUPPI data using the IFA calibration method. Right panels show the same observations calibrated using the method presented in Sect. 3, which is based on the METM method of PSRCHIVE (see van Straten 2013).

All NUPPI observations of the considered pulsars are shown in these figures, and were normalized to the maximum value of the total intensity. The individual 2048-bin profiles were plotted with a transparency factor proportional to s_i/s_max, with s_i the S/N of the i-th observation and s_max the maximum S/N value, and such that a factor of 0 would correspond to a fully transparent (i.e., invisible) profile. The profiles shown in Figs. A.1 to A.4 were not corrected for Faraday rotation caused by magnetic fields along the line-of-sight to the pulsars, across the frequency bandwidth of 512 MHz.

Fig. A.1 Polarimetric profiles for PSRs J0125-2327, J0613-0200, and J1022+1001, as measured with the NUPPI backend at 1.4 GHz. The black lines represent the total intensity (Stokes parameter I), the red lines indicate the linear polarization (L), and the blue lines represent the circular polarization (V). The left panels show the profiles obtained using the IFA calibration method, while polarimetric profiles in the right panels were formed using the calibration method presented in Sect. 3, that uses the METM method of PSRCHIVE. All NUPPI observations centered at 1.484 GHz are shown. The observations were normalized to the maximum value of the total intensity.

Fig. A.2 Same as Fig. A.1, but for PSRs J1024-0719, J1600-3053, and J1643-1224.

Fig. A.3 Same as Fig. A.1, but for PSRs J1730-2304, J1744-1134, and J1857+0943.

Fig. A.4 Same as Fig. A.1, but for PSRs J1909-3744, J2124-3358, and J2145-0750.

Appendix B Comparison of S/N values and TOA uncertainties

In Sect. 4.1 we presented a comparison of S/N values and TOA uncertainties fora selection of MSPs (see Table 3), for two different NUPPI datasets: one calibrated using observations of a noise diode conducted prior to each pulsar observations and using the IFA method, and the other one calibrated using the method presented in Sect. 3. In this Appendix we show the ratios of S/N and TOA uncertainty values plotted as a function of time, for all the MSPs used as test cases in our study. TOAs were determined using the Fourier-domain Monte Carlo method implemented in the “FDM” algorithm of pat. The ratios as shown in Figs. B.1 to B.3.

Fig. B.1 S/N and TOA uncertainty ratios as a function of time for PSRs J0125 - 2327, J0613-0200, J1022+1001, and J1024-0719. The S/N and TOA uncertainty values were derived from 1.4 GHz NUPPI observations calibrated using the IFA method and the calibration method presented in Sect. 3 (see Sect. 4.1 for details on the analysis).

Fig. B.2 Same as Fig. B.1, but for PSRs J1600-3053, J1643-1224, J1730- 2304, and J1744-1134.

Fig. B.3 Same as Fig. B.1, but for PSRs J1857+0943, J1909-3744, J2124-3358, and J2145 −0750.

Appendix C Results of the timing analysis

In Sect. 4.2 we investigated the influence of the new polarization calibration method on pulsar timing quality, by analyzing NUPPI timing data on a selection of millisecond pulsars. In Table C.1 we provide the full list of best-fit parameters found for each pulsar and TOA dataset. Full descriptions of the analysis and of datasets can be found in Sect. 4.2.

References

All Tables

All Figures

Fig. 1 Left: measured and modeled Stokes parameters as a function of time and position angle Φ′ = Φ + α (where Φ denotes the parallactic angle and α is the orientation of the horn) for PSR J0742-2822, observed on MJD 59368. The pulsar was observed for ~1 h with the NUPPI backend at the central frequency of 1484 MHz. Across the observation, the feed horn was made to rotate by 180° The data points represent measured Stokes parameters; the dashed lines represent the modeled Stokes parameters as determined by analyzing the data using the MEM model (see Sect. 2 for the description of the analysis). The data were normalized by the total invariant interval (Britton 2000) to limit apparent variations in the gain caused by scintillation during the observation. The results shown in this panel were obtained by using an expression for the parallactic angle that included a sign error. Right: same results as in the left panel, after correcting for the sign error in the expression for the parallactic angle in the case of NRT observations.

In the text

Fig. 2 Left: best-fit calibration parameters as a function of channel index, as determined from the analysis of observations of a selection of pulsars (conducted between MJDs 59368 and 59373; see Table 1 for details) and of observations of the noise diode conducted prior to the pulsar observations. During the pulsar observations, the feed horn was made to rotate by 180° between the beginning and the end of the observation. The absolute gain G is expressed in units of the square root of the noise diode flux density in the considered frequency channel, $c_0=\sqrt{C_0}$. The two top panels show the feed ellipticities, ε0 and ε1, and the feed orientation parameter θ1. Right: best-fit parameters from the analysis of the same observations, after correction of the expression for the parallatic angle to be used in the case of NRT observations.

In the text

	Fig. 3 Same as Fig. 2, for the Stokes parameters Q, U, and V of the noise diode. Stokes parameters are expressed in units of C0, the noise diode flux density in the considered channel.
In the text

Fig. 4 Best-fit calibration parameters as a function of channel index, from two analyses. The blue dots show the results from a combined analysis of the seven pulsar observations with horn rotation conducted between MJDs 59368 and 59373 listed in Table 1 under the 0/0 analysis configuration (see Sect. 2 for details). For comparison, the results from the analysis of one of the seven pulsar observations, in this case that of J0742-2822 on MJD 59368, are shown as orange crosses.

In the text

Fig. 5 Comparison of S/N values of normal mode 1.4 GHz observations of millisecond pulsars with NUPPI, calibrated using the ideal feed assumption (IFA) calibration model, and using improved calibration models obtained from 0/0, 0/1, and 1/1 analyses. The middle panel includes a zoomed-in image of the values near zero. MSP observations within 30 days of MJD 59370 are shown as blue circles, those conducted within 30 days of MJD 60240 are shown as orange squares, and those within 30 days of MJD 60461 are displayed as green triangles.

In the text

Fig. 6 Gain, differential gain and differential phase parameters as a function of time, as determined from the METM analyses of 1.4 GHz NUPPI data on PSR J1136+1551. The best-fit parameters displayed as blue dots are those for the 64th 4 MHz channel recorded by NUPPI, which has a central frequency of 1482 MHz. The solid vertical lines indicate PSR J0953+0755 observation epochs. The dashed vertical lines delimit intervals in which the instrumental response of the NRT is consistent across the bandwidth, between the epochs (see Sect. 3 for details).

In the text

Fig. 7 Top: gain parameter G as determined from the METM analysis of individual observations selected within each of the time intervals delimited in Fig. 6. Bottom: gain parameters divided by the value found for the 64th frequency channel. While discontinuities between intervals can be readily seen from the bottom panel, we note that in certain time intervals the discontinuities were more visible in other calibration parameters.

In the text

Fig. 8 Calibration parameters as a function of frequency as determined from METM analyses of PSR J0953+0755 observations, in two different time intervals. The best-fit calibration parameters were corrected for the scale and offset factors determined from the analysis described in Sect. 3, to match the archetypes for the corresponding parameter and time segment, displayed as dashed black lines. A different color was used for each of the observations made within the MJD range specified in the title.

In the text

Fig. 9 Scale and offset parameters (blue points; a and b terms from Eq. (9)) as a function of time, as determined from fits of the variations of the Stokes U of the reference noise diode with frequency, for PSR J0953+0755. The dashed red lines represent fits of the scale and offset parameters with Gaussian processes. The orange shaded areas represent 95% confidence intervals on the parameter values predicted by the Gaussian processes. Additional details on this analysis can be found in Sect. 3.

In the text

Fig. 10 NUPPI polarimetric pulse profiles for PSR J1730-2304 at 1.4 GHz. The left panel shows the profiles obtained when calibrating the NUPPI data with the ideal feed assumption (IFA) method, and the right panel shows polarimetric profiles obtained with the calibration method presented in Sect. 3, based on the METM method of PSRCHIVE. The black, red, and blue lines respectively represent the total intensity (Stokes parameter I), the linear polarization (L), and the circular polarization (V). All NUPPI observations of J1730-2304 are shown, and were normalized to the maximum value of the total intensity (see Appendix A for additional details on the construction of this figure, and for the polarimetric profiles for the other analyzed MSPs).

In the text

Fig. 11 Ratios of S/N and TOA uncertainty values for PSR J1730-2304 as a function of time. Signal-to-noise ratio (S/N) and TOA uncertainty values (σ) were derived from NUPPI observations at 1.4 GHz, calibrated using the IFA method and the calibration method presented in Sect. 3. Individual NUPPI observations were fully integrated in time and integrated in frequency so as to form eight sub-bands of 64 MHz each. TOAs whose uncertainties are shown in this figure were extracted from individual observations using the FDM algorithm (see Sect. 4.1 for additional details on the analysis, and Appendix B for equivalent figures for the other analyzed MSPs).

In the text

	Fig. 12 Weighted rms residual values and reduced χ2 values from the analysis of timing data for different MSPs (see Sect. 4.2 for details on the different datasets and TOA extraction methods used).
In the text

	Fig. 13 Best-fit EFAC (Ef) and EQUAD (Eq) values for different pulsars and data calibration methods (see Section 4.2 for descriptions of the EFAC and EQUAD parameters).
In the text

	Fig. 14 Best-fit achromatic red noise power-law amplitudes and spectral indices for different pulsars and data calibration methods.
In the text

	Fig. 15 Best-fit DM power-law amplitudes and spectral indices for different pulsars and data calibration methods.
In the text

Fig. A.1 Polarimetric profiles for PSRs J0125-2327, J0613-0200, and J1022+1001, as measured with the NUPPI backend at 1.4 GHz. The black lines represent the total intensity (Stokes parameter I), the red lines indicate the linear polarization (L), and the blue lines represent the circular polarization (V). The left panels show the profiles obtained using the IFA calibration method, while polarimetric profiles in the right panels were formed using the calibration method presented in Sect. 3, that uses the METM method of PSRCHIVE. All NUPPI observations centered at 1.484 GHz are shown. The observations were normalized to the maximum value of the total intensity.

In the text

	Fig. A.2 Same as Fig. A.1, but for PSRs J1024-0719, J1600-3053, and J1643-1224.
In the text

	Fig. A.3 Same as Fig. A.1, but for PSRs J1730-2304, J1744-1134, and J1857+0943.
In the text

	Fig. A.4 Same as Fig. A.1, but for PSRs J1909-3744, J2124-3358, and J2145-0750.
In the text

	Fig. B.1 S/N and TOA uncertainty ratios as a function of time for PSRs J0125 - 2327, J0613-0200, J1022+1001, and J1024-0719. The S/N and TOA uncertainty values were derived from 1.4 GHz NUPPI observations calibrated using the IFA method and the calibration method presented in Sect. 3 (see Sect. 4.1 for details on the analysis).
In the text

	Fig. B.2 Same as Fig. B.1, but for PSRs J1600-3053, J1643-1224, J1730- 2304, and J1744-1134.
In the text

	Fig. B.3 Same as Fig. B.1, but for PSRs J1857+0943, J1909-3744, J2124-3358, and J2145 −0750.
In the text