Free Access
Issue
A&A
Volume 587, March 2016
Article Number A99
Number of page(s) 11
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201526490
Published online 24 February 2016

© ESO, 2016

1. Introduction

Pulsars are rotating neutron stars with a magnetic field misaligned with respect to their rotation axis. The periodic pulse of light that we see is caused by a beam of radiation produced in the magnetosphere and pointing in our direction at each rotation. Such a rotating magnetic dipole loses energy at the expense of its rotational energy so decelerates in time. The mechanism by which this takes place has been investigated ever since the discovery of pulsars but has not yet been completely understood. The idea that the slowdown mechanism consists of a radiative torque on a rotating magnetic dipole and of the torque by which the pulsar drives the pulsar wind was proposed early (Goldwire & Michel 1969), together with the idea that the braking torque is described by a power-law relation between the frequency of rotation f and its derivative: = −Kfn.

The braking index n=f¨f/2\hbox{$n = f \ddot f/{\dot f}^2$} is expected to be three in the case of braking by pure dipole radiation and one in the case of a pulsar wind dominated torque (Michel & Tucker 1969). This proposal has stimulated a long-lasting effort to determine and classify pulsars by their braking indices. After 23 years of observations of the Crab pulsar Lyne et al. (1993) found that, between sudden jumps in rotational frequency (glitches), the rotational slowdown is described well by a power law with n = 2.51 ± 0.01. This braking index does not hint at any simple model for the braking mechanism. Long-term well-defined values for the braking index have been confirmed for four other pulsars (Livingstone et al. 2007), and they consistently give n< 3, suggesting that magnetic dipole radiation alone is not sufficient to account for the observed spin-down.

The study of the pulsar braking mechanism is made difficult by known irregularities in pulsar clocks, which take the form of sudden jumps in rotational frequency (glitches) or more gradual deviations from regular spin-down (timing noise). Recent systematic analyses of timing noise of the Crab and other pulsars show that it is not random, as was assumed until quite recently (Scott et al. 2003). The comparison of a large number of pulsars has demonstrated that the evolution of the rotational phase of young pulsars is dominated by long (~ 1 yr) relaxation periods following the occurrence of significant glitches, whereas older pulsars show quasi-periodic behavior with phase modulations on typical timescales between ~ 1 and 10 yr (Hobbs et al. 2010). Such quasi-periodic changes are sometimes correlated withdiscrete variations in slowdown rate and pulse profile (Lyne et al. 2010). These properties have been interpreted in terms of phenomena occurring in pulsar magnetospheres (Lyne et al. 2010). The periodically active pulsar PSR B1931+24, which exhibits on- and of-switching of the radio emission and drastic changes in braking torque, was proposed as a prototype of a pulsar-magnetosphere interaction system, in which a varying flow of charged particles drives the braking mechanism (Kramer et al. 2006). Fluctuations in the size of acceleration gaps were also considered as possible sources of variations in particle current flow and braking torque (Cheng 1987).

Pulsar glitches have not been considered by most authors as part of the same mechanism, but as a different phenomenon that is connected with internal pulsar dynamics (Reichley & Downs 1969; Baym et al. 1969; Anderson & Itoh 1975; Ho & Andersson 2012). Indeed, the persistent increase in slowdown rate of the Crab pulsar from 1969 through 1993 (Lyne et al. 1993) was interpreted as being caused by a decrease in the moment of inertia due to interaction of the internal superfluid core with the crust of the neutronstar. An alternative explanation in terms of torque increase was discarded, because it was expected to be accompanied by a change in the configuration of pulsar’s magnetic field that would likely induce a change of pulse profile and of the braking index. This was not found in measurements taken in 19691993 (Lyne et al. 1993).

In this paper we analyze the phase history of the Crab pulsar and find a very accurate mathematical description of its behavior. Such unifying description indicates that, in our opinion, glitches and timing noise are part of the same braking mechanism that undergoes sudden changes during glitches. Preliminary results of this analysis were presented by one of us (AČ) at the Prague Synergy 2013 Conference (“Interaction of a pulsar with the surrounding nebula: the case of Crab”1).

Very recently, Lyne et al. (2015) have published results of a similar analysis using 45 years of radio data on the rotational history of the Crab pulsar. We will compare their method with ours, pointing out similarities and differences and, more important, giving a different interpretation of the observed phenomenology. The plan of the paper is as follows. Section 2 touches on questions raised by optical timing observations of the Crab pulsar and their comparison with the Jodrell Bank radio ephemerides (JBE2). Section 3 deals with phase analysis of JBE data and discusses the changing braking law index. Section 4 contains an analysis of phase residuals and completes the description of the evolution of rotational phase of the Crab pulsar. Conclusions follow in Sect. 5. Some more technical details are presented in the Appendix.

2. Optical timing of the Crab pulsar

During the period from the end of 2008 through the end of 2009, which was characterized by the absence of significant glitches in the JBE, we obtained three sets of high signal-to-noise ratio optical observations of the Crab pulsar. The first set of data was obtained in October 2008 with the ultra-fast photon counter Aqueye (Barbieri et al. 2009), mounted on the Copernico telescope at Asiago, while the second and third sets were taken with a similar instrument, Iqueye (Naletto et al. 2009), mounted on the ESO New Technology Telescope at La Silla in January and December 2009.

We measured optical pulse arrival times during two-second intervals by correlating the incoming photon rate with the average pulse profile and define the starting point for the optical phase at the maximum of the main pulse as defined by the template. Optical phase residuals during a typical observation run are Gaussian-distributed with the width consistent with photon-counting noise. In this way a typical one-hour observation yields a local phase model with statistical phase uncertainty of ~0.3 μs at Le Silla and ~0.45 μs at Asiago (Germanà et al. 2012; Zampieri et al. 2014).

Since our data were taken with two different instruments at very different locations on Earth, we chose to analyze them independently. Our analysis of optical data is illustrated in Fig. 1, where JBE radio phase residuals (calculated as described in the next section) are also shown. It turned out that the Iqueye data fit a braking index model (BIM, see below) with an unexpectedly high precision, i.e. the frequency and frequency derivative determined from January and December 2009 data are so tightly constrained that essentially a single braking index law solution with n = 2.437 also fixes the relative phase between January and December (left panel in Fig. 1). This leads to a phase description with respect to which the optical data points deviate by less then 10 μs. Thus, the 2009 n = 2.437 braking law solution appears to be a good baseline for measuring the 2009 Crab phase noise, which is clearly detected at the microsecond level during this quiet period of Crab history. The phase predicted by this solution agrees (to within the well-documented 150−250 μs radio phase lag) with radio ephemerides phase on the dates of when we gathered our data and differs from radio ephemerides on the whole interval by less than 3 ms (Zampieri et al. 2014). The Asiago October 2008 data were added to the analysis after checking the consistency of our timing protocols on the two sites and the equality of timing response of the two instruments. These data are again consistent with radio ephemerides, but do not follow the 2009 n = 2.437 braking index law as well, because they miss the prediction by almost 8 ms. It turns out that our complete set of optical data can be fitted by a n = 2.476 braking index law with respect to which radio phase differs by less than 4 ms in the whole 14-month interval (right panel in Fig. 1). However, in this case the optical phase distinctly shows large excursions (up to 50μs per day) with respect to this braking index law. The good fit of our January and December data to the braking index law appears as a rare coincidence, but it stimulated us to ask the question whether there may be a natural baseline with respect to which pulsar phase noise should be measured and how well and for how long a braking index law can approximate the phase history of the Crab pulsar.

thumbnail Fig. 1

Residuals of optical and radio phase with respect to a BIM model. Left: best-fitting BIM to 2009 Iqueye data (n = 2.437). Right: best-fitting BIM to 20082009 Aqueye and Iqueye data (n = 2.476). Optical residuals are plotted as red data points clustered at the time of optical observations, while Jodrell Bank residuals (over the same interval of time) are represented by a gray line crossed by error bars at epoch dates of JBE. Insets: zoom of average optical phase residuals during observation nights with 1σ error bars. The scales of the y-axis in the insets are very different.

Our analysis leads to the following conclusions:

  • The optical phase always leads the radio one. For three sets of observations, the time lapse between optical measurement and JBE reported radio phase (the latter refers to the center of the main pulse) was between 160 and 260 μs (Germanà et al. 2012; Zampieri et al. 2014). Therefore, there is no evidence for any significant change in the delay between the three epochs. Accurate estimates of the X-ray-radio delay have been reported by Rots et al. (2004) (344 ± 40 μs) and show that it does not change on a timescale of several years. No other difference between optical and radio outside the quoted interval of uncertainty has been found. We note that this delay is also consistent with other recent measurements performed in the optical band (e.g., Oosterbroek et al. 2006, 2008).

  • A braking index law can reproduce the phase history during the studied 14-month period to within one turn with a narrow range of braking parameters, yet residuals with respect to such a description represent timing noise present on all wavebands and thus most likely reflect genuine variations in rotation of the Crab pulsar. This then motivates us to use JBE data for studying the long-term intrinsic rotation properties of the pulsar.

  • The small residuals shown in the inset of lefthand part of Fig. 1 are real and should be considered as typical of short, daily timescale phase noise during a quiet glitch-less period of pulsar history. They are shown again in Fig. 2, together with polynomial fits through data points on the left, and frequency and frequency derivative residuals corresponding to these polynomials on the right. Thus, these data suggest “typical” frequency noise on a daily scale at the level of ~ 10-8 s-1 and “typical” frequency derivative noise on a daily scale at the level of a few times 10-13 s-2. These estimates are in reasonable agreement with the difference in frequency and frequency derivative derived from optical and radio data (Zampieri et al. 2014) which are January 15. 2009: foptfradio = −5.64 × 10-9 s-1, optradio = 1.55 × 10-14 s-2; and December 15. 2009: foptfradio = −6.45 × 10-9 s-1, optradio = −6.86 × 10-15 s-2.

2.1. Braking index model implementation

The braking index model implies that the phase ϕ can be expressed as ϕ(t)=c+a(1+bt)s,\begin{eqnarray} \varphi(t) = c + a(1+b t)^s, \label{EQN1} \end{eqnarray}(1)where c is an integration constant that can be used to adjust the initial phase, a and b are parameters that are related to the age and frequency of the pulsar at t = 0, and s=n2n1\hbox{$s =\frac{n-2}{n-1}$} is the braking parameter. The frequency and its derivative are then the following functions of time \hbox{$f = \dot\varphi=s a b (1+b t)^{s-1}$} and =¨ϕ=s(s1)ab2(1+bt)s2\hbox{$\dot f=\ddot\varphi=s(s-1) a b^2 (1+b t)^{s-2}$}. The last equation leads to the familiar form of the braking index law: =[b(s1)(sab)1s1]fs2s1,\begin{equation} \dot f=\left [b(s-1)\left (s a b \right )^\frac{1}{s-1} \right]f^\frac{s-2}{s-1} \label{EQN2} , \end{equation}(2)therefore n=s2s1\hbox{$n=\frac{s-2}{s-1}$}.

To compare our optical data with radio data as documented in JBE, we devised a numerical procedure that expresses the phase as the number of turns that the pulsar has made since the first pulse arrival time on May 15. 1988; i.e., φ(t) = 0 at MJD = 47 296.0000003712. More details are given in Appendix A.

thumbnail Fig. 2

January (top) and December (bottom) 2009 optical residuals (Fig. 1) fitted to polynomials that reveal “typical” noise in pulsar frequency and frequency derivative, as shown on the right, in solid line for frequency residuals (left scale) and dashed for frequency derivative residuals (right scale).

3. Braking episodes

That the phase history of the Crab pulsar can be approximated by a BIM to within one turn during a 14-month period raises the question of how long such a description can go on? This idea can best be analyzed by using the publicly available JBE Crab pulsar ephemerides since 1988 because they represent the most complete and uniform description of Crab rotational history, which has been shown to be consistent with data in other wavebands. As a first step in finding periods during which the braking index may be sufficiently constant, we plot residuals of JBE-published frequency derivatives with respect to frequency derivatives predicted by a BIM with s=n2n1=1/3\hbox{$s=\frac{n-2}{n-1}=1/3$} (n = 2.5) and the parameters a and b in equation (1) chosen by a least squares fit minimizing residuals res = JBE. These residuals are shown in Fig. 3.

The graph clearly shows that on average, the n = 2.5 braking index law describes the 26-year phase history reasonably well. However, large systematic changes in slope immediately after some (major) glitches signal that during such periods, the braking index must be different from the average value. Therefore, we divide the 26-year interval into ten subintervals, henceforth called episodes (see Fig. 3), each corresponding to periods characterized by the absence of pronounced changes in the braking index. We expect that the phase history during an episode can be largely approximated by a constant braking index law, as in the case of the interval of comparison of optical and radio data. At this stage in our discussion, the exact dates of the ten episodes are not yet well defined. In the next sections we show how they can be better constrained.

thumbnail Fig. 3

Residuals of JBE frequency derivatives with respect to a constant braking index law fit with n = 2.5, calculated over the whole JBE interval (from 1988 May through 2015 June 15th). The vertical gray lines denote the occurrence of glitches as reported by Espinoza et al. (2011)2 (see Table A.1). Nine arrows delimit chosen episodes and are placed at times when (major) glitches appear to change the braking index. The two dashed vertical lines mark the beginning and the end of the data set. The gray broken line indicates the second derivative of the continuous phase function defined in the text. Points used for the fit with the braking law model ϕj(t) are displayed in red and green, while points in black are excluded from the fit, as explained in the text. Some post glitch residuals with the value below −0.5 × 10-12 s-2 go beyond the scale and are not shown.

To be able to properly categorize a glitch as a discontinuous change in frequency and frequency derivative, one would need densely distributed high resolution optical data from such an event. Unfortunately, they are not available yet. To get an idea of how the spin-down behaves before and after the glitch event, we used the public available data (Lyne et al. 1993)2 and built the phase history since May 15, 1988, assuming that the phase prediction by the suggested third-degree polynomial is accurate to a (small) fraction of a period during each ephemerides epoch and that the phase is continuous between successive periods. We built up a table of pulsar’s integer radio phase at the exact Julian dates of pulse arrival times given by JBE. This table T is a table of entries of main pulse arrival times and integer number of turns the pulsar has made from the starting date. This table is broken into ten episodes according to boundaries shown in Fig. 3 and determined as explained below (see also Appendix A).

Table 1

Parameters of local fits for different j episodes, expressed in the form of a Taylor series of the form ϕj(t)=ψj+νjt+12ν̇t2j+16¨νjt3+124...νjt4\hbox{$\varphi_j(t)=\psi_j+\nu_j t+\frac{1}{2}\dot \nu_j t^2+\frac{1}{6}\ddot \nu_j t^3 +\frac{1}{24}\dddot \nu_j t^4$} where the coefficients are constrained to fit the braking index law: i.e., this expression is a Taylor series expansion of Eq. (1).

Inside each episode we seek a braking index law ϕj(t) in the same form as in Eq. (1) with sj as an additional fitting parameter. Our goal is to obtain the sj for the braking law that fits the particular episode, in such a way that it gives the smallest residuals over the whole episode. It turns out that this goal can only be achieved by using data in the apparently quiet part of the episode, disregarding scattered data immediately after the glitch at the beginning of the episode. Data points whose phase is not considered in the fit are colored black in Fig. 3 (we comment on this choice below). In this way we obtain local fits ϕj(t) = cj + aj(1 + bjt)sj, valid on complete intervals of episodes {Tjb,Tj+1b}\hbox{$\lbrace T^b_j,T^b_{j+1}\rbrace$}, where Tb is the starting MJD for each episode.

These local fits are joined into a continuous curve Φ(t) by choosing cj + 1 in such a way that ϕj+1(Tj+1b)=ϕj(Tj+1b)\hbox{$\varphi_{j+1}(T^b_{j+1}) = \varphi_{j}(T^b_{j+1})$}. Residuals of the fit (R(t) = T−Φ(t)) are shown in Fig. 4 at all data points, and the parameters of ϕj(t),expressed as a Taylor series expansion of Eq. (1), are given in Table 1. When varying the parameters within the reported errors, the solution deviates from the best fits by less than 0.001 turn3. The fits to those residuals, discussed in the next section, are also shown as blue and cyan curves on top of the R(t) dots in Fig. 4 (dots also include black data points in Fig. 3). Residuals in Fig. 4 show that, during each episode, it is possible to find a local braking law with respect to which Rt)(\hbox{$R\left (t\right )$} is never more than a few turns, even if thepulsar makes billions of turns during the same period. It must be emphasized that residuals that are as small as those shown in Fig. 4 can only be obtained by fitting the phase on different episodes Tj with markedly different braking law indices, varying from ~ 2.1 to ~ 2.6.

As already mentioned in the discussion of fitting optical data (see Fig. 1), the precise value of the braking index on the episode may depend on the selection of data points after the glitch. Looking at Fig. 3, different choices for the black, green, and red points can possibly be made. However, the requirement that phase residuals be as small as only a few turns very clearly narrows down an already narrow range of acceptable values of the braking index for each episode. Finally, since Φ(t) is a continuous function, the residuals R(t) must also be continuous at the boundaries between episodes. This requirement narrows down the dates of boundaries between episodes to values Tjb\hbox{$T^b_j$} as listed in Table 1.

thumbnail Fig. 4

JBE phase residuals R(t) with respect to local fits ϕj(t) calculated for the 10 chosen episodes (see Fig. 3). Gray vertical lines are plotted at the dates of reported glitches (Espinoza et al. 2011). The braking index during episodes is shown in orange. Sinusoidal fits to residuals, discussed below, are shown in blue and cyan to distinguish episodes. Gray points at the bottom show ten times the difference between R(t) and sinusoidal fits Rj(t) at data points; the horizontal dashed lines bracket these final residuals with their standard deviation of 0.057 turns; the distribution of final residuals has wider wings than a Gaussian.

4. The rotational phase history of the Crab pulsar

It is remarkable that the evolution of the rotational phase of the Crab pulsar can be split into two parts: a regular phase Φ(t) consisting of a constant braking index during a given episode, but different for each episode, and a small residual phase R(t) that, during all this time, wanders by no more than 35 turns, to be compared with the ~ 2.5 × 1010 turns the pulsar has made during the 9800 days covered by the JBE (see Fig. 4).

According to the dates reported in JBE and in Espinoza et al. (2014), jumps in braking index are clearly related to large glitches (see Table A1), while small glitches do not leave a notable imprint on the curve R(t). To recover the significance of small glitches from JBE, we calculated the second derivative of the function Rs(t), which is constructed as a differentiable function from tabular values of R(t) by cubic spline interpolation. The result is shown in Fig. 5 (left panel). In this representation, small glitches, except possibly a few (n. 10, 11, 13), show up as spikes barely above the general noise in ¨Rs(t)\hbox{${\ddot R}_{\rm s}(t)$}. Only nine glitches (n. 6, 7, 9, 14, 16, 20, 23, 24, 25) show second-derivative ¨Rs(t)\hbox{${\ddot R}_{\rm s}(t)$} beyond 2 × 10-14 s-2, which is smaller than the second derivative of optical phase residuals derived from the January 2009 data (¨RJan10-13s-2\hbox{${\ddot R}_{\rm Jan}\approx 10^{-13}~{\rm s}^{-2}$} see Fig. 2). All of them are clearly related to a change in the braking index (Fig. 4).The righthand panel of Fig. 5 shows phase residuals Rs(t), together with phase residuals calculated as suggested by explanatory notes of JBE. Breaks (a few thousands of a turn) between the ephemerides dates, occur because the suggested expression for calculating the phase increment, ϕ(t)=ϕ(t0)+1P(tt0)12P2(tt0)2+132P3(tt0)3,\begin{equation} \varphi(t)=\varphi(t_0)+\frac{1}{P}(t-t_0)-\frac{1}{2}\frac{\dot P}{P^2}(t-t_0)^2+\frac{1}{3}\frac{{\dot P}^2}{P^3}(t-t_0)^3 \label{Eq3} , \end{equation}(3)(P is the nominal period) does not yield an exact integer number of turns between two entries in Table T. We suspect that the actual pulsar phase noise variation (Wong et al. 2001) is the dominant cause for those phase residuals breaks. The truncated Taylor series in Eq. (3) is probably not sufficient to take complete account of phase noise variation. Intrinsic timing noise has also been confirmed by optical observations (Fig. 2), and according to Zampieri et al. (2014), the difference in frequency derivative obtained from radio ephemerides and from optical observations is consistent with the difference between the red and blue curve in the righthand panel of Fig. 5.

thumbnail Fig. 5

Left: second time derivative of the residual phase Rs(t), a third-order spline fit through phase residuals from Table T. The orange vertical lines are drawn at dates of reported glitches and labeled according to Table A.1. Episodes are shaded intermittently as light gray and white. Right: phase residuals as a function of time for a short interval during episode 3. The red line shows Rs(t), while the blue line is calculated as suggested by explanatory notes of JBE (see text for details). Error bars show JBE quoted arrival-time uncertainties.

The curves of residuals in Fig. 4 all have a characteristic shape. It is customary to express functions describing sets of data with a linear combination of a certain Hilbert space basis (eigenvectors) and natural to choose such a basis so that the concrete experimental result can be described by the fewest components. Having tried different possibilities (including the widely used combination of sinusoids + decaying exponentials), we find that the residuals can be described by only two complex Fourier components: Rj(t)=Ψj+Ajeλj(tTjb)sin(ωj(tTjb)+δj)\begin{eqnarray} R_j(t)&=&\Psi_j+A_j {\rm e}^{-\lambda_j (t-T^b_j)}\sin\left (\omega_j (t-T^b_j )+\delta_j\right) \nonumber \\ &&\quad + {\overline A}_j {\rm e}^{-{\overline \lambda}_j (t-T^b_j )}\sin\left ({\overline \omega }_j(t-T^b_j)+\overline {\delta}_j\right). \label{EQN4} \end{eqnarray}(4)We believe that this mathematical approach leads to a simpler description of pulsar noise. The coefficients are given in Table 2 and are sufficient to produce a fit of R(t) without systematic trends in the residuals, as shown in gray at the bottom of Fig. 4. The standard deviation of the distribution of phase residuals is 0.057. It can be compared with the standard deviation of the set of fractional phase residuals between ephemerides epochs, which is 0.063. The declared arrival-time uncertainty in JBE is 0.01 or 300 μs. The details of these fits, together with some other pertinent information, are shown in Fig. 6.

thumbnail Fig. 6

Sinusoidal phase functions Rj(t) (red lines) that fit the JBE phase residuals (blue points) for episodes from j = 1 to 10 (the short episode 5 is included in the graph showing episode 6). Time on the abscissa is in MJD, and the scale is the same for all episodes. Green dots are numbered as in Espinoza et al. (2011) at dates of published glitches with the ordinate as calculated from Rj(t). Ordinate scales are different and adjusted to different amplitudes of oscillations. Each panel shows data of the complete episode and includes the first point of the next episode. The difference between horizontal dotted lines represents a measure of the strength of the perturbation Rj(t) caused by the glitch. Dashed gray curves are the plots of ϕj(t)−ϕj−1(t), to show the difference of phases between two contiguous episodes. Gray arrows point to dates of the > 100 MeV X-ray flares detected by Fermi (Abdo et al. 2011; Ojiha et al. 2013; Bühler et al. 2012) and AGILE (Tavani et al. 2011; Striani et al. 2013).

Table 2

Coefficients of the fitting function Rj(t) in Eq. (2).

The ansatz in Eq. (4) is purely mathematical and is an almost satisfactory description of phase residuals R(t). A posteriori, we note that the first sinusoid has very low values of ω so is essentially a decaying component, while in all episodes but one, the second sinusoid is not damped.

thumbnail Fig. 7

Left: dispersion measure DM(t) (blue dots), braking index n(tτ) for τ = 1100 days (continuous red line). The dashed red line shows n(t). Right: correlation between braking index (n) and dispersion measure DM. The correlation coefficient is 0.7.

It should be understood that glitches may not be modeled exactly by our analysis because the time resolution of JBE data is not sufficient to describe the exponential decay of a glitch that lasts a few days (Wong et al. 2001). However, because of all that has been said, the JBE data are reliable as to the global phase behaviour and on timescales longer than about a month. In this respect, the decaying component at the beginning of several episodes resembles the usual after-glitch recovery behavior frequently modeled in the literature through a single-parameter exponential (e.g., Lyne et al. 2015, with a characteristic timescale of 320 days).

Figures 4 and 6 clearly illustrate the meaning of the split of the phase behavior into the braking index part Φ(t) and residuals R(t). It is quite clear that toward the end of an episode, residuals relax to the perfect braking law solution and eventually oscillate for years by only a very small fraction of a turn. On the other hand, as hinted in Fig. 6, the braking-law part of an episode (ϕj(t)) generally differs quite considerably from the braking law of the previous episode (ϕj−1(t)). In fact, if the difference ϕj(t)−ϕj−1(t) was plotted to the end of episodes, it would reach values several tens or a hundred times the value ΔN. An exception is the short episode 8, which follows a relatively weak glitch. For this episode R(t) and ϕj(t)−ϕj−1(t) are comparable, and therefore, the split between the braking index part and residuals is not as clear cut. Figure 6 also confirms that, within limits allowed by the time resolution of available data, the breaks between episodes occur on the dates of reported glitches. It appears plausible to classify glitches into two groups: those for which the change in the phase Φ(t) dominates residuals R(t) by many factors of ten (6, 7, 9, 16, 20, 23, 24, 25) and those where the integral change in ϕj(t)−ϕj−1(t) is comparable or insignificant with respect to the amplitude of R(t) (all other glitches beyond #6).

5. Discussion and conclusions

Our analysis shows that the phase evolution of the Crab pulsar can be described as a series of constant braking-law episodes, with the braking index changing abruptly after each episode in the range of values between 2.1 and 2.6. Phase residuals with respect to such a smooth phase description amount to only a few turns in ~ 109 turns executed during an episode. The split between the smooth braking-law-dominated part and residuals is not mathematically unique, but requirements that phase residuals be as small as only a few turns and that the phase between ephemerides epochs clearly converges to the characteristic braking index solution of the episode narrow the choice of braking index parameters.

A similar conclusion concerning the behavior of the braking index has been recently obtained by Lyne et al. (2015) from an independent analysis of 45 years of radio data on the rotational history of the Crab pulsar. Results obtained from a fit with a single Taylor series returns a behavior of f and (their Figs. 1 and 2) very similar to the one shown in Fig. 3 (with a best fitting n = 2.34). Variations in the braking index in the period between 1996 and 2006, characterized by a high concentration of glitches, were also noted by Lyne et al. (2015). While they consider it a weak, unexplained effect on the background of the previous rotational history (described by a simple slowdown with braking index 2.519(2)), we offer here a different interpretation that can account for the overall timing irregularities of the Crab pulsar.

According to our interpretation, glitches and abrupt changes in the braking mechanism may be part of the same physical process that also drives semi-periodic timing noise between glitches (Fig. 4). In fact, JBE data provide an interesting correlation of the braking index and dispersion measure. Namely, the dispersion measure (as listed in JBE) follows the braking index with a time delay of 1100+450-250\hbox{$^{+450}_{-250}$} days as shown in Fig. 7. The delayed response of dispersion measure to the value of braking index lends support to the idea of Anderson & Itoh (1975) that change of the braking index has to do with a pulsar-wind-driving torque and is also consistent with the idea that eddy currents threading the pulsar nebula ionize it and thus inject varying amounts of free electrons.

From this new perspective, glitches, timing noise, changes in braking torque, and dispersion measure appear to be part of a common mechanism that jumps between different braking modes. Distinctly long timescales of timing noise oscillations suggest that the mechanism can hardly be connected to phenomena occurring within the pulsar. It appears plausible that the “instantaneous” change in braking torque is caused by some instability, which varies the configuration of the external electromagnetic field and currents in the nebular plasma through which the pulsar interacts with its nebula. Such an interaction is needed in order to understand the acceleration and energizing mechanisms of the pulsar nebula (Trimble 1968; Weyler & Panagia 1978), such as the highly dynamical flow of relativistic particles in the form of equatorial wind and polar jets, as seen in Hubble Space Telescope and Chandra images (Hester et al. 2002).

The possible occurrence of plasma instabilities causing “instantaneous” changes in the braking mechanism is expected to produce observable changes in radiation from the Crab nebula. It is tempting to consider the possibility (Cerutti et al. 2014) that plasma instabilities occurring through magnetic field line reconnection drive the recently observed gamma ray flares (Abdo et al. 2011; Tavani et al. 2011; Ojiha et al. 2013; Striani et al. 2013; Bühler et al. 2012). Arrows in Fig. 6 point to instants when the six observed flares occurred in the Crab nebula. The same mechanism, which also appears to be acting in the solar corona to produce gamma ray flares (Ajello et al. 2014), may be able to provide a sudden short release of braking torque by disconnecting the magnetic field lines from the pulsar from those threading the nebular plasma. The same mechanism may also be responsible for sharp increases in dispersion measure (Fig. 7, left) by emitting highly energetic particles into the neutral nebula, thus ionizing it.

The back reaction of plasma instabilities on the pulsar, possibly associated to these flares, has not been studied yet, but it does not seem to be simultaneous with the occurrence of the flare (Weisskopf et al. 2013; Zampieri et al. 2013). In view of the delayed correlation between dispersion measure and braking index, this appears understandable – perturbations caused by magnetic reconnection travel long distances before reaching the pulsar or before permeating the nebula. Nevertheless, the interflare time inferred from numerical simulations (hundreds of days) (Mignone et al. 2013; Porth et al. 2014; Cerutti et al. 2014) appears to be broadly consistent with the timescales observed in oscillations of JBE phase residuals. It is quite remarkable that the occurrence of the reported gamma-ray flares is also consistent with such a timescale.

Our ultra-fast optical observations of the Crab pulsar with Aqueye and Iqueye stimulated the line of research presented here. Future simultaneous radio and optical timing measurements, as well as optical imaging and gamma-ray observations, will be crucial for revealing the source of the braking mechanism, in particular if it is located in the external electromagnetic field through which the pulsar interacts with the surrounding plasma, as suggested by our results.


3

As noted before, the split T = Φ(t) + R(t) between Φ and R is not unique. The one presented here is the result of our attempts to find one, where the contribution of R(t) is as small and regular as we can find.

Acknowledgments

This work is based on observations made with ESO Telescopes at the La Silla Paranal Observatory under program IDs 082.D-0382 and 084.D-0328(A) and on observations collected at the Copernico Telescope (Asiago, Italy) of the INAF-Osservatorio Astronomico di Padova. We acknowledge the use of the Crab pulsar radio ephemerides available on the web site of the Jodrell Bank Radio Observatory (http://www.jb.man.ac.uk/~pulsar/crab.html, see Lyne et al. 1993). This research has been partly supported by the University of Padova under the Quantum Future Strategic Project, by the Italian Ministry of University MIUR through the program PRIN 2006, by the Project of Excellence 2006 Fondazione CARIPARO, and by INAF-Astronomical Observatory of Padova. One of us (A. Č) would like to express his gratitude to the relativity group at the Silesian University in Opava for their support and friendship. The authors wish to thank the referee, Patrick Weltevrede, for his valuable comments that helped to improve the paper.

References

  1. Abdo, A. A., Ckermann, M., Ajello, M., et al. 2011, Science, 331, 739 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  2. Ajello, M., Albert, A., Allafort, A., et al. 2014, ApJ, 789, 20 [NASA ADS] [CrossRef] [Google Scholar]
  3. Anderson, P. W., & Itoh, N. 1975, Nature, 256, 25 [NASA ADS] [CrossRef] [Google Scholar]
  4. Baym, G., Pethick, C., Pines, D., & Ruderman, M. 1969, Nature, 224, 872 [NASA ADS] [CrossRef] [Google Scholar]
  5. Barbieri, C., Naletto, G., Occhipinti, T., et al. 2009, J. Mod. Opt., 56, 261 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bühler, R., Scargle, J. D., Blandford, R. D., et al. 2012, ApJ, 749, 26 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2014, ApJ, 782, 104 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cheng, K. S. 1987, ApJ, 321, 799 [NASA ADS] [CrossRef] [Google Scholar]
  9. Espinoza, C. M., Lyne, A. G., Stappers, B. W., Kramer, M., et al. 2011, MNRAS, 414, 1679 [NASA ADS] [CrossRef] [Google Scholar]
  10. Espinoza, C. M., Antonopoulou, D., Stappers, B. W., Watts, A., Lyne, A. G., et al. 2014, MNRAS, 440, 2755 [NASA ADS] [CrossRef] [Google Scholar]
  11. Germanà, C., Zampieri, L., Barbieri, C., et al. 2012, A&A, 548, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Goldwire, H. C., & Michel, F. C. 1969, ApJ, 156, L111 [NASA ADS] [CrossRef] [Google Scholar]
  13. Hester, J. J., Mori, K., Burrows, D., Gallagher, J. S., et al. 2002, ApJ, 577, L49 [NASA ADS] [CrossRef] [Google Scholar]
  14. Ho, W. C., & Andersson, N. 2012, Nat. Phys., 8, 787 [NASA ADS] [CrossRef] [Google Scholar]
  15. Hobbs, G., Lyne, A. G., & Kramer, M. 2010, MNRAS, 402, 1027 [NASA ADS] [CrossRef] [Google Scholar]
  16. Kramer, M., Lyne, A. G., O’Brien, J. T., Jordan, C. A., & Lorimer, D. R. 2006, Science, 312, 549 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  17. Livingstone, M. A., et al. 2007, Ap&SS, 308, 317 [NASA ADS] [CrossRef] [Google Scholar]
  18. Lyne, A. G., Pritchard, R. S., & Graham-Smith, F. 1993, MNRAS, 265, 1003 [NASA ADS] [CrossRef] [Google Scholar]
  19. Lyne, A., Hobbs, G., Kramer, M., Stairs, I., & Stappers, B. 2010, Science, 329, 408 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  20. Lyne, A. G., Jordan, C. A., Graham-Smith, F., et al. 2015, MNRAS, 446, 857 [NASA ADS] [CrossRef] [Google Scholar]
  21. Michel, F. C., & Tucker, W. H. 1969, Nature, 223, 277 [NASA ADS] [CrossRef] [Google Scholar]
  22. Mignone, A., Striani, E., Tavani, M., & Ferrari, A. 2013, MNRAS, 436, 1102 [NASA ADS] [CrossRef] [Google Scholar]
  23. Naletto, G., Barbieri, C., Occhipinti, T., et al. 2009, A&A, 508, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Ojha, R., Buehle, R., Hays, E., & Dutka, M. 2013, Atel, 4855 [Google Scholar]
  25. Oosterbroek, T., de Bruijne, J. H. J., Martin, D., et al. 2006, A&A, 456, 283 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Oosterbroek, T., Cognard, I., Golden, A., et al. 2008, A&A, 488, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Porth, O., Komissarov, S. S., & Keppens, R. 2014, MNRAS, 438, 278 [NASA ADS] [CrossRef] [Google Scholar]
  28. Reichley, P. E., & Downs, G. S. 1969, Nature, 222, 229 [NASA ADS] [CrossRef] [Google Scholar]
  29. Rots, A. H., Jahoda, K. L., & Andrew, G. 2004, ApJ, 605, L129 [NASA ADS] [CrossRef] [Google Scholar]
  30. Scott, D. M., Finger, M. H., & Wilson, C. A. 2003, MNRAS, 344, 412 [NASA ADS] [CrossRef] [Google Scholar]
  31. Striani, E., Tavani, M., Vittorini, V., et al. 2013, ApJ, 765, 52 [NASA ADS] [CrossRef] [Google Scholar]
  32. Tavani, M., Bulgarelli, A., Vittorini, V., et al. 2011, Science, 331, 736 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  33. Trimble, V. 1968, AJ, 73, 535 [Google Scholar]
  34. Weisskopf, M. C., Tennant, A. F., Arons, J., et al. 2013, ApJ, 765, 56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Weyler, K. W., & Panagia, N. 1978, A&A, 70, 419 [NASA ADS] [Google Scholar]
  36. Wong, T., Backer, D. C., & Lyne, A. G. 2001, ApJ, 548, 447 [NASA ADS] [CrossRef] [Google Scholar]
  37. Zampieri, L., Verroi, E., Naletto, G., et al., 2013, Atel, 4878 [Google Scholar]
  38. Zampieri, L., Čadež, A., Barbieri, C., et al. 2014, MNRAS, 439, 2813 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Crab phase data

We construct the decadal timing solution of the Crab pulsar as a series of pairs Tj = { MJDj + (1 / 86400)TOAj,zj }, where MJDj is the mean Julian date of the jth entry of JB ephemerides, TOAj is the first arrival time of the pulse at MJDj (the tJPL entry in JBE), and zj is the integer number of turns the pulsar has made from the starting date May 15, 1988 (MJD1 = 47 296). Integers zj are calculated by summing the integer number of turns Nj that the pulsar has made between MJDj−1 and MJDj (zj=k=1jNk\hbox{$z_{j} = \sum{}^{j}_{k=1}N_{k}$}). Nk is calculated using the JB published values of frequency and frequency derivative. Using the formula suggested by explanatory notes to JBE, we calculate it as follows.

Let Δt=43200Tj(Tj1)\hbox{$\Delta t =43200\left ({\bf T}_j-{\bf T}_{j-1}\right )$}, and let νj and \hbox{${\dot \nu}_j$} be the frequency and frequency derivative listed at the jth ephemerides entry in JBE. Define Pj = 1 /νj and \hbox{${\dot P}_j=-{\dot \nu_j}{\nu_j}^{-2}$}, then ΔN=(1Pi1+1Pi)Δt(j1Pj12jPj2)Δt22+(j12Pj13+j2Pj3)Δt33\hbox{$\Delta N=\left(\frac{1}{P_{i-1}}+\frac{1}{P_i}\right)\Delta t-\left(\frac{\dot P_{j-1}}{{P_{j-1}}^2}-\frac{\dot P_{j}}{{P_{j}}^2}\right )\frac{\Delta t^2}{2}+\left(\frac{{\dot P_{j-1}}^2}{{P_{j-1}}^3}+\frac{{\dot P_{j}}^2}{{P_{j}}^3}\right )\frac{\Delta t^3}{3}$} and Nj = IntegerPartN]. The evaluation of ΔN in most cases yields an integer, as it should. However, in 25 cases, almost all of them occurring at the time of a glitch, the calculated number of turns between successive TOA’s has a fractional part that is inconsistent with the quoted TOA accuracy. The fractional phase errors are shown in Fig. A.1.

Table A.1 lists all glitches reported in Espinoza et al. (2011)4 and the starting dates of episodes 2 to 10.

thumbnail Fig. A.1

Fractional part of calculated number of turns between successive TOAs’ fractional phase error. Vertical lines denote the dates of reported glitches, starting from glitch No. 5 in Table A.1.Horizontal dashed lines are at ± 0.062, the variance of the (wide winged) phase error distribution.

Table A.1

Numbers and dates of glitches (from Espinoza et al. 2011).

All Tables

Table 1

Parameters of local fits for different j episodes, expressed in the form of a Taylor series of the form ϕj(t)=ψj+νjt+12ν̇t2j+16¨νjt3+124...νjt4\hbox{$\varphi_j(t)=\psi_j+\nu_j t+\frac{1}{2}\dot \nu_j t^2+\frac{1}{6}\ddot \nu_j t^3 +\frac{1}{24}\dddot \nu_j t^4$} where the coefficients are constrained to fit the braking index law: i.e., this expression is a Taylor series expansion of Eq. (1).

Table 2

Coefficients of the fitting function Rj(t) in Eq. (2).

Table A.1

Numbers and dates of glitches (from Espinoza et al. 2011).

All Figures

thumbnail Fig. 1

Residuals of optical and radio phase with respect to a BIM model. Left: best-fitting BIM to 2009 Iqueye data (n = 2.437). Right: best-fitting BIM to 20082009 Aqueye and Iqueye data (n = 2.476). Optical residuals are plotted as red data points clustered at the time of optical observations, while Jodrell Bank residuals (over the same interval of time) are represented by a gray line crossed by error bars at epoch dates of JBE. Insets: zoom of average optical phase residuals during observation nights with 1σ error bars. The scales of the y-axis in the insets are very different.

In the text
thumbnail Fig. 2

January (top) and December (bottom) 2009 optical residuals (Fig. 1) fitted to polynomials that reveal “typical” noise in pulsar frequency and frequency derivative, as shown on the right, in solid line for frequency residuals (left scale) and dashed for frequency derivative residuals (right scale).

In the text
thumbnail Fig. 3

Residuals of JBE frequency derivatives with respect to a constant braking index law fit with n = 2.5, calculated over the whole JBE interval (from 1988 May through 2015 June 15th). The vertical gray lines denote the occurrence of glitches as reported by Espinoza et al. (2011)2 (see Table A.1). Nine arrows delimit chosen episodes and are placed at times when (major) glitches appear to change the braking index. The two dashed vertical lines mark the beginning and the end of the data set. The gray broken line indicates the second derivative of the continuous phase function defined in the text. Points used for the fit with the braking law model ϕj(t) are displayed in red and green, while points in black are excluded from the fit, as explained in the text. Some post glitch residuals with the value below −0.5 × 10-12 s-2 go beyond the scale and are not shown.

In the text
thumbnail Fig. 4

JBE phase residuals R(t) with respect to local fits ϕj(t) calculated for the 10 chosen episodes (see Fig. 3). Gray vertical lines are plotted at the dates of reported glitches (Espinoza et al. 2011). The braking index during episodes is shown in orange. Sinusoidal fits to residuals, discussed below, are shown in blue and cyan to distinguish episodes. Gray points at the bottom show ten times the difference between R(t) and sinusoidal fits Rj(t) at data points; the horizontal dashed lines bracket these final residuals with their standard deviation of 0.057 turns; the distribution of final residuals has wider wings than a Gaussian.

In the text
thumbnail Fig. 5

Left: second time derivative of the residual phase Rs(t), a third-order spline fit through phase residuals from Table T. The orange vertical lines are drawn at dates of reported glitches and labeled according to Table A.1. Episodes are shaded intermittently as light gray and white. Right: phase residuals as a function of time for a short interval during episode 3. The red line shows Rs(t), while the blue line is calculated as suggested by explanatory notes of JBE (see text for details). Error bars show JBE quoted arrival-time uncertainties.

In the text
thumbnail Fig. 6

Sinusoidal phase functions Rj(t) (red lines) that fit the JBE phase residuals (blue points) for episodes from j = 1 to 10 (the short episode 5 is included in the graph showing episode 6). Time on the abscissa is in MJD, and the scale is the same for all episodes. Green dots are numbered as in Espinoza et al. (2011) at dates of published glitches with the ordinate as calculated from Rj(t). Ordinate scales are different and adjusted to different amplitudes of oscillations. Each panel shows data of the complete episode and includes the first point of the next episode. The difference between horizontal dotted lines represents a measure of the strength of the perturbation Rj(t) caused by the glitch. Dashed gray curves are the plots of ϕj(t)−ϕj−1(t), to show the difference of phases between two contiguous episodes. Gray arrows point to dates of the > 100 MeV X-ray flares detected by Fermi (Abdo et al. 2011; Ojiha et al. 2013; Bühler et al. 2012) and AGILE (Tavani et al. 2011; Striani et al. 2013).

In the text
thumbnail Fig. 7

Left: dispersion measure DM(t) (blue dots), braking index n(tτ) for τ = 1100 days (continuous red line). The dashed red line shows n(t). Right: correlation between braking index (n) and dispersion measure DM. The correlation coefficient is 0.7.

In the text
thumbnail Fig. A.1

Fractional part of calculated number of turns between successive TOAs’ fractional phase error. Vertical lines denote the dates of reported glitches, starting from glitch No. 5 in Table A.1.Horizontal dashed lines are at ± 0.062, the variance of the (wide winged) phase error distribution.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.