Evolution of Spin Period and Magnetic Field of the Crab Pulsar: Decay of the Braking Index by the Particle Wind Flow Torque

The evolutions of a neutron star's rotation and magnetic field (B-field) have remained unsolved puzzles for over half a century. We ascribe the rotational braking torques of pulsar to both components, the standard magnetic dipole radiation (MDR) and particle wind flow ( MDR + Wind, hereafter named MDRW), which we apply to the Crab pulsar (B0531 + 21), the only source with a known age and long-term continuous monitoring by radio telescope. Based on the above presumed simple spin-down torques, we obtain the exact analytic solution on the rotation evolution of the Crab pulsar, together with the related outcomes as described below: (1) unlike the constant characteristic B-field suggested by the MDR model, this value for the Crab pulsar increases by a hundred times in 50~kyr while its real B-field has no change; (2) the rotational braking index evolves from $\sim$3 to 1 in the long-term, however, it drops from 2.51 to 2.50 in $\sim$45 years at the present stage, while the particle flow contributes approximately 25% of the total rotational energy loss rate; (3) strikingly, the characteristic age has the maximum limit of $\sim$10 kyr, meaning that it is not always a good indicator of real age. Furthermore, we discussed the evolutionary path of the Crab pulsar from the MDR to the wind domination by comparing it with the possible wind braking candidate pulsar PSR J1734-3333.


Introduction
About 55 years has passed since the first pulsar was discovered in 1967 [1], identified as the beacon phenomena of rotating neutron stars (NSs) [2,3]. From then on, more than 3500 radio pulsars have been observed [4], including over 600 ones recently detected by fivehundred-meter aperture spherical radio telescope (FAST) [5][6][7][8], but the puzzle of how the rotation and magnetic field of pulsar evolves remains [9,10]. To answer these fundamental questions, the Crab pulsar (PSR B0531 + 21) is usually considered one of the best astronomical labs because it is the only pulsar with a known age. This famous pulsar was discovered in 1968 in the Crab nebula with a supernova remnant (SNR) [11][12][13][14] that was identified as a arXiv:2212.04674v1 [astro-ph.HE] 9 Dec 2022 product of a massive star explosion in 1054 AD with a clear historical record by ancient Chinese astronomers [15,16]. In addition, this young pulsar has been continuously monitored for half a century, yielding fruitful and accurate observational data [17,18].
Recently, the magnetic dipole radiation (MDR) model [2,19] was proposed and developed to account for the rotation slowdown of pulsar, in which the loss rate of the rotational kinetic energy (Ė) of pulsar is supposed to equate the emission power of a magnetic dipole in vacuum L d = K 1 Ω 4 . For example, a perpendicular rotator withĖ ≡ L d infers that −IΩΩ = K 1 Ω 4 , where the definitions of the conventional quantities are K 1 = 2B 2 R 6 /3c 3 ; I (10 45 g cm 2 )the moment of inertia; Ω (Ω)-the rotation angular frequency (its derivative); B-the NS surface magnetic field (B-field); R (10 6 cm)-the stelar radius; and c-the speed of light [15,16,20]. This standard MDR model has successfully predicted the B-field strength of the Crab pulsar and normal pulsars by the observed timing quantities as follows. In terms of spin period (P = 2π/Ω) and its derivativeṖ, the derived characteristic B-field is B = B ch ≡ [3c 3 I/(8π 2 R 6 )] 1/2 (PṖ) 1/2 10 12 (PṖ/10 −15 ) 1/2 G [15,16,20], which is very close to the direct measurement by the cyclotron absorption lines of X-rays for an accreting NS in the high mass X-ray binary (HMXB) [21,22]. Moreover, the evolution of the NS characteristic B-field has been widely investigated by previous studies [23][24][25][26]; and for a recent review, we refer to Ref. [27].
Without the direct measurement of B-field for a normal pulsar, the validity of the characteristic B-field as a replacement of the real B-field is often debated. To answer above doubts, astronomers can measure the braking index (n) of a pulsar that is defined by the spin timing parameters as n ≡ ΩΩ/Ω 2 , whereΩ is the second derivative of angular frequency, whereby the theoretical canonical constant value of n = 3 is obtained for the MDR model [28,29]. However, for the Crab pulsar, n = 2.515 was first reported in the 1970s [30,31], and its continuous monitoring followed the stable and accurate values of n = 2.51 in 1993 [17] and n = 2.50 in 2015 into the present, constituting an effort of 45 years [18]. Owing to the fact that the accurate measurement of braking, the index requires pulsars with a high-Ė or high-Ṗ. In general, only eight young radio pulsars have been measured, and to date, the stable values of braking indices have approximately ranged from 1 to 3 [18], which is deviated from the assumption of the basic MDR model. Apparently, the simple MDR model responsible for the pulsar spin-down torque should face a substantial modification [18,28,29,32,33]. One possibility is dedicated to the decoupling of the superfluidity vortex lines in the NS core that transfers the angular momentum into the NS crust; hence, through the variation of the moment of inertia I over time [34], the braking index can evolve to depart from the canonical value of 3 [35]; another possibility is that the B-field or magnetic angle between the spin and magnetic axis changes with time [36][37][38][39][40][41][42]; moreover, when introducing the plasma-filled and non-vacuum magnetosphere, the effects and some possible observed results of pulsar braking are also analyzed [43][44][45][46]. Besides the above, the particle wind flows responsible for the pulsar spin-down were also noticed [47][48][49], where the electric field accelerates the flow of electrons out of the magnetosphere and takes away the NS angular momentum [50,51], which can be also taken as the causes of a pulsar wind nebula (PWN) [52][53][54]. In addition, the on-off radio emission phenomena of intermittent pulsar PSR B1931 + 24 [55] and rotating radio transients (RRATs) [56] are interpreted as evidence for the switches of particle wind outflow [51]. The braking torques (labeled as T) of these particle flows can be described in the form of T ∝ Ω, or the loss rate of kinetic energy of pulsaṙ E ∼ L f = K 2 Ω 2 ∝ Ω 2 , with an undetermined parameter K 2 , which has been thoroughly discussed as a cause of the observed range of the pulsar braking index 1 < n < 3 [18,[57][58][59][60][61][62]. Furthermore, from the aspects of the Crab Nebula, the wind component should be a part of the pulsar braking [54,63]. According to the particle wind torque model for a pulsar spin-down, the parameter K 2 can be expressed as K 2 = πΦ 2 /(4c), where Φ is the magnetic flux of the particle flow. The above parameter K 2 was first studied by Michel [48], who extended the relativistic treatment from the solar-wind torque [64], an expression that we also apply in our work. Here, we consider the particle wind flow to modify the MDR model (hereafter referred as MDRW) by introducing the wind flow torque that can independently explain a braking index of n ∼ 1, which can gradually influence the constant braking index of n = 3, as predicted by the MDR itself.
In this article, the MDRW model is one of a thought experiment model that can account for the spin-down torque of the Crab pulsar and its analogues. Then, we acquire an analytical solution of the spin evolution for the Crab pulsar in Section 2, thereby deriving the evolution formula of the braking index, characteristic B-field, and characteristic age in Section 3. The tentative discussions and conclusions of the evolution evidence for the Crab pulsar are presented in Section 4.

Spin Evolution Model Furthermore, Results
In this section, we aim to give an exact solution when the additional wind torque is related to Ω 2 , and try to answer the possible evolution of the Crab pulsar. The reasons for which we chose particle wind outflow are not only due to the wind component being familiar with the astronomical issue and having been considered in the pulsar braking, as mentioned in the Introduction [48], but also because the Crab pulsar wind nebula has been observed, the luminosity of which is comparable to its rotational energy loss rateĖ [54]. Meanwhile, the Gamma ray luminosity by the FERMI-LAT observation is constrained to be approximately 13% ofĖ for the Crab pulsar [65]. Thus, we think that, without taking the particle flow into account, the simple MDR model is incomplete in understanding the spin-down evolution of the Crab pulsar.

Analytic Solution of the Pulsar Spin Evolution
For a pulsar, the loss rate of its kinetic rotational energy is assumed to equate the emission powers contributed by both the MDR (L d ) and particle wind flow (L f ), i.e.,Ė ≡ L d + L f , expressed below: where two undetermined parameters are defined by the MDR model [19] and particle flow model [48] for the spin-down of the pulsar, respectively, as mentioned above, or, equivalently, Equation (2) can be simplified as [18,62] −Ω = aΩ 3 + bΩ , with the undetermined dipolar parameter a = K 1 /I and flow parameter b = K 2 /I, where the condition b = 0 or K 2 = 0 corresponds to the conventional case of MDR. As expected, if the two presumed constants a and b were known, then the analytic solution of pulsar rotation might be achieved. By defining the fraction factors with respect toĖ contributed by the dipolar and flow components as which satisfy the condition of d + f = 1 that is the core assumption of the model. We proceed by submitting the spin derivative Equation (4) into the braking index Equation (1), and obtain a relation by which we have the following expressions These results demonstrate that the component fraction factors of pulsar emission are intimately related to the braking index. For the Crab pulsar, by employing the present observed value (denoted by subscript 'o') of the braking index n o = 2.50, we obtain d o = 3/4 and f o = 1/4, implying that the particle flow as a non-dipolar component has contributed to 25% of the total kinetic energy loss. Furthermore, by Equation (5), we have where τ c = Ω/Ω is defined as the characteristic age of pulsar. Necessarily, it is stressed that we do not use the definition of the characteristic age τ = Ω/2Ω as the magnetic dipole model does [16,66], because the braking index of our model evolves with time and departs from the canonical value 3 [28]. After arrangement, we obtained the two parameters of model a and b, respectively, which are calculated by the present spin period and its derivative of the Crab pulsar as, a = 2.7 × 10 −16 c.g.s and b = 3.2 × 10 −12 c.g.s, corresponding to K 1 = 2.7 × 10 29 c.g.s and K 2 = 3.2 × 10 33 c.g.s, respectively. Similarly, in terms of the fractional ratio ( ) of the particle flow relative to the magnetic dipole (hereafter referred as the flow-dipole ratio) defined by Ref. [18], ≡ L f /L d , we obtain = (Ω m /Ω) 2 = (P/P m ) 2 , where Therefore, combining Equations (4) and (11), the differential equation for the spin-down torque evolution can be transformed into the following form d /dt = 2b(1 + ), (12) and the exact analytic solution of above differential equation can be solved by the integral, then we have where i is an integral constant, or the initial value of that can be determined by the present spin parameters of the Crab pulsar at t o = 960 yr (see Table 1, the data from [18] and ATNF Pulsar Catalogue [4]), as i = (t = 0) = 0.1. Thus, the analytic solution of the Crab pulsar's spin evolution is settled down as below, which is plotted in Figure 1 and shown in Table 1, together with the other related parameters at different ages. To test the validity of this analytic solution Equation (13), we process it by the  Table 1.
Taylor expansion, then acquire the two special solutions, respectively, for the pure dipolar case b → 0 with Ω −2 = Ω −2 i + 2at and the pure non-dipolar case a → 0 with where Ω i denotes the initial angular frequency. For more details about the derivation of our model, please see Appendix D below.

Magnetic Field and Braking Index
The pulsar observable quantities such as the characteristic B-field, rotational braking index, and characteristic age are all determined by the spin periods and their derivatives. Therefore, their evolutions can be derived in terms of the spin period solution, which are described in the following.

Growth of Characteristic Magnetic Field
When the particle flow recedes to null ( = 0, f = 0 and d = 1), the MDRW model returns to the conventional dipolar case. If the particle flow switches on, then the loss rate of kinetic energyĖ = L d (1 + ) can derive the characteristic B-field as B ch ∝ √ PṖ, but the modification factor of (1 + ) to the dipolar case is introduced. Thus, the relation between the characteristic B-field (B ch ) and the real B-field strength (B) is given by which infers that B ch is not a constant and increases with the spin evolution. In other words, the real B-field of the Crab pulsar remains a constant and should be calculated by both the characteristic B-field and the flow-dipole ratio factor , expressed in terms of the pulsar parameters at the present time As shown in the diagrams of P −Ṗ (Figure 2), the B ch of the Crab pulsar increases by two orders of magnitudes to ∼ 10 14 G when the real time goes to ∼ 50 kyr. Thus, if our model is correct, then this result may indicate that the Crab pulsar may become a high characteristic B-field pulsar in the future, which could answer the previously proposed conjectures [32,67] that some pulsars may move to the high B-field range with evolution. In detail, the P −Ṗ evolutionary track of the Crab pulsar shows that it moves from its birth place with the initial spin period of P i = 18.3 ms, derivative ofṖ i = 6.4 × 10 −13 s s −1 , and the initial characteristic B-field of B chi = 3.5 × 10 12 G. Then, the path passes through the locations of the Vela pulsar (PSR B0833-45, n = 1.4, P = 89 ms, B ch = 3.38 × 10 12 G) [18] and the high B-field radio pulsar PSR J1734-3333 (n = 0.9 ± 0.2, P = 1.17 s, B ch = 5.2 × 10 13 G), arriving at those of the "magnetar" population, if there were no other factors influencing the pulsar spin down. It is noted that the real B-field of the Crab pulsar is still maintaining its same initial value of B = 3.3 × 10 12 G, and the characteristic field might appear much higher than the actual one if interpreted with the consideration of a smaller braking index. Here, it is clarified that we do not suspect the potential super-strong B-fields of magnetars, since the emission properties of both the Crab pulsar and most magnetars are very different; therefore, our simple but thoughtful model could not be automatically applied to those special sources of soft gamma-ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs) [68][69][70][71][72] with extremely intense high-energy outbursts (more explanations are noted in Appendix E). In Figure 2, the tendency of the B ch curve evolution can be understood by its evolutionary equation whereby the slop parameter k of the B ch curve is almost null at the early age, horizontally going along the constant B-field line (n = 3), and slop k gradually approaches unity while the B ch curve moves into the "magnetar" population (n = 1). Equivalently, in the P −Ṗ diagram, the slop of the B ch curve moves from k = −1 to 1, which is similar to the predicted route from the Crab pulsar via PSR J1734-3333 to the high B-field population [67], corresponding to the cases of the braking index of n = 3 and n = 1, respectively, [41]. For a further illustration, the evolution of the characteristic B-field B ch with the low values of the B-field and flow parameter b is also plotted, and its evolutionary path covers the vast population of normal pulsars in the P −Ṗ diagram.  [22,50]. The light-dashed lines stand for the different characteristic age and B-field, respectively. The data were taken from the ATNF Pulsar Catalogue [4]. The red stars represent particular pulsars, including the Crab pulsar, the Vela pulsar, and PSR J1734-3333. Below the symbol "Recycled pulsars" are those experienced the binary accretions [22,73,74].

Decay of Braking Index
Currently, approximately eight radio pulsars are measured by the stable braking index with long-term observation [18,75], and most of them are randomly and quite evenly spread in the expected range between 1 and 3. These phenomena are consistent with the predictions of the MDRW model [48]. Therefore, we can obtain the evolutionary equations of the braking index by solving Equation (10), that is With evolution, the flow-dipole ratio ranges from 0 to ∞, which is corresponding to the flow-total fraction factor from 0 to 1, implying an index ranging from 3 to 1 owing to the two extreme cases of dipolar or non-dipolar domination, which are plotted in Figure 3, where we notice the decay of the braking index with time or spin-down. Notably, the braking index of the Crab pulsar decreases from 2.51 to 2.50 over 45 years until the present stage, which is consistent with the observational results [18]. Meanwhile, some other evidence may also be consistent with the above relation, such as the specific pulsars mentioned in Figure 2, such as the Vela pulsar (n = 1.4) and PSR J1734-3333 (n = 0.9 ± 0.2), which are seemingly exhibited as the later evolutionary phase of the Crab pulsar [18]. If we assume that our thought experiment MDRW model has merit, then it can easily explain values such as these. Moreover, for the index n = 2, corresponding to = 1 as shown in Table 1, it represents a transitional point at which the dipolar and non-dipolar balance appears, with f=d=0.5 and P = P m , Ω = Ω m as listed in Table 1, implying that the two components account for the same proportion of total radiation.

Evolution of Characteristic Age
The characteristic age is usually used to estimate the approximate age of young, regular middle-aged, or old recycled pulsars, although it is known that this may not be accurate enough [66,[76][77][78]. It is interesting to examine the evolution of τ c in our MDRW model, as exhibited below: the evolution of which is shown in Figure 3. Strikingly, for the Crab pulsar, τ c has an upper limit of τ cmax = 1/b = 4τ o 10 kyr in our model, which corresponds to the flow-total fraction f = 1 or flow-dipole ratio → ∞, meaning that the characteristic age of a pulsar will not increase forever with evolution. In addition, the minimum value of the characteristic age τ cmin = 0.9 kyr is obtained by the initial condition. That is to say, τ c varies within a limited range of (0.9-10) kyr, therefore, it is a coincident event for the Crab pulsar that the characteristic age is close to the real age at the present stage. This result reminds us that the characteristic age cannot represent the real ages of pulsars in general, if the particle flows share a big portion in the total radiation. For a clear view, τ c values at various time stages are also shown in Table 1.

Measurement of Braking Index
For pulsar braking index measurement, stochastic timing variations by the NS spin instabilities exist which may be related to the spin-down torque, so it is not easy to acquire the precise braking index [79,80], and the elimination of noise is very important [81]. To date, only eight pulsars have reliable braking index values, because it is hard to obtain an accurate second derivative of the period. Generally speaking, there are three difficulties [82] which need long-term and continuous monitoring by radio telescopes: namely the Crab pulsar, which has a 45-year observation [18]; the pulsar timing irregularities including glitches, which may significantly affect the measurement accuracy, and the effect of the random walk related to the micro-glitches which may be supposed to be the main source for the impossible accurate braking index measurement; and some other aspects, which may have less potential influence but still need to be mentioned, including interstellar scatting, scintillation, etc. Wavelet analysis is a powerful time-frequency analysis method, like Fourier transform, which is particularly suitable for processing unstable signals such as complex noise or sudden changes in the data. Although this method is effective and has been applied in some fields of astronomy such as white dwarf [83] and pulsar [84,85], it has not yet been widely used in pulsar data processing from our knowledge. Perhaps the method of wavelet analysis will be helpful for the measurement of a braking index in the future [84]. Furthermore, due to the James Webb Space Telescope having been successfully launched and operated [86], more details of the Crab nebula will be discovered and analyzed. Since the magnetic field of NS has been estimated by its nebula [87,88] before the first pulsar was discovered, the surrounding environment is quite important for studying the characters of its central engine. Particularly for the pulsar wind, accurate observations and studies will enhance our knowledge of the torque of particle flow and additionally constrain and test our MDRW model.

Limitations and Assumptions of MDRW Model
Firstly, our MDRW model about pulsar spin-down has two parts, namely vacuum MDR component and wind flow. Vacuum around the pulsar should be a simple case (the modification may be due in the future, and more explanations are listed in Appendix E), and the wind flow is borrowed from a previous study [48]. Then, we assume constant a and b parameters, which mean that there is no change in the real B-field and magnetic inclination angle in our model. Therefore, variations in the above parameters will alter our conclusions, as shown in Figure 2 (various initial B-field and parameter b). Meanwhile, if the multipole B-field exists [62,89], the results of the paper should also be changed. After that, the Crab pulsar perhaps cannot represent all pulsars, because the braking index of PSR J1640-4631 is higher than 3 (n = 3.15 ± 0.03) [82], meaning that the other mechanism should be possible, and our model is probably not suitable for the apparent non-Crab-like group. While for most pulsars with accurate measurement [18], the braking index is less than 3, such as PSR B0540-69 (n = 2.14 ± 0.009) [90], the Vela pulsar PSR B0833-45 (n = 1.4 ± 0.2) [77], PSR J1119-6127 (n = 2.684 ± 0.002) [91], PSR B1509-58 (n = 2.839 ± 0.001) [90], PSR J1734-3333 (n = 0.9 ± 0.2) [67], PSR J1833-1034 (n = 1.8569 ± 0.001) [92], and PSR J1846-0258 (n = 2.65 ± 0.1) [90], our model may be suitable for these sources. Besides these, some other pulsars also have potential measurements [62] such as PSR J1208-6238 (n = 2.598) [93]. Thus, our model is limited to Crab-like pulsar with strong wind, while the population of pulsars can be diverse [78]. Further questions and problems about MDR and MDRW are addressed in Appendices A-E as listed below.

Theoretical Interpretation and Observation Evidence
Based on both the MDR and particle wind flow contributions to the NS spin-down torque, this study has obtained an exact analytic solution for a pulsar spin period evolution, thereby applied to the Crab pulsar to give rise to three primary outcomes, i.e., (1) the substantial enhancement of the characteristic B-field from ∼3.3 × 10 12 G to ∼10 14 G in 50 kyr, while the real B-field remains unchanged; (2) the decay of the braking index from 2.51 to 2.50 over 45 years to the present day and from 2.82 to 1 for a long-term evolution; and (3) the saturation of a characteristic age at 10 kyr with the continuation of time.
These tentative conclusions are consistent with the observations, parts of which were also noticed by astronomers [18,32]. Although the characteristic B-field increases approximately hundred times, the real B-field does not change at all, and hence, the Crab pulsar could possibly evolve to a high characteristic B-field pulsar such as PSR J1734-3333 [67,94]. Therefore, we do not think that the older Crab pulsar could evolve to a real magnetar with a super-strong true B-field [68], and our model is unsuitable for the magnetars due to their violent high energy emission outbursts (see Appendix E). Now, a question arises: do all pulsars evolve to high characteristic B-field pulsars like the Crab pulsar potentially does? Our answer is "No". If we consider the lower value of particle flow contribution (e.g., the low value of K 2 or b) than that of the Crab pulsar (see Curve-5 in Figure 2), we find that the P −Ṗ evolutionary path goes to the region of normal pulsars (∼10 12−13 G). In other words, most radio pulsars seem to follow almost constant B-field routes over millions of years. Thus, the different particle flow contributions may account for the distribution of pulsars in the P −Ṗ diagram, which reminds us to relate them to the two types of pulsars with and without the SNRs [78].
As known from the MDR model, the characteristic age is bigger than the real age, from which many debates and discussions have arisen [66,76]. Several observations have shown that the characteristic age of pulsar is far from the expansion age of its supernova remnant or proper motion age [75,94,95], which are usually taken as the indicators to estimate their real ages. In addition, it was noticed that the increment in the characteristic age and real age may be not synchronized [32]. However, our model pointed out that the discrepancy between the two ages may be due to the significant contribution of the particle flow component. Hence, we concluded that the characteristic age of pulsar is not a good indicator of the real age if sufficient particle wind flow exists. The wind component is approximately 25% of the rotational kinetic energy loss rate for the Crab pulsar (Ė = 5 × 10 38 erg/s) inferred by the braking index n = 2.5 [18]. In addition, the observed luminosity of the Crab nebula is approximately 1.3 × 10 38 erg/s, which is also noticed to be 26% of the total luminosity of the Crab pulsar [54]. Thus, the values of ratio (f) between the wind flow and total energy loss rate are consistent based on the two observations, the pulsar timing and nebula luminosity, which may favor our model to some extent.
Furthermore, we predicted the imaginable future parameters of the Crab pulsar based on the MDRW, while the torque of wind dominates its braking, e.g., after 20 kyr, its braking index is about 1, the spin period will slow down to ∼0.6 s, and the characteristic B-field will reach ∼3.5 × 10 13 G. Then, another question arises, namely that of whether there a source to proof this evolution path? The answer is a strong "maybe", considering the assumptions of our model. Currently, only a few radio pulsars have the reliable braking index measurements as mentioned above, among which only one pulsar has n very close to 1, namely PSR J1734-3333, holding with n ∼ 1, P = 1.17 s, and B ch ∼ 5.2 × 10 13 G [67], which are close to the future features of the Crab pulsar. For the pulsar PSR J1734-3333, both the estimated SNR age (23 kyr) [35] and proper motion age (45-100 kyr) [94] are much older than its characteristic age (8.1 kyr), which might hint at the fact that it may be saturated like the Crab pulsar probably is. If we calculate the ratio parameter between the wind flow and MDR with the age of 23 kyr and initial magnetic field of 3.3 × 10 12 G by Equation (16), we obtain ∼ 248, implying that the torque of the wind flow is higher than that of its MDR for PSR J1734-3333. Meanwhile, the corresponding braking index of 1.008 is obtained by Equation (18), which satisfies its observed constraint of n = 0.9 ± 0.2. Moreover, the X-ray luminosity of this pulsar is weak [67], which may be due to the fact that the high energy particles do not emit in our line of sight, or PSR J1734-3333 has insufficient overall energy (since itsĖ ∼ 5.6 × 10 34 erg/s is four orders of magnitude lower than the Crab pulsar). Thus, PSR J1734-3333 is a possible candidate for its contribution of wind flow in braking torque as it is significantly higher than that of MDR.

Acknowledgments:
We are grateful to Richard Strom for carefully reading the manuscript. Meanwhile, we sincerely thank all referees for the meaningful comments and suggestions, which have significantly improved the quality of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Basic Information for MDR and Wind Component
As one of the simplest theories to describe pulsar emission, the MDR model was developed by Gunn and Ostriker [19,96] based on pioneer NS proposals [2,3] soon after the discovery of the first pulsar [1], which was taken as a "standard version" for the evolution and emission mechanism in the pulsar and NS text books [15,16,20,[97][98][99][100][101]. It has been a successful model in pulsar astronomy for over 50 years, and is usually applied on the pulsars powered by their rotational energy loss. Moreover, the characteristic magnetic field of pulsar is estimated by the MDR, which has been popularly used by astronomers and astrophysicists worldwide and in the data base of the ATNF Pulsar Catalogue [4]. In the MDR model, the pulsar is regarded as a rapidly rotating rigid body with a strong magnetic field located in a vacuum environment, so the corresponding physical picture of MDR can be regarded as a simple magnetic dipole form. From the mathematical expression, the momentum is proportional to the angular velocity Ω 4 , and can be further expressed as IΩΩ. Therefore, if we want to know how the pulsar evolves with time under MDR, equating the differential equation related to Ω andΩ can be obtained as made by the "standard MDR model".
However, as illustrated in this work, the simple MDR cannot explain some phenomena, such as the braking index lower than 3 (n = 3 is expected by MDR). Thus, one of modified methods is introduced by the contribution of particle flow during the pulsar spin-down [28,29]. After introducing the flow, the pulsar braking is not only caused by MDR but also by the angular momentum loss due to particle wind. Mathematically, the latter effect can be considered to be proportional to Ω 2 [18,58]. Therefore, when considering the evolution under the MDRW model, we can also solve the differential equation between Ω andΩ (Equation (2) in the manuscript), and fortunately, acquire an exact analytic solution.

Appendix B. Characteristic B-Field, Braking Index, and Characteristic Age
These are important parameters for pulsar evolution, and there are some explanations for these three parameters. To begin with, the characteristic B-field is calculated based on the MDR model, which is derived by the period and period derivative. However, this value may not be the real B-field of the pulsar, and only some pulsars in the high-mass X-ray binaries have their measured value inferred by the cyclotron absorption [102]. Next, it is known that the rotation for normal pulsars are slowing down, so the braking index accurately describes the spin-down of the pulsar [103]. Under the MDR model, the braking index is obtained to be 3, while from the aforementioned observations, we know that none of the measured n is 3, and most of them are in the range of 1-3. Then, due to the braking of the pulsar, the characteristic age represents the time for the pulsar to evolve to its present spin period with a consistent rate. Therefore, astronomers used to take the ratio between the period and period derivative to estimate the approximate age, but now this needs to change. For example, the Crab pulsar, the real age of which is ∼960 yr, and whose characteristic age is ∼2500 yr (under MDR, this value changes to 1260 yr). The real age and characteristic age are not consistent.

Appendix C. Parameter List for the Crab Pulsar and MDRW Model
We added a parameter table for a better understanding on MDRW model that we used as follows. The coefficients in the model have already been mentioned in Introduction section, such as the moment of inertia (∼10 45 gcm 2 ) and radius (∼10 km) of the NS we applied. These coefficients are widely used in pulsar astronomy as well as in the ATNF Catalogue and Pulsar textbook. Based on these, the mass of NS can be estimated by I ≈ 2/5MR 2 , so the mass is approximate ∼ 1.3M (M ∼ 2 × 10 33 g). However, this result is under the assumption of a uniform sphere. If considering an inhomogeneous sphere, a modification coefficient of less than 1 needs to be multiplied (e.g., 0.7 or 0.8), so the calculated mass is ∼ 1M . The above values are consistent in our derivation and analysis, and even though they slightly change, they will not affect our results.

Appendix D. Derivation of MDRW Model
If we assume that the energy loss rate (Ė) of pulsars is not only dipole radiation (L d ), but may also contain other radiations, such as particle flow wind radiation (L f ), then we obtaiṅ The above equation can be written as where K 1 = 2B 2 R 6 /(3c 3 ) and K 2 = ηB are parameters described in Sections 1 and 2; I is the rotational inertia of a pulsar; Ω andΩ are the angular velocity and its first derivative, respectively; B is the real magnetic field strength of a pulsar; R is the pulsar's radius; c is the speed of light; and η is the coefficient of the particle flow radiation. Let a = K 1 /I and b = K 2 /I, then the above equation can be simplified aṡ and take the derivative of bothΩ Two sides multiplied by IΩ 2 can be obtained The expression for the breaking index is The numerator and denominator are multiplied by IΩ and take the Equation (A5) in it We can simplify Equation (A1) to obtain where d = L d /Ė and f = L f /Ė and Equation (A7) becomes Therefore, solving Equations (A8) and (A9) simultaneously can be obtained Then, we determine the coefficients K 1 and K 2 by applying the present observed values denoted by "o" for the Crab pulsar (B0531+21) breaking index (n o = 2.50), the energy loss rate (Ė o = 4.5 × 10 38 erg s −1 ), and angular velocity (Ω o = 188.2 rad s −1 ) We take I = 10 45 g cm 2 , R = 10 6 cm, and c = 3 × 10 10 cm s −1 into calculation. Thus, K 1 = 2.7 × 10 29 c.g.s and K 2 = 3.2 × 10 33 c.g.s, and a = 2.7 × 10 −16 c.g.s and b = 3.2 × 10 −12 c.g.s. Moreover, we know that d and f represent the proportion of L d and L f inĖ, respectively. Thus, the ratio factor between f and d is defined as (based on the definition of , as can also be seen in [18] where Ω m = 108.6 rad s −1 and P m = 57.8 ms are the mean values of Ω and P that correspond to n = 2, respectively. From Equations (A1) and (A2), we have two components of energy loss rate (Ė). We can simplify and solve the above differential equation and obtain the relationship between Ω and t.
Here, we write X = Ω 2 and the above equation can be rewritten as Solve the above differential equation Calculate the integral on both sides Therefore, we obtain We take the Crab pulsar's initial angular velocity (Ω i = 343.8 rad s −1 ) and t i = 0 yr into the above equation. Then, constant C = ( b X i + a) is obtained by us, which can give the observed angular velocity value of the Crab pulsar at t = 960 years, implying that the exact solution is consistent.
Then, we can examine our calculation with the conditions of b = 0 and a = 0, respectively. If b = 0, then energy loss rate goes back to the MDR (Ė = L d ), and we obtain the solution below Meanwhile, On the other hand, if a = 0, then the energy loss rate is occupied by particle flow radiation (Ė = L f ). Thus, the solution of the equation can be obtained, that is Correspondingly, Equation (A17) can be written as Then, if we take Equation (A17) into Equation (A12), we can obtain Thus, Equation (A17) can be understood as a general solution of Ω 2 , and Equations (A18) and (A20) can be understood as special solutions. Meanwhile, considering P = 2π/Ω, we can obtain the relation between the spin period (P) and time (t). When √ PṖ is the characteristic magnetic field strength. We takeĖ = L d into L d /Ė = (n − 1)/2, where L d is written by the real B, andĖ is written by P andṖ, and accordingly B ch , so we have the relation between B and B ch Then, the evolution of B can be written as As Equation (A9) told us, breaking index (n) can be written as and we take d = L d /Ė into it Furthermore, taking Equation (A12) into the above equation, we can obtain the evolution of the braking index (n) with the angular velocity (Ω), that is The characteristic age is defined as τ c = −Ω/Ω. Thus, τ c can be rewritten as Then, we take L f /Ė = (3 − n)/2 into the above equation, and we can obtain By Equation (A27), we have Therefore, τ cmax = 1/b = 3.15 × 10 11 s ≈ 10 4 yr.

Appendix E. More Explanation on MDRW Model
Although the MDR model is simple, we rely on it to make a modification and study pulsar evolution by adding the wind flow component, and we would like to explain it in some detail. To begin with, the MDR has been widely used in pulsar astronomy since it was created and it can describe the basic phenomena of pulsars. Then, many parameters are calculated under MDR, such as the characteristic magnetic field and energy loss rate, which were recorded in the mature database of the ATNF Pulsar Catalogue [4], as mentioned in the above text. Meanwhile, it is a mainstream method to revise MDR by considering magnetosphere and wind, which has been discussed by many researchers [10]. For example, researchers have analyzed the wind on pulsar spin-down [62,111], but they did not present the exact solutions for the Crab pulsar evolutions. What we achieved is that the exact solutions for the Crab pulsar evolutions of spin, characteristic magnetic field, and braking index are obtained, by which we can quantitatively study the Crab pulsar compared with its well-measured observational data. Therefore, if MDR is totally discarded, this means that the parameters in the pulsar database will not be available for model calibration, and we keep the basic MDR and add the wind component in our MDRW model. Furthermore, the vacuum environment assumed in the MDR is another widely concerning issue [111][112][113][114][115][116], but we do not take this into account in the current work because the plasma propagation effect in the magnetosphere is quite intricate, which will let us lose the analytical expression of the MDRW. Additionally, we do not consider the magnetic quadrupolar radiation in the present work, since the quadrupolar power should be related to the pulsar spin frequency in 5 power [62], while for most pulsars, the braking index is lower than 3, which is not consistent with the observed value of the Crab pulsar-namely 2.5. However, the magnetic quadrupolar structure is expected to be complicated [89,[117][118][119][120][121], which is far beyond our scope. Thus, in the present state, the MDRW is a concise and focused model to explain the braking index and evolutions of the Crab pulsar.
In the end, to avoid potential confusion, the issue of magnetars from the Crab pulsar should be clarified, because the simple MDRW model indicates that the characteristic B-field of the Crab pulsar grow from ∼10 12 G to ∼10 14 G, while the true B-field of the Crab pulsar has no change at all, which possibly hints at the potential existence of the "impostor" magnetars. However, this association is ambiguous and unrealistic, the reasons for which are considered below. To begin with, the conceptions of magnetars are based on soft gamma-ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs) that often exhibit the intense high-energy outbursts, which are well-explained by the assumed super-strong true B-field magnetars [68][69][70][71][72]. However, the outburst phenomena of SGRs/AXPs are also considered as other interactions between the wind flow materials or fall-back disk and the magnetosphere [122][123][124]. Then, the emission properties of the Crab pulsar are significantly different from those of the magnetars, and the persistent X-ray luminosity (L x ) of the magnetar is often higher than its spin-down kinetic energy loss rate, inferring that both sources should have different origins and properties [71,72,125]. It is remarked that perhaps the true B-field of some magnetars may be overestimated because the electron X-ray absorption cyclotron line has not yet measured them [21,102]-which is now also estimated and inferred by the spin period P and its derivative based on the MDR model. Conservatively and possibly, MDRW may be applied to several high characteristic B-field radio pulsars such as PSR J1734-3333 [67]; meanwhile, it shows a low braking index close to 1, which is consistent with the prediction of the wind flow torque itself. Thus, if the MDRW model dominates its spin-down torque, we just guess and assume its true B-field to be approximately ∼10 12 G. Nevertheless, we cannot assure this conclusion, since there is no direct measurement of the true B-field for this source. Additionally, the Galactic magnetar SGR 1935 + 2154 generated fast radio bursts (FRBs) [126,127], supporting the hypothesis of the magnetar origin of FRBs [128,129]. However, the properties of FRBs are far from the pulses of normal radio pulsars [130], such as energy, polarization, and narrow band, indicating that some basic differences between both should exist. In other words, we stress that it is too early to apply the simple MDRW model to all types of NSs, because their origins and formation mechanisms may be different [78,123,[131][132][133].