11institutetext:  Laboratoire de Physique et Chimie de l’Environnement et de
l’Espace, Université d’Orléans / CNRS, 45071 Orléans Cedex 02, France
11email<EMAIL_ADDRESS>22institutetext: Observatoire
Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université
d’Orléans, CNRS, 18330 Nançay, France 33institutetext: Institute for Radio
Astronomy & Space Research, Auckland University of Technology, Private Bag
92006, Auckland 1142, New Zealand 44institutetext: Manly Astrophysics,
15/41-42 East Esplanade, Manly, NSW 2095, Australia 55institutetext: LUTH,
Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92195
Meudon, France 66institutetext: GEPI, Observatoire de Paris, Université PSL,
Université Paris Cité, CNRS, 92195 Meudon, France

# Improving pulsar polarization and timing measurements with the Nançay Radio
Telescope

L. Guillemot 1122 I. Cognard 1122 W. van Straten 3344 G. Theureau 112255 E.
Gérard 66

(Received 26 May 2023 / Accepted 2 August 2023)

###### Abstract

Context. Accurate polarimetric calibration of the radio pulse profiles from
pulsars is crucial for studying their radiation properties at these
wavelengths. Additionally, inaccurate calibration can distort recorded pulse
profiles, introducing noise in time of arrival (TOA) data and thus degrading
pulsar timing analyses. One method for determining the full polarimetric
response of a given telescope is to conduct observations of bright polarized
pulsars over wide ranges of parallactic angles, to sample different
orientations of their polarization angle and in turn determine the cross-
couplings between polarization feeds.

Aims. The Nançay decimetric Radio Telescope (NRT) is a 94m equivalent meridian
telescope, capable of tracking a given pulsar for approximately one hour
around transit. The NRT therefore cannot sample wide ranges of parallactic
angles when observing a given pulsar, so until late 2019 the polarimetric
calibration of 1.4 GHz pulsar observations with the NRT was rudimentary. We
therefore aimed to develop a method for improving the calibration of NRT
observations, overcoming the above-mentioned limitation. Ultimately, our goal
was to improve the quality of NRT pulsar timing, with better calibrated pulsar
pulse profiles.

Methods. In November 2019, we began conducting regular observations of the
bright and highly linearly polarized pulsar PSR J0742$-$2822, in a special
observing mode in which the feed horn rotates by $\sim 180^{\circ}$ over the
course of the one hour observation, mimicking wide parallactic angle
variations and in principle enabling us to determine the polarimetric response
of the NRT at 1.4 GHz. In addition, we assessed the quality of the NRT timing
of a selection of millisecond pulsars (MSPs), namely, J1730$-$2304,
J1744$-$1134, and J1857+0953, with conventional TOAs extracted from total
intensity pulse profiles, and TOAs extracted with the Matrix Template Matching
(MTM) technique, designed to compensate for putative polarimetric calibration
errors.

Results. From the analysis of the rotating horn observations of PSR
J0742$-$2822 we could determine the cross-couplings between the polarization
feeds and also constrain the Stokes parameters of the noise diode signal,
which prior to this work was erroneously assumed to be ideal and was used as
the only reference source for the calibration of pulsar observations. The
improved polarimetric response of the NRT as determined from these
observations was applied to observations of a selection of MSPs with published
polarimetric properties. We find that the new polarimetric profiles and
polarization position angles are consistent with previous findings, unlike NRT
polarimetric results obtained with the previously used method of calibration.
The analysis of the timing data shows that the new calibration method improves
the quality of the timing, and the MTM method proves very effective at
reducing noise from imperfect calibration. For pulsars with sufficient degrees
of polarization, the MTM method appears to be the preferred method of
extracting TOAs from NRT observations.

###### Key Words.:

pulsars: general, pulsars: individual (PSRs J0613$-$0200, J0742$-$2822,
J1022+1001, J1024$-$0719, J1600$-$3053, J1643$-$1224, J1713+0747,
J1730$-$2304, J1744$-$1134, J1857+0943, J1909$-$3744, J2124$-$3358),
polarization

##  1 Introduction

Pulsars are rapidly rotating, highly magnetized neutron stars that emit beams
of electromagnetic radiation. These beams are swept across the sky as they
rotate, so that the emission appears to be “pulsed” to distant observers whose
lines of sight are crossed by the beams. Pulsars can be used as tools to probe
a variety of astrophysical topics, via the “pulsar timing” technique. Pulsar
timing consists of measuring pulse times of arrival (TOAs), which correspond
to epochs at which a fiducial phase of the pulsar’s periodic signal is
detected at a telescope. The TOAs of a given pulsar encode much information
about its spin properties, its orbital properties if the pulsar belongs to a
multiple system, information about the intervening interstellar medium, and
other characteristics of the pulsar and its environment (see, for instance,
Taylor, 1992; Stairs, 2003; Lorimer & Kramer, 2004, for reviews of pulsar
timing and its applications). These properties can be determined by minimizing
the differences between the TOAs and the predictions from a model for the
pulsar encoding the above-listed information, the so-called “timing
residuals”. For some pulsars with high rotational stability, timing residuals
with standard deviations on the order of 100 ns can be achieved (see e.g., Liu
et al., 2020).

The Nançay Radio Telescope (NRT) is a meridian transit-type telescope of the
Kraus/Ohio State design (see e.g., Kraus, 1960, 1966). It consists of a
tiltable flat primary mirror (with dimensions 200 $\times$ 40m), a fixed
spherical secondary mirror (300 $\times$ 35m) separated from the primary
mirror by 460m, and a receiver carriage that follows the focal point during
observations by moving on a 80m train track. The NRT is equivalent to a
parabolic dish with a diameter of about 94m, and can track objects with
declinations above $\sim-39^{\circ}$ for approximately one hour around
transit. It is equipped with two cryogenically cooled receivers that
respectively cover the $\left[1.1;1.8\right]$ GHz and $\left[1.7;3.5\right]$
GHz frequency ranges, allowing for measurement of the four Stokes parameters
over those ranges. Since August 2011, pulsar observations with the NRT have
been conducted using the NUPPI backend, a version of the Green Bank Ultimate
Pulsar Processing Instrument used at the Green Bank Telescope (DuPlain et al.,
2008), and designed for the NRT. The NUPPI instrument (for further details,
see e.g., Desvignes et al., 2011; Cognard et al., 2013) 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. The 512 MHz bandwidth is split into 128 channels of 4 MHz each, and
in practice the bulk of pulsar observations with the NUPPI backend are
conducted at a central frequency of 1484 MHz, to exploit as much of the
frequency interval covered by the 1.1–1.8 GHz receiver as possible. A smaller
fraction of the observations is conducted with the high-frequency receiver at
a central frequency of 2539 MHz, or at other central frequencies, depending on
the pulsar or science objective. Brief descriptions of the instruments used
for conducting pulsar timing observations with the NRT prior to NUPPI can be
found in Desvignes et al. (2016).

Pulsar timing data recorded with the NUPPI backend feature prominently in
several recent high-precision pulsar timing studies. Examples include the
timing of PSR J0737$-$3039A in the so-called “double pulsar” system, enabling
highly precise tests of General Relativity (GR) in the strong-field regime
(Kramer et al., 2021), tests of the strong equivalence principle of GR via the
timing of PSR J0337+1715 in a stellar triple system (Voisin et al., 2020), the
measurement of the mass of PSR J2222$-$0137 in a massive neutron star – white
dwarf binary system (Cognard et al., 2017; Guo et al., 2021), the timing of
the eclipsing black widow pulsar PSR J2055+3829 discovered at Nançay
(Guillemot et al., 2019), or searches for low-frequency gravitational waves
using Pulsar Timing Arrays (PTAs; see e.g., Chen et al., 2021). All of these
studies analyzed TOAs that had been extracted from NUPPI data by cross-
correlating a standard pulse profile for the pulsar of interest with
individual observations, partially averaged in frequency and/or time. A
comprehensive review of the TOA extraction process and of sources of
measurement uncertainties can be found in e.g., Verbiest & Shaifullah (2018),
and a recent comparison of TOA creation practices using data from the NRT and
other telescopes can be found in Wang et al. (2022). One crucial assumption
made when extracting TOAs from timing data by cross-correlating a standard (or
template) pulse profile with observations of a given pulsar is that the
observations resemble the standard profile; that is, they are scaled versions
of the template profile on top of varying white noise levels. Distortions of
the observed pulse profile caused by time-varying instrumental artifacts can
cause individual observations to differ from the template profile, affecting
the accuracy of the TOAs and introducing systematic timing errors in the TOA
dataset. As a consequence, the individual pulsar timing observations need to
be carefully calibrated in polarization. Moreover, accurate polarization
properties of a given pulsar can provide a useful insight into the physical
processes at work in its magnetosphere (see e.g., Philippov & Kramer, 2022).

As will be presented in further detail in Sect. 2, until recently the
procedure followed for calibrating pulsar observations with the NRT was
limited. With seemingly correct pulsar timing results with the NRT (as e.g.,
determined in the context of multi-telescope pulsar timing studies), and a
relative paucity of pulsar polarization studies with this telescope (see for
instance Theureau et al., 2011; Dyks et al., 2016; Desvignes et al., 2022, for
examples of NRT pulsar polarization measurements), the validity of the former
polarimetric calibration procedure had not been studied in detail prior to
this work. However, detailed comparisons of NRT polarization profiles with
reference results obtained at other telescopes, and significant modifications
in the observed polarization properties of pulsars concurrent with
instrumental changes in 2019, prompted us to revisit the procedure for
calibrating pulsar observations with the NRT, and the procedure for extracting
TOAs from NRT data. In Sect. 2 we present efforts to improve the polarimetric
calibration procedure. In Sect. 3 we present the 1.4 GHz NUPPI data archive,
and assess the impact of the improved calibration procedure on the quality of
the timing, investigating a novel use of the Matrix Template Matching
technique (MTM; see van Straten, 2006) with NRT data. This technique uses the
timing information in all four Stokes components of the pulsar signal while
compensating for residual polarimetric calibration errors when extracting TOAs
from observations. Finally, we present a summary of our findings and some
prospects from this study in Sect. 4.

##  2 Observations and polarimetric calibration

###  2.1 First order polarimetric calibration of NUPPI observations

In August 2011, the NUPPI backend became the primary pulsar instrument in
operation at the NRT, replacing the Berkeley-Orléans-Nançay (BON) pulsar
backend. As mentioned in the introduction, the majority of NRT pulsar timing
observations are conducted at a central frequency of 1.484 GHz with the low-
frequency receiver. In this work we focus on these observations, leaving
observations with previous backends or with the high-frequency receiver for
future work.

At the beginning of each pulsar observation with the NRT, a 3.33 Hz noise
diode injects a polarized reference signal into the receiver feed horn for ten
seconds. An example of a noise diode observation is displayed in Figure 1. As
can be seen from the figure, the noise diode used at the NRT has three
different states: an “OFF” state, in which the diode is inactive, and two “ON”
states. In the first “ON” state (i.e., the “ON1” state), the noise diode
theoretically injects a 100% linearly polarized signal that illuminates both
feeds equally, which corresponds to 100% Stokes $U$ polarization (note that in
Figure 1 the noise diode observation is not calibrated, so the “ON1” state
appears to contain little emission along Stokes $U$, and strong circularly
polarized emission). Under the latter assumption, the signal from the noise
diode can thus be used to perform a first-order calibration of the complex
gains of the two polarization feeds, as implemented in e.g., the SingleAxis
calibration method of the PSRCHIVE software library (Hotan et al., 2004) used
for manipulating the pulsar data presented throughout this article. The “ON2”
state, on the other hand, theoretically consists of a 100% linearly polarized
signal along Stokes $-U$ and is not used by PSRCHIVE in the calibration
procedure.

In ten seconds, the noise diode is typically detected in each 4 MHz channel
with a signal-to-noise ratio (S/N) of 2000. Throughout this article, S/N
values are calculated using Equation 7.1 of Lorimer & Kramer (2004), i.e.,

$\mathrm{S/N}=\dfrac{1}{\sigma_{\mathrm{Off}}\sqrt{W_{\mathrm{eq}}}}\sum_{i=1}^{N_{\mathrm{bins}}}\left(p_{i}-p_{\mathrm{Off}}\right),$
(1)

where $p_{i}$ denotes the amplitude of the $i^{\mathrm{th}}$ bin,
$\sigma_{\mathrm{Off}}$ and $p_{\mathrm{Off}}$ represent the off-pulse
standard deviation and mean amplitude, and $W_{\mathrm{eq}}$ is the equivalent
width of a top-hat pulse with the same area and peak height as the profile,
expressed in number of bins. With an S/N of 2000, and given the
characteristics of the NRT (see Table 6 of Wu et al., 2018), we estimate,
using Equation 7.10 of Lorimer & Kramer (2004), that the noise diode increases
the receiver temperature by $\sim 20$% of the system temperature. This
significant temperature increase enables accurate determination of the
calibration parameters.

A comprehensive description of the formalism employed by PSRCHIVE for the
calibration of pulsar observations can be found in van Straten (2004). In a
nutshell, the SingleAxis method of PSRCHIVE assumes that the polarization
feeds are perfectly orthogonal and that the reference signal from the noise
diode is 100% linearly polarized and illuminates the two feeds equally and in
phase, as should be the case under the assumption made above. Under these
hypotheses, the $2\times 2$ complex-valued Jones matrix
$\boldsymbol{J}_{\mathrm{SA}}$ describing the response of the instrument can
be written as:

$\boldsymbol{J}_{\mathrm{SA}}=G\
\boldsymbol{B}_{\boldsymbol{\hat{q}}}\left(\gamma\right)\
\boldsymbol{R}_{\boldsymbol{\hat{q}}}\left(\varphi\right),$ (2)

where $G$ is the absolute gain, $\gamma$ denotes the differential gain between
the feeds, $\varphi$ is the differential phase, and
$\boldsymbol{B}_{\boldsymbol{\hat{q}}}$ and
$\boldsymbol{R}_{\boldsymbol{\hat{q}}}$ denote Lorentz boost and rotation
transformations (Britton, 2000). Under the assumptions made here, the Jones
matrix has the following simple form:

$\boldsymbol{J}_{\mathrm{SA}}=G\ \pmatrix{e}^{\gamma+i\varphi}&0\\\
0e^{-\left(\gamma+i\varphi\right)}.$ (3)

The measured coherency matrix, $\boldsymbol{\rho^{\prime}}_{\mathrm{ref}}$, of
the reference signal from the noise diode is then given by the polarization
measurement equation:

$\boldsymbol{\rho^{\prime}}_{\mathrm{ref}}=\boldsymbol{J}_{\mathrm{SA}}\
\boldsymbol{\rho}_{\mathrm{ref}}\ \boldsymbol{J}_{\mathrm{SA}}^{\dagger},$ (4)

where $\boldsymbol{\rho}_{\mathrm{ref}}$ is the coherency matrix of the input
ideal reference signal, in our case:

$\boldsymbol{\rho}_{\mathrm{ref}}=\dfrac{C_{0}}{2}\pmatrix{1}&1\\\ 11\\\ ,$
(5)

with $C_{0}$ the flux density of the noise diode at the considered frequency.
The SingleAxis method uses Equation 4 and measurements of
$\boldsymbol{\rho^{\prime}}_{\mathrm{ref}}$ at different frequencies to
determine $G$, $\gamma$, and $\varphi$ across the bandwidth recorded by NUPPI.
Figure 2 shows the measured $G$, $\gamma$, and $\varphi$ parameters as a
function of frequency, as obtained from the analysis of the noise diode
observation plotted in Figure 1.

Figure 1: Example of a noise diode observation conducted prior to a pulsar
observation with the NRT. The black line represents the total intensity, $I$,
of the signal as a function of pulse phase (the noise diode signal is
periodic, with a frequency of 3.33 Hz), and the red, magenta, and blue lines
correspond respectively to the $Q$, $U$, and $V$ Stokes parameters. Note that
the $I$, $Q$, $U$, and $V$ Stokes parameters are uncalibrated in this figure.
The data were extracted from a 10-s observation of the noise diode with the
NUPPI backend, conducted on MJD 59368 prior to an observation of PSR
J0742$-$2822 at the central frequency of 1.484 GHz. The data shown in this
figure were taken from the 46th 4-MHz channel recorded by NUPPI, corresponding
to the frequency of 1.410 GHz. Figure 2: First order calibration parameters as
determined from the SingleAxis analysis of the noise diode observation shown
in Figure 1, as a function of frequency. (Top) Differential phase, $\varphi$.
(Middle) Differential gain, $\gamma$. (Bottom) Absolute gain, $G$, expressed
in units of the square root of the flux density of the noise diode in the
considered frequency channel, $c_{0}=\sqrt{C_{0}}$.

The SingleAxis Jones matrix characterizing the instrumental response
determined from the analysis of the noise diode observation can then be used
to perform the polarimetric calibration of the pulsar observation conducted
subsequently. The Jones matrix to be used for this transformation is
$\boldsymbol{J}_{\mathrm{SA}}\
\boldsymbol{R}_{\boldsymbol{\hat{v}}}\left(\Phi\right)$, where
$\boldsymbol{J}_{\mathrm{SA}}$ is the Jones matrix describing the instrumental
response and $\Phi$ is the time-varying parallactic angle, given in the case
of the NRT by:

$\sin\left(\Phi\right)=\sin\left(\delta\right)\times\sin\left(\mathrm{HA}\right),$
(6)

where $\delta$ and HA are the declination and the (time-varying) hour angle of
the source.

In the left-hand panels of Figures 3 to 6 we show polarimetric profiles at 1.4
GHz for a selection of MSPs observed with NUPPI and calibrated using the
SingleAxis method of PSRCHIVE. For each pulsar we selected a high S/N
observation, which was then calibrated using the first order calibration
parameters (gains, differential gains, and differential phases) determined
from the analysis of the noise diode observation conducted prior to it. The
data were cleaned of radio-frequency interference (RFI) using the Surgical
method of the CoastGuard pulsar analysis library (Lazarus et al., 2016). The
original number of 2048 profile bins was reduced by a factor of up to eight to
improve readability, and the data were corrected for Faraday rotation across
the band caused by magnetic fields along the line-of-sight to the pulsars,
using the Rotation Measure (RM) values from Dai et al. (2015). The average
baseline levels for the $I$, $Q$, $U$, and $V$ Stokes parameters were set to
zero. Stokes parameters in the figures are in accordance with the conventions
described in van Straten et al. (2010). The average parallactic angles of all
the considered observations were within $\pm 1^{\circ}$, except for that of
J1730$-$2304, for which the average parallactic angle was $-1.8^{\circ}$.
Additional details of the selected observations can be found in Table 1. The
right-hand panels of Figures 3 to 6 show “reference” polarimetric profiles
from Dai et al. (2015), obtained by analyzing 20 cm data from the Parkes
telescope in Australia. These profiles were obtained by summing up multiple
observations conducted at a central frequency of 1.369 GHz, and over a
bandwidth of 256 MHz, i.e., from 1.241 to 1.497 GHz. We discarded NUPPI
frequency channels outside of this frequency range to prevent potential pulse
profile evolution with frequency from altering the integrated NUPPI profiles,
and therefore facilitate a comparison with the Dai et al. (2015) results.
Unlike NUPPI data, the 20 cm Parkes data analyzed in Dai et al. (2015) were
not coherently dedispersed. As a result, the profiles are affected by
dispersion smearing, indicated by markers at the top right. Details on
dispersion smearing and the amount of smearing for each of the considered
pulsars are given in Dai et al. (2015). In a few cases, such as those of PSRs
J0613$-$0200, J1600$-$3053, and J1909$-$3744, narrow profile structures
visible in the NUPPI profiles can be smeared out in the Parkes profiles, due
to this effect. Finally, in Figures 7 and 8 we show plots of the position
angles (PAs), $\Psi$, of the linear polarization as a function of rotational
phase for the same pulsars, calculated as:

$\Psi=\dfrac{1}{2}\arctan\left(\dfrac{U}{Q}\right).$ (7)

We note that the PAs shown in Figures 7 and 8 were obtained by correcting
Faraday rotation to infinite frequency111With the “\--aux_rm” option of the
PSRCHIVE tool “pam”., hence the fact that the “reference” PAs are different
from those shown in Dai et al. (2015).

Table 1: Properties of the 1.4 GHz NUPPI observations presented in Figures 3 to 6. The first five columns give the pulsar names and discovery articles, their rotational periods (P), and Dispersion Measure (DM) values, and the Rotation Measure (RM) values from Dai et al. (2015) used for correcting the NUPPI data for Faraday rotation. The following columns give the observation epochs and durations, the numbers of pulse phase bins of the NUPPI profiles displayed in Figures 3 to 6, and the signal-to-noise (S/N) ratios of the profiles, when calibrating the data using the SingleAxis method of PSRCHIVE (“Orig. S/N”) and with the improved calibration method presented in Section 2.2 (“New S/N”). See Section 2.2 and Equation 1 for details on the calculation of the S/N values. Pulsar | Discovery article | P | DM | RM | Epoch | Duration | $N_{\mathrm{bins}}$ | Orig. S/N | New S/N
---|---|---|---|---|---|---|---|---|---
|  | (ms) |  (pc cm-3) |  (rad m-2) | (MJD) | (s) |  |  |
J0613$-$0200 | Lorimer et al. (1995) | 3.06 | 38.77 | 9.7 | 59346 | 3192.3 | 256 | 279.1 | 316.6 (+13.4%)
J1022+1001 | Camilo et al. (1996a) | 16.45 | 10.27 |  $-$0.6 | 59339 | 3931.7 | 2048 | 6898.1 |  6384.0 ($-$7.5%)
J1024$-$0719 | Bailes et al. (1997) | 5.16 | 6.49 |  $-$8.2 | 59323 | 3901.1 | 512 | 349.1 | 454.6 (+30.2%)
J1600$-$3053 | Jacoby et al. (2007) | 3.60 | 52.33 |  $-$15.5 | 59338 | 3467.7 | 1024 | 794.8 | 841.5 (+5.9%)
J1643$-$1224 | Lorimer et al. (1995) | 4.62 | 62.40 |  $-$308.1 | 59400 | 1897.0 | 512 | 561.1 | 579.3 (+3.2%)
J1713+0747 | Foster et al. (1993) | 4.57 | 15.99 | 8.4 | 59384 | 2830.2 | 2048 | 4736.7 | 4966.0 (+4.8%)
J1730$-$2304 | Lorimer et al. (1995) | 8.12 | 9.63 |  $-$7.2 | 59368 | 2080.6 | 512 | 484.4 | 509.3 (+5.1%)
J1744$-$1134 | Bailes et al. (1997) | 4.07 | 3.14 |  $-$1.6 | 59347 | 1902.1 | 1024 | 1236.4 | 1510.6 (+22.2%)
J1857+0943 | Segelstein et al. (1986) | 5.36 | 13.30 | 16.4 | 59401 | 3039.3 | 512 | 1025.0 |  1018.3 ($-$0.7%)
J1909$-$3744 | Jacoby et al. (2003) | 2.95 | 10.39 |  $-$6.6 | 59359 | 3707.4 | 2048 | 1242.1 | 1418.1 (+14.2%)
J2124$-$3358 | Bailes et al. (1997) | 4.93 | 4.60 |  $-$5.0 | 59365 | 2738.4 | 512 | 439.7 |  425.8 ($-$3.2%)
J2145$-$0750 | Bailes et al. (1994) | 16.05 | 9.00 |  $-$1.3 | 59338 | 3013.8 | 2048 | 3004.3 | 3133.6 (+4.3%)

Figure 3: Polarimetric profiles for PSRs J0613$-$0200, J1022+1001, and
J1024$-$0719\. In each panel, the black line represents the total intensity
(Stokes parameter $I$), the red and magenta lines correspond to the Stokes
parameters $Q$ and $U$ describing the linear polarization, and the blue line
is the circular polarization (Stokes parameter $V$). For each pulsar, the
left-hand panel shows NUPPI pulse profiles calibrated with the SingleAxis
method as described in Section 2.1. The middle panel shows NUPPI polarimetric
profiles calibrated with the improved scheme presented in Section 2.2.
Reference results from Dai et al. (2015) are displayed in the right-hand
panel. In the latter panel, the marker at the top right indicates the
dispersion smearing resulting from incoherent dedispersion of the Parkes data.
The displayed pulse phase ranges were restricted to intervals with significant
emission features, to facilitate comparisons. The profiles were normalized to
the maximum value of the total intensity. See Table 1 for details on the NUPPI
observations presented in this figure.

Figure 4: Same as Figure 3, for PSRs J1600$-$3053, J1643$-$1224, and
J1713+0747.

Figure 5: Same as Figure 3, for PSRs J1730$-$2304, J1744$-$1134, and
J1857+0943.

Figure 6: Same as Figure 3, for PSRs J1909$-$3744, J2124$-$3358, and
J2145$-$0750.

Figure 7: Position angle of the linear polarization as a function of pulse
phase, for PSRs J0613$-$0200, J1022+1001, J1024$-$0719, J1600$-$3053,
J1643$-$1224, and J1713+0747. The position angles, $\Psi$, were calculated
using the $Q$ and $U$ Stokes parameters plotted in Figures 3 to 6, as
$\Psi=\frac{1}{2}\arctan{\left(\frac{U}{Q}\right)}$. Results from NUPPI
observations calibrated with the SingleAxis method are represented as orange
squares. Those obtained with the improved calibration scheme are shown as
green circles. Reference results from Dai et al. (2015) are shown as blue
crosses. For each pulsar, the pulse profile is shown in light gray for
reference.

Figure 8: Same as Figure 7, for PSRs J1730$-$2304, J1744$-$1134, J1857+0943,
J1909$-$3744, J2124$-$3358, and J2145$-$0750.

It is clear, from Figures 3 to 8, that the polarimetric profiles and linear
polarization angles derived from NUPPI data calibrated using the SingleAxis
method are highly inconsistent with the reference results from Dai et al.
(2015). The Parkes Multibeam receiver, used to collect the data presented in
the latter article, is also equipped with a noise diode injecting a 100%
linearly polarized signal, which is used to determine differential gains and
phases prior to the pulsar observations. However, in addition to the noise
diode observations, observations of the nearby, highly polarized MSP
J0437$-$4715 (Johnston et al., 1993) are regularly conducted over wide ranges
of parallactic angles in order to determine cross-couplings and polarization
leakages between the feeds, as described in van Straten (2004). This model,
developed further in Section 2.2, provides a more accurate estimate of the
instrumental response. Polarimetric profiles and PA measurements from NUPPI
data calibrated using the SingleAxis method are therefore incorrect. This
polarimetric calibration procedure was used until an improved calibration
scheme (see next section) was developed in late 2019.

###  2.2 Improving the polarimetric calibration of NUPPI observations

As demonstrated in Section 2.1, a simple polarimetric calibration of NUPPI
data that assumes perfectly orthogonally polarized feeds and an ideal, 100%
linearly polarized noise diode signal illuminating the two feeds equally and
in phase does not produce satisfactory results. Cross-coupling between the two
feeds, and a potentially non-ideal noise diode signal, can alter polarimetric
results. In order to determine these second-order calibration corrections, a
number of studies have used observations of bright, broadband-emitting
polarized pulsars over wide ranges of parallactic angles (see e.g., Stinebring
et al., 1984; Xilouris, 1991; Johnston, 2002). Doing so, different
orientations of the angle of polarization of the pulsar are sampled during the
observation, and cross-coupling of the polarization feeds can in turn be
determined. The NRT, however, can only observe sources for approximatively one
hour around transit. From Equation 6, it can be seen that in the most
favorable situation, corresponding to $\delta=90^{\circ}$, the parallactic
angle varies by only $15^{\circ}$ in an hour of observation, insufficient for
appropriately sampling pulsar polarization angles. However, the discrepancies
between NUPPI polarimetric profiles obtained with a simple SingleAxis
calibration with published results (such as the comparisons made in Section
2.1) prompted us to explore strategies for circumventing the above-mentioned
limitation of the telescope. The solution found in late 2019 is to conduct
observations of bright, polarized pulsars in a special observational mode, in
which the feed horn, which is held fixed during normal observations, is made
to rotate by 180∘ during the $\sim 1$-hr observation, i.e., at a rate of $\sim
3^{\circ}$ per minute, close to the maximum horn rotation rate allowed by the
instrument. Doing so effectively mimics a variation of the parallactic angle,
thereby enabling the NRT to sample much wider ranges of pulsar polarization
angles during the observation.

The framework for simultaneously determining the full polarimetric response of
the receiver and the actual Stokes parameters of the reference signal from the
analysis of such observations over wide parallactic angle ranges, or, as in
our case, over wide horn orientation ranges, is presented in van Straten
(2004) and is implemented in the Reception calibration model of PSRCHIVE. In
this model, the Jones matrix representing the instrumental response is given
by

$\boldsymbol{J}_{\mathrm{R}}=G\
\boldsymbol{B}_{\boldsymbol{\hat{q}}}\left(\gamma\right)\
\boldsymbol{R}_{\boldsymbol{\hat{q}}}\left(\varphi\right)\ \boldsymbol{C}.$
(8)

In this equation, the term $\boldsymbol{C}$, which represents the response of
a receiver with non-ideal feeds, is given by

$\boldsymbol{C}=\boldsymbol{\delta}_{0}\ \boldsymbol{S}\left(\theta_{0},\
\epsilon_{0}\right)+\boldsymbol{\delta}_{1}\ \boldsymbol{S}\left(\theta_{1},\
\epsilon_{1}\right),$ (9)

where $\boldsymbol{\delta}_{i}$ is the 2 $\times$ 2 selection matrix:

$\boldsymbol{\delta}_{i}=\pmatrix{\delta}_{0i}&0\\\ 0\delta_{1i}\\\ ,$ (10)

with $\delta_{ij}$ the Kronecker delta. In Equation 9, $\epsilon$ and $\theta$
represent the ellipticities and orientations of the two feeds, and
$\boldsymbol{S}\left(\theta,\
\epsilon\right)=\boldsymbol{R}_{\boldsymbol{\hat{u}}}\left(\epsilon\right)\
\boldsymbol{R}_{\boldsymbol{\hat{v}}}\left(\theta\right)$. The model
describing the noise diode observation is then
$\boldsymbol{\rho^{\prime}}_{\mathrm{ref}}=\boldsymbol{J}_{\mathrm{R}}\
\boldsymbol{\rho}_{\mathrm{ref}}\ \boldsymbol{J}_{\mathrm{R}}^{\dagger}$,
where $\boldsymbol{\rho}_{\mathrm{ref}}$ no longer corresponds to the ideal
reference signal given in Equation 5, and the model representing the pulsar
observation is

$\boldsymbol{\rho}_{\mathrm{psr}}^{\prime}(\Phi^{\prime})=\boldsymbol{J}_{\mathrm{R}}\
\boldsymbol{R}_{\boldsymbol{\hat{v}}}(\Phi^{\prime})\
\boldsymbol{\rho}_{\mathrm{psr}}\
\boldsymbol{R}_{\boldsymbol{\hat{v}}}^{\boldsymbol{\dagger}}(\Phi^{\prime})\boldsymbol{J}_{\mathrm{R}}^{\dagger},$
(11)

with $\Phi^{\prime}=\Phi+\alpha$ where $\alpha$ is the orientation of the horn
at a given time, and $\boldsymbol{\rho}_{\mathrm{psr}}$ is the coherency
matrix that is intrinsic to the observed pulsar signal. Given
$\boldsymbol{\rho^{\prime}}_{\mathrm{ref}}$ and
$\boldsymbol{\rho^{\prime}}_{\mathrm{psr}}$ for different values of
$\Phi^{\prime}$, and assuming $\boldsymbol{\rho^{\prime}}_{\mathrm{psr}}$ is
constant in all time samples, the Reception model can solve for the actual
Stokes parameters of the reference noise diode signal, and for the values of
$G$, $\gamma$, and $\varphi$ and of the ellipticies, $\epsilon$, and
orientations, $\theta$, of the two feeds. In practice, multiple on-pulse phase
bins from a given pulsar can be included as input source polarizations, to
increase the number of constraints for the fit. Finally, as described in Ord
et al. (2004), the polarimetric response determined from the above-described
analysis at a reference epoch can be used to calibrate regular pulsar
observations using the following Jones matrix:

$\boldsymbol{J}=G^{\prime}\
\boldsymbol{B}_{\boldsymbol{\hat{q}}}\left(\gamma^{\prime}\right)\
\boldsymbol{R}_{\boldsymbol{\hat{q}}}\left(\varphi^{\prime}\right)\
\boldsymbol{J}_{0}.$ (12)

In the above expression, $\boldsymbol{J}_{0}$ is the Jones matrix representing
the instrumental response at the reference epoch, and $G^{\prime}$,
$\gamma^{\prime}$, and $\varphi^{\prime}$ denote the changes in absolute gain,
differential gain, and differential phase since that reference epoch, as
determined from a SingleAxis analysis of a new observation of the noise diode,
using the actual Stokes parameters of the diode determined from the Reception
analysis. This assumes that the polarization properties of the noise diode
have not changed since the reference epoch, and similarly, the model described
in Equation 12 assumes that the ellipticities and orientations also have not
varied.

We conducted the first observation of a bright polarized pulsar with feed horn
rotation on MJD 58816 (November 29, 2019). The observed pulsar was PSR
J0742$-$2822 (B0740$-$28, Bonsignori-Facondi et al., 1973), a relatively
bright (flux density at 1.4 GHz of $26\pm 2$ mJy, see Jankowski et al., 2018),
slow-period ($P\sim 0.167$ s) pulsar with a very strong linear polarization
component: Johnston & Kerr (2018) measured a percentage of linear polarization
($L/I$) of $\sim 90$% for this pulsar, making it an ideal target for our
calibration observations. Since MJD 58816, PSR J0742$-$2822 has been observed
with NUPPI in calibration mode on a regular basis, with one to two
observations per month on average. We here present the results of the analysis
of a calibration mode observation of PSR J0742$-$2822 conducted on MJD 59368
(June 3, 2021), close in time to the regular mode NUPPI observations presented
in Figures 3 to 6. We integrated the $\sim 1$-hr observation in time to form
one sub-integration per minute of observation. The original frequency
resolution of 128 channels of 4 MHz each was kept. The original number of
phase bins per rotation of 2048 was reduced by a factor of four in order to
increase the S/N of the pulsar in the individual on-pulse phase bins and time
sub-integrations. We then used the PSRCHIVE tool “pcm” to carry out the
Reception model analysis of the observation of PSR J0742$-$2822 and of the
noise diode conducted prior to the pulsar observation, selecting 32 on-pulse
phase bins as sources of input pulsar polarization data for the analysis. As
explained in detail in Appendix B of van Straten (2004), solutions to Equation
11 are degenerate under commutation, and two assumptions need to be made to
constrain the mixing between Stokes $I$ and $V$ on the one hand and between
$Q$ and $U$ on the other hand. In our analysis we first assumed that the feeds
have equal ellipticities, which implies that the mixing between $I$ and $V$ is
zero. In that case, the circular polarization of the reference signal from the
noise diode is allowed to vary in the fit. Second, we assumed that the
misorientation of the first polarization feed is zero. Under this hypothesis,
the $Q$ parameter of the noise diode signal is also allowed to vary in the
fit. Finally, the data were normalized by the total invariant interval
(Britton, 2000) in order to limit apparent variations of the gain caused by
scintillation of the pulsar across the observation.

Figure 9: Measured and modeled Stokes parameters for PSR J0742$-$2822 as a
function of time and position angle $\Phi^{\prime}=\Phi+\alpha$, where $\Phi$
is the parallactic angle and $\alpha$ is the orientation of the horn at a
given time. The pulsar was observed at 1484 MHz with the NUPPI backend, in an
observing mode where the feed horn is made to rotate by $3^{\circ}$ per
minute, during the $\sim$1-hr observation. The modeled Stokes parameters,
represented as dashed lines, were determined by analyzing the data with the
Reception model of PSRCHIVE, enabling the determination of the polarimetric
response of the telescope. The data points shown in this plot correspond to
the Stokes parameters in the phase bin in which the pulsar is brightest, with
512 bins per rotation in this case. See Section 2.2 for the description of the
analysis.

In Figure 9 we show a plot of the Stokes parameters as a function of time for
the MJD 58816 observation of J0742$-$2822, along with the modeled Stokes
parameters represented as dashed lines. The data shown in this figure were
taken from the phase bin in which the pulsar is brightest. Stokes parameters
$Q$ and $U$, which describe the linear component of the polarized signal,
display clear sinusoidal variations, as expected given that the feed horn
rotated during the observation. The modeled Stokes parameters generally
represent the actual measurements well, with slight discrepancies
(particularly in the $Q$ and $U$ parameters) observed during the first half of
the observation. These discrepancies, which suggest possible time-dependence
of the instrumental response during observations, will be investigated further
in future work.

The gain, differential gain, and differential phase values, as well as the
ellipticities and orientation parameters as determined from the Reception
analysis, are plotted in Figure 10. The Stokes parameters of the noise diode
plotted as a function of channel index are shown in Figure 11. The analysis
finds significant non-orthogonality in most frequency channels, confirming the
presence of cross-coupling between the polarization feeds. The gain and
differential phase parameters are generally similar to those shown in Figure
1, determined from a simple SingleAxis analysis of the noise diode
observation. Differential gains, however, strongly differ from those plotted
in the latter figure in most channels. As can be seen from Figure 11, the
analysis confirms that the noise diode signal mainly consists of linear
emission along the $U$ Stokes parameter, albeit with significant Stokes $Q$
and $V$ emission in most frequency channels, indicating that the noise diode
signal is also not ideal. The ellipticities and orientations of the feeds, as
well as the Stokes parameters of the noise diode plotted in Figures 10 and 11,
need to be taken with a grain of salt since, as mentioned in the previous
paragraph, strong assumptions about the non-orthogonality parameters needed to
be made to constrain the mixing between $I$ and $V$, and between $Q$ and $U$.
The analysis nevertheless demonstrates the existence of significant cross-
coupling between the feeds and that the noise diode is most likely non-ideal.

Figure 10: Best-fit calibration parameters as determined from a joint analysis
of an observation of PSR J0742$-$2822 with the feed horn rotating by 180∘
across the $\sim 1$-hr observation, and of an observation of the noise diode
conducted prior to the pulsar observation. The three bottom panels are the
same as in Figure 2. The blue points represent the best-fit results from this
analysis, and the yellow squares show the results from the SingleAxis analysis
of the noise diode observation only, as in Figure 2, for comparison. The two
top panels show the feed ellipticities (assumed to be equal in our analysis,
see Section 2.2 for details), and the orientations ($\theta_{0}$ is assumed to
be 0 in the analysis). Figure 11: $Q$, $U$, and $V$ Stokes parameters of the
noise diode signal as a function of frequency, expressed in units of the flux
density of the diode in the considered channel, $C_{0}$. The results plotted
in this figure were determined by analyzing the noise diode observation shown
in Figure 1, as part of the procedure for improving the calibration of NUPPI
data described in Section 2.2.

The improved calibration should in principle lead to more accurate
polarimetric results. We repeated the analysis of a selection of 1.4 GHz
observations of MSPs with NUPPI presented in Section 2.1, this time using
improved calibration parameters obtained from the analysis of observations of
PSR J0742$-$2822 with feed horn rotation. For each of the MSP observations
listed in Table 1, we used the PSRCHIVE command “pac” to calibrate the data
using the second order calibration parameters from the closest-in-time
observation of PSR J0742$-$2822 and the noise diode observation conducted
before the MSP observation. Equation 12 was then used by pac to perform the
calibration. The new polarimetric profiles are displayed in Figures 3 to 6,
and the new PA measurements are shown in Figures 7 and 8. Additionally, in
Table 1 we list the profile S/N values as obtained when calibrating the data
with the simple SingleAxis calibration scheme, and with the improved method.
The polarimetric profiles formed with the improved calibration method are much
closer to the Dai et al. (2015) profiles than those obtained with the generic
calibration method, as can be seen from Figures 3 to 6. Similarly, the new PA
values are very different from those obtained with the SingleAxis calibration
method (see, for instance, the results for PSR J1713+0747 in Figure 7 or PSR
J1857+0943 in Figure 8 for particularly discrepant results), and are also very
close to the Dai et al. (2015) measurements. The latter results and those
obtained with the improved calibration method are, in all cases highlighted in
Figures 7 and 8, qualitatively consistent. The results are, however, not
identical: some differences between the PA values are for example seen in the
leading edge of the first peak of PSR J1022+1001 or in the main peak of
J1713+0747. It should be noted that the NUPPI observation of PSR J1713+0747
was carried out about 60 days after the pulsar underwent an abrupt profile
change on MJD 59321 (Xu et al., 2021; Singha et al., 2021; Jennings et al.,
2022), while the reference data were taken several years before the event. As
a consequence, the NUPPI and Parkes profiles are not expected to be identical.
Residual calibration imperfections, or inaccurate RM values at the epochs of
the NUPPI observations, could also produce discrepancies. However, the new
NUPPI measurements generally agree well with the reference results, indicating
that NUPPI polarimetric profiles calibrated with the improved scheme can be
confidently used for studying the polarization emission properties of pulsars
or determining RM values along their lines-of-sight.

As can be seen from Table 1, the new S/N values are in almost all cases larger
than those of the profiles calibrated with the SingleAxis method, with
increases of about 20% in the case of PSR J1744$-$1134, or 30% for PSR
J1024$-$0719\. One notable exception appears to be that of PSR J1022+1001, for
which the new S/N value is significantly smaller with the improved calibration
method, and for which the total intensity profile as obtained with the
SingleAxis calibration method is closer to the reference profile than that
obtained with the improved calibration. As mentioned earlier, we assumed that
the mixing between Stokes $I$ and $V$ is zero when determining the
polarimetric response. In PSR J1022+1001, the maximum of the total intensity
profile coincides with the maximum of Stokes $V$, in absolute value. We could
thus expect significant differences in Stokes $I$ in parts of the profile
where $|V|$ is large. It should also be noted that the pulse profile of PSR
J1022+1001 is known to exhibit significant time variations (see e.g., Camilo
et al., 1996b; van Straten, 2013; Padmanabh et al., 2021) over timescales of
years, and could therefore differ in the NUPPI and Parkes observations
considered here. Besides, the Stokes $Q$, $U$, and $V$ parameters as
determined with the improved scheme are more consistent with the reference
results. We therefore conclude that the new polarimetric profiles for PSR
J1022+1001 shown in Figure 3 are a better representation of the true
polarimetric properties of this pulsar, despite the lower S/N value.

The new differential gain values shown in Figure 10 are significantly
different from those derived from the analysis of the noise diode observation
alone. S/N values are calculated using the total intensity profiles, which are
formed by adding data from the two polarization feeds. Under the assumption of
ideal feeds recording independent and identically distributed (i.i.d.) data,
the noise levels in each feed have equal variance, and the S/N is proportional
to the square root of the number of polarizations summed (two in this case).
However, if data are calibrated using inaccurate differential gains and
unaccounted cross-couplings then the i.i.d. assumption no longer applies,
thereby reducing the effective number of polarizations added to less than two.
We argue that the larger S/N values in Table 1 are a consequence of the fact
that the improved calibration, and particularly the substantial corrections to
the differential gain parameter over large portions of the observed bandwidth,
restored the i.i.d. assumption in most cases, resulting in larger effective
sample sizes. Foster et al. (2015) studied the influence of inaccurate data
calibration on pulsar observations and pulsar timing, based on simulations,
finding that inaccurate calibration results in lower observed S/N than the
ideal value. In their simulations, however, noise was added after the
instrumental polarization transformation, whereas in our case the majority of
the system noise was added before the differential gain and was thus subjected
to the dominant boost transformation. Our conclusions are therefore consistent
with and complementary to those of Foster et al. (2015).

In general, the improved calibration method led to larger S/N values for the
observations we analyzed. As a consequence, NUPPI TOAs extracted from these
observations should have smaller uncertainties. Furthermore, regular NRT
observations of PSR J0742$-$2822 with horn rotation enable us to determine the
polarimetric response of the instrument on a regular basis. Therefore, we also
expect NUPPI TOAs to be more accurate, in addition to being more precise, as
pulse profile distortions from inaccurate calibration should in principle be
mitigated. Effects of the improved polarimetric calibration on NUPPI timing
results are discussed in detail in the following section.

##  3 Pulsar timing with NUPPI

###  3.1 The 1.4 GHz NUPPI timing archive

To assess the influence of the improved calibration on the quality of NUPPI
timing data, we consider three of the MSPs analyzed in Section 2: PSRs
J1730$-$2304, J1744$-$1134, and J1857+0943, all three of them being well-timed
PTA pulsars observed with high cadence since the beginning of NUPPI
operations. The rationale for selecting these three MSPs is the following: of
the MSPs considered in Section 2, PSR J1744$-$1134 is the most polarized
pulsar and should thus be particularly susceptible to calibration
inaccuracies. It is also among the best targets used by PTAs to search for
gravitational waves (see e.g., Chen et al., 2021). Conversely, PSR J1857+0943
has the lowest degree of polarization in the sample according to Dai et al.
(2015), and thus represents an interesting test case. Finally, we chose to
consider PSR J1730$-$2304, for which van Straten (2013) predicted that the MTM
technique (see Section 3.2) should determine TOAs with much ($\sim 30$%) lower
uncertainties than techniques that extract TOAs from total intensity profiles.
In this subsection we present the 1.4 GHz NUPPI timing archive using PSR
J1744$-$1134 as an illustrative example, since this pulsar has been observed
with high cadence with the NRT since NUPPI was commissioned, and with previous
backends also.

In Figure 12 we show timing residuals for PSR J1744$-$1134 as a function of
time between August 2011 and December 2022, where the residuals correspond to
the differences between the measured TOAs and the predictions from an accurate
timing model for the pulsar. From the figure it can be seen that pulsar timing
observations have been conducted at a relatively stable cadence between August
2011 and December 2022, with intervals of weeks to months when no timing
observations were carried out. We here provide a list of the main
interruptions, and what prevented observations from being carried out during
those time intervals:

  1. a)

Painting activity in July 2013 and an intense observational campaign on the
magnetar at a similar right ascension J1745$-$2900 (Eatough et al., 2013)
prevented J1744$-$1134 from being monitored normally between April and July
2013.

  2. b)

A fire in the electrical cabinet of the focal carriage and subsequent repairs
prevented the NRT from observing between late May and early July 2015.

  3. c)

Due to the failure of an amplifier of the low-frequency receiver in June 2016,
no 1.4 GHz observations could be carried out for several weeks. Incidentally,
the failure of a pump in the hydraulic system of the focal carriage in
September 2016 made it unable to move and track sources. The latter issue
required heavy repairs, and normal observations resumed in January 2017.

  4. d)

Issues with the motion of the focal carriage along the track prevented NRT
observations between August and November 2017.

  5. e)

The NRT was inactive for a few weeks in August 2018 due to electrical
problems.

  6. f)

The failure of a pre-amplifier, which was replaced after a few weeks,
prevented NRT observations with the low-frequency receiver from happening
between late March and late May 2019.

  7. g)

A full repair of the focal track in early 2021 prevented the NRT from
operating normally for a few weeks.

  8. h)

Electrical problems on the site caused the NRT to be inactive from late
December 2022 to mid-February 2023.

The above-mentioned interruptions are highlighted in the top panel of Figure
12.

Figure 12: NUPPI timing data on PSR J1744$-$1134\. (Top) Timing residuals
(i.e., differences between measured pulse arrival times and arrival times
predicted by a model for the pulsar) as a function of time for the millisecond
pulsar J1744$-$1134 observed with the Nançay Radio Telescope at the central
frequency of 1.484 GHz using the NUPPI backend. The 512 MHz of frequency
bandwidth recorded by NUPPI were split into eight sub-bands of 64 MHz each,
and the observations were integrated in time. Pulse arrival times were
extracted from the timing data using the MTM method, as further described in
Section 3. Time intervals when the NRT was inactive or when the low-frequency
receiver was unavailable are highlighted with arrows (for interruptions
shorter than 40 days) and shaded areas (when intervals were longer than 40
days). See Section 3.1 for details on the individual interruptions. (Bottom)
Average number of 1.4 GHz NUPPI observations of PSR J1744$-$1134 per month as
a function of time. The average cadence was estimated by counting the number
of observations within a time window of 100-d centered on the considered date,
and normalizing the obtained number to 30 days.

###  3.2 Timing data preparation

We selected NUPPI observations of PSRs J1730$-$2304, J1744$-$1134, and
J1857+0943 at 1.4 GHz, and cleaned the data of RFI using the Surgical method
of the CoastGuard package. We then used the cleaned observations to prepare
three distinct datasets for each pulsar. In the first dataset, the
observations were left uncalibrated. In the second, observations were
calibrated with the SingleAxis method described in Section 2.1, i.e., using
calibration parameters obtained from the analysis of the noise diode
observations conducted prior to each of the pulsar observations, and assuming
that the noise diode signal and the polarization feeds are ideal. In the third
dataset, NUPPI observations were calibrated using the improved method
described in Section 2.2. In practice, the calibration parameters from the
closest-in-time rotating horn observation of PSR J0742$-$2822 were used,
setting a maximum threshold for the time separation of 100 days. The first
observation of J0742$-$2822 with feed horn rotation was conducted on MJD
58816; therefore, for NUPPI observations of PSRs J1730$-$2304, J1744$-$1134,
and J1857+0943 conducted before MJD 58716, a generic SingleAxis calibration
was applied, as in the second dataset. Finally, all observations in the three
datasets were integrated in time and frequency, forming eight frequency sub-
bands of 64 MHz each.

In Figure 13, we compare the S/N values of the total intensity pulse profiles
from the second and third datasets. As expected, between MJDs 55800 and $\sim
58700$ the S/N values are identical, the calibration procedure being the same.
However, it is clear from the figure that the improved calibration procedure
generally leads to significantly larger S/N values as a consequence of
accurate modeling of cross-couplings and polarization leakages, and
particularly so for the highly polarized pulsar J1744$-$1134, with median
increases of 4%, 20%, and 2%, respectively, for PSRs J1730$-$2304,
J1744$-$1134, and J1857+0943.

Figure 13: Ratios of S/N and TOA uncertainty values for PSRs J1730$-$2304,
J1744$-$1134, and J1857+0943 as a function of time. The values are derived
from NUPPI pulse profiles observed at 1.4 GHz and calibrated using the
SingleAxis method (“no corr.”) and the improved method described in Section
2.2 (“corr.”). For each pulsar, the top panel (blue points) shows the ratios
of the S/N values and the bottom panel (orange points) shows the ratios of the
TOA uncertainties. Individual observations were fully integrated in time, and
integrated in frequency to form eight sub-bands. TOAs were extracted from the
observations using the MTM technique. See Section 3.2 for additional details
about the data preparation.

In standard pulsar timing analyses, TOAs are extracted from the data by
comparing total intensity profiles and a template profile for the pulsar at
the considered frequency, using a cross-correlation algorithm. Wang et al.
(2022) compared different methods for generating the template profiles, and
found that for bright pulsars such as the ones considered here, templates
formed by smoothing the addition of a selection of high S/N observations of
the pulsar were preferred over other standard template preparation methods.
Wang et al. (2022) further recommended the usage of the “FDM” method, as
implemented in the PSRCHIVE tool “pat” for extracting TOAs from the
observations. In the FDM method, the phase-gradient of the cross-power
spectrum of the template profile and the observation are fitted in the Fourier
domain (Taylor, 1992), as in the PGS method of pat, to determine the TOAs, but
with TOA uncertainties estimated from a Monte Carlo Analysis. The FDM method
was found by Wang et al. (2022) to produce TOAs and TOA uncertainty estimates
that are generally more reliable than in the other tested methods. We adopted
the recommendations of this systematic study of TOA creation procedures, and
generated template profiles for each of the three datasets on the three
pulsars by adding the total intensity profiles for the ten highest S/N
observations made after MJD 58800 and then smoothing the summed profiles. The
time selection was used to ensure that the template profile for the third
dataset be formed by adding observations calibrated with the same procedure.
TOAs for the three datasets were then extracted by comparing the observations
with the respective template profiles, using the FDM method of pat.

As mentioned above, calibration observations of PSR J0742$-$2822 have been
conducted at Nançay since MJD 58816. We therefore do not have measurements of
the full polarimetric response of the NRT before that date. In addition, as
mentioned in Section 3.1, various instrumental changes, such as the
replacement of a pre-amplifier in 2019, occurred over the time span of the
NUPPI dataset. Some of these instrumental changes likely modified the
polarimetric response of the NRT, potentially leading to varying total
intensity pulse profiles at different epochs. In addition to extracting TOAs
from total intensity profiles using the FDM method, we thus also investigated
the use of the MTM technique presented in van Straten (2006) and implemented
in pat, to extract TOAs from the NUPPI data. Unlike the FDM method, the MTM
technique uses all four Stokes parameters and models a transformation between
the template profile and the considered observation, to reduce timing errors
caused by instrumental distortions of the pulse profile. In this approach,
observed polarization profiles represented by the coherency matrix,
$\boldsymbol{\rho^{\prime}}$, are related to the template profile for the
pulsar, $\boldsymbol{\rho}_{0}$, by the equation:

$\boldsymbol{\rho^{\prime}}\left(\Phi\right)=\boldsymbol{J}\
\boldsymbol{\rho}_{0}\left(\Phi-\Delta\Phi\right)\boldsymbol{J}^{\dagger}+\boldsymbol{\rho}_{\mathrm{DC}}+\boldsymbol{\rho}_{\mathrm{noise}}\left(\Phi\right),$
(13)

where $\Phi$ is the pulse phase, $\Delta\Phi$ is the phase shift,
$\boldsymbol{\rho}_{\mathrm{DC}}$ and $\boldsymbol{\rho}_{\mathrm{noise}}$ are
the coherency matrices for the DC offset between the profiles and for the
noise, and $\boldsymbol{J}$ is the polarimetric transformation between
$\boldsymbol{\rho}_{0}$ and $\boldsymbol{\rho^{\prime}}$. By comparing the
four Stokes parameters of the observations and of the template profile in the
Fourier domain, the MTM method determines the $\Delta\Phi$ parameter and the
Jones matrix, $\boldsymbol{J}$. We generated template profiles for each of the
three datasets on the three MSPs by summing the ten highest S/N observations
carried out after MJD 58800 and smoothing the added profiles, this time
keeping all four Stokes parameters. The PSRCHIVE command pat was finally used
to extract TOAs using these template profiles and the MTM algorithm. In the
bottom panel of Figure 13 we compare the uncertainties of TOAs extracted with
the MTM method from the observations from the second and third datasets.
Before MJD $\sim 58700$, the ratios are centered on one, but are not
necessarily equal to one due to the fact that, although the procedure followed
for calibrating the observations is the same for these epochs, the template
profiles used for extracting the TOAs are different between the two datasets.
After MJD $\sim 58700$, TOA uncertainties for the third dataset are generally
smaller, as expected, given that the S/N values are generally higher with the
improved calibration scheme.

After the above-described TOA preparation steps, we obtained six different TOA
datasets for each pulsar: TOAs extracted from uncalibrated NUPPI total
intensity profiles with the FDM algorithm (we will refer to these TOA datasets
as “Uncal. FDM”), TOAs derived from uncalibrated NUPPI polarimetric profiles
with the MTM algorithm (“Uncal. MTM”), TOAs obtained after a simple SingleAxis
calibration with the FDM algorithm (“Cal. FDM”) and with the MTM method (“Cal.
MTM”), and finally TOAs obtained after application of the improved calibration
method with the FDM algorithm (“Cal. + corr. FDM”) and with the MTM method
(“Cal. + corr. MTM”). Figure 14 presents histograms of the TOA uncertainties
for the six different TOA datasets on each MSP. One can see that, as expected,
the improved polarimetric calibration scheme leads to the smallest TOA
uncertainties among the MTM and FDM datasets, while the largest TOA
uncertainties are obtained for the uncalibrated data. Additionally, for the
“Cal. + corr.” datasets, we find that the ratios of the median TOA
uncertainties estimated by the MTM and FDM algorithms are respectively 0.68,
1.49, and 0.91 for PSRs J1730$-$2304, J1744$-$1134, and J1857+0943. These
values are very close to the theoretical expectations listed in Table 1 of van
Straten (2013), of 0.71, 1.56, and 0.89, respectively. For PSR J1744$-$1134,
which is almost 100% linearly polarized (see Dai et al., 2015), MTM finds
larger TOA uncertainties due to large correlations between the phase shift and
the parameters describing the polarization transformation. As can be seen from
Table 1 of van Straten (2013), MTM is expected to find smaller TOA
uncertainties in most pulsars, because of their lower fractional
polarizations.

Figure 14: Histograms of the TOA uncertainties, with TOAs extracted from
uncalibrated NUPPI data (“Uncal.”), data calibrated using the SingleAxis
method (“Cal.”) and data calibrated with the improved calibration method
described in Section 2.2 (“Cal. + corr.”). Three of the TOA datasets were
derived from the NUPPI data using the FDM algorithm; the other three were
extracted using the MTM technique.

###  3.3 Timing results

The “Cal. + corr. MTM” TOA datasets were built from better calibrated
observations than the other datasets, and the TOAs were extracted using the
MTM algorithm, which models putative profile distortions caused by calibration
artifacts at the same time as it determines the TOAs. Of the different TOA
datasets described in Section 3.2, they thus likely represent our best
datasets for timing analyses of PSRs J1730$-$2304, J1744$-$1134, and
J1857+0943. We analyzed the three “Cal. + corr. MTM” datasets using the TEMPO2
pulsar timing package (Hobbs et al., 2006). Topocentric TOAs were converted to
Barycentric Coordinate Time (TCB) using the DE440 solar system ephemeris
described in Park et al. (2021). The timing models were purposefully kept
simple: for each pulsar we fitted for the sky coordinates, the proper motion
and timing parallax, the pulsar’s spin frequency and its first two time
derivatives, the orbital parameters in the case of PSR J1857+0943, and for the
DM and its first two time derivatives. That is, we did not model any spin
frequency or DM red noise components that could absorb timing instabilities
caused by calibration issues. Similarly, we also did not model putative
systematic errors in TOAs and/or in their uncertainties using EFAC and EQUAD
parameters (see e.g., Lentati et al., 2014). The best-fit ephemerides are
consistent with previously published results on these pulsars (e.g., Desvignes
et al., 2016). In all cases the second spin frequency derivatives were
significant, indicating the presence of spin noise in these pulsars.

Figure 15: Differences between timing residuals of PSR J1744$-$1134 for the
TOA datasets described in Section 3.2, and those for the “Cal. + corr. MTM”
TOA dataset, as a function of time. In all cases, the residuals were
determined using an ephemeris for the pulsar obtained by analyzing the “Cal. +
corr. MTM” dataset with TEMPO2. Colored symbols represent the residual
differences, and black dots show averages over 30-d sliding windows.

The timing residuals for PSR J1744$-$1134 as a function of time for the best-
fit model and the “Cal. + corr. MTM” TOA dataset are plotted in Figure 12. The
weighted rms residual is $\sim 0.628$ $\mu$s, confirming PSR J1744$-$1134 as
one of the most precisely timed pulsars known currently. The reduced
$\chi^{2}$ value is $\sim 2$, indicating that the timing data could be
whitened further than what is done by our simple ephemeris for this pulsar. We
used the best-fit model for each pulsar to determine timing residuals for the
five other TOA datasets, and compared the timing residuals with those for the
“Cal. + corr. MTM” dataset. In Figure 15 we show the differences in the timing
residuals (i.e.,, the differences between the residuals for the various TOA
datasets, and those for the “Cal. + corr. MTM” dataset) for PSR J1744$-$1134\.
As expected, the largest differences are seen for the “Uncal. FDM” dataset
(see panel a), with significant structures in the residual differences and
many data points well away from zero. These differences reflect changes in the
properties of the instrument that strongly influence the total intensity pulse
profile of J1744$-$1134 at a given epoch, and can in particular make it
inconsistent with the template profile used for extracting the TOAs.
Interestingly, much less structure is seen for the “Uncal. MTM” TOAs (panel
b), demonstrating that the MTM algorithm manages to compensate for many of the
profile changes induced by instrumental artifacts, albeit not completely,
likely due to bandwidth depolarization (van Straten, 2002) in the 8 $\times$
64 MHz frequency sub-bands. Much smaller residual differences are observed in
the case of the “Cal. FDM” dataset (panel c) than for the “Uncal. FDM” TOAs.
However, a systematic offset is detected between the data taken before and
after the pre-amplifier change at the NRT in 2019, showing that the total
intensity pulse profiles are inconsistent in these two parts of the NUPPI
archive. This offset is compensated for by the MTM algorithm, as can be seen
from panel d of Figure 15. In the latter case, the residual differences are in
all cases close to zero, indicating that for this pulsar the MTM algorithm
finds similar results, independently of the calibration scheme. Finally, it
can be seen from e that the above-mentioned offset is not seen in data
prepared with the improved polarimetric calibration scheme. However, as in
panels a and c, significant structure is seen in the residual differences,
corresponding to cases where the MTM algorithm found differences between the
observed pulse profiles and the template profile, which it compensated for.
For the other two MSPs, J1730$-$1134 and J1857+0943, differences between
timing residuals display the same behavior, albeit with more noise than in
J1744$-$1134 due to the fact that they are less precisely timed objects.

Strictly speaking, although significant structure is seen in panel e, we
cannot tell from Figure 15 if one of the two TOA datasets (“Cal. + corr. FDM”
or “Cal. + corr. MTM”) is a better dataset than the other. We therefore
conducted independent TEMPO2 analyses of all TOA datasets for the three
pulsars to obtain best-fit timing models for each dataset and in turn
determine weighted rms and reduced $\chi^{2}$ values. These parameters are
listed in Table 2. As expected, of the three TOA datasets extracted from total
intensity pulse profiles with the FDM algorithm, the lowest weighted rms
residual and reduced $\chi^{2}$ values are obtained for the dataset that is
based on the best calibration method. The same conclusion also applies to TOA
datasets derived using the MTM algorithm, although in this case, the “Cal.
MTM” and “Cal. + corr. MTM” datasets have similar weighted rms residuals and
reduced $\chi^{2}$ values. In Figure 16 we show the distributions of the
reduced $\chi^{2}$ values of the residual profiles from template matching for
each pulsar and TOA dataset. The median goodness-of-fit values are in all
cases close to unity, indicating that template matching fits are generally
good. However, for all three pulsars, the ranges between the 5th and 95th
percentiles are significantly larger for the uncalibrated data, and in the
case of the calibrated data the ranges for MTM TOAs are systematically smaller
than for FDM TOAs. Profile residuals therefore seem to provide a metric for
assessing the correctness of the polarization calibration in this case. We
note finally that in all three MSPs the weighted rms residuals and the reduced
$\chi^{2}$ values from the timing analyses are higher for the “Cal. + corr.
FDM” TOA datasets than for their MTM counterparts. This analysis confirms our
initial assumption that the “Cal. + corr. MTM” datasets, which are based on
better calibrated observations and include TOAs extracted with the MTM
algorithm to account for potential changes of the instrumental response with
time, are our best datasets for conducting the timing analysis of PSRs
J1730$-$2304, J1744$-$1134, and J1857+0943 with NUPPI. Additionally, the
analysis independently confirms the conclusion from Foster et al. (2015),
based on simulations, that the MTM method produces better timing than methods
that use total intensity profiles. van Straten (2006) noted that the MTM
method may not perform well when the analyzed pulsar has a low degree of
polarization or the estimated phase shift is highly covariant with the MTM
model parameters that describe the unknown polarimetric transformation. In
these situations where the FDM method or another method using the total
intensity pulse profiles would be preferred over the MTM algorithm, our
analysis has shown that better calibrated data lead to better timing datasets,
as expected.

Table 2: Weighted rms residual values and reduced $\chi^{2}$ values from individual timing analyses of the six TOA datasets described in Section 3.2 with TEMPO2. The TOA datasets for PSRs J1730$-$2304, J1744$-$1134, and J1857+0943, respectively, comprised 2379, 2085, and 1612 TOAs. TOA dataset | Uncal. FDM | Cal. FDM | Cal. + corr. FDM | Uncal. MTM | Cal. MTM | Cal. + corr. MTM
---|---|---|---|---|---|---
J1730$-$2304
Weighted rms residual ($\mu$s) | 2.213 | 1.710 | 1.604 | 1.090 | 0.974 | 0.960
Reduced $\chi^{2}$ | 4.091 | 2.426 | 2.259 | 2.149 | 1.839 | 1.843
J1744$-$1134
Weighted rms residual ($\mu$s) | 1.612 | 1.042 | 0.705 | 0.905 | 0.635 | 0.628
Reduced $\chi^{2}$ | 27.951 | 10.034 | 5.648 | 4.480 | 1.895 | 2.041
J1857+0943
Weighted rms residual ($\mu$s) | 1.364 | 1.359 | 1.280 | 1.272 | 1.166 | 1.120
Reduced $\chi^{2}$ | 2.001 | 1.958 | 1.903 | 2.203 | 1.898 | 1.817
Figure 16: Distributions of the reduced $\chi^{2}$ values of the residual
profiles from template matching for each TOA, for the three pulsars and six
dataset types considered in Section 3.3. Orange lines represent the median
goodness-of-fit values, boxes encompass the first to third quartiles, whiskers
encompass the 5th and 95th percentiles, and gray dots represent outliers.

##  4 Summary and prospects

Until recently, pulsar data recorded with the NUPPI backend of the NRT at 1.4
GHz were polarization-calibrated using the observations of a noise diode that
injects a polarized reference signal into the feeds, assuming that both the
instrumental response and the noise diode signal are ideal. In late 2019, we
began observing the bright, highly linearly polarized pulsar J0742$-$2822 at a
regular cadence, in a special observing mode in which the feed horn, which is
fixed during normal pulsar observations, rotates by $\sim 180^{\circ}$ across
the typically one-hour observation. Doing so enables us to more fully
determine the polarimetric response of the telescope at the epochs of these
observations as well as the actual Stokes parameters of the noise diode
signal. We found significant cross-couplings between the polarization feeds
and found that the noise diode does not emit 100% linearly polarized emission
along the Stokes $U$ parameter, as was assumed prior to this work. We applied
these improved calibration parameters to the observations of a selection of
MSPs with published polarimetric results, and found our new results to be
consistent with those found by Dai et al. (2015) using the Parkes radio
telescope. Finally, we analyzed NUPPI timing data on the MSPs J1730$-$2304,
J1744$-$1134, and J1857+0943, and found that TOAs extracted from total
intensity pulse profiles obtained with the improved calibration procedure have
lower residual rms and reduced $\chi^{2}$ values than TOAs extracted from
uncalibrated data or from data calibrated with the generic method used prior
to this work. We also found that the MTM method for extracting TOAs from
polarimetric profiles while compensating for potential calibration errors
gives excellent results for these three MSPs, with TOA residuals that have
lower rms and $\chi^{2}$ values than those obtained with TOAs extracted from
total intensity profiles. It should be noted that, even though calibration
errors are expected to affect the pulse profiles of pulsars differently
depending on their polarization properties, they introduce time-correlated
signals in the timing data of different pulsars, and thus contribute noise to
searches for gravitational waves by PTAs. Correcting polarization calibration
errors is thus critical to these searches.

The polarimetric response of the NRT at the epochs of the rotating horn
observations and the properties of the reference source were determined by
solving Equation 11. As explained in detail in van Straten (2004), solutions
to this equation are degenerate under commutation, and assumptions need to be
made to constrain the mixing between $I$ and $V$ and between $Q$ and $U$. We
assumed that the two feeds have equal ellipticities and that the
misorientation of the first receptor is zero. A possibility for determining
the individual receptor ellipticities would be to combine the rotating horn
observations with observations of standard flux candles. In future work we
will carry out such observations of standard candles, close in time to
rotating horn observations, to determine the individual ellipticities across
the bandwidth covered by NUPPI. Additionally, some of the properties of the
NRT, such as the gain of the telescope, are known to vary with the declination
of the observed source (see e.g., Theureau et al., 2005). In this work we have
calibrated observations of MSPs at different sky locations using instrumental
responses determined from observations of J0742$-$2822, which may or may not
be optimal, given the potential dependency of some of the calibration
parameters on declination or other pulsar parameters. This, and the possible
variations of the instrumental response during observations (as noted in
Section 2.2), will also be investigated in future work. Pulsar observations
conducted in a mode in which appropriate feed horn rotation would maintain a
fixed parallactic angle across the observations may be informative in this
regard.

As mentioned in Section 2, the first pulsar observations with the NUPPI
backend were conducted in August 2011, and regular observations of
J0742$-$2822 with feed horn rotation began in late 2019. The improved
calibration solutions derived from these later data cannot be directly applied
to earlier observations because the response of the NRT likely evolved with
the various instrumental changes that occurred since 2011. However, using the
well-calibrated observations of pulsars observed since late 2019 as polarized
reference sources, it is possible to model the instrumental response over the
first eight years of NUPPI observations using Measurement Equation Template
Matching (METM; van Straten, 2013). The improved calibration solutions for the
NUPPI observation archive can be expected to lead to improved timing results,
which is of great interest for pulsar timing studies and PTA analyses.

Observations of PSR J0742$-$2822 with feed horn rotation are conducted at
Nançay to overcome the fact that the NRT cannot sample wide ranges of
parallactic angles when observing a given pulsar. The work presented in this
article is thus relevant for any telescope affected by a similar limitation.
The results obtained with the improved calibration procedure demonstrate that
accurate polarization calibration is essential for any pulsar polarimetric or
timing study, and the need for accurate calibration will get more and more
acute for the increasingly sensitive radio telescopes of the future.

###### Acknowledgements.

We would like to thank the anonymous referee for very constructive and
valuable comments that helped to improve the manuscript. 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. 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/. Part of this research has
also made use of the EPN Database of Pulsar Profiles maintained by the
University of Manchester, available at:
http://www.jodrellbank.manchester.ac.uk/research/pulsar/Resources/epn/. The
published data displayed in Figures 3 to 6 were obtained from the EPN database
under the Creative Commons Attribution 4.0 International licence.

## References

  * Bailes et al. (1994) Bailes, M., Harrison, P. A., Lorimer, D. R., et al. 1994, ApJ, 425, L41
  * Bailes et al. (1997) Bailes, M., Johnston, S., Bell, J. F., et al. 1997, ApJ, 481, 386
  * Bonsignori-Facondi et al. (1973) Bonsignori-Facondi, S. R., Salter, C. J., & Sutton, J. M. 1973, A&A, 27, 67
  * Britton (2000) Britton, M. C. 2000, ApJ, 532, 1240
  * Camilo et al. (1996a) Camilo, F., Nice, D. J., Shrauner, J. A., & Taylor, J. H. 1996a, ApJ, 469, 819
  * Camilo et al. (1996b) Camilo, F., Nice, D. J., Shrauner, J. A., & Taylor, J. H. 1996b, ApJ, 469, 819
  * Chen et al. (2021) Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970
  * Cognard et al. (2017) Cognard, I., Freire, P. C. C., Guillemot, L., et al. 2017, ApJ, 844, 128
  * Cognard et al. (2013) Cognard, I., Theureau, G., Guillemot, L., et al. 2013, in SF2A-2013: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. L. Cambresy, F. Martins, E. Nuss, & A. Palacios, 327–330
  * Dai et al. (2015) Dai, S., Hobbs, G., Manchester, R. N., et al. 2015, MNRAS, 449, 3223
  * Desvignes et al. (2011) Desvignes, G., Barott, W. C., Cognard, I., Lespagnol, P., & Theureau, G. 2011, in American Institute of Physics Conference Series, Vol. 1357, Radio Pulsars: An Astrophysical Key to Unlock the Secrets of the Universe, ed. M. Burgay, N. D’Amico, P. Esposito, A. Pellizzoni, & A. Possenti, 349–350
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341
  * Desvignes et al. (2022) Desvignes, G., Cognard, I., Smith, D. A., et al. 2022, A&A, 667, A79
  * DuPlain et al. (2008) DuPlain, R., Ransom, S., Demorest, P., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7019, Advanced Software and Control for Astronomy II, ed. A. Bridger & N. M. Radziwill, 70191D
  * Dyks et al. (2016) Dyks, J., Serylak, M., Osłowski, S., et al. 2016, A&A, 593, A83
  * Eatough et al. (2013) Eatough, R. P., Falcke, H., Karuppusamy, R., et al. 2013, Nature, 501, 391
  * Foster et al. (2015) Foster, G., Karastergiou, A., Paulin, R., et al. 2015, MNRAS, 453, 1489
  * Foster et al. (1993) Foster, R. S., Wolszczan, A., & Camilo, F. 1993, ApJ, 410, L91
  * Guillemot et al. (2019) Guillemot, L., Octau, F., Cognard, I., et al. 2019, A&A, 629, A92
  * Guo et al. (2021) Guo, Y. J., Freire, P. C. C., Guillemot, L., et al. 2021, A&A, 654, A16
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302
  * Jacoby et al. (2007) Jacoby, B. A., Bailes, M., Ord, S. M., Knight, H. S., & Hotan, A. W. 2007, ApJ, 656, 408
  * Jacoby et al. (2003) Jacoby, B. A., Bailes, M., van Kerkwijk, M. H., et al. 2003, ApJ, 599, L99
  * Jankowski et al. (2018) Jankowski, F., van Straten, W., Keane, E. F., et al. 2018, MNRAS, 473, 4436
  * Jennings et al. (2022) Jennings, R. J., Cordes, J. M., Chatterjee, S., et al. 2022, arXiv e-prints, arXiv:2210.12266
  * Johnston (2002) Johnston, S. 2002, PASA, 19, 277
  * Johnston & Kerr (2018) Johnston, S. & Kerr, M. 2018, MNRAS, 474, 4629
  * Johnston et al. (1993) Johnston, S., Lorimer, D. R., Harrison, P. A., et al. 1993, Nature, 361, 613
  * Kramer et al. (2021) Kramer, M., Stairs, I. H., Manchester, R. N., et al. 2021, Physical Review X, 11, 041050
  * Kraus (1960) Kraus, J. D. 1960, AJ, 65, 54
  * Kraus (1966) Kraus, J. D. 1966, Radio astronomy
  * Lazarus et al. (2016) Lazarus, P., Karuppusamy, R., Graikou, E., et al. 2016, MNRAS, 458, 868
  * Lentati et al. (2014) Lentati, L., Alexander, P., Hobson, M. P., et al. 2014, MNRAS, 437, 3004
  * Liu et al. (2020) Liu, K., Guillemot, L., Istrate, A. G., et al. 2020, MNRAS, 499, 2276
  * Lorimer & Kramer (2004) Lorimer, D. R. & Kramer, M. 2004, Handbook of Pulsar Astronomy, Vol. 4
  * Lorimer et al. (1995) Lorimer, D. R., Nicastro, L., Lyne, A. G., et al. 1995, ApJ, 439, 933
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993
  * Ord et al. (2004) Ord, S. M., van Straten, W., Hotan, A. W., & Bailes, M. 2004, MNRAS, 352, 804
  * Padmanabh et al. (2021) Padmanabh, P. V., Barr, E. D., Champion, D. J., et al. 2021, MNRAS, 500, 1178
  * Park et al. (2021) Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, AJ, 161, 105
  * Philippov & Kramer (2022) Philippov, A. & Kramer, M. 2022, ARA&A, 60, 495
  * Segelstein et al. (1986) Segelstein, D. J., Rawley, L. A., Stinebring, D. R., Fruchter, A. S., & Taylor, J. H. 1986, Nature, 322, 714
  * Singha et al. (2021) Singha, J., Surnis, M. P., Joshi, B. C., et al. 2021, MNRAS, 507, L57
  * Stairs (2003) Stairs, I. H. 2003, Living Reviews in Relativity, 6, 5
  * Stinebring et al. (1984) Stinebring, D. R., Cordes, J. M., Rankin, J. M., Weisberg, J. M., & Boriakoff, V. 1984, ApJS, 55, 247
  * Taylor (1992) Taylor, J. H. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 117
  * Theureau et al. (2005) Theureau, G., Coudreau, N., Hallet, N., et al. 2005, A&A, 430, 373
  * Theureau et al. (2011) Theureau, G., Parent, D., Cognard, I., et al. 2011, A&A, 525, A94
  * van Straten (2002) van Straten, W. 2002, ApJ, 568, 436
  * van Straten (2004) van Straten, W. 2004, ApJS, 152, 129
  * van Straten (2006) van Straten, W. 2006, ApJ, 642, 1004
  * van Straten (2013) van Straten, W. 2013, ApJS, 204, 13
  * van Straten et al. (2010) van Straten, W., Manchester, R. N., Johnston, S., & Reynolds, J. E. 2010, PASA, 27, 104
  * Verbiest & Shaifullah (2018) Verbiest, J. P. W. & Shaifullah, G. M. 2018, Classical and Quantum Gravity, 35, 133001
  * Voisin et al. (2020) Voisin, G., Cognard, I., Freire, P. C. C., et al. 2020, A&A, 638, A24
  * Wang et al. (2022) Wang, J., Shaifullah, G. M., Verbiest, J. P. W., et al. 2022, A&A, 658, A181
  * Wu et al. (2018) Wu, J., Clark, C. J., Pletsch, H. J., et al. 2018, ApJ, 854, 99
  * Xilouris (1991) Xilouris, K. M. 1991, A&A, 248, 323
  * Xu et al. (2021) Xu, H., Huang, Y. X., Burgay, M., et al. 2021, The Astronomer’s Telegram, 14642, 1