What brakes the Crab pulsar?
^{1} Faculty of Mathematics and Physics, University of Ljubljana, 1000 Ljubljana, Slovenia
email: Andrej.Cadez@fmf.unilj.si
^{2} INAF, Astronomical Observatory of Padova, 35122 Padova, Italy
^{3} Department of Physics and Astronomy, University of Padova, 35122 Padova, Italy
^{4} Department of Information Engineering, University of Padova, 35122 Padova, Italy
^{5} CNRIFN UOS Padova LUXOR, 35131 Padova, Italy
^{6} Department of Physics, University of Atacama, 5 Copiapo, Chile
Received: 7 May 2015
Accepted: 3 October 2015
Context. Optical observations provide convincing evidence that the optical phase of the Crab pulsar follows the radio one closely. Since optical data do not depend on dispersion measure variations, they provide a robust and independent confirmation of the radio timing solution.
Aims. The aim of this paper is to find a global mathematical description of Crab pulsar’s phase as a function of time for the complete set of published Jodrell Bank radio ephemerides (JBE) in the period 1988−2014.
Methods. We apply the mathematical techniques developed for analyzing optical observations to the analysis of JBE. We break the whole period into a series of episodes and express the phase of the pulsar in each episode as the sum of two analytical functions. The first function is the bestfitting local braking index law, and the second function represents small residuals from this law with an amplitude of only a few turns, which rapidly relaxes to the local braking index law.
Results. From our analysis, we demonstrate that the power law index undergoes “instantaneous” changes at the time of observed jumps in rotational frequency (glitches). We find that the phase evolution of the Crab pulsar is dominated by 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. Deviations from such a regular phase description behave as oscillations triggered by glitches and amount to fewer than 40 turns during the above period, in which the pulsar has made more than 2 × 10^{10} turns.
Conclusions. Our analysis does not favor the explanation that glitches are connected to phenomena occurring in the interior of the pulsar. On the contrary, timing irregularities and changes in slow down rate seem to point to electromagnetic interaction of the pulsar with the surrounding environment.
Key words: pulsars: general / pulsars: individual: PSR B0531+21 / radiation mechanisms: general / stars: magnetic field
© 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 powerlaw relation between the frequency of rotation f and its derivative: ḟ = −Kf^{n}.
The braking index 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 longlasting 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. Longterm welldefined 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 spindown.
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 spindown (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 quasiperiodic behavior with phase modulations on typical timescales between ~ 1 and 10 yr (Hobbs et al. 2010). Such quasiperiodic 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 ofswitching of the radio emission and drastic changes in braking torque, was proposed as a prototype of a pulsarmagnetosphere 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 1969−1993 (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 (JBE^{2}). 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 signaltonoise ratio optical observations of the Crab pulsar. The first set of data was obtained in October 2008 with the ultrafast 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 twosecond 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 Gaussiandistributed with the width consistent with photoncounting noise. In this way a typical onehour 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 welldocumented 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 14month 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.
Fig. 1 Residuals of optical and radio phase with respect to a BIM model. Left: bestfitting BIM to 2009 Iqueye data (n = 2.437). Right: bestfitting BIM to 2008−2009 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 yaxis in the insets are very different. 

Open with DEXTER 
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 Xrayradio 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 14month 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 longterm 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 glitchless 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: f_{opt}−f_{radio} = −5.64 × 10^{9} s^{1}, ḟ_{opt}−ḟ_{radio} = 1.55 × 10^{14} s^{2}; and December 15. 2009: f_{opt}−f_{radio} = −6.45 × 10^{9} s^{1}, ḟ_{opt}−ḟ_{radio} = −6.86 × 10^{15} s^{2}.
2.1. Braking index model implementation
The braking index model implies that the phase ϕ can be expressed as (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 is the braking parameter. The frequency and its derivative are then the following functions of time and . The last equation leads to the familiar form of the braking index law: (2)therefore .
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.
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). 

Open with DEXTER 
3. Braking episodes
That the phase history of the Crab pulsar can be approximated by a BIM to within one turn during a 14month 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 JBEpublished frequency derivatives with respect to frequency derivatives predicted by a BIM with (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 26year 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 26year 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.
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. 

Open with DEXTER 
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 spindown 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 thirddegree 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).
Inside each episode we seek a braking index law ϕ_{j}(t) in the same form as in Eq. (1) with s_{j} as an additional fitting parameter. Our goal is to obtain the s_{j} 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) = c_{j} + a_{j}(1 + b_{j}t)^{sj}, valid on complete intervals of episodes , where T^{b} is the starting MJD for each episode.
These local fits are joined into a continuous curve Φ(t) by choosing c_{j + 1} in such a way that . 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 turn^{3}. 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 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 T_{j} 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 as listed in Table 1.
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 R_{j}(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. 

Open with DEXTER 
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 × 10^{10} 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 R_{s}(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 . Only nine glitches (n. 6, 7, 9, 14, 16, 20, 23, 24, 25) show secondderivative beyond 2 × 10^{14} s^{2}, which is smaller than the second derivative of optical phase residuals derived from the January 2009 data ( 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 R_{s}(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, (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.
Fig. 5 Left: second time derivative of the residual phase R_{s}(t), a thirdorder 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 R_{s}(t), while the blue line is calculated as suggested by explanatory notes of JBE (see text for details). Error bars show JBE quoted arrivaltime uncertainties. 

Open with DEXTER 
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: (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 arrivaltime 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.
Fig. 6 Sinusoidal phase functions R_{j}(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 R_{j}(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 R_{j}(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 Xray 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). 

Open with DEXTER 
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.
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. 

Open with DEXTER 
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 afterglitch recovery behavior frequently modeled in the literature through a singleparameter 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 brakinglaw 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 brakinglaw 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 ~ 10^{9} turns executed during an episode. The split between the smooth brakinglawdominated 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 semiperiodic 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 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 pulsarwinddriving 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 gammaray flares is also consistent with such a timescale.
Our ultrafast 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 gammaray 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.
Acknowledgments
This work is based on observations made with ESO Telescopes at the La Silla Paranal Observatory under program IDs 082.D0382 and 084.D0328(A) and on observations collected at the Copernico Telescope (Asiago, Italy) of the INAFOsservatorio 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 INAFAstronomical 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
 Abdo, A. A., Ckermann, M., Ajello, M., et al. 2011, Science, 331, 739 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Ajello, M., Albert, A., Allafort, A., et al. 2014, ApJ, 789, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, P. W., & Itoh, N. 1975, Nature, 256, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Baym, G., Pethick, C., Pines, D., & Ruderman, M. 1969, Nature, 224, 872 [NASA ADS] [CrossRef] [Google Scholar]
 Barbieri, C., Naletto, G., Occhipinti, T., et al. 2009, J. Mod. Opt., 56, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Bühler, R., Scargle, J. D., Blandford, R. D., et al. 2012, ApJ, 749, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2014, ApJ, 782, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Cheng, K. S. 1987, ApJ, 321, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Espinoza, C. M., Lyne, A. G., Stappers, B. W., Kramer, M., et al. 2011, MNRAS, 414, 1679 [NASA ADS] [CrossRef] [Google Scholar]
 Espinoza, C. M., Antonopoulou, D., Stappers, B. W., Watts, A., Lyne, A. G., et al. 2014, MNRAS, 440, 2755 [NASA ADS] [CrossRef] [Google Scholar]
 Germanà, C., Zampieri, L., Barbieri, C., et al. 2012, A&A, 548, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldwire, H. C., & Michel, F. C. 1969, ApJ, 156, L111 [NASA ADS] [CrossRef] [Google Scholar]
 Hester, J. J., Mori, K., Burrows, D., Gallagher, J. S., et al. 2002, ApJ, 577, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Ho, W. C., & Andersson, N. 2012, Nat. Phys., 8, 787 [NASA ADS] [CrossRef] [Google Scholar]
 Hobbs, G., Lyne, A. G., & Kramer, M. 2010, MNRAS, 402, 1027 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Livingstone, M. A., et al. 2007, Ap&SS, 308, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Lyne, A. G., Pritchard, R. S., & GrahamSmith, F. 1993, MNRAS, 265, 1003 [NASA ADS] [CrossRef] [Google Scholar]
 Lyne, A., Hobbs, G., Kramer, M., Stairs, I., & Stappers, B. 2010, Science, 329, 408 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lyne, A. G., Jordan, C. A., GrahamSmith, F., et al. 2015, MNRAS, 446, 857 [NASA ADS] [CrossRef] [Google Scholar]
 Michel, F. C., & Tucker, W. H. 1969, Nature, 223, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Striani, E., Tavani, M., & Ferrari, A. 2013, MNRAS, 436, 1102 [NASA ADS] [CrossRef] [Google Scholar]
 Naletto, G., Barbieri, C., Occhipinti, T., et al. 2009, A&A, 508, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ojha, R., Buehle, R., Hays, E., & Dutka, M. 2013, Atel, 4855 [Google Scholar]
 Oosterbroek, T., de Bruijne, J. H. J., Martin, D., et al. 2006, A&A, 456, 283 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oosterbroek, T., Cognard, I., Golden, A., et al. 2008, A&A, 488, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Porth, O., Komissarov, S. S., & Keppens, R. 2014, MNRAS, 438, 278 [NASA ADS] [CrossRef] [Google Scholar]
 Reichley, P. E., & Downs, G. S. 1969, Nature, 222, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Rots, A. H., Jahoda, K. L., & Andrew, G. 2004, ApJ, 605, L129 [NASA ADS] [CrossRef] [Google Scholar]
 Scott, D. M., Finger, M. H., & Wilson, C. A. 2003, MNRAS, 344, 412 [NASA ADS] [CrossRef] [Google Scholar]
 Striani, E., Tavani, M., Vittorini, V., et al. 2013, ApJ, 765, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Tavani, M., Bulgarelli, A., Vittorini, V., et al. 2011, Science, 331, 736 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Trimble, V. 1968, AJ, 73, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Weisskopf, M. C., Tennant, A. F., Arons, J., et al. 2013, ApJ, 765, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Weyler, K. W., & Panagia, N. 1978, A&A, 70, 419 [NASA ADS] [Google Scholar]
 Wong, T., Backer, D. C., & Lyne, A. G. 2001, ApJ, 548, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Zampieri, L., Verroi, E., Naletto, G., et al., 2013, Atel, 4878 [Google Scholar]
 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 T_{j} = { MJD_{j} + (1 / 86400)TOA_{j},z_{j} }, where MJD_{j} is the mean Julian date of the jth entry of JB ephemerides, TOA_{j} is the first arrival time of the pulse at MJD_{j} (the t_{JPL} entry in JBE), and z_{j} is the integer number of turns the pulsar has made from the starting date May 15, 1988 (MJD_{1} = 47 296). Integers z_{j} are calculated by summing the integer number of turns N_{j} that the pulsar has made between MJD_{j−1} and MJD_{j} (). N_{k} 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 , and let ν_{j} and be the frequency and frequency derivative listed at the jth ephemerides entry in JBE. Define P_{j} = 1 /ν_{j} and , then and N_{j} = IntegerPart[ΔN]. 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.
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. 

Open with DEXTER 
Numbers and dates of glitches (from Espinoza et al. 2011).
All Tables
All Figures
Fig. 1 Residuals of optical and radio phase with respect to a BIM model. Left: bestfitting BIM to 2009 Iqueye data (n = 2.437). Right: bestfitting BIM to 2008−2009 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 yaxis in the insets are very different. 

Open with DEXTER  
In the text 
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). 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 
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 R_{j}(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. 

Open with DEXTER  
In the text 
Fig. 5 Left: second time derivative of the residual phase R_{s}(t), a thirdorder 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 R_{s}(t), while the blue line is calculated as suggested by explanatory notes of JBE (see text for details). Error bars show JBE quoted arrivaltime uncertainties. 

Open with DEXTER  
In the text 
Fig. 6 Sinusoidal phase functions R_{j}(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 R_{j}(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 R_{j}(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 Xray 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). 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 