HIGH-FIDELITY RADIO ASTRONOMICAL POLARIMETRY USING A MILLISECOND PULSAR AS A POLARIZED REFERENCE SOURCE

High-precision pulsar timing exploits the exceptional rotational stability of millisecond pulsars to measure orbital dynamics in binary systems and perform unique tests of general relativity in the strong field regime (e.g., Stairs 2006). The basic measurement in timing analyses is the pulse time-of-arrival (TOA), the epoch at which a fiducial phase of the pulsar's periodic signal is received at the observatory. The difference between the measured TOA and the arrival time predicted by a best-fit model is known as the timing residual. Timing residuals with a standard deviation of the order of 100 ns have been achieved for a growing number of millisecond pulsars (e.g., Lommen 2001; Hotan 2007), which has renewed both theoretical and experimental interest in the detection of gravitational radiation using a pulsar timing array (PTA; e.g., Foster & Backer 1990; Jaffe & Backer 2003; Wyithe & Loeb 2003; Jenet et al. 2006; Sesana et al. 2008; Yardley et al. 2011; van Haasteren et al. 2011).

The sensitivity of a PTA and the confidence limits that may be placed on any gravitational wave detection directly depend upon the precision and accuracy with which pulse arrival times can be estimated. Therefore, an important part of the ambitious PTA effort is the quantification and reduction of systematic error through the development of improved methods of arrival time estimation (e.g., Hotan et al. 2005; van Straten 2006), instrumental calibration (Jenet & Anderson 1998; van Straten 2004), radio-frequency interference (RFI) mitigation (Kocz et al. 2010; Nita & Gary 2010), compensation for the effects of propagation through the interstellar medium (e.g., You et al. 2007a; Hemberger & Stinebring 2008; Coles et al. 2010), statistical analysis of the stochastic nature of the pulsar signal (Manchester et al. 1975; Cordes & Downs 1985; Rathnasree & Rankin 1995; Osłowski et al. 2011), and the characterization of timing noise (Cordes & Downs 1985; Shannon & Cordes 2010).

At the Parkes Observatory, the majority of pulsar timing data are obtained from observations of PSR J0437−4715, the closest and brightest millisecond pulsar known. Discovered in the Parkes 70 cm survey (Johnston et al. 1993), PSR J0437−4715 has a spin period of ∼5.7 ms; at 20 cm, it has a pulse width at half-maximum of about 130 μs (Navarro et al. 1997) and an average flux of 140 mJy (Kramer et al. 1998). Owing to its sharp peak and large flux density, it is an excellent target for high-precision pulsar timing studies. It is also fortuitously located away from the Galactic plane (l = 245°, b = −42°), which permits long observing sessions during periods of local sidereal time when competition for the telescope is relatively low.

Early timing of this pulsar highlighted the need for more accurate polarimetric calibration (Sandhu et al. 1997). The systematic timing error due to instrumental polarization artifacts is readily observed in the dramatic variation of arrival time residuals as a function of parallactic angle (see Figure 1). This issue motivated the discovery and development of the polarimetric invariant profile (Britton 2000). Arrival times derived from the invariant profile proved to be significantly more accurate than those of the total intensity profile, a breakthrough that led to the first detection of annual-orbital parallax and a new test of general relativity (van Straten et al. 2001).

Zoom In Zoom Out Reset image size

Despite its demonstrated advantages over the total intensity profile, there remain a number of drawbacks to the use of the invariant profile. First, the signal-to-noise ratio (S/N) of the invariant profile is lower than that of the total intensity, reaching zero for completely polarized sources. Second, the noise in the invariant profile is neither normally distributed nor homoscedastic with respect to pulse longitude (see Appendix A); these two properties are assumed, either explicitly or implicitly, by most template-matching algorithms commonly used for pulsar timing. Finally, the computation of the invariant profile requires estimation and correction of the bias due to noise, and the relative precision of the bias estimate is inversely proportional to the square root of the number of discrete samples within off-pulse regions of pulse longitude. For a typical observation of PSR J0437−4715, the relative error in the bias correction is of the order of 8%, an unacceptably large uncertainty when aiming for arrival time uncertainty that is four orders of magnitude smaller than the pulse width.

These concerns motivated the development of matrix template matching (MTM; van Straten 2006), a technique that exploits the additional timing information in the polarized component of the pulsar signal while simultaneously eliminating residual polarimetric calibration errors. For some pulsars, the mean polarization profile exhibits sharper features than observed in the total intensity; these features generate more power at higher harmonics of the pulsar spin frequency, yielding greater arrival time precision than might be expected from the polarized flux alone (e.g., see Figure 1 of van Straten 2006). Furthermore, MTM improves arrival time accuracy by mitigating systematic timing errors due to instrumental polarization artifacts.

Although the development of MTM was primarily driven by the high-precision timing of PSR J0437−4715, van Straten (2006) predicted that PSR J1022+1001 would benefit most from this technique, including a ∼33% increase in timing precision and a significant decrease in systematic timing error. PSR J1022+1001 was contemporaneously discovered in both the Princeton-Arecibo Declination-Strip Survey (Camilo et al. 1996) and the Green Bank Northern Sky Survey (Sayer et al. 1997). Camilo et al. (1996) observed variations in the shape of its mean pulse profile on a timescale of the order of minutes and noted that the accuracy of arrival time estimates was limited by this effect. Mean profile variations over such timescales would be distinct from pulse-to-pulse shape variations such as giant pulses, microstructure, and drifting sub-pulses (e.g., Hankins 1971; Lyne et al. 1971; Backer 1973), which when averaged over many pulses can be described as phase jitter (Cordes & Downs 1985) or stochastic wide-band impulse-modulated self-noise (Osłowski et al. 2011). As it is typically assumed that the mean pulse profile does not change with time, detection of systematic profile variations in a millisecond pulsar would have a serious impact on the field of high-precision timing.

Accordingly, the stability of the PSR J1022+1001 pulse profile has been the subject of a number of detailed studies (Kramer et al. 1999; Ramachandran & Kramer 2003; Hotan et al. 2004a). These are reviewed in Section 2, following a mathematical treatment of the systematic timing error due to instrumental distortion of the total intensity profile. Section 3 describes the observing system and the novel procedures used to calibrate the PSR J1022+1001 data for this study, including a new scattered power correction algorithm that is described in more detail in Appendix B. As part of this analysis, a 7.2 year history of the polarimetric response at Parkes is presented and its impact on systematic timing error is discussed. In Section 4, the results of this study are compared and contrasted with previous work, the relevance of the calibration method to other pulsars is discussed, and potential future directions are considered.

In the following discussion, the polarization of electromagnetic radiation is described by the second-order statistics of the transverse electric field vector e, as represented by the four Stokes parameters:

$$\begin{equation} S_k = \langle{\bm e}^\dagger {\bm \sigma }_{k} {\bm e}\rangle. \end{equation} \tag{ 1 }$$

The angular brackets denote an ensemble average; e† is the conjugate transpose of e; ${\bm \sigma }_{0}$ is the 2 × 2 identity matrix; ${\bm \sigma }_{1}$, ${\bm \sigma }_{2}$, and ${\bm \sigma }_{3}$ are the Pauli spin matrices; S₀ is the total intensity, and S = (S₁, S₂, S₃) is the polarization vector (e.g., Britton 2000). The reception and propagation of polarized radiation is described by linear transformations of e, as represented using complex 2 × 2 Jones matrices. Any non-singular matrix can be decomposed into the product of a unitary matrix and a self-adjoint matrix. Unitary transformations result in a three-dimensional Euclidean rotation of the Stokes polarization vector, leaving the total intensity unchanged and preserving the degree of polarization. Self-adjoint Jones matrices effect a Lorentz boost of the Stokes four-vector, mixing total intensity and polarized flux (Britton 2000). The boost component, also known as the poldistortion (Hamaker 2000), of the instrumental response distorts the shape of the total intensity pulse profile and introduces systematic arrival time errors. Physically, instrumental boost transformations are caused by differential amplification of the signals from the two orthogonally polarized receptors and by cross-coupling between these receptors.

In van Straten (2006), the systematic timing error induced by a given level of polarization distortion is predicted for a number of millisecond pulsars using a numerical simulation in which the template profile for each pulsar is copied and subjected to a boost transformation. The orientation of the boost axis in Poincaré space is varied, and the minimum and maximum induced phase shifts between the distorted total intensity profile and the template profile are recorded. It is also possible to predict the systematic timing error due to polarization calibration errors by computing the rate at which the best estimate of the phase shift varies with respect to the instrumental boost transformation parameters. Beginning with Equation (12) of van Straten (2006), the best estimate of the phase shift φ derived using only the total intensity is given by

$$\begin{equation} \varphi = { { \sum \phi _{0,m} |S_{0,m}| \nu _m } \over { 2\pi \sum |S_{0,m}| \nu _m^2 } }. \end{equation} \tag{ 2 }$$

Here, S_{0, m} = S'*₀(ν_m)S₀(ν_m) is the cross-spectral power of the observed total intensity S'₀ and the template total intensity S₀ in the Fourier domain, ϕ_{0, m} is the argument of S_{0, m}, and the summations in the numerator and denominator are performed over m = 1 to m = N, where N is the maximum harmonic at which the fluctuation power spectrum exhibits significant power.

Now consider an observed total intensity profile that is a copy of the template subjected to an instrumental boost transformation, as parameterized in Section 4.1 of van Straten (2004) such that sinh²β = b · b. Referring to Equation (10) of Britton (2000), the transformed total intensity is given by

$$\begin{equation} S_0^\prime = (1+2\,{\bm b\cdot {\bm b}}) S_0 + 2 b_0 \, {\bm b\cdot {\bm S}}, \end{equation} \tag{ 3 }$$

where $b_0 = \cosh \beta = \sqrt{1+{\bm b\cdot {\bm b}}}$. The partial derivatives of S'₀ with respect to the boost parameters b_k are used to compute the partial derivatives of the phase shift with respect to these parameters when the boost parameter β is zero,

$$\begin{equation} \dot{\varphi }_k \equiv {\partial \varphi \over \partial b_k }\bigg |_{\beta =0} = {{ \sum {\rm Im}[S_k^*(\nu _m) S_0(\nu _m)] \nu _m } \over { \pi \sum |S_0(\nu _m)|^2 \nu _m^2 }}. \end{equation} \tag{ 4 }$$

The above expression defines the components of a three-dimensional gradient ${\dot{\bm {\varphi }}}$. The magnitude of the gradient $\dot{\varphi }_\beta \equiv |\mbox{$\dot{\bm \varphi }$}|$ provides a measure of the susceptibility of arrival time estimates from a given pulsar to instrumental distortion.

To aid in the physical interpretation of the susceptibility $\dot{\varphi }_\beta$, the first-order approximations to Equations (25) and (32) of van Straten (2006) are used to define the systematic timing error induced by either a 1% differential gain error or receptor non-orthogonality of 0.01 rad (∼06),

$$\begin{equation} \tau _\beta \equiv 5\times 10^{-3} \dot{\varphi }_\beta P, \end{equation} \tag{ 5 }$$

where P is the pulsar spin period. Values of τ_β for each of the 20 millisecond pulsars that are regularly observed as part of the Parkes Pulsar Timing Array (PPTA; Manchester et al. 2012) project are reported in Table 1. The quantities in this table are derived from the high S/N polarization profiles presented by Yan et al. (2011). For the majority of pulsars in the PPTA, a modest degree of instrumental distortion induces systematic timing error of the order of 100 ns, which is sufficient to seriously impede progress toward the detection of the stochastic gravitational wave background (Jenet et al. 2005).

Table 1. Relative Arrival Time Uncertainties and Systematic Timing Error

Pulsar	$\hat{\sigma }_{\varphi }$	$\hat{\sigma }_{\tilde{\varphi }}$	τβ
			(ns)
J0437−4715	0.85	1.43	207
J0613−0200	0.92	1.46	59
J0711−6830	0.88	1.54	81
J1022+1001	0.68	1.65	282
J1024−0719	0.74	2.11	34
J1045−4509	0.88	1.48	338
J1600−3053	0.90	1.39	115
J1603−7202	0.85	1.55	142
J1643−1224	0.91	1.40	266
J1713+0747	0.85	1.58	6
J1730−2304	0.71	1.70	198
J1732−5049	0.96	1.38	185
J1744−1134	1.56	6.43	105
J1824−2452A	0.88	2.56	18
J1857+0943	0.89	1.43	124
J1909−3744	1.02	1.51	22
J1939+2134	0.95	1.49	44
J2124−3358	0.85	1.45	127
J2129−5721	1.15	1.61	211
J2145−0750	0.95	1.44	147

Download table as:  ASCIITypeset image

To first order, the systematic timing error due to instrumental distortion is given by the inner product, $\mbox{${\bm b \cdot \dot{\bm \varphi }}$}$. In the polar coordinate system best suited to describing linearly polarized receptors, $\mbox{ ${\bm b}$}\simeq (\gamma ,\delta _\theta ,\delta _\epsilon)$, where γ parameterizes the differential receptor gain, and δ and δ_θ describe the non-orthogonality of the receptors (see Section 3.3 of van Straten 2006). These properties of the instrument may vary as a function of time for a variety of reasons. Diurnal variations are introduced by the parallactic rotation of the receiver with respect to the sky, as shown in Figure 1. Furthermore, step changes in differential gain are introduced when the levels of attenuation applied to each of the two orthogonal polarizations are reset at the start of each observation and/or periodically updated during the observation. In contrast, the cross-coupling of the receptors is typically assumed to remain relatively stable over timescales of the order of months.

Of particular interest to a long-term timing experiment such as a PTA are any step changes in b due to modifications of the instrumentation and/or slow variations in b due to gradual degradation of one or more system components. The temporal variations in distortion-induced error due to the long-term evolution of a given instrument will be correlated between the different pulsars observed using that system. In an array of N pulsars, the spectral structure of the correlated systematic error due to instrumental distortion is described by

$$\begin{equation} {\bf C}_\tau (f) = {\bm \Upsilon} {\bf C}_\beta (f) {\bm \Upsilon }^T, \end{equation} \tag{ 6 }$$

where C_τ is the N × N matrix of the auto- and cross-power spectra of the timing residuals (e.g., Yardley et al. 2011), $\mbox{ ${ \bm \Upsilon }$}$ is the N × 3 Jacobian matrix,

$$\begin{equation} \Upsilon _j^k \equiv {\partial \tau _j\over \partial b_k} = P_j {\partial \varphi _j\over \partial b_k} \end{equation} \tag{ 7 }$$

(P_j is the spin period of the jth pulsar) and C_β(f) is the 3 × 3 matrix of the auto- and cross-power spectra of the instrumental distortion parameters (the components of b).

In practice, an optimally weighted sum of the cross-power spectral harmonics is used to define a detection statistic with maximum sensitivity to a stochastic background of low-frequency gravitational waves (e.g., Anholm et al. 2009; Yardley et al. 2011). Given only the mean polarization profile of each pulsar in the timing array, it is possible to predict the impact of correlated systematic timing error due to instrumental distortion on such a detection statistic. First, only the dominant harmonic of the cross-power spectrum of two pulsars with distortion gradients $\mbox{ ${ \dot{\bm \varphi }}$}_A$ and $\mbox{ ${ \dot{\bm \varphi }}$}_B$ is considered. Next, it is assumed that the instrumental distortion parameters are uncorrelated and homoscedastic, such that C_τ = σ²_τI, where I is the 3 × 3 identity matrix. Following these simplifications, a first-order approximation of the coefficient of correlated systematic timing error is given by

$$\begin{equation} c_{AB} = { \mbox{${ \dot{\bm \varphi }}$}_A\, \mbox{${\bm \cdot }$}\, \mbox{ ${ \dot{\bm \varphi }}$}_B \over |\mbox{ ${ \dot{\bm \varphi }}$}_A| |\mbox{ ${ \dot{\bm \varphi }}$}_B| }. \end{equation} \tag{ 8 }$$

This correlation coefficient depends only on the shape of the mean polarization profile and has no modal structure on the celestial sphere. That is, the correlated timing error due to instrumental distortion will corrupt all of the moments in a multipole expansion of timing delay, which will impact on all of the primary goals of any PTA experiment (Foster & Backer 1990), including the long-term measurement of time (monopole moment; Guinot & Petit 1991), corrections to the solar system ephemeris (dipole moment; Champion et al. 2010), and detection of the gravitational wave background (quadrupole moment; Hellings & Downs 1983).

To illustrate the potential importance of correlated systematic timing error, the high S/N polarization profiles presented by Yan et al. (2011) are used to compute the cross-correlation coefficients for each of the 15 pulsar pairs presented in Figure 4 of Yardley et al. (2011). The cross-spectral power measurements of Yardley et al. (2011) are compared with the predicted values of c_AB in Figure 2; the coefficient of correlation between predicted and measured values is ρ ∼ 0.35. The agreement between predicted and experimental measures of correlation indicates that the systematic error due to instrumental distortion is a plausible contributor to the apparent anti-correlation between the results of Yardley et al. (2011) and the expected quadrupolar signature predicted by Hellings & Downs (1983).

Zoom In Zoom Out Reset image size

Also shown in Table 1 are the relative arrival time uncertainties $\hat{\sigma }_{\varphi }$ and $\hat{\sigma }_{\tilde{\varphi }}$ defined in van Straten (2006). The former is the ratio between the theoretical error in arrival time estimates derived using MTM and the error in estimates derived using only the total intensity; the latter is the ratio between the uncertainties in arrival times derived from the invariant interval and from the total intensity. The statistical uncertainty in the invariant profile arrival times is underestimated because the effects of heteroscedasticity and bias correction error are not considered. Regardless, for every pulsar in the PPTA, precision is lost through use of the invariant interval. The increase in uncertainty is most dramatic for PSR J1744−1134, which emits almost 100% linearly polarized radiation. For this pulsar, MTM is also predicted to yield arrival times with greater error than those derived from the total intensity owing to the large coefficient of multiple correlation (∼0.9) between the phase shift used to compute the TOA and the model parameters that describe the polarization transformation.

2.1. PSR J1022+1001

Table 1 also predicts that MTM will yield the greatest relative improvement in arrival time precision for PSR J1022+1001. This pulsar also has the second-highest value for predicted systematic timing error owing to the susceptibility of the total intensity profile to polarization distortion. This may explain the profile shape variations and excess arrival time error first noted by Camilo et al. (1996). In a detailed study of the spectrum and polarization of these variations using the Effelsberg 100 m Radio Telescope, Kramer et al. (1999) assert that the observed fluctuations in pulse shape cannot be explained by instrumental effects. However, the data used in this analysis were calibrated using the "ideal feed" assumption that there is no cross-coupling between receptors and that the reference source used for gain calibration produces an equal and in-phase response in each receptor. Given that receptor cross-coupling of the order of 10%–20% is commonly observed when more accurate calibration is performed (e.g., Stinebring et al. 1984; Xilouris 1991; van Straten 2004), an instrumental origin of the variability deserves consideration.

For instance, the temporal variation of the total intensity profile observed by Kramer et al. (1999) is reproducible using a simple model in which the polarized signal is rotated by the parallactic angle of the receiver and then subjected to an instrumental boost that mixes ∼10% of the linearly polarized flux with the total intensity. The similarity between the simulated observations plotted in Figure 3 and the total intensity profiles presented in Figure 2 of Kramer et al. (1999) supports the argument that the variations are the result of uncalibrated instrumental distortion.

Zoom In Zoom Out Reset image size

Similarly, the small frequency-scale variations observed by Kramer et al. (1999) may be readily explained by spectral variation of the instrumental response (e.g., Figure 3 of van Straten 2004). The observations used to demonstrate the narrowband variations and presented in Figure 8 of Kramer et al. (1999) were recorded using the Effelsberg–Berkeley–Pulsar–Processor with a bandwidth of 56 MHz. In this mode, only the fluxes in each of the circularly polarized signals (i.e., two out of four Stokes parameters) are retained; consequently, only the mixing between circularly polarized flux and total intensity can be calibrated in these data. Therefore, frequency-dependent mixing between linearly polarized flux and total intensity is a plausible explanation for the observed profile shape variations.

Ramachandran & Kramer (2003) analyzed the variability of PSR J1022+1001 using multi-frequency observations made at the Westerbork Synthesis Radio Telescope (WSRT). They argue that, because the WSRT is composed of equatorially mounted antennas, it is possible to rule out pulse shape variations due to instrumental distortion combined with parallactic rotation of the feed. However, the polarimetric response of a phased array can vary in both time and frequency for a wide variety of reasons, including imperfect phase calibration and (as for any single-dish antenna) instrumental gain variations. Indeed, the automatic gain controllers used to prevent instrumental saturation (Vote et al. 2002) also make it impossible to accurately calibrate the polarimetric response of the WSRT phased array (Edwards & Stappers 2004).

The high degree of susceptibility of the PSR J1022+1001 pulse profile to polarization distortion and the inadequate instrumental calibration employed in these previous works cast doubt on the conclusion that the profile variations are intrinsic to the pulsar. Furthermore, Hotan et al. (2004a) find no evidence of significant profile shape variations in PSR J1022+1001 observations made at the Parkes Observatory. This study achieved arrival time residuals with a standard deviation of the order of 2.3 μs and χ² per degree of freedom ∼1.4, whereas Kramer et al. (1999) obtained a reduced χ² close to unity only after systematic errors of the order of 10 μs were added in quadrature to the formal uncertainties.

The analysis by Hotan et al. (2004a) is based on observations integrated over 5 minutes and 64 MHz and spanning 1.3 years. Increasing the integration length to 60 minutes and the time span to 2.3 years yielded arrival time residuals with a standard deviation of the order of 1.5 μs and reduced χ² ∼ 3 (Hotan et al. 2006). With the same instrumentation and similar integration lengths, Verbiest et al. (2009) achieved 1.6 μs residuals in data spanning 5.1 years. In the absence of systematic error and other sources of noise (e.g., dispersion variations), the 1 hr integrations used in these two studies should have yielded timing residuals of the order of 660 ns. This is of the same order as the systematic timing error introduced by a differential gain error of ∼2.3% or by receptor non-orthogonality of ∼14 (cf. Table 1). However, the arrival times used in these studies were derived from either uncalibrated (Hotan et al. 2006) or inaccurately calibrated (Verbiest et al. 2009) total intensity pulse profiles, which are highly susceptible to instrumental distortion. Therefore, it is reasonable to anticipate that high-fidelity polarimetric calibration will improve the accuracy of PSR J1022+1001 arrival time estimates.

This conjecture is verified in the following section, which describes a new method of deriving the instrumental response by using PSR J0437−4715 as a polarized reference source. The method is used to calibrate high-precision timing observations of PSR J1022+1001 and demonstrated to reduce the systematic error due to instrumental polarization distortion.

Dual-polarization observations of PSR J0437−4715 and PSR J1022+1001 were made at the Parkes Observatory using the H-OH receiver and the center element of the Parkes 21 cm Multibeam receiver (Staveley-Smith et al. 1996). Two 64 MHz bands, centered at 1341 and 1405 MHz, were digitized using two bits per sample and processed using the second generation of the Caltech–Parkes–Swinburne Recorder (CPSR2; Bailes 2003; Hotan 2007). To maintain optimal linear response during digitization, the detected power was monitored and the sampling thresholds were updated approximately every 30 s. The baseband data were reduced using dspsr (van Straten & Bailes 2011), which corrects quantization distortions to the voltage waveform using the dynamic level setting technique (Jenet & Anderson 1998, hereafter JA98), then performs phase-coherent dispersion removal (Hankins 1971; Hankins & Rickett 1975) while synthesizing a 128 channel filterbank (van Straten 2003). The Stokes parameters are then formed and integrated as a function of topocentric pulse phase, producing uncalibrated polarization profiles with 1024 discrete pulse longitudes.

3.1. Minimizing Quantization Distortions

During analog-to-digital conversion, the radio signal is subjected to nonlinear distortions that significantly alter the observed pulse profile (JA98; Kouwenhoven & Vote 2001). JA98 introduce a dynamic output level setting technique that is employed by dspsr to correct the underestimation of undigitized power. To implement this technique, the digitized data are divided into consecutive segments of L samples and, for each segment, the number of low-voltage states M is counted. A histogram of occurrences of M is archived with the pulsar data and used offline to evaluate the quality of the digital recording.

When the voltage input to the digitizer is normally distributed, the ratio Φ = M/L has a binomial distribution as in Equation (A6) of JA98. The difference between this theoretical expectation and the recorded histogram of M provides a measure of the degree to which either the input signal is not well described by a normal distribution or the sampling thresholds diverge from optimality (e.g., at the start of an observation). This difference, called the two-bit distortion, is given by

$$\begin{equation} D = \sum _{M=0}^L [{\mathcal {P}}(M/L) - {\mathcal {H}}(M)]^2, \end{equation} \tag{ 9 }$$

where ${\mathcal {P}}(\Phi)$ is the expected binomial distribution and ${\mathcal {H}}(M)$ is the recorded distribution of M. Separate histograms of M are maintained for each polarization, and the reported distortion is simply the sum of the distortion in each polarization.

For each 16.8 s integration output by CPSR2, the two-bit distortion D was computed and all data with D > 3.5 × 10⁻⁴ were discarded. This threshold was chosen using the histograms plotted in Figure 4; beyond the range of plotted values, a sparse tail of outliers extends to a maximum value of D ∼ 0.98. Approximately 3.3% of the observations were rejected; the remaining data were corrected for scattered power using the method described in Appendix B and implemented in the psrchive software package for pulsar data archival and analysis¹ (Hotan et al. 2004b). The corrected data were then summed until the integration length was of the order of 5 minutes.

Zoom In Zoom Out Reset image size

3.2. Polarimetric Template Profile

To use a pulsar as a polarized reference source requires a high S/N, well-calibrated, template polarization profile. From over 700 hr of PSR J0437−4715 observations were selected only those sessions with low levels of two-bit distortion and of sufficient duration to achieve a precise estimate of the instrumental response. About 50 sessions, each ∼8 hr in duration and with corresponding pulsed noise diode observations, were calibrated using measurement equation modeling (MEM; van Straten 2004). The observed Stokes parameters were normalized to account for pulsar flux variations as described in Appendix A. As in van Straten (2004), the two degenerate model parameters (a boost along the Stokes V axis and a rotation about this axis) were constrained by assuming that observations of 3C 218 (Hydra A) have negligible circular polarization and that the orientation of the receiver is known. Although these assumptions are not necessarily valid, absolute certainty in the Stokes parameters is not required for the purposes of a high-precision timing experiment. It is sufficient that the observed Stokes parameters can be mapped to the template profile by a single Jones transformation.

This prerequisite may be evaluated using the objective merit function of the MTM least-squares fit; the reduced χ² statistic will be greater than unity whenever the transformation between the template and the observed pulse profile is not well described by a single Jones matrix. Various phenomena may contribute to a poor model fit. For example, if the mean polarized emission from the pulsar evolves over time, then the difference between the template and the observed profile will vary as a function of pulse longitude in a manner that cannot be adequately modeled by a single Jones transformation. Furthermore, as Jones matrices describe only linear transformations of the electric field, MTM will fail to model any nonlinear component of the instrument.

Analog-to-digital conversion using only two bits per sample is an intrinsically nonlinear process; therefore, variations in the response of the digitizer adversely affect the MTM goodness of fit, as shown in Figure 5. Here, 5 of the ∼50 MEM-calibrated PSR J0437−4715 observations were integrated to produce a template profile and the mean polarization profiles from each of the original ∼50 observations were fit to this template using MTM. The reduced χ² of each MTM fit (averaged over all frequency channels) is highly correlated with the median value of the two-bit distortion D in each session; i.e., greater distortion reduces the merit of the MTM fits. For the five observations used to create the template profile, χ²/N_free < 1 because covariance between the observation and the template is not included in the definition of the MTM objective merit function.

Zoom In Zoom Out Reset image size

Figure 5 justifies the experimental design decision to regularly update the two-bit sampling thresholds during each observation. Independent sampling thresholds were applied to each of the two orthogonal polarizations; consequently, the astronomical signal was subjected to an unknown differential gain that varies with time. Such differential gain variations can be modeled using MTM, whereas the nonlinear response of the digitizer cannot. Therefore, when recording the baseband signal with a two-bit digitizer, it is more important to maintain optimal linear response than it is to maintain a constant flux scale.

Out of the ∼50 sessions that were calibrated using MEM, 16 were selected with low median two-bit distortion and MTM reduced χ² within 4% of unity. The MTM-corrected mean profiles from these 16 sessions were integrated to form a polarimetric template profile with a total integration length of ∼100 hr, shown in Figure 6. The overpolarization seen from pulse phase 0.88 to 0.94 in this figure is an instrumental artifact that is currently not understood. It may be that the unpolarized quantization noise is overestimated by the two-bit scattered power correction algorithm described in Appendix B. However, a very similar artifact is observed in polarization data obtained with the second generation of the Parkes Digital Filter Bank (PDFB2; Yan et al. 2011), which employs an eight-bit digitizer. These facts might point to the existence of a nonlinear component in the signal path shared by CPSR2 and PDFB2. However, the third and fourth generations of the PDFB also regularly exhibit ∼5% overpolarization when it is not observed using the latest generation of baseband recording instrumentation at Parkes. Therefore, it may be only coincidence that a similar artifact is observed in both PDFB2 and CPSR2 data.

Zoom In Zoom Out Reset image size

3.3. Measurement Equation Template Matching

MTM has been previously utilized to calibrate and derive arrival time estimates from observations of PSR J0437−4715 (e.g., Verbiest et al. 2008). However, this method cannot be applied directly to PSR J1022+1001 because the mean flux of this pulsar at 20 cm (∼6 mJy; Manchester et al. 2012) is much lower and the precision of calibrator solutions derived from similar integration lengths is inadequate. It is not possible to integrate over longer intervals because the instrumental response varies with time, primarily due to the parallactic rotation of the receiver. Therefore, to calibrate the PSR J1022+1001 data, a new technique was developed that combines MTM with MEM. The new method, called Measurement Equation Template Matching (METM), matches multiple pulsar observations to the template profile, thereby increasing the precision of the solution. As with MEM, variations of the instrumental response (e.g., the parallactic rotation of the receiver) are included in the model. METM also incorporates measurements of an artificial reference source (e.g., a noise diode coupled to the receptors), enabling the backend component of the calibrator solution to be later updated as in Ord et al. (2004) so that the solution may be applied to other astronomical sources. In contrast with MEM, the METM method does not require any additional observations of a source of known circular polarization to constrain the boost along the Stokes V axis because the template profile provides this information. Using the new METM technique, a detailed history of the instrumental response at Parkes was derived from the PSR J0437−4715 data and used to calibrate the PSR J1022+1001 data.

First, the high S/N template profile (created as described in the previous section) and the new METM method were used to derive a calibrator solution from every observation of PSR J0437−4715 with more than 3.5 hr of integration. A total of ∼350 solutions were produced, spanning 2003 April to 2010 June. The instrumental response was modeled using Equation (19) of Britton (2000), which includes three independent Lorentz boost transformations described by the differential gain and the non-orthogonality parameters, δ and δ_θ. The derived non-orthogonality parameter estimates are plotted in Figure 7. In addition to smooth variations in δ_θ and δ, there are a number of distinct steps that correspond to physical changes in the 21 cm Multibeam receiver and downconversion system at Parkes. The gaps in the data correspond to maintenance periods during which the Multibeam receiver was unavailable (2003 November 9 to 2004 September 7, and 2006 December 16 to 2007 May 17). Immediately following the first maintenance period, the receptor cross-coupling remains stable for a period of approximately one year. The distinct changes in state that occur around MJD 53660 (2005 October 17) and MJD 53800 (2006 March 6) are due to changes in the downconversion system (the signal paths for the two CPSR2 bands were swapped). Following the second maintenance period, the cross-coupling parameters are notably less stable in both time and frequency. Figure 8 provides a closer look at the variations of these parameters over a period of one week and a bandwidth of 20 MHz. In addition to significant quasi-periodic variations in these parameters as a function of radio frequency, the spectral structure appears to drift as a function of time. For example, from 2010 February 11 to 12, the main features in the spectrum have shifted toward lower frequencies by ∼1.5 MHz. After MJD ∼ 54300, it is invalid to assume that the receptor cross-coupling remains stable for any period of time longer than a day.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 plots the non-orthogonality parameters in units of centiradians; these values can be multiplied by the estimates of τ_β listed in Table 1 to approximate the variation in systematic arrival time error introduced by receptor cross-coupling. For example, for PSR J1022+1001, the peak-to-peak variations of ∼0.06 rad in both δ_θ and δ are predicted to induce systematic timing errors of the order of 1.7 μs. However, after integrating over all frequency channels, the net effect of the cross-coupling variations will be diminished. Furthermore, receptor cross-coupling describes only two of the three forms of instrumental distortion. In addition to the presumption of zero receptor cross-coupling, calibration techniques based on the ideal feed assumption (IFA) typically also incorporate the false premise that the artificial reference source illuminates both receptors identically. Under this assumption, any intrinsic imbalance in the induced amplitude of the reference source signal in each receptor will be misinterpreted as differential gain. In contrast, the METM technique includes direct estimation of the Stokes parameters of the reference source (as in Figure 2 of van Straten 2004). Using these estimates, the systematic error in differential gain δ_γ due to unbalanced reference source amplitudes may be derived.

Given estimates of δ_γ, δ_θ, and δ as a function of frequency, it is possible to directly compute the systematic arrival time error at each epoch by applying the instrumental distortion transformation to a copy of the template profile and estimating the induced phase shift between the template and its distorted copy. Application of the METM-derived instrumental distortion parameters to the PSR J1022+1001 template profile yields the systematic timing error shown in Figure 9. The instrumental distortion transformation is applied independently in each frequency channel before integrating over frequency. The induced phase shift is computed using only the total intensity profile and the Fourier-domain template-matching algorithm described by Taylor (1992). During the periods of maintenance on the Parkes 21 cm Multibeam receiver, the solutions derived from observations made with the H-OH receiver have been used. Uncalibrated instrumental distortion introduces peak-to-peak variations in systematic timing error of the order of 600 ns; the largest jumps in timing error occur when switching between receivers. These systematic errors are present in arrival times derived from either uncalibrated or inaccurately calibrated PSR J1022+1001 observations. As shown in the following section, these errors are significantly reduced after calibration with the METM-derived solutions.

Zoom In Zoom Out Reset image size

3.4. Arrival Time Estimation and Analysis

The solutions produced using the new METM method were applied to calibrate the 5 minute PSR J1022+1001 integrations as in Ord et al. (2004). The calibrated data were then integrated in frequency and time to form approximately 64 minute integrations. The observations with the greatest S/N were added to form a full-polarization template profile with an integration length of ∼67 hr, shown in Figure 10.

Zoom In Zoom Out Reset image size

To enable a quantitative comparison between the analysis described in this paper and previous work, the PSR J1022+1001 observations were also calibrated using the IFA employed by Hotan et al. (2004a). Arrival times for the two data sets were then derived using only the total intensity profile (Taylor 1992) and MTM (van Straten 2006). These variations produced four sets of arrival time estimates, calibrated using either full measurement equation template matching (METM) or the IFA and timed using either MTM or standard total intensity (STI) methods. The post-fit arrival time residuals for these four sets are compared in Table 2.

Table 2. Arrival Time Residual Standard Deviation and Merit

Method	στ	χ2/Nfree
	(μs)
IFA–STI	1.9	2.2
METM–STI	1.4	2.1
IFA–MTM	0.92	1.1
METM–MTM	0.88	1.0

Download table as:  ASCIITypeset image

After calibration based on the IFA and arrival time estimation with the standard method of template matching using only the total intensity profile (STI), the residuals are of similar quality to those presented by Hotan et al. (2006). That is, compared to uncalibrated data, the IFA-calibrated data are equally adversely affected by distortions to the total intensity profile. This is not surprising, given that the IFA does not apply to the 21 cm Multibeam receiver at Parkes (van Straten 2004). Calibration via METM reduces both the standard deviation of the arrival time residuals and the reduced χ² of the model fit. Assuming that the systematic error removed by METM is uncorrelated with the remaining noise in the METM–STI residuals, instrumental distortion of the IFA-calibrated total intensity profile gives rise to timing errors with a standard deviation of approximately 1.3 μs.

For both METM- and IFA-calibrated data sets, MTM yields arrival times of similar quality. This indicates that MTM is able to correct the remaining instrumental calibration errors in the IFA data and produce arrival time estimates with reduced systematic error. It is not possible to test MTM on uncalibrated data because the variation of instrumental response with frequency leads to bandwidth depolarization, which cannot be modeled using a single Jones matrix. That is, although inaccurate, calibration based on the IFA may be sufficient to avoid depolarization.

Focusing on only the METM-calibrated data, comparison of arrival times generated using STI and MTM demonstrates that MTM reduces the standard deviation of arrival time residuals by approximately 35%, which is consistent with the relative uncertainty predicted in Table 1 and by van Straten (2006). When compared to the results of the typical data analysis employed in most high-precision timing experiments to date (IFA–STI), MTM reduces arrival time residuals by over 50%. The corresponding improvement in the experimental sensitivity of this data set is equivalent to that of quadrupling the integration length of each observation.

3.5. Results

The results presented in this section are based on the data calibrated using METM and arrival time estimates derived using MTM. Table 3 lists the best-fit timing model parameters obtained using the tempo2 analysis software (Hobbs et al. 2006; Edwards et al. 2006). These results include a significant detection of the Shapiro delay; after subtracting two degrees of freedom from the least-squares fit with the addition of the shape s = sin i and range r = m₂ parameters, χ² is reduced by ∼20% (from ∼330 to ∼260). Although the detection is significant, the elongated contours of constant Δχ² shown in Figure 11 demonstrate that s and r are highly covariant; the coefficient of correlation between these two model parameters is −0.99. Furthermore, the constraint on the companion mass is not very informative; the 95% confidence interval given by the projection of the Δχ² = 4 contour onto the m₂ axis ranges from 0.25 to >3 solar masses.

Zoom In Zoom Out Reset image size

Table 3. Best-fit Model Parameters for PSR J1022+1001

Parameter Name	Value
Reference clock	TT(TAI)
Planetary Ephemeris	DE414
P epoch (MJD)	53869.00016629324189
Pulse period, P (ms)	16.4529299518323(2)
Period derivative, $\dot{P}$ (10−20)	4.33399(6)
Ecliptic longitude, λ (°)	153.865881(3)
Ecliptic latitude, β (°)	−0.06(4)
Proper motion in λ (mas yr−1)	−15.97(2)
Proper motion in β (mas yr−1)	38(19)
Annual parallax, π (mas)	1.4(2)
Binary Model	DDH
Orbital period, Pb (days)	7.805135(1)
Projected semimajor axis, x (lt-s)	16.76541(4)
Orbital eccentricity, e (10−5)	9.702(7)
Epoch of periastron, T0 (MJD)	53876.1038(2)
Longitude of periastron, ω (°)	97.757(9)
Orthometric amplitude, h3 (μs)	0.52(7)
Orthometric ratio, ς	0.5(1)
Periastron advance, $\dot{\omega }$ (° yr−1)	0.009(3)
Secular variation of x, $\dot{x}$ (10−15)	14(1)

Note. For each parameter, the formal uncertainty (one standard deviation) in the last digit quoted is given in parentheses.

Download table as:  ASCIITypeset image

The shape and range parameters are highly covariant because, in nearly circular orbits that are only moderately inclined with respect to the line of sight, the Shapiro delay is readily absorbed in the Roemer delay by modification of the Keplerian orbital parameters. The unabsorbed remnant of the Shapiro delay is best described as a weighted sum of harmonics of the orbital frequency using the orthometric parameterization introduced by Freire & Wex (2010). Accordingly, Table 3 lists the best-fit values of the ratio of the amplitudes of successive harmonics:

$$\varsigma = { \sin i \over 1 + |\cos i| }$$

and the amplitude of the third harmonic

$$h_3 = m_2 \varsigma ^3.$$

The values of s and r derived from ς and h₃ are nearly identical to the values yielded by directly modeling the shape and range of the Shapiro delay. However, the coefficient of correlation between h₃ and ς is only −0.38. There is close agreement between the orbital inclination angle, i ∼ 57°, constrained by the Shapiro delay and previous estimates of the inclination of the pulsar spin axis, ζ ∼ 61° (Xilouris et al. 1998) and ζ ∼ 55° (Ramachandran & Kramer 2003), constrained using the rotating vector model (Radhakrishnan & Cooke 1969; Everett & Weisberg 2001). This agreement is compatible with the expectation that the pulsar spin and orbital angular momentum vectors are aligned during the period of mass accretion that led to the formation of the millisecond pulsar (e.g., Thorsett & Chakrabarty 1999).

The best-fit model parameters also include the first significant detections of the precession of periastron $\dot{\omega }$ and the secular variation of the projected semimajor axis of the orbit, $\dot{x}$. The estimate of $\dot{\omega }\sim (9\mbox{$^\circ $}\,{\pm}\, 3\mbox{$^\circ $})\times 10^{-3}$ yr⁻¹ is consistent with the value predicted by General Relativity, $\dot{\omega }^{\rm GR}\sim 8\mbox{$^\circ $}\times 10^{-3}$ yr⁻¹, based on the masses of the pulsar and its companion yielded by the Shapiro delay detection. The measured value of $\dot{x}$ is seven orders of magnitude larger than the general relativistic prediction, $\dot{x}^\mathrm{GR} \sim 5 \times 10^{-21}$, because it is dominated by the secular variation due to proper motion (Kopeikin 1996). This geometric effect can be used to place an upper limit on the orbital inclination angle; i.e., $\tan i < \mu { x / \dot{x} }$, where μ = 42 ± 18 mas yr⁻¹ is the total proper motion. The resulting constraint, i < 83°, rules out only 12% of possible orientations and is consistent with the stronger constraint yielded by measurement of the Shapiro delay.

The proper motion of the system also induces a quadratic Doppler shift (Shklovskii 1970) that contributes Tμ²d/c to observed period derivatives, where T is the period (e.g., spin or orbital), d is the distance to the pulsar, and c is the speed of light. The magnitude of the Shklovskii effect is estimated using the distance d_π = 750⁺¹⁰⁰₋₉₀ pc derived from the annual trigonometric parallax after correction for Lutz–Kelker bias, π = 1.340^+0.175_−0.165 mas (Verbiest et al. 2010). The Shklovskii contribution to the spin period derivative (∼4.9 × 10⁻²⁰) is slightly larger than its measured value, which would imply that the intrinsic $\dot{P}$ of the pulsar is negligible. However, the predicted contribution to the orbital period derivative (∼2 × 10⁻¹²) is around 20 times larger than the upper limit derived by including $\dot{P}_\mathrm{b}$ in the timing model. The relative uncertainties in the above estimates of the Shklovskii effect are almost 100%, primarily due to the large uncertainty in the estimate of the proper motion in ecliptic latitude. Using only the proper motion in ecliptic longitude, the predicted kinematic contribution to $\dot{P}_\mathrm{b}$ is only three times larger than the timing-derived upper limit.

To explain the remaining discrepancy, three other contributions to the observed $\dot{P}_\mathrm{b}$ are considered (as in Damour & Taylor 1991). First, the loss of energy due to gravitational radiation as predicted by General Relativity contributes around −3 × 10⁻¹⁶, approximately four orders of magnitude smaller than the Shklovskii effect. Second, the contribution due to differential rotation about the Galactic center (approximately −7 × 10⁻¹⁵) can also be neglected. Finally, the acceleration of the PSR J1022+1001 system in the Galactic gravitational potential is considered. Using the parallax distance d_π and the Galactic latitude b = 51°, the component of gravitational acceleration perpendicular to the Galactic disk, K_z ∼ 6.6 × 10⁻¹¹ m s⁻² (Bahcall 1984), contributes approximately −1.2 × 10⁻¹³ to the observed $\dot{P}_\mathrm{b}$. At large heights above the Galactic disk, the relative uncertainty in K_z is more than a factor of two; therefore, this contribution is large enough to negate the Shklovskii effect, resulting in an observed $\dot{P}_\mathrm{b}$ that is consistent with zero.

Instrumental distortion of the total intensity profile introduces significant systematic timing errors that are correlated between pulsars observed with the same instrument. Therefore, high-fidelity polarimetry is inextricably linked with the long-term objectives of high-precision timing and the primary goals of PTA experiments. This paper presents a novel method of calibrating the instrumental response by exploiting the long-term stability of the polarized emission from a millisecond pulsar. Full-polarization timing observations of PSR J0437−4715 are used to model variations in the 21 cm Multibeam and H-OH receivers and downconversion system at Parkes. Over a period of ∼7.2 years, temporal and spectral variations of the receptor cross-coupling introduce systematic error of the order of 1 μs in arrival time estimates derived from observations of PSR J1022+1001. High-fidelity calibration of instrumental polarization and arrival time estimation using MTM is demonstrated to correct systematic error and double the sensitivity of the experiment, yielding significant detections of the Shapiro delay, the precession of periastron, and the secular variation of the projected semimajor axis of the orbit. The improvements in both timing accuracy and precision are consistent with the hypothesis that the average polarized emission from these millisecond pulsars is stable over the timescales of relevance to this experiment.

With a modest instrumental bandwidth of 64 MHz and median integration length of 64 minutes, PSR J1022+1001 yields arrival time residuals with an uncertainty-weighted standard deviation of only σ_τ = 880 ns, roughly five orders of magnitude smaller than the pulsar spin period P ∼ 16.5 s. The remarkable timing precision of these observations supports the conclusions of Hotan et al. (2004a), who argue that the PSR J1022+1001 profile variations observed by Kramer et al. (1999) and Ramachandran & Kramer (2003) are not intrinsic to the pulsar and are most likely due to instrumental calibration errors.

The precision of the arrival time data presented in this paper places PSR J1022+1001 among the top 10 sources regularly observed as part of the PPTA project. Other PPTA sources may also benefit from the application of METM; e.g., six of the sources listed in Table 1 have instrumental distortion susceptibility factors ranging from 200 ns to 340 ns. Referring to the measured values of the non-orthogonality parameters plotted in Figure 7, these factors roughly correspond to systematic timing error variations of the order of 1 μs. Therefore, the timing precision of these pulsars would likely be improved by application of the methods developed for this study. For PSR J1744−1134 and PSR J2129−5721, MTM is predicted to yield arrival times with greater uncertainty than those derived from the total intensity profile. In these two cases, the best results would be obtained by calibrating using METM and deriving arrival times using STI.

In a wider analysis of all of the pulsars in the PPTA, Manchester et al. (2012) calibrated approximately 5.9 years of PSR J1022+1001 timing observations using MEM-derived solutions. Using typical integration lengths of 64 minutes and an instrumental bandwidth of 256 MHz, this study yielded arrival time residuals with an uncertainty-weighted standard deviation of 1.7 μs and reduced χ² ∼ 9.3. The residual noise is approximately 2.4 times greater than expected based on simple extrapolation of the METM–STI results presented in Table 2. This discrepancy may be explained by the fact that the PDFB instruments have a nonlinear response that introduces overpolarization, which cannot be calibrated using MEM. The accuracy of MEM-based calibration is also limited by the fact that there is no unique solution to the measurement equation when the only constraining transformation is the geometric rotation of the receptors with respect to the sky (van Straten 2004). To constrain the otherwise degenerate boost along the Stokes V axis, it is necessary to include observations of a source that is assumed to have negligible circular polarization (e.g., for the MEM fits performed in Section 3.2, observations of Hydra A were used to constrain δ). However, the use of an unpulsed source of radiation as a constraint is intrinsically susceptible to variability in other contributions to the system temperature, including receiver noise and ground spillover. That is, it is safer to assume that the polarized emission from PSR J0437−4715 is constant than it is to rely on a source of unpulsed radiation during calibration. The results presented by Manchester et al. (2012) may also be limited by the instability of the receptor cross-coupling, as shown in Figure 7. The temporal variations of the non-orthogonality parameters are sufficiently resolved only via application of the METM method presented in this paper, which yields a seven-fold increase in the number of available calibrator solutions.

As noted in previous studies (Kramer et al. 1999; Ramachandran & Kramer 2003; Hotan et al. 2004a; Manchester et al. 2012), other phenomena may also currently limit the timing precision of PSR J1022+1001. For example, owing to its low ecliptic latitude, the radio signal from this pulsar is subject to significant dispersive delays in the solar wind (You et al. 2007b, 2012). No corrections for dispersion variations were applied to the data presented in this work; however, observations made when the line of sight to PSR J1022+1001 passes near the Sun (around late August of each year) were excluded from the data set. In addition to fluctuations in dispersion, the observed flux of the pulsar varies as a function of time and radio frequency due to scintillation in the interstellar medium. The average profile of PSR J1022+1001 varies significantly as a function of radio frequency (e.g., see Figure 1 of Ramachandran & Kramer 2003) and, when modulated by interstellar scintillation, the frequency-integrated mean profile may fluctuate with time. The potentially significant arrival time estimation errors induced by this effect could be mitigated through the use of a frequency-dependent template profile.

It is reasonable to expect that PSR J1022+1001 timing will be improved by the current generation of instrumentation at the Parkes Observatory. The data presented in this paper were observed using a system with low dynamic range that performed analog-to-digital conversion of the radio signal with only two bits per sample. Two-bit quantization is an intrinsically nonlinear process and the techniques applied to restore linearity (dynamic output level setting) and mitigate quantization noise (scattered power correction) are based on a linear approximation to the response of the digitizer (Jenet & Anderson 1998). As discussed in Appendix B, the scattered power correction algorithm applied in this analysis is accurate only to first order and it is possible that overestimation of the unpolarized scattered power contributes to the overpolarization noted in Figure 6. To overcome the limitations of two-bit sampling, a new baseband recording and processing system with greater dynamic range was commissioned at the Parkes Observatory in 2010 April. Designed in collaboration with the Center for Astronomy Signal Processing and Electronics Research (CASPER) at Berkeley, the CASPER–Parkes–Swinburne Recorder (CASPSR) digitizes a dual-polarization 400 MHz band with eight bits per sample and performs real-time RFI excision based on the spectral kurtosis estimator (Nita & Gary 2010). This system currently operates in parallel with the third and fourth generations of the Parkes Digital Filter Bank (PDFB) instruments as part of the PPTA program.

The calibration techniques applied in this paper could potentially be improved by reducing the number of degrees of freedom in the estimates of the receptor cross-coupling parameters. For example, this could be achieved by fitting an analytic model to the temporal and spectral variations of the calibrator solutions or by employing a lossy compression algorithm that can be applied to irregularly sampled data. The smoothing effected by such a transformation would reduce the instantaneous noise in the applied calibrator solutions and might also provide a robust means of interpolating between solutions. However, further refinements of the calibration technique may yield only marginal improvements to arrival time accuracy and precision. Table 2 demonstrates that, even when the data are calibrated using the inaccurate IFA, MTM yields arrival time estimates that are nearly as good as those derived from data calibrated using the new METM technique described in this paper. This seems to suggest that, as long as the signal is not depolarized by integrating over time and/or frequency, the fidelity of the calibration technique has negligible impact on arrival times derived using MTM.

This conclusion may have a significant impact on the design of the next generation of instrumentation for high-precision timing. For example, it has been asserted that in order to achieve timing precision of the order of 100 ns, the Square Kilometre Array (SKA) must achieve net polarization purity corresponding to β < 10⁻⁴ (Cordes et al. 2004). However, even when non-orthogonality as large as β ∼ 10⁻² is ignored, as is the case when the 21 cm Multibeam receiver at Parkes is assumed to be ideal, MTM yields accurate arrival time estimates. That is, by relaxing the requirement for polarization purity, the application of MTM has the potential to reduce the cost of SKA development for high-precision timing.

All of the software required to perform MEM, polarimetric calibration, and MTM is freely available as part of psrchive (Hotan et al. 2004b); the use of this software is demonstrated by van Straten et al. (2012) and more fully documented online.³

I am grateful to Matthew Bailes, Paul Demorest, Mike Keith, Michael Kramer, Stefan Osłowski, and John Reynolds for helpful discussions during this research project and for insightful comments that greatly improved this report. Joris Verbiest implemented and assisted with tempo2 support of ecliptic coordinates and the orthometric parameterization of the Shapiro delay. The Parkes Observatory is part of the Australia Telescope which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO.

The invariant profile is formed by computing the Lorentz interval S as a function of pulse longitude, where $S^2\equiv S_0^2-|\mbox{${\bm S}$}|^2$. Assuming that the noise in each of the Stokes parameters is independent and normally distributed with standard deviation ς, define the normalized Stokes parameters S'_k = S_k/ς. The noise in S'²₀ has a noncentral χ² distribution with one degree of freedom and noncentrality parameter λ = 〈S'₀〉². Likewise, the noise in $|\mbox{${\bm S}$}^\prime |^2$ has a noncentral χ² distribution with three degrees of freedom and $\lambda ={|}\langle \mbox{${\bm S}$}^\prime \rangle |^2$. The distribution of the noise in S² is given by the cross-correlation of the distributions of S'²₀ and |S'|²; examples are shown in Figure 12. The variance of S² is equal to 4||〈S'〉||² + 8, where ||S|| is the Euclidean norm of the Stokes four-vector, defined by ||S||² ≡ S²₀ + |S|². As all four Stokes parameters vary as a function of pulse longitude, the noise in S'² is heteroscedastic, which violates one of the basic premises of most template-matching algorithms.

Zoom In Zoom Out Reset image size

Both the distribution of the invariant interval and its tendency toward zero for highly polarized radiation adversely affect its usefulness in template matching (as a replacement for the total intensity) and as a normalization factor during MEM (van Straten 2004). Gain normalization is necessitated by interstellar scintillation, which causes the received flux of the pulsar to vary as a function of time and radio frequency. To address the problems with the invariant interval during the application of MEM in Section 3.2, the Stokes parameters were normalized by the sum of the invariant interval over all on-pulse longitudes. For the typical observation of PSR J0437−4715, this summation reduces the relative error of the normalization factor by ∼97%, thereby yielding more normally distributed normalized Stokes parameters.

The data presented in this paper were processed using a new scattered power correction algorithm that can be used to correct digitization distortions to folded pulsar profiles. In the case of two-bit sampled voltages, the scattered power (or quantization noise) is ∼12% of the total power. The derivation of the method starts with the mean digitized power $\hat{\sigma }^2$ given by Equation (A5) of Jenet & Anderson (1998, hereafter JA98),

$$\begin{equation} \hat{\sigma }^2 = \sum _\Phi {\mathcal {P}}(\Phi)f(\Phi), \end{equation} \tag{ B1 }$$

where Φ is the fraction of samples that fall between the chosen thresholds, f(Φ) is the digitized power as a function of Φ, given by Equation (A4) of JA98, and ${\mathcal {P}}(\Phi)$ is the discrete probability distribution for Φ, given by Equation (A6) of JA98. The parameter Φ is eliminated by the summation in the above equation; however, both f(Φ) and ${\mathcal {P}}(\Phi)$ are also parameterized by the expectation value 〈Φ〉. Therefore, Equation (B1) represents the relationship between the mean digitized power and the mean value of Φ. That is, given the mean digitized power $\hat{\sigma }^2$, Equation (B1) can be inverted to compute 〈Φ〉, which can in turn be used to estimate the mean undigitized power σ² and the mean scattered power 1 − A via Equations (45) and (43) of JA98, respectively.

Equation (A5) of JA98 can be inverted using the Newton-Raphson method and the partial derivatives of Equations (A4) through (A6) of JA98 with respect to 〈Φ〉. These are simplified in the case of two-bit sampling by noting that Equation (A4) of JA98 reduces to

$$f(\Phi) = \langle \Phi \rangle y_3^2(\Phi) + (1-\langle \Phi \rangle) y_4^2(\Phi),$$

where y₃(Φ) and y₄(Φ) are given by Equations (41) and (40) of JA98.

The folded profile scattered power correction algorithm is based on the following assumptions and approximations. First, at the time of recording the astronomical signal, it is necessary that the baseband voltages are sampled using the optimum two-bit input thresholds defined in Table 1 of JA98. To first order, this condition is satisfied by excluding data with excessive two-bit distortion as defined by Equation (9) of this paper. Second, Equations (45) and (43) of JA98 only approximately relate the expectation values of Φ, σ², and A. For example, the relationship between A and σ² defined by Equation (43) of JA98 is concave down (see Figure 3 of JA98); therefore, Jensen's inequality dictates that the value of A estimated from the mean undigitized power will be greater than the expectation value of A. Consequently, the mean fractional scattered power may be systematically underestimated. This limitation cannot be overcome without prior knowledge of the distribution of total intensity fluctuations intrinsic to the pulsar signal. Finally, it is assumed that the signal in each frequency channel has not been significantly altered, such that there is a well-defined relationship between the mean undigitized power and the mean digitized power. For example, after phase-coherent dispersion removal, the flux density over a given time interval no longer represents the voltage fluctuations in the digitizer over that interval; the previously smeared pulsar signal will be recovered and the scattered power will be smeared. To estimate the scattered power using coherently dedispersed digitized power, the dispersion smearing in each frequency channel must be less than or of the order of the time resolution of the folded profile. This condition is satisfied in the 20 cm timing observations of PSR J0437−4715 and PSR J1022+1001 presented in this paper.