Confronting the neutron star population with inverse cascades

The origin and evolution of magnetic fields of neutron stars from birth has long been a source of debate. Here, motivated by recent simulations of the Hall cascade with magnetic helicity, we invoke a model where the large-scale magnetic field of neutron stars grows as a product of small-scale turbulence through an inverse cascade. We apply this model to a simulated population of neutron stars at birth and show how this model can account for the evolution of such objects across the $P\dot{P}$ diagram, explaining both pulsar and magnetar observations. Under the assumption that small-scale turbulence is responsible for large-scale magnetic fields, we place a lower limit on the spherical harmonic degree of the energy-carrying magnetic eddies of $\approx 40$. Our results favor the presence of a highly resistive pasta layer at the base of the neutron star crust. We further discuss the implications of this paradigm on direct observables, such as the nominal age and braking index of pulsars.


INTRODUCTION
Neutron stars harbor the strongest known magnetic fields in the universe, with the large-scale poloidal field strength B p in the so-called magnetars exceeding values of B p ≈ 10 15 G. The magnetic field strength plays a crucial role in determining the observational properties of a neutron star and whether, for example, it will be seen as a standard radio pulsar or a magnetar (e.g., Borghese 2023). Furthermore, the large-scale dipolar field is thought to be primarily responsible for the spindown of neutron stars (Ostriker & Gunn 1969), leading to a spindown law of the formΩ ≈ Ω n , where Ω = 2π/P is the spin frequency of the star with P being the rotation period, and n is the so-called braking index that can be observationally constrained from the combination n =ΩΩ/Ω 2 . For dipole magnetic spindown, one has n = 3, with deviations from this value expected if other mechanisms are driving the spindown (e.g., n = 1 for winds, and n = 5 for gravitational wave emission from quadrupolar 'mountains').

nikhil.sarin@su.se
In this picture, neutron stars are assumed to be born rapidly rotating, in the upper left corner of the PṖ diagram, and evolve towards the bottom right corner and the so-called 'death' line as they spin down with a constant magnetic field. Observed braking indices are, however, often different from n = 3, and there is significant observational and theoretical evidence for magnetic field evolution during the lifetime of a neutron star; see Igoshev et al. (2021b) for a recent review. Therefore, understanding neutron star formation and evolution mechanisms is essential to understanding the observed neutron star population.
The origin of neutron star magnetic fields is, however, still debated. It is generally agreed that the fossil field inherited from the progenitor star is insufficient to explain the strongest neutron star fields; therefore, some field amplification is necessary. Most models invoke a large-scale dynamo at birth, powered by convection or differential rotation (Thompson & Duncan 1993). The source of the turbulence driving the dynamo comes from neutrino-driven convection. An additional source of turbulence is the magnetorotational instability Reboul-Salze et al. 2022), which draws energy from the differential rotation, i.e., ulti-mately from potential energy. This instability requires the presence of a magnetic field, which is accomplished through a positive feedback loop whereby magnetic energy can be amplified further by dynamo action (Brandenburg et al. 1995). Other possible dynamo mechanisms tap energy through fallback accretion (Barrère et al. 2022) or through precession (Lander 2021). In all those cases, small-scale and large-scale dynamo action appear together. In this Letter, we define a largescale field as the part contained in the dipole component, while the small-scale field is the remainder.
Magnetic field amplification can also occur later in the life of a neutron star due to the re-emergence of a buried magnetic field due to Ohmic dissipation (Muslimov & Page 1996;Ho 2011) or Hall drift (Gourgouliatos & Cumming 2015. Most of these models of magnetic field evolution in neutron stars assume initially low-order multipoles (Gourgouliatos & Hollerbach 2018;Gourgouliatos et al. 2020), although the true magnetic field structure is a matter of considerable debate (Igoshev et al. 2021b). Alternatively, the magnetic field could be predominantly of small scales at a very young age of the neutron star, but the field could then undergo what is known as inverse cascading (Wareing & Hollerbach 2009, 2010Igoshev et al. 2021a). The spectral magnetic energy at small wavenumbers or low multiples would then increase with time rather than decrease. This idea was advanced in a recent paper (Brandenburg 2020, hereafter B20), where it was shown that the resulting large-scale field, B LS , increases approximately linearly with time, typically by three orders of magnitude in the models of B20, while the thermal emission continues to decrease.
In this Letter, we investigate the implications of the scenario proposed in B20 on a simulated population of neutron stars and observables such as the braking index, characteristic age τ char , and the viability of this idea to explain the neutron star population. We begin in Section 2 discussing the motivation for a model where the magnetic field grows from a small-scale turbulent field at birth to a large-scale field. In Section 3, we describe the spin evolution of a neutron star under this hypothesis and perform simulations showcasing the evolution of a population of neutron stars for various observables in Section 4. We discuss implications of our results and conclude in Section 5.

GROWTH OF LARGE-SCALE FIELD
A proto-neutron star is typically fully convective at birth and remains so for tens of thousands of turnover times (Epstein 1979), corresponding to a time of O(1 min) (Burrows & Lattimer 1986). It is com-monly assumed that, within this time, the progenitor magnetic field is destroyed and then regenerated by a dynamo (Thompson & Duncan 1993). However, given the limited time between the collapse of the progenitor star and the end of the convective period, it is conceivable that only a small-scale magnetic field with scales comparable to the size of the turbulent eddies is generated. At the neutron star surface, the scale of such a field can be estimated by the density scale height H, which we expect to coincide roughly with the thickness of the crust of ≈ 5% of the neutron star radius (Cumming et al. 2004). This corresponds to a typical wavenumber k = 1/H, and thus to a typical spherical harmonics degree of ℓ = (H/R) −1 ≈ 20, where R is the radius of the neutron star. Note, however, that the thickness of the crust will depend on the star's mass, the details of the equation of state, and the scale height will be some fraction of this thickness. Therefore, we will explore a range 10 ≤ ℓ max ≤ 40.
Eventually, when the neutron star solidifies, the magnetic field in the crust is still randomly distributed in space. It then decays in a way analogous to a turbulent decay, but here the dynamics are governed by the Hall cascade (Goldreich & Reisenegger 1992). In this case, similar to the case of ordinary magnetohydrodynamic (MHD) turbulence, the magnetic field undergoes an inverse cascade, where not only the peak of the magnetic energy spectrum moves toward smaller k, but also the spectral energy at the smallest wavenumber increases in time. Numerical simulations in B20 found the large-scale field B LS to increase with decreasing rms magnetic field strength, B rms , like B LS ∝ B −5 rms ; see his Figure 12(b).
In the Hall cascade, B rms decays with time like B rms ∝ t −p/2 , where p = 2/5 is the decay exponent of the magnetic energy for a helical field (B20). Thus, so the large-scale magnetic field increases linearly with time; see Sarin et al. (2023) for supplemental material showing that the growth of B LS is actually closer to t 1.2 , although this minor departure from the expected linear growth will be ignored here. This linear growth applies to the case when the magnetic field is fully helical (see Runs E and F in Figure 12 of B20). However, it also applies to the case of nearly nonhelical fields (see Runs C and D in Figure 12 of B20) because the magnetic helicity is nearly unchanged, but the magnetic energy decays, so the ratio increases, making the magnetic field in the end nearly fully helical; see Tevzadze et al. (2012) for similar behavior in the MHD case. The increase of B LS continues until the largest scale in the system has been reached. Mathematically, the scale increases with t like ξ M ∝ t q , where q = 2/5 = 0.4 in the helical case, and 4/13 ≈ 0.31 in the nonhelical case (Brandenburg 2023). When ξ M becomes comparable to R, we must adopt global spherical geometry.
So far, inverse cascading has not yet been seen in global simulations. There are multiple reasons for this. First, it helps to initialize the simulations with a magnetic field having a broken power law energy spectrum with a proper inertial range and a sufficiently steep subinertial range. Second, large numerical resolution is required to obtain inverse cascading, at least in the nonhelical case (B20 used 1024 3 mesh points). The global simulations of Gourgouliatos et al. (2020), for example, did not have power laws, extended only to spherical harmonic degrees of eighty, and did not include magnetic helicity, which is the case considered here. Dehman et al. (2023) also presented global simulations with initially complex magnetic field configurations, but again not with power-law initial fields or with helicity.
In the absence of suitable global simulations, we argue here that we can substitute the spherical harmonics representation of the Laplacian, ℓ(ℓ + 1)/R 2 , by just k 2 , as in Cartesian geometry. Simplifying this further to k ≈ ℓ/R, and using k ≈ ξ −1 M ∝ t −q with q = 2/5 in the helical case, we have ℓ ≈ kR ∝ Rt −q . (2) As argued above, the initial value of ℓ is expected to be somewhere around ℓ max ≈ 20, and will then decrease to where τ 1 marks the start of the decay. Thus, we have ℓ = 1 at t ≡ τ 2 = τ 1 ℓ 1/q max . After that time, the inverse cascade stops, and then even the large-scale field can only decay. This decay is expected to be either exponential or in a power-law fashion. We assume here the latter and thus propose for the time evolution of the large-scale magnetic field, where B p,0 is the initial large-scale poloidal (dipolar ℓ = 1) component of the neutron star magnetic field.
Here, m dictates the slope of the large-scale field decay, which is not well constrained. Finally, let us note that the value of τ 1 is linked, in the simulations of B20, to the diffusive timescale t d ≈ ξ 2 M /η, where η is the magnetic resistivity. Low values of τ 1 generally correspond to a high resistivity, as we will discuss in the following sections.

SPIN EVOLUTION
With the behavior of the magnetic field described above, we can continue to describe the evolution of the neutron star spin. For the sake of simplicity, we start with the standard vacuum spin-down formula (Ostriker & Gunn 1969) ignoring other spin-down mechanisms such as gravitational-wave emission. Using Equation (4), our neutron star spin-evolution is then governed byΩ Here, c is the speed of light, and I is the moment of inertia of the neutron star. We emphasize that the above model is effectively the common vacuum-dipole spin-down equation (Ostriker & Gunn 1969) except for accounting of the evolution of the large-scale magnetic field following Eq. (4). Under the assumption that the magnetic field evolves from small scales to large scales in a fully helical manner, the timescales τ 1 and τ 2 can be related to each other via Here, ℓ max can be estimated as the inverse ratio of the electron density scale height to the radius of the neutron star, i.e., (H e /R) −1 . As discussed above, the magnetic field grows until ℓ → 1, i.e., until the wavenumber approaches a dipole from the initial small-scale wavenumber ℓ max . For later reference, values of ℓ max = 40, 20, and 10 correspond to τ 2 /τ 1 ≈ 10 4 , 1800, and 300. Owing to the linear scaling of B LS in Eq. (1), the ratios τ 2 /τ 1 also correspond to the amplification factor of the large-scale field.

SIMULATION
With the model for the spin and magnetic field evolution of neutron stars described above, we now consider the effect of such a model on several neutron star properties and observables through a series of simulations. For each of our simulations, we draw the initial spin-period P 0 and initial magnetic field B p,0 from astrophysically motivated distributions p (Igoshev et al. 2022;Pagliaro et al. 2023), in particular, p(ln P 0 /s) = N (−1.25, 0.99), p(log 10 B p,0 /G) = U(8, 11), where N and U denote normal and uniform distributions in given ranges.
We also fix m = 2.3 throughout this work, consistent with the decay of the magnetic field under ohmic dissipation (e.g., Pons & Geppert 2007). We note here that a precise value of m is not important to this work as we focus on the growth of the magnetic field rather than the decay. We have verified that the results are robust to our choices of m, the initial spin, and large-scale magnetic field. Our motivation for choosing these ranges were to show 1) that one does not need large scale magnetic fields at birth to make magnetars and 2) that a regular neutron star spin population can make magnetars.
We start by investigating the evolution with time; we set τ 1 = 10 2 yr to τ 2 = 10 6 yr, i.e., ℓ max ≈ 40, and evolve 100 neutron stars drawn with initial periods and magnetic fields from the distributions above forward in time. Figure 1 shows the evolution for different observables. In particular, we show the evolution of neutron stars across the PṖ diagram, the period and period derivative as functions of time, the observable braking index, dipole component of the magnetic field and the characteristic age. Figure 1 helps illustrate the model's main features. For example, large-scale initial fields as small as 10 11 G can be amplified to magnetar strengths (10 14 G, alleviating the concern of needing a large-scale poloidal field at birth. In particular, the age at such field strengths depends on the choice of τ 1 and ℓ max . Similarly, the model can naturally produce a range of braking indices despite being built only on the assumption of vacuum-dipole emission, and naturally returns values of the braking index n ≲ 2 for young pulsars, as measured in the observed population (Espinoza et al. 2017). We also show the evolution of the characteristic age with time. We find that at young ages, there is a huge discrepancy between the characteristic age and the actual age of the neutron star, while after some time these ages tend to agree. This is expected as our spin evolution is only different from n = 3 at early times. The evolution across the PṖ diagram also demonstrates that the model can accommodate neutron stars in various parts of the PṖ diagram despite a narrow starting position and explain the different locations we observe for real neutron stars.
To illustrate this last point, we now simulate a hypothetical PṖ diagram with τ 1 = 1 yr and ℓ max = 40 corresponding to τ 2 = 10 4 yr, ignoring selection effects and drawing ages randomly throughout the lifespan of the neutron star. In Figure 2, the left panel shows a real PṖ diagram of pulsars obtained from the ATNF catalog (Manchester et al. 2005) through the psrqpy software package (Pitkin 2018). The middle panel shows a 1000 simulated neutron stars with our model drawn ran- Figure 2. Real PṖ diagram of pulsars obtained from the ATNF catalog while the middle panel shows a simulated PṖ diagram with our model for a 1000 different neutron star initial conditions using τ1 = 10 2 yr and τ2 = 10 6 yr, corresponding to ℓmax = 40. Both panels include lines for constant magnetic fields and characteristic ages and the pulsar death line, calculated using the psrqpy software package. The color of the dots in the middle panel indicates the age of the selected neutron star, with neutron stars older than 10 9 yr in black, red for neutron stars between 10 5 -10 9 yr, blue for neutron stars between 10 2 -10 5 yr, and green for neutron stars younger than 10 2 yr. The last panel shows the histogram of the large-scale magnetic field of the 1000 neutron stars at their selected ages. domly from the above distributions and selected at different ages. Both panels include lines for constant magnetic fields and characteristic ages and the pulsar death line (Zhang et al. 2000), calculated through psrqpy with the region below shaded in yellow. The last panel shows the histogram of the large-scale magnetic field of the 1000 simulated neutron stars at their selected ages. We emphasize that we do not model processes such as mass transfer that is likely responsible for neutron stars observed in the lower left of the PṖ diagram, i.e., the millisecond pulsars. We caution the reader against a direct comparison of the observed population and our simulated population. We do not model the radio luminosity of any of the neutron stars and therefore cannot account for any selection effects that may limit our ability to see neutron stars in certain parts of the PṖ diagram (Faucher-Giguère & Kaspi 2006;Szary et al. 2014). Figure 2 also demonstrates that this evolutionary model for the magnetic field can produce neutron stars in the "magnetar" region of the PṖ diagram while also accounting for the distribution of neutron stars across the PṖ diagram for a small range of initial conditions. Furthermore, it helps illustrate the significant differences in neutron stars' actual ages and magnetic fields compared to the characteristic spin-down age and constant magnetic-field models. This may be important to reconcile the often significant discrepancies in inferred ages of neutron stars in supernova remnants with their location in the PṖ diagram.
The above PṖ diagram represented one single case of τ 1 and ℓ max . However, as we mentioned in Section 2, these values are unknown. In Figure 3, we show simu-lated PṖ diagrams for a range of these values. We note that each panel has the same initial conditions (i.e., P 0 and B p,0 ) and all neutron stars are selected at the same ages, and the difference in their "observed" locations is only due to different τ 1 and ℓ max . Figure 3 can be interpreted in two primary ways. First, for almost any choice of τ 1 and ℓ max , this model can accommodate a large portion of the PṖ diagram for a small range of initial conditions and provide a natural explanation for the often discrepant age estimates from spindown and a supernova remnant. Second, to explain magnetars with this model, we need large ℓ max and small τ 1 , i.e., we need the inverse cascades to begin early and need the initial field to be confined to small scales.

CONCLUSIONS
Our results demonstrate that it is indeed possible to explain the production of large-scale magnetic fields through an inverse cascade. This means that an initial small-scale magnetic field gets gradually converted into a large-scale one. The efficiency of this process increases with increasing dynamic range in space and time, i.e., it becomes more efficient, the smaller the scale of the initial field (larger ℓ max ), and the earlier the inverse cascade begins (smaller τ 1 ).
Our models for the growth of the magnetic field early on in the life of the star can explain the low (n ≲ 2) braking indices measured in young pulsars, and their observed evolution in the PṖ diagram (Espinoza et al. 2017). This is consistent with numerical results presented in Gourgouliatos & Cumming (2015), where the low braking indices of Vela and other young pulsars were interpreted as a byproduct of the increase in Figure 3. Simulated PṖ diagrams for the same initial conditions in all panels and ages but different τ1 and ℓmax. The colors correspond to the ages of the neutron stars, with the same ranges as described in Fig. 2. large-scale magnetic fields in Hall MHD. However, note that Gourgouliatos & Cumming (2015), started with initially strong large-scale magnetic fields, while we can reproduce magnetars with much smaller initial large-scale fields.
The most extreme case studied in this Letter is ℓ max = 40 and τ 1 = 1 yr. We estimated ℓ max = 20 based on the scale height, but this provided a rough estimate, because the convective eddies during the first minute of the neutron star's life can well be smaller than the local scale height as discussed in Section 2. Somewhat larger values of ℓ max = 40 appear therefore feasible. However, the situation with the value of τ 1 is less obvious. Once the crust has been established, the relevant timescale is the magnetic diffusion time, which was also used as the natural time unit in B20 and Brandenburg (2023). In their Cartesian models, it was defined as τ d = ξ 2 M /η. Using η = 4 × 10 −8 m 2 s −1 and ξ M = R/ℓ max ≈ 700 m with ℓ max = 20, we have τ d = 0.4 Myr. Somewhat smaller values are expected as we increase ℓ, for example to 100, in which case we can have ξ M = 150 m, and therefore τ d = 20 kyr. In addition, the time τ 1 when the inverse cascade begins is a certain fraction of τ d , for example τ 1 /τ d = 0.01, corresponds to the time when a certain characteristic quantity (the Hosking integral) reached a maximum; see Figure 4 of Brandenburg (2023). Thus, the smallest conceivable value of τ 1 is 200 yr for a standard crustal composition. It should be noted, however, that recent calculations of the conductivity in the presence of a nuclear pasta phase (Pelicer et al. 2023) have yielded significantly smaller conductivities, which could accommodate values of τ 1 ≈ 1 yr.
Our results show that low values of τ 1 (together with high values of ℓ max ) will create stars in the magnetar range in most of our models, and therefore point towards the presence of a thin, highly resistive pasta layer at the base of the crust, as suggested also by Pons et al. (2013), to explain the absence of long period isolated X-ray pulsars. We note that all models produce neutron stars below the death line, but those would not be observable as pulsars owing to their low luminosity.
Another implication of this model is on the ages of magnetars. If the magnetic field is confined to smallscales at birth and the initial dipole component is small, then there needs to be a minimum time such that the field can grow to magnetar strengths. For the choice of τ 1 = 10 yr and ℓ max = 40, this is ≥ 500 yr for our simulation, i.e., there should not exist any neutron stars with B p ≳ 10 14 G, younger than ≈ 500 yr. A comparison between characteristic age and true age is presented in Figure 1, where we see that at early times most of our systems are much younger than they appear (we caution, however, that the systems we produce with large characteristic ages have very small values ofṖ at birth, which would make them challenging to observe), while in the range between 10 3 and 10 5 years, a small number of systems, which correspond to the magnetars, are slightly older than they appear. However, we note that we do not capture the true diversity of initial conditions such as spin or considerations of the equation of state which could produce magnetar-like field strength earlier in time. We will explore the full effect of different initial conditions, equations of state, and compare directly to observations in future work.

ACKNOWLEDGMENTS
We are grateful to Clara Dehman, Cristobal Espinoza, and Kostas Gourgouliatos for helpful comments. We also thank the referee for their helpful feedback that has improved the presentation of this work. N.S. is supported by a Nordita fellowship. A.B. acknowledges support from the Swedish Research Council (Vetenskapsrådet, grant number 2019-04234). Nordita is supported in part by NordForsk. B.H. acknowledges support from the National Science Centre Poland (NCN) via grant OPUS 2018/29/B/ST9/02013.