A HIGH BRAKING INDEX FOR A PULSAR

Pulsars emit light by extracting energy from their rotational kinetic stores. As such, their spin-down is expected to follow the form (Manchester et al. 1985)

$$\begin{eqnarray}&&\dot{\nu }=-K{\nu }^{n},\end{eqnarray} \tag{ 1 }$$

where ν is the spin frequency of the pulsar, $\dot{\nu }$ the frequency derivative, and K a constant of proportionality related to the pulsar's moment of inertia and magnetic field structure (Gunn & Ostriker 1969). Here, n is the braking index. In the standard pulsar model of an unchanging magnetic dipole in a vacuum, electrodynamics predicts a value of three (e.g., Manchester & Taylor 1977). In more realistic models of a pulsar and its magnetosphere, the braking index is predicted to always lie between 1.8 and 3 (Melatos 1997). Values less than this can be obtained by relaxing the various assumptions of the model—e.g., allowing magnetic field evolution (Blandford & Romani 1988), momentum loss due to a particle wind (Harding et al. 1999), or a varying angle between the spin and magnetic poles (Lyne et al. 2013).

Taking the time derivative of Equation (1) gives us the following fundamental equation that contains only the braking index and observable quantities:

$$\begin{eqnarray}&&n=\displaystyle \frac{\nu \ddot{\nu }}{{\dot{\nu }}^{2}},\end{eqnarray} \tag{ 2 }$$

where $\ddot{\nu }$ is the second derivative of the spin frequency. A measurement of n can be made only for the youngest pulsars for which $\ddot{\nu }$ is large enough to detect on human timescales. As such, only 8 pulsars of the ∼2400 known have measured braking indices, with values ranging from 0.9 ± 0.2 to 2.839 ± 0.001 (Lyne et al. 2015 and references therein).

PSR J1640–4631 was discovered as a pulsating X-ray source in a NuSTAR survey of the Norma region of the Galactic plane (Lemiere et al. 2009; Gotthelf et al. 2014). The pulsar is located in the center of the supernova remnant G338.1−0.0 and powers the pulsar wind nebula (PWN) HESS J1640−465, first detected in very high energy gamma-rays and thought to be the most luminous TeV source in our Galaxy (Gotthelf et al. 2014). We undertook X-ray timing observations of PSR J1640–4631 starting shortly after discovery with the aim of measuring its braking index.

All X-ray observations presented in this work were taken using NuSTAR, which consists of two co-aligned X-ray telescopes sensitive to photons with energies from 3 to 79 keV (Harrison et al. 2013). NuSTAR observations of PSR J1640–4631 were typically 20–50 ks, and the observation cadence can be seen in Figure 1. Level 1 data products were obtained from HEASARC and reduced using nupipeline v0.4.4. Photons from a circular region having a 30'' radius centered on the source were extracted. To maximize the signal-to-noise ratio of the pulse, we used only photons with energies in the 3.0–55 keV range.

Zoom In Zoom Out Reset image size

Photon arrival times were corrected to the solar system barycenter using the Chandra position of PSR J1640–4631, R.A. = 16^h40^m4352, decl. = −46°31'354 (Lemiere et al. 2009) using barycorr from HEASOFT v6.17 and v052 of the NuSTAR clock file.

Photon arrival times were used to derive an average pulse time-of-arrival (TOA) for each observation. The TOAs were extracted using a maximum likelihood (ML) method. The ML method compares a continuous model template of the pulse profile to the computed phases of the photon arrival times from an observation (Livingstone et al. 2009). To create the template, all observations were folded into a high signal-to-noise profile. This high signal-to-noise profile was then fitted to a Fourier model using the first two harmonics. Two harmonics were chosen to optimally describe the pulse shape, as determined by use of the H-test (de Jager et al. 1989). We verified that TOAs extracted using a cross-correlation method give consistent results. NuSTAR's absolute timing calibration is accurate to ±3 ms (Madsen et al. 2015), smaller than our measurement uncertainties.

The TOAs were fitted to a standard timing model in which the phase as a function of time t is described by a Taylor expansion:

$$\begin{eqnarray}\phi (t) & = & {\phi }_{0}+{\nu }_{0}(t-{t}_{0})+\displaystyle \frac{1}{2}\dot{{\nu }_{0}}{(t-{t}_{0})}^{2}\\ & & +\displaystyle \frac{1}{6}\ddot{{\nu }_{0}}{(t-{t}_{0})}^{3}+\ldots \end{eqnarray} \tag{ 3 }$$

This was done using the TEMPO2 (Hobbs et al. 2006) pulsar timing software package. In Table 1, we present a fully phase-coherent timing solution for PSR J1640–4631 over the NuSTAR observation campaign. This is the only solution that provides a statistically acceptable fit, i.e., we have verified there are no pulse counting ambiguities. The residuals, the difference between our timing model and the observed pulse phases, can be seen in Figure 1. In the top panel, we show these residuals accounting only for ν and $\dot{\nu }$. The bottom panel shows the residuals for the full timing solution, accounting for $\ddot{\nu }$. Fitting for an extra frequency derivative, $\dddot{\nu }$ does not significantly improve the fit with the F-test indicating a 52% probability of the improvement of χ² being due to chance. We measure $\ddot{\nu }=(3.38\pm 0.03)\times {10}^{-22}$ s⁻³ corresponding to a braking index of n = 3.15 ± 0.03, where the uncertainty represents the 68% confidence interval.

Table 1.  Phase-coherent Timing Parameters for PSR J1640–4631

Dates (MJD)	56463.0–57298.8
Dates	2013 Sep 29–2015 Oct 3
Epoch (MJD)	56741.00000
ν (s−1)	4.843 410 287 0(5)
$\dot{\nu }\$ (s−2)	−2.280 830(4) × 10−11
$\ddot{\nu }\$ (s−3)	3.38(3) × 10−22
$| \dddot{\nu }| \$(s−4)	<1.4 × 10−30
rms residual (ms)	6.17
rms residual (phase)	0.030
${\chi }_{\nu }^{2}$/dof	0.98/46
Braking index, n	3.15(3)

Note. Figures in parentheses are the nominal 1σ tempo2 uncertainties in the least-significant digits quoted. Upper limits are quoted at the 2σ level. The source position was held fixed at the Chandra position.

Download table as:  ASCIITypeset image

This measured braking index is 5σ higher than that expected in the standard magnetic dipole scenario. In Figure 2, we show all braking indices where the long-term electrodynamic spin-down is believed to be dominant; note how PSR J1640–4631 is the only measurement greater than the canonical n = 3 magnetic dipole line.

Zoom In Zoom Out Reset image size

2.1. Timing Noise Simulations

A possible way to explain such a large measured braking index is contamination from timing noise. Timing noise refers to unexplained low-frequency modulations found in the timing residuals of many, particularly young, pulsars (Arzoumanian et al. 1994). Timing noise in radio pulsars has been observed to be spectrally "red," that is, having most power at low frequencies (Arzoumanian et al. 1994). As such, it can contaminate measurements of $\ddot{\nu }$, and thus n (Hobbs et al. 2010). The power spectral density (PSD) of timing noise is often modeled as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{TN}}(f)=A{\left(1+\displaystyle \frac{{f}^{2}}{{f}_{c}^{2}}\right)}^{-q/2},\end{eqnarray} \tag{ 4 }$$

where A is the spectral density amplitude, f_c the corner frequency, and q the power-law index (Lasky et al. 2015). A q of 0 represents a white noise power spectrum, whereas indices of 2, 4, and 6 represent random walks in pulse arrival phase, the pulse frequency, and frequency derivative, respectively.

To quantify the probability of timing noise biasing our measurement of the braking index, we conducted a series of simulations that aimed to determine whether any reasonable form of red noise artificially results in a measurement of n > 3, given the observed (i.e., white) noise properties of the resulting timing residuals. As timing noise in real pulsars can have PSDs with many indices, we created 10⁵ realizations of red noise and injected them into simulated pulsar TOAs, with TOA uncertainties, ν, and $\dot{\nu }$ identical to that of the pulsar, but with a braking index of n = 3.

To do this, we used the simRedNoise plug-in of TEMPO2. We simulated parameters on a grid, drawing 10 iterations from each set of parameters. We simulated A ranging from 10⁻²⁰ to 10⁻¹⁸ s² yr⁻¹ with 20 log-spaced steps, q between 0 and 6 with 25 linear steps, and f_c from 10⁻³ to 10^0.5 yr⁻¹ in 20 log-space steps. The upper bound on A was chosen to ensure phase connection was possible, as larger values of A precluded phase connection more that 50% of the time, and hence are ruled out by our observations.

After the noise injection, the new TOAs were fitted to the full timing model to measure n, allowing ν, $\dot{\nu }$, and $\ddot{\nu }$ to vary. We considered any iterations where the χ² value indicated that a probability of less than 1% was ruled out by our measured residuals, as this would indicate unaccounted for noise. After this, only 0.01% of all the simulations could produce a braking index greater than three with the measured significance. The parameter regime that gave this highest probability of artificially producing a high braking index was when f_c was of the order the observing length with the highest values of A. There is only a weak dependence on q, with larger values producing more false positives.

Another way to check for the possible contamination of our measured value of n by timing noise is by considering the third frequency derivative. In our data, the third frequency derivative is consistent with zero at the 1σ level. We ran another suite of simulations within the parameter space described above wherein we fitted up to a third frequency derivative. We note that there is a covariance, with the second and third frequency derivatives being anti-correlated. Only in 0.008% of simulated sets of residuals can we reproduce a braking index significantly greater than three without having a detectable third frequency derivative.

These simulations indicate that only a very low level of timing noise can be present in our data, and the measured braking index of n = 3.15 ± 0.03 is highly unlikely to be due to timing noise. We note that the assumed value of n = 3 in our simulations is conservative, since when assuming n < 3, it is even less likely for timing noise to result in a measured n > 3.

2.2. Parkes Observations

In order to search for radio pulsations from PSR J1640–4631, we undertook observations with the 64 m Parkes Telescope. Observations were performed in two sessions totaling 14.96 hr at the position of PSR J1640–4631 (Lemiere et al. 2009). Data were taken by observing with the central beam of the 21 cm multi-beam receiver, using the BPSR pulsar backend. These observations were taken on 2014 April 18 and 2015 April 27 (MJDs 56765 and 56775, respectively) at center frequency 1382 MHz over 400 MHz of bandwidth, divided into 1024 channels. The data from each channel were detected and the two polarizations summed to form a time series with 64 μs samples.

We searched these data using the known pulsar spin frequency and frequency derivative from the phase-coherent timing analysis found in this article. We searched over 4704 dispersion measures between 0 and 1600 pc cm⁻³. No signal was found in these data, and so we quote an upper limit to the pulsar flux at this frequency, using the radiometer equation for the rms noise from the observing system:

$$\begin{eqnarray}&&{\sigma }_{{\rm{rms}}}=\displaystyle \frac{{T}_{{\rm{sys}}}}{G\sqrt{{n}_{p}{t}_{{\rm{obs}}}{\rm{\Delta }}f}},\end{eqnarray} \tag{ 5 }$$

where T_sys is the system temperature in Kelvin, G is the receiver gain in K Jy⁻¹, n_p is the number of polarizations, t_obs is the total integration time in seconds, and Δf is the observing bandwidth in Hz. From this, we calculate a 3σ upper flux limit of 0.018 mJy at 1.4 GHz. This upper limit assumes a 50% duty cycle, with the upper limit scaling as $\sqrt{\mathrm{DC}/(1-\mathrm{DC})}$, where DC is the duty cycle. This flux limit is low, but not unusually so, especially when one considers that the estimated distance to the source is ∼12 kpc.

A possible contaminant to the measured braking index is a long-term recovery from an unseen glitch prior to our monitoring. If such a glitch occurred, a typical exponential recovery would, in general, cause an artificially high value for $\ddot{\nu }$ to be measured (e.g., Johnston & Galloway 1999; Hobbs et al. 2010). Indeed, for pulsars with τ_c < 10⁵ years, there is a clear preference for positive $\ddot{\nu }$, compared to older pulsars that are equally likely to have positive or negative $\ddot{\nu }$ (e.g., Hobbs et al. 2010). If we assume an exponentially recovering glitch contaminating a constant braking index, we can, from our upper limit of $\dddot{\nu }$, place a lower limit on the decay timescale for such a glitch to be $\tau =250\sqrt[3]{{\rm{\Delta }}{\nu }_{d}/{10}^{-8}\;\mathrm{Hz}}$ days, where ν_d is the size of the unseen decaying glitch. In Yu et al. (2013), of the 107 glitches detected, 27 had a detected exponential recovery. Of those, only 3 were longer than 200 days. We thus cannot rule out contamination due to a prior unseen glitch recovery as an explanation for the high measured braking index, although the recovery timescale would have to be longer than most yet observed. Ultimately, this hypothesis can be tested by continued monitoring.

Measurements of braking indices could also be contaminated by uncertainties in the pulsar's position, or proper motion (Bisnovatyi-Kogan & Postnov 1993). The position of PSR J1640–4631 is well determined by Chandra to a 3σ error radius of 06 (Lemiere et al. 2009). This positional error at the ecliptic latitude of PSR J1640–4631 would add a 1.0 ms rms signal to our timing residuals, far smaller than our measurement uncertainties. For PSR J1640–4631's estimated distance of 8–13 kpc (Lemiere et al. 2009), and a typical pulsar kick velocity of 300 km s⁻¹ (Hansen & Phinney 1997), an unmodeled proper motion would change the measured braking index by less than one part in a million. Thus, neither of these effects can account for our measured braking index.

The pulsar's high-luminosity PWN has been argued to be powered by a relativistic outflow or wind from the neutron star. Under this assumption, predictions for the pulsar's spin-down history, and hence braking index, have been made (Gotthelf et al. 2014). The 0.2–10 TeV luminosity of the PWN powered by the source represents ∼6% of the pulsar's current spin-down luminosity (Gotthelf et al. 2014). A pulsar whose spin-down is driven solely by a particle wind would result in a braking index of one (Michel 1969; Manchester et al. 1985). Furthermore, a combination of magnetic dipole radiation and wind braking would result in a braking index with a value between one and three. In this case, the braking index as a function of , the fraction of spin-down power due to a particle wind, is given by (Lyne et al. 2015)

$$\begin{eqnarray}&&n=\displaystyle \frac{2}{\epsilon +1}+1.\end{eqnarray} \tag{ 6 }$$

This implies a maximal expected braking index of n = 2.89 for PSR J1640–4631. Indeed, a more thorough modeling of the pulsar and PWN system suggested an even smaller braking index, n ≈ 1.9 (Gotthelf et al. 2014), clearly at odds with our result.

A changing magnetic field has also been put forth as a possibility for a braking index that is different from three by the growth or decay of the field (Blandford & Romani 1988; Gourgouliatos & Cumming 2015), or a change in the angle between the magnetic and rotation axes (Lyne et al. 2013). In this case, n is given by

$$\begin{eqnarray}&&n=3+2\displaystyle \frac{\nu }{\dot{\nu }}\left(\displaystyle \frac{\dot{B}}{B}+\displaystyle \frac{\dot{\alpha }}{\mathrm{tan}\alpha }\right).\end{eqnarray} \tag{ 7 }$$

For the decaying field case, this would imply a magnetic field decay rate of ∼200 MG per century. This decay rate is very close to that predicted in some magneto-thermal evolutionary models (Viganò et al. 2013), and in this interpretation might be providing direct observational evidence of a decaying magnetic field. On the other hand, population synthesis studies find no strong evidence for field decay in the radio pulsar population as a whole (Faucher-Giguère & Kaspi 2006).

If a change in the alignment angle between the magnetic and rotational axes is the cause of the anomalous braking index, then either α is less than π/2 and the rotation and magnetic axes are moving toward alignment, or α is greater than π/2 and the rotation and magnetic axes are counter-aligning. If α is changing on order of the rate of the Crab pulsar (Lyne et al. 2013), the only pulsar for which such a change has been measured, at ∼1° per century, this would imply that α is ∼5° away from being an orthogonal rotator. This is at odds with the pulse profile being single peaked, since an orthogonal rotator would typically be seen to have emission as each pole enters our line of sight. In general, the value of α for pulsars can be independently determined by modeling of the gamma-ray or radio pulse profiles. PSR J1640–4631, however, is radio quiet; see Section 2.2. There are also no detected gamma-ray pulsations from PSR J1640–4631 (Gotthelf et al. 2014), and so the value of α is unknown for this source.

Another possibility to explain a braking index greater than three is to invoke higher order multipoles (Pétri 2015). A pure quadrupole, either a magnetic quadrupole, or a mass quadrupole leading to gravitational radiation (Blandford & Romani 1988), would yield a braking index of 5 and could coexist with the magnetic dipole to give a braking index between 3 and 5. Analogous to the case of a wind, the fraction of spin-down due to a quadrupole versus a dipole, _Q is (Palomba 2000):

$$\begin{eqnarray}&&{\epsilon }_{Q}=\displaystyle \frac{n-3}{5-n}.\end{eqnarray} \tag{ 8 }$$

In our case, this implies ∼8% of the spin-down is due to quadrupolar radiation. In the case of a mass quadrupole, this would imply that the pulsar has an ellipticity of ∼0.005, which cannot be reproduced by theoretically proposed dense matter equations of state, for a neutron star rotating at 4.84 Hz (Owen 2005). If such an ellipticity did exist, it would produce gravitational waves having a maximum strain of $\sim 4\times {10}^{-26}\left(\frac{12\;{\rm{kpc}}}{d}\right)$ at twice the spin period of the pulsar (Palomba 2000), which is far below the detection sensitivity of current technology.

The existence of a magnetic quadrupole is in principle testable with future X-ray polarimeter missions. X-ray polarization measurements of neutron stars are in principle sensitive to the magnetospheric configuration (van Adelsberg & Lai 2006; Taverna et al. 2014), be it a quadrupolar field structure or a change in the alignment of the spin and magnetic poles. The specific magnetic field structure of pulsars has a strong impact on the inferred magnetic field strength, as well as predicted radio and gamma-ray pulse profiles (Pétri 2015). Thus, X-ray polarimetric observations of PSR J1640–4631 could help us understand the origin of the pulsar's high n and shed light on the range of possibilities of neutron-star magnetic field structure.

Since the first measurement of the Crab's braking index in 1972 (Boynton et al. 1972), we have known that various physical mechanisms, such as angular momentum loss due to a wind, can result in a pulsar braking index less than the canonical dipole value. Our results for PSR J1640–4631 now show that other physics, such as the quadrupole moment of the magnetic field, affect the evolution of this source, and likely rotation-powered pulsars in general. Given that two other young, high-magnetic field pulsars have experienced glitches that resulted in a significant drop in the braking index (Antonopoulou et al. 2015; Archibald et al. 2015), it is clear that continuous study of braking indices provides an important window into additional physical processes at work in the youngest and most energetic of neutron stars.

This work made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. Parkes radio telescope is part of the Australia Telescope National Facility, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. We also thank an anonymous referee for helpful comments that improved the manuscript. R.F.A. acknowledges support from an NSERC Alexander Graham Bell Canada Graduate Scholarship. E.V.G. received support from the National Aeronautics and Space Administration through Chandra Award Number GO5-16061X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. V.M.K. receives support from an NSERC Discovery Grant and Accelerator Supplement, Centre de Recherche en Astrophysique du Quebéc, an R. Howard Webster Foundation Fellowship from the Canadian Institute for Advanced Study, the Canada Research Chairs Program, and the Lorne Trottier Chair in Astrophysics and Cosmology. Part of this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.