On the Relaxation Behaviors of Slow and Classical Glitches: Observational Biases and Their Opposite Recovery Trends
aa r X i v : . [ a s t r o - ph . H E ] S e p On the Relaxation Behaviors of Slow and Classical Glitches:Observational Biases and Their Opposite Recovery Trends
Yi Xie , , Shuang-Nan Zhang , ABSTRACT
We study the pulsar timing properties and the data analysis methods dur-ing glitch recoveries. In some cases one first fits the time-of-arrivals (TOAs)to obtain the “time-averaged” frequency ν and its first derivative ˙ ν , and thenfits models to them. However, our simulations show that ν and ˙ ν obtained thisway are systematically biased, unless the time intervals between the nearby datapoints of TOAs are smaller than about 10 s, which is much shorter than typ-ical observation intervals. Alternatively, glitch parameters can be obtained byfitting the phases directly with relatively smaller biases; but the initial recoverytimescale is usually chosen by eyes, which may introduce a strong bias. We alsoconstruct a phenomenological model by assuming a pulsar’s spin-down law of˙ νν − = − H G ( t ) with G ( t ) = 1 + κe − t/τ for a glitch recovery, where H is aconstant and κ and τ are the glitch parameters to be found. This model canreproduce the observed data of slow glitches from B1822–09 and a giant classicalglitch of B2334+61, with κ < κ >
0, respectively. We then use this model tosimulate TOA data and test several fitting procedures for a glitch recovery. Thebest procedure is: 1) use a very high order polynomial (e.g. to 50th order) toprecisely describe the phase; 2) then obtain ν ( t ) and ˙ ν ( t ) from the polynomial;and 3) the glitch parameters are obtained from ν ( t ) or ˙ ν ( t ). Finally, the uncer-tainty in the starting time t of a classical glitch causes uncertainties to someglitch parameters, but less so to a slow glitch and t of which can be determinedfrom data. Subject headings: stars: neutron - pulsars - individuals: B1822–09, B2334+61 -general - magnetic fields National Astronomical Observatories, Chinese Academy of Sciences, Beijing, 100012, China University of Chinese Academy of Sciences, Beijing, 100049, China Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,Beijing 100049, China; [email protected]
1. Introduction
Pulsars are very stable rotators. However, many pulsars exhibit significant timing ir-regularities, i.e., unpredicted arrival times of pulses. There are two main types of timingirregularities, namely ‘timing noise’ which is consisted of low-frequency quasi-periodic struc-tures, and ‘glitches’ which are abrupt increases in their spin rates followed by relaxations.Glitch activities are more frequent in relatively young pulsars with a characteristic ageof 10 − yr (Shemar & Lyne 1996; Wang et al. 2000). For the hundreds of glitchesobserved , their typical fractional jumps in spin frequency ν are in the range of ∆ ν/ν ≈ − − − , and their relative increment in frequency derivative is ∆ ˙ ν/ ˙ ν ∼ − . Despitethe abundance of observational data accumulated for over 40 years, we are still far fromsatisfactory understanding of glitch events. Traditional models mainly involve the expectedsuperfluid nature of part of the neutron star interior (Anderson & Itoh 1975; Ruderman1976), and the angular momentum is carried in the form of microscopic, quantized vortices,whose density determines the rotation rate of a pulsar. Mostly, these vortices are pinnedto the crust and the charged matter in the core of the star, thus their outward driftingmotions are prevented (Anderson & Itoh 1975; Alpar 1977; Pines et al. 1980; Alpar et al.1981; Anderson et al. 1982). However, as the crust spins down due to the electromagneticbraking, a rotational lag and stress (Magnus force) gradually builds up. A glitch occurswhen the stress reaches some critical value and the pinning breaks, vortices suddenly moveoutward and impart their angular momentum to the crust. Immediately after the glitch, thevortices are pinned to other parts again and the superfluid is effectively decoupled from thecrust.Following the seminal work of Baym, Pethick & Pines (1969), there are two classes ofmodels that have been developed to explore the dynamical evolution of pinned superfluidduring the post-glitch recovery. One kind of models involve a weak coupling between thesuperfluid and the crust due to the interaction between free vortices and the coulomb lattice ofnuclei (Jones 1990, 1992, 1998). Another kind of models assume that the vortices creep rate ishighly temperature-dependent. As the vortices creep through the crust, angular momentumis gradually transferred (Alpar 1984a, 1984b; Link, Epstein & Baym 1993; Larson & Link2002). Superfluid vortex dynamics can model the relaxation well; however, there are stillmany significant problems unsolved. For instance, the mechanism that triggers the glitch inthe first place and the detailed processes of angular momentum transfer during the recoveryare still controversial. It has been suggested that such an event may be triggered by largetemperature perturbations (Link & Epstein 1996), or caused by starquakes (Baym & Pines ν ,accompanied by a rapid decrease in | ˙ ν | and subsequent exponential increase back to its initialvalue (Zou et al. 2004; Shabanova 2005). That is the so-called ‘slow glitch’. Currently, thereis still no convincing theoretical understanding for slow glitches. Peng & Xu (2008) proposedthat, after a collapse or a small star-quake, the solid superficial layer of a rigid quark starmay be heated and becomes a viscous fluid, which will eventually produce a gradual increasein ν . However, Hobbs et al. (2010) and Lyne et al. (2010) argued that the slow glitcheshave the same origin as the timing noise of many pulsars.For the recovery processes of both glitch and slow glitch events, the variations of spin 4 –frequency ν and its first derivative ˙ ν of pulsars are obtained from polynomial fit results ofarriving time epochs of pulses. The local TOAs of the mean pulses for individual observingsessions are determined from the maximum cross correlation between the observed meanpulses and a Gaussian profile template. The profile template is a mean pulse with high signal-to-noise ratio, obtained by summing the best-quality mean pulses over several observingsessions. Correction of TOAs to the solar system barycenter can be done using TEMPO2 program with the Jet Propulsion Laboratory DE405 ephemeris (Standish 1998). These TOAsare then weighted by the inverse squares of their estimated uncertainty. Since the rotationalperiod is nearly constant, these observable quantities, ν , ˙ ν and ¨ ν can be obtained by fittingthe phases to the third order of its Taylor expansion over a time span t s ,Φ i = Φ + ν ( t i − t ) + 12 ˙ ν ( t i − t ) + 16 ¨ ν ( t i − t ) . (1)One can thus get the values of ν , ˙ ν and ¨ ν at t from fitting to Equation (1) for independent N data blocks around t , i.e. i = 1 , ..., N . Apparently, these observational quantities obtainedthis way are not instantaneous results, rather, the “averaged” results over each data block (i.e.over each t s ) and extrapolated to t , which are not necessarily the same as the instantaneousvalues (denoted as ν I and ˙ ν I ). Thus, they are called “averaged” values (denoted as ν A and˙ ν A ) in this work. Usually, t s is much less than pulsar’s spin-down age τ c , thus the differencesbetween instantaneous values and “averaged” values are not significant, and consequently ν A and ˙ ν A are good approximations for ν I and ˙ ν I in most cases. However, it has been foundrecently that oscillations of the “apparent” magnetic fields of neutron stars are responsiblefor the observed signs and magnitudes of ¨ ν , the second derivative of frequency, and brakingindices (Biryukov et al. 2012; Pons et al. 2012; Zhang & Xie 2012a, 2012b). We furthersuggested that the oscillation time scales are between 10-100 yr, comparable to t s , thusmaking the fitted spin-down parameters different from the true and instantaneous spin-downparameters. Similarly, considerable biases may also exist when fitting the glitch recoverydata, since the glitch recovery time scales are also comparable with t s .In section 2, we simulate several pulsar timing data analysis procedures for glitch re-coveries, and find that the glitch parameters, obtained from the averaged ν A and ˙ ν A , havesignificant systematic biases compared with that obtained with the instantaneous ν I and ˙ ν I .In order to get the true glitch parameters with the reported, yet averaged glitch recoverydata ν O and ˙ ν O , a phenomenological or physical glitch model is needed to be combined withsimulations. We thus present a phenomenological spin-down model during a glitch recovery,and model several slow glitch recovery events and the recovery of a giant classical glitch ν ( t ) and ˙ ν ( t ). In Section 5, we discusshow to obtain the model parameters of glitch recoveries more accurately. The results aresummarized in section 6.
2. Simulating Data Analysis of Glitch Recoveries2.1. Simulation for ˙ ν -fitting procedure By fitting the TOA set { Φ( t i ) } to Equation (1), one can get { ν ( t ) } and { ˙ ν ( t ) } . When { ν ( t ) } and { ˙ ν ( t ) } show exponential relaxations, their variations following the jump at epoch t can be described as the following empirical functions (e.g. Yuan et al. 2010, Roy et al.2012), ν ( t ) = ν ( t ) + ∆ ν p + ∆ ˙ ν p ∆ t + 12 ∆¨ ν p ∆ t + X j ∆ ν d j e − ∆ t/τ j , (2)and ˙ ν ( t ) = ˙ ν ( t ) + ∆ ˙ ν p + ∆¨ ν p ∆ t + X j ∆ ˙ ν d j e − ∆ t/τ j , (3)where ∆ t = t − t , ∆ ν p and ∆ ˙ ν p are permanent changes in ν and ˙ ν relative to the pre-glitchsolution ν ( t ) and ˙ ν ( t ), ∆ ν d j is the amplitude of the j th decaying component with a timeconstant τ j , and ∆ ˙ ν d j = − ∆ ν d j /τ j . One can get the glitch parameters ∆ ν p , ∆ ˙ ν p , ∆¨ ν p ,∆ ν d j , ∆ ˙ ν d j and τ by fitting ν ( t ) and ˙ ν ( t ) to Equations (2) and (3), respectively. The twofunctions describe the post-glitch behaviors fairly well, especially for the case of a long termrecovery, and usually multiple decay terms with different decay time constants can be fitted(e.g. there are up to five exponentials are fitted for Vela 2000 and 2004 glitches; Dodson etal. 2002, 2007). For simplicity, the cases that ν varies as one exponential decay term or twoexponential decay terms are assumed in the following simulations. Slow glitches are characterized by a gradual increase in ν with a long timescale of several months, accompanied by a rapid decrease in | ˙ ν | by a few percent,which is sometimes even shorter than the observation interval and thus cannotbe seen. Then | ˙ ν | experiences an exponential increase back to its initial valuewith the same time scale as that of ν increase (Shabanova 2005). Analogous tothe classical glitches, we suggest that the slow glitches can be described by the following twofunctions: ν ( t ) = ν ( t ) + ∆ ν p + ∆ ˙ ν p ∆ t + 12 ∆¨ ν p ∆ t + X j ∆ ν d j (1 − e − ∆ t/τ j ) , (4) 6 –and ˙ ν ( t ) = ˙ ν ( t ) + ∆ ˙ ν p + ∆¨ ν p ∆ t + X j ( − ∆ ˙ ν d j ) e − ∆ t/τ j , (5)where the parameters are the same as those in Equations (2) and (3). Since the glitch or slow glitch recoveries can be described by Equations (2)-(5), somesimple models can also be derived from them. For a classical glitch, we simply assume ν ( t ) = ∆ ν d exp ( − ∆ t/τ ) , (6)i.e., ν = ∆ ν p = ∆ ˙ ν p = ∆¨ ν p = 0. We will use this equation to produce simulated data,and obtain the “instantaneous” ν I = ν ( t ) and ˙ ν I = dν ( t ) /dt , with the parameters (∆ ν d and τ ) given later. On the other hand, the “averaged” values are obtained by the followingprocedure. Firstly, we get the phase by Φ( t ) = R tt ν ( t ′ ) dt ′ . For convenience we take t = 0.However, in practice t cannot be known precisely due to discontinuous observations; we willshow later that this will cause some uncertainty in estimating the parameters of a classicalglitch, but not so for slow glitches. We assume a certain time interval ∆ T int between eachtwo nearby TOAs, i.e. ∆ T int ≡ t i +1 − t i is a constant. We set ten adjacent TOAs in oneblock (i.e. N = 10 in Equation (1)), and the latter five TOAs are used as the first five TOAsin the next block. We then fit the TOA blocks to Equation (1) to obtain ν A and ˙ ν A , whichare the fitted coefficients of ν term and ˙ ν term of the equation, respectively. The time t for ν A and ˙ ν A is taken as the middle epoch of each block, i.e., t = ( t + t ) /
2, and is also“averaged” (e.g. Yuan et al. 2010).In Figure 1, we show these instantaneous values and averaged values with different ∆ T int for a glitch with ∆ ν d = 0 . µ Hz and τ = 50 days. One can see that both ν A and ˙ ν A haveremarkably different decay profiles from ν I and ˙ ν I during the recovery process, respectively.This systematic biases are independent of ∆ T int , and it seems that the recovery time-scale τ is the key parameter that is mainly biased. By fitting ν A and ˙ ν A to Equation (2) andEquation (3), respectively, we find that all the recovery time scales of ν A and ˙ ν A are muchlonger than the time scale of 50 days (e.g. τ ≈
95 day for ∆ T int = 10 s). The systematicdifferences between the decay profiles of ν A or ˙ ν A and the profile of ν I or ˙ ν I are considerable,and apparently caused by the procedure of fitting TOAs to Equation (1); thus for higherorder fits, one cannot consider the first order coefficient to be the “frequency”. This procedureis thus abandoned for glitch data analysis in the following. -0.020.000.020.040.060.080.100.12 ( H z ) ( - H z s - ) Instantaneous
Averaged: s Averaged: s Averaged: s Instantaneous
Averaged: s Averaged: s Averaged: s Time (days) . Fig. 1.— Variations of ν (upper panel) and ˙ ν (bottom panel) for a simulated classical glitchwith one decay term . Solid lines represent their instantaneous values; circles, triangles, andcrosses represent the averaged values obtained by fitting the simulated TOAs to Equation (1)for the cases of ∆ T int = 10 , 10 , and 10 s, respectively. 8 –However with the TEMPO2 software, ν may be obtained from the TOAs by fitting toΦ i = Φ + ν ( t i − t ) , (7)and ˙ ν may be obtained by fitting toΦ i = Φ + ν ( t i − t ) + 12 ˙ ν ( t i − t ) , (8)i.e. the first two or three terms of Equations (1), respectively (Yu 2013). Here, we first fitthe TOA blocks to Equation (7) to obtain ν A , which is the fitted coefficients of ν term ofEquation (7). We then separately fit the TOA blocks to Equation (8) to obtain ˙ ν A , which isthe fitted coefficients of ˙ ν term of Equation (8). In the left panels of Figure 2, we show theinstantaneous and averaged values obtained this way, with different ∆ T int for a glitch withthe same ∆ ν d and τ . Clearly now the profiles of both ν A and ˙ ν A follow that of ν I and ˙ ν I without obvious distortions. By fitting to Equation (2) or Equation (3), we find that all therecovery time scales of ν A or ˙ ν A equal the time scale of 50 days, i.e. τ has not been biased.We can then obtain the normally reported glitch parameters ∆ ν d and ∆ ˙ ν d , as listed inTable 1, by fitting { ν A } or { ˙ ν A } to Equation (2) or Equation (3) with different ∆ T int ; forcomparison we also list ∆ ν d and ∆ ˙ ν d obtained from ν I and ˙ ν I . One can see that “averaged”∆ ν d or ∆ ˙ ν d (denoted as ∆ ν dA or ∆ ˙ ν dA hereafter) have systematic differences from instan-taneous ∆ ν d or ∆ ˙ ν d (denoted as ∆ ν dI or ∆ ˙ ν dI hereafter). For ∆ T int = 10 s, the differencesare tiny and the glitch parameters can be restored satisfactorily; however, for ∆ T int > s,both the “averaged” ∆ ν d and ∆ ˙ ν d may be considerably smaller than the instantaneous ∆ ν d and ∆ ˙ ν d , respectively.For a slow glitch, we assume ν ( t ) = ∆ ν d (1 − exp ( − ∆ t/τ )) , (9)where t = 0 and ∆ ν d = 0 . µ Hz and τ = 50 days. We show the averaged glitch parametersand profiles with different ∆ T int in Table 1 and the right panels of Figure 2, as well as thoseinstantaneous ones. One can see that we always have ∆ ν dA = ∆ ν dI for any ∆ T int , since ∆ ν dA is determined by the differences of | ˙ ν A | between the data points slightly beforethe starting point of the glitch and the data points at the end of the recovery,and both of them are always available for slow glitch observations. However, ∆ ˙ ν dA is biased in the same way as for the simulated classical glitch. We simply assume ν ( t ) = ∆ ν d1 exp ( − t/τ ) + ∆ ν d2 exp ( − t/τ ) for a classical glitchwith two decay terms, where ∆ ν d1 = 0 . µ Hz, τ = 21 . ν d2 = 0 . µ Hz, 9 – -0.020.000.020.040.060.080.100.12 ( H z ) -0.020.000.020.040.060.080.100.12 Averaged: s Averaged: s Averaged: s Instantaneous
Averaged: s Averaged: s Averaged: s Instantaneous
Averaged: s Averaged: s Averaged: s Instantaneous
Averaged: s Averaged: s Averaged: s Instantaneous ( - H z s - ) -24-20-16-12-8-40 Time (days) . Time (days)
Fig. 2.— Variations of ν and ˙ ν for a simulated classical glitch (left panels) and a simulatedslow glitch (right panels) with one decay term . Solid lines represent their instantaneousvalues; circles, triangles, and crosses represent the averaged values obtained by fitting thesimulated TOAs to Equations (7) (to get ν A ) and (8) (to get ˙ ν A ) for the cases of ∆ T int = 10 ,10 , and 10 s, respectively. 10 – τ = 147 days (the parameters are adopted from pulsar B2334+61 for its very large glitchbetween 2005 August 26 and September 8, Yuan et al. 2010). We also assume ν ( t ) =∆ ν d1 (1 − exp ( − t/τ )) + ∆ ν d2 (1 − exp ( − t/τ )) for a slow glitch with two decay terms. Theinstantaneous values and averaged values are obtained with the same methods describedabove and the main results are presented in Table 1 and Figure 3, in which the similarresults can be found with the case of one decay term. For ∆ T int = 10 s, the differencesfor all of τ , ∆ ν d and ∆ ˙ ν d are tiny and the glitch parameters can be restored satisfactorily.However, things for two decay terms are a little more complicated. For ∆ T int & s,though the data points still converged to the instantaneous values as shown in Figure 3 (i.e.variation trends are the same), the fitted glitch parameters (including τ ) for each componentsare still somewhat biased, and it seems that larger ∆ T int corresponds to a smaller τ for shorttime-scale component. The biases are probably due to the fact that the data are too sparsefor ∆ T int . Actually if τ of the short term decay component is comparable to orshorter than the interval between the observations, then τ of this componentwould be difficult to determine and can only be set as the internal. Similar resultscan be found for a slow glitch, but P ∆ ν dI = P ∆ ν dA is always kept.The above simulations unveil significant biases caused by the averaging procedures (i.e.fitting to Equation (7) and Equation (8)) for ν and ˙ ν during glitch recoveries. Thus, ν A and˙ ν A obtained this way and ν O and ˙ ν O (the subscript “o” means observed values) reported inliterature should not be used directly to test physical models. It should be noted that, forone-decay-term cases, the reported amplitudes of ∆ ν and ∆ ˙ ν of a classical glitch are usuallyunderestimated; the reported amplitude of ∆ ˙ ν of a slow glitch is also underestimated, butthat ∆ ν is not. However, these biases were never noticed in almost all previous theoreticalworks modeling glitch recoveries, and { ˙ ν ( t ) } are usually directly modeled, e.g. the post-glitch fits for Vela pulsar, Crab pulsar and PSR 0525+21 with vortex creep model (Alpar etal. 1984b; Alpar, Nandkumar, & Pines 1985; Alpar et al. 1993; Chau et al. 1993; Alpar, etal. 1996; Larson & Link 2002), and the two-component hydrodynamic model for Vela (vanEysden & Melatos 2010). In these works, the observed data { ˙ ν O ( t ) } (i.e. { ˙ ν A ( t ) } ) are shownin ˙ ν - t diagram and fitted directly by theoretical models. In order to make optimum use of all available data (Shemar & Lyne 1996), the pulsephase induced by a glitch is usually fitted to the following equation, which can give τ and∆ ν d (e.g. Yu et al. 2013): 11 – ( H z ) Averaged: s Averaged: s Averaged: s Instantaneous
Averaged: s Averaged: s Averaged: s Instantaneous -120-100-80-60-40-200 Averaged: s Averaged: s Averaged: s Instantaneous
Averaged: s Averaged: s Averaged: s Instantaneous
Time (days) . Time (days) ( - H z s - ) Fig. 3.— Variations of ν and ˙ ν for a simulated classical glitch (left panels) and a simulatedslow glitch (right panels) with two decay terms . Solid lines represent their instantaneousvalues; circles, triangles, and crosses represent the averaged values obtained by fitting thesimulated TOAs to Equations (7) (to get ν A ) and (8) (to get ˙ ν A ) for the cases of ∆ T int =1 . × , 10 , and 10 s, respectively. 12 – φ g = ∆ φ + ν p ∆ t + 12 ˙ ν p ∆ t + [1 − e − (∆ t/τ i ) ]∆ ν d i τ i , (10)where the ν p and ˙ ν p are the permanent increments in ν and ˙ ν , respectively. However, it isdifficult to get τ i directly by fitting to Equation (10). Actually, TEMPO2 implements onlya linear fitting algorithm, and one thus needs to have a good initial estimate for τ i , whichis estimated from post-glitch ˙ ν variation by eye inspecting. Then the estimated value wasintroduced into Equation (10) fits. By increasing or decreasing τ i , a best estimated τ i canbe eventually found via minimum post-fit χ (Yu et al. 2013). This procedure is widely usedfor classical glitches, but not applied to slow glitches.
We simulate the fitting procedure of Equation (10) as described above, and find thatboth τ and ∆ ν d can be obtained with high precision for one-decay-term case, if a goodinitial estimate for τ is taken, as shown in Table 2. For two-decay-term case, we also assume ν ( t ) = ∆ ν d1 exp ( − t/τ ) + ∆ ν d2 exp ( − t/τ ), where ∆ ν d1 = 0 . µ Hz, τ = 147 days and∆ ν d2 = 0 . µ Hz, τ = 21 . t ) = R ν ( t ) dt . Firstly, we estimate τ for the long term one, and get the best-fit τ by fitting { Φ( t i ) } to Equation (10). Thenwe fix τ and get the timescale of the short term τ the same way. This process is widelyadopted in the data analysis of glitch recovery with TEMPO2 software (Yu et al. 2003).However, we find that the glitch parameters of the long decay term in the two term cases,i.e τ and ∆ ν d1 obtained by this way are already biased, as shown in Table 2.These biases are probably caused by the procedure that fitting the long-decay-term andthe short-decay-term in different steps; the short-decay-term may slightly interfere the firstfitting for τ and ∆ ν d1 of the long-decay-term, thus the results are biased. If the biased τ is fixed, one will also get a biased τ to fit { Φ( t i ) } again to Equation (10), since a localminimum χ will obtained, as shown in Figure 4. Therefore, we suggest that the two termsshould be fitted simultaneously .If ∆ T int . s, the simultaneous fits can be realized by the following steps:(i) Get { ˙ ν } series by fitting { Φ( t i ) } to Equation (8);(ii) Estimate τ and τ by fitting { ˙ ν } to Equation (3) (the calculation cost needed bythis fit is much lower than fitting to Equation (10));(iii) Use the estimated τ and τ as initial values and fit { Φ( t i ) } again to Equation (10),then the best fitted τ i and ∆ ν d i will be obtained.The results of the simultaneous fit are ∆ ν d1 = 0 . µ Hz, τ = 147 . ν d2 =0 . µ Hz, τ = 21 . T int . We also simulate the fitting processes with different 13 –Table 1. ∆ ν d ( µ Hz), ∆ ˙ ν d (10 − Hz s − ) and τ (days) for classical (upper part) and slow(bottom part) glitch simulations. The data in left part and right part are results from onedecay term and two decay term simulations, respectively. “Instantaneous” represents theinstantaneous values, and the different values of ∆ T represent the time intervals betweeneach TOA for the “averaged” values. ∆ ν d ∆ ˙ ν d τ ∆ ν d ∆ ˙ ν d τ Instantaneous 0.100 -23.15 50.00 (0.19, 0.119) (-102.8, -9.4) (21.4, 147.0)∆ T int = 10 s 0.099 -22.88 50.00 (0.18, 0.119) (-100.2, -9.4) (21.3, 145.8)∆ T int = 10 s 0.089 -20.67 50.00 (0.13, 0.139) (-100.4, -16.8) (14.5, 95.6)∆ T int = 10 s 0.041 -9.53 49.52 (2.73, 0.081) (-2960.9, -6.4) (10.7, 146.8)Instantaneous 0.100 23.15 50.00 (0.19, 0.119) (102.8, 9.4) (21.4, 147.0)∆ T int = 10 s 0.100 23.04 49.95 (0.19, 0.120) (100.2, 9.4) (21.3, 145.8)∆ T int = 10 s 0.100 20.67 49.99 (0.19, 0.122) (100.4, 16.8) (14.5, 95.7)∆ T int = 10 s 0.100 9.88 49.36 (0.25, 0.055) (2959.0, 6.4) (10.7, 146.8) Table 2. ∆ ν d ( µ Hz), τ (days) for classical glitch simulations with direct phase-fitprocedure. The data in left part and right part are the results from one decay term andtwo decay term simulations, respectively. Subscript “1” indicates the long decay term intwo term cases. ∆ ν d τ ∆ ν d1 τ Instantaneous 0.100 50.00 0.119 147.0∆ T int = 10 s 0.109 50.00 0.150 128.0∆ T int = 10 s 0.108 50.00 0.151 127.6∆ T int = 10 s 0.106 50.00 0.158 124.0
14 –
18 19 20 21 22 23 240.000.040.080.120.160.20 (days) =145.0 days =145.5 days =146.0 days =146.5 days =147.0 days =147.5 days Fig. 4.— The local minimum of χ produced by fitting { Φ( t i ) } to Equation (10). Duringeach fitting, τ is fixed to a certain value around 147 days. 15 –values of τ i and ∆ ν d i , and all the glitch parameters are restored with relatively small biases,some of which are even better than the previous procedures for ∆ T int . s. Here, wewant emphasize that the fitting procedures described in literature are in chaos. Many authorsadopted the procedure of fitting the TOAs to Equation (10) to obtain the pulsars parameters,and using Equation (8) to get { ˙ ν } , but only Equation (1) is mentioned in their papers (Yu2013). Thus, we suggest that the exact fitting procedure should be described in detail.
3. A Phenomenological Spin-down Model
In this section, we develop a phenomenological spin-down model to describe the glitchand slow glitch recoveries, so that the model can be a tool to simulate data to test the dataanalysis procedures for the recoveries in the next section. Classically, a magnetic dipolewith a magnetic moment M = BR , rotating in vacuum with angular velocity Ω, emitselectromagnetic radiation with a total power 2 M Ω / c . Assuming the pure magneticdipole radiation as the braking mechanism for a pulsar’s spin-down, the energy loss rate isthen given by ˙ E = I Ω ˙Ω = − BR sin χ ) c Ω , (11)where B is its dipole magnetic field at its magnetic pole, R is its radius, I is its moment ofinertia. Equation (11) is modified slightly in order to describe a glitch event,˙ΩΩ − = − BR sin χ ) c I G ( t ) , (12)in which G ( t ) represents very small changes in the effective strength of dipole magneticfield B sin χ , or the effective moment of inertia I of both the pulsar and its magnetosphereduring a glitch recovery. The left hand are observable quantities, and the right hand are alltheoretical quantities. In the following we assume G ( t ) = 1 + κe − ∆ t/τ . Then Equation (12)can be written as ˙ νν − = − H (1 + κe − ∆ t/τ ) , (13)where ˙ ν = ˙Ω / π and H = π ( BR sin χ ) c I = 1 / τ c ν , and τ c = − ν/ ν is the characteristic ageof a pulsar.Integrating and solving Equation (13), we have ν ( t ) = ν p t + κτ (1 − e − (∆ t ) /τ )) /τ c . (14)The derivative of ν is ˙ ν ( t ) = − (1 + κe − ∆ t/τ ) ν / τ c (1 + (∆ t + κτ (1 − e − ∆ t/τ )) /τ c ) / . (15) 16 –We know ∆ t ∼ τ ∼
100 days and generally τ c & years, and the term | (∆ t + κτ (1 − e − ∆ t/τ )) /τ c | ≪ κ ≪
1, the expression of ν and ˙ ν can be approximately written inthe same forms of Equations (2) and (3), which give ∆ ν d = ν κτ / τ c and ∆ ˙ ν d = − ν κ/ τ c .Numerical calculations show that Equations (2) and (3) with these parameters give identicalresults as Equations (14) and (15) for all known ranges of glitch parameters. The expressionof ∆ ν p and ∆ ˙ ν p that relate to the initial jumps of ν and ˙ ν , are not given by the model,since the glitch relaxation processes are only considered here. It has been suggested thatthese non-recoverable jumps are the consequence of permanent dipole magnetic field increaseduring the glitch event (Lin and Zhang 2004). ∆¨ ν p is the jump of timing residual, which isbeyond the scope of the present work. In the following we attempt to apply this phenomenological model to fit thereported data of several slow glitches of B1822–09 and one classical glitch ofB2334+61. Since the reported data points of ν and ˙ ν are too sparse (about onepoint per 150 days for B1822–09 or 50 days for B2334+61) and the TOAs ofthese glitches are not available in literature, we cannot apply our model to fitthe reported to obtain both τ and κ simultaneously, as that done in the abovesimulations. As a compromise, we focus only on determining κ by applyingour phenomenological model and simply take (the inevitably biased) τ obtainedby fitting directly the reported ν and ˙ ν . Therefore κ remains the only glitchrecovery parameter to be determined from observations in the following. Ourmain purpose here is to show the applicability of the our phenomenological modelto describe glitch observations.3.1. Modeling several slow glitches of B1822–09 We first model slow glitches because they are simpler than the classical glitches, espe-cially they have no jumps in ν and ˙ ν , i.e. ∆ ν p = 0 and ∆ ˙ ν p = 0. Shabanova (2005) reportedthree slow glitches of B1822–09 (J1825-0935) over the 1995-2004 interval. The pulsar has ν ≃ .
30 s − , a relatively large | ˙ ν | ≃ . × − s − (note ˙ ν < τ c ≃
232 kyrand B ≃ . × G. As shown in Figure 5, the pulsar experienced three slow glitchesfrom 1995 to 2005. A gradual increase in ν is well modeled by an exponential function withtimescales of 235, 80 and 110 days, respectively. For ˙ ν , the fractional decreases of | ˙ ν | (i.e.increases of ˙ ν ) are about 0 .
7, 2 . . | ˙ ν | (i.e. decreases of ˙ ν ) back to the previous values with the same time scales are also welldescribed by exponential functions. The third slow glitch was separately detected by Zou etal. (2004). 17 – SimulatedSimulated
Fig. 5.— Slow glitches of Pulsar B1822–09. Observational results are taken from Shabanova(2005). Upper panels: variations of ∆ ν relative to the pre-glitch solution. Bottom panels:variations of ˙ ν . Left panels: comparison between the reported and simulated (both are alsotime-averaged) ∆ ν and ˙ ν . Right panels: comparison between the reported and restored (i.e.model-predicted) instantaneous ∆ ν and ˙ ν . 18 –Since the detailed data on ∆ T int are not reported in literature, we assume an uniformTOA distribution with ∆ T int = 3 . × s. We take the following steps in modeling theobserved data for each slow glitch event:(i) We get our model-predicted TOAs with ∆ T int by integrating Equations (2) or (14),with a κ for each slow glitch event.(ii) We simulate the data analysis process by fitting every block of ten adjacent TOAsto Equations (7) or Equations (8) to obtain one set of ν A and ˙ ν A ; and the latter five TOAsare also used in the next TOA block.(iii) The above simulated ν A and ˙ ν A are compared with the reported glitch profile ν O and ˙ ν O ; κ is adjusted until reasonable agreements between them are reached.With the above steps, we confirm that the slow glitch behavior can be explained byour phenomenological model with κ <
0. Our modeling results are shown in Figure 5. Thefit parameter κ is − . − .
06 and − .
04 for the three slow glitch events, respectively.In Table 3 we show the relative magnitudes of ∆ ν and ∆ ˙ ν for the three slow glitches; forcomparison we also list in Table 3 the results for the giant classical glitch from B2334+61obtained in the next section. It is found that the relative magnitudes of ∆ ν A , ∆ ν O and ∆ ν I are identical, i.e. ∆ ν I = ∆ ν A = ∆ ν O , as expected from the above simulations. It is also clearthat the instantaneous values of ∆ ˙ ν , which are calculated directly from the modelwith the parameters determined above, are much larger than the reported results inliterature, e.g. the ∆ ˙ ν I are larger than two times of ∆ ˙ ν O for the second and third slowglitches.Table 3. The relative values of ∆ ν and ∆ ˙ ν for slow glitches of B1822–09 and the giantclassical glitch from B2334+61. In the classical glitch part, the superscripts ‘i’ and ‘ii’represent for results of one-term fit and two-term fit, respectively. Slow Glitches of B1822–09 Classical Glitch of B2334+61∆ ν /ν ∆ ˙ ν / ˙ ν ∆ ν /ν ∆ ˙ ν / ˙ ν ∆ ν /ν ∆ ˙ ν / ˙ ν ∆ ν i /ν ∆ ˙ ν i / ˙ ν ∆ ν ii /ν ∆ ˙ ν ii / ˙ ν (10 − ) (%) (10 − ) (%) (10 − ) (%) (10 − ) (%) (10 − ) (%)Reported 12.9 0.7 28.6 2.7 25.2 1.7 75.8 -2.96 75.8 -2.96Simulated 13.2 0.7 29.7 2.7 25.4 2.0 35.6 -3.15 54.5 -2.87Instantaneous 13.2 0.94 29.7 6.0 25.4 4.0 64.4 -3.85 79.8 -3.98
19 –
The pulsar PSR B2334+61 (PSR J2337+6151) was discovered in the Princeton-NRAOsurvey using the 92 m radio telescope at Green Bank in 1985 (Dewey et al. 1985). It has ν ≃ .
019 s − , ˙ ν ≃ − . × − s − , τ c ≃ . × yr, and B ≃ . × G. It islocated very close to the center of the supernova remnant G114.3+0.3. Yuan et al. (2010)reported the timing observations of PSR B2334+61 for seven years with the Nanshan 25 mtelescope at Urumqi Observatory. A very large glitch occurred between 2005 August 26and September 8 (MJDs 53608 and 53621), the largest known glitch ever observed, with afractional frequency increase of ∆ ν/ν ∼ . × − . Yuan et al. (2010) obtained each ν , ˙ ν and ¨ ν by fitting ten adjacent TOAs to Equation (1), and the latter five TOAs had also beenused as the first five TOAs in the next fit. The rotational behavior during this glitch event isshown in Figure 6. A large jump of rotational frequency could be seen in the top panel with∆ ν ≈ × − Hz. The bottom panel shows a very significant long-term increase in | ˙ ν | afterthe time of jump, and the corresponding braking indices are 10 . ± . . ± . ∼ . ∼
147 days (Yuan et al. 2010).We follow almost the same steps as for the slow glitches above to model the reported,yet time-averaged glitch recovery data of this classical glitch, with ∆ T int = 1 . × s; theonly difference is that a slope of ∆ ˙ ν p = − . × − s − is taken in Equation (2), followingLyne et al. (2000).In the left panels of Figure 6, we show the fits with one exponential term G ( t ) =(1 + κ exp ( − ∆ t/τ )) for a comparison with the “realistic” simulation of two termsbelow . The best parameters for this glitch event are κ = 0 .
038 and τ = 50 days. We thenmodel the glitch recovery process with G ( t ) = (1 + κ exp ( − ∆ t/τ ) + κ exp ( − ∆ t/τ )), asshown in the right panels of Figure 6. The best parameters for this glitch event are κ = 0 . κ = 0 .
012 ( τ = 21 . τ = 147 days are fixed by the observed values). Table3 gives the relative magnitudes of ∆ ν and ∆ ˙ ν for both fits. In order to distinguish betweenthe two fits, we show them in the logarithmic coordinates in Figure 7. Clearly the simulatedprofiles of the two term fit match the reported ones better than that of the one term fit. Onecan see that | ∆ ˙ ν I | are also slightly larger than the reported | ∆ ˙ ν O | for both the one-term fitand two-term fit.In Figure 8 we show ∆ ν with the slope of ∆ ˙ ν p removed, and ˙ ν O − ˙ ν I . It is clearlyshown that one exponential term cannot fit the observed data at the end of decay profile,and this is also the reason why ∆ ν I is smaller than ∆ ν O for this fit, as given in Table 3.Thus, the one-term decay is ruled out , and we focus on the two-term fit below. Using 20 – Simulated Simulated
Fig. 6.— The giant glitch of Pulsar B2334+61. Observational results are taken from Yuanet al. 2010. ∆ T int = 1 . × s is adopted in the fit. Upper panels: variations of ∆ ν .Bottom panels: variations of ˙ ν . The left and right panels represent for models with one andtwo decay components, respectively. 21 – Simulated Simulated
Fig. 7.— The giant glitch of Pulsar B2334+61. Observational results are taken from Yuan etal. 2010. ∆ T int = 1 . × s is adopted in the fit. Upper panels: variations of ∆ ν . Bottompanels: variations of ˙ ν . The left and right panels represent for models with one and twodecay components, respectively. This figure is the same with Fig. 6 except for a logarithmicabscissa. 22 – Simulated Simulated
Fig. 8.— The giant glitch of Pulsar B2334+61. Upper panels: variations of ∆ ν from whichthe slope ∆ ˙ ν p is removed. Bottom panels: residuals of ˙ ν O − ˙ ν I . The left and right panelsrepresent for models with one and two decay components, respectively. 23 – ν , τ c and the determined κ and κ , we obtain these glitch parameters: ∆ ν d1 = 0 . µ Hz,∆ ν d2 = 0 . µ Hz, ∆ ˙ ν d1 = − . × − s − , ∆ ˙ ν d2 = − . × − s − ; some of themhave significant differences from the reported results of Yuan et al. (2010), in which ∆ ν d1 =0 . µ Hz, ∆ ν d2 = 0 . µ Hz, ∆ ˙ ν d1 = − . × − s − , ∆ ˙ ν d2 = − . × − s − .In Figure 8, one can also see an exponential increase of ν after the glitch recovery, whichis a very common, but not well understood behavior (Lyne 1992, see an example of a Crabglitch). We suggest that the exponential increase component is probably a slow glitch, andthe fact that a slow glitch following a classical glitch recovery may be an important clue tothe enigmas of glitch phenomena.
4. Testing Several Fitting Procedures Based on the PhenomenologicalSpin-down Model
In section 3, we showed that the recovery processes of glitches and slow glitches canbe well modeled by the phenomenological model, which can also be used to simulate a realglitch recoveries, in order to fully test different fitting procedures. We get Φ( t ) by integratingEquation (13) for a certain τ i and κ i , and get the TOA set { Φ( t i ) } by assuming a certain∆ T int . We test the biases produced by the following four fitting procedures in the section; three of them are discussed in section 2 for the simplified model of classical andslow glitches. Here all four procedures are examined with a more “realistic”model, i.e., our phenomenological spin-down model. Fitting Procedure I: obtain { ˙ ν ( t i ) } by fitting { Φ( t i ) } to Equation (1), and get τ and∆ ν d by fitting { ˙ ν ( t i ) } to Equation (3) (e.g. Roy et al. 2012, Espinoza et al. 2011, Yuan etal. 2010, Zou et al. 2008, Zou et al. 2004, Dall’Osso et al. 2003, Urama 2002, Dodson etal. 2002, McCulloch et al. 1990). We take τ = 50 days, κ = 0 .
03 in the model, and showinstantaneous results and the fitted results for different ∆ T int in Table 4. The instantaneous∆ ν d is given by ∆ ν d = ν κτ / τ c . It is found that both τ and ∆ ν d are seriously biased. It isalso noticed that the total time span T s (i.e. the time span from the begin of the glitch to theend of the recovery assumed), has considerable impact on the fitting results else. We takehigher order polynomials to fit the phase and find that τ is also seriously biased, even worsethan the lower order one. Thus, if one takes a higher order polynomial and calls the first(linear) term “frequency” then this is clearly not a good approximation, since a higher orderpolynomial would lead to this not being the “frequency”, given that part of it is reabsorbedinto other coefficients. Fitting Procedure II: obtain { ˙ ν ( t i ) } by fitting { Φ( t i ) } to Equation (8), and get τ and 24 –∆ ν d by fitting { ˙ ν ( t i ) } to Equation (3) (e.g. Shabanova 2005, Shabanova 1998). Firstly, weconduct the one decay component case and take τ = 50 days, κ = 0 .
03 in the model. Themain results are shown in Table 5. One can see that the instantaneous values can be wellrestored for ∆ T int = 10 s, and the results with ∆ T int = 10 s are also good approximations.Then, we conduct the two-component case and take τ = 21 . τ = 147 days, and κ = 0 . κ = 0 .
012 in the model. One can see that the instantaneous values can onlybe restored for ∆ T int = 10 s. It is noticed that T s should be long enough for both theone-component and two-component cases. Thus, this procedure is a good approximation onvery small ∆ T int . Fitting Procedure III: get τ and ∆ ν d directly by fitting { Φ( t i ) } to Equation (10) (Yu etal. 2013, Edwards et al. 2006, Shemar & Lyne 1996). Firstly, we also conduct the one decaycomponent case and take τ = 50 days, κ = 0 .
03 in the model. The main results are shownin upper part of Table 6. It is found that the instantaneous values can be well restored andthe fit results is nearly independent of ∆ T int . Then, we conduct the two-component caseand take τ = 21 . τ = 147 days, and κ = 0 . κ = 0 .
012 in the model. We fit thetwo decay terms simultaneously (at a high computing cost) and the main results are shownin the bottom part of Table 6. It is also found that the instantaneous values can be restoredsatisfactorily and the results are independent of ∆ T int . However, T s should not be too longfor both the one component and two components cases, which is opposite to procedure II. Fitting Procedure IV: the phase { Φ( t i ) } is fitted by a very high order polynomial, suchas Φ( t ) = Φ + ν ( t − t i ) + 12 ˙ ν ( t − t i ) + 16 ¨ ν ( t − t i ) + · · · + 150! ν (50) ( t − t i ) . (16)The fitted polynomial (Φ( t ), a continuous function) can very precisely describe the TOAseries { Φ( t i ) } . One can then take its first or second derivative to obtain ν or ˙ ν , i.e. ν = Φ ′ ( t )or ˙ ν = Φ ′′ ( t ). This procedure is suggested by the anonymous referee.We also simulate the one component and two component cases, respectively. The resultsTable 4. Classical glitch simulations with the phenomenological model. ∆ ν d (0.1 µ Hz)and τ (days) obtained by fitting procedure I. T s = 1 year T s = 3 years T s = 5 years τ ∆ ν d τ ∆ ν d τ ∆ ν d Instantaneous 50 1.01 50 1.01 50 1.01∆ T int = 10 s 213.38 1.15 96.45 2.35 92.34 2.22∆ T int = 10 s 161.99 6.49 90.74 2.51 87.52 2.11∆ T int = 10 s 27.02 3.68 37.75 3.07 38.75 3.04
25 –Table 5. Classical glitch simulations with the phenomenological model. ∆ ν d (0.1 µ Hz)and τ (days) obtained by fitting procedure II. T s = 1 year T s = 3 years T s = 5 years τ ∆ ν d τ ∆ ν d τ ∆ ν d Instantaneous 50 1.01 50 1.01 50 1.01∆ T int = 10 s 49.95 1.00 49.97 1.00 49.98 1.00∆ T int = 10 s 45.84 0.95 47.82 0.96 48.01 0.96∆ T int = 10 s 22.23 5.02 27.37 3.15 27.85 2.42Instantaneous (21.40, 147.00) (1.90, 1.19) (21.40, 147.00) (1.90, 1.19) (21.40, 147.00) (1.90, 1.19)∆ T int = 10 s (21.14, 113.58) (1.85, 0.86) (21.26, 144.47) (1.87, 1.18) (21.27, 145.45) (2.33, 1.20)∆ T int = 10 s (2.26, 23.86) (0.59, 1.70) (9.58, 42.54) (0.76, 1.57) (14.40, 81.41) (1.90, 1.45) Table 6. Classical glitch simulations with the phenomenological model. ∆ ν d (0.1 µ Hz)and τ (days) obtained by fitting procedure III. T s = 1 year T s = 3 years T s = 5 years τ ∆ ν d τ ∆ ν d τ ∆ ν d Instantaneous 50.00 1.01 50 1.01 50 1.01∆ T int = 10 s 50.16 1.02 53.13 0.99 67.23 0.83∆ T int = 10 s 50.16 1.02 53.11 0.99 67.01 0.84∆ T int = 10 s 50.15 1.02 52.86 1.00 65.29 0.87Instantaneous (21.40, 147.00) (1.90, 1.19) (21.40, 147.00) (1.90, 1.19) (21.40, 147.00) (1.90, 1.19)∆ T int = 10 s (21.38, 147.92) (1.92, 1.20) (21.73, 152.13) (1.92, 1.20) (18.98, 160.90) (2.33, 1.20)∆ T int = 10 s (21.40, 147.93) (1.92, 1.20) (21.76, 152.17) (1.92, 1.20) (19.04, 160.87) (2.31, 1.20)∆ T int = 10 s (21.40, 147.96) (1.92, 1.20) (21.77, 152.19) (1.93, 1.20) (19.35, 160.57) (2.20, 1.20)
26 – ( H z s - ) Time (days)
Instantaneous-1 Restored-1, 50 Instantaneous-2 Restored-2, 50 Restored-2, 30 Restored-2, 20 Restored-2, 15 . Fig. 9.— Classical glitch simulations with the phenomenological model. ˙ ν ( t ) for both one andtwo component cases (denoted as 1 and 2 in the legend, respectively) are obtained by fittingprocedure IV. Note that restored results are also continuous functions; here we take theirdiscrete values in order to compare them more conveniently with the instantaneous values.The numbers 50, 30, 20 and 15 in the legend denote the orders of the fitting polynomials. 27 –are shown in Figure 9. One can see that the instantaneous values of ν ( t ) or ˙ ν ( t ) can berestored with very high precision for both the cases. In the figure, ∆ T int = 10 s is taken,and it is checked that the results are almost independent of both ∆ T int and T s . Then,one can get τ i and ∆ ν d i by fitting the restored ˙ ν ( t ) to Equation (3). The fitted glitchparameters for one component case are τ = 50 .
00 days and ∆ ν d = 1 . × − Hz, and fortwo component case are τ = 21 .
37 days, τ = 146 .
93 days, and ∆ ν d1 = 1 . × − Hz,∆ ν d2 = 1 . × − Hz. They are all consistent with instantaneous values (see e.g. Table6) with very high precisions. In Figure 9, we show the fitting results of two component casewith different order polynomials else. It is found that the order of the polynomial must bevery high, e.g. &
35, which requires that the TOA data points should not be too sparse. Wealso test the fit procedure with different values of τ i and ∆ ν d i , and all the glitch parametersare restored satisfactorily.In conclusion, procedure III is a reasonable choice to get τ and ∆ ν d ; however, the twocomponents should be fit simultaneously (in order to avoid some local minimum of χ ),and T s should not be too long. Procedure IV seems to be the best choice for pulsar glitchdata analysis, which gives { ν ( t ) } and { ˙ ν ( t ) } with very high precision, and then the glitchparameters τ i and ∆ ν d i can be satisfactorily estimated by fitting the restored ˙ ν ( t ) to Equation(3). We thus suggest that theorists should always use the full timing solution, rather thantry to compare models to individual parameters of fits, as these may be highly inaccurate.Furthermore working in phase seems to be the most accurate and reliable method.
5. Discussions5.1. How to obtain the correct model parameters of pulsars?
We have shown recently that fitting the observed TOAs of a pulsar to Equation (1) willresult in biased (i.e., averaged) spin-down parameters, if its spin-down is non-secular andthe variation time scale is comparable to or shorter than the time span of the fitting (Zhang& Xie 2012a, 2012b). In particular we predicted that the reported braking index should bea function of time span and approaches to a small and positive value when the time span ismuch longer than the oscillation period of its spin-down process, which can be tested withthe existing data (Zhang & Xie 2012b).We notice that in some of the literature (e.g. Roy et al. 2012, Espinoza et al. 2011,Yuan et al. 2010) only Equation (1) is referred to, when describing the fitting process, evenfor the glitch data analysis. However we have shown in Figure 1 that this will produce asignificantly distorted glitch profile. Instead, one can fit to Equations (7) and (8) to obtain 28 –the un-distorted (but still averaged) glitch profile; probably this is usually done in practice,though not explicitly described in the literature (Yu 2003). We suggest that the exact fittingprocedures should be described when reporting the analysis results of the observed glitchdata.However, neither Equation (1) nor Equations (7) and (8) describe exactly the physicalspin-down processes of all pulsars. The pulsar parameters fitted by Equations (10) are alsoslightly biased especially if T s is not properly taken. Ideally, a spin-down model (empirical,phenomenological or physical) should be used to the fit the observed TOA data, in orderto obtain the model parameters. To serve this purpose, the observed TOAs of each pulsarshould be made available, and the exact fitting procedure should be described along withthe reported spin-down parameters of a pulsar. As shown in Figure 10, we simulate a glitchrecovery with parameters τ = 21 . τ = 147 days, and ∆ ν d1 = 1 . × − Hz,∆ ν d2 = 1 . × − Hz, and ∆ T int = 10 s. Then we have TOAs from the phenomenologicalmodel, and the simulated “reported” { ˙ ν } is obtained by fitting TOAs to Equation (8) andis represented by the solid line. By fitting TOAs to Equation (10), we have the “reported”glitch parameters τ = 21 . τ = 152 days. With the timescales, we simulate { ˙ ν } againwith ∆ T int = 10 . The model parameters ∆ ν d1 and ∆ ν d2 can be adjusted until simulatedfits match the “reported” ones, and the best fit model parameters ∆ ν d1 = 1 . × − Hz,∆ ν d2 = 1 . × − Hz, which are agree well with original parameters. We show the restored { ˙ ν } with circles. With the same parameters, the restored { ˙ ν } for ∆ T int = 5 × and 2 × s are represented by diamonds and triangles, respectively. One can see that { ˙ ν } can be wellrestored if TOAs are known. If ∆ T int taken in simulation is not the right one, the { ˙ ν } profilesare apparently different from the “reported” one, even though the model parameters are allcorrect.When TOAs are available, one can then follow the steps we used above to combine amodel with simulations to obtain model parameters. Alternatively, Φ( t ), ν ( t ) and ˙ ν ( t ) givenby procedure IV can also be fitted directly by physical models. In the above analysis, we have assumed that t is known; however, t is usually takenas the averaged time of the last reported TOA just before the glitch and the first reportedTOA of the glitch. This means we have an uncertainty in t : σ t = ∆ T int /
2. Then fromEquation (6) for a classical glitch, we find 29 –
10 100 1000-880-860-840-820-800-780 ( - H z s - ) Time (days) "Reported", 10 s Restored, 10 s Restored, 5 10 s Restored, 2 10 s . Fig. 10.— The effects of ∆ T int on variation of ˙ ν . The “reported” { ˙ ν } is represented by thesolid line, the restored { ˙ ν } with ∆ T int = 10 , 5 × and 2 × s are represented by circles,diamonds and triangles, respectively. 30 – σ ∆ ν ∆ ν = σ ∆ ˙ ν ∆ ˙ ν = σ t τ , (17)where σ ∆ ν and σ ∆ ˙ ν are the uncertainties of the restored ∆ ν and ∆ ˙ ν , respectively. For theclassical glitch of B2334+61, σ t ∼ . τ ∼ . σ ∆ ν / ∆ ν = σ ∆ ˙ ν / ∆ ˙ ν ∼ σ ∆ ν ≈ ν d is determined by the dataat the end of the recovery, i.e. ν ∼ ∆ ν d for ∆ t ≫ τ from Equation (9). However, fromthe derivative of Equation (9) , we have ˙ ν = ∆ ν d τ e − ( t − t ) /τ and ∆ ˙ ν d ( ≡ ∆ ν d /τ ) is closelyrelated to t , which resembles the case of a classical glitch. However, for the slow glitch wecan fit for t of the glitch by calculating where the rise and the pre-glitch solutions intersect,which will cause a much smaller uncertainty. This is a major difference from analyzing thedata of a classical glitch. Unfortunately, this has not been realized previously and thus t wasnot determined from the reported ν O with this method for slow glitch data analysis. Thiscauses an uncertainty to ∆ ˙ ν in the same way as in Equation (17), i.e., the bias is related to∆ T int . For instance, in Figure 2 of Zou et al. (2004), the observed results for a slow glitchevent of B1822–09 are ∆ ν a = (40 . ±
26) nHz and ∆ ˙ ν a ≃ . × − s − ; and for the sameevent, the results in Shabanova (2005) are ∆ ν b = 40 . ν b ≃ . × − s − . Asexpected above, ∆ ν a = ∆ ν b , but ∆ ˙ ν a = ∆ ˙ ν b . For the event, τ ∼
110 d, and σ at ∼ . σ bt ∼ . σ a ∆ ˙ ν = 1 . × − s − , σ b ∆ ˙ ν = 6 . × − s − ,and σ ∆ ˙ ν = p ( σ a ∆ ˙ ν ) + ( σ b ∆ ˙ ν ) = 6 . × − s − . Then we have (∆ ˙ ν a − ∆ ˙ ν b ) /σ ∆ ˙ ν ≃ . σ ∆ ˙ ν . Based on observational results, we generalize the variations of ν and ˙ ν for slow andclassical glitch recoveries, as shown in Figure 11. The pre-glitch tracks are representedby dotted line. After the jump, the classical glitch recoveries (represented by solid line)generally have ν variation that tends to restore its initial values, and usually the restorationis composed by a exponential decay and a permanent linear decrease with slope ∆ ˙ ν p ; however,for slow glitches (represented by dashed line), ν monotonically increases, as shown in panel(1). In panel (2), ˙ ν of classical glitch recoveries that tends to restore its initial values, butcannot completely recover for ∆ ˙ ν p = 0; ˙ ν of slow glitch recoveries almost completely recoverto its initial value, corresponding to the increase of ν .In sections 3 and 4, we have shown that the classical and slow glitch recoveries can bewell modeled by a simple function, G ( t ) = 1 + κ exp ( − ∆ t/τ ), with positive or negative κ , asshown in panel (3), respectively. However, it is should be noticed that the model only have 31 – . d p d p p =0 . (3)(2) .. p . dd . -2 0 2 4 6 8 . Classical Glitch Slow Glitch Before Glitch G ( t ) Time ( t ) (1)
Fig. 11.— Schematic depictions of ν , ˙ ν and G ( t ) for the slow and classical glitch recoveries.The pre-glitch tracks are represented by dotted line. The classical glitch recoveries arerepresented by solid lines. The slow glitches are represented by dashed lines. 32 –two parameters, κ and τ , from which we can obtain ∆ ν d and ∆ ˙ ν d , but not ∆ ν p and ∆ ˙ ν p ,which are not modelled. Nevertheless, we conclude that the major difference between slowglitch and classical glitch recoveries are that they show opposite trends with opposite signsof κ , in our phenomenological model.
6. Summary
In this work we studied the data analysis procedures of pulsar’s glitch observationsand found the conventionally used methods produce biases to the true glitch parameterswith varying degrees. We presented a phenomenological model for the recovery processes ofclassical and slow glitches, which is used to model successfully the observed slow and classicalglitch events from pulsars B1822–09 and PSR B2334+61, respectively. Based on the model,we tested four different data analysis procedures. Our main results are summarized asfollows:1. The timing analysis method of fitting the observed TOAs with Equation (7) or Equa-tion (8) results in significant biases to glitch parameters of variation magnitude, asshown in Figures 2 and 3, and Table 1. The biases can be ignored only when ∆ T int ≤ s; otherwise biases still exist to some extend.2. With Equation (10), one can obtain the glitch parameters by fitting the phase di-rectly, which produce relatively smaller biases. However, for the case with multipledecay terms, the timescales are usually fixed by eye for their initial values, which mayintroduce strong biases.3. We propose a phenomenological model of glitch recovery (Equation (13)), which canreproduce the commonly observed exponential glitch recovery profiles. The recoveryprocesses of both slow and classical glitches can be explained as the G ( t ) = 1 + κ exp ( − ∆ t/τ ) with κ < κ > ν ( t ) and ˙ ν ( t ). Then the glitch parameters can be obtainedfrom ν ( t ) and ˙ ν ( t ) (e.g. fitting ˙ ν ( t ) to Equation (3)). We suggest that this procedureshould be used in pulsar timing analysis.5. The uncertainty in the starting time ( t ) of a classical glitch causes uncertainties to 33 –the glitch parameters ∆ ν d and ∆ ˙ ν d (Equation 17), but less so to a slow glitch and t of a slow glitch can be determined from data.However our phenomenological model cannot account for the non-recoverable jumps in ν and ˙ ν , which are observed for some classical glitches and may be due to the permanentincrease of a pulsar’s dipole magnetic field due to glitches (Lin & Zhang 2004). In the work,we also assumed uniform TOA distributions to simulate both the slow and classical glitchrecoveries, since the observed TOAs are not reported in literature. The glitch parameterscan be better restored, if the observed TOAs are available and fitted directly with a glitchmodel; this is actually generally desired for pulsar timing studies. Thus we suggest thatTOAs should be made available to the community when possible or that the full fittingprocedure and fit parameters for different epochs made available. Also theorists could tryto calculate phase as an output, thus making the comparison more accurate.We thank Jianping Yuan and Meng Yu for valuable discussions. We would like to thankthe anonymous referee for his/her comments and suggestions that led to a significant improve-ment of this paper. SNZ acknowledges partial funding support by 973 Program of Chinaunder grant 2009CB824800, by the National Natural Science Foundation of China undergrant Nos. 11133002 and 10725313, and by the Qianren start-up grant 292012312D1117210. REFERENCES
Anderson, P. W., & Itoh, N. 1975, Nature, 256, 25Anderson, P. W., Alpar, M. A., Pines, D., & Shaham, J., 1982, Philos. Magazine A, 45, 227Andersson, N., Comer, G. L., & Prix, R. 2003, Phys. Rev. Lett., 90, 091101Andersson, N., Glampedakis, K., Ho, W. C. G., & Espinoza, C. M. 2012, Phys. Rev. Lett.,109, 241103Alpar, M. A. 1977, ApJ, 213, 527Alpar, M. A., Anderson, P. W., Pines, D., & Shaham, J., 1981, ApJ, 249, L29Alpar, M. A., Anderson, P. W., Pines, D., & Shaham, J. 1984b, ApJ, 278, 791Alpar, M. A., Chau, H. F., Cheng, K. S., & Pines, D. 1993, ApJ, 409, 345Alpar, M. A., Chau, H. F., Cheng, K. S., & Pines, D. 1996, ApJ, 459, 706 34 –Alpar, M. A., Nandkumar, R., & Pines, D. 1985, ApJ, 288, 191Alpar, M. A., Pines, D., Anderson, P. W., & Shaham, J. 1984a, ApJ, 276, 325Baym, G., Pethick, C., & Pines, D., 1969, Nature, 224, 872Baym, G., & Pines, D., 1971, Ann. Phys., 66, 816Biryukov, A., Beskin, G., & Karpov, S. 2012, MNRAS, 420, 103Chamel, N., 2013, Phys. Rev. Lett., 110, 011101Cheng, K. S., Chau, W. Y., Zhang, J. L., & Chau, H. F., 1992, ApJ, 396, 135Chau, H. F., McCulloch, P. M., Nandkumar, R., & Pines, D., 1993, ApJ, 413, L113Dall’Osso, S., Israel, G. L., Stella, L., Possenti, A., Perozzi, E., 2003, ApJ, 599, 485Dewey, R. J., Taylor, J. H., Weisberg, J. M., & Stokes, G. H. 1985, ApJ, 294, L25Dodson, R. G., McCulloch, P. M., & Lewis, D. R., 2002, ApJ, 564, L85Dodson, R. G., Lewis, D. R., & McCulloch, P. M. 2007, Ap&SS, 308, 585Edwards, R. T., Hobbs, G. B., & Manchester, R. N., 2006, MNRAS, 372, 1549Espinoza, C. M., Lyne, A. G., Stappers, B. W., & Kramer, M. 2011, MNRAS, 414, 1679Glampedakis, K., & Andersson, N., 2009, Phys. Rev. Lett., 102, 141101Haskell, B., & Antonopoulou, D., 2013, arXiv:1306.5214Haskell, B., Pizzochero, P. M., & Seveso, S. 2013, ApJ, 764, L25Haskell, B., Pizzochero, P. M., & Sidery, T. 2012, MNRAS, 420, 658Hobbs, G., Lyne, A. G., & Kramer, M. 2010, MNRAS, 402, 1027Jones, P. B., 1990, MNRAS, 243, 257Jones, P. B., 1992, MNRAS, 257, 501Jones, P. B., 1998, Phys. Rev. Lett., 81, 4560Larson, M. B., & Link, B., 2002, MNRAS, 333, 613Lin, J. R. & Zhang, S. N. 2004, ApJ, 615, L133 35 –Link, B., Epstein, R. I., & Baym G. 1992, ApJ, 390, L21Link, B., & Epstein, R. I., 1996, ApJ, 457, 844Link, B., Epstein, R. I., & Baym G. 1993, ApJ, 403, 2Link, B., Epstein, R. I., & van Riper, K. A. 1992, Nature, 359, 616Link, B., Franco L. M. & Epstein, R. I. 1998, ApJ, 508, 838Lyne, A. G., Smith, F. G., Pritchard, R. S. 1992, Nature, 359, 706Lyne, A. G., Hobbs, G., Kramer, M., Stairs, I., & Stappers, B. 2010, Science, 329, 408Lyne, A. G., Shemar, S. L, & Smith, F. G, 2000, MNRAS, 315, 534McCulloch, P. M., Hamilton, P. A., McConnell, D., & King, E. A., 1990, Nature, 346, 822Peng, C., & Xu, R. X. 2008, MNRAS, 384, 1034Pines, D., Shaham, J., Alpar, M. A., & Anderson, P. W., 1980, Progress Theor. Phys. Suppl.,69, 376Pons, J. A., Vigan`o, D., & Geppert, U. 2012, A&A, 547, A9Pizzochero, P. M., 2011, ApJ, 743, L20Roy, J., Gupta, Y., Lewandowski, W., 2012, MNRAS, 424, 2213Ruderman, M. 1972, ARA&A, 10, 427Ruderman, M. 1976, ApJ, 203, 213Seveso, S., Pizzochero, P. M., & Haskell, B., 2012, MNRAS, 427, 1089Shemar, S. L., & Lyne, A. G., 1996, MNRAS, 282, 677Sedrakian A. D., & Cordes J. M., 1999, MNRAS, 307, 365Shabanova, T. V. 1998, A&A, 337, 723Shabanova, T. V. 2005, MNRAS, 356, 1435Shemar, S. L., & Lyne, A. G. 1996, MNRAS, 282, 677Standish, E. M. 1998, A&A, 336, 381 36 –Urama, J. O., 2002, MNRAS, 330, 58van Eysden, C. A., & Melatos, A. 2010, MNRAS, 409, 1253Wang, N., Manchester, R. N., Pace, R. T., Bailes, M., Kaspi, V. M., Stappers, B. W. &Lyne, A. G. 2000, MNRAS, 317, 843Yuan, J. P., Manchester, R. N., Wang, N., Zhou, X., Liu, Z. Y., & Gao, Z. F. 2010, ApJ,719, L111Yu, M., et al. 2013, MNRAS, 229, 688Yu, M. 2013, personal communicationZou, W. Z., Wang, N., Wang, H. X., Manchester, R. N., Wu, X. J., & Zhang, J. 2004,MNRAS, 354, 811Zou, W. Z., Wang, N., Manchester, R. N., Urama, J. O., Hobbs, G., Liu, Z. Y., & Yuan, J.P., 2008, MNRAS, 384, 1063Zhang, S.-N., & Xie, Y. 2012, ApJ, 757, 153Zhang, S.-N., & Xie, Y. 2012, ApJ, 761, 102