A connection between proto-neutron-star Tayler–Spruit dynamos and low-field magnetars

Magnetars play a special role in modern high-energy astrophysics. They have been suggested as central engines for superluminous supernovae¹ and ultralong γ-ray bursts². They produce at least a fraction of the mysterious fast radio bursts^3,4. While Galactic magnetars are scarce due to their short life—with 30 known magnetars, compared with 3,500 radio pulsars according to Australia Telescope National Facility Pulsar Catalogue v.2.1.1⁵—it is estimated that around 10% of all neutron stars (NSs) are observed as magnetars at some point in their evolution⁶.

The standard magnetar model explains quiescent X-ray emission, spin period, bursts, outbursts and giant flares observed from anomalous X-ray pulsars and soft gamma repeaters (SGRs) by assuming that these NSs have strong dipolar magnetic fields, approximately 10¹⁴–10¹⁵ G (refs. ^7,8). However, a substantial fraction of magnetars (5 out of 30 known objects) in fact have dipolar magnetic fields well below 10¹⁴ G and have therefore been named low-field magnetars^9,10,11,12. It has been suggested that low-field magnetars are old NSs primarily powered by crust-confined toroidal magnetic fields with strength of approximately 10¹⁴ G (refs. ^10,13). Rea et al.¹¹ suggested that low-field magnetars are born with both poloidal and toroidal magnetic fields greater than 10¹⁴ G, but that the poloidal component decays by a factor of six in approximately 500 kyr. Phase-resolved X-ray observations show that in two cases low-field magnetars host small-scale magnetic fields which are 10–100 times stronger than their dipolar fields^14,15.

The origin of magnetar magnetic fields is a subject of debate¹⁶. Different dynamo mechanisms have been proposed to explain the formation of the strongest magnetic fields, including proto-NS convection^7,17,18,19, magnetorotational instability (MRI)²⁰ and more recently the Tayler–Spruit dynamo^21,22,23. The Tayler–Spruit dynamo is a particularly promising mechanism for generating magnetars’ magnetic fields in cases when the progenitor core is slowly rotating and the proto-NS is spun up by fallback accretion²². In cases of rotation periods larger than 5 ms, a normal core-collapse supernova is expected to occur, in agreement with observational constraints for the majority of magnetars^24,25. After the first minute, the proto-NS cools down, it shrinks in radius, crust solidification begins and the remnant becomes an NS. After this time, the initially complicated crustal magnetic field slowly relaxes due to Ohmic decay and Hall evolution on a timescale of 10⁵–10⁶ yr (ref. ²⁶).

Previous simulations of magnetothermal evolution have assumed idealized initial conditions (dipole, quadrupole)¹³ or random small-scale fields^27,28,29,30 rather than detailed magnetic configurations generated by a specific dynamo mechanism. However, the study of more realistic initial conditions is of key importance to obtain realistic predictions of magnetar properties. Indeed, Hall evolution has been shown to preserve certain aspects of the initial conditions^31,32. Hence, the observational properties of magnetars, and low-field magnetars in particular, should contain information about the proto-NS magnetic field.

The proto-NS dynamo and NS crust stages are modelled separately because of their very different timescales and physical conditions. While the dynamo is formulated as a magnetohydrodynamics (MHD) problem for a stably stratified fluid with shear caused by fallback accretion (which is mimicked through the boundary condition, see Methods for more details) over a timescale of a few tens of seconds, the magnetothermal evolution of the NS crust occurs on a much longer timescale of 1 Myr and is formulated as electron MHD (eMHD).

The initial condition for our NS simulation is a magnetic-field configuration taken from the late stages of model Ro0.75s from ref. ²³ (Table 1 in their supplemental materials) and corresponds to a Tayler–Spruit dynamo branch recently discovered in direct numerical simulations and characterized by a dipolar symmetry (that is, equatorially symmetric)²³. The radius of the proto-NS at this stage is 12 km, which is very similar to the NS radius. This magnetic field is obtained using the three-dimensional spherical MHD code MagIC³³ for rotation frequencies of the outer and inner spheres respectively Ω_o = 4Ω_i = 628 rad s⁻¹ (which corresponds to a surface rotation period of 10 ms; see Methods for a more detailed description).

In the transition from MHD to eMHD we preserve the angular structures up to ℓ = 30, where ℓ is the spherical harmonic degree, and downsample the radial part to take into account the differences in simulation set-up such as different geometric aspect ratios and presence of crust in eMHD simulations; full details can be found in Methods and Extended Data Figure 1. We also run a separate shorter simulation keeping angular structure up to ℓ ≤ 60, which evolves similarly to our basic one (see Extended Data Figure 2 for comparison).

The magnetic field is predominantly toroidal and reaches values up to 3 × 10¹⁵ G inside the volume, but the field at the outer boundary is much weaker. Assuming a scenario in which the core magnetic field is expelled to a crust-confined configuration, we extract the magnetic field in the top 10% of the simulation volume and adapt it to our code to model crust-confined NS magnetothermal evolution (Methods). Figure 1 shows the initial configuration of the magnetic field inside the NS crust.

We then use the PARODY code to integrate the coupled magnetic induction and thermal diffusion equations for 1 Myr before analysing the NS magnetic characteristics (Methods).

Figure 2 shows the dipolar and quadrupolar poloidal magnetic-field intensities, which are the only components that could contribute substantially to electromagnetic spin-down. The surface dipolar magnetic field increases by a factor of only three during the first million years, reaching a maximum value of 1.5 × 10¹² G, and the quadrupole component remains similarly small, with a maximum of around 6 × 10¹² G. These values are two to three orders of magnitude smaller than the internal magnetic-field strength in the crust.

Figure 3 shows a complex surface magnetic-field topology featuring individual arches elongated in the north–south direction. The local field strength at the footpoints of these arches on the NS surface reaches 10¹⁴ G, 100 times stronger than the surface dipolar magnetic field. Surface small-scale magnetic fields remain dominant at all times from the beginning of the evolution until 1 Myr later (Fig. 2). Our numerical simulation therefore successfully reproduces two crucial properties of low-field magnetars: (1) weak dipolar magnetic field and (2) presence of very strong (50–100 times stronger) small-scale magnetic fields, similar to those found in SGR 0418+5729¹⁴ and Swift J1882.3-1606¹⁵.

The X-ray observations of low-field magnetars are consistent with thermal emission from isolated hotspots with sizes of ≤1 km (ref. ³⁴) and black-body temperatures reaching T_bb = 0.12–0.6 keV. The bulk NS emission is not detected with typical upper limits less than 10³¹ erg s⁻¹. SGR 0418+5729 has a pulsed fraction of 62 ± 10% in the [0.3–1.2] keV range³⁴. CXOU J164710.2-455216 and Swift J1822.3-1606 have quiescent pulsed fractions 80 ± 3% and 38 ± 3% in the [0.5, 10] keV range respectively³⁵. The upper limit on the bulk thermal emission indicates that low-field magnetars are at least approximately 200 kyr old because the bulk X-ray emission drops below 10³¹ erg s⁻¹ after 200 kyr (ref. ³⁶) for strongly magnetized NSs (with internal field strengths of approximately 10¹⁵ G).

Strong magnetic fields can create large variations in surface temperature, as illustrated in Fig. 3 (see also Extended Data Figure 3). We see variations of an order of magnitude between the hottest (T ≈ 4.8 × 10⁵ K) and coldest (T ≈ 4.3 × 10⁴ K) regions. These variations could cause up to 20% pulsed fraction but would stay undetectable because of the low bulk X-ray luminosity, 10³¹ erg s⁻¹, and small effective T_bb = 0.028 keV.

We suggest that the observational properties of low-field magnetars can be explained by magnetospheric heating at the footpoints of small-scale magnetic arches visible in Fig. 3. Magnetospheric currents flow along the field lines and heat the surface, forming hotspots. The size of individual footpoints is a fraction of a kilometre, thus emission generated from these footpoints would have properties of emission seen from low-field magnetars, that is, very high temperature and small emission area.

Models of the magnetosphere generally assume the force-free condition J ∝ ∇ × B = μB (ref. ³⁷). This means that the strongest radial currents in the magnetosphere near the NS surface coincide with the strongest surface radial magnetic fields. These radial currents heat the surface, producing magnetospheric hotspots. Here we assume that only footpoints with radial magnetic field ∣B_r∣ > 7 × 10¹³ G are heated. Under this assumption, it is possible to form up to ten independent hotspots (Fig. 3), which, if heated to 3 × 10⁶ K, produce luminosity 2 × 10³² erg s⁻¹ and an emission area with radius ≈ 0.9 km. The lightcurve is sine-like with a pulsed fraction reaching 92% for a favourable orientation even without beaming, in agreement with X-ray observations of low-field magnetars (Extended Data Figs. 4 and 5; see also Table 1. Increasing the critical ∣B_r∣ leads to fewer hotspots with smaller areas, while decreasing the critical ∣B_r∣ results in a larger heated area. If the X-ray thermal emission is indeed generated close to the footpoints of these arches, the arches themselves provide natural sites where Compton scattering occurs and absorption features are formed¹⁴.

To assess whether this magnetic-field configuration can power the X-ray activity characteristic of magnetars, we examine the magnetic stresses inside the crust. Bursts and outbursts of magnetars are indeed thought to be caused by crust failure or plastic deformation due to the magnetic stresses^8,38,39. We apply the Lander–Gourgouliatos model³⁹ and compare the crustal magnetic stresses with the von Mises criterion for crust yielding (Methods). To obtain a conservative estimate, the crust is assumed to have completely relaxed only after 2 kyr. Figure 4 shows the average depth of crust failure regions developed at the age of 200 kyr. All the failing regions are located close to the original north and south magnetic poles, coinciding with the regions of strongest magnetic field generated by the proto-NS dynamo. The crust failure regions are much larger in the northern than in the southern hemisphere, due to the properties of the initial magnetic field. This is very different from earlier simulations with simple dipolar initial conditions⁴⁰, in which the crust failure occurred around the original magnetic equator.

The electromagnetic energy that can potentially be released in such a crustal failure is⁸

where l ≈ 1 km is the typical angular extent of the failing region and ∣B∣ ≈ 2 × 10¹⁴ G as extracted from simulations (Fig. 4). This value is actually well above the typical burst energy of approximately 10³⁷ erg of two low-field magnetars: SGR J0418+5729⁹ and CXOU J164710.2-455216⁴¹. Our modelling provides an upper limit on the extent of crust failure because it maps all the regions that could fail by a certain age.

By electromagnetic braking alone, NSs with dipolar fields of a few times 10¹² G cannot reach the spin periods of 8–11 s typical for low-field magnetars on a timescale of 1 Myr. However, it is essential to also take accretion into account, since the Tayler–Spruit dynamo can only develop if the proto-NS accretes fallback material, and this accretion will continue even after the NS is formed. Using the formalism of Ronchi et al.⁴² to model torques from the fallback disk, we naturally obtained periods of 8–11 s after 170 kyr for NSs with a dipolar magnetic field similar to that in our simulations (Fig. 5). More details of these calculations are summarized in Methods.

Most of the spin-down occurs during the propeller stage when the NS decelerates due to the interaction of its magnetosphere with the fallback disk (Extended Data Fig. 6). After 200 kyr, this propeller phase has spun the NS down to a rotation period of P = 8.5 s and a period derivative (dot{P}=8.5times 1{0}^{-13}) s s⁻¹. According to the standard magnetic dipole spin-down formula, the inferred surface magnetic dipole should then be B_dip ≈ 3.8 × 10¹³ G, which overestimates the true surface magnetic dipole in our simulation by a factor of approximately 40. This inferred value of B_dip is comparable to that of Swift J1822.3-1606 (1.4 × 10¹³ G) and below the upper limit measured for CXOU J164710.2-455216 (<6.6 × 10¹³ G) as well as 3XMM J185246.6+003317 (<4 × 10¹³ G).

The apparent magnetic field estimated using the instantaneous period and period derivative might be smaller if the disk is partially depleted and provides less torque. Depending on the exact amount of material left in the disk, the period derivative (dot{P}) could range from approximately 10⁻¹⁵ (electromagnetic spin-down only) to approximately 10⁻¹² (non-depleted disk; green area in Fig. 5). All low-field magnetars with a measured period derivative fall within this area.

Previous magnetothermal simulations have considered idealized, large-scale magnetic fields¹³. Some of these simulations can be made more similar to low-field magnetars by assuming a magnetar-strength dipolar magnetic field which is then dissipated by an increased crust resistivity¹¹. Moreover, these simulations remain highly axially symmetric because of the symmetries of the initial conditions. Although some previous studies have considered more complicated field structures, as expected from proto-NS evolution^29,30, our study directly implements the field from a self-consistent dynamo simulation. The crucial properties of our new magnetic-field configuration are that this field is predominantly toroidal and is initially localized near the polar regions of the crust. As a consequence, we find that crustal fractures are most likely to occur in these regions.

In comparison with the model suggested by Rea et al.¹¹, with an initial dipolar field as strong as 1.5 × 10¹⁴ G, the dipolar magnetic field stays very low, approximately 10¹² G, and does not decay substantially in our simulations.

The following mechanisms have been suggested to be responsible for formation of magnetar magnetic fields: (1) fossil fields^43,44, (2) convective dynamo^7,17, (3) dynamo due to MRI²⁰ and (4) Tayler–Spruit dynamo²³. There are also more exotic mechanisms such as the chiral magnetic instability in proto-NSs⁴⁵ leading to the formation of magnetic fields with sub-millimetre to centimetre scales.

Despite its simplicity, the fossil-field hypothesis is among the most uncertain ones since the magnetic-field configuration would be subject to multiple instabilities during the proto-NS stage. If the fossil field does survive the proto-NS stage, the strengths of poloidal and toroidal components should be roughly the same⁴³, which would point towards classical magnetars where the dipolar component is large and comparable to the strength of toroidal fields.

Raynaud et al.¹⁷ show that the convective dynamo can lead to classical magnetar formation if the proto-NS rotates very rapidly (spin period < 10 ms). Rapidly rotating proto-NSs are expected to be rare, and most proto-NSs should have longer spin periods. In the slow-rotation case, the convective dynamo works in the turbulent regime and would be expected to form mostly small-scale magnetic fields. It was suggested⁴⁶ that, due to an inverse Hall cascade, small-scale fields decay while large-scale fields grow. Other papers^27,47 were inspired by this turbulent dynamo and tried to test the inverse Hall cascade. The initial conditions were created by injecting energy for 10 ≤ ℓ ≤ 20. The outcome of these simulations was central compact objects because the average field was designed to be smaller than the one required to cause crustal yielding. Moreover, only a very limited inverse cascade was found. These papers^27,47 did not reproduce any actual dynamo simulations.

In a recent paper, Dehman et al.³⁰ performed a simulation of magnetothermal evolution starting with an initial magnetic configuration similar to one produced by MRI⁴⁸. They focused on energy spectra and did not try to precisely reproduce the spatial distribution of electric currents, similarly to refs. ^27,47. Note also that, unlike the Tayler–Spruit dynamo, the MRI develops when the proto-NS has a radius R ≈ 40 km, which is notably larger than that of a typical NS, which forced ref. ³⁰ to rather roughly mimic this contraction as no simulation has ever been made of the relaxation of MRI-generated magnetic fields. Similar to our research, they³⁰ found no substantial growth of the dipolar component, and most of the total magnetic energy was stored in small-scale magnetic fields. They interpreted observational properties of their objects as similar to those of central compact objects. It is unknown if their configuration leads to the formation of surface patches with strong radial magnetic field or to crust yielding. In their scenario, there is no need to consider a fallback disk as the rapid rotation originates from the progenitor; these NSs should thus have much shorter spin periods, which indeed makes them more similar to central compact objects.

The following aspects of the Tayler–Spruit dynamo scenario (considering proto-NSs with rotation periods longer than 6 ms) make it especially suitable for the formation of low-B magnetars: (1) formation of strong toroidal magnetic fields necessary for magnetar activity, (2) weak dipole and quadrupole at the surface, which naturally explain low-B magnetars, (3) presence of small-scale fields, which form patches of very strong radial magnetic field as a result of the Hall evolution, and (4) presence of a fallback disk, which allows the NS to gain a long spin period. Note that the Tayler–Spruit dynamo may also generate a stronger magnetic dipole and lead to the formation of a classical magnetar for rotation periods shorter than 6 ms (P.B., J.G., R.R. & A. Reboul-Salze, manuscript under review). Note that such a rotation period is close to the minimum initial proto-NS period derived from the kinetic energy of supernova remnants associated with magnetars²⁴.

Our results also suggest an important connection between low-field magnetars and recently discovered long-period radio pulsars, such as PSR J0901-4046⁴⁹. If the NS continues to operate in the propeller phase, it will ultimately reach periods comparable to 75 s by 10 Myr (Fig. 5). The external magnetic-field configuration remains complex, with large open field-line curvature near the NS surface facilitating radio pulsar operation. Thus, pulsar radio emission could occur if the disk is depleted.

Mahlmann et al.⁵⁰ performed numerical simulations for X-ray outbursts with energies up to 10⁴³ erg produced by a twisted magnetar magnetosphere. In our simulations, we see the development of individual magnetic arcs and the evolution of their footpoints. Thus, our results can be used as the initial magnetic field for future relativistic magnetosphere simulations. An additional open question is how magnetospheric currents flowing along small-scale arches will be sustained by strong-field quantum electrodynamic reactions and plasma dynamics in the magnetosphere.

Our work opens new perspectives for testing extreme dynamos operating in proto-NSs. We suggest that different dynamos leave their unique imprint on magnetic-field configurations, thus allowing us to identify different magnetic amplification processes using the magnetothermal properties of young isolated NSs. While we suggest that the formation of low-field magnetars is linked to the Tayler–Spruit dynamo, the formation of classical magnetars as well as the internal structure of their magnetic fields remains an open question.

Simulation of the proto-NS dynamo

We simulate a proto-NS with a mass of 1.4 M_⊙ and a radius R_NS = 12 km. Its interior is modelled as a stably stratified fluid enclosed between two spherical shells. To control the differential rotation, we impose constant rotation frequencies on both shells (spherical Taylor–Couette configuration), with the outer shell rotating faster than the inner shell to be consistent with the fallback formation scenario. We solve the Boussinesq MHD equations by using the pseudospectral code MagIC^51,52. In this code, the different lengths r, the time t, the temperature T and the magnetic field B are scaled as follows:

with the gap between the two spheres d = r_o − r_i = 9 km, the kinematic viscosity ν = 3.5 × 10⁹ cm² s⁻¹, the temperatures of the outer and inner spheres T_o and T_i respectively, the constant density ρ = 4.1 × 10¹⁴ g cm⁻³, the resistivity η = 3.5 × 10⁹ cm² s⁻¹ and Ω_o = 628 rad s⁻¹. So, the dimensionless equations solved using MagIC read

where v and B are the velocity and magnetic field. The potential Π includes all gradient forces and functions as a Lagrange multiplier to ensure the solenoidality of the velocity field (equation (4)). D/Dt ≡ ∂/∂t + v·∇ is the Lagrangian derivative. E, Ra, Pr and Pm are dimensionless numbers, which depend on the fluid properties. The Ekman number E is defined as the ratio of the rotation period to the viscous timescale,

The thermal and magnetic Prandtl numbers are defined by

where κ = 3.5 × 10¹⁰ cm² s⁻¹ is the thermal diffusivity. Finally, the Rayleigh number Ra measures the ratio between the timescales of thermal transport by diffusion to the thermal transport by convection,

where

is the Brunt–Väisälä frequency. The gravitational acceleration is assumed to be purely radial, g = g₀r/r₀e_r. Y_e and S are the electron fraction and the entropy, respectively.

We apply no-slip, electrically insulating and fixed-temperature boundary conditions on both shells. The spin-up torque exerted at the outer boundary by the choice of a fixed angular velocity mimics in an approximate way the effect of fallback accretion. The resolution used is (n_r, n_θ, n_ϕ) = (257, 256, 512), where n_r, n_θ and n_ϕ are numbers of grid points in the radial, latitudinal and longitudinal direction. For more information on the numerical methods, see the supplementary materials of ref. ²³.

Conversion between MagIC and PARODY codes

There are substantial differences between MagIC and PARODY codes. In particular, in MagIC the radial part of the potentials is expanded using Chebyshev polynomials, while in PARODY the radial part is represented using finite differences. While PARODY covers only the NS crust, the MagIC simulations include a core. Therefore, we downsample the result of MagIC simulations in the radial direction.

The poloidal–toroidal decompositions and thus the magnetic potentials are defined differently in the MagIC and PARODY codes. Specifically,

where the superscript M/P refers to MagIC/PARODY, respectively.

Moreover, the codes use different normalization factors C_lm for the spherical harmonics ({Y}_{l}^{,m}(theta ,phi )). The spherical harmonics are normalized in the PARODY code as

while the normalization in the MagIC code reads

where δ_m,0 is the Kronecker delta.

Thus, for the radial magnetic field, we have

Making this comparison for each (l, m) separately we thus obtain

Expanding and simplifying this expression we obtain two different equations for axisymmetric and non-axisymmetric poloidal potentials,

Similarly, we can proceed with the θ component of the magnetic field computed using only the toroidal potential

Thus, the normalization is

which simplifies to

In this work, we preserve the angular structure obtained in dynamo simulations at the surface and in the middle of the crust up to ℓ = 30, which corresponds to surface structures of approximately 1 km. In addition, we ran shorter simulations where we preserve the angular structure up to ℓ = 30. We compare the results in Extended Data Figure 1 and show that they are very similar to our basic simulation.

Crust-confined magnetic-field configurations

In addition to the technical details in the previous section, the proto-NS dynamo set-up and the magnetothermal crust evolution set-up differ in their geometry, having geometric aspect ratios χ_pNS = 0.25 and χ_NS = 0.9, respectively. Thus, to create a magnetic-field configuration that is similar to proto-NS results but is also crust confined, we should extract only the top 10% of the proto-NS simulation.

Our approach for importing the results of the dynamo simulations is to require all components of the magnetic field to coincide between MHD and eMHD simulations at certain points within the crust. We consider the poloidal and toroidal potentials for each individual spherical harmonic, and require both these potentials to exactly coincide with our numerical fits at the following points: r₁ = 0.93 and r₂ = 0.96. We require our fit for the poloidal potential to coincide at the surface. We also require our poloidal and toroidal potentials to satisfy the potential boundary condition at the surface and the ‘no-currents’ boundary condition at the core–crust interface.

Similarly to recent work⁵³, we represent the radial part of the poloidal and toroidal potentials as a polynomial expansion

Overall, all conditions for the radial part of the poloidal potential can be written as

Here β_p(r) are coefficients of the spectral expansion for poloidal magnetic field extracted from the proto-NS MagIC simulations. These conditions are individually satisfied for each l and m, and translate into the following system of linear equations:

For the toroidal potential we use the following conditions:

Similarly, β_t(r) here are the coefficients of the spectral expansion for the toroidal magnetic field extracted from the proto-NS simulations. These conditions then translate into the linear system

Simulation of NS magnetothermal evolution

Code PARODY^54,55,56 was modified to solve the magnetothermal evolution of NSs. The test cases (benchmarks) for the code are available in the literature⁵⁷. The pseudospectral code PARODY was modified to solve the following system of dimensionless partial differential equations for B and T:

where the first equation is the magnetic induction equation and the second is the thermal diffusion equation. The terms on the right-hand side of the first equation correspond to the Hall effect, Ohmic decay and the Biermann battery effect. The terms on the right-hand side of the second equation correspond to anisotropic thermal diffusion, Ohmic heating and entropy carried by electrons. The derivation of the above equations is summarized in ref. ⁵⁷. The same code was also used to compute the evolution of off-centred dipole configurations⁵³.

To ensure the solenoidality of B, we write the magnetic field as a sum of poloidal and toroidal parts,

The scalar potentials b_pol and b_tor are expanded in spherical harmonics.

The electron chemical potential varies within the crust as

The electrical resistivity varies within the crust, but is not sensitive to temperature in our simulations. Our chosen parameters roughly mimic resistivity caused by impurity with parameter Q_imp ≈ 40. In a real NS the resistivity is larger during the first 10 kyr, which means that the role of Hall evolution is slightly smaller during this time. In our simulations, we are interested in longer timescales where the Hall evolution dominates. The tensor (hat{chi }) describing the anisotropy of the heat transport is written as

where δ_ij is the Kronecker symbol and ϵ_ijk is the Levi-Civita symbol.

The dimensionless Hall (Ha), Seebeck (Se), Péclet (Pe) and Roberts (Ro) parameters depend on the chosen scales for the magnetic field and temperature, which we take to be B₀ = 10¹⁴ G and T₀ = 1.0 × 10⁸ K. The Hall number is defined by

where e is the electron charge, c is the speed of light, τ₀ = 9.9 × 10¹⁹ s is the electron scattering relaxation time⁵⁸ and μ₀ = 2.9 × 10⁻⁵ erg is the electron chemical potential at the top of the crust. The Seebeck number is defined by

where k_B is the Boltzmann constant and n₀ = 2.603 × 10³⁴ cm⁻³ is the electron number density at the top of the crust. Finally, the Péclet and Roberts numbers are

and

We model the core as a perfect conductor, which implies the following inner boundary conditions at r = 0.9:

We model the region outside the NS as a vacuum, which implies the following outer boundary conditions at r = 1:

where ({b}_{,text{pol},}^{lm}) is the coefficient of degree l and order m in the spherical harmonic expansion of the poloidal potential b_pol.

The temperature is fixed to its initial value at the core–crust boundary (see more details about modelling cooling at the end of the section). The outer boundary condition for the temperature is

where the (dimensional) surface temperature T_s is related to the crustal temperature T_b as

using a simplified relation⁵⁹.

The numerical resolution is n_r = 96 grid points in the radial direction and spherical harmonic degrees up to l_max = 128.

To take into account the NS cooling, we restart calculations at 200 kyr, changing the core temperature to 10⁶ K. We run calculations for 1 kyr to allow the simulation to relax, that is, crust temperatures stop evolving on short timescales, creating a stable surface thermal pattern. The timescale of 1 kyr is estimated numerically and is related to assumptions of how we model the thermal capacity of the crust.

Properties of thermal emission

We use the open-source code Magpies to model X-ray thermal lightcurves. We show these results in Extended Data Figure 4. The maximum pulsed fraction reaches 93% for the most favourable orientation of the rotational axis with respect to the original dipole axis.

Similarly to ref. ¹³, we try to fit the soft-X-ray lightcurve in the range 0.3–2 keV. We show the results in Extended Data Figure 5. We summarize the obliquity angle as well as inclination angles in Table 1. While SGR 0418+5729 and Swift J1822.3-1606 are fitted relatively well, the two remaining magnetars have more features in the lightcurves.

Crust failure

We use here a model developed in ref. ³⁹ based on earlier work in ref. ³⁸. Essentially, we use the von Mises criterion for crust yielding, following equation (14) of ref. ³⁹:

Here B₀ is the relaxed (initial) state of the magnetic field, which we assume to coincide with our first simulation snapshot at 2 kyr. τ_el is the scalar yield stress. B is computed at 200 kyr. We compute the critical strain following the procedure in ref. ³⁹ with a correction (S. Lander, personal communication)

where ({mathcal{R}}) is computed as

R_cc = 0.9 is the location of the crust–core interface and R_nd = 1 is the location of the neutron drip point. Thus, our critical strain varies from approximately 8 × 10²⁶ g cm⁻¹ s⁻² close to the neutron drip boundary to 4.6 × 10²⁹ g cm⁻¹ s⁻² at the core–crust boundary. Following our normalization, the stress caused by Lorentz forces (right-hand side of equation (42)) is multiplied by a numerical factor, (10¹⁴ G)². This von Mises criterion is written assuming that failure occurs in the form of shearing motion³⁹.

Accretion-driven spin-down

To explain the NS spin-down to the regime of low-field magnetars, we invoke the propeller mechanism⁶⁰ due to the interaction between the NS magnetic field and the remaining fallback disk. The evolution of the NS fallback depends on the three different radii: (1) the light-cylinder radius, (2) the magnetospheric radius and (3) the corotation radius, which are defined by the respective expressions

M_NS and Ω_NS are the NS mass and rotation rate; μ = B_NSR_NS³ is its magnetic dipole moment. (dot{M}) is the accretion rate. Strictly speaking, (dot{M}) is the material loss rate from the accretion disk. In the propeller regime this quantity remains positive even though the material is not accreted onto the NS.

If the disk penetrates the magnetosphere (r_lc > r_mag), it can either spin up the NS by accreting matter if r_cor > r_mag, or spin down the NS in a propeller phase if r_cor < r_mag. In this propeller phase, the magnetic field accelerates the inner disk to super-Keplerian speeds, which produces a centrifugal outflow. Angular momentum is therefore transported from the NS toward the disk, which can efficiently spin down the NS.

The modelling of the NS fallback evolution we use is strongly inspired by ref. ⁶¹ except for the mass accretion rate, which reads⁴²

where t_ν ≈ 30 s is the viscous timescale and ({dot{M}}_{0}={M}_{{rm{d}},0}/{t}_{nu }approx 6.5times 1{0}^{29}) g s⁻¹ is the initial accretion rate, and M_d,0 = 0.01 M_⊙ is the initial fallback disk mass. The torques exerted on the NS by the accretion disk are given by

where ({varOmega }_{{rm{K}}}=sqrt{G{M}_{{rm{NS}}}/{R}_{{rm{NS}}}^{3}}) is the Keplerian angular velocity. The dipole spins the NS down as follows:

Therefore, the NS angular velocity evolves as

where I_NS = 1.45 × 10⁴⁵ g cm² is the NS moment of inertia. Extended Data Figure 6 shows the time series of the characteristic radii and NS rotation period that result from the solution of equation (51) for B_NS = 10¹² G, M_d,0 = 0.01 M_⊙ and an initial rotation period of 10 ms. We clearly find that the NS is strongly spun down during the propeller phase and reaches the period range of the observed low-field magnetars at approximately 170 kyr. This timescale varies up to approximately 550 kyr for B_NS = 5 × 10¹¹ G.