The Dipole Magnetic Field and Spin-down Evolutions of The High Braking Index Pulsar PSR J1640-4631

In this work, we interpreted the high braking index of PSR J1640$-$4631 with a combination of the magneto-dipole radiation and dipole magnetic field decay models. By introducing a mean rotation energy conversion coefficient $\overline{\zeta}$, the ratio of the total high-energy photon energy to the total rotation energy loss in the whole life of the pulsar, and combining the pulsar's high-energy and timing observations with reliable nuclear equation of state, we estimate the pulsar's initial spin period, $P_{0}\sim (17-44)$ ms, corresponding to the moment of inertia $I\sim (0.8-2.1)\times 10^{45}$ g cm$^{2}$. Assuming that PSR J1640$-$4631 has experienced a long-term exponential decay of the dipole magnetic field, we calculate the true age $t_{\rm age}$, the effective magnetic field decay timescale $\tau_{D}$, and the initial surface dipole magnetic field at the pole $B_{p}(0)$ of the pulsar to be $(2900-3100)$ yrs, $1.07(2)\times10^{5}$ yrs, and $(1.84-4.20)\times10^{13}$ G, respectively. The measured braking index of $n=3.15(3)$ for PSR J1640$-$4631 is attributed to its long-term dipole magnetic field decay and a low magnetic field decay rate, $dB_{\rm p}/dt\sim -(1.66-3.85)\times10^{8}$ G yr$^{-1}$. Our model can be applied to both the high braking index ($n>3$) and low braking index ($n<3$) pulsars, tested by the future polarization, timing, and high-energy observations of PSR J1640$-$4631.


INTRODUCTION
Pulsars are commonly recognized as highly magnetized and rapidly rotating neutron stars (NSs). A pulsar's secular spindown is mainly caused by its rotational energy losses (Lyne et al. 2015). An important and measurable quantity closely related to a pulsar's rotational evolution is the braking index n, defined by assuming that the star spins down in the light of a power law (Lyne et al. 1993) where Ω andΩ are the angular velocity and its derivative of the star, respectively, and K is a proportionality constant. A standard way to define the braking index is whereΩ is the second derivative of Ω, ν = Ω/2π the spin frequency, and P = 1/ν the spin period (Manchester & Taylor 1977, and references therein). When the magneto-dipole radiation (MDR) solely causes the pulsar to spin down, the braking index is predicted to be n = 3. Since the braking index of a pulsar can provide some information about the pulsar's energy loss mechanisms, it has been investigated by many au-thors (for recent works, see, e.g., Espinoza et al. 2011, Ferdman et al. 2015Franzon et al. 2016;Rogers & Safi-Harb 2017). Until now, only 9 of the ∼ 2500 known pulsars 6 have reliably measured braking indices. Most of them are remarkably lower than 3, demonstrating that the spin-down mechanism is not pure MDR. Recently, Dupays et al. (2008Dupays et al. ( , 2012 put forward an energy loss mechanism, called quantum vacuum friction (QVF), which results from the interaction between the magnetic dipole moment of a pulsar and its induced quantum vacuum. Taking into account the QVF effect, the measured braking indices of n < 3 can be well explained . To account for the observed braking indices, several other interpretations have been put forward (see Gao et al. 2016 for a brief summary).
Recently, Gotthelf et al. (2014) reported the discovery of PSR J1640−4631, associated with the TeV γ−ray source HESS J1640−465 using the N uST AR X-ray observatory.
The possible origin of the braking index of PSR J1640−4631 has been discussed in the literature. The decrease in the magnetic inclination angle α (the angle between the rotational and the magnetic axes) can result in a braking index higher than 3 (Ekşi et al. 2016). A combination of gravitational wave (GW) and magnetic energy dipole loss mechanism radiation could give rise to a braking index 3 < n < 5 (de Araujo et al. 2016 a, b). Alternatively, such a braking index can be accounted for if the magnetic dipole braking still dominates but the star is experiencing a longterm dipole magnetic field decay.
The surface dipole magnetic field has a significant effect on the spin evolution of a pulsar, as well as on its braking index (e.g., Goldreich & Reisenegger 1992;Tauris & Konar 2001;Jones 2009;Marchant et al. 2014). In the case of a varying dipole magnetic field at the pole B p , the braking index n can be simply expressed as follows: where τ c = P/2Ṗ is the characteristic age and τ B = B p /Ḃ p is the timescale on which the dipole component evolves. The above equation shows that any variation in B can lead to a deviation from 3. For a decreasing B p , we will always obtain n > 3 owing to a decaying dipole braking torque (e.g., Rea & Esposito 2011;Kaminker et al. 2014;Potekhin et al. 2015), whereas n < 3 for an increasing B p (e.g., Mereghetti 2008;Gao et al. 2016). The motivation of this work is to account for the high n of PSR J1640−4631 within a theoretical model of a longterm dipole magnetic field decay. The rotational evolution of PSR J1640−4631 is assumed to be dominated by the magnetic field anchored to the solid crust of the star. Since the core field evolves on much longer timescales, the crustal field evolves mainly through ohmic dissipation and the Hall drift. The Hall drift mainly occurs in the neutron-drip solid or out of the core, conserving the magnetic energy (e.g., Haensel et al. 1990;Hertman et al. 1997;Page et al. 2000), while ohmic diffusion occurs in the range of density below the neutron-drip threshold and dissipates magnetic energy through Joule heating (e.g., Pons et al. 2009Pons et al. , 2012Gourgouliatos & Cumming 2015).
The remainder of this work is organized as follows: The key problems in estimating the true age of PSR J1640−4631 are presented in Section 2, and the initial spin period is constrained in Section 3, The initial magnetic field strength, the true age, and the effective dipole field evolution timescale of the pulsar are determined in Section 4, the dipole magnetic field and spin-down evolutions are numerically simulated in Section 5, and alternative dipole magnetic field decay models are presented in Section 6. Comparisons with other works and discussions are given in Section 7.

KEY PROBLEMS IN ESTIMATING T AGE
To investigate the long-term spin evolution of a pulsar, it is best to know the true age of the pulsar. It is well known that τ c is a poor age approximation, and it coincides with the true age only if the current spin period is much longer than the initial value and the braking torque and the moment of inertia I are constant during the entire pulsar life. A pulsar's true age t age can be estimated by the age of its associated supernova remanant (SNR), because it is universally considered that pulsars originate from supernova explosions (e.g., Kouveliotou et al. 1994;Gaensler et al. 2001;Vink & Kuiper 2006).
The SNR ages are currently estimated based on the measurements of the shock velocities, the X-ray temperatures, and/or other quantities. Assuming that an SNR is in the Sedov phase 7 (Sedov 1959), the SNR age can be estimated by a 7 Sedov (1959) showed that there is a self-similar solution for the adia-convenient expression: where the SNR radius R SNR is in units of pc, and the postshock plasma temperature T is in units of keV (for details see Gao et al. 2016). G338.3−0.0 is a shell-type, 8 ′ diameter SNR (Shaver & Goss 1970;Whiteoak & Green 1996), and spatially relates with HESS J1640−465, which is considered to be the the most luminous γ-ray source in the Galaxy. The X-ray pulsar PSR J1640−4631 was recently discovered within the shell of SNR G338.3−0.0 (Gotthelf et al. 2014). Unfortunately, no X-ray emission was detected from the shell of the SNR (probably due to its large intervening column density, N H ∼ 1.4 × 10 23 cm −2 and low temperature; Lemiere et al. 2009;Castelletti et al. 2011); thus, the true age of PSR J1640−4631 cannot be estimated from Equation (4).
The pulsar's true age can be inferred from Equation (1), where P 0 is the initial spin period of the star and the braking index is assumed to be a constant since its birth (e.g., Manchester et al. 1985;Espinoza et al. 2011;Rogers & Safi-Harb 2017). As discussed in Gao et al. (2016), one indeed cannot derive an analytical expression such as Equation (5) to estimate t age if we consider that either K or n in Equation (1) is a time-dependent quantity. Since the braking index is influenced by the variation of braking torque during the observational intervals, it is necessary to assume a constant mean braking index to calculate the age of a young pulsar by using the above equation ). The mean braking index n of a pulsar is defined as n = n(t)dt/ dt.
As seen from Equation (5), t age is less than τ c if n = 3. In addition, Equation (5) implies an upper limit of the true age, t age = 2τ c /(n − 1), corresponding to (P 0 /P ) n−1 = 0. However, it is difficult to obtain the values of P 0 and n of the pulsar; t age cannot be directly derived from Equation (5).
3. INITIAL SPIN PERIOD OF PSR J1640−4631 The initial spin period P 0 of a pulsar is also very important in our ability to investigate the pulsar's long-term spin evolution. Since one cannot know the value of P 0 of PSR J1640−4631, theoretically estimating its initial spin period becomes urgent. This section is composed of the following two subsections.
3.1. X-ray and γ-ray luminosities As mentioned above, PSR J1640−4631 is located inside G338.1−0.0, the pulsar wind nebula (PWN), which is primarily responsible for the γ-ray emission of HESS J1640−465. batic expansion phase, and firstly applied this solution to estimate the age of SNRs in this phase, i.e., the Sedov age. The Sedov age estimates rely on many assumptions, but ultimately on the SNR size and the SNR temperature. If a SNR is too young (when the expansion is free), or too old (when radiation energy losses brake the expansion), the Sedov expansion is no longer applicable. The detection of the pulsar in the SNR provides longawaited evidence that a PWN powers the γ-ray source HESS J1640−465, and its properties can be used to test the radiation mechanisms (e.g., Reynolds 2008;Slane et al. 2014).
Recently, employing Markov Chain Monte Carlo algorithm, Gotthelf et al. (2014) fitted the parameters of the PWN in G338.3−0.0 and obtained the best-fitted values of free parameters, including the high-energy photon ("photon " in short) fluxes. For the purpose of illustration, in Table 1 we list the photon fluxes of PSR J1634−4731 and its PWN, respectively.
In Table 1, the spectrum flux of F X [2−25 keV] includes the X-ray fluxes of PSR J1634−4731 and its PWN. All the given fluxes are in units of erg cm −2 s −1 . Then, the total luminosity (the currently measured photon luminosity) is calculated as where d 12 is a dimensionless distance in units of 12 kpc.

Mean rotational energy transformation coefficient
In order to estimate P 0 of PSR J1634−4731, we introduce an energy transformation coefficient, ς, which is the ratio of the total photon luminosity, L X,γ , to the spin-down luminosity of the pulsar L spin (L spin = IΩ(t)Ω(t)) From the values of current P andṖ , we get the current luminosity L spin ∼ 4.36 × 10 36 I 45 erg s −1 , where I 45 is the moment of inertia in units of 10 45 g cm 2 . Inserting the values of L X,γ and L spin into Equation (8) yields a general solution for the current coefficient, ζ = 0.067d 2 12 I −1 45 . Recently, Getthelf et al. (2014) obtained the best-fit source distance to G338.3−0.0, d ∼12 kpc. This distance agrees well with the typical supernova explosion scenario (e.g., the order of magnitude of explosion energy is 10 51 erg). In this paper, we adopt the distance of d 12 = 1.
Since ζ is variable with time, in order to estimate the initial spin parameters including P 0 , a constant ζ should be assumed in a specific model. We define a mean rotation energy transformation coefficient ζ, where E X,γ is the total photon energy and E spin is the total rotation energy loss, which can be estimated by for P 0 ≪ P (Tanaka Shuta 2016). For a constant braking index n, the spin-down luminosity L spin (t) evolves in the light of a simple power-law form, where L spin (0) is the initial spin-down luminosity and τ 0 is the initial spin-down timescale of the star, which is defined as with the initial period derivativeṖ 0 . Since PSR J1640−4631 is a pulsar powered by rotational energy loss, the total photon luminosity, L X,γ (t), as well as L spin (t), should decrease with time. For the sake of simplicity, we assume that L X,γ and L spin (t) have the same evolution form. Then the total photon energy from t = 0 to t = t age is calculated as where L X,γ (0) is the initial total photon luminosity and the relations of τ 0 + t age = 2τc n−1 and τ 0 = 2τc n−1 ( P0 P ) n−1 are utilized.
Inserting Equations (10) and (13) and τ c ≃ 3360 yrs into Eq.(9), we obtain the initial spin period It is interesting to note that, although n has been used in the derivation above, the ultimate value of P 0 is free of n. Recently, Abdo et al. (2010Abdo et al. ( , 2013 presented catalogs of highenergy gamma-ray pulsars detected by the Large Area Telescope on the F ermi satellite and estimated the values of ζ for some sources, (here ζ is the ratio of the gamma-ray luminosity to rotational energy loss rate of a pulsar). According to Abdo et al. (2010Abdo et al. ( , 2013, the range of ζ is about (0.01 − 1.0)I −1 45 , but its average value is about 0.08I −1 45 , close to our estimate of ζ = 0.067d 2 12 I −1 45 . In Equation (14), we have adopted fiducial NS parameters: mass M = 1.4 M , radius R=10 km, and the moment of inertia I 45 = 1. Strictly speaking, the true value of ζ is unknown and may be constrained by specific NS nuclear equation of state (EOS).
From Equation (14), it is clear that P 0 is a function of moment of inertia I. In order to investigate how different values of I (or the mass [M ]-equatorial radius [R] relation) modify P 0 , we need to take into account the M − R relation of the uniformly rotating NS and feasible EOSs.

Nuclear EOS and moment of inertia
Very recently, to investigate the possibility that some of soft gamma-ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs) family could be canonical rotation-powered pulsars, Coelho et al. (2017) adopted realistic NS structure parameters (e.g., mass, radius, and moment of inertia) instead of fiducial values and considered the issue related to the high-energy luminosity and the spin-down luminosity of SGRs and AXPs in detail. Their results provide useful constraints on the initial spin period P 0 of PSR J1640−4631.
As we know, the maximum NS mass predicted by EoSs is model dependent (e.g., Dutra et al. 2014;Li et al. 2016;Xia et al. 2016;Zhou et al. 2017). The largest sample of measured NS masses available for analysis is publicly accessible online at http://www.stellarcollapse.org/, and has been updated by J. Lattimer. A similar analysis can be found in Valentim et al. (2011), from which one can get a range of about 1 − 2 M for the observational NS masses. The minimum NS mass, as well as the maximum NS mass, is still a matter of debate. Though the NS masses could be could be less than 1.0 M from the viewpoint of the EOS, it is difficult to explain their formations from the supernova explosion mechanism. Pulsars having mass of ∼ 2.0 M are confirmed by Demorest et al. (2010) and Antoniadis et al. (2013). The typical nuclear matter EOSs obtained fromÖzel & Freire (2016) predict that pulsars have maximum mass larger than 2.0 M , which prefers the APR (Akmal et al. 1998) model and the RMF (relativistic mean-field) model (e.g., Lalazissis et al. 1997;Toki et al. 1995). In the APR model, properties of dense nucleon matter and the structure of NSs are studied using variational chain summation methods and the new Argonne v18 two-nucleon interaction (Av18). However, in this work, we adopt the APR3 model (one version of APR EOS) instead of the RMF model. The main reasons are as follows: 1. The RMF theory is based on effective coupling constants that take the correction effect into consideration. However, these coupling constants are density dependent, and a microscopic theory is needed to calculate them.
2. Akmal et al. (1998) investigated the properties of dense nucleon matter and the structure of NSs and provided an excellent fit to all of the nucleon-nucleon scattering data in the Nijmegen database. The authors not only considered the nonrelativistic calculations with Av18 and Av18+UIX (Urbana IX three-nucleon interaction) models for nuclear forces but also described the relativistic boost interaction model (denoted as δv) with and without three-nucleon interaction (UIX * ).
3. In this work, we prefer the APR3 model, e.g., Av18+δv+UIX * model, which provides a constraint on the maximum NS mass M ≤ 2.2 M . This model has also been preferred and included in recent works (e.g., Ozel & Freire 2016;Zhou et al. 2017), which limit the range of NS masses to < 2.2 M because of the absence of any data to constrain the relation at higher mass.
In order to calculate the moment of inertia I, it is of importance for us to adopt a better M − R − I relation for NSs. In this work, we prefer an improved approximation expression of M − R − I with M ≥ 1 M , provided by Steiner et al. (2016),

. Hall drift and Ohmic decay
Here we restrict ourselves to the dipole magnetic field evolution in the crusts of NSs, ignoring possible effects involving the fluid core. Assuming that ions are locked into a column lattice, the only freely moving charged species in the crust are electrons. In order to reconcile the measured braking index of PSR J1640−4631 with the dipole magnetic field decay, it is necessary to introduce the Hall induction equation describing the evolution of the crust magnetic field: where σ is the electric conductivity parallel to the magnetic field, e ν is the relativistic redshift correction, n e is the electron number density, and e the electron charge. This equation contains two different effects that act on two distinct timescales, which can be estimated as where t Ohm is the ohmic dissipation timescale with a typical value of ∼ 10 6 yr or more (Viganò et al. 2013), t Hall is the Hall drift timescale, and L is a characteristic length scale of variation, which can be taken to be the thickness of the neutron star crust (e.g., Haensel et al. 1990;Pons & Gepport 2007;Pons et al. 2009;Aguilera et al. 2008;Ho 2011). There is a typical value of several ×(10 4 − 10 5 ) yr for t Hall , which is constrained by the magnetic field strength of high magnetic field pulsars and magnetars (Viganò et al. 2013). There are two field configurations: one in which the field is confined in the crust, and another in which the field extends into the core (there could be an electron current flowing through a superconducting core). For the purposes of illustration, we choose a dipole magnetic field B p constrained in the crust and assume a simple exponential decay of B p over time, as done in other works (e.g., Pons et al. 2009;Viganò et al. 2013). This assumption requires that the pulsar has an initial magnetic field at the pole B p (0) (corresponding to t = 0) and that the magnetic field decays at a rate proportional to its strength. τ D is an effective dipole magnetic field timescale, which is defined as The Hall timescale t Hall is universally one order of magnitude lower than the ohmic timescale (e.g., Pons & Gepport 2007;Pons et al. 2009). This would imply that the effective timescale we define in Equation (19) is basically dominated by the Hall timescale. By integrating Equation (18), we get a time-dependent field With respect to two parameters of τ D and B p (0), there are three aspects that need to be addressed explicitly: 1. In general, the effective timescales of pulsars are in the range of about 10 4 − 10 7 yrs (see Muslimov & Page 1995Geppert & Rheinhardt 2002;Lyutikov et al. 2015;Mereghetti et al. 2015 and references therein). However, τ D of PSR J1640−4631 is determined mainly by t Hall . To exactly constrain τ D of the star requires both observational and theoretical investigations.
2. Though there are several methods for roughly measuring the magnetic field strength of a pulsar, such as magneto-hydrodynamic pumping, Zeeman splitting, cyclotron lines, magnetar bursts and etc., a direct estimate of the strength B p (0) from the observations is unavailable.
3. For PSR J1640−4631, the initial dipole magnetic field was inherited from its progenitor star and cannot be estimated from its initial spin parameters via a simple MDR model. In this work, we have assumed a variable dipole magnetic field for the pulsar.

Equations of P (t),Ṗ (t) andP (t)
Assuming that the MDR of a pulsar with a decaying field can be responsible for the high braking index of PSR J1640−4631, we here constrain three parameters B p (0), t age and τ D of the pulsar.
If the magnetic field evolution of PSR J1640−4631 cannot be ignored and the dipole braking still dominates, according to Blandford & Romani (1988), the braking law of the pulsar is reformulated asν where a constant inclination angle α = 90 • is assumed for the sake of simplicity. Integrating Equation (21) gives the spin frequency, where ν 0 is the initial spin frequency of the pulsar. Then we get the spin period, Utilizing the differential method, we get the first-order derivative of the spin periodṖ (t), and the second-order derivative of the spin periodP (t), In the above three equations there are unknown variables of B p (0), t age and τ D , which will be constrained by numerically simulating in the next subsection. (14) and (15) with the APR3 EOS, take t = t age , and insert the current values of P = 206.4 mṡ P = 9.7228 × 10 −13 s s −1 andP = −5.27(13) × 10 −24 s s −2 into Equations (23)− (25), then we obtain the values of τ D , t age and B p (0) given a specific value of P 0 (corresponding to a specific I of the pulsar).

The results of numerical simulations If we combine Equations
By numerically simulating, we make a plot of B p (0) versus I for PSR J1640−4631, as shown in Figure 2. In Figure 2 EoS. Also, we obtain the initial spin period P 0 ∼ (17 − 44) ms, the true age t age ∼ 2800 − 3100 yr, the initial dipole magnetic field strength B p (0) ∼ (1.84 − 4.20) × 10 13 G, and the effective timescale τ D ∼ 1.07(2) × 10 5 yr for the pulsar. As we know, the estimates of the crustal ohmic and Hall timescales represent the choice of microscopic parameters. Recently, Gourgouliatos & Cumming (2014) estimate the average initial value of τ Hall for young pulsars about 100 − 200 kyr. This timescale is close to that given by us.

NUMERICAL SIMULATIONS OF THE DIPOLE MAGNETIC FIELD AND SPIN-DOWN EVOLUTIONS
5.1. The evolution of the braking index After obtaining B p (0), t age and τ D of the pulsar, we then use them in numerically simulating the spin-down evolution, especially accounting for its high braking index.
Combining Equations (23)-(25) with n = 2 − PṖ P 2 , we get From the above equation, it is easily seen that n = 3 if τ D → ∞. In order to investigate the evolution of n, we plot the diagrams of n versus t for the pulsar. In Figure 3, the red solid line, yellow solid line, black dashed line, and blue solid line stand for the predictions of the dipole magnetic field decay model with P 0 = 17, 26, 34 and 44 ms, respectively. The blue dot-dashed line stands for the prediction of the MDR model, in which we assume the NS mass M = 1.4 M , radius R = 10 km and momentum of inertial, I = 10 45 g cm 2 , corresponding to P 0 ∼ 21 ms, t age ∼ 3220 yr, τ D ∼ 1.07×10 5 yr and B p (0) ∼ 2.97×10 13 G. Notice that the signs in Figure 3 are universally adopted in other figures in the subsequent sections of this work. The horizontal blue dotted line and the surrounding shaded region denote, respectively, the measured braking index of n = 3.15 and its possible range given by the uncertainty 0.03 of PSR J1640−4631. From Figure 3, it is obvious that n increases with t owing to the decay of the dipole magnetic field.

The evolution of the dipole magnetic field
By using Eq.(20) and the constrained parameters of t age , τ D and B p (0), we estimate the present value of the dipole magnetic field strength, B p (t age ) ∼ (1.77 − 4.10) × 10 13 G. This range is different from that estimated by the characteristic magnetic field at the polar of a pulsar B c ≃ 6.4 × 10 19 × (PṖ I 45 ) 1/2 = 2.87 × 10 13 G, because the latter assumes a constant dipole braking torque and a constant momentum of inertial I 45 = 1. Then we obtain a mean magnetic field decay rate of the pulsar, ∆B p /∆t = (B p (t age ) − B p (0))/t age ≈ −(2.41 − 3.51) × 10 7 G yr −1 , and plot B p versus t of PSR J1634−4631 in Figure 4(a).
The decay rateḂ p of PSR J1634−4631 is also an important issue. From Equation (20), we obtain the expression ofḂ p and t for the pulsar, The magnetic field decay rateḂ p of the pulsar declines with t, as shown in Figure 4(b). When t = t age , we get the present value of dB p /dt ∼ −(1.66 − 3.85) × 10 8 G yr −1 . The constraints on three parameters B p (0), t age and τ D ) of the pulsar also can be useful in accounting for the age difference between τ c and t age . Inserting Equation (23) If τ D → ∞, then Equation (28) is approximated as From Equations (28)-(29), we plot the diagrams of τ c vs. t for the pulsar in Figure 5. In Figure 5, the red dotted line denotes τ c = t age . It is easy to see that τ c increases with the age, but the fitted curves given by Equation (28) are always above the boundary of τ c = t age . This is because the characteristic age of a pulsar is always higher than its inferred age if n ≥ 3, which can be seen from Equation (5) in Section 2.
5.4. Diagrams of N − t, P − t,Ṗ − t andP − t. Considering a pulsar with rotational angular momentum L (L = IΩ(t)), the braking torque acting on the pulsar is given by where I is assumed to be constant in time. The decay of B p inevitably causes a decrease in the dipole braking torque N , which is described by If τ D → ∞, Equation (31) is rewritten as (32) Then we plot the diagrams of N − t for PSR J1634−4631, as shown in Figure 6(a). It is obvious that N decreases as t increases.
As a comparison, below we give the relations of P −t,Ṗ −t, andP − t in the case of the MDR model: (34) and (35) Based on Equations (22)-(24) and (31)-(33), we numerically simulate the relations of P − t,Ṗ − t andP − t for PSR J1634−4631. In Figure 6, the measured values of N , P ,Ṗ andP are shown with the red dots, and the errors in data point are smaller than the sizes of the symbols. In the terms of the dipole magnetic field decay model, as t increases, P increases rapidly at the earlier evolution stage and then increases slowly at the later evolution stage, whereasṖ andP decrease slowly at the earlier evolution stage and then decrease rapidly at the later evolution stage. Unlike PSR J1734−3333, which could evolve into a magnetar , PSR J1634−4631 will eventually come down to the death valley, due to the pulsar death effect.

ALTERNATIVE MODELS FOR THE DIPOLE
MAGNETIC FIELD DECAY It is worth emphasizing that our results and conclusions strongly depend on Equation (20) in which an exponential form of dipole magnetic field decay is assumed. In order to investigate how different magnetic field decay laws would affect our results, we assume that the dipole magnetic is deeply buried in the inner crust and extends into the core (see, e.g., Muslimov & Page 1995;Gourgouliatos & Cumming 2015). If the magnetic field decays via the Hall drift and ohmic decay with a nonlinear form: As a comparison, we give the relations of P −t,Ṗ −t andP − t in this nonlinear decay form. Inserting Equation (36) into Equation (22), we obtain the expression of the spin period, Utilizing the differential method, we get the first-order derivative of the spin periodṖ (t),  and the second-order derivative of the spin periodP (t), If we assume M ∼ 1.0-2.2 M and use the same APR3 EOS, then we obtain the effective magnetic field decay timescale τ D ∼ 0.91(8)×10 5 yr and the true age t age ∼ 3100−3200 yr, the initial dipole magnetic field B p (0) ∼ (1.72 − 3.83) × 10 13 G, and the dipole magnetic field decay rate dB p /dt ∼ −(1.76 − 3.94) × 10 8 ) G yr −1 . Alternatively, the magnetic field is assumed to decay with a power-law decay form: as done in Muslimov & Page (1995, where ε < 0 is the magnetic field index. In this case, Equations (37)-(39) are replaced by respectively. For the sake of simplicity, we assume that the range of t age obtained from the power-law decay form (Equation (40)) is very close to those obtained from the exponential and nonlinear decay forms (Equation (20) and Equation (36)), e.g., t age ∼ 2800 − 3200 yr. Utilizing the same NS mass range and fitting method, we get the magnetic field index ε = −0.034(1), the effective field decay timescale τ D ∼ (0.4 − 2.5) × 10 6 yr, the initial dipole magnetic field strength B p (0) ∼ (1.44 − 3.03) × 10 13 G, and the magnetic field decay rate dB p /dt ∼ −(1.78 − 4.42) × 10 8 ) G yr −1 . It is worth addressing that for ordinary radio pulsars the effect of field decay may not be very pronounced, and so it is difficult to uncover it. For example, Mukherjee & Kembhavi (1997) synthesized a population of pulsars using assumed theoretical distributions of age, initial magnetic field, period of rotation, position, luminosity, and dispersion measure and concluded that the timescale for the field decay must be greater than 160 Myr. Other studies (Bhattacharya et al. 1992;Hartman et al. 1997;Regimbau & de Freitas Pacheco 2001) also confirm these long timescales. We recall that the pulsars considered by these authors are older radio pulsars having the kinetic ages and/or characteristic ages τ ∼ 10 − 100 Myr, the dipole magnetic field strengthes B ∼ (2 − 5) × 10 12 G, which are different from young X−ray/γ−ray pulsar PSR J1640−4631 with τ c ∼ 3360 yr (the true age t age ∼ 2800 − 3200 yr). Furthermore, using the selection effect of pulsars' distance distribution and velocity distribution, Mukherjee & Kembhavi (1997) obtained the data for radio pulsars from the catalog of pulsars in the Princeton database (Taylor et al. 1995) and dropped binaries, X-ray pulsars and gamma-ray pulsars from the simulated pulsar samples. Thus, the magnetic field decay timescale obtained from Mukherjee & Kembhavi (1997) is larger than the one we give. Nevertheless, the magnetic field decay timescale of pulsars is an interesting subject worth long-term study.

COMPARISONS AND CONCLUSIONS
In this work, the high braking index of PSR J1640−4631 is well interpreted within a combination of MDR and dipole magnetic field decay. The related comparisons and conclusions are given as follows.

Comparisons with our previous work
In our previous work , we investigated the braking indices of eight magnetars. We attributed the larger braking indice (n > 3) of three magnetars to the decay of external braking torque, which might be caused by magnetic field decay, magnetospheric current decay, or the decay of magnetic inclination angle, and also calculated the dipolar field decay rates, which are compatible with the measured values of n for the following three magnetars: Compared with PSR J1640−4631, the three magnetars in Table 2 have higher magnetic field decay rates and larger braking indices. The main reasons are twofold:(i) magnetars are a kind of special pulsar powered by their magnetic energy rather than their rotational energy (Duncan & Thompson 1992;Gao et al. 2013;Liu 2017);and (ii) to account for the high-value braking indices of magnetars, we adopted the updated magnetothermal evolution model of Viganò et al. (2013) (see Gao et al. 2016 for details).
In addition, we show the long-term rotational evolution of PSR J1634−4631 in Figure 7. Five fitted lines in Figure 7 represent the evolution paths of the pulsar in the context of dipole magnetic field decay, and the arrows represent the spin- down evolution directions with different initial spin period in the next 100 kyr. PSR J1634−4631 would be placed in the top left region of the P −Ṗ diagram at its birth. As it spins down, the pulsar moves toward the bottom right of the P −Ṗ diagram in ∼ 10 5 − 10 6 yr and will eventually come down to the death valley, based on the dipole magnetic field decay model.

Comparsions with other models
Recently, to understand better the existence of both pulsar braking indices both larger than 3 and smaller than 3, different braking mechanisms have been proposed (e.g., Ekşi et al. 2016;Chen 2016;Clark et al. 2016;Coelho et al. 2016;de Araujo et al. 2016 a, b;Magalhaes et al. 2016;Tong & Kou 2017). For the purpose of illustrating the feasibility of this theoretical model, it is worth comparing our results with those of other models.
1. Initial spin period estimate. Recently, by numerically simulating, Gotthelf et al. (2014) predicted that PSR J1640−4631 has a short initial spin period P 0 ∼ 15 ms. Such a short P 0 requires a smaller braking index (n < 3) and a higher true age (t age > τ c ) (Gotthelf et al. 2014), which departures from the observed n and τ c of the pulsar obviously. In this paper, we firstly introduce a mean braking index n and a mean rotation energy conversion coefficient ζ and then estimate the initial spin period P 0 = 17 − 44 ms. The largest uncertainty of our model comes from ζ, because the true value of ζ is unable to be obtained. However, if P 0 ≪ P , the uncertainty of ζ will have little effect on the results of the three parameter of t age , τ D and B p (0), according to our model. In this work, the initial spin period (P 0 ∼ 17 − 44 ms) is obtained by combining pulsar's high-energy observations with the EOS.
2. True age estimate. Recently, using the leptonic PWN model 8 , Slane et al. (2010) assumed the age of G338.3−0.0 t SNR ∼10 kyr, and Gotthelf et al. (2014) gave an age t SNR ∼ 6.8 kyr. The former requires a constant braking index n = 3, while the latter requires a smaller braking index n ≈ 2.0 (Gotthelf et al.2014). However, their methods are universally based on supernova explosion mechanisms, due to the unknown explosion energy E and the ambient mass density ρ 0 . In this work, the true age estimate (t age ∼ 2900 − 3100 yr) is obtained by solving the pulsar spin-down evolution equation set.
3. Inclination angle estimate. Based on the MDR model with corotating plasma, Ekşi et al. (2016) found that the high braking index is consistent with two different inclination angles, α ∼ 18.5(3) • and 56 (4) • . Very recently, based on the vacuum gap model and low mass (M ∼ 0.1 M ) neutron star candidate, employing the vacuum gap model, Chen (2016) proposed that a lowmass neutron star (M ∼ 0.1M ) with a large inclination angle (close to the perpendicular case) could interpret the high braking index, as well as the radio-quiet nature of the pulsar. Note that in these two models a constant magnetic field has been uniformly assumed for PSR J1640−4631, and the discrepancy of α they obtained is obvious. As we know, it is the most meaningful to determine a pulsar's inclination angle from its timing and polarization observations. We expect that PSR J1640−4631 has a reliable measurement of α from the future timing and polarization observations. 7.3. No detected bursts in PSR J1640−4631 Since the higher-order magnetic moments decay rapidly with the distance to the stellar center, the influence of multipole magnetic moments of the star on the braking index can be ignored (if multipole magnetic fields exist in the pulsar). Interestingly, Chen (2016) supposed that there could be superhigh multipole magnetic fields inside PSR J1640−4631 (e.g., quadrupole magnetic field) resulting in a large stellar deformation and strong GW emission. If this is true, the burst behaviors should have been detected. However, up to date, no burst has been observed. If our model is correct, the undetected bursts in PSR J1640−4631 will be well explained as follows: by using the currently estimated values of B p andḂ p , we estimate the magnetic field energy decay rate, dE B /dt = d dt ( B 2 8π ) × 4 3 πR 3 = (1.2 − 2.8) × 10 32 erg s −1 , where B 2 p /8π is the magnetic field energy density. This magnetic field energy decay rate is about two orders of magnitude lower than the X-ray luminosity of PSR J1640−4631. Nevertheless, our prediction could be tested by the future highenergy observations of the pulsar.
In this work, we interpret the high braking index of PSR J1640−4631 with a combination of the dipole magnetic field decay and MDR models, ignoring the influences of all other possible braking torques. It is suggested that the decay of dipole magnetic field could occur universally in pulsars with n < 3. If so, why are the measured braking indices of the eight pulsars lower than 3 (see, e.g., Gao et al. 2016;Tong & Kou 2017)? Below we suggest two possible reasons. 1. The outflowing particle winds (the relativistic outflow mainly composed of the electron-positron pairs) luminosities in the eight rotationally powered pulsars with low braking indices may be so high that additional braking torques provided by winds can be comparable to the magneto-dipole braking torques in magnitude. Thus the stars' braking indices are less than 3. A PWN is dominated by a bright collimated feature, which is interpreted as a relativistic jet directed along the pulsar spin axis (e.g., Gaensler et al. 2003;Grondin et al. 2013). The PWNs of five low braking index pulsars have been detected. 9 Especially, the spin-down luminosity of the Crab pulsar ∼ 5 × 10 38 erg s −1 is converted into radiation with a remarkable efficiency, approaching ∼ 30% ( 2. Another possibility is that the eight pulsars are experiencing toroidal magnetic field (e.g., octupole field) decay via Hall drift and ohmic dissipation. For example, magnetar-like outbursts from three pulsar of PSR J1846−0258, PSR J1734−3333 and PSR J1119-6127 were reported (e.g., Gavriil et al. 2008;Göǧüş et al.2016). These outbursts could be caused by the decay of initial toroidal multipole magnetic fields. The multipole magnetic fields are merging through crustal tectonics to form a dipole magnetic field causing a low braking index.
3. For PSR J1640−4631, though it has an associated PWN with very low rotational energy transformation coefficient (see Section 3 of this work), the pulsar may be at an inactive wind epoch. Maybe this pulsar has experienced the stage of multipole fields merging, i.e., its dipole magnetic field stops increasing. At the current epoch, the emission and braking properties of the pulsar are dominated by the dipole magnetic field decaying, which causes a high mean braking index. In summary, attempts to explain the lower braking indices (n < 3) of pulsars within the dipole magnetic field decay model coupled with wind braking and/or toroidal field decaying will be considered in our future studies.

Conclusions
In this work, we interpreted the high braking index of PSR J1640−4631 with a combination of the MDR and dipole magnetic field decay. By introducing a mean rotation energy con-version coefficient ζ and combining the EOS with the highenergy and timing observations of the pulsar, we give the constrained values of P 0 , t age , τ D and B p (0), and numerically simulat the dipole magnetic field and spin-down evolutions of the pulsar. The high braking index of 3.15(3) of PSR J1640−4631 is attributed to its long-term dipole magnetic field decay at a low rate dB p /dt ∼ −(1.66−3.85)×10 8 G yr −1 . Considering the uncertainties and assumptions in this work, thus our theoretical model needs to be tested and modified by the future polarization, timing, and high-energy observations of PSR J1640−4631. Our results may apply to other pulsars with higher braking indices.