OBSERVATIONS AND MODELING OF RELATIVISTIC SPIN PRECESSION IN PSR J1141−6545

Binary pulsars with short orbital periods and massive companions exhibit a range of relativistic effects, the most prominent of which are relativistic periastron precession, gravitational time dilation, transverse Doppler effects, orbit decay due to loss of energy to gravitational waves, and the Shapiro delay resulting from the passage of the pulsar signal through the gravitational field of the companion. These effects may be described by the so-called post-Keplerian parameters (Damour & Deruelle 1985, 1986) and are most readily observed in double-neutron-star systems such as the Hulse–Taylor binary pulsar, PSR B1913+16 (Weisberg & Taylor 2005) and the Double Pulsar, PSR J0737−3039A/B (Kramer et al. 2006). Another relativistic effect that is potentially observable in such systems is the precession of the pulsar spin axis resulting from coupling between the spin and orbital angular momenta (Damour & Ruffini 1974; Barker & O'Connell 1975). The angular rate of precession, Ω_p, is given by

$$\begin{equation} \Omega _p = \frac{G^{2/3}}{c^2} \left(\frac{P_b}{2\pi }\right)^{-5/3} \frac{m_c(4m_p + 3m_c)}{(1-e^2)(m_p + m_c)^{4/3}}, \end{equation} \tag{ 1 }$$

where G is the gravitational constant, c is the velocity of light, P_b is the binary orbital period, e is the orbital eccentricity, and m_p and m_c are the pulsar and companion masses, respectively. For PSR B1913+16, the expected rate from Equation (1) is 121 yr⁻¹, corresponding to a precessional period of about 300 years, whereas for PSR J0737−3039A/B the rate is about a factor of 4 larger and the precessional period is about 75 yr.

Since the emission from pulsars is believed to be beamed, probably along the open field lines associated with the magnetic poles on the neutron star, precession of the pulsar spin axis will result in changes in the beam aspect as viewed from the Earth. One would expect this to lead to changes in the observed pulse profile and changes attributed to this effect have indeed been observed in PSR B1913+16 (Weisberg et al. 1989; Kramer 1998) and PSR B1534+12 (Arzoumanian 1995; Stairs et al. 2004). For PSR B1913+16, a significant fraction of the expected beam diameter has been traversed in the time since its discovery and this has been used to map the two-dimensional structure of the emission beam (Weisberg & Taylor 2002; Clifton & Weisberg 2008). Because of its relatively short precessional period, similar effects were expected to be observable for the Double Pulsar, but surprisingly, they have not so far been observed in the A pulsar (Manchester et al. 2005; Ferdman et al. 2008). A possible reason for this is that A's spin–orbit misalignment angle (the angle between the spin and orbital angular momenta) is small (see also Willems et al. 2004; Stairs et al. 2006). However, large variations in the pulse profile and orbital visibility of pulsar B have been observed (Burgay et al. 2005) and these are certainly due to some combination of relativistic periastron precession and precession of the B spin axis. There is now direct evidence for the relativistic spin precession of pulsar B. The emission from pulsar A is eclipsed for about 30 s when it passes behind pulsar B (Lyne et al. 2004; Kaspi et al. 2004) with the eclipse profile showing modulation related to the spin period of pulsar B (McLaughlin et al. 2004). Modelings of the long-term variations in the eclipse profile are consistent with the relativistic precession of B's spin axis with a spin–orbit misalignment angle of about 130° (Breton et al. 2008).

PSR J1141−6545 is a 394 ms pulsar in an eccentric 4.7 hr binary orbit with a relatively massive companion (m_c ∼ 1.0 M_☉) discovered in the Parkes Multibeam Pulsar Survey (Kaspi et al. 2000). With these characteristics, relativistic effects are expected to be detectable and indeed, measurement of the relativistic advance of periastron was reported in the discovery paper. The expected rate of relativistic precession of the pulsar spin axis is 136 yr⁻¹, corresponding to a precessional period of about 265 yr. Unusually for such a binary system, the pulsar characteristic age is relatively low, about 1.4 × 10⁶ yr, suggesting that the initially more massive primary star evolved to form a white dwarf, the present companion star, and in the process transferred mass to the secondary, which subsequently exploded, leaving behind a neutron star, the present pulsar (Dewey & Cordes 1987; Kaspi et al. 2000; Tauris & Sennels 2000). Continued timing observations (Bailes et al. 2003; Bhat et al. 2008) reinforced these ideas and gave detections of relativistic time dilation and orbital decay which were consistent with the predictions of general relativity. A clear orbital modulation in the timescale for interstellar scintillation was detected by Ord et al. (2002), giving estimates of the orbital inclination angle i (or 180 − i) =760 ± 25 and the system transverse space velocity, 115 ± 10 km s⁻¹. From measurements of the relativistic Shapiro delay, Bhat et al. (2008) derive a value for the inclination angle of i = 730 ±  15, consistent with the scintillation measurement.

Observations at frequencies around 1.4 GHz by Hotan et al. (2005a) showed that the mean pulse profile for PSR J1141−6545 changed significantly over the five years from 1999 July to 2004 May (MJD range 51381–53134), with an approximately linear increase in W₁₀, the pulse width at 10% of the pulse peak. Furthermore, based on two measurements separated by about 0.7 yr, there appeared to be steepening of the gradient of polarization position angle (P.A.) dψ/dϕ, where ψ is the P.A. and ϕ is the pulse phase (measured in the same units as the P.A.) near the profile center, at a rate of 2.3 ± 0.4 yr⁻¹. Variations of P.A. across pulse profiles are often well described by the "rotating-vector model" (RVM) in which the emission is assumed to be polarized parallel (or perpendicular) to the projected direction of magnetic fields in the vicinity of a magnetic pole on the rotating neutron star (Radhakrishnan & Cooke 1969). In the RVM, ψ varies as

$$\begin{equation} \tan (\psi - \psi _0) = \frac{\sin \alpha \sin (\phi - \phi _0)}{\sin \zeta \cos \alpha - \cos \zeta \sin \alpha \cos (\phi - \phi _0)}, \end{equation} \tag{ 2 }$$

ψ₀ is the P.A. at ϕ₀, the center of symmetry of the P.A. variation and the projected direction of the rotation axis for a dipole field, ζ = α + β is the inclination of the observer's line of sight relative to the rotation axis of the star, α is the inclination of the magnetic axis relative to the rotation axis, and β is the impact parameter of the observer's line of sight, that is, the minimum angle between the magnetic axis and the observer's line of sight which occurs at ϕ = ϕ₀.⁹ From Equation (2),

$$\begin{equation} (d\psi / d\phi)_{\rm m} = \sin \alpha /\sin \beta, \end{equation} \tag{ 3 }$$

where (dψ/dϕ)_m is the maximum value of dψ/dϕ at ϕ₀. Therefore, the increase in P.A. gradient observed by Hotan et al. (2005a) for PSR J1141−6545 was interpreted as a decrease in |β|, implying that our line of sight was moving closer to the beam axis. This and the changing pulse width were interpreted in terms of relativistic precession of the pulsar spin axis with a spin–orbit misalignment angle greater than 15° and probably less than 30°.

In this paper, we report on observations of PSR J1141−6545 made over 9.3 yr from 1999 August to 2008 November with the Parkes radio telescope at frequencies close to 1.4 GHz. Since 2004 April, we have observed with systems recording full polarization data but before that, only total-intensity data were obtained. The observations and analysis procedures are described in Section 2 and results are presented in Section 3. In Section 4, we describe the interpretation of the observed variations in terms of relativistic precession of the pulsar spin axis and the implied shape of the emission beam is described in Section 5. The implications of our precessional model for the formation of the system are described in Section 6. In Section 7, we summarize results and give our conclusions.

We have observed PSR J1141−6545 using the Parkes 64 m radio telescope between 1999 August and 2008 November (MJDs 51411–54785) using (at different times) two receiving systems, either the center beam of the Parkes 20 cm multibeam (MB) receiver (Staveley-Smith et al. 1996) or the "H–OH" receiver. The MB receiver has a bandwidth of about 300 MHz centered at about 1.4 GHz and a system equivalent flux density (S_sys) of approximately 29 Jy. The H–OH receiver has a wider bandwidth, approximately 600 MHz from 1.2 to 1.8 GHz, but a somewhat higher S_sys of about 35 Jy. Both receivers have orthogonal linearly polarized feeds and a linearly polarized, broadband and pulsed calibration signal which can be injected into the feed at 45° to the two signal probes. Four different back-end systems were used, all of which recorded data across a total bandwidth of 256 MHz: an analog filter bank (AFB) having 512 frequency channels on each polarization and recording total-intensity data using a one-bit digitizer system (see Manchester et al. 2001), a wideband correlator (WBC) which recorded all four Stokes parameters with 1024 frequency channels, and two digital filterbank systems (PDFB1 and PDFB3) which also recorded full polarization data with 512 channels across the band. For the WBC and the PDFBs, each observation was preceded by a short (2 minutes) observation of the pulsed calibration signal. Observations of Hydra A, assumed to have a flux density of 43.1 Jy at 1.4 GHz, were used to set the flux-density scale. Off-line processing made use of the psrchive pulsar data analysis system (Hotan et al. 2004).

A log of the observations is given in Table 1. AFB observations, which commenced soon after the discovery of the pulsar, were recorded on tape and processed off-line to form mean total-intensity profiles. Frequency channels containing known interference were given zero weight before data were summed across the band. For the purposes of this paper, observations made over intervals ranging from a few hours to a few weeks have been grouped as shown in Table 1. Pulse widths and flux densities were determined for each observation and mean values and their uncertainties for the group determined from these. The flux-density scale for the AFB data was established by comparison with contemporaneous WBC observations.

WBC and PDFB observations were typically of 1 hr duration and were folded online to form 1 minute sub-integrations. Frequency channels containing strong interference, affected by resonances in the feed or with low system gain (normally 5% of the band at the band edges) and sub-integrations affected by strong impulsive interference were given zero weight. Using the pulsed calibration signal as a reference, the variations in instrumental gain and phase across the band were removed, Stokes parameters formed and the data placed on a flux-density scale. Data were summed to form 15 minute sub-integrations and mean profile widths and flux densities computed as for the AFB data. All profiles were recorded with 1024 bins across the pulse period.

Observed P.A.s given in this paper follow the astronomical convention, with P.A. measured from celestial north and increasing toward east (counterclockwise looking toward the source). The IEEE definition of circular polarization is adopted (i.e., at a given point in space and looking at the source, the E-vector for a right-circular wave rotates in a counterclockwise direction). Both of these definitions are in accord with IAU recommendations (Trans. IAU, 15B, 166, 1973). Stokes V is defined in the sense left-circular minus right-circular to conform with established pulsar polarization conventions which pre-date the IAU recommendations. File header parameters defined the polarization setup for each receiver/back-end combination and these were used by psrchive to correct the Stokes parameters so that they conformed to these definitions. For further details on polarization conventions and their implementation, see van Straten et al. (2010). The strong millisecond pulsar, PSR J0437−4715, was observed almost every session and used to check the sign of the calibrated Stokes parameters and the absolute value of P.A. by comparison with the results given by Navarro et al. (1997) which are known to be in accordance with these definitions. The linearly polarized intensity L = (Q² + U²)^1/2 is a positive definite quantity which is biased by noise. Plotted values have been corrected for this bias using the relations given by Lorimer & Kramer (2005).

Before summing in frequency, the observed pulse profiles must be corrected for Faraday rotation across the band. The nominal rotation measure (RM) for PSR J1141−6545 is −84 ± 2 rad m⁻² (Han et al. 2006). RMs were measured for each epoch having polarization data by summing the upper and lower halves of the bandpass separately using the nominal RM, taking a weighted mean of the P.A. differences between the upper and lower bands across the profile, recomputing the RM and then iterating until convergence. We show below that the observed P.A. variations across the pulse profile are complex, both in time and in frequency, with quite different behavior in the central and outer parts of the profile. Consequently, RMs were determined separately for the central and outer parts; we believe that the outer parts represent the true interstellar RM.

The MB receiver suffers from significant coupling between the nominally orthogonal feed probes. Observed polarization variations resulting from this coupling are a function of parallactic angle and the coupling parameters can be measured by analysis of a series of observations of a polarized source covering a wide range of parallactic angles (van Straten 2004). PSR J0437−4715 was observed every few months, typically for 10 minutes each hour during a 10 hr transit, for each of the WBC and PDFB configurations used with the MB receiver (Table 1). These data were analyzed using the psrchive program pcm to determine the feed ellipticities and their relative orientation as a function of frequency across the band. Observations of PSR J1141−6545 were then calibrated to remove the effects of the cross-coupling. Cross-coupling in the H–OH receiver is an order of magnitude less and this calibration step was not required. Figure 1 shows (truncated) PDFB1 polarization profiles for PSR J0437−4715 taken with the MB and H–OH receivers using the observational and processing methods described above. The results are essentially identical within the noise uncertainties, confirming that our calibration procedures are robust.¹⁰

Zoom In Zoom Out Reset image size

While not the primary objective of this work, it was necessary to monitor period variations of the pulsar in order make accurate predictions of the topocentric period for use during the observations and to sum observations in off-line processing. Initially, the timing model of Bailes et al. (2003) was used for this purpose. However, while processing the 2007 July 18 observations, it was realized that the pulsar had suffered a sizable glitch. Glitch parameters resulting from a fit of a model including post-glitch exponential decay (see, e.g., Wang et al. 2000) to pulse time-of-arrival (ToA) data from 2006 April to 2008 November are given in Table 2. ToAs were obtained from average pulse profiles for data segments of 15 minute duration and the program tempo2 (Hobbs et al. 2006) was used to determine the timing parameters. There was only a small relaxation of the pulse frequency toward the extrapolated pre-glitch solution following the glitch, i.e., the Q parameter was small. Significant timing noise remains after the fit and the uncertainties have been multiplied by five (${\sim} \sqrt{(}\chi ^2_r)$) to allow for this.

Table 2. PSR J1141−6545 Glitch Parameters

Parameter	Value
Glitch epoch (MJD)	54277 ± 20
Δνg/ν	(5.890 ± 0.006) × 10−7
$\Delta \dot{\nu }_g/\dot{\nu }$	(5.0 ± 0.9) × 10−3
Q	0.0040 ± 0.0007
τd (d)	495 ± 140
Data span (MJD)	53834–54785
Number of ToAs	254
Rms residual (μs)	220

Download table as:  ASCIITypeset image

Table 3 gives the pulsar timing parameters from a fit to post-glitch data from 2007 December to 2008 November; the ∼160 days following the glitch were omitted to avoid biasing the result by the post-glitch relaxation. Over the fitted data span there is significant timing noise and the fit uncertainties have been multiplied by a factor of 3. Only the pulse frequency ν and $\dot{\nu }$ were fitted for; the position and the binary parameters were held fixed at the values given by Bhat et al. (2008). Parameters are quoted in TDB-compatible units, and the "DD" binary model (Damour & Deruelle 1986) and the Jet Propulsion Laboratory Solar-System ephemeris DE405 (Standish 1998) were used.

Table 3. PSR J1141−6545 Post-glitch Timing Parameters

Parameter	Valuea
R.A. (J2000)	11h41m070140b
Decl. (J2000)	−65°45'191131b
Pulse frequency (ν) (Hz)	2.538723048486(4)
Pulse frequency time derivative ($\dot{\nu }$) (s−2)	−2.77615(8) × 10−14
Epoch (MJD)	54637.00
Dispersion measure (cm−3 pc)	116.080b
Binary period (days)	0.1976509593b
Binary period time derivative	−4.3 × 10−13b
Orbit semimajor axis (s)	1.858922b
Eccentricity	0.171884b
Periastron time (MJD)	51369.8545515b
Longitude of periastron (°)	42.4561b
Longitude time derivative (° yr−1)	5.3096b
Grav. redshift/time dilation (γ) (s)	0.000773b
Data span (MJD)	54435–54785
Rms timing residual (μs)	56

Notes. ^aUncertainties in the last quoted digit are given in parentheses. ^bFrom Bhat et al. (2008).

Download table as:  ASCIITypeset image

Figure 2 shows the observed time variations of the 50% mean pulse width (W₅₀) and mean flux density at 1400 MHz over the 9.3 yr data span. Such dramatic long-term variations in pulse width (more than a factor of 3) are unprecedented. The most straight-forward explanation is that they result from precession of the pulsar spin axis changing our view of the pulsar beam. The variations are complex with a ∼400 day interval around MJD 52800 (2003 June) where the width clearly decreases with time before increasing again. Recent data show that a maximum width was reached around MJD 54150 (2007 February) and now the width is decreasing. There is good agreement between the widths obtained with the different back-end systems, although the AFB widths tend to be a few percent less than those measured with the WBC or PDFB1 at the same or similar time owing to differences in the instrumental impulse response. The time variation in W₁₀ is similar, but with a smaller relative change over the data span.

Zoom In Zoom Out Reset image size

The mean pulsed flux density has also changed dramatically over the data span with a broad peak around MJD 52500 (mid-2002). There is a significant day-to-day variation in the measured flux densities. This is indicated by the large error bars on the AFB averages (which typically cover a week or more) and by the scatter in the WBC/PDFB measurements (which typically are made on a single day), both much larger than the uncertainty in an individual measurement. These short-term variations can be attributed to refractive scintillation for which the expected timescale is about 7 days (Ord et al. 2002). However, the broad rise and decay of the flux density over the nine-year data span is unlikely to be a scintillation effect and we attribute this to the changing aspect of a complex beam pattern sweeping across the Earth as the pulsar spin axis precesses.

Average total-intensity profiles as a function of time are shown in Figure 3 illustrating the profile evolution¹¹. At early times, the profile was dominated by a strong trailing peak with a broad ramp of emission leading up to it, whereas at late times, the profile is more symmetric with an approximately Gaussian shape. Figure 4 shows the profile evolution in the form of grayscale plots with both linear (left) and logarithmic (right) intensity scales. The logarithmic plot clearly illustrates a striking property of the profile evolution, that the width at a very low level (the lowest contour is at 1% of the peak) is remarkably constant, suggesting that the low-level flux-density contours represent the overall beam extent. As will be further described in Section 4, the outer parts of the profile are relatively stable, both in shape and flux density. The profiles in Figures 3 and 4 have therefore been aligned using the midpoint at a constant flux density of 2% of the maximum as a reference point. This is an important assumption as it affects not only the interpretation of the various profile components, but also the parameters derived from any timing analysis. With this assumption, the peak of the profile has moved in a somewhat step-wise fashion from near its trailing edge to near its center. This evolution will be interpreted in terms of the growth and decay of components representing bright regions of the beam in Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the mean pulse polarization profiles for PSR J1141−6545 at two epochs separated by about three years. Plotted P.A.s refer to the band-center frequency as marked on the plots.¹² Dramatic changes are observed over the three-year span in the (total intensity) pulse profile and in the polarization parameters. Once the sign differences are allowed for and ignoring the absolute value of P.A., the 2004 profile in Figure 5 agrees well with the profile given by Hotan et al. (2005a) taken about 2 months earlier. In particular, the P.A. changes rapidly near the profile center, indicating that our line of sight traverses the beam relatively close to the magnetic axis. However, there are clear and significant departures from the P.A. variation expected for the simple RVM.

Zoom In Zoom Out Reset image size

The evolution of the profiles of linearly polarized intensity L and of Stokes V across the pulse over the 4.5 yr span of the polarization data is shown in Figure 6. Both show systematic changes as a function of time. In the case of the L profiles, components at different pulse phases appear to get stronger or weaker as a function of time, whereas the variation of Stokes V is more consistent across the whole profile. It is striking that the variation of V across the pulse changes smoothly from a positive-to-negative sign change at early epochs to a negative-to-positive change at late epochs. Especially at later times, there are clear dips in the linearly polarized intensity at pulse phases near ±0.01. As will be discussed further below, these are attributed to overlapping of approximately orthogonally polarized components.

Zoom In Zoom Out Reset image size

Figure 7 shows the evolution of the P.A. variation across the pulse. The observed P.A. variations are clearly not well described by the RVM and there are clear systematic trends in the P.A. variations as a function of time. While there are some differences between adjacent P.A. profiles which may in part result from residual calibration errors or low-level radio-frequency interference, in general the trends in polarization properties with time are smooth and consistent despite the use of different receivers and back-end systems at different times, again showing that our calibration procedures are effective.

Zoom In Zoom Out Reset image size

Most striking is the clear difference in P.A. evolution between the central and the outer parts of the profile. Figure 8 shows the P.A. variations over the 4.5 yr data span averaged over the inner part (phase −0.0078 to 0.0068) and the outer parts (−0.0253 to −0.0117 and 0.0136–0.0253) of the profile. Variations of the mean P.A. for the inner part of the pulse profile are very significant with an initial decrease followed by a rapid increase which slows at later epochs. For the outer parts of the profile, the mean P.A. at early times shows some fluctuations but overall is consistent with a slow decrease, flattening at later times. Especially at early times, the polarization in the outer parts of the profile is quite weak and hence the fluctuations in P.A. can be attributed to remaining systematic errors and/or contributions from the wings of the stronger emission in the center of the profile.

Zoom In Zoom Out Reset image size

These P.A. changes could be due to variations in RM or to aspect changes resulting from the precessional motion. In Figure 9, we show the RMs separately for the central and outer parts of the profile computed using the iterative method described in Section 2. For the outer parts of the profile, there is no significant time variation. The weighted mean RM value is −93.4 ± 2.2 rad m⁻² with a reduced χ² of 3.1; the quoted uncertainty is the weighted rms deviation multiplied by $\sqrt{(}\chi ^2_r)$.¹³ In contrast, for the central part, there is a clear systematic variation of measured RM with time. We believe that this is not a true variation in the interstellar RM—it would be much larger than RM variations in other pulsars that are attributed to a changing path through the interstellar medium, e.g., for the Vela pulsar (Hamilton et al. 1985). Furthermore, it is inconsistent with the changes in the central P.A. shown in Figure 8. Rather, we believe the observed changes in both P.A. and apparent RM for the central part of the PSR J1141−6545 profile are due to variations in the relative amplitudes of observed profile components with different P.A.s and spectral indices resulting from the changing aspect of the line of sight relative to the beam axis as the pulsar spin axis precesses. This is somewhat analogous to the apparent RM variations (as a function of pulse phase) seen in PSR B2016+28 by Ramachandran et al. (2004) which are attributed to overlapping non-orthogonal components.

Zoom In Zoom Out Reset image size

We adopt P.A.s in the outer zones as representing the underlying magnetic-field structure. Consequently, these P.A. variations should be consistent with the RVM. The two sides are consistently at approximately the same P.A., as expected in the RVM. Also, in some cases, e.g., the 2007 August 4 profile shown in Figure 5, there is evidence for the P.A. variations expected for the RVM at the inner edges of the outer P.A.s. Since there is no significant RM variation for these parts, we interpret the observed slow decrease in outer-zone P.A.s shown in the right panel of Figure 8 as intrinsic and resulting from precession of the pulsar spin axis.

By modeling the observed P.A. variations in terms of precessional motion of the pulsar's spin axis, we can infer additional properties of the system, for example, the spin–orbit misalignment angle and the precessional phase. As a first step in the modeling, we fit the observed total-intensity profiles (Figure 3) with Gaussian components. A total of six components are fitted, although for later epochs only four are required. We have used fixed central phases for the components; only the amplitude and width of each component is allowed to vary as a function of epoch. The component phases, given in Table 4, were chosen using an iterative process by fitting the first and last observed profiles. Fitted components and their sum are shown in Figure 10 for a selected set of profiles which cover most of the observed data span, showing that an excellent fit to the observed profiles is obtained. This retroactively justifies our decision to fix the component phases. Figure 11 gives the widths and mean flux densities of the fitted components. The mean flux densities are scaled so that their sum equals the smoothed flux-density variation shown in Figure 2. For components 3 and 5, the fitted amplitudes dropped to zero after MJDs 52100 and 52500, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is striking that the amplitude of each component varies smoothly from one epoch to the next and that the widths of the components are relatively stable. This suggests that these components represent real physical zones of emission that are fixed in longitude on the star. We note that most of the observed flux-density variation comes from the central and trailing components. Although the leading component 1 is relatively weak at early epochs and dominates the profile at late epochs, its absolute amplitude remains approximately constant across the whole data span.

As discussed above, the observed polarization P.A. and L variations (Figures 7 and 6) strongly suggest that these components may be grouped into two overlapping and quasi-independent emission zones, one dominant in the wings of the observed pulse and the other dominant in the central regions. We note that similar non-RVM variations (not accountable for by overlapping orthogonal polarization modes) are seen in other short-period pulsars, for example, PSR B1913+16 (Blaskiewicz et al. 1991) and PSR B1534+12 (Arzoumanian et al. 1996; Stairs et al. 2004). The outer components (1 and 6) may be identified with the outer zones in the P.A. and L plots. As discussed above, the P.A.s in these zones together are well fitted by a single RVM. Note that, as discussed in Sections 1 and 2, for modeling purposes, we have adopted a right-hand coordinate system in which the sign of P.A. is reversed compared to the astronomical convention used in Section 3. We also choose ϕ₀ = 0.0 and hold this fixed for all epochs.

Most previous discussions of relativistic spin precession in this and other pulsars (e.g., Weisberg et al. 1989; Kramer 1998; Stairs et al. 2004; Hotan et al. 2005a) have concentrated on measurement and interpretation of the rate of change of the impact parameter. For PSR J1141−6545, the impact parameter is not very well determined because of the presence of the evidently independent emission in the central part of the profile. However, since we have measured absolute P.A.s, we have another, potentially more sensitive, observable: the central P.A. of the RVM, ψ₀ (cf. Stairs et al. 2004; Kramer & Wex 2009). In the RVM, this is the direction of the spin axis of the star projected on the sky plane (see Figure 12). This of course varies as the pulsar precesses:

$$\begin{equation} \psi _0 = \Omega _{\rm asc} + \eta, \end{equation} \tag{ 4 }$$

where Ω_asc is the longitude of the ascending node and η, the precessional longitude, is the angle of the pulsar spin axis relative to the ascending node projected on the sky plane. In this equation, ψ₀ is the observed P.A. corrected to infinite frequency using the measured RM of −93.4 rad m⁻². The polar angle λ varies with the precession according to

$$\begin{equation} \cos \lambda = \cos \delta \cos i - \sin \delta \sin i \cos \Phi, \end{equation} \tag{ 5 }$$

$$\begin{equation} \cos \eta \sin \lambda = \sin \delta \sin \Phi, \end{equation} \tag{ 6 }$$

where δ is the spin–orbit misalignment angle, Φ is the precessional angle measured in the orbit plane from $\bf -j$, and Φ₀ is the value of Φ at t = t₀ (DT92) and

$$\begin{equation} \Phi = \Omega _p (t-t_0) + \Phi _0. \end{equation} \tag{ 7 }$$

From Equation (4), since Ω_asc is effectively constant,¹⁴ the precessional variation in ψ₀ is determined by the variation of η. From Equation (5), we can determine cos λ as a function of time and the fixed angles δ, i, and Φ₀. We take i = 730 from Bhat et al. (2008), take t = t₀ at MJD 53000.0, near the center of the data span and the beginning of our polarization monitoring, and take the GR prediction for the spin precession rate Ω_p, approximately 136 yr⁻¹. Then, using Equation (6) (where there is no sign ambiguity in sin λ since 0 ⩽ λ < π), we obtain cos η. From the same spherical triangle used to derive Equation (6), we also have

$$\begin{equation} \cos \delta = \cos \lambda \cos i - \sin i \sin \lambda \sin \eta \end{equation} \tag{ 8 }$$

giving us sin η and hence η as a function of time (and the fixed but unknown angles Ω_asc, δ, and Φ₀) without further ambiguity (cf. Kramer & Wex 2009).

Consequently, the set of three parameters, δ, Φ₀, Δψ₀, (where Δψ₀ absorbs the contribution to ψ₀ by Ω_asc and any other intrinsic constant offset to the P.A. values) describes the variation of ψ₀ as a function of time. At the same time, the P.A. swing at each epoch is ideally described by the RVM, which uses two further parameters: the magnetic inclination α and the impact parameter β. (We have held the central longitude of the RVM, ϕ₀, at zero as discussed above.) However, since

$$\begin{equation} \lambda = \pi - \zeta = \pi - \alpha - \beta, \end{equation} \tag{ 9 }$$

the angle β can be computed from Equation (5) for a given α. As a result, only four parameters should describe the behavior of the P.A. variations at all epochs.

Zoom In Zoom Out Reset image size

As described above, the central range of the P.A. swing shows clear deviations from an RVM. However, in a blind search for the best-fit solution, it nevertheless adds valuable information as it clearly indicates (by its slope) that β < 0 (in the DT92 convention). It can also be noted that the slope steepens and then flattens again toward the end of our data set (see Figure 7). In order to use the geometric information provided by the central part but to still tie the results to the outer wings, we decided to decrease the relative weight of the central part in the fitting process (by multiplying its errors by a factor of 5¹⁵). We allowed a P.A. offset between the central part relative to the RVM determined by the outer wings, and this was fitted for independently for each epoch. We take longitudes in the range of |ϕ| ⩽ 3° (0.00833 in phase) to define the central region, whereas P.A.s in the transition zone 3° < |ϕ| < 5° (0.00833–0.0139 in phase) were ignored in the fit.

As a first step of the fitting procedure, we constrain the allowable values of δ and Φ₀ using just the outer-zone P.A.s. For this purpose, we compute the mean P.A. of the outer wings and fit its behavior as a function of time (see Figure 13). We use a grid search in the δ–Φ₀ plane, where at each grid point we compute the χ² for the fit of the model values (based on Equations (4)–(8)) to the computed values of ψ₀ (Figure 13). Figure 14 is a plot of χ² in the δ–Φ₀ plane, showing two clear and well-defined solutions. For this plot, we have taken i = 730; the plot for 180 − i is mirror symmetric.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The solution with Φ₀ ∼ 0 corresponds to λ>π/2 while the other solution implies λ < π/2. We can use the information provided by the P.A. swing to identify the appropriate solution in the following way. From the P.A. swing (in the DT92 convention), the slope is negative and hence β < 0. Also, the fact that the P.A. becomes flat at the outer wings suggests an outer line of sight (that is, the line of sight is on the equatorial side of the magnetic axis when ϕ = ϕ₀; see, e.g., Lorimer & Kramer 2005) rather than an inner one. The combination of both is only possible in a solution with ζ = α + β>π/2, corresponding to a solution where λ = π − ζ < π/2. This uniquely identifies the region around δ ∼ 90° and Φ₀ ∼ 175° as the correct solution.

The best solution from this fitting procedure corresponds to λ being very close to zero. We also know that β is close to zero from the rapid swing of P.A. in the central region and the flat P.A.s in the outer zones which are only a few degrees from the pulse center. Therefore, from Equation (9), we have that α ∼ π. However, the total-intensity profiles and P.A. swing show that the magnetic axis cannot be exactly aligned (or counteraligned) with the rotation axis. The fact that the sign of the central P.A. swing remains the same at all epochs also shows that β does not change sign as a result of the precessional motion. This suggests that the fit to the P.A. data as a whole will lie somewhat to the right of the best-fit point shown in Figure 14. This idea is confirmed by a blind fit to all P.A. data using only the four parameters δ, Φ₀, Δψ, and α. The search is again performed on a δ–Φ₀ grid, where at each grid point we execute a simplex algorithm to minimize for the other two parameters. The results (collapsed in Δψ and α space) are shown in Figure 15. We show in grayscale the previous results of the P.A.-offset fit and overlay contour levels for the global RVM fit to the P.A. data. The model variation of ψ₀, shown in Figure 13, has a reduced χ² of 7.4 with 920 degrees of freedom. Assumed and derived parameters for the system are given in Table 5. Figure 16 shows the joint probability density function for the spin–orbit misalignment angle δ derived from the fits to the P.A. variations and the global RVM fit (Figure 15).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Model Parameters from the RVM Fit to the Observed Position-angle Variations

Parameter	Valuea
Assumed parameters
Orbit inclination i	73°
Precession rate Ωp	136 yr−1
RVM reference phase ϕ0	0.0
Reference time for precession	MJD 53000.0
Derived parameters
Spin–orbit misalignment angle δ	93° (−9°, +16°)
Precession angle Φ0	1755 (−28, +18)
Offset in P.A. Δψ0	−1187 ± 02
Magnetic inclination α	160° (−16°, +8°)

Note. ^aMost probable value and 68% confidence limits.

Download table as:  ASCIITypeset image

Given the geometry of the system we can obtain the model time variation of impact parameter β using Equations (5)–(7) and (9). The result, given in Figure 17, shows that the observer's line of sight has approached the magnetic axis direction through most of our observed data span, reaching its closest approach around MJD 54000, and is now receding again. This reversal in the variation of β is supported by the fact that the pulse width reached a maximum at about MJD 54200 and is now decreasing (Figure 2). Note that we have obtained these results by fitting only the polarization data, i.e., we did not make use of the total-intensity information. Although the overall variation of pulse width is consistent with the precessional model, the detailed variations are more complex. This shows that the assumption of a uniform circular beam is not appropriate.

Zoom In Zoom Out Reset image size

Given the variation of β over the data span (Figure 17) and the time variation of component shapes and amplitudes (Figure 11), we can compute the two-dimensional intensity profile of the emission beam over the traversed region. For each epoch, we compute the traverse of our line of sight across the polar region and accumulate the pulse intensity in a two-dimensional grid centered on the magnetic axis. The resulting beam pattern is shown in Figure 18. Because of the reversal in the time derivative of β when it was close to zero, we only see one half of the polar-cap region. Despite this we can clearly see that the beam is quite asymmetric with no evidence for a core-cone or ring structure that is symmetric about the magnetic pole. The partially filled beam can be described as "patchy," albeit with just one major patch in the region scanned so far. Although this is the first two-dimensional map of an emission beam to clearly show such patchy structure, there is good evidence that the emission beam from most pulsars is best described in this way (Lyne & Manchester 1988; Han & Manchester 2001).

Zoom In Zoom Out Reset image size

Although the observed pulse width is about average (the median pulse width for all pulsars is about 10°), the small value of 180° − α ∼ 20° implies that the intrinsic beamwidth is small. Figure 18 shows that the emitting region fits within a circle of radius about 4° centered on the magnetic axis. Excluding millisecond pulsars, observed pulse widths are generally consistent with the relation

$$\begin{equation} \rho = 6\mbox{.\!\!^\circ $}5 P^{-0.5} \end{equation} \tag{ 10 }$$

(Gil et al. 1993; Kramer et al. 1994). For PSR J1141−6545 the predicted value of ρ is about 103, much larger than the emitting zone traversed so far.

PSR J1141−6545 was not detected in the Parkes 70 cm survey (Manchester et al. 1996) although, even at the present relatively low flux-density levels, a detection with signal/noise of the order of 50 would have been expected. Observations within half a beamwidth of the pulsar position were made 1993 July 14 (MJD 49182). Figure 17 shows that the impact parameter β at that time was about −8°. This non-detection therefore suggests that the beam half-width in latitude is ≲8°, although this must be qualified because of the patchy beam structure. There are good arguments (e.g., Narayan & Vivekanand 1983; Manchester 1996) that beams in young pulsars are elongated in the latitude direction and hence more fan-like than circular. An elongated beam has also been suggested for PSR B1913+16 (Weisberg & Taylor 2002; Clifton & Weisberg 2008). Observations over the next decade or two will establish whether or not this is the case for PSR J1141−6545 as β returns to large (negative) values.

There is good and increasing evidence (e.g., Johnston et al. 2005; Ng & Romani 2008) that the velocity and spin vectors of young pulsars are nearly aligned. This suggests that the kick from the supernova (SN) explosion, the dominant factor in determining the pulsar space motion, may also determine the spin direction (see also Spruit & Phinney 1998). Since PSR J1141−6545 is a young pulsar, its spin direction is therefore likely to be aligned with the SN kick. While precession will have caused the azimuthal angle of the pulsar spin relative to the (post-SN) orbital angular momentum to wrap many times over the lifetime of the pulsar, the polar angle (δ in this paper) should still reflect the direction between the kick and the post-SN orbit normal.

We have explored if this expectation can set any interesting constraints on the kick imparted to the forming neutron star, using the binary orbit and kick description in Kalogera (1996) and Wex et al. (2000). In terms of their variables, we find

$$\begin{equation} \cos \delta = \frac{ V_{\rm kz} (2 V_{\rm ky} +V_0)}{|V_{\rm k}| \big[V_{\rm kz}^2+(V_{\rm ky}+V_0)^2\big]^{1/2}}. \end{equation} \tag{ 11 }$$

We follow the procedure set out in Thorsett et al. (2005) and the Bayesian version in Stairs et al. (2006), which follows the principles of Wex et al. (2000) but which puts priors on the system velocity, the angle θ between the pre- and post-SN orbit normals¹⁶ and the orientation Ω of the binary system on the sky. Each sampling of the priors determines the coefficients of a quadratic equation in m_2i, the mass of the pre-SN star (Equation (8) of Thorsett et al. 2005), which may or may not have a solution that falls within acceptable ranges for the pre-SN mass and the size of the (assumed circular) pre-SN orbit. As this is a young system, we did not evolve the system's motion back in time through the Galaxy to identify plausible birth sites (Wex et al. 2000), nor did we evolve the orbital size and eccentricity to account for gravitational radiation losses (Peters 1964), but simply assumed the current velocity and orbital parameters reflect the birth properties. Note that the results are therefore independent of the assumed distance; we assumed a distance of 3.7 kpc (Ord et al. 2002). We put uniform priors on θ (sampling from 0° to 180° and testing both positive and negative values in each trial) and Ω (0°–360°). We sampled the companion mass over (1.02 ± 0.01) M_☉ and the total system mass over (2.28911 ± 0.00026) M_☉, based on the orbital parameters reported in Bhat et al. (2008). We subtracted these to get the pulsar mass for the given trial, and used the mass function computed from the orbit presented in Bhat et al. (2008), 0.176550265 M_☉, to derive sin i. We considered only the cases where cos i>0 in order to match the modeling in the previous section.

While there is a scintillation velocity measurement for this system (Ord et al. 2002), which could in principle provide some constraints on the angle between the proper motion and Ω, similar measurements have been shown to be quite unreliable in the Double Pulsar, probably because of the effects of anisotropies in the interstellar medium (Ransom et al. 2004; Coles et al. 2005; Kramer et al. 2006) and hence we do not make use of this information here. Instead, we investigate Maxwellian velocity distributions with one-dimensional dispersions of 50 km s⁻¹ and 100 km s⁻¹ (these dispersions are of order the transverse velocity derived from the scintillation measurements), combined with a second Maxwellian with one-dimensional dispersion of 12 km s⁻¹ representing the pre-SN peculiar velocity. We then derive Bayesian probability distributions for the pre-SN and kick parameters as described in Stairs et al. (2006), considering pre-SN masses in the range 1.4–10.0 M_☉ to allow for the binding energy of the neutron star (Lattimer & Prakash 2001). For each acceptable solution of the quadratic equation, the relevant set of parameters is assigned a uniform likelihood (as in Stairs et al. 2006) for a "natural" weighting. We then make use of the probability distribution for the post-SN pulsar spin–orbit misalignment angle δ obtained from the P.A. fitting (Figure 16) to give a constrained weighting on the range of acceptable pre-SN parameters and kick velocity. The derived probability density functions for the pre-SN stellar mass, the kick velocity, and the spin–orbit misalignment angle after constraining with the allowed misalignment-angle distribution (Figure 16) are shown in Figure 19. Table 6 gives the median values of the parameters and the limits at 68% and 95% confidence. We note that the natural weighting tends to disfavor spin–orbit misalignment angles around 90°, so the constrained distribution is biased toward slightly larger angles compared to Figure 16.

Zoom In Zoom Out Reset image size

Not unexpectedly, the assumed distribution of post-SN system velocities has a significant effect on the likely kick velocity, with larger values requiring a larger kick. However, the most probable kick velocity is relatively small in both cases, about 100 km s⁻¹ for a dispersion of 50 km s⁻¹ and about 180 km s⁻¹ for a dispersion of 100 km s⁻¹. Median velocities (Table 6) are somewhat higher since the distributions have a high-velocity tail. Nevertheless, they are still small compared with mean or median of the inferred three-dimensional pulsar velocity distribution (Lyne & Lorimer 1994; Hobbs et al. 2005; Faucher-Giguère & Kaspi 2006).

Relatively low progenitor masses are strongly favored with median values of about 2 M_☉. The pulsar mass implied by the timing fits for companion mass and total mass (Bhat et al. 2008) is just 1.27 M_☉, lower than average (Thorsett & Chakrabarty 1999), but very similar to the mass of PSR J0737−3039B, the second-born and slow pulsar in the Double Pulsar system. This system also has a relatively low implied kick velocity from the B pulsar formation (e.g., Piran & Shaviv 2005; Willems et al. 2006; Stairs et al. 2006; but see also Kalogera et al. 2008). These low masses and kick velocities may imply formation of the neutron star in an electron-capture collapse of the core of an ONeMg white dwarf (Podsiadlowski et al. 2005) but other evolutionary histories are also possible (e.g., Dewi & Pols 2003).

The inferred spin–orbit misalignment angle δ in the PSR J1141−6545 system is large, with a most probable value of between 100° and 110°. That δ is not small (in contrast to pulsar A of the Double Pulsar system; Stairs et al. 2006; Ferdman et al. 2008), is consistent with the major changes that we see in the pulse profile as a result of the precession and in fact results from the fitting of the precessional model to these changes.

The dramatic long-term changes in pulse shape and amplitude observed for PSR J1141−6545 (Figures 2 and 3) are unprecedented in pulsar astronomy. Pulse profile changes are observed in mode-changing pulsars (e.g., Wang et al. 2007) and in radio emission from magnetars (e.g., Camilo et al. 2007) but these are of a quite different character to the variations seen in PSR J1141−6545 and are certainly caused by fluctuations in the emitted pulsar beam. Slow and systematic variations similar to those observed in PSR J1141−6545 are seen in the pulse profiles of PSR B1913+16 and PSR B1534+12, but these are of much smaller amplitude. All three of these pulsars are in close binary orbits with massive companions. Consequently, the observed profile variations most probably result from precession of the pulsar spin axis resulting from spin–orbit coupling changing our view of the pulsar beam.

Precession of the pulsar's spin axis has resulted in large and systematic changes in both the amplitude and the shape of the observed pulse profile as our line of sight scanned across the polar region. There are clearly large variations in emissivity across the nominal beam area, consistent with the idea of a partially filled or patchy beam. The large variations in the observed mean pulse profile raise the issue of profile alignment when comparing different epochs. Figure 4 shows that the pulse width is relatively stable at the lowest contour levels. We have chosen to align the profiles according to the midpoint of a constant flux-density level near their extreme wings. This choice is motivated by the fact that the outer components are relatively stable, both in flux density and shape (Figure 11). With this choice, alignment of the polarization features over the data span is also very stable, giving added support to our method of defining the absolute pulse phase relative to the star. We note that this choice also defines a reference phase for pulse timing. It is interesting to note that, for PSR B1913+16 also, the pulse width defined by the lower beam contours is more constant than, say, the pulse width at the 50% level (Weisberg & Taylor 2005).

Major changes are also observed in the pulse polarization parameters over the 4.5 yr in which the polarization has been monitored. There has been a steady evolution of the circular polarization with the sign of the sense reversal near the profile center changing from positive-to-negative to negative-to-positive over the data span. Variations in the linear polarization are complex with quite different behavior in the central and outer zones of the pulse profile (Figure 7). Our choice of taking the outer-zone P.A.s to represent the underlying magnetic-field structure is supported by a number of factors. First, the outer-zone P.A.s are well fitted by the RVM whereas this is less true of the central-zone P.A.s. More importantly, there is a large and variable offset between the central P.A.s and the outer-zone P.A.s. The time variation of the mean central-zone P.A. is rapid and inconsistent with precessional motion. Furthermore, there is an apparent time variation of RM for the central zone which is not present for the outer zones. We attribute the complex behavior of the central-zone polarization to the overlapping of independent components that have different P.A.s and different spectral indices. It may be possible to model the polarization variations based on the variations of the total-intensity components shown in Figure 11 but this is beyond the scope of the present work. We note that similar but much less dramatic non-RVM variations are observed in the central parts of the profiles for PSR B1913+16 and PSR B1534+12.

Fitting of a precessional model to the observed P.A. variations leads to the conclusion that the spin–orbit misalignment angle δ is very large in this system, with the pulsar spin axis nearly orthogonal to the orbital angular momentum vector. This is consistent with the large profile shape and polarization changes as the spin axis precesses. Unfortunately, because of the complex and poorly understood emission physics, it is not possible to turn the argument around and use these observations as a test of the precessional predictions of general relativity. It is interesting to note though that a large spin–orbit misalignment angle is quite possible despite the low progenitor mass and modest kick velocity. This result is not in conflict with the conclusion of Hotan et al. (2005a), that the spin–orbit misalignment angle is probably less than 30°, since that was on the proviso that the pulsar was not at a precessional phase when the angle between the observer's line of sight and the spin axis (λ in our terminology) was changing rapidly. In our solution, λ is close to zero and is changing rapidly.

Our observational data span encompassed the time when the precessional longitude passed through 180°. This means that the impact parameter β reached an extremum during our data span. In fact, |β| reached a minimum value very close to zero at about MJD 54000 (Figure 17). This solution is consistent with the observed maximum in the pulse width at around the same time and implies that we will only ever traverse one side of the emission beam. The slope of the β time variation has now reversed and so we are now retracing our earlier path across the beam. We therefore predict that over the coming decade we will see a time-reversed "replay" of the recent profile amplitude and shape evolution.

We thank our colleagues at the Parkes Multibeam Pulsar Survey team for assistance with the observations reported in this paper. We also thank Steve Thorsett and Rachel Dewey for their contributions to the modeling of the effects of the pulsar kick velocity (Section 6). Pulsar research at UBC is supported by an NSERC Discovery Grant and a CFI New Opportunities Grant to I.H.S. and M. Berciu. C.A.M. held an NSERC USRA and I.H.S. an NSERC UFA for part of this work. I.H.S. further acknowledges support from the ATNF Distinguished Visitors Program and the Swinburne University Visiting Distinguished Researcher Scheme. The Parkes radio telescope is part of the Australia Telescope which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO.