Spin evolution of neutron stars

In this paper we review the basics of magneto-rotational properties of neutron stars focusing on spin-up/spin-down behavior at different evolutionary stages. The main goal is to provide equations for the spin frequency changes in various regimes (radio pulsar, propeller, accretor, etc.). Since presently spin behavior of neutron stars at all stages remains a subject of many uncertainties, we review different suggestions made over the years in the literature.


Introduction
In many respects, observational appearances of neutron stars (NSs) are determined by their spin period and magnetic field.Thus, on the one hand, an understanding of magneto-rotational evolution is necessary to explain observed phenomena and to predict properties of yet undiscovered types of sources.On the other hand, it allows us to derive parameters of NSs (e.g., magnetic field) from observations and understand the formation and evolution of various sources.
Magneto-rotational evolution is partly driven by the internal physics of NSs (e.g., field decay), and partly by the interaction of an NS with the surrounding medium.Typically, NS activity (particles wind, electromagnetic radiation, magnetic field) works against a tendency of external matter to penetrate down to the NS surface either due to gravity or just ram pressure.This counterbalance of external and internal forces determines the stage of evolution of an NS.Here we use the following classification (see e.g., [1]): ejector, propeller, accretor, and georotator.Transitions from stage to stage often can be specified in terms of the equality of some specific radii, see Fig. 1.Let us define the most important of them.
The equality of kinetic energy and the absolute value of potential energy of the matter surrounding a compact object defines the gravitational capture radius (aka Bondi radius): Here G is the Newton constant, M is the compact object mass, and v is velocity relative to the medium accounting for the sound velocity c s : Magnetic field lines corotate with the NS.Thus, the linear velocity of a field line grows (in the equatorial plane) as ΩR, where R is the distance from the star's center and Ω = 2π/P is spin frequency.As this velocity is limited by the speed of light, there is a critical distance called the light cylinder radius: where c is the velocity of light.Plasma frozen in the magnetosphere corotates with it with the maximum velocity ΩR.If this value is larger than the local Keplerian velocity √ GM/R then a centrifugal barrier prevents accretion down to the NS surface.Equality of the linear and Keplerian velocity defines the corotation radius: Recently, Lyutikov [2] (see also [3], Sec.2.3) reconsidered the penetration of matter in the cases of spherical and disc accretion regarding centrifugal force.He derived a value of the centrifugal barrier, R br , that is smaller than R co .In particular, for the spherical accretion for an aligned dipole in the equatorial plane: R cb = (2GM/3Ω 2 ) 1/3 .
Magnetic pressure (energy density) of a dipolar magnetosphere behaves as P mag ∝ R −6 .Thus, it rapidly grows toward the compact object.The pressure of gravitationally captured free-falling matter grows as P ext ∝ R −5/2 under the standard assumption that the sound velocity is about the free-fall velocity, i.e.T ∝ r −1 , and density is ∝ r −3/2 from the continuity relation.Then at some distance R m the matter can be stopped by the magnetic field (for low fields it can happen that R m is smaller than the NS radius R NS , then the influence of the field can be neglected).The magnetospheric radius R m might be calculated differently in the case of disc and spherical accretion.Also, the radius is calculated differently depending on the relative value of R m and R G .For spherical accretion and R m < R G we have: This critical distance is also called the Alfven radius, R A .Here µ = BR 3 NS is a magnetic moment, where B is the equatorial surface dipolar field.As it is an approximate estimate, the numerical coefficient in the denominator can be different for different assumptions.The accretion rate can be estimated as: where ρ = m p n is the density of the accreted medium at R G (m p is proton mass).The coefficient ξ acc depends on details of accretion (for example, eq. 4 is obtained for ξ acc = 4π), and below we discuss it in every particular case.Note, that Ṁ can appear in equations even if accretion is not possible at that stage; then this value just characterizes the properties of the external medium.Shvartsman radius R Sh is determined by the balance between the pressure of relativistic particles wind and external medium.For R Sh > R G and standard magneto-dipole rate of losses we have: However, the coefficient in the numerator can differ for other models of rotation energy losses (see Sec. 3 below).R Sh must be larger than the light cylinder radius as it is assumed that the wind is generated in the vicinity of the light cylinder boundary.If R Sh < R G then the boundary between the wind and external plasma is rapidly shrinking as the wind pressure P w ∝ R −2 and the pressure of the gravitationally captured matter grows as P ext ∝ R −5/2 .Thus, at the ejector stage, we have R Sh > max (R G , R l ).
The ejector stage is over when the relativistic particle wind is not produced anymore and the Poynting flux becomes sufficiently low.Usually, it is assumed that this happens when external plasma can penetrate the light cylinder.This can happen due to spin-down, Ejector Propeller magnetic field decay, or modifications of external conditions.Inside the light cylinder, matter starts to interact with the magnetosphere.Unless some very exotic situation is realized, immediately after the transition from the ejector stage, R m is just slightly smaller than R l .This means that R m > R co as the Keplerian velocity cannot be larger than c.This corresponds to the propeller regime.Finally, when R m < R co and R m < R G matter can reach the surface of the NS.So, we have accretion.Some other more exotic situations can be realized.We discuss them in Sec. 6.
The basic picture of magneto-rotational evolution looks like this (see Fig. 2).An NS is born rapidly rotating (P ≲ 1 s) with a sufficient magnetic field (B ∼ 10 12 G) at the ejector stage (often as a radio pulsar).Then it slows down (and, maybe, the magnetic field decays) coming to the propeller stage at some t = t E .The transition happens at a critical period P E .Consequent spin-down can bring the NS to the stage of accretion, especially in the case of an isolated object, when another critical period, P A , is reached at t = t P This route can be followed by an isolated NS as well as by an NS in a binary system.At the accretion stage, the spin period can change due to changes in the angular momentum of the incoming matter.For example, in the case of an isolated NS, this can be due to the interstellar turbulence (see Sec. 5.2).Turbulence becomes important at the period P cr .At the stage of turbulent accretion, the spin period fluctuates around the value P turb .
This simple picture can be perplexed by an initial stage of fallback, non-trivial field evolution, variations of properties of the surrounding medium, etc. Mainly, in this review we focus on the spin evolution.Magnetic field evolution was recently reviewed, e.g. in [4].
In the following sections, we discuss the main evolutionary stages and transitions between them in different astrophysical scenarios.We describe the spin evolution of an NS in 'chronological' order, starting from the early stage of fallback of a part of supernova ejecta onto a newborn NS (Sec.2).Then we describe the ejector stage paying special attention to radio pulsars (sec.3).This stage is followed by the propeller (Sec.4) and accretor stages (Sec.5).In Sec. 6 we discuss some hypothetical regimes of interaction of an NS with surrounding medium.Our conclusions are presented in Sec. 7.
In the text we use convention that X x = X/10 x .

Fallback
In this section, we discuss the so-called fallback -an episode of very intense accretion of the ejecta onto a newborn compact object.The total amount of accreted matter can reach  ≲ 0.1M ⊙ and the accretion rate at early stages can be from a ∼few up to ∼ 10 4 M ⊙ yr −1 .For the first time this process was discussed more than half a century ago in [5,6].Modern approaches to modeling fallback are considerably influenced by the seminal paper by Roger Chevalier [7].
In [7] the author analyzed in detail different regimes of fallback accounting for heating due to radioactive decay and changes in opacity (in time and at different distances from the compact object).Asymmetries due to the kick velocity, rotation, and pulsar emission were also briefly discussed.The author used mainly a semianalytical approach so that all involved physical processes could be clearly visible and understood.Strong initial fallback is due to a reverse shock and interaction of the ejecta with the envelope.At the initial stages, the accretion rate dependence on time roughly follows the law Ṁ ∝ t −3/2 .Later, when the ballistic regime is established, the standard law of fallback Ṁ ∝ t −5/3 is realized.Duration of a quasi-spherical fallback depends on the properties of the progenitor and for red supergiants can reach ∼ 1 year.
Fallback properties at early stages significantly depend on the structure of the pre-SN.Thus, the total amount of the accreted matter also depends on the parameters of the progenitor.Emission properties are dependent on the trapping of photons by the envelope.Only at a low accretion rate, do they start to diffuse out.Before this moment, the effects of the central compact object are hidden.
Fallback can be important also to explain some properties of gamma-ray bursts and superluminous supernovae, see e.g.[8] and references therein.However, here we will not discuss these issues focusing on consequences related to the magneto-rotational evolution of NSs as the fallback can drastically modify the initial properties of these compact objects.For example, it can reduce the surface magnetic field at the early phases of evolution, can significantly change the spin period, and finally, can even change the evolutionary stage of a young NS (e.g., preventing the appearance of a radio pulsar).We discuss all these possibilities below.
Significant fallback (here the accretion rate Ṁ is more important than the total accreted mass ∆M, however, these two parameters are interlinked as ∆M = Ṁ(t)dt) can strongly influence the observational appearance of a young NS.In particular, the radio pulsar activity can be completely switched off.This possibility was suggested and analyzed in [9].If the pressure of the external medium falling onto the compact object is larger than both pressure of the relativistic wind and magnetic field pressure for all radii larger than R NS then the field is submerged at depth ∼ a few tens or hundred meters below the surface.Thus, the external field that determines the radio pulsar activity and many other manifestations of an NS is orders of magnitude lower than the initial value.For a range of parameters used in [9] it was found that the field diffuses out on a time scale ∼ 10 3 yrs.This value mostly depends on the depth of submergence which is determined by the amount of accreted matter.Later this process was studied in several papers in more realistic frameworks.
In [10] the authors modeled the field diffusion in more detail.In particular, it was demonstrated that due to the field submergence, the characteristic ages of young pulsars, defined as P/2 Ṗ, can be significantly different (larger) in comparison with their actual ages.The authors notice that due to compression the field in the crust is amplified by ∼ 2 orders of magnitude.And as the spatial scale for this field is relatively small it is a subject of relatively rapid Ohmic decay.With the parameters used in this paper, the field diffuses out in ∼ 10 3 − 10 4 yrs.Also, the authors notice that submergence can explain the absence of radio pulsars in some supernova remnants (SNRs).
A new interest in the fallback field submergence appeared in the 2010s.This was related to the necessity to explain the properties of central compact objects in SN remnants (CCOs).Wynn Ho was probably the first who used the field submergence to explain the observational appearance of CCOs [11].The author shows that accretion of ∆M ≲ 10 −4 M ⊙ is enough for significant submergence with the time of re-emergence ≲ 10 4 yrs.This is sufficient to explain the properties of the most studied CCOs.
A 2D model of submergence and re-emergence coupled with detailed calculations of field decay and thermal evolution was presented in [12].The authors obtained results qualitatively in correspondence with the study [11].Corresponding tracks in the P − Ṗ diagram can be found in [13].In 2D it was possible to analyze changes in the topology of the field during re-emergence.In more detail, this process was studied in [14].The authors presented a 2D model of field re-emergence after a fallback episode with ∆M = 10 −3 M ⊙ .It is shown that higher multipoles, which are strongly suppressed during fallback, re-emerge faster than the dipole.Thus, they can dominate at ages ∼ 10 4 − 10 5 yrs.This is important for the pulsar switch-on.
Details of field submergence were calculated numerically in several papers by Bernal et al.: [15], [16], [17].In [16] the authors considered the interaction of the fallback matter with a magnetic field loop anchored in the crust.They demonstrated with 2D and 3D modeling that for sufficiently strong accretion on a timescale ∼ 100 ms after the reverse shock starts to interact with the loop, the field (which is always below the accretion shock) is submerged.Convection in the accretion envelope is crucially important for the fate of the field.If the envelope remains convective then the field is not submerged completely.If the reverse shock is not developed at all the field is not modified significantly.
Field submergence can modify the evolutionary status of a newborn NS.In a 1D approach (i.e., for spherical symmetry) it was analyzed in [18] and [19].The authors aim to account for the back reaction of the NS (due to a relativistic particle wind and magnetic pressure) on the fallback flow.This allows them to explain different types of young neutron stars (radio pulsars, CCOs, magnetars, etc.) by different outputs of such interaction.In [18] the authors presented self-similar analytical considerations and in [19] -numerical modeling.
In [19] the model is mostly determined by two parameters: characteristic fallback time, t fb , and the fallback mass, M fb .The latter parameter can be calculated from the dependence of the fallback rate on time (notice, that it is not the rate of fallback on the NS surface as accretion can be prevented by relativistic wind or the magnetosphere): The fallback time was set to t fb,1 = (t fb /10 s) = 1.Among other parameters, it determines the initial encounter radius R enc at which the relativistic wind from the NS collides with the fallback matter: The luminosity of the wind, L sd , is calculated either with the usual magneto-dipole formula, or it can be enhanced due to accretion of some fraction of the fallback matter through a disc as it was suggested by [20].
The authors determine three possible regimes: significant field submergence, disturbed magnetosphere, and matter expulsion.If (R enc L sd )/(GM Ṁfb,ini ) < 1 then the fallback matter reaches the surface and submerges the magnetic field or significantly disturbs the magnetosphere (in the latter case the authors speculate about a magnetar formation due to field enhancement).
Fallback can result in the formation of a disc around the NS if the surrounding matter has enough angular momentum.These can be viscous accretion discs or passive discs like the one observed around the anomalous X-ray pulsar 4U 0142+61 [21].In this century, such discs have been extensively studied by many authors, starting with the paper by Chatterjee et al. [22], especially as an alternative to explain the properties of anomalous X-ray pulsars.See, e.g.[23,24] and references therein.Here we do not pretend to present an extensive review of fallback disc studies.Instead, we focus on a few papers mostly relevant to the subject of the present paper.
Ideally, fallback properties might be obtained from realistic calculations of SN explosions.At the moment, numerical models of SN are far from being complete.Still, new important results provide some optimism.In particular, recent calculations presented in [25] provide a realistic setup for analysis of the fallback influence on magneto-rotational evolution of NSs.
In this model, the newborn NS is not fixed in the center of the grid.The movement of the compact object allows accounting for new effects.In particular, the authors demonstrate that an NS can be spun up by the fallback and the obtained spin is correlated with the kick velocity.Vorticity in the fallback matter is not related to the spin of the progenitor.Instead, it is generated by non-radial hydrodynamic instabilities.Late fallback is associated with the highest angular momentum.Accretion of just ∼ 10 −2 M ⊙ can result in the spin period ∼ 0.05 s which is even slightly shorter than a usual value for young NSs [26].Typically in this scenario, the angle between the spin axis and velocity vector is ≲ 30 • .Alignment is more pronounced for large kicks as in this case, the NS can move farther from the center of the explosion, so accretion is less symmetric.For low kicks, an NS accretes from all directions and this results in averaging of the accreted angular momentum.According to the proposed model, fallback discs might be widely spread among young NSs with not very low spatial velocities.
Joint effects of fallback disc formation and field emergence on the appearance of radio pulsars are considered in [27].The authors perform population synthesis modeling accounting for spin up/down due to the NS interaction with the disc.In their model, the magnetic field is not submerged by the initial strong fallback.Instead, the authors use the equation identical to the one used for field behavior in an accreting binary system: Here B 0 is the field before accretion starts.In the population synthesis it is assumed to have a log-normal distribution with ⟨log B 0 /G⟩ = 12.35 and σ logB 0 = 0.4.And ∆M is the total amount of the accreted matter.After accretion stops, the field re-emerges with the rate: An NS can pass stages of accretion, propeller, and finally ejector.At the ejector stage, the NS can appear as a radio pulsar.The authors of [27], in the first place, intended to explain statistics of the braking indices distribution.By construction, in their model braking indices might be ≤ 3. Thus, the results are more or less compatible with the data on pulsars in SNRs.Still, it is very important that the authors used in their population model complicated early spin evolution due to NS-disc interaction and sub/re-emergence of the magnetic field.
Discovery of several long-period radio sources attracted new interest to the model of fallback discs as the interaction between an NS and the disc can result in a drastic spin-down of the compact object.The first new source with a peculiar period was the 76-second pulsar PSR J0901-4046 [28].In addition, two radio sources with periods ∼ 1000 s were found: GLEAM-X J1627-52 [29] and GPM J1839-10 [30].As their properties point towards magnetar-scale magnetic field, all three objects might be young.So, it is impossible to explain their spin periods simply by a magneto-dipole-like spin-down.Ronchi et al. [31] successfully explain long spin periods of freshly discovered sources in the framework of a fallback disc around a young strongly magnetized NS.
As specific angular momentum of the fallback matter is ∼ 10 16 − 10 17 cm 2 s −1 , see [25], the circularization radius where the accreted matter settles at a circular orbit forming a disc is R c ∼ 10 6 − 10 8 cm.The viscous time scale in the disc is t ν ≈ 2.1 × 10 3 T −1 c,6 R 1/2 c,8 s.Here T c,6 is the disc's central temperature in units 10 6 K and R c,8 = R c /10 8 cm.The mass flow rate in the outer parts of the disc is: The coefficient α depends on the opacity and in [31] the authors assume it to be equal to 1.2.
The accretion rate at the inner edge of the disc is limited by the Eddington rate.The authors do not account for the field submergence, but allow for the decay as they are interested in times scales up to ≳ 10 5 yrs.The field decay is calculated according to the analytical expression from [32]: Here τ Ohm is the Ohmic dissipation time scale which is assumed to be 4.4 × 10 6 yrs.And τ Hall,0 = 6.4 × 10 4 (10 14 G/B 0 ) is the initial Hall time scale which depends on the initial magnetic field B 0 .With such a setup and realistic fallback rates, it is possible to obtain spin periods ∼ 100 s for B 0 ∼ 10 13 G within 10 5 yrs and ∼ 1000 s for B 0 ∼ 10 14 − 10 15 G.

Ejector stage and radio pulsar activity
With more than 3000 known rotation-powered pulsars 1 , ejectors are the most common type of Galactic NSs observed.The rotational evolution of these sources is not affected by the surrounding interstellar medium but is driven by the complex electrodynamical processes within their magnetospheres.
For now, it is fairly clear, that most of the rotational energy of the ejectors is carried away by the wind of energetic charged particles (electrons and positrons) accelerated by strong electric fields within the magnetosphere of the star.At the same time, the physics of this process is still not well understood (there are detailed reviews on this topic [34,35]).Typical observed ejectors have spin periods in the range P ∼ 10 −3 − 10 s while their surface magnetic fields are of the order ∼ 10 9 − 10 13 G.This is consistent with observed spin-down rates of ejectors Ṗ ∼ 10 −20 − 10 −9 s/s (e.g.[36]).Therefore, their spin-down luminocities L sd = 4π 2 I ṖP −3 are of order 10 30 − 10 38 erg/s (see Fig. 3).
In classical theory, single neutron stars are viewed as rotating, highly magnetized, conducting bodies with a primarily dipolar magnetic field.Indeed, higher-order multipoles decrease very rapidly with the distance from the star.However, multipolar magnetic fields have been also considered in the literature (see e.g.[39][40][41][42]).
The magnetosphere of the star is created by its magnetic field, electric fields, and electric currents.If there is an electric field E ⊥ perpendicular to the local magnetic field B, a non-zero Poynting flux S ∝ (E × B) exists, which carries away the rotational energy of the star.This is true for both vacuum and plasma-filled rotating magnetospheres, which are two fundamental models for the physics of the ejector.
In the vacuum approach, the spin-down of a neutron star is due to the surface electric currents induced by the corresponding electric fields.So the deceleration of the star is driven by the Ampère's force.On the other hand, the environment of real neutron stars is not vacuum, but rather filled with plasma [43,44].In general, the magnetosphere of a neutron star consists of two parts, containing closed and open magnetic field lines respectively.The plasma can freely escape through the region of open lines, giving rise to the so-called 'pulsar wind'.Its back-reaction will also be responsible for the loss of the rotational energy [45].Notice, that magnetospheric currents define local electric and magnetic fields in a self-consistent manner.
So, finally, the spin-down torque acting on the neutron star can be described as the sum of two components: the vacuum and the wind.Hereafter we assume that the spin-down torque is positive and that the neutron star is spherical with a moment of inertia I.This means that for the ejector stage: For the vacuum component, the surface quadrupolar electric field E ∼ (ΩR/c)B is induced by the rotation, and so the corresponding spin-down torque is where µ ⊥ = µ sin χ is the orthogonal component of the magnetic moment, χ is the angle between the magnetic and spin axes, while f ⊥ is a structure factor.The latter is of the order of unity to a first approximation, although its details depend on certain assumptions.The vacuum losses are inefficient for a perfectly aligned rotator: K sd,vac = 0 when χ = 0.This is similar to the behavior of the emission from a rotating magnetic dipole.Furthermore, if f ⊥ = 2/3, then the equation ( 14) coincides with the classical one for the electromagnetic loss of a point-sized dipole.However, it is important to point out that there is no real magneto-dipolar emission from a neutron star (see the discussion in Section 2.3 of [35]).So the above equation is nothing more than a law for the spin-down of a magnetized, conducting, rotating sphere in a vacuum.As for the wind component of the ejector's spin down, its power L sd,w is based on the electric potential drop ∆φ near the star's surface, so that where r p = R NS √ ΩR NS /c is the star's polar cap radius, ρ GJ ≈ ωB ∥ /c is the Goldreich-Julian charge density [46], while κ ∼ 10 3 is the coefficient describing the primary particle density [47].Taking the maximum potential drop for a rotating dipole ∆Φ = µ ∥ /R 2 l [44] as Bottom panel: Sources with negative Ṗ are shown in the same coordinates.Most of these pulsars are in globular clusters and their "negative spin-down" is due to the Shklovskii effect [37].The pulsar graveyard border is plotted according to [38].Upper panel: Vertical lines show the spin period values for pulsars for which Ṗ has not yet been estimated.Grey lines are for all such objects, green ones mark the values for radio-quiet sources, red ones are for those exhibiting high-energy emission, while blue lines are for pulsars in binary systems.
an essential scale factor, the wind braking torque can be formally expressed in a manner similar to (14) as where µ ∥ = µ cos χ and f ∥ = 2κ∆φ/∆Φ is another structure factor that depends ultimately on the microphysics of the acceleration gap.The wind torque calculated as shown above disappears for a perfectly orthogonal rotator (χ = 90 • ) but is present for a perfectly aligned one.
At the turn of the millennium Harding et al. [48] and then Xu & Qiao [49] suggested that in fact, both spin-down mechanisms have a right to exist and that the total spin-down is driven by their combination.However, as the pulsar spins down, the potential drop may not sustain the need to accelerate particles [44].Thus, the maximum period P death exists where wind losses are efficient.This is the so-called 'death line' (or 'death period') of radiopulsars, beyond which these objects cannot be observed as radio-loud sources (since plasma acceleration and secondary particle production are thought to be related to pulsar radio emission).After crossing the death line, the wind component become negligible, so the pulsar's spin-down is expected to follow the vacuum approximation [50].
In this way, a complete torque that is responsible for the loss of rotational energy of a single ejector can be written as follows: The values of the these factors and their dependence on the parameters of the star are of particular interest for an accurate understanding of the rotation history of ejectors.Both quantities depend on the details of the physics of neutron star magnetospheres -the topology of the field and currents, the interaction of the currents with the stellar surface, and the physical structure and composition of the surface itself.And, in principle, there are two approaches to unravelling these factors theoretically: analytical and numerical ones.The purely analytical approach works well when considering a neutron star with a vacuum magnetosphere (and hence in determining f ⊥ , see the section 3.2), whereas the complex structure of the filled magnetosphere (taking into account both f ⊥ and f ∥ ) can only be revealed numerically as described in the section 3.3.But first, a brief review of the pulsar 'death condition' must be given.

Pulsar 'death line'
There has been much discussion regarding the 'death period' P death since the discovery of active radio pulsars.It is widely believed that pulsar radio emission is generated by the secondary plasma near the polar regions of neutron stars [43,44].The 'death' of pulsars is therefore closely linked to the possibility of the formation of electron-positron pairs.As shown by Ruderman and Sutherland [44], the critical condition for this process is the equality of the height of the vacuum gap (where ρ c is the radius of curvature of the local magnetic fields) to the polar cap radius r p .Corresponding 'death' condition can be written as where a ≈ 0.6 s, b = 4/11 and Ṗ−15 = Ṗ/(10 −15 s/s).This equation represents a monotonous curve (in particular, a straight line) on the P − Ṗ diagram for neutron stars when its plotted in logarithmic scales.This is why the pulsar 'death condition' is often referred to as the 'death line'.However, equation ( 18) only partially explains the properties of the observed population of radio pulsars since there are many sources with P > P death if the 'death period' is calculated using eq.( 18).Then Chen and Ruderman [51] considered several variants of the magnetic field structure near the polar cap and obtained refined parameters for the death condition.In particular, they found (a, b) =(3.2 s, 4/9) and (11 s, 1/2) as possible solutions.Accounting for relativistic frame dragging affects the condition of the electron-positron pairs production due to additional induced electric fields [52].Thus, Zhang et al. [53] have considered two polar cap acceleration models and found P death ≈ 0.67 Ṗ2/5 −15 s for a realistic death line in the case of dipolar and P death ≈ 5.75 Ṗ1/2 ρ −1/2 c,6 s for multipolar field configurations, respectively.
In the influential paper by Faucher-Giguère and Kaspi [54] describing the population synthesis of isolated radiopulsars, the death line condition was used in the form which corresponds to P death ≈ 3.28 Ṗ1/3 −15 s.This relation was earlier established empirically by Rawley et al. [55] from observations of a newly discovered long-period pulsar.In spite of its purely empirical nature, it was successful in reproducing the observed population of active pulsars.However, many authors (like Gullón et al. [56] and Graber et al. [57]) have not assumed the existence of a death line at all, since pulsars become undetectable earlier due to small luminosities and beaming fractions.
At the same time, Beskin et al. recently presented a more sophisticated analysis of the processes near the polar cap of active pulsars [38,58] and found that the classical equation ( 18) with b = 4/11 is still a good approximation to the death line condition but either a = 8.27 or a = 4.6 s depending on the spin-down model.In fact, they have considered both a numerical MHD model based on Spitkovsky's work [59] and their own analytical model [60].They argued, that the analytical model (with a = 8.27 s) fits the observations much better since it allows many long-period pulsars to exist.However, the pulsars with the longest known periods (J0250+5854 with P ≈ 23.5 s, J0901-4046 with P ≈ 76 s) are still contradict this model.So the physics of pulsar 'death' is far from complete.
In addition, it is worth mentioning that the 'death line' is likely to depend not only on its magnetic field (or spin-down rate), but also on the magnetic angle (see e.g.[35]).This fact is often ignored in the analysis.A reasonable form for the 'death' condition can then be approximately expressed as where A ≈ 0.4 − 0.5 from the analysis of more than 150 pulsars with known χ [61] (see also [40]), B 12 = B/(10 12 G) -is the equatorial magnetic field of a pulsar and period P is in seconds.This means that nearly orthogonal pulsars shut down much earlier than aligned ones, which affects their statistics.

Spin-down of a star with vacuum magnetosphere
The braking torque resulting from the non-zero Poynting flux of a magnetized star rotating in a vacuum can be calculated analytically and in great detail.
As early as November 1967 (a few months before the official discovery of radiopulsars by Hewish and Bell [62]) Pacini published a note briefly discussing the energy budget of a rotating neutron star with a dipolar magnetic field [63].He suggested that such a star could lose its rotational energy with luminosity L sd ∼ µ 2 Ω/R 3 l .But more than a decade earlier, Deutsch had derived [64] the first vacuum solution for the electromagnetic field of a perfectly conducting, rigidly rotating spherical star. 2 In general, the braking torque 2 See also an interesting pedagogical note by Satherley and Gordon [65] regarding Deutsch's classic work.
vector K sd,vac acting on such a star consists of three orthogonal terms.If e Ω and e µ are unit vectors collinear with the spin and magnetic axes of the star, respectively, then Here K sd,0 = µ 2 /R 3 l , while η 1,2,3 (r NS ) are dimensionless functions, depending on the assumptions made about the structure of the fields inside and outside the star [66][67][68][69].These functions depend on the small parameter r NS = R NS /R l .
In terms of the spin-down equation ( 17) one then has f ⊥ ≡ − η 3 .As the third term of ( 21), the second term represents the braking torque arising due to radiation in the far (wave) zone.It is responsible for the evolution of the magnetic angle χ.In the vacuum solution η 2 > 0, therefore the magnetic angle decreases systematically ( χ < 0) on the timescale of the pulsar spin-down.Moreover, typically The Euler equations describing the evolution of a spherical star with moment of inertia I, are therefore as follows: for the evolution of the spin period and for the evolution of the magnetic angle.As for the first term in ( 21), proportional to η 1 (r NS ), it does not directly affect the observed spin evolution of the star.Since it is always perpendicular to both the rotational and magnetic axes, it causes a precession around the magnetic axis with angular velocity or, in terms of period, T 1 ≈ 6 × 10 4 PB −2 12 cos −1 χ years.Here we assumed I = 10 45 g cm 2 , R NS = 10 6 cm, and the spin period P is taken in seconds [68,70].The value of T 1 is so short relative to the spin-down time scale at the ejector stage because η 3 (r NS ) ∼ 1 to a first approximation, while η 1 ∝ 1/r NS .Indeed, the first term in ( 21) can be interpreted either as a result of the inertia of the magnetic dipole itself [71], or as a near-field radiation torque [69], or as the appearance of an internal field momentum [35].This torque is often referred to as the 'anomalous braking torque' (see e.g.[72]).Melatos [69] has given an accurate calculation of η 1,2,3 , assuming the internal magnetic field to be a point dipole field: and respectively.Thus it is clear that if only leading order terms are taken the account, then As a result, unlike the 'conventional torques' which are proportional to 1/R 3 l , the 'anomalous' one scales as ∝ 1/R 2 l , making it up to four orders of magnitude larger, since r NS ∼ 10 −4 − 10 −2 for most real pulsars.Extensive analysis made in [69] has shown that this torque could be of high importance for a non-spherical neutron star, causing significant quasi-periodic variations in the spin-down rate and magnetic angle of a star with a typical timescale T 1 .The first term (∝ r −1 NS ) of the equation ( 26) has also been considered by Beskin & Zheltoukhov in [72].They have shown, that if the angular momentum of the internal electric field is also taken into account, then the anomalous torque depends on the structure of the magnetic field inside the star.Thus, if there is a uniform magnetic field up to radius R in < R NS , while the dipole configuration is relevant for R > R in , then It is interesting, that R in = 0 leads to the divergence of the η 1 and, therefore the divergence of the anomalous braking torque.However, this situation seems to be unphysical and would never be realized.

Spin-down of a star with plasma-filled magnetosphere
Much work has been done in the last two decades to obtain self-consistent solutions of pulsar magnetospheres in both axisymmetric [73][74][75][76][77] and oblique [59,[78][79][80][81][82] configurations.Significant progress has therefore been made in this area.Most models have assumed the force-free approximation, but radiative magnetospheres have also been considered (e.g.[83]).It has been shown that a magnetized spherical rotator with a plasma-filled magnetosphere is indeed affected by a torque proportional to µ 2 /R 3 l .However, the corresponding dimensionless coefficient is slightly different from that found in the vacuum approximation.
In his seminal paper Spitkovsky obtained a general numerical solution for the spindown torque of an oblique rotator [59]: where k 0 ∼ k 1 ∼ 1.Later Philippov et al [81] refined the parameters of the spin-down law: That is, in terms of (17) we get f ⊥ = k 0 + k 1 and f ∥ = k 0 .Since for most real pulsars r NS ∼ 10 −4 − 10 −2 is small, a good approximation is to keep k 0 and k 1 constant.Further PIC modelling [82] also produced to a result that is consistent with the results of MHD simulations.Namely, it was found that k 0 = 1.0 ± 0.1 and k 1 = 1.1 ± 0.1.If the radiative magnetosphere is taken into account then k 0 = 1.36 − 1.42 and k 1 = 1.73 − 1.76 [83].
The solution by Philippov et al. also contains the component responsible for the magnetic angle evolution.It is with the second term of (21) (or with scalar equation ( 23)) assuming η 2 ≡ 1.So, numerical consideration of pulsar magnetospheres predicts magnetic alignment of radio pulsars on spin-down timescales.
This result is the most complete and accurate theoretical description to date of the spin evolution of a single neutron star at the ejection stage.If one does not care about the value of the magnetic angle and its evolution, then spin-down losses can be estimated quantitatively simply by assuming an isotropic distribution of magnetic angles such that ⟨sin 2 χ⟩ = 2/3.Thus, assuming (30), one finally obtains where ξ ≈ 1.93.For actual numerical analysis, ξ = 2 remains a good approximation.

Observational verification
In general, one has to accept that the theory of pulsar spin-down described above allows one to understand and predict the properties of the actual ejectors in much detail.On the other hand, the large number of neutron stars observed as active radio pulsars and the high precision of measurements of their timing parameters (P, Ṗ, P, etc.) suggest that this theory can indeed be well verified in observations.But the reality is much more complicated.Despite many indirect arguments, there is no direct observational evidence to support the theory.
In this subsection, without pretending to provide a complete review, we briefly describe three relevant issues that remain unresolved and deserve to be discussed.

Magnetic angle evolution
The evolution of the magnetic angle χ during the spin-down of an isolated neutron star is an important part of its evolution.In a more general approach, the Euler equations ( 22) and ( 23) can be rewritten in terms of the so-called symmetric and anti-symmetric components of the magnetospheric currents normalized to the Goldreich-Julian current: i s = j s /j GJ and i a = j a /j GJ respectively (see e.g.[35]).In this approach one gets [60]: and It is clear from these equations, that the magnetic angle tends to evolve in such a way that the loss of rotational energy decreases with time.This means that if a perfectly aligned rotator (χ = 0) will spin down less quickly relative to the perfectly orthogonal one (χ = 90 • ), then one should expect a magnetic alignment ( χ < 0).On the other hand, if the physics of pulsar spin-down suggests that an orthogonal rotator loses rotational energy less efficiently, then magnetic orthogonalisation should occur ( χ > 0).As shown in Section 3.2, spin down is suppressed in vacuum magnetospheres when χ = 0, and hence χ < 0 there.The same result has also been obtained in numerical simula- tions of plasma-filled magnetospheres.However, in the 1980s a theory was developed by Beskin, Gurevich and Istomin [60,[84][85][86] which assumes a non-vacuum magnetosphere of a neutron star, but predicts magnetic counter-alignment.The MHD-based theory presented in the paper by Philippov [81] and the analytical description by Beskin et al. (BGI theory) give clear predictions about the evolution of the magnetic angle from fundamental physics.To our knowledge, these are the only two such considerations published to date.And their predictions are not consistent with each other.On the other hand, there is still no observational evidence in favor of one of any of these approaches, although many analyses have been carried out [54,61,[87][88][89][90][91][92].So this part of the rotational evolution of an isolated NS is not entirely clear.

Pulsar timing irregularities and braking indices
The spin-down equation (17) assumes that an isolated ejector loses its rotational energy monotonically.However, the actual observed spin evolution of such objects is not like this.At least about 6 per cent of known radio pulsars exhibit sudden and rapid spin-up events, known as glitches [93][94][95].However, the spin-down rate dP/dt of almost every pulsar appears to be significantly contaminated by additional irregular, almost stochastic variations on short time scales of months to years, typically referred to as 'timing noise' [96][97][98][99][100][101].
The strength of this unmodelled component of the pulsar spin-down can be characterized numerically in various ways.The standard approach is to use either the second derivative of the pulsar spin period [102] or the variance of the timing residuals (e.g.[97]).These two measures are correlated.There is also a widely used dimensionless combination: -the so-called 'braking index'.It has a simple and clear physical meaning.In particular, if the loss of rotational energy follows the expression L sd = −KΩ n+1 with K = const, then the combination (34) is nothing other than n.
Since the spin-down law for a spherical star is where the parameters of a neutron star are assumed to be variable and characteristic age of a pulsar.Classical vacuum losses ( f ⊥ = 2/3, f ∥ = 0) with constant magnetic fields and moment of inertia give n br = 3 + 2 cot 2 χ ≲ 103 for real pulsars.Only the third term in the brackets is non-zero in this case due to χ < 0. The theory by Beskin et al. [84] gives a similar result with n br ≈ 1.93 + 1.5 tan 2 χ.On the other hand, Spitkovsky's model [59] predicts a narrower interval for this quantity under the same assumptions (see also [103]).Namely, But a big problem is that the estimated values of n br for hundreds of pulsars are surprisingly far from all these predictions.The values have been found to range from ∼ −10 6 to ∼ 10 6 , and are negative for about half of the objects (e.g.[70,104,105]).They are therefore unlikely to represent the secular spin evolution of radio pulsars.In addition, the braking indices of most pulsars do not seem to be stable from observation to observation 3 There are only about a dozen sources that are accepted as having meaningful and stable n br ∈ 0.03-3.15[107,108].
Although the physics of pulsar timing irregularities remains generally unclear, the proposed solutions to the 'anomalous braking index' problem can be qualitatively divided into two categories.On the one hand, there are theoretical models that assume the relatively slow variability of the NS parameters affecting equation (35).The authors of these ideas assume either the variability of the magnetic moment [13,105,109,110] or the magnetic angle [69,91,92], or the effective moment of inertia [111][112][113][114][115] of the star.On the other hand, some researchers suggest the existence of an additional either quasi-periodic or purely stochastic external component in the spin-down torque Ωext .The underlying physics has been proposed to be related either to magnetospheric perturbations [116][117][118][119], or the imprint of the anomalous vacuum torque discussed above [70,106,120].
Biryukov et al. [70] attempted to estimate Ωext from the statistics of the observed second derivatives of pulsar spin frequencies.Assuming that anomalous braking indices are due to long-term variations in pulsar spin-down (with typical timescales of hundreds and thousands of years), they found 0.5 < ε < 0.8, where ε = | Ωext / Ω| -is the relative strength of the effect.
On the other hand, it has been shown in many papers that the observable parameters of Galactic pulsars can be well reproduced within a population synthesis adopting (17) with various forms of f ⊥ and f ∥ [54,56,121] but neglecting the corrections due to ε ̸ = 0.Moreover, the model-independent analysis of pulsar kinetics in the P − Ṗ diagram even led to a reasonable value of their birth rate of ∼ 1 per century, simply assuming f ⊥ = f ∥ ≡ 1 and ε ≡ 0 [122,123].This ultimately means that the simple equation for the isolated ejector spin-down torque K sd = ξµ 2 Ω 3 /c 3 , where ξ ∼ 1 is a 'universal constant', is still meaningful within a numerical analysis.

Isolated neutron stars birthrate
Knowledge of the birth rate ṄINS of Galactic isolated neutron stars (INSs) is also important in the context of neutron star population analysis and, in particular, details of their spin down.In general, several types of objects contribute to the total value of ṄINS : (i) ordinary rotation-powered pulsars (RPPs), (ii) rotating radio transients (RRATs), and (iii) young radio-quiet neutron stars (XDINS or the "Magnificent Seven") 4 , (iv) central compact objects in supernova remnants (CCOs) and (v) magnetars.
The birthrate of each of them, however, is not simply proportional to the number of observed sources but must be revealed assuming models of their spin-down, activity, magnetic and thermal evolution.For instance, despite the overwhelming number of ordinary pulsars, a significant fraction of neutron stars appear to be born as magnetars [127,128].
However, the total INSs birthrate is not simply the sum of birthrates of their different types.Keane and Kramer have shown that being simply summarized they give the value ≳ 5 − 10 century −1 [122] which is in tension with a proposed rate of Galactic core-collapse supernova explosions as it has been constrained by different methods [122,[129][130][131]. Within these calculations, the birthrate of RPPs, as the most numerous class of observed sources, can be estimated with formally small uncertainties.And understanding of their spin-down evolution is crucial, since ṄRPP ∼ N RPP /t ej .
So the analysis of pulsar kinetic current by Vranešević et al. [132] initially gave ṄRPP ∼ 1 − 2 century −1 .Later Lorimer et al. [133] used the same approach and estimated ṄRPP to be 1.4 ± 0.2 century −1 .And finally, Vranešević & Merlose [123] perform an even more accurate consideration, taking into account the pulsar spin-down law in the more general form P n−2 Ṗ =const.They found ṄRPP = 0.8 ± 0.3 century −1 for n = 2.9 ± 0.1.All these results are more or less consistent with each other.
Another way of measuring the birth rate of RPPs -a population synthesis approachtends to give even higher values (probably due to selection effects being taken into account).For example, Faucher-Giguère and Kaspi [54] reported ṄRPP = 2.8 ± 0.5, Ridley and Lorimer found ṄRPP from ∼ 2 to 5.3 century −1 depending on the details of the assumed spin-down model [121].Although a spin-down law of the form P Ṗ ∝ µ 2 (k 0 + k 1 sin 2 χ) was included in the calculations in the latter paper, the exponential decay of a pulsar obliquity was also assumed.which is inconsistent with the MHD model of pulsar rotational evolution.Finally, Gullón et al. found ṄRPP = 2.5 ± 0.1 century −1 including pulsar magneto-thermal evolution in their simulations [56].Therefore, the pulsar population synthesis generally predicts values of ṄRPP within 2-3 century −1 -about one pulsar per century more than the kinetic current approach and which is close to the rate of core-collapse supernovae themselves.
The obvious solution to the paradox of the overestimated birthrates of Galactic neutron stars is to assume that their populations are connected.Their observable properties depend on their current spin period, age, magnetic field strength and surface temperature, while the transition between the populations is driven by the magneto-rotational evolution [134][135][136][137].The latter is therefore crucial for understanding the properties and fate of INS.They need to be considered as a connected population of sources, unified by the common origin, rather than a set of separate classes5

Propeller stage
The propeller stage is, probably, the least understood one.On one hand, from the observational point of view, it is not easy to study this phase of evolution.This is so, partly because an NS luminosity can be rather low at this stage.In addition, the stage can be relatively short in many cases.On the other hand, the theoretical analysis of propellers is complicated and results of modeling often are not conclusive.We start by presenting different theoretical approaches to model the propeller stage.And then, observational results are reviewed and discussed.

Theory
The idea that a rapidly rotating magnetosphere can prevent penetration of matter down to the surface of an NS was initially proposed by Shvartzman [138].However, as the paper by Shvartsman has been published in a not-so-well-known journal, this idea became popular only later, mainly due to the paper by Illarionov and Sunyaev [139] who also proposed the name for the evolutionary stage: propeller.The existence of the centrifugal barrier for accretion onto NSs was also mentioned by [140], and discussed in more detail by [141], who also provided an estimate for the spin-down torque: Later on, a more detailed study of the propeller stage was presented by [142], [143].
Realistic 3D modeling of the interaction between the magnetosphere and surrounding matter is a very challenging task.Not many attempts have been made.That's why, early studies were analytical and typically limited by some simplified but physically motivated approaches.
At the propeller stage, a rapidly rotating magnetosphere transfers rotational energy and angular momentum to the medium.Thus, the matter is heated and moves away carrying some angular momentum, so that the NS spins down.Depending on the relative efficiency of heating and cooling, the properties of the shell around the magnetosphere can be different.External pressure at the magnetospheric boundary determines its radius, R m , which can significantly deviate from the Alfven radius, R A , see eq. ( 4).Without a disc formation, typically R m > R A .
In the simplest approach, e.g.[144], it is assumed that the external matter is gravitationally captured and cools effectively.Thus, . This provides a large density at the magnetospheric boundary.The magnetospheric boundary is determined from So, R m ∝ B 4/13 P 4/13 ρ −2/13 ∞ v 6/13 .Then, it is assumed that the matter is expelled with the rotational velocity of the magnetosphere: ΩR m .Thus, it carries away a huge amount of angular momentum.Altogether, this results in a very effective spin-down: So, the spin period increases rapidly.
Both assumptions (effective cooling and large velocity of the expelled matter) in this approach were strongly criticized.At first, due to significant heating, the density at R m might be much lower than Then, the velocity ΩR m is too high for the outflow.Davidson and Ostriker [141] assumed that the angular momentum is carried away with the free-fall velocity: In the seminal paper [139] the authors define the magnetospheric boundary as the usual Alfven radius and it is assumed that rotational energy loss is (1/2) Ṁv ff (R m ) 2 .Then the spin-down torque is: Fabian [145] obtained the radius of the magnetospheric boundary by equating the magnetic pressure and stellar wind pressure.To calculate the spin-down torque he used energetic considerations.In this approach, the rotational energy losses 1  2 I(Ω 2 1 − Ω 2 2 ) are equal to the kinetic energy of escaping matter In this model, the spin-down rate is identical to the one by Illarionov, Sunyayev [139].It can be written in an alternative form: Note, that R m can be defined in different ways.What is important, is that K sd ∝ Ω −1 .Then the characteristic time scale of spin-down τ sd ≡ P/ Ṗ ∝ (P −2 1 − P −2 2 ).This means the time scale is defined by the initial period.
In their very influential paper Davies and Pringle [143] discussed several regimes of the propeller.In their notation, the most standard case, typically discussed in the literature, is dubbed "supersonic propeller".In this regime, R m > R co and also typically R m ≪ R l .At the subsonic propeller stage R m < R co , accretion is still not possible as matter cannot cool efficiently.Finally, "a very rapid rotator" regime can be realized in the situation when R m ≲ R l , i.e. when the magnetospheric radius is just slightly below the light cylinder radius.
One of the main points, used and discussed by Davies and Pringle is related to the thermal balance in the envelope around the magnetosphere.The existence of a hot envelope results in a different pressure outside the magnetosphere with respect to the case of an efficiently cooled inflow assumed for calculation of R A .
At the supersonic propeller stage, the magnetospheric radius is: It is assumed that at this stage the energy released at the magnetospheric boundary is carried away by convection or turbulent motions.The turbulent velocity at any radius is assumed to be equal to the sound velocity.These assumptions result in the polytropic index n = 1/2 and ρ ∝ r −1/2 .Thus, we obtain: Matter is expelled from the magnetospheric boundary with the velocity significantly lower than the free-fall velocity.
At the phase of a very rapid rotator, the linear velocity of the rotating magnetosphere, ΩR m , is very high.This heats up the surrounding matter significantly.Effectively, this means that pressure is constant in the envelope up to "infinity", i.e. up to R G where properties are determined by the external medium (e.g., ISM in case of an isolated NS).This stage can be realized at low accretion rates or with high magnetic fields when the outer boundary of the heated envelope extends beyond R G .Otherwise, the pressure gradient becomes important and so, from the ejector stage the NS immediately come to the supersonic propeller stage.
At the very rapid rotator stage matter is expelled with the rotational velocity of the magnetosphere, but the density is low, so the spin-down rate is low, too: In addition to the propeller stage for which R m > R co , Davies and Pringle [143] discussed also a "subsonic propeller" stage with R m < R co (the stage was later re-considered e.g., in [146]).Here, despite the centrifugal barrier does not exist, matter cannot effectively accrete due to its high temperature.For low accretion rates, the subsonic stage can be quite long, and accretion can be postponed.However, instabilities at the magnetospheric boundary might allow the matter to accrete.This situation was recently analyzed in the framework of the so-called settling accretion model [147].Now, let us say a few words about the numerical modeling of the propeller stage.Wang and Robertson [148] used numerical simulations to model spin-down at the propeller stage.Specifically, they accounted for instabilities at the magnetospheric boundary.They obtained a very rapid spin-down: The pre-factor (2ζ WR )/(πβ) is expected to be > 1.Notice the similarity with the equation by Shakura [144], see eq. (39).What is also important, in their model R m is a function of the spin period, which is quite natural at the propeller stage, but this is typically not taken into account explicitly in many analytical models.
Numerical modeling of the propeller stage in the case of an isolated NS interacting with the interstellar medium was performed in [149], see also their earlier study in [150].Fitting of the numerical results provided the following scaling: K sd ∝ B 0.8 Ω 1.3 .This is a very rapid spin-down, similar to the one proposed in [144], where K sd ∝ Ω. Curiously, the dependence ∝ B 4/5 can also be obtained analytically, as discussed in [149].At first, the authors define the radius of the magnetospheric boundary from the condition And then, the spin-down torques is defined as Note, that here again R m depends on the spin period.
Francischelli and Wijers [151] compared several variants of the outflow velocity at the propeller stage and application to fossil discs around young NSs and to the source SMC X-1.They concluded that the highest velocity -ΩR m , -is not compatible with observations of SMC X-1.

Observations
Observational data about NSs at the propeller stage are not very abundant.Mainly, they are limited to compact objects in low-mass (LMXBs) or high-mass (HMXBs) X-ray binaries with a disc formation.Unfortunately, measurements of the spin-down rate at the propeller stage are not done, yet.Still, it is possible to determine the critical luminosity at which the transition to this stage occurs.
The first identification (not very secure at that time) of the propeller stage was probably made for the X-ray pulsar V0332+53 in [152].Some other candidates have been also identified in this paper, including 4U 0115+63 which was later studied by [153].This source demonstrated a decrease of the X-ray flux by more than two orders of magnitude in nearly one-half of a day.Another example -Aql X-1, -was presented in [154].This object is of special interest as it has a millisecond spin period.
In [158] the authors present observations of GRO J1008-57.The initial goal was to detect the transition to the propeller stage.Surprisingly, at some value (larger than the expected critical number) the flux decrease stopped.Thus, the authors conclude that the source instead of reaching the propeller stage, changed the mode of accretion.In this case, the disc at R m has temperature T < 6500 K. Interestingly, this regime is valid only for pulsars with relatively long spin periods: P > 36.6 B 0. 49  12 s.For smaller periods a usual propeller regime exists.
In [159] the authors consider various systems (LMXBs, HMXBs, cataclysmic variables, and young stars) for which transitions to/from the propeller stage have been proposed.Among the systems with NSs the authors discuss three LMXBs -SAX J1808.4-3658,IGR J18245-2452, and GRO J1744-28, -and three HMXBs -4U 0115+63, V 0332+53, and SMC X-2.In addition, three cataclysmic variables with accreting white dwarf and one young stellar object have been considered.
The authors tested the basic equation for the critical luminosity at which the transition to the propeller stage happens: L lim ∝ µ 2 P −7/3 R −1 .The analyzed systems provide the possibility to check this equation for a wide range of magnetic moments, spin periods, and radii of accretors.Surprisingly, the observed data perfectly confirm the dependences.In particular, the authors parametrized the equation as: L lim ∝ µ α P −β R −γ .And the fit gives: α = 1.9 +0.6 −0.3 , β = 2.3 +0.8 −0.4 , γ = 0.9 +0.3 −03 .More recently, the source GRO J1750-27 was proposed as a system with accretion/propeller transition [160].The authors observed a drastic drop in the X-ray flux accompanied by a decrease in the spin-up rate of the X-ray pulsar.Interpretation of this decrease in luminosity as transition to the propeller stage allows (for the known distance) determination of the magnetic field of the NS: ∼ 4 × 10 12 G.
An X-ray binary SWIFT J0850.8-4219 with a red supergiant donor is proposed to be at the propeller stage due to its low luminosity in comparison with expectations based on typical stellar wind mass losses from red supergiants [161].Still, this hypothesis deserves a more detailed analysis as the present-day data are not sufficient to make a firm conclusion.
Future instruments, like Athena [162], with large collecting area might allow detecting pulsations at low fluxes corresponding to the propeller stage.This will give an opportunity to test directly theoretical models of spin-down.

General properties of accretion flows onto magnetized NSs
NSs may accrete mass under various circumstances and from different sources, including ISM, gravitationally bound supernova ejecta (the case of fallback, see section 2), or a donor star.Interaction with the accreting matter is determined by the set of characteristic radii introduced in section 1: if both Shvartsman and magnetospheric radii are smaller than R G , mass is gravitationally captured.The flows inside the gravitational capture radius depend on the amount of linear and angular momentum in the matter captured by the star.If both are negligibly small, the star accretes spherically in the regime of Bondi accretion [163].Such a regime would be realized if the star is immersed in a static uniform gas.The gravity of the star creates a radial transonic inflow of matter that in its inner parts is modified by the presence of the solid surface and magnetosphere.Though unrealistic, Bondi accretion model allows to estimate the mass accretion rate in the case of subsonic motion through a uniform ambient medium according to formula (5).
Accretion from a supersonically moving uniform media may be thought of as gravitational focussing of the flow: gravity changes the direction of a streamline with a shooting parameter b by an angle ∼ R G /b, resulting in a shock wave of a roughly conical shape that dissipates some of the kinetic energy of the flow and allows for mass accretion onto the star from the downstream direction.This is an idealized case of Bondi-Hoyle-Lyttleton accretion [164].
An important factor that modifies the two scenarios mentioned above is the net angular momentum of the captured matter.For the net angular momentum of j the size of the disc is determined by the circularization radius If R c ≫ R m , the dynamics of the magnetosphere are determined by its interaction with the disc.In the case of Roche lobe overflow, j ∼ √ GM tot a, where a is the binary separation in the system, and M tot is the total mass of the system, hence R c ∼ a (though usually several times smaller).If the NS captures the stellar wind of its companion, j is much smaller, and the accretion flow is much better described by the Bondi-Hoyle-Lyttleton scenario (see however section 5.3.6).
All the accretion flows described above exist only outside the magnetosphere.Within the magnetosphere, magnetic stresses redirect the plasma along the field lines, creating funnel flows that start at R ∼ R m in the regions where the field lines interact with the inflow.Making assumptions about the shape of the magnetic field lines allows us to predict the shape of the funnel flow and the regions of the NS surface where the accreting matter finally falls and releases its gravitational energy.For a dipolar magnetic field, these accretion sites are associated with the polar caps or rings or arcs around them.If a dipolar magnetosphere acquires mass in its equatorial regions at a distance of R m , the matter, following a field line, hits the NS surface at a polar angle sin θ ∼ √ R NS /R m ∼ 0.01 (for the values typical for an accreting strongly magnetized NS, see section 5.3.1).
Accretion onto magnetized stars is a subject of many numerical studies, that are often designed for young stellar objects but are also applicable for NS accretion with an identical set of dimensionless parameters.To describe the accretion flow geometry and torques acting on an NS with a dipolar magnetic field accreting from a thin disc, the required dimensionless parameters are the so-called fastness parameter ω = Ω/Ω K (R m ) (where Ω K (R) = GM/R 3 is the local Keplerian frequency), R m /R NS , and magnetic angle χ.The rotation of the NS is assumed aligned with the rotation in the disc.We will talk about such simulations in more detail later in Sec.5.3.4.Simulations of Bondi and Bondi-Hoyle-Lyttleton accretion onto magnetized stars are much less numerous (see however the following section 5.2).

Isolated neutron stars
Accretion onto isolated NSs from the ISM was proposed more than half a century ago by Shvartsman [165,166] and Ostriker, Rees, and Silk [167] (the papers by Shvartsman were submitted slightly earlier).All of these authors correctly noticed that due to gradual spin-down or/and significant magnetic field decay an NS finally can start accreting the interstellar gas.In [167] it is estimated that for a typical pulsar magnetic field accretion onto a low-velocity NS can start in a few hundred million years due to spin-down and corresponding decrease of the pulsar wind power which prevents accretion at earlier times.
These authors applied Bondi accretion to estimate the luminosities of accreting isolated neutron stars (AINSs), see eq. ( 5).As they mainly considered low-velocity NSs (in 1969-1970 it was not known that compact objects could obtain large kicks), the estimated luminosity appeared to be ∼ 10 31 -10 32 erg s −1 .Shvartsman also accounted for a 'negative feedback' related to the heating of the external medium by the NS emission.This somehow lowered the expected luminosity and shifted the spectrum from X-rays to UV.
Interestingly, in the original papers [165,167] the authors did not discuss in detail the possibility of detecting AINSs as individual X-rays sources but mainly focused on the influence (heating and ionization) of the accretion luminosity on the surrounding medium.Later in the 1990s different authors intensively discussed observations of individual accreting isolated sources with space X-ray observatories.
Blaes et al. analyzed the observability of AINSs in several papers [168][169][170][171].In particular, in [169] the authors calculated the number of AINSs potentially detectable by ROSAT, and obtained that about several thousand such sources can be found.This very optimistic prediction is mainly due to the fact that the initial velocities of NSs were underestimated at that time.Despite the authors discussed some evolutionary aspects, finally, for their estimates they assumed that all old NSs accrete at the Bondi rate.This is, of course, an oversimplification.Reasonable accounting for spin evolution for realistic velocities of NSs has been done only in [172] and the realistic magnetic field distribution was added later by [173].In [171] the authors considered the influence of the AINS emission onto the surrounding medium and 'negative feedback' on a deeper level than it has been done in early studies.Naturally, this process reduces the number of potentially observable AINSs.
3D hydrodynamical simulations of accretion onto an INS were presented in [174].Magnetic field was not explicitly included in calculations but the authors modeled accretors of different sizes comparable with magnetospheric radius.The accretion rate was found to be reduced by a factor ∼ 2 in comparison with the Bondi-Hoyle-Littleton rate (eq. 5 with ξ acc = π).
Accretion onto isolated magnetized NSs was studied numerically in a series of papers by Romanova and her colleagues [149,150,175,176].The main result presented in [176] is related to the dependence of the accretion rate on the magnetic field of the NS.For magnetic field ∼ 10 12 G the authors obtain a reduction in the accretion rate by a factor ∼ 2 with approximate dependence Ṁ ∝ µ −0.4 .
An important question related to accretion onto INSs is about the spin properties of these objects.Naively, one can expect that an AINS is just spinning down, and the braking torque is ∝ µ 2 /R 3 co as the external medium does not bring, on average, angular momentum.However, the picture is more complicated due to turbulence.Accreting matter brings fluctuating angular momentum which is equal to zero only after averaging over a relatively long time ∼ R G /v. Up to our knowledge, for the first time, this situation was briefly considered in [177].In this paper the authors proposed that an AINS does not spin down continuously, but can reach some quasi-equilibrium period P turb ∝ µ 2/3 around which the spin fluctuates.The situation was considered in more detail in a brief note [178].
The evolutionary equation for the spin frequency can be written as: Here Φ is the turbulent torque and ⟨Φ⟩ = 0. Specific angular momentum j is determined as min , where v K is Keplerian velocity and v t is turbulent velocity taken at the gravitational capture radius.The coefficient k t is usually taken to be equal to 1/3.The evolution of the spin frequency can be represented as diffusion.Then for Ω turb = 2π/P turb we have: Here For j = v t( R G )R G we obtain: Here R t is the characteristic scale of interstellar turbulence for which v t (R t ) = 10 km s −1 and the turbulent velocity depends on the spatial scale as v t ∝ r 1/3 .Normalized units are I 45 = I/10 45 g cm 2 and µ 30 = µ/10 30 G cm 3 .
Numerically, the spin evolution of AINSs was studied in [179].The authors studied the spin evolution at the stage of accretion of a population of isolated NSs with realistic velocity and magnetic field distributions, but for a constant ISM density.After the onset of accretion, an NS is spinning down: dΩ/dt ≈ µ 2 /R 3 co .The influence of turbulence starts to be significant when the spin period reaches the value Typically this happens in ∼ 10 5 -10 7 yrs.After that, the period on average continues to increase, finally reaching P turb on a time scale ≳ 10 9 yrs, and fluctuating around this value.A typical time scale for Ṗ fluctuations is R G /v.The spin period distribution of AINSs appears to be very wide with typical values in the range ∼ 10 5 -10 7 s.Note, that the authors did not discuss the influence of shocks around AINSs on the angular momentum transfer by turbulent accreted matter.
As accretion onto INSs proceeds at low rates, the settling accretion regime might be applicable.We describe this mode of accretion in Sec.5.3.6.To the problem of AINSs this model was applied in [180].It was demonstrated that in the settling accretion regime, low-velocity NSs are expected to be transient sources with a typical duration of bursts of about a few hours and luminosity ∼ 10 31 erg s −1 .Persistent luminosities are expected to be three orders of magnitude lower, and so out of the reach of the present-day X-ray surveys.In addition, the authors proposed that inside convective envelopes AINSs can spin down much faster, and so they can reach quasi-equilibrium in a shorter time than it has been discussed in [179].
In the early 1990s, high hopes to detect AINSs were related to ROSAT observations [181].
Step by step, it became clear that due to different effects, the number of observable AINSs cannot be very high (see e.g., [182] and references therein).In the late 1990s, early enthusiasm was changed by pessimism related to the non-detection of any AINS by ROSAT (see a review in Treves et al. [183]).Finally, calculations presented in [172,184] demonstrated that it is impossible to expect even a few AINSs in the ROSAT data.
One of the most detailed searches for INSs in the ROSAT data was presented in [185].Up to now there are no known AINSs or good candidates.Still, the discovery of such sources will be a major breakthrough in understanding low-rate accretion and long-term evolution of isolated NSs.

Accretion in binaries 5.3.1. Variety of NSs in binary systems
Among the massive stars producing NSs, most are components of binary or multiple systems [186].Existence of a binary companion may affect the evolution of the NS in multiple ways, from tidal torques adjusting the rotation of the progenitor to mass exchange.Observationally, NSs in binaries manifest themselves as X-ray binaries (XRBs) and as millisecond pulsars (MSP).In both cases, the presence of a binary companion is crucial for the rotational properties of the NS.MSPs are radiopulsars thought to have accretion episodes in their past [187].The connection between the two classes of systems is confirmed by the existence of accreting millisecond X-ray pulsars (AMXPs, Archibald et al. [188]).We will cover these classes of objects later in section 5.3.7.
Classification of XRBs may be done in different ways, using such criteria as mass transfer type (Roche-lobe overflow or wind capture) or magnetic field of the accretor, but the traditionally adopted distinction is between high-mass and low-mass XRBs (HMXBs and LMXBs, respectively), with the mass of the donor being the crucial parameter.Higher donor mass (several Solar masses or larger) implies a relatively young age of the NS, a large magnetic field, and often mass transfer through the wind.Low donor mass makes the formation of the binary an unlikely scenario, as progenitor systems should be disrupted by supernova explosions.Thus, many LMXBs, especially in dense stellar environments like globular clusters, should have been formed by exotic processes such as binary companion exchange [189], or require a finely tuned natal kick.In terms of the number of active systems, the low formation probability is compensated by the abundance of low-mass stars and their longevity.The NSs in such systems are on average much older, which implies low magnetic field strengths.Also, accretion from a low-mass donor is usually associated with Roche-lobe-overflow mass transfer.
Low magnetic fields of LMXBs make it difficult to study their rotational evolution.With the exception of AMXPs, their spin periods are inferred indirectly from the other variability patterns such as post-burst oscillations in X-ray bursters [190] or quasi-periodic oscillations in the sources lacking coherent periodic pulsations.The latter includes the so-called 'Z' and 'Atoll' sources introduced by [191], containing most of the known LMXBs (see [189] for a review.)It is however unclear if the QPO frequencies are reliable indicators of the spins of these objects, see Méndez and Belloni [192]).
Because the matter accreted by a non-magnetized NS through an accretion disc has the net angular momentum equal to Keplerian at the inner disc edge R in ≃ 10km (determined either by the surface of the star or by its gravity field), angular momentum is gained at the time scales only several times smaller than the time scales of mass growth.For accretion at the Eddington limit ( Ṁ ≃ 10 18 g s −1 ), spinning up to a several millisecond period would take millions of years.
Strong magnetic field makes the spin-up process much more efficient, as a large magnetosphere intercepts accreting matter with a larger (∝ √ R m , if accretion runs through a disc) net angular momentum.On the other hand, rapid rotation in combination with strong magnetic fields suggests pulsar torques and pulsar wind formation.Even more importantly, rapid rotation of the NS may lead to propeller regime and thus quenching of mass transfer (see section 4).Hence, the spin evolution in HMXBs is usually thought of as deviations from an equilibrium between spin-up and spin-down torques.The equilibrium point shifts with the changing mass accretion rate and some other parameters, likely related to the geometry of the accretion flows.Many HMXBs are observed as X-ray pulsars (XRPs, i. e., sources of coherent periodic variability in the X-ray range, clearly related to rotation).Studying X-ray pulsations allows to probe their rotational evolution in real time and on long archival data series.One of the most important projects of this kind is the monitoring program held out with Fermi/GBM in the hard X-ray range [193].The series obtained by Fermi/GBM span more than 10 years of observations on 39 X-ray pulsars.
The majority of X-ray pulsars are members of binaries with Be-class donors (Be/X-ray systems, Raguzova and Popov [194], Reig [195]).The population of such systems is diverse, with the orbital periods in the range 20 − 200d and NS spins 1 − 10 3 s.Orbital and spin periods are known to correlate with each other [196], supporting the idea that that their rotation is close to equilibrium between different torques acting on the star (see section 5.3.4).The exact origin of this correlation is however unclear.Mass transfer rate in Be/X-ray systems is strongly modulated by orbital eccentricity and other factors such as the mass accumulation in the decretion disc of the star.This results in a complex phenomenology of periodic and aperiodic flares, during which the NS is spun up by accreting matter, and quiescent states where the rotation slows down.This picture is not universal and allows for objects that are, for example, spun up even during the quiescent state, suggesting that the equilibrium is reached beyond the time span of observations.For a given object, however, Ṗ is well correlated with the luminosity [193].
Another important category of XRPs are objects accreting from the winds of earlytype supergiants.The prototype and the best studied object of this class is Vela X-1 [197].Accretion disc may form in such systems but not necessarily, as the circularization radius is not necessarily larger than the radius of the magnetosphere.
In a homogeneous wind, the net angular momentum is (see for instance Illarionov and Sunyaev [139]) This angular momentum is aligned with the rotation of the binary.Strong inhomogeneity of the wind would create an additional randomly oriented net angular momentum component j w ∼ R G v w ≫ j 0 , where v w ≫ R G Ω orb is the velocity of the wind.Hence, for such systems, the angular momentum of the accreting matter is likely to be randomly aligned.

Spin-up by an aligned disc
First, let us consider a planar thin conducting accretion disc interacting with the dipolar magnetic field of a star.Even in its simplified version, the problem of disc-magnetosphere interaction does not have a general analytical solution.A thorough review of the existing models and approaches to disc-magnetosphere interaction was made, for instance, by Lai [198].
The general picture is dictated by the existence of a magnetosphere and the strong radial decrease in magnetic tensions.One might expect magnetic field lines threading the disc within some depth ∆R that is determined by intermediate-scale processes such as turbulence and reconnections between the stellar magnetic fields and the fields frozen into the accreting matter.Large ∆R ≳ R models are usually referred to as 'magnetically threaded disc', or MTD [199,200] models.Within the penetration region, the internal angular momentum flow in the disc, normally governed by effective viscosity forces [201], is contributed by the torque from the distorted magnetospheric lines.Purely poloidal magnetic field of an aligned dipole does not produce any torque, but the relative motion of the magnetosphere with respect to the disc creates toroidal field components B ± φ , having different signs on the different surfaces of the disc.Toroidal fields should also change sign at the corotation radius and are likely of the same order as poloidal fields.Magnetic field torque acting on the NS has two contributions that are not clearly distinguishable: the angular momentum coming with the accreting matter itself and the torque produced by Maxwell stresses B z B ± φ /4π (where B p is the poloidal field).As the flow inside the magnetosphere is bound to corotate with the NS, its angular momentum is also primarily transmitted by magnetic stresses.While the net angular momentum at the magnetospheric boundary is ΩR 2 m , its value near the surface of the NS is ∼ ΩR 2 NS , that is normally several orders of magnitude smaller.The torque acting on the surface of the NS from the disc may be calculated as the star-disc interaction (SDI) torque, containing both spin-up and spin-down parts where integration is performed over the surfaces of the disc threaded by magnetic field lines (usually, upper and lower surfaces).Straightforward integration of this expression requires knowledge about both field components.In general, the SDI torque may contain a spin-up part created by the field lines threading the disc inside the corotation radius, and a spin-down part from the lines extending outside (their tensions marked with green and red arrows in Fig. 4).Magnetic field opening and screening (see below) tend to decrease the contribution of the outer parts of the disc, making SDI torque essentially a spin-up torque.
The MTD approach assumes that the poloidal field interacting with the disc is equal or at least close to the original dipolar field in a broad range of radii ∆R ≳ R, and the toroidal fields are generated by the difference in the rotation velocity between the star and the disc.The expression used for the toroidal fields has the form The last condition is related to the stability of predominantly toroidal force-free magnetic field configurations.In these assumptions, the torque related to disc-magnetosphere interaction is convenient to normalize with the accretion torque spinning the NS up, This kind of normalization may also be used in a more general case when many different aligned torques affect the rotation of the star at the same time.The simplified picture where all the matter accreting through the disc gives its angular momentum at R m to the star corresponds to n(ω) = 1.In general, the dimensionless factor n(ω) accounts for all the processes of angular momentum exchange between the disc and the star [200,202].The classic MTD model assumes that the magnetic field retains its poloidal component well outside the radius of the magnetosphere, as one would expect if magnetic diffusivity in the disc were high.Rotation of the disc creates additional azimuthal components scaling with the poloidal dipolar field components B φ ∼ B p .Integrating the torques over the surface of the disc, one yields for the torque (assuming B φ = ±B z ) Such a configuration allows for a rotational equilibrium, where part of the disc (inside the corotation radius) spins the star up, while the outer regions spin the star down (see Figure 4).At a constant mass accretion rate, this kind of equilibrium keeps the corotation radius close to the magnetospheric radius.
Decreasing the penetration depth restricts the integration range in (52) to a narrow strip where the rotation frequency of the disc is likely to be faster than Ω.The spin-up torque should still be close to the accretion torque (53), because of the angular momentum advected with the matter.However, as we show later in section 5.3.4,this spin-up is opposed by other forces, not directly related to star-disc interaction.

Spin-up of an inclined dipole
The very existence of an XRP phenomenon requires an NS with misaligned rotational and magnetic axes.Observational estimates of magnetic angles χ (see section 3) in X-ray binaries cover a broad range, from nearly aligned objects towards nearly orthogonal (see in particular the recent results from IXPE Doroshenko et al. [203]).The observational estimates for magnetic angles in MSPs are strongly model-dependent [204] but also favour a broad distribution.For all these reasons, it is important to be able to relax the alignment condition for the dipolar magnetic field of the star.At the same time, its rotational axis is likely to be aligned with the spin-up torque, with a couple of important exceptions we consider later in section 5.3.5.Relaxing the magnetic obliquity leads to several consequences.First, the condition for the magnetospheric boundary is now satisfied at different radii for different azimuthal angles.According to Bozzo et al. [205], this leads to a correction factor related to the magnetic field strength dependence on the polar angle.Even more importantly, magnetic inclination breaks the axial symmetry of the flow, and the funnel flow becomes concentrated to a couple of streams restricted in the azimuthal direction.
An MTD approach allows calculating the spin-up torque, assuming that the magnetic field interaction with the plasma of the disc creates a non-dipolar component (generalization of B ± φ in the aligned case) and consequently angular momentum exchange between the disc and the star.Accurate calculation in the MTD approach was performed by [206] and involved three surfaces: upper, lower, and inner faces of the disc.The unperturbed magnetic field line of an inclined dipole was complemented by an additional field component directed along the relative velocity difference between the disc and the magnetosphere.Integration over all the surfaces leads to a torque that changes sign at some critical fastness parameter values close to unity.

Spin-down torques and equilibrium
The MTD picture described above has a major caveat related to the opening of the field lines due to asynchronous rotation between the disc and the star [198,207].A more realistic picture involves field lines loaded with mass and extending to infinity, which creates possibilities for mass ejections.A lot of inspiration here comes from the studies initially aimed at magnetized young stellar objects that evolve under a similar combination of torques.In Matt and Pudritz [208], the existence of magnetized winds draining the star of angular momentum and carrying away some part of the accreting matter was proposed as an explanation for the spin period distribution in T Tau stars.Many detailed 3D MHD simulations [209][210][211] designed for young star magnetospheres are scalable for the smaller and rapidly rotating magnetospheres of NSs with identical dimensionless parameters (fastness ω, magnetic angle, and disc thickness), that allows to make implications about XRPs as long as relativistic effects are ignored.All the simulations reproduce a moderate penetration depth of the order disc thickness, field line opening and magnetized ejections.The estimated torque related to the ejections is of the order K sd ∼ ΩR 2 m , similar to the braking torque given by the propeller theory K sd ∼ µ 2 /R 3 co .The exact value contains correction factors depending on the fastness parameter [211].Fully relativistic MHD simulations [212,213] confirm the general picture, but are so far limited only to very compact magnetospheres.
In an accreting magnetosphere, some of the field lines remain open and free of the matter supplied from the disc.As these lines tend to corotate with the star and extend beyond the light cylinder, they create a braking torque similar to the pulsar torque considered in section 3. Particle creation due to pair cascades is in this case probably suppressed, but the contribution from the vacuum losses remains as long as there are open magnetic field lines.On the other hand, interaction with the disc changes the field shape outside the radius of the magnetosphere, and the resulting electromagnetic torque is likely to be higher than in the case of an ejector.In the aligned case, Parfrey et al. [20] estimate the value of the torque as where ζ ∼ 1 is a dimensionless parameter, and K sd corresponding to the pulsar braking torque given by equation (29).Three-dimensional MHD simulations confirm the picture drawn by Parfrey et al. [20] and also show that the open field lines of the accreting magnetosphere behave like Poynting-dominated jets in accreting black hole systems powered by the rotation of the compact object [212,213].
It is unclear how applicable these results are to a more general geometry.In particular, is this spin-down torque aligned with the rotation of the NS or with the disc that opens up the field lines?How does this particular spin-down torque depend on the magnetic inclination of the rotator?Answering these questions would require full-scale simulations with non-zero magnetic and disc inclination angles.
Combined, two or more torques of different signs, dependent on mass accretion rate and the spin on the NS, result in a rotational evolution equation which has a single equilibrium point at n = 0. Different models using different assumptions about the spin-down torques differ in the dimensionless factor n(ω).The solution n(ω) = 0 is usually close to the propeller limit ω = 1, or, in terms of spin period P eq ≃ 2π(GM) 5/7 µ 6/7 Ṁ−3/7 ≃ 5 µ 10 30 G cm 3 6/7 Ṁ 10 16 g s −1 −3/7 s. (57) Applying this expression to an XRP with a known spin period and a reliable estimate of the mass accretion rate (that is proportional to the bolometric accretion luminosity) allows to estimate the magnetic moment in the assumption of equilibrium.In a strong flare, when the spin-up torque dominates over the spin-down one, the magnetic moment may be probed by the observed spin-up rate.Several methods of magnetic field estimation using the spin-up and spin-down torques are described in Chashkina and Popov [214].
Depending on the variability of the mass supply rate and potentially other parameters, some objects might oscillate between accretor and propeller regimes, while others proceed to accrete in the low-accretion-rate state.There is observational evidence that some Be/Xray systems enter the propeller stage (see section 4), while others retain a small persistent mass accretion rate during quiescence [158].The reason for such a dichotomy is unclear and may be connected to additional parameters such as magnetic angle, or with the mass accretion rate variations in the disc.
The time required to reach the equilibrium state may be estimated as eq Ṁ(GM) 2/3 ∼ 10 3 I 10 45 g cm 2 Ṁ 10 18 g s −1 −3/7 yr, (58) where the radius of the magnetosphere was assumed R m ∼ R co (that allows to exclude the value of Ω eq = 2π/P eq from the equation), that is a reasonable assumption near the equilibrium.Cases, when the spin period changes considerably during the observational period, are scarce and related to unique and probably catastrophic events.A bona fide example of such a case is the super-Eddington XRP NGC300 ULX1 [215] which was spun up by about a factor of five during the period of observations.Detailed timing studies of XRPs allow testing of the models of magnetospheric accretion flows.In particular, the flares of Be/X-ray systems are often thought of as the case when spin-up torque dominates over all the spin-down mechanisms, hence , is generally supported by observations [216].There are indications for a steeper dependence at higher luminosities [217] that might be a signature for a more elaborate disc-magnetosphere interaction [158,218].

Misaligned spin-up torques
If the spin-up torque does not change its direction (in particular, if its direction is fixed by the orbital plane of the system), one would expect the rotation axis of the NS to align with the torque within the spin-up time (58).But on some of the evolutionary stages and in certain accretion regimes in X-ray binaries the rotational axis of the star may deviate strongly from the external spin-up torque.First, before the mass exchange, the accreting NS is likely to have a rotational axis misaligned with the orbital plane due to the kick received during the supernova explosion [219] or a three-body process responsible for the creation of the system.
Another important case is accretion from an inhomogeneous wind when the density and/or velocity gradients within the wind make the direction of the torque essentially decoupled from the orbital plane [220].Numerical simulations of such an accretion regime were done, for instance, by El Mellah et al. [221].The spin of the NS in this case is lower than the equilibrium value, and its direction is involved in a random walk adjusting to the rapidly changing external torque.The observed spin derivative in this case should be weakly correlated with the flux, because of the relative orientation of the torque and the angular momentum.Such a behavior is observed for Vela X-1 [197] and other XRPs fed by supergiant winds.
Misaligned torques are also important as a possible trigger for magnetic angle evolution.Adding the spin-up torque to Euler equations of rotational dynamics [222] reveals an additional term proportional to sin χ cos χ sin α cos α, where α is the misalignment angle between the rotational axis of the NS and the external spin-up torque, that increases the magnetic angle during the alignment stage.As a result, the alignment of the NS rotational axis increases the magnetic angle, increasing the possibility of observing the NS as an XRP.Wind accretion with a randomly aligned spin-up torque results in a random walk evolution of χ, increasing the number of strongly oblique pulsars.

Quasi-spherical accretion regime
If the circularization radius becomes smaller than R m , the Keplerian disc can not form, and the torques acting on the NS should be revised.The net angular momentum of the matter interacting with the magnetosphere is now much smaller than Keplerian at the magnetosphere radius, meaning that the spin-up torque is severely diminished.For a homogeneous wind, j ∼ Ω orb R 2 G [139].The matter near the magnetospheric boundary is heated by shock waves to temperatures close to virial, which implies a quasi-spherical envelope around the magnetosphere.
The rotation of the NS adjusts to the rotation of the envelope via magnetic stresses at the edge of the magnetosphere.
The structure of the envelope depends on the presence of cooling and angular momentum transfer within it.In the case of weak cooling, an elaborate theoretical model was developed by Shakura et al. [147].As the main process responsible for angular momentum transfer, the authors considered an anisotropic turbulent velocity field.Anisotropy allows to change the direction of the angular momentum transfer [223].In particular, predominantly radial motions lead to a flat distribution in angular momentum, meaning that the rotation frequency in the envelope Ω env ∝ R −2 .Shakura et al. [147] adopt a more general power-law scaling for the rotation frequency.
Equilibrium spin periods predicted by the quasi-spherical model are much longer than for disc accretors, roughly corresponding to the equal net angular momenta of the outer magnetosphere rim and the captured matter ΩR 2 m ≃ Ω orb R 2 G , implying, up to a coefficient of the order unity, P eq, sph ∼ 10 3 v 10 8 cm s −1 For identical magnetic moments and mass accretion rates, this value is about two orders of magnitude longer than the equilibrium spin for disc accretion (compare equation 57).Such an accretion scenario may be applied to wind-fed HMXBs like Vela X-1, to symbiotic wind-fed LMXBs like GX 1+4 [224], and potentially to some of the Be/X-ray systems having very long spin periods that require suspiciously large magnetic fields if interpreted as disc accretors [214].

Recycled and transitional pulsars
About 10 percent of the RPP population belongs to millisecond pulsars (MSPs), a category of short-period, weakly magnetized NSs with the estimated ages of about Gyrs, well beyond the expected period of pulsar lifetime [187].This population is understood as recycled, i. e. spun up by accretion in a binary system to a period short enough to support a pulsar cascade.The majority of MSPs have binary companions, and some of them have masses in excess of 1.8M ⊙ .
Among the NSs with reliably measured rotation periods, the fastest rotator is PSR J1728-2446ad [225] having a rotation frequency of more than 700Hz.Several other MSPs are rotating at similar rates ∼ 500 − 600Hz.There are indications that some non-pulsating NSs in LMXBs have frequencies around 500-600Hz [226].These values are several times smaller than the Keplerian rotation rate (∼ 2.5kHz for a 2M ⊙ star with a radius of 10km), which is sometimes interpreted as a signature of an additional braking mechanism.If the excess ∼ 0.3 − 0.5M ⊙ of mass observed in MSPs was gained via disc accretion, the expected rotation frequencies should be pretty close to critical.For I ≃ 10 45 g cm 2 and total angular momentum gained via accretion ∆M √ GMR NS ∼ 10 49 g cm 2 s −1 , this implies a spin period of P s ≲ 0.6ms, or a frequency ∼ 1.5kHz.One factor limiting neutron star rotation frequency is gravitational waves [227,228].The upper spin limit may be also related to the electromagnetic braking term (see Parfrey et al. [20] and sec.5.3.4), or simply to the decrease in spin-up efficiency with the shrinkage of the magnetosphere [229].Population synthesis simulations run by Tauris [230] shows that in certain assumptions the observed spin limit may be interpreted as a consequence of the braking torques slowing the NS down on the propeller stage.
Our understanding of MSPs as a product of binary evolution is supported by the existence of LMXBs containing NSs with properties similar to MSPs (accretion-powered millisecond pulsars, AMXPs, see Patruno and Watts [231]) and, most importantly, transitional MSPs (TMSPs) that combine the radio emission typical for RPPs with an evidence for an accretion disc [232].Observations show TMSPs switch between different activity states.The prototypical object PSR J1023+0038, in particular, is observed either as a radiopulsar with an X-ray luminosity L X < 10 33 erg s −1 , or in an X-ray bright state during which it has pulsations in the X-ray range, but is not observed in radio [233].X-ray pulsations are observed in both phases, but the pulse shapes are different.In the AMXP state, the pulse profiles are dominated by the first overtone.An obvious interpretation of this behavior is the transition between ejector and accretor states (see sec. 1).Pulsating X-ray emission is interpreted as a signature of a funnel flow within the magnetosphere.An alternative explanation proposed by Veledina et al. [234] considers that the source is always rotationpowered, and the X-ray emission observed during the bright states is produced by the disc existing outside the light cylinder and reprocessing part of the pulsar wind power.Anisotropy of the pulsar wind leads to a modulation of the X-ray emission with a doubled spin frequency of the pulsar.

Magnetic field burial by accretion
Ohmic magnetic field decay given by equation ( 12) predicts a late-time exponential decay on millions of years, but is definitely inapplicable to objects as old as hundreds of millions of years, such as MSPs and AMXPs.The situation is unclear, and most studies assume that Ohmic losses are negligibly small during late evolution (see discussion in [54]).
The proposed long accretion histories and large accreted masses in MSPs imply that most of the crust material in these objects is renewed.This should also affect the magnetic field strength and structure.The idea of the field being screened by a layer of accreted matter was proposed already in Bisnovatyi-Kogan and Komberg [235].It is known to be an important effect in other circumstances, such as fallback accretion (see sec. 2).
Due to the presence of the magnetosphere, an XRP is primarily accreting onto a relatively small fraction of its surface, about the size of the two polar caps (∼ 2πR 4  NS /R 2 m ) or smaller.Magnetic field tensions are able to support only a small layer with a surface density of where g is surface gravity.This led to the suggestion that the newly accreted matter spreads over the surface of the NS from the magnetic poles toward the equator, screening the field and distorting the shapes of the field lines [236].In an idealized axisymmetric approach, the dipolar component of a buried magnetic field decays approximately quadratically with the accreted mass.For instance, Payne and Melatos [237] approximate their results with a power-law relation Accretion of only about 10 −4 M ⊙ would lead to a magnetic moment decrease by a factor of several.Local magnetic fields near the magnetic equator in this scenario are not buried and are even expected to grow due to flux conservation, forming a "magnetic tutu".The increase in field gradients, both vertical and latitudinal, increases Ohmic losses, which was also supposed to contribute to the long-term evolution of the field [238].
It should be noted that the spreading of the accreted matter over the surface of the NS is accompanied by a number of instabilities acting on different scales and caused by different processes [187].Local interchange instabilities [239,240] operate on very small time scales (several local dynamics) and effectively result in magnetic field diffusion in the favorable parts of the flow.Also, the strongly magnetized layer buried under a lessmagnetized layer of accreted matter is susceptible to Parker instability [241].The resulting magnetic fields are suppressed, decayed, and have a much more complex topology than the initial dipole.

Exotic and hypothetical stages
In this section, we briefly discuss several stages of magneto-rotational evolution that remain hypothetical, i.e. no direct evidence in favor of their existence is known.

Georotator
The name of this stage is related to the fact that the shape of the magnetosphere and some processes are similar to the case of the Earth magnetosphere in the Solar wind.The stage with R m < R co but R m > R G was mentioned already in [139].The authors just note that in this case matter flows around the magnetosphere as the gravitational influence of the compact object is too small.Such a situation can be realized if the velocity of matter relative to the NS is too high (e.g., due to very fast stellar wind in a wide binary, or due to a large spatial velocity of the compact object in the case of an INS), if the magnetic field is too large, or if matter density is too low (this can happen e.g., for an INS far above the Galactic plane).
It is worth to note, that the condition R m < R co is important.Correspondingly, the georotator stage is an alternative to the stage of accretion, not to the propeller stage.Thus, one can consider evolution as Ejector→Propeller→Georotator, or switching from accretor to georotator.The latter scenario can be realized e.g., while an INS moves in the Galactic potential rising high above the Galactic plane and then coming back.
Also, we note that the magnetospheric radius is defined by different equations in the case of accretion, at the propellers stage, and in the case when R m > R G .That's why the critical condition for transition to/from the georotator stage might be carefully calculated accounting for necessary effects.If we consider a typical INS with µ ≈ 10 30 G cm 3 and the ISM density n ≈ 1 cm −3 , then the critical velocity for transition to the georotator stage is ∼ 400 − 500 km s −1 .Note, that for such parameters duration of the ejector plus propeller stages is longer than the Galactic age.Still, e.g.INSs with large magnetic fields can become georotators.
To our knowledge, spin-down at the georotator stage was never analyzed in detail.From general considerations, it can have some similarities with the subsonic propeller stage [143] or/and with the settling accretion [147].When a long spin period is reached, the influence of turbulence on the rotational frequency might become important [178,179], similar to the case of accreting INSs.
The case R m > R G was studied numerically for non-rotating INSs in [175].The authors dub this stage a 'magnetic plow'.In this study, just the case of aligned magnetic axis and spatial velocity vector was analyzed.It is demonstrated that a tiny amount of matter still can accrete onto the NS surface, but the energy release is small and hardly can be detected.Instead, the authors focus on the energy dissipation in the bow shock and in the magnetic tail.
In the bow shock, the energy release can be estimated as: It is expected that the energy is emitted in X-ray and, maybe, in softer bands.The case of energy dissipation in the magnetic tail is more interesting as different regimes can be realized.If the tail is not expanding, i.e. even far from the NS its cross-section has the size ∼ R m , then the total energy is: Here l is the length of the tail.This energy can be dissipated due to reconnection on the dynamical time scale Energy can be also dissipated in bursts in a volume ∼ πR 3 m , so the emitted amount is: Then the power is about E rec /t A , where t A ≈ R m /v A and the Alfven velocity v A = B(R m )/ 4πρ.In (66) But notice the dependence on n which is low.
In the case of reconnection in the tail, the energy might go into accelerating electrons and then it can be emitted in the radio band.
Potentially, the tail can even reach the light cylinder.So, a region of opened field lines can be formed.It is tempting to call such an object a geoejector.However, a particleproducing cascade in the magnetosphere might be switched off as the NS left the ejector stage long before it reached the georotator stage.Still, elongated structures due to large spatial velocities can accompany radio pulsar activity, see a brief discussion in [175], sec.6.5.

Pulsating magnetospheres and other exotic regimes
Two exotic transient regimes with a pulsating magnetosphere can be realized at the propeller and the ejector stages.
It was hypothesized (see [1] and references therein) that for effective cooling, a dense envelope can grow around the magnetosphere at the propeller stage.Equilibrium at the magnetospheric boundary can be described as: Here M sh is the mass of the accumulating shell.As it grows, the magnetosphere shrinks until it reaches R co .Then the whole envelope collapses on the free-fall time scale.A similar regime of transient accretion can be realized also in the framework of the settling accretion model [180].In this case, episodes of accretion are due to the development of the Rayleigh-Taylor instability at the magnetospheric boundary.
Another transient regime proposed in [242], can be related to the stage of ejector.It can be realized e.g., in a dense stellar wind of the companion in a binary system.Relativistic wind produced by a pulsar can be 'captured' by the stellar wind.In this case, a cavern is formed around the NS.With time the pressure inside the cavern grows and so it is expanding.Pressure in the surrounding medium is not uniform and is decreasing outwards.Thus, at some point, the cavern opens.This might produce a modest radio flare.This scenario was proposed to explain properties of the Galactic center radio transient GCRT 1745-3009 [243].Such pulsating caverns can be also formed around isolated NSs.
Finally, we want to mention again a regime of enhanced rotational energy losses at the radio pulsar stage when the configuration of magnetic field lines is modified by an accretion flow.This model was proposed by Parfrey et al. [20] and we described it in some detail above (Sec.5).This is again a hybrid stage as a relativistic wind from the compact object and an accretion disc inside a light cylinder around coexist.It might be possible to identify this stage (if it is realized in nature) in binary systems.

Discussion and conclusions
Spin properties of NSs -period, period derivatives, their correlations with other parameters e.g., flux, etc. -are often well-measured and provide a model-independent source of information about these compact objects.Thus, it is important to know how these properties are related to other parameters (e.g., magnetic field, accretion rate, etc.) and how they evolve in time.
Unfortunately now, even in a clearer case of radio pulsars, we do not have a complete understanding of how spin parameters are related to other properties.Thus, above we provided a review of different models of spin behavior of NSs at various stages.
Lack of understanding is related partly to complicated physical processes at work, and partly to the absence of detailed observational data.Theoretical questions are related to the unknown equation of state (EOS) of NS interiors (see e.g., a review in [244] and references therein), to poorly known topology and evolution of the magnetic field (see reviews [4,245]), and to many uncertainties in plasma-magnetic field interaction (see [246,247]).From the observational side, for example, there is no final description of the radio pulsar stage.On the other hand, there are no detected isolated accreting NSs or NSs at the stage of georotator.This situation results in many different approaches to describe the spin evolution of NSs.
It is tempting to say that in the near future, thanks to advances in 3D MHD modeling of processes around NSs and with the appearance of new observational facilities many puzzles of magneto-rotational behavior and evolution of NSs will be solved.However, more than half a century of studies in this field of research suggests that we have to be more realistic.It seems that many approaches suggested already in the 1970s despite all the uncertainties and simplifications will be used in the forthcoming years to derive properties of NSs from spin measurements.
We conclude by stating that despite our understanding of the spin behavior of NSs is far from being complete, there are many useful approaches suggested by many authors, and step by step we approach a better description of the magneto-rotational evolution of NSs.

Figure 1 .
Figure 1.A schematic representation of relation between critical radii at different stages of magnetorotational evolution.

Figure 2 .
Figure 2. A schematic representation of spin evolution of an isolated NS.The compact object consequently passes the stages of ejection, propeller, and accretor.Finally, it can reach the stage when the spin behavior is determined by the accretion of turbulent angular momentum.

Figure 3 .
Figure 3. Spin periods and their derivatives for 3473 pulsating neutron stars collected in the ATNF Pulsar Catalogue.Middle panel: Classical single rotation-powered pulsars are indicated by grey dots; red circles mark sources with high-energy pulsed emission; green crosses indicate radio quietness, while horizontal blue stripes are for pulsars in binary systems.Lines of constant surface magnetic fields B sd ∝ √ P Ṗ Gs, ages τ sd = P/2 Ṗ and luminosities L sd ∝ ṖP −3 are also drawn assuming R NS = 10 km and I = 10 45 g cm 2 .Bottom panel: Sources with negative Ṗ are shown in the same coordinates.Most of these pulsars are in globular clusters and their "negative spin-down" is due to the Shklovskii effect[37].The pulsar graveyard border is plotted according to[38].Upper panel: Vertical lines show the spin period values for pulsars for which Ṗ has not yet been estimated.Grey lines are for all such objects, green ones mark the values for radio-quiet sources, red ones are for those exhibiting high-energy emission, while blue lines are for pulsars in binary systems.

Figure 4 .
Figure 4. Three-dimensional sketch illustrating the MTD model of field penetration into an aligned disc.The blue circle in the disc plane marks the corotation radius.The field lines have dipolar and toroidal components, the latter related to the motion of the disc and limited by the absolute value of the seed poloidal field.Arrows show the rotation of the disc and the tension along the field lines.