title: Braking indices of young radio pulsars: theoretical perspective
authors: A. P. Igoshev S. B. Popov
text: 
  
  firstpage–lastpage
  Countdown games, and simulation on (succinct) one-counter nets
      [
      
  ==============================================================
  
  
  
  
  
  Recently, Parthsarathy et al.  analysed long-term timing observations of 85 young radio pulsars. They found that 11
  objects have braking indices ranging ∼ 10-100, far from the classical value n=3. They also noted a mild correlation between measured value of n and characteristic age of a radio pulsar. In this article we systematically analyse possible physical origin of large braking indices. We find that a small fraction of these measurements could be caused by gravitational acceleration from an unseen ultra-wide companion of a pulsar or by precession. Remaining braking indices cannot be explained neither by pulsar obliquity angle evolution, nor by complex high-order multipole structure of the poloidal magnetic field. The most plausible explanation is a decay of the poloidal dipole magnetic field which operates on a time scale ∼ 10^4-10^5 years in some young objects, but has significantly longer time scale in other radio pulsars. This decay can explain both amplitude of measured n and some correlation between n and characteristic age.
  The decay can be caused by either enhanced 
   crystal impurities in the crust of some isolated radio pulsars, or more likely, by enhanced 
  resistivity related to electron scattering off phonons  due to slow cooling of low-mass neutron stars. 
  If this effect is indeed the main cause of the rapid magnetic field decay manifesting as large braking indices, we predict that pulsars with large braking indices are hotter in comparison to those with n≈ 3.
  
  
  
  
  
  stars: neutron – pulsars: general – stars: magnetic field
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Rotational evolution of young radio pulsars is usually described in a simplistic way using so-called magneto-dipole formula <cit.>. This description is a zeroth-order approximation at best because it lacks multiple physical effects: (1) a radio pulsar is not a point source (see  for details), (2) angle between the magnetic dipole  and the rotational axes (so-called obliquity angle α) seems to evolve with time <cit.>, (3) the pulsar magnetosphere is not in vacuum state  and contains charges, so the slowing down proceeds mainly due to emission of pulsar wind <cit.>. Moreover, the poloidal magnetic field, which causes this slowing down cannot be constant due to finite crust conductivity of the neutron star (NS) and should evolve over time.
  
  In response to challenges based on observational data, many models with evolving magnetic field (e.g., see reviews and references in ) and angle <cit.> were developed, as well as several modifications for the equation describing rotational energy losses were proposed (e.g.,  and references therein). These scenarios predict different timing evolution for single radio pulsars. Precise measurements of spin period (P) and its time derivatives (Ṗ, P̈)  provide one of the best tests for various theoretical scenarios.  
  
  Recently, <cit.> measured braking indices, n=2-(P P̈)/(Ṗ^2), for 19 young radio pulsars. 
  <cit.> argued that this result is robust against uncertainties in the irregular (stochastic) component of pulsars' timing noise and reflects a regular pulsar evolution at least at the timescales of tens of years.
  The obtained values for eleven of them are in the range ∼ 10 – 100, two are negative, two are above 100 and only four are compatible with the classical value n=3 predicted by the magneto-dipole formula for constant field and obliquity angle. 
   The original sample of young radio pulsars with decade-long timing measurements used by <cit.> included 85 radio pulsars.
  For most of them only upper limits  on the braking indices (in the range n≲10 – 2000) were obtained. 
  
  
  <cit.> mentioned several  hypotheses to explain large values of the braking index.
  In the present article we test these hypotheses and aim at identifying the physical cause for these measurements.
  First, we examine what fraction of the sample presented by <cit.> could be in fact wide binaries with orbital periods of tens or hundreds years. Second, we systematically analyse physical phenomena occurring in truly isolated radio pulsars which affect the process of slowing down and increase the value of braking index. Such phenomena are the magnetic field decay analysed in Section <ref>, obliquity angle evolution considered in Section <ref>, rigid body precession calculated in Section <ref> or complicated multipolar configuration of the poloidal magnetic field discussed in Section <ref>.  We conclude that only magnetic field decay provides plausible explanation for the majority of large braking indices.
  
  
  
  
  
  
  § WIDE BINARIES
  
  
  A large fraction 
  of NS progenitors (massive OB stars) are expected to be born 
  in binaries <cit.>. Many of these binaries are disrupted at the moment of supernova explosion due to sudden mass loss and pulsar natal kick <cit.>.
  
  <cit.> analysed population of compact objects in supernova remnants to probe the fraction of binaries among them. These authors provide the following estimates: <10 per cent of non-interacting binaries with NSs survive after core-collapse SNae, and so do ≲ 10 per cent of interacting ones. In the case of objects under study, we are definitely not dealing with very close binaries, as otherwise the orbital motion could be directly recognized via pulsar timing <cit.>. Thus, we are possibly left with wide pairs having orbital periods ≳ tens of years. 
  
  NSs born with weak natal kicks could stay bound even in wide binaries with relatively massive companions. 
  If an orbital period exceeds 50-100 years, periodic fluctuations are not seen in the timing measurements of radio pulsar because undisturbed timing series continues for a few decades at best for majority of radio pulsars. Also, many systems might be highly eccentric after an SN explosion, thus a pulsar might spend a long time in apastron, where orbital velocity is very low.
  
  Unseen as periodic modulation in time series, the companion gravity still accelerates the pulsar contributing to measured P̈ and increasing or decreasing the braking index.
  According to <cit.> we are left with small (but non-negligible) fraction of such systems, few per cent at most.  
  <cit.> discussed the possibility that PSR J0857–4424 is a member of a wide binary system with a NS or a WD as a companion. 
  For two pulsars in the sample (PSR  J1637–4642, PSR J1830–105) presence of planetary companions have been discussed <cit.>, but detailed analysis did not confirm this hypothesis. 
  
  
  
  Another way of looking at the <cit.> data is to count number of objects with positive and negative braking indices. If all large braking indices are caused by presence of a companion on a wide orbit, we would expect  random orientation of this companion's orbit in respect to the line of sight and random orbital phase, thus the gravitational acceleration (and braking index) should be positive in approximately half of the cases (p=0.5) and negative in the other half (see eq. 19 in  which depends on the orbital phase ϕ as cosϕ). We plot oscillations of braking index for parameters of PSR J0857-4424 estimated by <cit.> (mass of unseen companion M_c = 1 M_⊙, orbital period 130 years, inclination i=30^∘). As clearly seen from the figure the braking index is indeed negative half of the time.
  
  
  
  We can quantify a probability of finding k=2 negative braking indices in n=19 measurements. This is done using the binomial distribution. 
  
  The resulting probability is:
  
      f(k,n,p) = n!/k! (n-k)! p^k (1-p)^n-k≈ 3× 10^-4,
  
  which is very low. Therefore, it is highly unlikely that all non-standard braking indices could by explained by presence of massive unseen companion on a wide orbit.  
  
  
  Additionally, to even better constrain possible fraction of wide companions, we search for optical counterparts to 19 pulsars with measured braking indices using the second Gaia data release <cit.>. We search for all stars located at an angular separation up to 7 from the catalogue position of the radio pulsar. We find that 12 out of 19 pulsars with measured  braking indices have an optical counterpart in the Gaia catalogue. Most of these stars have G band magnitude ∼ 18–20 and  angular separation  ≈ 2. Neither  parallax measurements for such dim stars are precise nor many pulsars have  parallax measurement, so further work is required to check if any of these optical counterparts are in fact associated with any pulsar. 
  
  With the modern data it is possible to statistically check if a part or all of these counterparts can be a simple projection  of otherwise physically unrelated sources. To do so we add half a degree to the right ascensions of all real pulsars and search for counterparts to these “fake” sources. We find 10 counterparts. If we subtract half a degree from the right ascension of the real pulsars and repeat our search for counterparts we find 6 counterparts. It means that a few counterparts to real pulsars could be the pulsars themselves or their distant companions visible in optical band. Of course the Gaia could still miss multiple sources fainter than 20th magnitude.
  
  It is important to discuss the case of PSR J0857-4424 as it has the largest measured braking index of n=2890± 30. In our search for optical counterparts we find a star Gaia DR2 5331775184393659264 located at angular separation of 044 from the pulsar sky location. The parallax of this star is 3.28± 0.40 mas, which results in distance estimate of ≈ 300 pc while the pulsar J0857-4424 seems to be located at distance d_DM=2.83 kpc or d_DM=1.94 kpc using the electron density model  YMW2017 <cit.> and NE2001 <cit.> respectively. According to <cit.> the pulsar distance should be in range of 0.1-1.9 d_DM in 95 per cent of cases, so 0.1 × 2.83 ≈ 300 pc is compatible with the Gaia distance. In reality, the distance estimated by the dispersion measure falls more frequently outside of this wide interval <cit.>.
  Unfortunately, it is impossible to conclude without further observations if the star seen in the Gaia is a pulsar itself, its companion or a completely unrelated object. A wide orbit companion to PSR J0857-4424 was discussed in <cit.> with no solid conclusion. 
  
  
  Overall, following these multiple lines of reasoning described above we derive that a few objects (our guess is from two to six objects) in the sample by <cit.> of pulsars with large braking index could be objects with a distant companion. This fraction of 4/19≈ 21 per cent to some extent exceeds the stellar statistics expectations <cit.>. Four pulsars with braking indices caused by gravitational acceleration nicely agrees with two negative braking indices and with an excess of optical counterparts seen in the Gaia data. We must conclude that remaining ≈ 9-13 radio pulsars have large braking indices which cannot be explained by their gravitational acceleration. 
  
  
  
  
  
  
  § MAGNETIC FIELD DECAY
  
  
  The process of slowing-down of NSs rotation is typically controlled by large scale poloidal magnetic field with an exception of magnetars in an active state, where the crustal motions produce magnetospheric twist; i.e. external field is not purely poloidal anymore, and contains toroidal component. The twist leads to enhanced spin-down of the NS <cit.>. Such a twist is unstable and the magnetosphere tends to untwist itself quasi-steadily <cit.> or through a magnetic reconnection episodes <cit.>. These reconnection events in overtwisted magnetospheres produce high-energy outbursts. Before a strong flare as the twist grows, spin-down rate  increases and X-ray spectrum becomes harder. This type of behaviour was clearly observed in the case of the hyperflare of SGR 1806-20 in 2004 (see  and references therein). 
  Such activity is  not typical for normal radio pulsars. These sources might have moderate toroidal field in the crust, but this does not significantly influence their rotational evolution. Due to these factors we consider from here on only poloidal magnetic field. 
  
  Magnetic field decay in isolated radio pulsars can be described as an action of three distinct  processes (e.g.,  and references therein) driven by: (1) scatter of electrons off crystal phonons, which we refer to as “phonon resistivity”, (2) Hall evolution, and (3)  resistivity due to crystal impurities in the crust. Correspondingly, one can introduce three timescales: τ_ph, τ_Hall, τ_Q.
  
  The expected timescale for the decay due to phonon resistivity is related to cooling of NS and for the  so-called minimal scenario of thermal evolution (i.e., neglecting direct URCA processes, exotic phases, additional heating due to field decay, etc., see e.g. ; we use in our estimates a particular realization presented by )  the mean timescale seems to be around 400 kyr  <cit.>. 
  
  The Hall time scale depends on the magnetic field of a NS as  <cit.>:
  
      τ_Hall = 4π n_e e L^2/c B_p,
  
  where  c is the speed of light, e is the elementary charge,  B_p is the poloidal dipolar magnetic field at the pole, n_e is the electron number density, and L is length scale related to electric currents structure in the crust. The Hall evolution is dominant in magnetars due to their high magnetic fields.
   For a normal pulsar the Hall time scale is typically about a few Myrs.
   
  The timescale of decay due to crustal impurities seems to be the longest in normal pulsars. It exceeds at least 8 Myr and is compatible with the value even larger than 20 Myr <cit.>.  Note that stable long-term magnetic field decay is not clearly visible in lifelong evolution of radio pulsars <cit.>,[Several authors <cit.> came to different conclusions, as they obtained better fits with inclusion of magnetic field decay. What we want to highlight, is that the evolution of the whole population of normal radio pulsars on time scale ∼ few tens of Myrs can be described without significant continuous field decay distinguishable in present day observations already on a time scale ∼ few ×10^5 yrs.] which probably agrees with weak  dissipation due to impurities, since it operates along the whole life span of an NS. Therefore,  crystal impurities in the crust cannot contribute significantly to field decay in all young pulsars. 
   
   Dissipation due to phonons and  impurities can be joined together to define the Ohmic time scale:
  
  
      τ_Ohm^-1=τ_ph^-1 + τ_Q^-1.
  
  The Hall process is not a dissipative one, only changing the field structure in the crust. Thus, due to this reconfiguration which results in decrease of the length scale L, the rate of dissipation can be significantly enhanced, as τ_Ohm∝ L^2 <cit.>.
  
  Magnetic field decay affects the first and the second period derivatives and therefore, affects both spin-down age and braking index and makes them correlate with each other. Below, we first estimate the timescale of decay and then construct a more sophisticated model comparing it particularly to the measurement presented by <cit.>. 
  
  
  
   §.§ Simple timescale estimate
  
  
  
  The braking index is a well-known characteristic of radio pulsars (see, e.g.,  and references therein). It can be defined as:
  
      n = ΩΩ̈/Ω̇^2=2 - P P̈/Ṗ^2.
  
  
  If for a moment we assume no evolution of the obliquity angle and no changes in inertia moment, we can absorb these values in a constant. So, the pulsar evolution is described as:
  
      P Ṗ = 2/3β B_p^2
  
  where β is:
  
      β = π^2 R_NS^6/Ic^3,
  
  where R_NS is the neutron star radius and I is its moment of inertia. For typical parameters of a NS — R_NS = 10 km, I = 10^45 g cm^2 — the value of β is equal to β=3.7× 10^-40 s G^-2.
  
  If we combine eqs. (<ref>) and (<ref>) together we can get a simple estimate for braking index when the magnetic field evolves:
  
      n = 3 - 2 P/ṖḂ_p/B_p.
  
  
  Note that for increasing magnetic field the braking index can become significantly smaller than 3, and even reach large negative values. This corresponds to an upward movement in the P – Ṗ diagram. We do not consider such scenarios here (see  and  for a description of this possibility and references to earlier papers), as we focus on data from <cit.> for pulsars with n∼10-100.
  
  If we assume that the magnetic field exponentially decays due to Ohmic losses in the crust on timescale τ_Ohm we obtain:
  
      n = 3 + 2 P/Ṗ1/τ_Ohm.
  
  It is interesting to note that as τ_ch = P / 2Ṗ is the spin-down age, then eq. (<ref>) is simplified even further:
  
      n = 3 + 4 τ_ch/τ_Ohm.
  
  If we take typical values for n and τ_ch from <cit.>, we obtain:
  
      τ_Ohm = 10^4 - 10^5 yrs.
  
  It is important to note that despite such rapid evolution is possible in young NSs, it can not proceed for a long time, otherwise the field dissipates too quickly and the number of observed radio pulsar might be significantly reduced.
  
  
  
  
  
   §.§ Approximate model
  
  
  An approximate description for field evolution with the Hall term was presented by <cit.>:
  
      dB_p/dt = - B_p/τ_Ohm (T_crust) - 1/B_0B_p^2/τ_Hall_0.
  
  In this equation τ_Ohm is temperature dependent and τ_Hall_0 corresponds to the initial value of the Hall time scale which depends on the initial magnetic field: τ_Hall_0∝ B_p, 0^-1. Sometimes an analytical solution of this equation B(t) is given for 
  τ_Ohm being constant. Such solution can also be used for varying 
  τ_Ohm, but in this case the solution should be used iteratively on small timescales. Here, instead, we numerically integrate eq. (<ref>). Moreover, we add the Hall attractor <cit.>, so the Hall evolution is essentially turned off after three initial Hall timescales.
  
  
  As mentioned earlier, the actual τ_Ohm is a combination of the following terms: (i) phonon resistivity and (ii) resistivity due to  impurities, see eq. (<ref>).
  The phonon resistivity is sensitive to the crustal temperature, T_crust, as τ_ph∝τ_ph, 0 / T_crust^-2. 
   The τ_ph, 0 also depends on the depth inside the NS crust where  most of the electric current flows. 
  
  In observations of young isolated radio pulsar usually it is only possible to probe the large scale, poloidal magnetic field of NSs. However, the total magnetic field configuration can include small scale poloidal magnetic field components and toroidal crust-confined magnetic field <cit.>. The magnetic field configuration is formed in convective liquid which solidifies at the proto-neutron stage. So the initial conditions (exact depth of electric current) can differ even for different neutron stars. Therefore, in this work we study: (1) whole range of possible resistivities and (2) we use effective quantities instead of exact distributions across the whole crust.  
  
  For example, <cit.> suggested the following equation to estimate  τ_ph:
  
      τ_ph = 2.2 Myr ρ^15/6_14/T_8^2(Y_e/0.05)^5/3(Y_n/0.8)^10/3(f/0.5)^2 (g_14/2.45)^-2=
  
  
      =τ_ph, 0/T_8^2
  
  where ρ_14 is the density inside the NS crust where the most electric current flows in units of 10^14 g cm^-3; T_8 is the deep crust temperature 
  
  in units of 10^8 K; Y_e, Y_n are the number fractions of electrons and neutrons, respectively; f is a factor  accounting for interactions between neutrons; finally, g_14 is the local gravity in units of 10^14 cm s^-2.  
  This estimate is based on the assumption that typical height for the crust layer with current is proportional to typical pressure height. Applying eq. (<ref>) for density range of 2× 10^13 g cm^-3 to 10^14 g cm^-3 we get the estimate for typical timescale of the phonon decay  range from 4× 10^4 years to 2.2 Myr for T_8 = 1. In contrast to our previous works, here we do not turn-off the phonon resistivity when a critical temperature is reached, in agreement with results by <cit.>. 
  
  The value of τ_Q also depends on the crust density. Recent analysis by <cit.> showed that this value is above 
  8 Myr. So it cannot affect our model on timescales of ≈ 1 Myr, unless some NSs have significantly different value of impurity, Q. At  densities around 10^14 g cm^-3 we 
  can deal with the highly-resistive pasta layer <cit.>. This layer can have impurity parameter of Q≈ 100 which allows the magnetic field to decay significantly by the NS age of ≈ 1 Myr. This layer is important in explaining the absence of long-period magnetars. The past layer cannot play significant role in majority of normal radio pulsars evolution because the fast magnetic field evolution is not seen in these objects <cit.>. In our simulations we assume τ_Q = 200 Myr.      
  
  
  
  
  Different magnetic field decay curves are shown in Figure <ref>.
  
  
  
  
  
  
  
  NSs cool down differently depending on their initial mass. Here we use our fits to cooling curves calculated by <cit.> for several masses.   In our calculations we neglect additional heating due to decaying magnetic field, as we are dealing with normal radio pulsars, so the decay time scale is relatively long and magnetic field energy is relatively low. 
  For NSs with the mass 1.32 M_⊙ we use the following fit:
  
      T_crust (t) = 7.65 × 10^8  K(t/1  year)^-0.182exp(-t/857 kyr).
  
  Here T_crust is the temperature below the heat blanket.
  For NS with mass 1.25 M_⊙ we use:
  
      T_crust (t) = 1.78 × 10^9  K(t/1  year)^-0.26exp(-t/1143 kyr)    .
  
  
  For NSs with the lowest mass of 1.1 M_⊙  the internal temperature drops significantly around the age ≈ 30 kyr. This behaviour is not well described by a fit in a similar form.  To obtain a correct description of cooling in this case we interpolate the numerical curve with a cubic spline. 
  
  
  
  
   
   
  
  
  In order to describe behaviour of braking index and spin-down age, we solve the differential eq. (<ref>) together with eq. (<ref>) 
  with the initial condition P(0) = P_0. We estimate the spin-down age and braking index at each time interval using eq. (<ref>). 
  
  
  
  
  
  We show the dependence of braking indices on the spin-down ages in Figures <ref> and <ref>. In the first plot we variate the initial period and typical phonon decay timescale, but keep the neutron star mass fixed at M=1.32 M_⊙. In Fig. <ref> we variate the NS mass (i.e., we modify the cooling curve, respectively) but keep the phonon decay timescale fixed at τ_ph,0 = 10^5 yrs and initial period fixed at P_0 = 0.04 s. 
  
  We see that n(τ_ch) curves reach values smaller than n≈ 70 before τ_ch≈1 Myr if M_NS = 1.32 M_⊙ even if unrealistically small τ_ph,0 are selected. On the other hand, in Figure <ref> we see that the curve for M=1.1 M_⊙ can explain large braking indices even if τ_ph,0 = 10^5 yrs which is close to the realistic estimate. This result is in correspondence with findings in <cit.>, where we analysed pulsars with few×10^4≲τ_ch<10^6 yrs and demonstrated that P-Ṗ data require field decay with  τ∼few×10^5 yrs which is switched-off for older objects. 
  
  In the set of cooling curves that we use, the low-mass NSs cool down much slower than massive NSs. This is a general feature of many models of thermal evolution of NSs <cit.>. It happens because of much weaker neutrino emission as  bremsstrahlung in neutron-neutron collisions dominates in such stars down to their cores. Other much more efficient cooling mechanisms, even modified URCA process, are damped by proton superfluidity, and they become available only at higher densities, i.e., in more massive NSs <cit.>. Therefore, internal temperature stays above 10^8 K for a  much longer period of time which enhances the phonon resistivity. Phonon resistivity  results in significant decay of the magnetic field (up to a few orders of magnitude) which strongly increases the braking index. Thus, we can potentially “weigh” NSs by timing technique.
  Of course, precise “weighting” requires much better knowledge of thermal and field evolution, as well as a better model of spin-down is necessary. Without such information we cannot determine precisely the masses necessary to explain a large observed value of the braking index.  Still, an optimistic conclusion is that we can obtain a model-dependent  guess of the fraction of low-mass NSs.
  
  It is unclear whether there are enough low-mass young pulsars in the total radio pulsar population to explain observed braking indices.  Formation of such low-mass NSs was discussed in different scenarios, for example, in the ultra-stripped supernova model <cit.>.  <cit.> analysed masses of millisecond radio pulsars which went through an episode of accretion.
  
  If we directly use their bimodal model, we find that 1.1 M_⊙ is 4.5σ away from the peak, so a negligible number of radio pulsars could have this or lower mass. If we  use their unimodal model, we find that 1.1 M_⊙ is 1.7σ away from the peak, so up to 10 per cent of radio pulsars could be born with this or smaller mass. 
  If the model of <cit.> is applicable only for millisecond pulsars,  the isolated radio pulsar with M_NS≈ 1.1 M_⊙ might be more frequent. In particular, if a typical millisecond pulsar accretes ≈ 0.1 M_⊙, the peak of the mass distribution by <cit.> will shift to smaller masses. 
  
  Note that the population of millisecond pulsars might be different from isolated NSs due to selection effects as well. Thus, considerations about the mass distribution of isolated NSs may be just partly based on the data on recycled pulsars.
  
  The main test for the enhanced phonon resistivity hypothesis is to measure NSs temperatures. This test is challenging for young radio pulsars since their X-ray emission is often dominated by bright non-thermal component. NSs with large braking indices are the ones which experience fast magnetic field decay, momentarily, so  even without accounting for additional heating, which, as we stated above, can be neglected in our model, their bulk surface thermal luminosity (and blackbody temperatures respectively) should be systematically larger than thermal luminosities and temperatures of pulsars with braking index ≈ 3 and comparable spin-down ages,  as the later have systematically large masses. If fast magnetic field decay is indeed associated with the lowest mass neutron stars, these NSs typically have largest radii, additionally enhancing their surface thermal luminosity.
  
  
  
   §.§ Crystal impurities
  
  
  If  impurities in the crust are large in some normal radio pulsars and reach values of Q = 20-200 (similar to magnetars ), it also can cause braking indices n≈ 10-100. In comparison to the phonon decay, decay due to impurities is not terminated when the NS is getting cold enough. Therefore  radio pulsar activity of such NSs might cease on a  timescale not much longer than the one of the Ohmic decay, ∼a few × 10^5 years. It seems impossible for a large group of radio pulsars to have impurities Q≈ 100, because the magnetic field would decay before pulsars reach P≈ 1 s and this decay should be easy to notice in analysis similar to <cit.>.  
  
  
  If one assumes that large values of Q are responsible for large braking indices of some NSs, then these compact objects must rapidly loose their fields, and on time scale ∼ 10^6 yrs we expect B∼ 10^8 G.  This possibility can also be tested with high-mass X-ray binaries (HMXBs). NSs in such systems spend ∼ few Myrs at the Ejector and Propeller stages prior to the moment when accretion from the secondary companion is started  (see  for the stages description). In addition, accretion rate from stellar wind is low, thus, in ∼1 Myr just ∼10^-4  M_⊙ is accreted, which is smaller than the mass of the crust. Thus, magnetic field of such NSs have enough time to evolve down to low values in the regime equivalent or similar to the one in normal radio pulsars, if decay is driven by large impurities.  HMXBs with low-field NSs can be selected by short spin periods, as accreting compact objects are expected to reach approximate equilibrium between spin-up and spin-down processes. Analysis of known sources of this type demonstrates that there is no significant fraction of low-field NSs in HMXBs (see  and references therein). We conclude that  impurities in the crust are not a viable explanation for large values of braking index of noticeable fraction
  of young pulsars.
  
  
  
  
  
  § EVOLUTION OF OBLIQUITY ANGLE
  
  
  
  The angle between directions of the magnetic dipole and rotational axis of the radio pulsar is expected to evolve due to  electromagnetic forces caused by the magnetosphere (see, e.g,  for a recent discussion and references to previous studies).  Earlier analytical studies, which considered vacuum magnetospheres <cit.> predicted a fast alignment of the rotational axis with the dipole one. Such a configuration does not slow its rotation due to the electromagnetic emission and can reach arbitrary large values of the braking index. 
  
  First analytical approaches in analysis of the plasma-filled magnetosphere <cit.> predicted that the evolution proceeded towards the counter-alignment.  
  Modern detailed MHD studies of the filled magnetosphere show alignment of the magnetic dipole orientation and the rotational axis <cit.>. This alignment proceeds much more gradually than in the case of vacuum magnetosphere. The configuration with small obliquity angle α  still slows down effectively, which agrees with observations. Therefore, in plasma filled magnetosphere it is hard to reach large braking indices. Below we perform basic calculations to demonstrate quantitatively validity of these statements.
  
  
  
   §.§ Vacuum magnetosphere
  
  
  In a theoretical limit of empty magnetosphere <cit.> the pulsar slows down as:
  
      P Ṗ = 2/3β B_p^2 sin^2 α    .
  
  The obliquity angle evolves with time as:
  
      sinα (t) = sinα_0 exp( - t/τ_v),
  
  where α_0 is the initial obliquity angle and τ_v is the evolution timescale which is equal to:
  
      τ_v = 3/2P_0^2/cos^2α_01/β B_p^2   .
  
  Note, that in this model a priori there are no specific restrictions on the value of τ_v because it is very sensitive to cosα_0, and from observations the distribution of initial obliquity angle is not well constrained. 
  If we follow the procedure described in Section <ref> we obtain the following estimate:
  
      n = 3 + 4 τ_ch/τ_v    .
  
  To explain the large braking indices we require τ_v≈ 10^4-10^5 yrs, which is not consistent with observational data. Observations indicate that the aligment proceeds on timescale of ∼ 7× 10^7 yrs <cit.>. Additionally, the vacuum magnetosphere struggles to explain formation of radio pulsars with periods of ≈ 1 sec because aligned rotators in this model do not spin down significantly <cit.>. 
  
  
  
  
  
   §.§ Plasma-filled magnetosphere
  
  
  To compute the spin evolution of a pulsar with plasma-filled magnetosphere we solve simultaneously two equations:
  
      PṖ = (κ_0 + κ_1 sin^2 α) β B_p^2
  
  and
  
      dα/dt = -κ_2 βB_p^2/P^2sinαcosα   .
  
  Here values of constants κ_0 ≈ 1, κ_1 ≈ 1.2 and κ_2 ≈ 1  are estimated via MHD simulations by <cit.> and α is the obliquity angle, as above. These equations are solved using the fourth-order Runge-Kutta integration method with initial conditions P(0) = P_0 and α (0) = α_0. 
  As soon as the period, the period derivative, and the obliquity angle are computed we calculate the braking index using:
  
      n = 3 - 2 κ_1 sin 2 α/(κ_0 + κ_1 sin^2 α)τ_chdα/dt    .
  
  This can be simplified further:
  
      n = 3 + 2 κ_1κ_2 sin^2αcos^2α/(κ_0 + κ_1sin^2α)^2
  
  From this equation we can see <cit.> that the braking index in plasma filled magnetosphere varies between 3.0 and 3.25. The results of our calculations for both vacuum magnetosphere and plasma-filled magnetosphere are shown in Figure <ref> for P_0 = 0.01 s and α_0 = 70^∘ (left panel). We see that the obliquity angle reaches values which are very close to 0^∘ for vacuum magnetosphere. In this case the braking behaves as sin^2 α which means that pulsars with obliquity angle close to 0 are not slowed down and the braking index can reach arbitrary large values, see Figure <ref> (right panel; for vacuum magnetosphere we use P_0 = 0.04 s). The presence of plasma in the pulsar magnetosphere changes this situation. In this case the angle evolves much slower. Also a pulsar with obliquity angle exactly 0^∘ is still slowed down and has braking index close to 3. The maximum braking index for the plasma filled magnetoshpere showed in Figure <ref> has a value of 3.125.
  
  
  
  
  Overall, a pulsar with vacuum magnetosphere can reach arbitrary large braking indices at any moment of its evolution. The exact value is determined mostly by the initial obliquity angle. Vacuum magnetosphere is not considered to be physical because it contradicts observations. On the other hand the evolution of obliquity angle following the detailed MHD simulations shows that the braking index is always in range n=3 – 3.25 independent of the initial obliquity angle. Therefore, we do not think that obliquity angle evolution can explain braking indices measured by <cit.>.
  
  
  
  
  
  § PRECESSION IN THE CASE OF PLASMA FILLED MAGNETOSPHERE
  
  
  
  If the tensor describing the moment of inertia for radio pulsar Î is not isotropic i.e. the mass distribution is not spherically symmetric and the rotation axis does not coincide with the axis of symmetry, the pulsar is expected to precess <cit.>. Precession means that the obliquity angle α evolves towards alignments as in the plasma filled magnetosphere model, but also oscillates <cit.>. 
  These oscillations change Ṗ and P̈ altering the braking index significantly on timescales of 100-1000 years.
  
  Here we use following prescription <cit.>: we project the angular rotational velocity Ω = 2π / P onto principal axes of the radio pulsar e⃗_1, e⃗_2, e⃗_3 using angles χ and θ. Angle χ is the angle between dipole moment and axis e⃗_3. The angle θ is between the rotational axis and e⃗_3. In this representation, the inertia moment has components I_1, I_2, I_3. We simplify the problem and assume that I_1 = I, I_2 = I (1 + ϵ_12) and  I_3 = I (1 + ϵ_13). We simplify the problem further and consider a bi-axial star with ϵ_12 = 0.
   As soon as we choose initial conditions α_0, χ_0 and θ_0 we calculate the component of the angular velocity as:
  
      ω_3 = 2π/Pcosθ_0 ,
  
  
      ω_1 = 2π/Pcosα_0 - cosθ_0 cosχ_0/sinχ_0,
  
  
      ω_2 = √((2π/P)^2 - ω_3^2 - ω_1^2).
  
  
  These values evolve further under influence of magnetospheric torques K⃗ following Euler's equations:
  
      Î dω⃗/dt + ω⃗× (Îω⃗) = K⃗    .
  
  This equation can be written in components as:
  
      dω_1/dt = - ϵω_2 ω_3 + K_1/I,
  
  
      dω_2/dt = ϵω_3 ω_1 + K_2/I ,
  
  
      dω_3/dt =  K_3/I (1+ϵ)   .
  
  The torques are taken from the MHD simulations by <cit.>:
  
      K_align = B_p^2 Ω^3 R^6_NS/0.25 c^3    ,
  
  
      K_x = κ_2 K_alignsinαcosα   ,
  
  
      K_y = κ_3 K_alignc/Ω R_NSsinαcosα    ,
  
  
      K_z = - K_align (κ_0 + κ_1 sin^2 α)
  
  with numerical coeffitients κ_0 = 1, κ_1 = 1.2, κ_2 = 1 and κ_3 = 0.1. The vector of the torque is transformed using matrix K⃗_123 = ÂK⃗_xyz to the coordinate system of the rotating body, see appendix 1 in <cit.>.
  We solve the system of ordinary differential eqs. (<ref>-<ref>) using the fourth order Runge-Kutta integrator. At each numerical step we compute the obliquity angle as:
  
      α = arccos(ω_1 sinχ + ω_3 cosχ/Ω)  .
  
  We perform a numerical simulation with θ_0 = 60^∘, χ = 4^∘, α=60^∘ and ϵ = 8× 10^-12.  We select parameters similar to <cit.> which they used to describe timing behaviour of the Crab radio pulsar with exception of χ, which we increased slightly.
  
  
  
  We show results in Figure <ref>. It is easy to see that small oscillations of the obliquity angle causes the braking index to deviate significantly from n=3. These deviations are periodic with approximate period of:
  
      T_prec = P/ϵ≈ 130 yrs    .
  
  Half of the time the braking index is expected to be negative. It is seen in Figure <ref> (right panel) that a correlation is present between spin-down age and braking index as it requires some time for the precession to fully develop. Increase in the initial angle χ could produce even stronger modulation in the braking index.
  
  Overall, the hypothesis that measured braking indices could be explained by rigid body precession has the same problems as the wide binary hypothesis: there are not enough negative braking indices. 
  Hence, a fraction of the total sample might have non-standard braking indices because of precession.
  If the extreme measured braking indices are indeed caused by pulsar precession this is easy to check in radio observations. The precession causes changes in the pulsar profile and the inclination geometry (probed though radio polarisation) should evolve with time. 
  
  
  
  
  
  § NON-DIPOLAR MAGNETIC FIELD CONFIGURATIONS
  
  
  
  The braking index depends on configuration of the poloidal magnetic field of a NS. The braking index is most sensitive to the component of the field which stays dominant at the light-cylinder radius. <cit.> showed that the braking index depends on the dominant magnetic multipole moment l as n=2l + 1.
  He also provided an equation to estimate the mean surface strength of non-dipolar component if it dominates:
  
      ⟨|B_lm|⟩ = S_lm(2l)!/2^l l!√(Ic^2l+1/2^2l-3π^2l-1)1/R_NS^l+2√(Ṗ(P/m)^2l-1),
  
  where S_lm is a numerical coefficient which is close to 1, exact values can be found in table 1 of <cit.>.
  
  At first, the dominant multipolar structure seems a plausible explanation 
  for certain pulsars with not too large braking index. But if we develop this hypothesis further we notice the following.
   If we assume that pulsar J1815-1738 with braking index n=9± 3 has dominant octupole poloidal magnetic field l=4, it means that surface strength needs to be ≈ 10^25 G to explain the combination of period and period derivatives, which is clearly nonphysical. The global magnetic field configuration 
   with pure hexapolar structure  (l=3) would require mean surface magnetic field ≈ 10^20 G and could produce braking index of n=7. This strength is still nonphysical because the magnetic pressure would exceed gravity. 
   
  It is worth to notice that pulsars with large braking indices do not look as outliers in the P – Ṗ diagram and do not demonstrate specific (for magnetars) types of activity. So it is hard to expect that these pulsars have much stronger fields. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSIONS
  
  
  In this article we analyse physical grounds for non-standard 
  braking indices (), n∼ 10 – 100, of some young radio pulsars. At first, we restrict a number of possible wide binaries in the sample by <cit.> based on the stellar statistics, a number of negative braking indices, and excess of optical counterparts seen in the Gaia second data release. We think that the majority of large braking indices (11/15) are intrinsic, thus they are not caused by  gravitational acceleration due to a stellar companion on a wide orbit. For these pulsars we investigated four hypotheses explaining  origin of large braking indices: (1) magnetic field decay, (2) evolution of the obliquity angle, (3) precession and (4) complicated multipole structure of the poloidal magnetic field.
  
  We find that magnetic field decay is the only plausible explanation for the majority of large braking indices. The evolution of obliquity angle can cause large n for certain initial obliquity angles only in the nonphysical case of vacuum magnetosphere. Plasma-filled magnetosphere gives n ranges 3 – 3.25 for all initial obliquity angles. Although  large multipoles l=3,4 can explain unusual braking indices of some objects, these surface fields need to have strength well in excess of physical limits for a NS stability. The precession can explain measured braking indices, but predicts equal number of positive and negative braking indices which is not observed.
  
  Magnetic field decay can proceed with different speed in different NSs depending on crust composition (crust impurity parameter Q) or cooling of NSs. 
  
  If in a large fraction of NSs the URCA processes are suppressed, the poloidal magnetic field could decay fast due to large phonon resistivity causing the braking index to reach values n=10-100. 
  
  
  Overall, if the method suggested by <cit.> does not discriminate against negative braking indices, the excess of large positive braking indices is statistically significant. Large positive braking indices can only be explained by magnetic field decay (possibly an episode of decay). Hence, it is necessary to study the phonon decay in details. Necessary prediction in this case is that large braking indices pulsars are hotter than normal braking indices objects. 
  If the method of <cit.> discriminates against negative braking indices, measurements of <cit.> can be explained by a mixture of (1) presence of unseen companion (only a fraction of whole sample because of the stellar statistics) and (2) precession. Precession has simple observational signatures: evolution of pulsar profile and polarisation angle with time. Therefore, further observations of pulsars from the list of <cit.> will be able to discriminate unambigiosly from these possibilities.
  
  
  
  
  § ACKNOWLEDGEMENTS
  
  We are grateful to the anonymous referee for comments and suggestions that helped to improve the manuscript.
  SP thanks Peter Shternin for comments on cooling curves.
  AI thanks support from the Science and Technology Facilities Council for research grant ST/S000275/1. SP acknowledges the support from the Program of development of M.V.Lomonosov Moscow State University (Leading Scientific School `Physics of stars, relativistic objects, and galaxies'). AI thanks A. Frantsuzova for help with proofreading of the manuscript.
  
  This work has made use of data from the European Space Agency (ESA) mission
  Gaia (<https://www.cosmos.esa.int/gaia>), processed by the Gaia
  Data Processing and Analysis Consortium (DPAC,
  <https://www.cosmos.esa.int/web/gaia/dpac/consortium>). Funding for the DPAC
  has been provided by national institutions, in particular the institutions
  participating in the Gaia Multilateral Agreement.   
  
  
  
  
  § DATA AVAILABILITY
  
  The data underlying this article will be shared on reasonable request to the corresponding author.
  
  
  
  
  
  
  
  
  
  
  
  mnras
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  