The effective magnetic field decay of radio pulsars: insights from the statistical properties of their spin frequency's second derivatives

TL;DR

This study introduces a new statistical method to analyze the decay of magnetic fields in radio pulsars, revealing that the decay timescale exceeds several million years and is independent of specific theoretical models.

Contribution

The paper presents a novel approach using the distribution of second derivatives of pulsar spin frequencies to estimate magnetic field decay timescales without relying on specific evolution theories.

Findings

Magnetic field decay timescale exceeds ~5 Myr for young pulsars.

Decay timescale exceeds ~100 Myr for older pulsars.

The extent of the closed magnetic line region is close to the light cylinder radius.

Abstract

We present a new method to investigate the effective magnetic field decay of isolated neutron stars, from the analysis of the long-term timing data of a large sample of radio pulsars \citep{2010MNRAS.402.1027H}. There are some differences between the distributions of frequency's second derivatives of the pulsar spins with different effective field decay timescales. Kolmogorov-Smirnov tests are performed to reexamine the consistency of distributions of the simulated and reported data for a series values of decay timescales. \textbf{We show that the timescale of the effective field decay exceeds $\sim 5~\rm {Myr}$ for pulsars with spin-down age $\tau_{\rm C}< 10^{7}~{\rm yr}$ or $\sim 100~\rm {Myr}$ for pulsars with $10^{7}<\tau_{\rm C}< 10^{9}~{\rm yr}$ in the sample. The result does not depend on any specific theories of the field evolution, the inclination decay or the variation in the moment of inertia. It is also found that the extent of the closed line region of the magnetic field is close to the light cylinder $r_{\rm lc}$, i.e., the corotating radius $r_{\rm c}\approx r_{\rm lc}$ is a good approximation for the observed pulsar population.

Figures

Equations

$$\dot{\nu }=-K{\nu }^3,$$

$$K=\frac{4{\pi }^2{B}^2{R}^6}{c^{3}I}(1+{\sin }^2\theta ).$$

$$\ddot{\nu }$$

$$\ddot{\nu }={\ddot{\nu }}_{\mathrm{O}}+{\ddot{\nu }}_{\mathrm{SD}}-2\dot{\nu }/{\tau }_{\mathrm{B}}.$$

$$B(t)=B_d(t)(1+\sum k\sin (\varphi +2\pi \frac{t}{T})),$$

$$\dot{\nu }=-AB(t{)}^2{\nu }^3,$$

$$\Phi (t)={\int }_{t_0}^t\nu (t^{\prime })\mathrm{d}{t}^{\prime }.$$

$$\Phi (t_i)={\Phi }_0+\nu (t_i-t_0)+\frac{1}{2}\dot{\nu }(t_i-t_0{)}^2+\frac{1}{6}\ddot{\nu }(t_i-t_0{)}^3.$$

$$4\log B_{\mathrm{p}}-6.5\log P=\mathrm{45.7.}$$

$$\dot{\Omega }=\frac{B^2{R}^6\Omega }{4cI{r}_{\mathrm{c}}^2}[{\sin }^2\theta +(1-\frac{{\Omega }_{\mathrm{death}}}{\Omega }){\cos }^2\theta ],$$

$$\begin{array}{rl}{P}_{\mathrm{death}}=& 8.1^{1/(2-\alpha )}\mathrm{s}(\frac{B}{10^{12}\mathrm{G}}{)}^{1/(2-\alpha )}\\ & \times (\frac{V_{\mathrm{gap}}}{10^{12}\mathrm{V}}{)}^{-1/(2-\alpha )}(\frac{P_0}{1\mathrm{s}}{)}^{-1/(2-\alpha )},\end{array}$$

$${\tau }_{\mathrm{Ohm}}\sim \frac{4\pi \sigma L^2}{c^2}=13.5(\frac{\sigma }{3\times 10^{24}{s}^{-1}})(\frac{L}{\mathrm{km}}{)}^2\mathrm{Myr},$$

$${\tau }_{\mathrm{Hall}}\sim \frac{4\pi e{L}^2{n}_{\mathrm{e}}}{cB}=\frac{16.8}{B_{13}}(\frac{n_{\mathrm{e}}}{2.5\times 10^{36}{\mathrm{cm}}^{-3}})(\frac{L}{\mathrm{km}}{)}^2\mathrm{Myr},$$

$${\tau }_{\mathrm{A}}\equiv \frac{\sin \theta }{\frac{d}{dt}\sin \theta }=2\frac{{\sin }^2\theta }{{\cos }^2\theta }{\tau }_{\mathrm{c}}.$$

Full text

stars: neutron — pulsars: general — stars: magnetic field

††journal: ApJ

1 Introduction

The magnetic field is probably one of the most important physical quantities affecting the evolution and the observational behaviours of radio pulsars. The field strength determines the loss rate of rotational energy, the luminosity of pulses, and thus the spin evolution and observability of a pulsar. The primary method used to determine the magnetic field is by measuring each pulsar’s spin parameters, which have actually provided a surprising amount of information on the nature of the pulsed radio sources. As such, the knowledge of the spin parameters is particularly valuable in elucidating whether magnetic field decay occurs in isolated neutron stars. Many impressive studies have been done on this issue during the past few decades (e.g. Ostriker & Gunn, 1969; Gunn & Ostriker, 1970; Lyne et al., 1982, 1985; Stollman, 1987; Narayan & Ostriker, 1990; Bhattacharya et al., 1992; Harrison et al., 1993; Lorimer et al., 1997; Han, 1997; Cordes & Chernoff, 1998; Tauris & Manchester, 1998; Tauris & Konar, 2001; Gonthier et al., 2002; Faucher-Giguère& Kaspi, 2006; Igoshev & Popov, 2015; Johnston & Karastergiou, 2017). Unfortunately, the conclusions of pulsar population investigations have often been conflicting (See e.g. Harding & Lai (2006); Ridley & Lorimer (2010); Lorimer (2011) for reviews). The lack of conclusive evidence on magnetic field decay is mainly due to the fact that the true age of a pulsar is unavailable, and the characteristic (spin-down) age $\tau_{\rm C}$ is normally significantly different from its true age (e.g. Zhang & Xie, 2011). This makes the evolution of the magnetic field remaining as one of the most important unresolved issues of the physics of neutron stars.

Hobbs et al. (2010) studied the timing noise in the residuals of 366 pulsars that had been regularly observed for 10 to 36 years, and showed that the magnitude of the frequency’s second derivatives of the pulsar spins, i.e. $|\ddot{\nu}|$, is much larger than that caused by magnetic braking of the neutron star, and the numbers of negative and positive $\ddot{\nu}$ are almost equal in the sample. It had also been noticed that the distributions between the positive and negative signs in $\ddot{\nu}-\tau_{\rm C}$ diagram show an approximate symmetry.

In this paper we present a new method to investigate the magnetic field decay of radio pulsars, from the analysis of the long-term timing data. We find that the decay timescale can be constrained by measuring the difference between the distributions of $\ddot{\nu}$ with different decay timescales in $\ddot{\nu}-\tau_{\rm C}$ diagram. This method can effectively avoid the “true age problem”. The method and its validity are described in section 2, the revealed restrictions on the timescale of effective magnetic field decay are shown in section 3, and the magnetospheric effects are also tested in the section. The physical implications of the decay timescales are discussed in section 4, and the results are summarized and discussed in section 5.

2 The method

2.1 The Spin-down Models

A basic model for a pulsar’s spin-down is the magnetic dipole radiation model (Pacini, 1968; Ostriker & Gunn, 1969; Gunn & Ostriker, 1970; Ruderman & Sutherland, 1975). The standard dipole (SD) radiation model assumes that the pure magnetic dipole radiation in vacuum as the braking mechanism, i.e.

[TABLE]

where $\nu$ and $\dot{\nu}$ are the spin frequency and its first derivative, respectively. The parameter $K=8\pi^{2}R^{6}B^{2}\sin^{2}\theta/3c^{3}I$ is a constant, $B=3.2\times 10^{19}\sqrt{-\dot{\nu}/\nu^{3}}/\sin\theta$ is the effective dipole magnetic field at equator, $R~{}(\simeq 10^{6}~{}{\rm cm})$ and $I~{}(\simeq 10^{45}~{}{\rm g~{}cm^{2}})$ are the radius and moment of inertia, respectively. Xu & Qiao (2001) improved the model by incorporating the effect of a longitudinal current outflow (relativistic particle wind) powered by a unipolar generator into the rotation energy-loss rate. This effect was confirmed by a few intermittent pulsars, whose rotation slows down faster when the pulsar is on than when it is off (e.g. Kramer et al., 2006). Further, Spitkovsky (2006) developed a numerical method for evolving time dependent force-free MHD equations and applied it to solving a dynamic pulsar magnetosphere. Similarly, the dynamic magnetosphere (DM) model also included both the dipole radiation and the unipolar generator mechanisms to contribute the total braking torque of an oblique pulsar, and they found a formula that gives a very good fit to the oblique spin-down for all inclinations. The formula have the same form with Eq.(1) but with a different parameter $K$,

[TABLE]

The inferred effective magnetic field at the magnetic equator is then $B=2.6\times 10^{19}\sqrt{-\dot{\nu}/\nu^{3}}(1+\sin^{2}\theta)^{-1/2}$, which can be up to $1.7$ times smaller than the estimate from the SD formula.

2.2 The Frequency’s Second Derivatives and The Revised Spin-down Model

Both the SD and the DM model predict that the frequency second derivative of a pulsar spin $\ddot{\nu}_{\rm SD}=\ddot{\nu}_{\rm DM}=3\dot{\nu}^{2}/\nu>0$. However, it is widely known that the observed $\ddot{\nu}$ for the majority of pulsars cannot be explained by these model with a constant field strength $B$. Particularly, the recent large-sample analysis showed (Hobbs et al., 2010; Zhang & Xie, 2012a, b) that $|\ddot{\nu}|\gg\ddot{\nu}_{\rm SD}$, $\ddot{\nu}^{+}(\tau_{\rm C})\approx-\ddot{\nu}^{-}(\tau_{\rm C})$ and $N(\ddot{\nu}^{+})\approx N(\ddot{\nu}^{-})$, where $N$ indicates the total number, the superscripts ‘+’ and ‘-’ indicate positive and negative signs of $\ddot{\nu}$ and $\tau_{\rm C}\equiv-\nu/2\dot{\nu}$ is the characteristic age of a pulsar. All the pulsars in the sample are shown in the $|\ddot{\nu}|-\tau_{\rm c}$ diagram in panel (1) of Fig. 1, in which 193 pulsars have $\ddot{\nu}>0$ and the remainder 173 pulsars have $\ddot{\nu}<0$.

Following Blandford & Romani (1988), we re-formulate the braking law of a pulsar as $\dot{\nu}=-K(t)\nu^{3}$, which assumes that the SD or DM model is responsible for the instantaneous spin-down of a pulsar, but $K(t)$ is time-dependent. Generically and without depending-upon any specific model for the time-dependence, $K(t)$ can be decomposed into a long-term monotonic term plus a short-term perturbation to the monotonic term. Again assuming $R$ and $I$ are constants, the decomposition is equivalent to $B(t)=B_{\rm M}(t)+B_{\rm O}(t)$, where $B_{\rm M}(t)$ is the long-term monotonic component and $B_{\rm O}(t)$ is the short-term perturbation around $B_{\rm M}$. The quasi-periodic oscillation structures, which is widespread in pulsar timing behaviours (Hobbs et al., 2004, 2010), can be phenomenologically described by $B_{\rm O}(t)$ (Zhang & Xie, 2013; Xie, Zhang, & Liao, 2015). One possible source of the perturbation may be their magnetospheric activities, which are known to influence their timing behaviours significantly (Lyne et al., 2010). Then after some simple algebra, we get

[TABLE]

where $\ddot{\nu}_{\rm M}=2\dot{\nu}\dot{B}_{\rm M}/B_{\rm M}$ and $\ddot{\nu}_{\rm O}=2\dot{\nu}\dot{B}_{\rm O}/B_{\rm M}$. For relatively old pulsars without significant glitch activities, $\dot{\nu}<0$ and $\dot{B}_{\rm M}\leqslant 0$, therefore $\ddot{\nu}_{\rm M}\geqslant 0$. Given the case of a quasi-periodic oscillation, $\ddot{\nu}_{\rm O}$ has a positive or a negative value with almost equal chances (Zhang & Xie, 2012a, b). Observationally, since generically and statistically $\ddot{\nu}^{+}(\tau_{\rm C})\approx-\ddot{\nu}^{-}(\tau_{\rm C})$ and $N(\ddot{\nu}^{+})\approx N(\ddot{\nu}^{-})$ for large number of pulsars, clearly $\ddot{\nu}_{\rm O}$ in Eq. (2.2) dominates the observed statistical properties of $\ddot{\nu}$. On the other hand, $\ddot{\nu}_{\rm SD}>0$ and $\ddot{\nu}_{\rm M}>0$ will cause some differences on the distributions of $\ddot{\nu}$ with different decay timescales, and some asymmetry between the observed $\ddot{\nu}^{+}$ and $\ddot{\nu}^{-}$. These effects might in turn provide clues of long-term magnetic field decay of pulsars, since $\ddot{\nu}_{\rm SD}$ can be calculated from the observed $\nu$ and $\dot{\nu}$.

3 Simulations and Tests

3.1 Monte-Carlo Simulations

We assume that the magnetic fields of pulsars in the sample have a typical decay timescale. Define the timescale as $\tau_{B}\equiv-B_{\rm M}/\dot{B}_{\rm M}$, we have $\ddot{\nu}_{\rm M}=-2\dot{\nu}/\tau_{\rm B}$ and

[TABLE]

We can simulate the distributions of $\ddot{\nu}$ with different $\tau_{B}$, as described below. Our strategy is then to search for the typical $\tau_{B}$ that can maximize the p-value of the Kolmogorov-Smirnov test against the hypothesis that the reported distribution is the same as the simulated distribution in $|\ddot{\nu}|-\tau_{\rm C}$ diagram for pulsars in the sample.

We construct a phenomenological model for the dipole magnetic field evolution of pulsars with a long-term decay modulated by short-term oscillations,

[TABLE]

where $t$ is the pulsar’s age, and $k$, $\phi$, $T$ are the amplitude, phase and period of the oscillation,respectively. $B_{d}(t)=B_{0}\exp(-t/\tau_{\rm B})$, in which $B_{0}$ is the field strength at the age $t_{0}$. Substituting Equation (5) into Equation (1), we get the differential equation describing the the spin frequency evolution of a pulsar as follows

[TABLE]

in which $A=\frac{8\pi^{2}R^{6}\sin\theta^{2}}{3c^{3}I}$ is a constant, $R~{}(\simeq 10^{6}~{}{\rm cm})$, $I~{}(\simeq 10^{45}~{}{\rm g~{}cm^{2}})$, and $\theta~{}(\simeq\pi/2)$ is the radius, moment of inertia, and angle of magnetic inclination of the neutron star, respectively. The constant $A$ and $B(t)$ are actually inseparable in Equation (6), which means that the decay timescales of the effective magnetic fields can be attributed not only to the magnetic field evolution, by also the changes of the inclination angle, or the small changes in the moment of inertia.

In order to model the $\ddot{\nu}$ distributions in $\tau_{\rm C}-|\ddot{\nu}|$ diagram, we first obtain $\nu(t)$ by integrating the pulsar spin-down law described as Equation (6), and the phase

[TABLE]

Then, these observable quantities, $\nu$, $\dot{\nu}$ and $\ddot{\nu}$ can be obtained by fitting the phases to the third order of its Taylor expansion over a time span $T_{\rm s}$,

[TABLE]

We thus get $\nu$, $\dot{\nu}$ and $\ddot{\nu}$ from fitting to Equation (8), with a certain time interval of phases $\Delta T_{\rm int}=10^{6}~{}{\rm s}$.

We assume that $k=10^{-4.6}\sim 10^{-1.9}$ and $T$ follows a uniform random distribution in the range from $0.1\sim 10~{}{\rm yr}$, and the reasons will be shown in the next subsection. It is also assumed that the sample of the phase $\phi$ of the field oscillations uniformly distributed in the range from [math] to $2\pi$. Drawing randomly a data set $\{\nu,\dot{\nu},T_{\rm s}\}$ from the reported sample space (i.e. from Table 1 of Hobbs et al. (2010)), and calculating a corresponding start time $t_{0}$, we can obtain a rotation phase set $\{\Phi(t_{i})\}$ using Equation (7). Then the values of $\nu$, $\dot{\nu}$ and $\ddot{\nu}$ can be obtained by fitting $\{\Phi(t_{i})\}$ to Equation (8). Hence one has each $\{\tau_{\rm c},|\ddot{\nu}|\}$. Repeat this procedure for $N$ times, we will have $N$ ($=283$) data points in the $|\ddot{\nu}|$-$\tau_{\rm c}$ diagram. We exclude all the pulsars which have glitch records from the sample, since a small variation in the moment of inertial due to glitches may have an impact on the overall spin-down evolution, for instance, a change on the amount of superfluid content of the star may impact the braking index (i.e. Ho and Andersson (2012)), and actually the timing noise of younger pulsars can be mainly attributed to glitch recovery (Hobbs et al., 2010). Meanwhile, millisecond pulsars ($\tau_{\rm C}\gtrsim 10^{9}~{}{\rm yr}$ or $\nu>100~{}{\rm s^{-1}}$) are also excluded from the sample, since the characteristic age of millisecond pulsars is highly deceptive, and the millisecond period in these systems reflects the recycling spin-up mechanism rather than secular spin-down evolution.

As examples, we show the simulated $\ddot{\nu}$ distribution with $\tau_{\rm B}=10^{6}~{}{\rm yr}$ in left panel of Fig. 1, the simulated distribution without field decay ($\tau_{\rm B}=10^{15}~{}{\rm yr}$ is taken, which is longer than the age of the universe) in the right panel of Fig. 1, and the reported $\ddot{\nu}$ distribution in both panels. In the left panel, one can see that the simulated data are much more sparse in the lower part than the reported data, especially inside the triangular area surrounded by dashed lines. In the right panel, the two distributions are completely consistent.

Some of the simulated pulsars in the bottom right part of the sample may have turned off as a radio pulsar and have crossed the death line in the $P-\dot{P}$ diagram. In most models, the period $P$ at turnoff depends upon the structure and the magnitude of the neutron star’s surface magnetic field (Chen and Ruderman, 1993; Zhang et al., 2000). Assuming a multipole magnetic field configuration, i.e. the polar cap area is similar to that of the pure dipole field, but with very curved field lines at the surface, and the radius of the curvature $r_{\rm c}\sim R=10^{6}~{}{\rm cm}$, and this field configuration will be discussed in section 4. The theoretical death line of the pulsar is then (Chen and Ruderman, 1993),

[TABLE]

After taking into account this observational effect, about 30 simulated pulsars have been excluded from the simulated sample in the right panel of Fig. 1.

The analysis on the scatter of $\ddot{\nu}$ versus $\tau_{\rm C}$ in Fig. 1 is potentially misleading, since the two quantities may not be entirely independent. We carry out a similar analysis as Lyne et al. (1975) to assure the reader that inherent correlations could not be found by plotting random pairings of pulsars, i.e. the value of $\nu$, $\dot{\nu}$, and $\ddot{\nu}$ are randomly taken from different pulsars in the sample. In this case, there is no clear correlations between $\ddot{\nu}$ and $\tau_{\rm C}$.

3.2 Kolmogorov-Smirnov Tests

We perform two-dimensional Kolmogorov-Smirnov (2DKS) test to reexamine the consistency of distributions of the simulated and reported $\ddot{\nu}$ for a series values of $\tau_{B}$. The 2DKS package111http://www.downloadplex.com/Scripts/Matlab/Development-Tools/two-sample-two-diensional-kolmogorov-smirnov-test¯432625.html (Peacock, 1983) is adopted for the test. Our strategy is then to search for a typical $\tau_{B}$ that can maximize the p-value of the 2DKS test against the hypothesis that the two distributions are consistent for the pulsars in the sample. We let $\tau_{\rm B}$ vary from $10^{5}$ to $10^{8}$ yr. The returned p-values are shown with solid lines in the upper panel of Fig. 2. The p-value $0.1$ is considered as the threshold level with probability $95\%$. From the panel, one can see that the decay timescale can be well constrained with p-values larger than $\sim 0.1$. It is shown apparently that the decay timescale $\tau_{B}\gtrsim 5\times 10^{6}~{}\rm{yr}$. In the bottom panel, we show the $N(\ddot{\nu}^{+})$ and $N(\ddot{\nu}^{-})$ as functions of $\tau_{\rm B}$, giving a constraint $\tau_{\rm B}\gtrsim 10^{5}$ yr for $95\%$ probability, which is much loose than the constraint from 2DKS tests. It should be noticed that the 2DKS test has only an approximate and stochastic p-value in each simulation, thus very intensive tests (with $\log\Delta\tau_{B}=0.01$) were performed, as shown in the upper panel. We also checked the validity of 2DKS tests with one-dimensional Kolmogorov-Smirnov (1DKS) test. A very similar result is obtained with 1DKS but $\tau_{B}\gtrsim 10^{6}\rm{yr}$.

We also performed the 2DKS test only for young pulsars with $\tau_{\rm C}<10^{7}~{}{\rm yr}$, the returned values also give $\tau_{B}\gtrsim 5\times 10^{6}~{}\rm{yr}$, as shown in the upper panel of Fig. 3. However, for the sample of middle-age pulsars with $10^{7}<\tau_{\rm C}<10^{9}~{}{\rm yr}$, the returned values give $\tau_{B}\gtrsim 10^{8}~{}\rm{yr}$, as shown in the bottom panel of Fig. 3. In addition, 23 millisecond recycled pulsars are excluded from our samples, and the number is too small to be tested independently with 2DKS.

It is very important to explore the parameter space of the simulations, especially regarding the dependence on the oscillation period $T$ and magnitude $k$. However, it is found that the method cannot place effective restrictions on $T$. For instance, there is no significant difference between the returned p-values for $T$ distributing uniformly from $5$ to $25~{}{\rm yr}$ or from $0.1$ to $100~{}{\rm yr}$. There are extreme examples of magnetars with torque variations which could be interpreted as a change in the spin-down magnetic field by a generous fraction within months (i.e. Archibald et al. (2015) in 1E 1048.1-5937). Thus, for a wider coverage the latter ($0.1\sim 100~{}{\rm yr}$) is chosen for all the simulations in this paper. For the magnitude $k$, the returned p-values for the lower limit and the upper limit ($k=10^{-k_{2}}\sim 10^{-k_{1}}$) of the power index are shown in panel (a) of Fig. 4. It can easily be seen that $3.9\lesssim k_{1}\lesssim 4.6$ and $1.9\lesssim k_{2}\lesssim 2.5$. Therefore, $k=10^{-4.6}\sim 10^{-1.9}$ is taken for the wider coverage. As an example, the young pulsar PSR ${\rm B1828-11}$ shows correlated shape and spin-down changes (Stairs et al., 2019), and the observed $0.7\%$ variation in $\dot{P}$ implies a fractional change of similar magnitude in the oscillation magnitude $k\simeq 3.5\times 10^{-3}$, which falls well within the limits.

3.3 The magnetospheric effects

Aside from the two prevailing models (SD and DM model), Contopoulos and Spitkovsky (2006) proposed a spin-down formula that takes into account the magnetospheric particle acceleration gaps and the misalignment of magnetic and rotation axes, as well as the mechanism of the magnetic field reconnection around the equatorial extent $r_{\rm c}$ of the closed-line region. This formula can be simply expressed as

[TABLE]

where the corotating region follows the light cylinder as $r_{\rm c}=r_{\rm lc}(\Omega/\Omega_{0})^{\alpha}$, $\Omega_{0}$ is the value of the angular velocity $\Omega$ at pulsar birth, and the parameter $\alpha$ ($0<\alpha<1$) depends on the efficiency of the reconnection around $r_{\rm c}$. If the reconnection is very efficient, $r_{\rm c}\approx r_{\rm lc}$, i.e. $\alpha=0$. However, if the reconnection is very inefficient, then the closed-line region cannot grow, thus $r_{\rm c}\approx{\rm const}.$ (Contopoulos and Spitkovsky, 2006). $\Omega_{\rm death}$($=2\pi/P_{\rm death}$) describes a pulsar is “death”, i.e., the cessation of pulsar emission. $P_{\rm death}$ can be written as

[TABLE]

in which $V_{\rm gap}$ is the gap potential.

We perform Monte Carlo simulations to confront the model with observations. The main procedures of the simulations are the same as in the previous subsections. We assume a lognormal distribution of polar magnetic fields with mean value $\langle\log(B/{\rm G})\rangle$ and standard derivation $\sigma_{\log B}$. Following Contopoulos and Spitkovsky (2006), $V_{\rm gap}=10^{13}~{}{\rm V}$ is taken, and the initial period $P_{0}$ is uniformly distributed between $10~{}{\rm ms}$ and $0.2~{}{\rm s}$. We assume the distribution of inclination angle $\theta$ is also uniform from [math] to $\pi/2$. The returned p-values for $\langle\log(B/{\rm G})\rangle$, $\sigma_{\log B}$ and $\alpha$ are shown in the panels (b), (c) and (d) of Fig. 4, respectively. The results show that $12.05\lesssim\langle\log(B/{\rm G})\rangle\lesssim 12.35$, $\sigma_{\log B}\lesssim 0.50$, and $\alpha\lesssim 0.11$.

The result of no floor for $\sigma_{\log B}$ implies that the distribution width of $\ddot{\nu}$ is determined by the oscillation magnitude $k$, rather than by the distribution width of the magnetic field $B$. The parameter $\alpha\lesssim 0.11$ means that our method cannot prove or rule out the spin-down law, but suggests that the reconnection of the north-south poloidal magnetic field around $r_{\rm lc}$ is very efficient, and the extent of the closed line region is close to the light cylinder. We thus propose that $\alpha\sim 0$ or $r_{\rm c}\approx r_{\rm lc}$, is a good approximation for the observed pulsar population.

4 Physical Implications

The time-dependent behaviors of $B(t)$, and thus the decay timescales of the effective magnetic fields, can be attributed not only to the magnetic field evolution, by also the changes of the inclination angle, or the small changes in the moment of inertia. However, since these tests only prescribe lower limits on the evolution timescales, the results are valid for all the three mechanisms.

Magnetic field are crucial for neutron stars’s activities. Understanding the long-term evolution of neutron stars’ magnetic fields might be key to unifying the observational diversity of isolated neutron stars (Viganò et al., 2013). The magnetic field in slow-rotating ultra-magnetized neutron stars, so-called magnetars (AXPs and SGRs), is believed to be decay on timescales of $10^{3}-10^{5}$ years (Thompson & Duncan, 1996), since their rotational energy is not sufficient to power the observed emission. The isolated X-ray pulsars with spin periods longer than $12~{}\rm s$ are still rarely observed. However, they are not subject to physical limits to the emission mechanism nor observational biases against longer periods. This puzzle could be well understood if their magnetic field is dissipated by one or even two orders of magnitude for $1~{}{\rm Myr}$, which is probably due to a highly resistive layer in the innermost part of the crust of neutron stars (Pons, Viganò, & Rea, 2013). For normal radio pulsars, some population synthesis studies suggest that $\tau_{\rm B}$ must be longer than $10~{}{\rm Myr}$ (Hartman et al., 1997; Regimbau & de Freitas Pacheco, 2001). However, there are also some other studies claimed short decay timescales, i.e. $0.1\lesssim\tau_{\rm B}\lesssim 10~{}{\rm Myr}$ (Lyne et al., 1985; Narayan & Ostriker, 1990; Gonthier, Van Guilder, & Harding, 2004; Popov et al., 2010; Gullón et al., 2014; Igoshev & Popov, 2014). The present method imposes a piecewise restriction on the decay timescale, i.e. $\tau_{\rm B}\gtrsim 5~{}{\rm Myr}$ for young pulsars ($\tau_{\rm C}<10^{7}~{}{\rm yr}$), and $\tau_{\rm B}\gtrsim 100~{}{\rm Myr}$ for middle-age pulsars ($10^{7}<\tau_{\rm C}<10^{9}~{}{\rm yr}$), and may contribute to our understanding of actual mechanisms of the field decay and magnetic configurations in neutron stars.

Three avenues for the magnetic field decay in isolated neutron stars have been intensively studied, i.e. Ohmic decay, ambipolar diffusion, and Hall drift (e.g. Goldreich & Reisenegger, 1992; Urpin & Shalybkov, 1999; Geppert & Rheinhardt, 2002; Rheinhardt & Geppert, 2002; Hollerbach & Rüdiger G., 2002; Cumming et al., 2004; Pons & Geppert, 2007, 2010; Pons, Viganò,& Geppert, 2012; Kojima & Kisaka, 2012; Geppert et al., 2013; Gourgouliatos & Cumming, 2014). Depending on the strength of the magnetic fields, each of these processes may dominate the evolution. Ohmic decay occurs in both the fluid core and solid crust. It is inversely proportional to the electric conductivity and independent of the field strength. The Hall drift is non-dissipative and thus cannot be a direct cause of magnetic field decay. However, it can enhance the rate of ohmic dissipation, since only electrons are mobile in the solid crust, and their Hall angle is large. This causes that the evolution of magnetic fields resembles that of vorticity, and then the fields undergo a turbulent cascade terminated by ohmic dissipation at small scales (Goldreich & Reisenegger, 1992; Cumming et al., 2004). Compared with the Hall drift, the timescale of the ambipolar diffusion is much longer for normal pulsars, however, it may be very important for magnetars (Thompson & Duncan, 1996). For a typical density and conductivity profile in the crustal region (e.g. Pons & Geppert, 2007; Gourgouliatos & Cumming, 2014), the Ohmic timescale is

[TABLE]

where $\sigma$ is the electric conductivity, and $L$ is the characteristic length scale of magnetic field in the crust. For the Hall timescale, one reads,

[TABLE]

in which $B_{13}\equiv B/(10^{13}~{}\rm G)$. The combined effect, i.e. Hall cascade, could cause a fast field evolution on a timescale of the order of $10~{}{\rm Myr}$ (Graber et al., 2015).

All these contradictory facts can be well understood by the natural assumption that the magnetic field is maintained by two current systems. The large scale dipolar field which is responsible for the pulsar spin down are supported by long living currents in the superconducting core. Currents in the crust support the small scale multipolar fields which decay on timescale that are comparable to the pulsar spin-down ages (Pons & Geppert, 2007). The two current systems and the corresponding field configurations are particularly demonstrated in the burst activities of a low dipole magnetic field magnetar, SGR 0418+5729, which is expected to harbor a sufficiently intense internal toroidal component (Rea et al., 2010). The present result for middle-age pulsars, i.e. $\tau_{B}\gtrsim 100~{}\rm{Myr}$, suggests that the dipole component that anchored in the crust are relatively low, and thus its decay has no observable influence on the spin frequency’s second derivatives of pulsars in the sample. In addition, the core-anchored field could be expelled and subsequently dissipated in the crust, and our result also implys that the timescale exceeds $5~{}\rm{Myr}$. This may be helpful to understand the poorly known physics at the crust-core boundary.

Our results are also suitable for changes of the inclination angle, which could be either due to rotation-magnetic axis alignment or three-dimensional magnetic field evolution (Philippov et al., 2014; Gourgouliatos and Hollerbach, 2018). Using polarization data for a large number of isolated pulsars, Tauris & Manchester (1998) found that the magnetic beam axis align with the spin axis on a timescale of $\sim 10~{}{\rm Myr}$. With new data, Young et al. (2010) found a shorter alignment timescale of $\sim 1~{}{\rm Myr}$. Theoretically, the electromagnetic torque which brakes the rotation of a pulsar also tends to align the magnetic axis with the rotation axis (Davis and Goldstein, 1970; Goldreich, 1970). The electromagnetic alignment timescale is related to the spin-down age as (Lander and Jones, 2018),

[TABLE]

For young or middle-age pulsars, our results imply the alignment timescale is most likely longer than $\sim 5~{}{\rm Myr}$ or $\sim 100~{}{\rm Myr}$, which is roughly consistent with the relation.

A small variation in the moment of inertia may have an impact on the overall spin-down evolution, for instance, a decrease in the effective moment of inertia due to an increase on the amount of superfluid content as the star cools through neutrino emission may impact the braking index (Ho and Andersson, 2012). However, most of the stars are typically young and glitching pulsars, which have been excluded from our sample. Our results imply that the populations without glitch record shows no long-term variation in the moment of inertia with timescale short than $5~{}{\rm Myr}$.

5 Summary and Discussion

The perturbation from the long-term dipole magnetic field decay will produce some differences on the distributions for the second derivatives of pulsars’ spin frequency with different decay timescales. This in turn provides a new method to investigate the magnetic field decay of radio pulsars, which does not depend on any specific theories of field evolution or inclination decay. We made use of the published large-sample timing data of radio pulsars to find evidence of their magnetic field decay with 2DKS tests. The method impose a piecewise restriction on the decay timescale, i.e. $\tau_{\rm B}\gtrsim 5~{}{\rm Myr}$ for young pulsars with $\tau_{\rm C}<10^{7}~{}{\rm yr}$, and $\tau_{\rm B}\gtrsim 100~{}{\rm Myr}$ for middle-age pulsars with $10^{7}<\tau_{\rm C}<10^{9}~{}{\rm yr}$. It is also proposed that the corotating radius $r_{\rm c}\approx r_{\rm lc}$ is a good approximation for the observed pulsar population. Though pulsars with major glitches have been excluded from the data, tiny glitch activities and other types of timing irregularities may still have some influences on the observed $\ddot{\nu}$, which may cause a small deviation.

We expect to gain much deeper understanding of pulsars from future larger sample of radio pulsars with higher precision data on $\ddot{\nu}$, to be brought by China’s soon-to-be operating Five-hundred-meter Aperture Spherical radio-Telescope (FAST) and the future Square Kilometer Array (SKA).

Acknowledgments

We thank J. Y. Liao for discussions. We thank the anonymous referee for comments and suggestions that led to a significant improvement in this manuscript. This work is supported by National Natural Science Foundation of China under grant Nos. 11603009, 11803009, 11373036,and 11133002, by the National Program on Key Research and Development Project under grant Nos. 2016YFA0400802, by the Key Research Program of Frontier Sciences, CAS, Grant No. QYZDY-SSW-SLH008, and by the Natural Science Foundation of Fujian Province under grant Nos. 2016J05013 and 2018J05006.