Monotonic and cyclic components of radio pulsars spindown
Abstract
In this article we revise the problem of anomalous values of pulsars’ braking indices and frequency second derivatives arising in observations. The intrinsic evolutionary braking is buried deep under superimposed irregular processes, that prevent direct estimations of its parameters for the majority of pulsars. We reanalyze the distribution of “ordinary” radio pulsars on a , , and diagrams assuming their spindown to be the superposition of a “true” monotonous term and a symmetric oscillatory term. We demonstrate that their effects may be clearly separated using simple ad hoc arguments. Using maximum likelihood estimator we derive the parameters of both components. We find characteristic timescales of such oscillations to be of the order of years, while its amplitudes are large enough to modulate the observed spindown rate up to 0.55 times and completely dominate the second frequency derivatives. On the other hand, pulsars’ secular evolution is consistent with classical magnetodipolar model with braking index
So, observed pulsars’ characteristic ages (and similar estimators that depend on the observed ) are also affected by long term cyclic process and differ up to 0.55 times from their monotonous values. This fact naturally resolves the discrepancy of characteristic and independently estimated physical ages of several objects, as well as explains very large, up to years, characteristic ages of some pulsars.
We discuss the possible physical connection of long term oscillation with a complex neutron star rotation relative its magnetic axis due to influence of the nearfield part of magnetodipolar torque.
keywords:
stars: pulsars: general – methods: statistical PACS 04.40.Dg – 95.75.Pq – 97.60.Gb – 98.62.Ve1 Introduction
Radio pulsars are highly periodic variable astrophysical objects, powered by neutron stars’ rotation with periods evolving with time. Observational determination of their timing parameters – instantaneous period and its derivatives – is an extremely complex task (Edwards et al., 2006; Cordes & Shannon, 2010), but after the correction of pulse times of arrivals (TOAs) for the effects of radio wave propagation in the interstellar medium, motion and position of the Earth, gravitational delays in Solar system potential etc., the resulting phase of the light curve may be well described by an (infinite) Taylor series
(1) 
dominated by the lower order terms. This may also be expressed in terms of observed rotational frequency as
(2) 
where values of and its derivatives are attributed to the epoch .
At the same time, the majority of published theoretical models of pulsar braking predict the spindown described by a differential equation of the form (Manchester & Taylor, 1977; Beskin et al., 1993)
(3) 
where is a (constant) braking index and depends on the individual physical properties of the pulsar. For the vacuum magnetodipolar model the canonical value is , pulsar wind decreases it to and multipole magnetic field increases it to (Manchester & Taylor, 1977).
For such a spindown with constant the braking index is equal to a simple combination of frequency and two its derivatives ^{1}^{1}1Note that the observed ones, estimated by fitting the TOAs with a model of (1) broken down to the third order term, are always biased even if the spindown is exactly according to equation (3) due to influence of nonzero higherorder terms. For actual pulsars, however, this bias is negligible and will be ignored below.
(4) 
That is why the determination of spin frequency second derivative via pulsar timing is important for understanding the pulsar’ spindown.
Several decades of very detailed timing of hundreds of pulsars (D’Alessandro et al., 1993; Baykal et al., 1999; Chukwude, 2003; Hobbs et al., 2004, 2010; Livingstone et al., 2005), including the best studied one – Crab pulsar (Scott et al., 2003) – demonstrated, however, that the measured rotational phase (1) evolution is not generally consistent with the one expected for a braking law (3) with physically reasonable parameters. Observed rotational phase (as well as , ) includes components with a very complex, irregular behaviour. The first kind of such irregularities is “glitches” – sporadic fast changes of pulsars’ periods and spindown rate with up to few months relaxation timescale (e.g. Manchester & Taylor (1977)), while the other one – quasirandom phase variations with typically red Fourier spectrum and characteristic timescales of a few months to few years – “timing noise” (e.g Cordes & Helfand (1980); Lyne (1999); Hobbs et al. (2010)).
Both of them affect the measurements of , and , and often make them dependent on the duration and epoch of the observations’ time span, when it is shorter or comparable to timescales of these processes. But, if observations are performed over sufficiently long intervals of about a decade or longer, the observed coefficients of series (1) become more stable – the influence of these irregularities is mostly suppressed (see Sec. 2.1 for additional discussion). To date, there are more than 200 pulsars observed for longer than 15 years (Baykal et al., 1999; Chukwude, 2003; Hobbs et al., 2004, 2010), and still the majority of values, measured over such long intervals, lead to extremely high, up to (and even more), braking indices. Moreover, their values are negative for nearly a half of all objects. In general, measured braking indices drastically differ from physically reasonable values, and this raises serious problem on the pulsars’ long term spin evolution.
While glitches are widely believed to be caused by a rapid transfer of momentum from the neutron stars interior (e.g. Espinoza et al. (2011) and references therein), the nature of the timing noise is still unclear and no selfconsistent and widely accepted model of pulsars’ phase residuals irregularities is constructed. There have been numerous attempts to attribute it to various stochastic phenomena in the neutron star magnetosphere (Cheng et al., 1987; Contopoulos, 2007; Lyne et al., 2010), in the interior (Cordes & Downs, 1985), to spindown torque variations (Urama et al., 2006), and to the existence of a number of different spindown regimes for a single pulsar (Lyne et al., 2010). Purely phenomenological attempts to describe observed noise as random walks of different orders (in phase, frequency or its derivative) (Cordes & Greenstein, 1981) have also not succeeded as with the increase of observations time spans the noise appeared to be more complex than these simple models predict. Some characteristics of this essentially irregular process, however, have been parametrized through several noise strength parameters and their correlations. So, Arzoumanian et al. (1994) parametrized pulsars’ timing noise through quantity and found a good correlation, where is a pulsar period.
In turn, Cordes & Helfand (1980), Cordes & Downs (1985) and then Chukwude (2003) have investigated activity parameter , which is a direct measure of timing residuals variance. A good correlations of versus and was also found.
However, the quantities used in these papers to parametrize timing noise were in fact just the measures of observed frequency second derivative. Even the activity parameter was calculated using timing residuals after subtraction of a second order polynomial in series (1). Therefore, the correlations obtained simply reflect the quite strong dependence.
At the same time, only a few longterm mechanisms have been offered to explain observed anomalous braking indices. Demiańsky & Proszyński (1979) studied the variations of pulsar timing parameters due to existence of possible massive partner – most of radio pulsars are, however, believed to be isolated objects. Alpar & Baykal (2006) investigated the impact of unresolved or missed glitches on the observed timing parameters, but it is very difficult to explain some extremely high braking indices by such mechanism. Gullahorn & Rankin (1977) considered variations of pulsar braking torque with timescales of years, probably due to interaction of the pulsar with interstellar medium in its vicinity. The idea of spindown torque variations itself is quite promising and reasonable, but the nature of variations on such timescale was not pointed out in their note and remained unclear.
In turn, in Beskin et al. (2006) and Biryukov et al. (2007) we considered the existence of a longterm cyclic variational process affecting pulsars’ spindown on a timescale of thousands of years. One may call such a process a third type of timing irregularities. Later, Barsukov & Tsygan (2010) suggested that the physical nature of such a process may be related to the “anomalous” magnetodipolar braking torque, which may produce a forced precession of pulsar’s rotational axis around its magnetic moment with the period of years, and demonstrated the possibility of significant variations of observed rotational parameters on such timescale. Even earlier, Melatos (2000) also investigated the triaxial neutron star rotation taking into account this torque and demonstrated that the rotation of such star may be very complex and even induce variations of magnetic inclination angle, which will also affect the neutron star spindown.
As a further development of this concept of a longterm variational process, in the current work we analyse the ensemble of 297 pulsars with published data on second frequency derivatives and demonstrate that these values may be used to estimate the parameters of both secular (evolutionary, monotonous) and additional, cyclic, components of a pulsar spindown under simple and reasonable assumptions. We formulate such a twocomponent model of pulsar braking and derive its parameters using a maximum likelihood estimator. The results are quite reasonable, as the parameters of monotonous component are in a good agreement with a standard magnetodipolar spindown model.
The article is organized as follows. In Section 2 we justify that measured values can indeed be used to analyse pulsars’ spindown. In Section 3, statistical analysis of a 297 objects is performed, phenomenology of a twocomponent pulsars’ spindown is presented and parameters of its irregular term are roughly estimated using modelindependent arguments. In Section 4 this twocomponent model is used for the determination of the secular spindown parameters also. In Section 5 we discuss the results and its astrophysical implications. The conclusions are given in Section 6.
2 Pulsars’ measurements
2.1 What ’s can tell us on the pulsars’ physics?
Anomalous values of (and braking indices) for most of ordinary pulsars raise at least two principal questions, critical for any serious analysis of these quantities. Namely:

Are values measured over long timespans (decades) really related to pulsars’ or interstellar medium^{2}^{2}2The effects of propagation of pulsar radiation through its “near” and “far” medium can introduce an additional delays to the times of pulses arrivals and therefore affect the timing solution (Cordes & Lazio, 2002; Cordes & Shannon, 2010). physics or they are just artifacts of somewhat incorrect observations and/or data reduction?
and

If they are not artifacts, then is the anomaly induced by the same phenomena responsible for the short time scale (months) quasiperiodic timing irregularities with redlike spectra, observed in timing residuals?
The answer to the first question above is clear. Pulsar timing is a complex, but welldefined and thoroughly described procedure, and the sources of observational uncertainties arising in it are well known and are taken into account during the reduction, as well as various effects that affect the timing systematically – e.g. motion of the Earth, relativistic effects, etc (see, for example, Edwards et al. (2006)). There is no reason to believe the significant cubic trends, clearly seen in phase solutions for hundreds of pulsars (Hobbs et al., 2004), are artifacts.
There is also a method for indirect measurements of a pulsar’ second derivative by comparing values acquired over sufficiently distant, separated by years, time spans (Johnston & Galloway, 1999). Results of such estimation are anomalous too, and are generally consistent with directly measured values.
Obviously, all accurately measured values, even being anomalous, are physically meaningful – they are indeed due to peculiarities of pulsars’ spin down, while their anomalousness is just an indication of insufficiency of our models of pulsars’ braking.
On the other hand, the answer to the second question is not so obvious, as stochastic timing irregularities, seen in a large fraction of pulsars, still are not well understood and are complex phenomena.
Phenomenologically, timing residuals’ irregularities and unexpected anomalous parameters of timing solution have to be separated in the analysis. Indeed, the former ones are seen directly within observational datasets and their stochastic nature is clear while the origin of the latter ones can only be speculated about. Below we argue in favour of indeed different origins of these phenomena.
The phase residuals after all appropriate corrections and quadratic trend subtraction (i.e. after extraction of and only) form a very plural zoo over the pulsar population. There are at least several different classes of typical pulsars’ timing series (see Figure 3 in Hobbs et al. (2010)) – the ones with dominant influence of a cubic term (e.g. PSRs B095954, B0943+10, B165713 etc.), the ones with more complex but generally smooth behaviour (PSRs B162026, B170616, B172747 etc.), and the ones that are purely white noise supposedly due to measurement errors. The latter class typically contains pulsars with inaccurate or unmeasured second frequency derivatives, and will not be analysed here. The second one is often used to depict the timing deviations from an expected spindown law as a process with red noiselike spectrum; there are, however, some evidences that this behaviour in fact is a combination of a cubic order term and a quasiperiodic component (Hobbs et al., 2010). Moreover, young pulsars (e.g. Crab, B150958, B2011+38) show strong shortterm quasiperiodic timing noise, along with cubic trends corresponding to a spindown with physically meaningful braking indices (Lyne et al., 1993; Scott et al., 2003; Livingstone et al., 2005; Hobbs et al., 2010).
Also, the behaviour of ’s with the increase of observational span is not consistent with their origin due to stochastic or short timescale noise processes. For example, for PSR B170616, variations of with an amplitude of s have been detected on a timescale of several years (see Figure 7 in Hobbs et al. (2004)) – its values, however, are always around the one revealed by the fit over the entire 25 year time span ( s), which still leads to a braking index .
All this strongly suggest that timing noise is a process distinct from the one producing large (and anomalous) cubic trends in (1); this noise acts mostly as a randomizing factor decreasing the accuracy of measurements, but not defining their properties in general.
As a result, we believe the second derivative values to be very promising tool for studying long timescale features of pulsars’ spindown. Although the values of individual pulsars bring little information on the spindown of neutron stars, the properties of its distribution over an ensemble of all pulsars seems to be appropriate for spindown study under a simple and reasonable assumptions on the properties of process producing these anomalous cubic trends.
2.2 The subset
The set of pulsars under investigation is similar to the one used in our previous works (Beskin et al., 2006; Biryukov et al., 2007). From the 393 objects of the ATNF^{3}^{3}3http://www.atnf.csiro.au/research/pulsar/psrcat/, revision from 6th Oct 2007 catalogue (Manchester et al., 2005) with known we have compiled a list of “ordinary” radio pulsars that

have ms and s/s (i.e. had not been recycled);

have relative accuracy of second derivative measurements better than 75%;

are not components of binary systems and not in a list of known anomalous xray pulsars.
19 supplementary pulsars from other sources (D’Alessandro et al., 1993; Chukwude, 2003) have been added. The final set consists of 297 objects including 247 from (Hobbs et al., 2004). 18 of them are associated with young supernova remnants^{4}^{4}4Data on PSRSNR associations was also taken from ATNF.
3 Statistical analysis of pulsars’ timing parameters
3.1 dependency. Cyclic evolution of pulsars’ spindown?
Our work is based on the study of relations between quantities derived from pulsars’ timing: spin frequency, its two derivatives and , observed braking index and characteristic age .
As we previously demonstrated in (Biryukov et al., 2007), the logarithms of show significant correlation with the ones of for both positive (172 pulsars, correlation coefficient ) and negative (125 pulsars, ) branches of the diagram (Figure 1).
The apparent separation of branches on the diagram is primarily due to the logarithmic scale of the plot, and actually objects continuously cover the full range of values, excluding only the gap on small ones due to limited accuracy of measurements (which is not better than ).
Young pulsars confidently associated with supernova remnants are systematically shifted to the left on the diagram (open symbols in Figure 1) and are absent on the right. Hence, pulsars seem to evolve towards lower values of . Moreover, if is an appropriate pulsar age estimator, then is the same too, as they are highly correlated ( in logarithmic scale, see Fig. 5)
On the diagram, shown in Figure 3, pulsars’ behaviour is similar: they are born with higher values of and, since , evolve toward lower values. The direction of the evolution is the same for positive and negative branches of . So all older pulsars have systematically lower values of which is consistent with evolutionary interpretation of diagram.
So, we conclude that each pulsar during its evolution moves along the branches of the and diagrams while its value increases and value decreases. However, on the negative branch of , the value of the first derivative, being negative too, can only decrease with time, since is a formal derivative of , and both of them are regular components of the observed rotational phase (1). So, pulsar motion along the branch may only be backward, which clearly contradicts its evolutionary interpretation suggested earlier. The solution we offer is to assume a cyclic behaviour of pulsars on this diagram. As pulsars evolve, they repeatedly change sign of , in a spirallike motion from branch to branch, and spend roughly half their lifetime on each one. Such behaviour is sketched in Figure 2.
The timescale of such variations has to be much shorter than the pulsar life time, and at the same time significantly larger than the one of the observations to systematically affect the timing solution. Of course, such cyclic behaviour should also affect other spindown parameters – and .
Then, on plot, shown in Figure 4, pulsars move from bottom right to the upper left, from high to low and , forming quasievolutionary sequence, with linear regression coefficient in logarithmic scale. Unfortunately, this diagram itself can’t be used to study the spindown evolution of pulsars, as it is distorted by various selection effects (death line crossing, emission beam width dependency on pulsar frequency (Tauris & Manchester, 1998), etc) and correlation of pulsars’ initial and at birth.
3.2 Phenomenology of observed pulsars’ spindown
Now we introduce an appropriate phenomenological description of the complex timing evolution of pulsars including the cyclic variations we proposed above.
The exact form of a pulsar’s trajectory on the diagram is unknown, but in general its rotational evolution may be described as the superposition of the regular component of the observed rotational phase , and an irregular, cyclic one :
(5) 
where
(6) 
over a long enough time interval comparable to a pulsar’s lifetime. After term by term differentiation of equation (5) one gets:
(7) 
for the rotational frequency and
(8) 
for the spindown rate . Obviously, due to secular loss of rotational energy, the value of should be always negative for an isolated pulsar. Quantity is a spindown rate relative deviation, which is not necessarily small, and is likely greater than 1:
(9) 
due to the absence of “ordinary” pulsars with in the subset under investigation. It is clear that deviation is also cyclic; in the simplest case it is a harmonic function of variational phase :
(10) 
where and is the relative amplitude of the oscillations, related to the absolute amplitude :
(11) 
By second differentiation of equation (5) for a rotational frequency second derivative we get:
(12) 
where (i.e. ) for most pulsars, and this is the reason behind anomalously high values of observed and braking indices. In other words, the absolute amplitude of the variations is much greater than the regular term . Note also that cyclic terms and are not independent – they are connected through the relation:
(13) 
Since pulsars secularly evolve to the higher values of and , the sign of should always be positive:
(14) 
Such positive contribution should introduce a small asymmetry of the observed values of in respect to , even if the variations of are intrinsically symmetric. Such asymmetry drives the average motion to the right on the diagram and affects the times that pulsars spend with positive and negative .
3.3 Asymmetry in observed ’s
Indeed, numbers of pulsars with positive () and negative () values of within the subset are significantly different. If this difference is just accidental, its probability is too small, only, assuming binomial distribution of 297 objects with chance to be on a positive branch. Therefore, null hypothesis of exactly symmetric branches is rejected on a % significance level, and the branches are indeed asymmetric in number of objects^{5}^{5}5At the same time, nonparametric KolmogorovSmirnov test rejects the hypothesis of the common distribution of of positive and negative branches on a 2.5% significance level..
We investigated in detail the behaviour of such significance levels for a number of subsets of pulsars with greater than and less then a given value. Corresponding dependencies are plotted in the Fig. 6. It is clearly seen that branches of the and diagrams are significantly different only when relatively young pulsars are taken into account, while behaviour of the older pulsars seems to be the same for both signs of . It is consistent with the existence of a positive evolutionary trend , which is negligible for older pulsars, but is large enough to affect the second derivatives of the younger ones.
In other words, for the initial stages of pulsars life, values of are less than (i.e. ), and the second derivatives are always positive (see the upper left region of Fig. 1). Later on, as decreases, the becomes quite small and the relative deviations grow up extremely over unity; as a result, observed values of turn out to be significantly different from an intrinsic , and it corresponds to the great increase of their absolute braking indices (see Fig. 7). These different stages of pulsars’ evolution are illustrated in the Fig. 2.
At the same time, the absence of a significant difference between the branches for relatively old pulsars – their symmetry – implies an important fact that variations are indeed approximately symmetric in respect to evolutionary trend: a pulsar deviates to higher and lower values of with nearly the same amplitude and spends an equal amounts of time with greater and less than . This fact will be used in Section 4 for the estimation of the parameters of pulsars’ twocomponent spindown model.
However, some generic characteristics of the cyclic, irregular component of behaviour – its amplitude and timescale – may be easily derived from an exploratory analysis of the and scatter plots using simple, modelindependent arguments. Such analysis is presented in the next two subsections.
3.4 Simple limits on the timescales and amplitudes of the cyclic process
As stated above, for the majority of pulsars the amplitude of the second derivative variations is much greater than the positive evolutionary trend . All the pulsars on the negative branch of diagram are due to subtraction of this amplitude from the trend. Therefore, the values of of the lower envelope of this branch (in terms of ) are sensible estimation of a lower limit on . At the same time, the values of of pulsars near upper envelope of positive branch may be considered an upper limit of .
Indeed, due to existence of the secular trend , the amplitude of the variational term can not exceed the values of ’s of the upper envelope of the positive branch (for a given ). For pulsars there , while for ones on the lower envelope of the negative branch .
We used the BoundFit generalized leastsquares code developed by Cardiel (2009) to acquire an analytical expressions for the mentioned envelopes of the negative and positive branches on the diagram, and found them to be
(15) 
and
(16) 
respectively. They are shown as dashed lines in Figure 1.
Since we consider pulsars on fits (15) and (16) to be observed with amplitude values of frequency second derivative, their ’s are close to turning points, where corresponding change signs during the oscillations (see also the sketch in Figure 2), and are obviously small (). Hence, the assumption of in equations (15) and (16) is justified.
Due to the absence of pulsars with positive , the amplitudes of variations of frequency first derivatives are likely to be less than their secular values: and . Also, for any more or less stable cyclic process, the following relations should be true:
(17) 
where is a characteristic frequency of cyclic variations. (For the special case (10) these relations are exact). For the upper limit while for the lower one . Hence, one can estimate:
(18) 
which leads to upper limit on the timescale :
(19) 
where is normalized to s. The value of changes from years for oldest to years for youngest pulsars. Then, the value of lower limit should depend on the upper limit of and typical minimal relative amplitude : . Hence
(20) 
or
(21) 
However, should be more than few times of typical observational span length for the pulsars under investigation, which is years. Hence,
(22) 
and is unlikely less than . Moreover, we will show below that is even greater than .
At the same time, assuming that fit to negative branch of diagram
(23) 
represents in the same way a typical variational amplitude (shown as solid line in Figure 1), one may, in a similar way, estimate a typical timescale able to provide observed spread. So, or
(24) 
with spread from to years for oldest and youngest objects, respectively.
Finally, according to equation (17) the amplitudes of the pulsars’ frequency variations should be roughly times the ones of its derivative and times – its second one. For typical maximal timescale (24):
(25) 
These values do not exceed a few Hertz even for young pulsars and significantly less than observed frequencies for most of pulsars. It justifies the neglection of variational term in the rotational frequency decomposition (7):
(26) 
3.5 Pulsar observed braking indices and characteristic ages. An analysis of the diagram.
There are number of widely accepted estimators of pulsars’ physical properties based on results of timing solution and powerlaw (3) spindown model with constant and : characteristic age , rotational energy loss rate and surface magnetic field
(27) 
But, if the amplitude of frequency derivative variations, , is large, these quantities may be significantly biased, as they are computed using observed values, and not actual evolutionary ones .
Using phenomenology introduced earlier and assuming powerlaw spindown defined by equation (3) with constant and evolutionary braking index , observed quantities and may be expressed as:
(28) 
and
(29) 
where is pulsar initial period, – its real age, and and are relative and deviations according to equations (8) and (12). If , then years for the typical pulsar and for the older objects.
So, if varies in range, then half of pulsars will be observed as up to times older than their true unbiased characteristic ages, while another half as times younger. The effect is therefore seemingly asymmetric, as characteristic ages will mostly appear overestimated.
The same bias will also appear in the rotational energy loss rate, while magnetic fields will be times underestimated for pulsars with .
Moreover, since observed braking indices , their absolute values will be additionally increased up to order of magnitude for a half of pulsars.
Since and depend on and , the pulsar behaviour on the plane, shown in Figure 7, is similar to that in the one. The positive and negative branches are not identical: the numbers of pulsars and shapes of branches are significantly different.
The most interesting regions of the diagram are the areas I, II and III with the relatively young pulsars. There are 19 pulsars with positive and only 2 with negative in the area III^{6}^{6}6Moreover, these two pulsars with negative (PSR B133862 and PSR B161050), while satisfying to the formal selection criteria described in Sec. 2.2, have relatively bad data sets – B133862 shows significant glitch activity, and time spans for both of them are only about 200 days long.. Binomial test with rejects the hypothesis of an accidental origin of such asymmetry on a 0.01% significance level. The pulsars with positive have measured , which are larger than any reasonable secular value. Moreover, they slightly deviate from the rest of positive branch of the diagram.
Is it possible to explain these as a result of decay of coefficient in spindown law (3), probably due to evolution of magnetic inclination angle (e.g. Davis & Goldstein (1970))? To yield pulsars with both and years, corresponding timescale should be as short as few thousands years only. However, e.g. FaucherGiguère & Kaspi (2006) have shown that there is no significant decay of coefficient on a timescale years for isolated pulsars.
We plotted pulsar evolutionary tracks corresponding to spindown law (3) with and and years on the diagram (blue dashed lines in Figure 7). It is clearly seen that first of them is insufficient to explain of pulsars in area III, while second one is unable to affect observed distribution of pulsars’ at all.
So, we conclude that such asymmetry results from combination of a small enough amplitudes of variations (i.e. , so the negative values of can not be reached yet) and high amplitudes of ones. Assuming and , for these pulsars with from relation (28) one get:
(30) 
Therefore, the existence of such a remarkable group of pulsars on the diagram argues in favour of significantly high relative amplitudes of variations – in decomposition (8).
Young pulsars in the area III are mostly the ones with negative . Since the variational process is likely symmetric relative to evolutionary trend , young pulsars with positive values of should exist with and that is likely less than .
The pulsars of area I seem to be exactly such objects. They all have and their observed characteristic ages are most likely shifted towards lower values. Also, there is a deficit of pulsars in the area II^{7}^{7}7Note, however, that at least six “ordinary” pulsars without measurement are known in this area. Their ages concentrate near 3 and 5 kyr as seen from the set of red vertical strokes in Figure 7 (with kyr) which may be a result of pulsars escaping from there to lower and higher due to positive and negative , correspondingly. So, the variations can be a reason for measured for the youngest pulsars known. We will return to the analysis of the youngest pulsars within our concept of cyclic variations in the Discussion.
Finally, all remaining pulsars on the diagram – objects of the area IV – form two nearly symmetric branches with positive and negative . The variational amplitudes of of these pulsars are much larger than corresponding evolutionary values: . They are also affected by large amplitudes of variations that may produce the oldest (with years) observed pulsars. Indeed, the largest pulsars’ seem to be physically unreasonable, as pulsars should cross their “death lines” in a few millions of years only (Bhattacharya et al., 1992; FaucherGiguère & Kaspi, 2006).
4 Pulsars’ spindown model
In the previous section we performed a statistical analysis of pulsars rotational parameters and concluded that observed pulsar spindown is consistent with an idea of the existence of some cyclic variational process. Modelindependent estimations show that this process operates on timescales as long as few hundreds or thousands of years, while relative amplitudes of variations are comparable to secular, monotonous component of pulsar spindown .
Also, we found some evidence that variations are symmetric in respect to evolutionary term . This fact may be useful to reveal the parameters of the evolutionary component, common for all pulsars.
In the Section below we construct a twocomponent model of observed pulsars’ rotational evolution which consists of evolutionary monotonous and cyclic components according to decomposition (8), and derive its parameters.
4.1 Monotonous component of observed spindown
For the secular spindown term of our model we will assume canonical powerlaw expression:
(31) 
with unique for each pulsar. On the other hand, we presume the constant evolutionary braking index to be the same for all pulsars. The observed quantity and model parameter are not equivalent and are connected with relation (28) for each pulsar.
Combining equations (7), (8) and (31), the evolutionary value of second frequency derivative for each pulsar with observed and may be written either as
(32) 
or as
(33) 
These expressions correspond to projections of an evolutionary trend onto either or planes, respectively. Moreover, when computed using actually measured timing parameters, these quantities may be treated as a statistically independent ones, as they arise from combinations of different observables.
4.2 Cyclic term of observed spindown
As a toy model, not necessarily exact but still representing all the important properties, we assume below a simple harmonic form for an oscillating term of sum (8):
(34) 
where is a constant relative amplitude, and is a phase of variations.
Actual values of for individual pulsars, as well as values of coefficient in the law (31), are unknown. Therefore, we will analyse our subset in terms of its average, ensemble characteristics only. The parameters of ’s distribution will be included in the model.
Namely, to take into account physical diversity of individual pulsars, we assume to be normally distributed with mean and variance :
(35) 
At the same time, the phase is obviously distributed uniformly,
(36) 
as it depends (as a modulus of quantity over ) both on an unknown value of initial phase and randomly selected moment of observations – pulsar’ age . It is also safe to assume the phases of individual pulsars to be uncorrelated both over the ensemble and with pulsar’ parameters.
4.3 Evolution of an ensemble of pulsars
Since we assume to be the same for all pulsars, their trends on the diagram are characterized by individual coefficients and initial frequencies only. These trends are straight lines in logarithmic scale (see Figure 4) and do not intersect each other. Hence, the average behaviour of a pulsars’ ensemble may be characterized by the line separating the set of their evolutionary trajectories in half. Such trend corresponds to the median value of the coefficient :
(38) 
where is a median over an ensemble and relations (8) and (26) are used. The value of can not be estimated from observables only since values of for individual pulsars are unknown, but one can build a distribution of assuming statistical properties of within the ensemble. This distribution is a function of three parameters:
(39) 
It will be used below to build a maximum likelihood estimator for the parameters of the pulsars in the ensemble.
4.4 Criterion for the estimation of model parameters
If the distribution is symmetric with respect to , then the numbers of objects with observed less and greater than corresponding should be nearly equal, both for the whole ensemble and for any subset defined by quantities uncorrelated with phase (i.e. or ). It provides a simple criterion for a model goodnessoffit for an appropriately chosen (defined by equations (32) and (33)) the numbers of pulsars observed with should simultaneously, for both be distributed binomially as
(40) 
where is the total number of pulsars in the ensemble. The same should obviously be valid for any subinterval along or : if the number of pulsars in th subinterval on th plane is , then the probability to have of pulsars with is
(41) 
Using this criterion and the dependence of on the parameters of pulsars’ spindown model, an obvious maximum likelihood estimator can be constructed to derive these parameters.
4.5 Maximum likelihood estimator for the model parameters
The estimation of for any interval of or requires exact values of and, hence, known specific and for each pulsar. These values can not be derived directly from observational data and therefore we have to hypothesize about them.
We adopted the following routine for its estimation. Since no ad hoc information is available on , we assumed it to be a random value distributed according to equations (34)(36). The specific has been, in turn, replaced with median value computed according to relation (38) with these values in mind. This way, the corresponding estimates for the monotonous component of the second frequency derivatives for each pulsar from equations (32)(33),
(42) 
and
(43) 
are in turn random quantities, uncorrelated to each other as they arise from combinations of different observables. They are, however, median estimates of its exact values due to functional form of equations, and therefore should roughly divide an ensemble of pulsars in half for adequately chosen model parameters – , and .
So, to construct maximum likelihood estimator we divided all pulsars into intervals of objects each () along and on the and planes, respectively. Then we calculated a likelihood function as follows:

On the first step specific values of , and have been chosen from the parameter space.

have been determined.

The number of pulsars with in interval have been calculated and

used to compute binomial probabilities defined (41).

Then, the likelihood functions has been constructed:
(44) 
We have also taken into account an absence of pulsars with positive in an ensemble of pulsars. This fact constraints the parameters of distribution, as probability for the pulsar to be observed with positive unlikely exceeds and corresponding probability
(45) 
for the pulsar to be observed with negative therefore has to be large enough.
The binomial probability to observe all pulsars with negative values of has therefore been added as a multiplicative term to the likelihood function, and

finally we calculated logarithmic likelihood functions as
(46) 
The indices = 1 and 2 here mark (statistically independent) likelihoods for the planes and , respectively.
We performed simulations of relative deviations for 297 pulsars times in a row^{8}^{8}8All simulations have been performed using the “Chebyshev” supercomputer of Moscow State University. for each set of given , and . In this way, the distributions (but not exact values) of for every combination of parameters have been obtained. We have computed these distributions over a grid of values , and in the intervals , and with steps , and , respectively.
4.6 Likelihood for the null hypothesis
Classical maximum likelihood technique presuppose a search of extremum of obtained within parameters space and then constraining the confidence regions assuming that distribution of is known if the null hypothesis is valid. Typically is a singlevalued function of model parameters, and the % confidence region consists of a values of for that , where is a random quantity distributed as a logarithmic likelihood within null hypothesis. The minimum of is obviously within this confidence region.
This method may be generalized in a straightforward way for the case when is not a deterministic function of parameters but a random quantity with known distribution.
To do it, we simulated the distribution of by computing its values in a way similar to one used earlier for , but using simulated according to binomial distribution with , and simulated uniformly distributed on the interval. Therefore, the null hypothesis corresponds to the exact symmetry of numbers of pulsars with greater and less than in each of intervals described above.
Then, using known distribution of , we independently estimated the probabilities and conflated them into the resulting united value using the method described by Hill & Miller (2010) and Hill (2011):
(47) 
This relation represents the cumulative probability that both likelihood functions have an extremum in the same point of a model parameters space. Thus, the isocontours of mark the bounds of a percent confidence regions.
4.7 Results of maximum likelihood estimation
The maps of final conflated probabilities are shown in the Fig. 10 as a collection of slices of the parameter space for fixed (upper plots) or (lower plots). The 99% confidence interval covers the range of :
(48) 
The corresponding range of depends on the accepted value of . The decision about the appropriate amplitude and its dispersion can be made from analysis of diagram provided above. Indeed, according to estimation (30), typical relative deviations are as high as , which corresponds to the in Figure 10. So, we believe that the most adequate choice is:
(49) 
and
(50) 
Obtained values of are in a good agreement with the expected from simple magnetodipolar pulsar braking theory with constant (or slowly evolving) (Manchester & Taylor, 1977).
The median evolutionary trend, assuming , relative amplitude and , is shown in Figure 4.
4.8 Timescales of cyclic process
The model of observed pulsars’ spindown does not include explicit relation between variational phase , pulsar age and timescale of variations . To estimate the distribution of , we will formulate it separately. Namely, let
(51) 
where .
Since functions and defined earlier are not fully independent and related through equation (13), the timescales (or frequencies) of a variational process for our model may be estimated directly: assuming secular spindown law (31), , , for each pulsar:
(52) 
To get the distribution of all possible ’s for the observed set we simulated uniformly distributed phases , as well as values of and – according to probability distributions within the confidence regions derived in Section 4.7. Negative values were rejected during this simulation as corresponding to physically impossible combinations of parameters. Distribution of simulated for is shown in Figure 11. Corresponding timescales are clustered inside a thousands of years region.
Also, the Figure 11 displays theoretical distribution for the timescales of NS precession caused by an anomalous braking torque (see Section 5.2 and equation (57)) for comparison. These timescale values depend on a number of physical parameters of NS, such as surface magnetic field, initial period etc. We used distributions of and derived by FaucherGiguère & Kaspi (2006): and distributed normally with , , s and s. Masses and radii of NS have been taken distributed uniformly within and km intervals, respectively. Initial magnetic dipole angles have been chosen randomly from interval.
Both obtained distributions are in a good agreement with each other and consistent with the idea of a longterm thousandsofyears variations of a pulsars’ timing parameters.
5 Discussion
In this work we revise the problem of anomalous values of isolated radiopulsars braking indices. Since individual values of are strongly dominated by nonmonotonous component of the spindown, we performed a statistical analysis of the ensemble of 297 objects. As a result, we concluded that pulsars’ behaviour on the and diagrams suggest the existence of longterm cyclic mechanism affecting observed timing parameters. We constructed a semiphenomenological model of observed pulsar spindown and estimated it parameters.
We found that the timescale of the longterm variations is likely close to few thousands of years – it follows both from simple modelindependent estimation (Section 3.4) and from accurate maximum likelihood analysis (Section 4.8).
In general, our analysis inevitably suggests that observed of radiopulsars vary with large enough relative amplitude. As a result, basic pulsar parameters depending on may be significantly over or underestimated. We already said some words about this effect in Section 3.5, and below we’ll continue its discussion.
5.1 Ages of young pulsars revised
PSR  , yr  B, G  , yr  , G  SNR  SNR age, yr  , s  
B1853+01  0.50  W44  [1]  0.190  
B2334+61  0.46  G114.3+0.3  [2]  0.346  
B175823  0.89  W28  [3,4]  0.391  
[1] Rho et al. (1994)  
[2] Fich (1986)  
[3] Long et al. (1991)  
[4] Kaspi et al. (1993) 
Let’s return to the diagram (Fig. 7). There are eight pulsars in area III with , associated with young supernova remnants. Assuming their are biased mostly due to variations, i.e. , corresponding instant can be derived as
(53) 
Therefore, their evolutionary characteristic ages
(54) 
and magnetic fields
(55) 
where is according to equation (27), can be estimated, which are expected to be closer to the physical ones.
We analysed data on three of these eight supernova remnants that have more or less independent and confident estimations of their ages, and summarized the results in Table 1. In general, corrected pulsars’ ages become more consistent with that of corresponding SNRs in all of these cases. While the estimations of magnetic fields in all cases exceed G. At the same time, if one assume initial periods for these pulsars, then the same corrections of characteristic ages can be provided by the proximity of the and the current period.
A similar correction for the pulsars of area I is not so evident. Whether the measured of the youngest Crablike pulsars are indeed biased from evolutionary due to the same variational process is not clear. However, according to equation (29), pulsars born with nonzero characteristic ages years. If so, then all pulsars in area I likely have greater than at least 2 thousands of years while their measured kyr imply positive .
We believe that understanding of a possibility of significant biasing of standard estimators of pulsars’ ages will induce new insights on the PSRSNR connections, evolution and kinematics.
5.2 Long timescale precession of a neutron star?
And finally we will discuss possible physical nature of a longterm variations of observed pulsars’ spindown suggested in this work.
The unusual timescale of these variations – thousands of years – is much larger than the pulsar rotational period and, on the other hand, much shorter than its typical lifetime. At the same time, even the canonical magnetodipolar braking model seems to already include an essential mechanism for longterm variations.
More than half a century ago Deutsch (1955) derived the vacuum solution for the electromagnetic field of a perfectly conducting, rigidly rotating spherical star. The full braking torque affecting such a star with magnetic moment and spin angular velocity may be described as a superposition of three orthogonal terms (in observer’s frame):
(56) 
where