Exploring pulsar timing precision: A comparative study of polarization calibration methods for NANOGrav data from the Green Bank Telescope

Abstract

Pulsar timing array experiments have recently uncovered evidence for a nanohertz gravitational wave background by precisely timing an ensemble of millisecond pulsars. The next significant milestones for these experiments include characterizing the detected background with greater precision, identifying its source(s), and detecting continuous gravitational waves from individual supermassive black hole binaries. To achieve these objectives, generating accurate and precise times of arrival of pulses from pulsar observations is crucial. Incorrect polarization calibration of the observed pulsar profiles may introduce errors in the measured times of arrival. Further, previous studies (e.g., van Straten, 2013; Manchester et al., 2013) have demonstrated that robust polarization calibration of pulsar profiles can reduce noise in the pulsar timing data and improve timing solutions. In this paper, we investigate and compare the impact of different polarization calibration methods on pulsar timing precision using three distinct calibration techniques: the Ideal Feed Assumption (IFA), Measurement Equation Modeling (MEM), and Measurement Equation Template Matching (METM). Three NANOGrav pulsars—PSRs J1643$-$1224, J1744$-$1134, and J1909$-$3744—observed with the 800 MHz and 1.5 GHz receivers at the Green Bank Telescope (GBT) are utilized for our analysis. Our findings reveal that all three calibration methods enhance timing precision compared to scenarios where no polarization calibration is performed. Additionally, among the three calibration methods, the IFA approach generally provides the best results for timing analysis of pulsars observed with the GBT receiver system. We attribute the comparatively poorer performance of the MEM and METM methods to potential instabilities in the reference noise diode coupled to the receiver and temporal variations in the profile of the reference pulsar, respectively.

1 Introduction

Pulsars are highly magnetized rapidly rotating neutron stars that emit beams of electromagnetic radiation along their magnetic axes. As the pulsar rotates, the beam(s) sweep(s) across the observer’s line of sight and pulses are seen at regular intervals. The remarkable rotational stability of pulsars, especially millisecond pulsars (MSPs), makes them invaluable tools in astrophysics. MSPs can serve as very accurate celestial clocks and are used for testing theories of gravity, probing the interstellar medium (ISM), and, perhaps most notably, detecting low-frequency gravitational waves (GWs) (Lorimer & Kramer, 2004). At the heart of these scientific endeavors employing MSPs lies the crucial concept of pulsar timing.

In pulsar timing, the rotation of a pulsar is accurately tracked by measuring the times of arrival (TOAs) of its pulses and comparing them to the TOAs predicted from a pulsar timing model. Deviations of the observed TOAs from the predicted ones, i.e., the timing residuals, can reveal important information about the pulsar, its environment, and GW signals present between the Earth and the pulsar. GWs, tidal ripples in the fabric of spacetime, induce minute changes in the TOAs from the pulsars, and pulsar timing array (PTA) experiments observe an ensemble of MSPs to measure those changes in order to detect GWs. Recently, evidence for the presence of a low-frequency stochastic GW background (GWB) in their respective PTA data sets was reported independently by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav: Agazie et al., 2023), the European Pulsar Timing Array (EPTA) + Indian Pulsar Timing Array (InPTA: Antoniadis et al., 2023), the Parkes Pulsar Timing Array (PPTA: Reardon et al., 2023), and the Chinese Pulsar Timing Array (CPTA: Xu et al., 2023) collaborations.

The next significant milestones for the PTAs to accomplish are characterizing the observed GWB more precisely, determining its source(s), and detecting continuous GWs from individual supermassive black hole binaries (SMBHBs). The PTA sensitivity directly depends upon the precision and accuracy with which pulse arrival times can be estimated and therefore generating accurate and precise TOAs from pulsar observations is a critical aspect of the PTA experiments. Observing pulsars using bigger and better telescopes with larger bandwidths and longer integration times is one of the possible ways to improve the precision of TOAs. Additionally, better precision and accuracy in TOAs can be achieved by developing improved methods of TOA estimation, radio frequency interference (RFI) mitigation, and instrumental calibration. It is noteworthy that inaccurate polarization calibration of the observed pulsar profiles may introduce errors in measured TOAs (see e.g., Rogers, 2020). Hence, adopting a robust and accurate polarization calibration procedure to calibrate pulsar profiles holds the potential to elevate the precision and accuracy of measured TOAs, thereby facilitating the realization of the current scientific goals set by PTA collaborations.

Pulsars are among the most polarized of all known radio sources and polarization measurement can provide additional insights into the pulsar emission process and the medium through which the radiation propagates. The polarization state of a pulsar signal can be described by the four Stokes parameters $I,Q,U$, and $V$, and can be represented by the Stokes vector (Stokes, 1851)

$$S=\begin{bmatrix}I\\ Q\\ U\\ V\\ \end{bmatrix}\,,$$ (1)

where $I$ is the total intensity, $Q$ and $U$ form the linear polarization $L=\sqrt{Q^{2}+U^{2}}$, and $V$ represents the circular polarization intensity. Using the International Astronomical Union’s (IAU’s) convention, right-handed circular polarization is positive and left-handed circular polarization is negative (Stokes, 1851). Further, the position angle of the linear polarization that represents the orientation of the linearly polarized radio waves emanating from the pulsar is given by

$$\Psi=0.5\,\tan^{-1}\left({\frac{U}{Q}}\right)\,.$$ (2)

The polarization angle is quoted using the IAU convention with the polarization angle increasing in the counterclockwise direction. As the radio waves from a pulsar travel through the ISM, they experience Faraday rotation, a frequency-dependent rotation of the polarization position angle caused by the Galactic magnetic field. Interstellar Faraday rotation ($\beta$) can be given by

$$\beta=RM\,\,\lambda^{2}\,,$$ (3)

where $\lambda$ is the wavelength of the radio waves and the overall strength of the effect is characterized by the rotation measure ($RM$). The $RM$ depends on the interstellar magnetic field component ($B_{||}$) parallel to the line of sight and free-electron density ($n_{e}$) as (in cgs units)

$$RM=\frac{e^{3}}{2\pi m_{e}^{2}c^{4}}\int_{0}^{d}n_{e}(l)B_{||}(l)\,dl\,,$$ (4)

where $e$ and $m_{e}$ are the charge and mass of electron, respectively, $c$ is the speed of light in vacuum, and $d$ is the distance to the pulsar. Measurements of Faraday rotations in linearly polarized pulsars are used to study the ISM and the large-scale Galactic magnetic field in the Milky Way (Han et al., 2018; Sobey et al., 2019).

When a pulsar is observed with a radio telescope, the processes of reception and detection of the signal introduce instrumental artifacts that make the measured Stokes vector $S_{\text{m}}$ differ from the intrinsic Stokes vector $S_{\text{i}}$. We can relate $S_{\text{m}}$ to $S_{\text{i}}$ with the Mueller matrix $M$ such that

$$S_{\text{m}}=M\,S_{\text{i}}\,,$$ (5)

where $M$ depends on differential gain and differential phase of the receiver, and ellipticity and non-orthogonality of the feeds (Heiles et al., 2001). During the polarization calibration, we determine $M$ by calibrating the observing system and thereafter solve Equation (5) to obtain the true Stokes vector $S_{\text{i}}$ from the observed $S_{\text{m}}$. However, there are different methods for performing the pulsar polarization calibration proposed in the literature.

In the NANOGrav data releases, the pulsar profiles are calibrated using Ideal Feed Assumption (IFA) which assumes the feeds to be perfectly linear and orthogonal. However, in reality, this assumption may not be valid, and therefore using IFA-calibrated pulsar profiles to generate TOAs can introduce systematic errors due to polarization miscalibration. Better polarization calibration methods, where a full polarimetric response (PR) of the observing system is used to calibrate the pulsar profiles, have been proposed to overcome the shortcomings of the IFA approach. van Straten (2004) developed Measurement Equation Modeling (MEM) that uses a pulsar with strong linear polarization observed over a wide range of parallactic angles to estimate the full PR of the observing system. Additionally, Measurement Equation Template Matching (METM) was introduced by van Straten (2013) as a polarization calibration method that matches multiple observations of a reference pulsar to a well-calibrated template profile of that pulsar to generate precise PR solutions at different epochs. In both the MEM and METM methods, the calculated receiver solutions are used to perform polarization calibration of the observed pulsar profiles, and then accurate and precise TOAs can be generated from the calibrated profiles.

The effects of these robust polarization calibration techniques on the accuracy of the timing analysis for different pulsars have been explored in many studies. Manchester et al. (2013) have found that, for 9 of the 20 pulsars observed with the Parkes radio telescope at 20-cm band, MEM-calibrated profiles provided timing residuals with reduced root-mean-square (RMS) and chi-square values compared to uncalibrated profiles. However, for rest of the pulsars, the MEM calibration made little difference to the reduced chi-square of the timing solution. van Straten (2013) compared the timing accuracy of PSR J1022+1001 for different combinations of polarization calibration and TOA generation methods using the Parkes telescope data. Two different methods for TOA generation were used: matrix template matching (MTM) which uses the full polarization profiles for generating TOAs (van Straten, 2006), and standard total intensity (STI) which uses only the total intensity profile. They found that calibration with the METM method reduces both the standard deviation of the arrival time residuals and the reduced chi-square of the model fit compared to IFA calibration, for both the STI and MTM methods of TOA generation. A detailed analysis of the effects of different polarization calibration methods on pulsar TOAs was also conducted by Rogers (2020) for five pulsars observed with the Parkes radio telescope. The study found that METM polarization calibration combined with MTM for TOA generation gives better timing accuracy compared to traditional IFA calibration followed by STI for TOA generation for all five pulsars.

Furthermore, Gentile et al. (2018) performed METM polarization calibration on a subset of NANOGrav data observed with the Arecibo Telescope to obtain some of the most sensitive polarimetric MSP profiles. This was repeated for a subset of Green Bank Telescope (GBT) pulsar data in Wahl et al. (2022) where they used MEM calibration method for the polarization calibration. However, a detailed timing analysis with those calibrated profiles have not been done yet (see Wahl (2022) for a precursor work). These results motivated us to explore the effects of different polarization calibration methods on the timing analysis for a subset of NANOGrav pulsars observed at the GBT.

In this paper, we present the timing analysis of three pulsars: PSRs J1643$-$1224, J1744$-$1134, and J1909$-$3744, with different polarization calibration methods. We have used three different methods for performing polarization calibration of the pulsar profiles: (i) IFA, (ii) MEM, and (iii) MEM+METM, and TOAs were generated from the calibrated profiles using the STI. The sets of TOAs generated for different polarization calibration methods are then individually used to perform timing analysis to compare the influence of different polarization calibrations on timing analysis for these pulsars.

The paper is structured as follows. In Section 2, we briefly describe the data we used in our paper. Different polarization calibration methods are discussed in Section 3 along with a brief overview of the timing analysis procedure used in this paper. Section 4 contains the results of our analyses and the results are summarized and discussed in detail in Section 5.

2 Data

In this paper, we focus on three pulsars observed with the 100-meter Green Bank Telescope in the NANOGrav program: PSRs J1643$-$1224, J1744$-$1134, and J1909$-$3744. We only use a subset of the data taken with the Green Bank Ultimate Pulsar Processing Instrument (GUPPI; DuPlain et al., 2008) backend system at both 820 MHz (with Rcvr_800) and 1500 MHz (with Rcvr1_2) frequencies with bandwidths of 200 MHz and 800 MHz, respectively. PSR J1643$-$1224 is a bright pulsar with average flux densities of 12.9 and 4.7 mJy at 820 MHz and 1500 MHz frequencies, respectively, and has a moderate polarization fraction ($\sim 21\%$ at both 820 and 1500 MHz frequencies). This pulsar was examined in Rogers (2020), which allows a direct comparison of our results. PSR J1744$-$1134 has average flux densities of 6.2/2.6 mJy at 820/1500 MHz frequencies and is highly polarized (polarization fractions at 820/1500 MHz are $\sim 78\%/88\%$). PSR J1909$-$3744, one of the best pulsars in NANOGrav, also has a fairly high polarization fraction (the polarization fractions at 820/1500 MHz are $\sim 51\%/45\%$ and the average flux densities are 3.6/1.3 mJy).

Although the GUPPI data acquisition instrument was used at GBT from 2010 March to 2020 April for NANOGrav observations, we only use the data from 2010 March (MJD 55265) to 2014 March (MJD 56739) in this paper. This is because of a technical problem, which arose in 2014 March, that made the time alignment of the digitizers for the X and Y polarization of the telescope signals unstable, thus corrupting the polarization cross products (Wahl et al., 2022). Therefore, it is impossible to recover the correct full Stokes parameters from the data taken after March 2014. However, this instability only affects the polarization information of the observed pulsar profiles and should not affect the total intensity profiles and therefore timing after 2014 March.

The pulsars were observed with an approximately monthly cadence and the data were coherently de-dispersed, with frequency resolution of $\sim 1.56$ MHz. The resulting time series were folded in real time using a nominal pulsar timing model to obtain folded pulse profiles as functions of time, radio frequency, and polarization. The folded profiles have 2048 phase bins and subintegrations of 10 seconds. We perform polarization calibrations on these folded profiles to get accurately calibrated pulsar profiles and thereafter use them to generate TOAs for timing analysis. In addition, we used three long-track observations of PSRs B1929+10 (on MJDs 56244, 56419, and 56608; Kramer et al., 2021) at 820 MHz and one long-track observation of J1022+1001 (MJD 55671) at 1500 MHz with the GUPPI backend system to generate full PRs of the observing system for the MEM calibration method. Additionally, observations of PSR B1937+21 at GBT during MJD $55265-56739$ are used as reference pulsar profiles for performing the METM polarization calibration.

3 Methods

In this section, we discuss various methods employed in this paper. Different polarization calibration methods, namely the IFA, MEM, and METM, are described in detail in Subsection 3.1. Subsection 3.2 outlines our methods for TOA generation and timing analysis.

3.1 Polarization Calibration Methods

As discussed in Section 1, the measured polarization profiles or the Stokes vector ($S_{\text{m}}$) of a pulsar observed with a radio telescope differs from the intrinsic one ($S_{\text{i}}$) due to instrumental effects. The measured and intrinsic Stokes vectors are related by the Mueller matrix $M$ (see 5) that characterizes the PR of the telescope. The Mueller matrix for a dual linear feeds can be written in the form (following Lorimer & Kramer, 2004; Heiles et al., 2001; Gentile et al., 2018)

$$M=\begin{bmatrix}1&E&A+EC&B+ED\\ E&1&AE+C&BE+D\\ AF-BG&DG-CF&F&-G\\ AG+BF&-CG-DF&G&F\\ \end{bmatrix}\,,$$ (6)

where

$$A = \epsilon_{1}\cos{\theta_{1}}+\epsilon_{2}\cos{\theta_{2}}\,,$$

$$B = \epsilon_{1}\sin{\theta_{1}}+\epsilon_{2}\sin{\theta_{2}}\,,$$

$$C = \epsilon_{1}\cos{\theta_{1}}-\epsilon_{2}\cos{\theta_{2}}\,,$$

$$D = \epsilon_{1}\sin{\theta_{1}}-\epsilon_{2}\sin{\theta_{2}}\,,$$

$$E = \gamma/2\,,$$

$$F = \cos{\phi}\,,$$

$$G = \sin{\phi}\,,$$

and and $\epsilon_{1}$ and $\epsilon_{2}$ represent the magnitude of the cross-coupling of the two respective feeds, $\theta_{1}$ and $\theta_{2}$ represent the phase of this cross-coupling, $\gamma$ represents the differential gain of the receiver, and $\phi$ represents the differential phase of the receiver. By calibrating the observing system, we can determine $M$ and solve Equation (5) for the true Stokes vector $S_{\text{i}}$. In different polarization calibration methods, the instrumental PR is calculated differently and therefore the calibrated profiles (i.e., the $S_{\text{i}}$s generated by solving Equation (5)) could be different. We now describe the different polarization calibration methods used in this paper.

3.1.1 Ideal Feed Assumption

Ideal Feed Assumption (IFA) is the most basic polarization calibration method and has been employed for all NANOGrav data releases. As the signals from the two feeds during the observation pass through slightly different amplifier chains, it introduces different gains and phases on the signal. In this method, observation of a reference source is used to determine the differential gain and phase of the receiver system. For NANOGrav observations, the reference source is an artificial noise diode coupled to the receptors, emitting a square wave with a 50% duty cycle and a period of 40 ms. The noise diode is expected to be 100% linearly polarized and should illuminate both receptors equally and in phase. This noise diode is observed before each pulsar observation and the observations are used to determine the complex gains of the instrumental response, described by the absolute gain $G$, differential gain $\gamma$, and differential phase $\phi$ of the receiver system, as a function of observing frequency. The noise diode signal amplitude is further calibrated by observing it on and off a bright and unpolarized standard continuum calibrator radio source (flux calibrator). Using flux calibrator observations during the calibration eliminates the assumption that the reference source illuminates both receptors equally. For NANOGrav GBT observation with the GUPPI system, quasar B1442+101 was used as the continuum calibrator and observed approximately once per month. One example IFA polarimetric response for calibrating the 800 MHz data is shown in panel (a) of Figure 1. These instrumental responses for each pulsar observations are then used to calibrate the pulsar profiles.

(a) IFA

(b) MEM	(c) METM

Figure 1: Example polarimetric response solutions used for different polarization calibration methods to calibrate the observed pulsar profiles. In all the panels $G$, $\gamma$, and $\phi$ represent the absolute gain, differential gain, and differential phase of the observing system, respectively. Panel (a): IFA PR solution obtained from the reference noise diode observation for PSR J1643$-$1124 on MJD 56614. Panel (b): MEM PR solution calculated from a long track observation of PSR B1929$+$10 on MJD 56608. Here, $\theta_{1}$ represents the orientation of receptor 1 with respect to receptor 0, and $\epsilon_{k}$ is the ellipticities of the two receptors (denoted by black and red points). Panel (c): METM-generated PR correction, calculated using B1937$+$21 observation on MJD 56614, to the MEM-generated PR on MJD 56608. Here, $\theta_{k}$ and $\epsilon_{k}$ represent the orientations and ellipticities of the two receptors (black and red points), respectively. Note that the absolute gain is not calculated for the METM-generated PR correction as we used the total invariant interval to normalize the Stokes parameters. See Section 3.1 for more details.

It is important to note that this model relies on the assumption that the receptors are perfectly orthogonal and the noise diode signal is 100% linearly polarized. We have used the pac command in PSRCHIVE (Hotan et al., 2004) to perform IFA calibrations of the observed pulsar profiles.

3.1.2 Measurement Equation Modeling

Developed in van Straten (2004), Measurement Equation Modeling (MEM) utilizes a long-track observation (or multiple observations) of a pulsar with strong linear polarization over a wide range of parallactic angles to generate a complete PR of the observing system at a fiducial epoch. This method is based on the polarization measurement equation (Hamaker, 2000), which relates the measured Stokes parameters to the intrinsic pulsar polarization and is employed to determine the unknown instrumental response.

At first, a preliminary IFA polarization calibration is performed on the long-track observation to create a preliminary total integrated calibrated profile. Thereafter, a specified set of pulse phase bins from this profile is chosen to be used as model constraints. The Stokes parameters at these phase bins are fitted as a function of parallactic angles using the polarization measurement equation, and a best-fit solution for the complete instrumental PR at the fiducial epoch is calculated. The complete PR is parameterized by the absolute gain $G$, differential gain $\gamma$, differential phase $\phi$, ellipticities of the two receptors $\epsilon_{k}$, and orientation of the receptor-1 with respect to receptor-0 $\theta_{1}$ as functions of observing frequency, when we use the van04e18 parametrization (equation 18 of van Straten, 2004); see panel (b) of Figure 1 for an example PR. During this process, a flux calibrator observation is used to constrain the mixing of Stokes I and V, and this breaks the degeneracy described in van Straten (2004).

The pcm command in PSRCHIVE is used to generate the MEM PRs of the observing system from pulsar observations over a wide range of parallactic angles. For calibrating Rcvr_800 data, we used three long-track observations of PSR B1929+10 at 820 MHz on MJDs 56244, 56419, and 56608 to generate three separate MEM PR solutions at those three epochs. This bright pulsar has well-known polarization characteristics and the data used here were acquired by Kramer et al. (2021) for the double pulsar. A single long-track observation of PSR J1022+1001 on MJD 55671 is used to derive the MEM solution at 1500 MHz for calibrating the Rcvr1_2 data. Although PSR J1022+1001 has been found to exhibit long-term pulsar profile variation, we do not expect that to affect our result as the MEM solution is calculated from an observation of the pulsar on a single day (Wahl et al., 2022).

To perform the polarization calibration of a pulsar observation using the MEM method, the MEM solution from the fiducial epoch closest to the pulsar observation is first used to correct for the complete system response at that fiducial epoch. Thereafter, the noise diode observation taken prior to the pulsar observation is employed to correct for any variations in differential gain or phase since the reference epoch. The MEM polarization calibration method is based on the assumption that the receptor orientations and ellipticities, as well as the polarization of the reference source, do not vary significantly with time.

3.1.3 Measurement Equation Template Matching

Measurement Equation Template Matching (METM) was developed by van Straten (2013), combining MTM with the MEM. In this method, the observation of a reference pulsar is matched with a well-calibrated template profile of that pulsar to obtain the best-fit METM model for the transformation between the template profile and the observation. The METM model can be used to fully represent the PR of the observing system at the observation epoch. This process can be repeated for multiple observations of the reference pulsar to obtain a complete PR of the observing system for each observation epoch. We can then use the METM models to calibrate other pulsar observations at those or nearby epochs. This method is valid under the assumption that the polarized emission from the reference pulsar remains constant over time.

In this paper, we use PSR B1937$+$21 as the reference pulsar for METM calibration due to its high brightness and well-known polarization characteristics. Additionally, B1937$+$21 is observed in the same session as the three pulsars used in our analysis by NANOGrav, thus providing METM solutions for the exact epochs we require. To generate METM solutions from the B1937+21 observations and perform polarization calibration on other pulsar observations, we followed the method outlined in Gentile et al. (2018). Initially, all B1937$+$21 observations are calibrated using the MEM method, following the procedure described in Section 3.1.2. Among the calibrated profiles, two profiles are selected as template profiles—one for the 820 MHz band and another for the 1500 MHz band. The template profiles are chosen to match the previously published polarization profiles of B1937$+$21 at the respective observation frequencies (Dai et al., 2015). Thereafter, METM PRs are calculated for each observation epoch by comparing the MEM-calibrated profiles to the template profiles. The pcm command from PSRCHIVE is used to obtain the METM PRs, and we used the option to normalize the Stokes parameters by the total invariant interval (Britton, 2000) instead of calculating the absolute gain. This decision was made because the reference pulsar B1937$+$21 exhibits pronounced scintillation (Turner et al., 2024), which could significantly bias the absolute gain due to variations in pulse intensity caused by scintillation across different frequencies.

Given that the observations used to obtain the METM PRs have already undergone calibration with the MEM-generated PR, it is appropriate to consider the METM-generated PRs as per-epoch corrections to the MEM-generated PR. Accordingly, we apply these corrections to other pulsar observations that are already calibrated using MEM PRs. This is done by selecting the METM PR correction whose epoch is closest to the epoch of the observation that is to be calibrated. In our paper, this method of polarization calibration is denoted as MEM + METM, signifying the combination of both MEM and METM methods. One example MEM+METM PR correction calculated from 800 MHz B1937+21 observation on MJD 56614 is shown in panel (c) of Figure 1. In the following section, we discuss the process of generating TOAs from the calibrated profiles and performing timing analysis.

3.2 TOA Generation and Timing Analysis

We started with the raw folded pulsar profiles and removed artifacts arising from the interleaved analog-to-digital converter scheme used by the GUPPI receiver system, as described in Section 2.3.1 of Arzoumanian et al. (2020). Thereafter, we performed the standard radio frequency interference (RFI) excision and calibrated the profiles using the three methods described above. After calibrating the pulsar profiles, we performed additional steps of excising RFI from the calibrated profiles, following Agazie et al. (2023). For MEM and MEM+METM calibrated profiles, we utilized rmfit to determine the optimal fit for the RM value associated with each profile and subsequently apply this calculated RM value to the profile header. The clean calibrated profiles are then frequency-averaged into 64 channels (for both Rcvr_800 and Rcvr1_2 data) and time-averaged into subintegrations of up to 30 minutes. Thereafter, narrowband template profiles are generated for each receiver band using the following steps. Profiles associated with a specific pulsar and receiver band are aligned, weighted based on signal-to-noise ratio, and then summed to create a final averaged profile. The resulting average profile is then ‘denoised’ using wavelet decomposition and thresholding, and the whole process is iterated multiple times to converge on the final template (for more details, see Section 3.1 of Demorest et al., 2013). Finally, TOAs are obtained from folded pulse profile data by measuring the time shift of the observed profile relative to the template, following the standard approach used in pulsar timing analyses for decades (e.g., Taylor, 1992). We have used only the total intensity profiles of the pulsars while generating the TOAs. The calculated TOAs are narrowband in nature, i.e, a separate TOA is measured for each frequency channel of the final profiles. We used the nanopipe (Demorest, 2018) data processing pipeline, which in turn uses the PSRCHIVE pulsar data analysis software package, to perform these steps. For more details on the TOA generation procedure, see Section 3 of Agazie et al. (2023) and references therein. It is crucial to emphasize that distinct templates have been generated for each polarization calibration method, and the corresponding template is employed when calculating the TOAs from profiles calibrated using different methods.

After obtaining the TOAs, we performed outlier removal on the files containing TOAs (tim files) for each polarization calibration method, following the procedures outlined in Section 3.3 of Agazie et al. (2023). At first, an initial timing solution (timing model parameter file or par file) was derived using the tim file, which contains both Rcvr_800 and Rcvr1_2 TOAs. Subsequently, the initial par and tim files were passed to an automated outlier analysis pipeline, which removed TOAs with outlier probabilities exceeding $0.1$ from the tim file. Furthermore, TOAs from a specific folded profile were excluded if a substantial percentage of them were flagged as outliers. In addition, any TOA associated with profile data that did not meet the signal-to-noise ratio threshold (S/N $>8$) was removed from the tim fil (see (Agazie et al., 2023) for more details). The resulting excised tim file was then used for the subsequent timing analysis.

We employed the PINT pulsar timing software (Luo et al., 2021) and adopted the procedure outlined in the NANOGrav 15-year data release (Agazie et al., 2023) for our timing analysis. We started with the NANOGrav 15-year .par files and the human-readable configuration (.yaml) files, and used the standardized Jupyter notebooks to automate our timing procedure. The timing analysis was conducted using the JPL DE440 solar system ephemeris (Park et al., 2021) and the TT(BIPM2021) timescale. Pulsar timing involves comparing observed TOAs with TOAs predicted by a timing model, yielding timing residuals. This model encompasses various timing parameters representing different physical effects influencing pulse arrival times. During the pulsar timing, best fit values of the timing model parameters are calculated to minimize the RMS of the timing residual. Here we briefly describe the different timing parameters used in our analysis.

For each pulsar, we fit for two spin parameters (rotational frequency and frequency derivative) and five astrometry parameters (two-dimensional sky position and proper motion, and parallax). For PSRs J1643$-$1224 and J1909$-$3744, we also fit binary parameters describing orbit with a companion star. For PSR J1643$-$1224, we employ the DD binary model (Damour & Deruelle, 1985), incorporating the following six parameters: orbital period $P_{b}$, projected semi-major axis $x$ and its time derivative $\dot{x}$, orbital eccentricity $e$, longitude of periastron $\omega$, and epoch of periastron passage $T_{0}$. The ELL1 binary model (Lange et al., 2001) is used for PSR J1909$-$3744, where in addition to $P_{b}$, $x$, and $\dot{x}$, we have incorporated $\dot{P}_{b}$, the companion mass $m_{2}$, orbital inclination parameter $\sin{i}$, two Laplace–Lagrange parameters ($\epsilon_{1}$, $\epsilon_{2}$) and the epoch of the ascending node $T_{asc}$.

The variation in the interstellar dispersion measure (DM) is mitigated using the DMX model, a piecewise constant function, with each DMX parameter describing the offset from a nominal fixed value. We also fit for additional time-independent but frequency-dependent delays on a per-pulsar basis using “FD” parameters (see NANOGrav Collaboration et al. (2015)). The FD parameters are thought to account for time offsets resulting from disparities between the observed pulse shape at a specific frequency and the template shape used in timing, and the number of these parameters included is determined via the F-test procedure discussed in Section 4.1.2 of Agazie et al. (2023). Furthermore, we incorporate “JUMP”s to address unknown phase offsets between data observed by different receivers.

The scripts used for performing the polarization calibration and TOA generation, along with all the MEM and METM PR solutions can be found in https://github.com/lanky441/psrcal_scripts. The PINT-based timing and outlier analysis packages, along with relevant example Jupyter notebooks used in our analysis, are available in PINT_pal.

4 Results

Figure 2: Uncalibrated and calibrated (using different methods) polarization profiles of PSR J1909$-$3744 in the Rcvr1_2 band (1500 MHz). The calibration method used to obtain the profiles are denoted in the plot titles. In each panel, the black, red, and blue lines indicate the total intensity (I), linear polarization (L), and circular polarization (V), respectively. Using the IAU’s circular polarization sign convention, right-handed circular polarization is positive and left-handed circular polarization is negative.

In this section, we present the outcomes of our diverse polarization calibration processes and the corresponding timing analysis for PSRs J1643$-$1224, J1744$-$1134, and J1909$-$3744. To illustrate the varied results from different polarization calibration procedures, Figure 2 displays all the Rcvr1_2 profiles of PSR J1909$-$3744 from different observation epochs. Each panel in the figure represents profiles obtained through distinct polarization calibration methods, with the black, red, and blue lines indicating the total intensity (I), linear polarization (L), and circular polarization (V), respectively. Additionally, we include the original uncalibrated profiles (after RFI excision) in the top-left panel of the figure for comparison. We observe from Figure 2 that, in the case of uncalibrated profiles, although the total intensity profiles closely align across different epochs, significant variability exists in both the linear and circular polarization profiles from epoch to epoch. Polarization calibration proves effective in mitigating variations in linear polarization profiles, with only a very few epochs showing exceptions, across all three calibration methods. The same holds true for circular polarization, although the MEM and MEM+METM methods seem to do a better job compared to the IFA polarization calibration. Similar trends are seen for Rcvr_800 observing frequency and for other pulsars, in general. However, the calibration processes are not always able to reasonably mitigate the epoch-to-epoch variations in the linear and circular polarization profiles. This is particularly prominent for PSR J1744$-$1134 where we see variations (at different levels) in the circular polarization even after performing the polarization calibration. The uncalibrated and calibrated profiles for J1909$-$3744 Rcvr_800 observations and of PSRs J1643$-$1224 and J1744$-$1134 are shown in the Appendix A.

Table 1: Timing analysis statistics for different polarization calibration methods

Pulsar	Method	XN${}_{\text{TOA}}$X	$\sigma_{\text{med}}$	Median	XRMS ($\mu$s)X	XWRMS ($\mu$s)X	Reduced
			($\mu$s)	SNR	All TOAs/Epoch-avg.	All TOAs/Epoch-avg.	chi-square
J1643$-$1224	Uncalibrated	7352	2.02	111.62	4.126 / 0.763	3.464 / 0.780	3.664
	IFA	7172	1.99	114.10	3.033 / 0.822	2.332 / 0.782	1.682
	MEM	7019	2.02	111.85	3.686 / 0.830	2.909 / 0.809	2.520
	MEM+METM	6851	2.13	104.66	3.771 / 0.806	2.918 / 0.793	2.163
J1744$-$1134	Uncalibrated	7467	0.997	76.93	2.823 / 0.749	0.655 / 0.277	2.892
	IFA	7433	0.996	77.41	2.786 / 0.544	0.454 / 0.169	1.475
	MEM	7190	1.07	72.08	2.959 / 0.470	0.572 / 0.179	2.030
	MEM+METM	6726	1.11	68.80	2.962 / 0.575	0.669 / 0.195	2.397
J1909$-$3744	Uncalibrated	9643	0.498	55.06	1.187 / 0.139	0.176 / 0.051	1.464
	IFA	9488	0.502	54.54	1.193 / 0.089	0.151 / 0.034	1.108
	MEM	9163	0.513	53.10	1.208 / 0.101	0.164 / 0.035	1.229
	MEM+METM	8732	0.544	50.06	1.247 / 0.112	0.183 / 0.037	1.330

In Table 1, various quantities and statistics related to the timing analysis for each pulsar are presented. The columns, from left to right, display the pulsar name, polarization calibration method, number of TOAs (N${}_{\text{TOA}}$), median uncertainty of the TOAs ($\sigma_{\text{med}}$), median signal-to-noise ratio (SNR) of each sub-band of the pulsar profiles from which the TOAs are calculated, RMS (all TOAs/epoch averaged) and weighted RMS (all TOAs/epoch averaged) of the timing residuals, and the reduced chi-square of the timing solution, respectively. For each pulsar, the first row shows the results when we do not perform any polarization calibration before generating the TOAs, while the subsequent three rows represent results for IFA, MEM, and MEM+METM calibration.

For all three pulsars, we observe that the number of TOAs utilized in the timing analysis (post-outlier analysis) is highest when no polarization calibration is performed. This is primarily because polarization calibration can occasionally corrupt a few frequency channels, typically due to presence of RFI in the reference noise diode or the reference pulsar observation. These corrupted channels are subsequently either zapped during RFI excision or the corresponding TOAs are excised during the outlier analysis. Additionally, the calibration process can fail for a few channels due to the absence of a solution for those channels in the polarization response calculated using the MEM or MEM+METM methods. Further, the MEM+METM calibration process failed for five J1744$-$1134 Rcvr_800 observations due to a mismatch in the center frequency and observation bandwidth between the observed profiles and METM PRs. We also see from the Table 1 that the calculated median TOA uncertainties for the .tim files generated from the uncalibrated and IFA-calibrated profiles are very similar and have the lowest values. Conversely, the $\sigma_{\text{med}}$ values are highest for the MEM+METM calibration, while the MEM calibration falls in the middle in this regard. This trend aligns with the observed variations in the SNR between the uncalibrated and different calibrated profiles. Typically, the SNR is highest for uncalibrated and IFA-calibrated profiles, lower for MEM-calibrated profiles, and lowest for MEM+METM-calibrated profiles.

We now shift our focus to examining various statistics of the timing solutions obtained for different calibration methods applied to each pulsar. Upon inspecting the reduced chi-square values, it becomes evident that the TOAs generated from all three calibration methods exhibit a superior fit to the timing model compared to the uncalibrated TOAs. Among the different calibration methods, the IFA calibration yields reduced chi-square values closest to unity for all three pulsars. However, there are variations observed while comparing the MEM and MEM+METM calibration methods. For PSR J1643$-$1224, the reduced $\chi^{2}$ value is closer to unity for the MEM+METM calibration method compared to the MEM method, whereas for PSRs J1744$-$1134 and J1909$-$3744, the opposite holds true.

In Table 1, we also present two sets of values for the RMS and Weighted RMS (WRMS) of the timing solutions: one labeled as All TOAs, and the other as Epoch-averaged. As their names suggest, the All TOAs RMS and WRMS values are computed using all the sub-banded TOAs. Conversely, the Epoch-avg. values are derived by averaging the TOAs from all sub-bands for a specific epoch to generate a single TOA. Upon comparing the RMS and WRMS values across different calibration methods, we observe variations in the trend among different pulsars. However, in most cases, the IFA calibration method consistently provides the lowest RMS and WRMS values compared to other methods. Therefore, it is evident from Table 1 that overall, the IFA polarization calibration works best for the GBT data taken with the GUPPI receiver system in terms of timing analysis. The other calibration methods, such as MEM and MEM+METM, offer improvements over performing no calibration, yet they demonstrate inferior performance compared to the IFA calibration method for our dataset.

5 Discussions

When a pulsar is observed with a radio telescope, instrumental artifacts can distort its total intensity profile, resulting in significant systematic timing errors. Therefore, ensuring accurate polarization calibration of the observed pulsar profiles is crucial for achieving high-precision timing, which is fundamental for PTA experiments. In this paper, we compare the performance of three different polarization calibration methods—IFA, MEM, and MEM+METM—using GBT observations of PSRs J1643$-$1224, J1744$-$1134, and J1909$-$3744 with the GUPPI receiver system. Our findings indicate that all three calibration methods improve timing precision compared to scenario where no polarization calibration is conducted. This improvement is expected as polarization calibration corrects for instrumental response, resulting in more stable intrinsic total intensity pulse profiles. Consequently, this leads to more accurate and precise estimation of TOAs, thereby enhancing timing precision.

Based on previous studies (van Straten, 2013; Manchester et al., 2013; Rogers, 2020), we anticipated that the MEM and MEM+METM calibration methods would yield better timing performance compared to the IFA calibration. This expectation also arises from the fact that, unlike MEM and MEM+METM methods, the IFA calibration does not correct for the ellipticities and non-orthogonality of the receiver feeds. Contrary to our expectations, however, it was found that the IFA-calibrated data produced TOAs with the smallest errors. In order to understand these results, it is important to recall that each of the polarization calibration methods operates based on distinct assumptions regarding the receiver and reference noise diode systems. For IFA, the assumptions are that the receptors are perfectly orthogonally polarized and the noise diode is 100% linearly polarized. The MEM calibration method assumes that the orientations and ellipticities of the receptors, as well as the polarization of the reference noise diode, do not significantly vary over time. Lastly, the METM calibration requires the polarization profiles of the reference pulsar to be stable and not subject to variation over time.

Upon examining the MEM PR solutions (one of which is depicted in panel (b) of Figure 1), we observe that the ellipticities and non-orthogonality of the receptors (represented by $\epsilon_{k}$ and $\theta_{1}$) are very small across the observing frequency range. Therefore, the assumption in IFA approach of perfectly orthogonally polarized receptors appears to be fairly accurate for the receiver systems at GBT and is unlikely to significantly degrade the accuracy of the polarization calibration. Determining whether the reference noise diode is inherently 100% linearly polarized presents a challenge as the instrumental response must be decoupled from the noise diode observations. PRs computed in the IFA approach assume the noise diode to be 100% linearly polarized, rendering them unsuitable for correcting the instrumental response in the noise diode observations. However, if MEM PR solutions are available from long-track observation of a pulsar, they can be employed to independently correct for the instrumental response of the associated noise diode observation and model the intrinsic Stokes parameters of the noise diode signal.

In Figure 3, we show the intrinsic Stokes parameters of the noise diode signal in the Rcvr_800 band by presenting the fractional Stokes $Q$, $U$, and $V$, modeled using the MEM PR solutions on MJDs 56244, 56419, and 56608. For an ideal reference source illuminating both receptors equally, $U$ should register at 100%, while $Q$ and $V$ should remain at zero across the entire frequency range. However, since we incorporate a flux calibrator observation during our IFA calibration, the assumption of equal illumination on both receptors is eliminated, and therefore, the noise diode signal only needs to be 100% linearly polarized. Inspection of Figure 3 reveals that across all three epochs, the reference signal is $\sim 95-100$% linearly polarized for the majority of the band. It is also evident from the figure that the intrinsic polarization characteristics of the reference noise diode exhibit variability from epoch to epoch. At the Rcvr1_2 observing frequencies, we observe from the single MEM solution available on MJD 55671 that the reference noise diode signal consists of $\sim 90$% linear polarization. However, due to the limited availability of MEM solutions for only one epoch, we lack information regarding the temporal variations in the polarization of the reference noise diode for the Rcvr1_2 system.

Figure 3: Intrinsic Stokes parameters of the noise diode reference signal for Rcvr_800, plotted as a function of observing frequency for three distinct epochs: MJDs 56244, 56419, and 56608. The modeled values of Stokes $Q$, $U$, and $V$ are expressed as percentages of the total intensity of the reference source.

Therefore, we see that for both the IFA and MEM approaches for polarization calibration, the assumptions made regarding the feeds and the reference noise diode signal are not entirely applicable to the GBT receivers. Our results suggest that the TOAs generated from IFA-calibrated pulse profiles have fewer systematic errors than those from MEM calibration. This leads us to believe that the very small ellpticities and non-orthogonality of the receiver feeds and the reference signal being $\gtrsim 90$% linearly polarized do not significantly affect the quality of the IFA calibration. Comparatively, the temporal variations in the polarization of noise diode signal is significant for the MEM calibration and affects the accuracy of the polarization calibration more adversely. However, further in-depth investigations are necessary to validate these assertions and quantify how the deviations from the underlying assumptions for different polarization calibration methods would affect the accuracy of the calibration.

For MEM+METM polarization calibration, we have used the PSR B1937$+$21 as the reference pulsar to generate the per-epoch corrections to the MEM-generated PRs. However, we observed that the MEM+METM calibration typically performs worse than both the IFA and MEM calibrations. We suspect that this occurs because the polarization profiles of PSR B1937$+$21 do not remain stable over time. One potential cause of this instability could be variable scatter broadening, leading to temporal profile variations, particularly at lower frequencies (see Brook et al., 2018).

In previous NANOGrav data releases (e.g., NANOGrav Collaboration et al., 2015; Arzoumanian et al., 2020; Agazie et al., 2023), the IFA approach has been employed to calibrate observed pulsar profiles. Our analysis of three pulsars with data obtained with GBT Rcvr_800 and Rcvr1_2 receivers, coupled with the GUPPI backend system, indicates that IFA polarization calibration yields the best results. Therefore, it is recommended to continue using IFA polarization calibration for future NANOGrav data sets until any potential changes in the GBT receiver systems. However, conducting similar analyses for additional pulsars and also using data obtained from other telescopes, such as Arecibo and the Very Large Array (VLA), and with other receiver systems would be valuable to validate these findings. Further, using some other pulsar as a reference pulsar instead of B1937+21 can potentially help achieving better results for the MEM+METM calibration method. Additionally, we plan to investigate whether using matrix template matching instead of the standard total intensity to generate TOAs from calibrated profiles alters our results.

Author Contributions

L.D. carried out the analyses and prepared the text, figures, and tables. M.A.M. and H.M.W. helped with the development of the framework. P.B.D., S.B.S., W.F., J.G., and R.J.J. helped with useful inputs and discussions. Rest of the authors (including M.A.M., H.M.W., and P.B.D.) contributed to the collection and analysis of the NANOGrav 12.5 yr data set.

Appendix A Calibrated and uncalibrated profiles

In this appendix, we present the uncalibrated profiles as well as profiles obtained by different calibration methods for PSRs J1909$-$3744, J1643$-$1224, and J1744$-$1134. Figure 4 shows the uncalibrated and calibrated polarization profiles for J1909$-$3744 observed with Rcvr_800 and GUPPI backend system at GBT. Full polarization profiles for J1643$-$1224 and J1744$-$1134, for both Rcvr_800 and Rcvr1_2 observation, are shown in Figures 5 and 6, respectively.

Figure 4: Polarization profiles for PSR J1909$-$3744 observed with the GUPPI 800 MHz receiver system at GBT. Both uncalibrated profiles and profiles obtained by different calibration methods are shown.

Figure 5: Uncalibrated and different calibrated profiles for J1643$-$1224 observed with Rcvr_800 (800 MHz) and Rcvr1_2 (1500 MHz) receivers and GUPPI backend system at the GBT.

Figure 6: Uncalibrated and different calibrated profiles for J1744$-$1134 observed with Rcvr_800 (800 MHz) and Rcvr1_2 (1500 MHz) receivers and GUPPI backend system at the GBT.