Magnetic field evolution in magnetar crusts through three dimensional simulations

Current models of magnetars require extremely strong magnetic fields to explain their observed quiescent and bursting emission, implying that the field strength within the star's outer crust is orders of magnitude larger than the dipole component inferred from spin-down measurements. This presents a serious challenge to theories of magnetic field generation in a proto-neutron star. Here, we present detailed modelling of the evolution of the magnetic field in the crust of a neutron star through 3-D simulations. We find that, in the plausible scenario of equipartition of energy between global-scale poloidal and toroidal magnetic components, magnetic instabilities transfer energy to non-axisymmetric, kilometre-sized magnetic features, in which the local field strength can greatly exceed that of the global-scale field. These intense small-scale magnetic features can induce high energy bursts through local crust yielding, and the localised enhancement of Ohmic heating can power the star's persistent emission. Thus, the observed diversity in magnetar behaviour can be explained with mixed poloidal-toroidal fields of comparable energies.

Current models of magnetars require extremely strong magnetic fields to explain their observed quiescent and bursting emission, implying that the field strength within the star's outer crust is orders of magnitude larger than the dipole component inferred from spin-down measurements.This presents a serious challenge to theories of magnetic field generation in a proto-neutron star.Here, we present detailed modelling of the evolution of the magnetic field in the crust of a neutron star through 3-D simulations.We find that, in the plausible scenario of equipartition of energy between global-scale poloidal and toroidal magnetic components, magnetic instabilities transfer energy to non-axisymmetric, kilometre-sized magnetic features, in which the local field strength can greatly exceed that of the global-scale field.These intense small-scale magnetic features can induce high energy bursts through local crust yielding, and the localised enhancement of Ohmic heating can power the star's persistent emission.Thus, the observed diversity in magnetar behaviour can be explained with mixed poloidal-toroidal fields of comparable energies.A n estimate of the magnetic field intensity in a neutron star can be obtained by assuming that the observed spin-down is caused by electromagnetic radiation from a dipolar magnetic field [1].Neutron stars for which this estimate exceeds the QED magnetic field 4.4×10 13 G are conventionally called magnetars, and typically exhibit highly energetic behaviour, as in the case of Anomalous X-ray Pulsars and Soft γ-ray Repeaters [2,3].Puzzlingly, not all high magnetic field neutron stars exhibit energetic behaviour [4], and conversely, "magnetar-like" activity has been observed in pulsars for which this field estimate is below [5] or only marginally above the QED magnetic field [6,7,8].Furthermore, despite the high thermal conductivity of neutron stars' solid outer crusts [9], observations of their thermal emission indicates that in some cases the surface is highly anisothermal [10], with kilometre-sized "hot spots" thought to be produced by small-scale magnetic features [11].Phase resolved spectroscopy has revealed that some magnetars have small-scale magnetic fields whose strength exceeds their large-scale component by at least an order of magnitude [12,13], which can be correlated with outbursting events [14].These observations all imply that the magnetic field structure in magnetars is more complicated, and varied, than the traditional picture of a simple inclined dipole.
The origin of the extreme magnetic fields in these objects, which are the strongest found in nature, is uncertain.Even if magnetic flux were exactly conserved during the star's formation, the resulting field would not exceed 10 13 G.It seems likely then that the strong fields in magnetars must result from a combination of differential rotation and dynamo action prior to the formation of the crust [15].Dynamo models generally predict magnetic fields with features over a wide range of scales, and poloidal and toroidal components of comparable strength [16].However, once the crust forms, small-scale surface features in the magnetic field would decay by Ohmic dissipation much faster than the global field.
The evolution of the crustal magnetic field is mediated by the Hall effect, corresponding to advection by free electrons [17].The Hall effect in magnetars typically operates on a shorter timescale than Ohmic dissipation, and might therefore explain the formation of small-scale magnetic features, whose dissipation could then power the star's ther-mal radiation.However, axisymmetric simulations of the crustal magnetic field only generate such features if a toroidal magnetic field exceeding 10 16 G is assumed to reside within the crust [18,19].The toroidal field strength must significantly exceed that of the poloidal dipole, which is the only component that is measured from spin-down observations, presenting additional challenges to models of magnetic field generation.A possible resolution is the growth of localised patches of strong magnetic field via magnetic instabilities.Such instabilities have been demonstrated in local, plane-parallel models [20,21,22], but not yet in a realistic global model.

Hall Evolution in Neutron Star Crusts
Motivated by this puzzle we study the fully non-linear 3-D problem of the magnetic field evolution in a magnetar crust through numerical simulations, using a modified version [23] of the PARODY code [24,25], to determine the necessary conditions for the spontaneous generation of strong localised magnetic field out of a largescale weaker one.The crust of the neutron star is treated as a solid Coulomb lattice, in which free electrons carry the electric current and the magnetic field evolution is described by Hall-MHD with subdominant Ohmic dissipation [17].This process is described by the Hall-MHD induction equation: where B is the magnetic induction, n e the electron number density, σ the electric conductivity, c the speed of light and e the elementary

Significance
The observed diversity of magnetars indicates that their magnetic topology is more complicated than a simple dipole.Current models of their radiative emission, based on axially symmetric simulations, require the presence of a concealed toroidal magnetic field having up to 100 times more energy than the observed dipole component, but the physical origin of such field is unclear.Our fully 3-D simulations of the crustal magnetic field demonstrate that magnetic instabilities operate under a range of plausible conditions, and generate small-scale field structures that are an order of magnitude stronger than the large-scale field.The Maxwell stress and Ohmic heating from these structures can explain magnetar bursts and surface hotspots, using comparable poloidal and toroidal magnetic fields.
Reserved for Publication Footnotes   3).The snapshot is at t = 15kyr overplotted with the magnetic energy density e m = B 2 /(8π) on part of the surface.The surface field is highly anisotropic, with small regions in which the magnetic energy density exceeds by at least an order of magnitude the average surface value.
charge.The first term in the right-hand-side of equation ( 1) describes the Hall effect and the second one Ohmic dissipation.We assume a neutron star radius R * = 10km, and a crust thickness of 1km.We shall express the quantities in terms of the normalised radial distance r = R/R * , where R is the distance from the centre.The electron number density is taken to be n e = 2.5 × 10 34 cm −3 1. .Such profiles are good analytical fits of the form σ ∝ n 2/3 e of more precise crust models [26] at temperatures ≈ 10 8 K.
The code uses spherical harmonic expansions in latitude and longitude, and a discrete grid in radius.The linear Ohmic terms are evaluated using a Crank-Nicolson scheme, while for the non-linear Hall terms an Adams-Bashforth scheme is used.We have used the Meissner superconducting condition assuming that no magnetic field penetrates the crust-core surface, serving as the inner boundary of the domain [27].This is a simplifying assumption, as the core may be a Type-II superconductor and magnetically coupled to the crust.Nevertheless, the exact structure of the field in the core is still uncertain with some works suggesting that the crustal magnetic field is practically disconnected from the core [c.f.Figs. 3 and 4 of [28], and Fig. 9 of [29]].With respect to the outer boundary we matched the magnetic field to an external vacuum field satisfying ∇ × B = 0. Given the lack of any other 3-D code or analytical solution of the Hall evolution, the code has been benchmarked against axially symmetric calculations [30,31] in various radial and angular resolutions with excellent quantitative agreement, and the free decay of Ohmic eigenmodes, in the absence of the Hall term, reproducing the analytical results.The results that we present here have a resolution of 128 grid-points in radius, and spherical harmonics Y m up to max , m max = 60 − 80.We have repeated some of the runs with the strongest magnetic fields, using twice the angular resolution, with max = m max = 120.We found good agreement in the overall magnetic field evolution of those runs with the lower resolution ones.

Simulations: Initial Conditions and Results
Given the complexity of the processes taking place during the formation of a neutron star, our knowledge of its initial magnetic field is limited.We have therefore explored a wide range of initial conditions, based on three plausible scenarios regarding the magnetic field formation: the fossil field, the large-scale dynamo, and the small-scale dynamo.In the scenario of a fossil field [32], the magnetic field has relaxed dynamically to an axially symmetric twisted torus with comparable energy in its poloidal and toroidal components.The poloidal magnetic field has dipolar geometry and the toroidal magnetic field is strongest near the dipole equator; both components have = 1 symmetry, where is the spherical harmonic degree.The second scenario corresponds to a magnetic field formed through the combination of differential rotation and dynamo action in the proto-neutron star [15,16]; in this case the most likely geometry is an = 2 toroidal and a dipolar poloidal field.In these two scenarios the field is predominantly large-scale and axially symmetric.In the third scenario, turbulent convection in the proto-neutron star generates a small-scale dynamo magnetic field [33].In this case, the field comprises small loops with size comparable to that of the convection cells.
We have performed 76 simulations integrating the Hall-MHD induction equation (1) in each of the three scenarios for the initial magnetic field in the star's crust (see S.I.Tables 1-3).In the fossil field and large-scale dynamo scenarios we have used mixed poloidal and toroidal fields, varying the ratio of the toroidal to total magnetic energy and the overall normalisation of the magnetic field intensity.The intensity of the poloidal field on the surface takes values between 10 13 G to 5 × 10 14 G, except for the purely toroidal initial conditions.
The magnetic field is expressed in terms of two functions V p (r, θ, φ) and V t (r, θ, φ): We express the fraction of the energy in the toroidal component of the magnetic field e t , allowing the following values: 0 (purely poloidal), 0.1, 0.5, 0.9 and 1 (purely toroidal).B 0 governs the overall amplitude and takes the values 0.5, 1, 2 and 4 expressed in units of 10 14 G; the normalisation is chosen so that the initial magnetic field energy for a simulation with amplitude B 0 is 2.25 × 10 46 B 2 0 erg.For the radial structure of the field, we have used two basic profiles: one where the latitudinal and azimuthal components decrease monotonically with radius; and a second one with a local maximum in the middle of the crust.We present them below, providing the exact expressions to allow reproducibility of our results.
The first family of potentials corresponding to dipole fields ( = 1), is given by: r 2 cos θ(−913.2837868+ 16749.59757r 3 −137700.2916r 5+ 292459.7609r6 − 269535.2857r7 +120217.0825r8 − 21277.28441r 9) , [ 3 ] The poloidal part of this potential gives a dipole magnetic field supported by a uniformly rotating electron fluid corresponding to a Hall equilibrium [34], with the magnetic field strength at the pole being equal to B 0 for e t = 0.Under this profile B θ and B φ change monotonically with radius.The maximum value of the B φ component under this profile occurs at the base of the crust on the equatorial plane, and is B φ,max = 15B 0 e 1/2 t .Thus, for the strongest toroidal field simulated (B 0 = 4, e t = 1), the maximum azimuthal field is 6 × 10 15 G.
The second family of = 1 profiles is the following: r 2 cos θ(734.5987631− 2333.649604r+2465.151852r 2 − 865.6887776r 3 ) , [ 5 ] i i "pnas˙arxiv" -2018/2/12 -18:16 -page 3 -#3 These profiles have B θ and B φ components with an extremum close to the centre of the crust.The choice of normalisation is such that the magnetic energy in the crust is the same as the corresponding V ps and V ts profiles.The maximum value of the B φ component under this profile occurs in the middle of the crust (r = 0.95) on the equatorial plane and is B φ,max = 12B 0 e 1/2 t .Thus, for the strongest toroidal field simulated (B 0 = 4, e t = 1), the azimuthal field has a maximum value of 4.8 × 10 15 G.
We also combined the poloidal dipole fields with = 2 toroidal fields.In such cases we have utilised the following profiles: which is used in conjunction with V ps , and which is used in conjunction with V pu .In the simulations with the quadrupolar toroidal field we assumed equipartition in energy between the poloidal and the toroidal field (e t = 0.5) and varied B 0 = 0.5, 1, 2, 4 for the case of the positive sign, while we run simulations with B 0 = 2 only for the case of the negative sign.The maximum value of the azimuthal field assumed, using the V ts,q profile is 5×10 15 G and occurs at the base of the crust at latitudes θ = π/4, 3π/4; while for the case of of V tu,q is 3.8 × 10 15 G and occurs at r = 0.95 and θ = π/4, 3π/4.In all these axially symmetric initial conditions, we have superimposed a non-axisymmetric random noise, containing 10 −4 − 10 −8 of the total energy.Axially symmetric initial conditions correspond to equilibria with respect to the φ coordinate, so any growth of the non-axisymmetric part of the field may indicate the presence of a non-axisymmetric instability.Finally, we have simulated a highly non-axisymmetric magnetic field.In this case we have populated the crust with a fully confined non-axisymmetric magnetic field, where ∼ 0.9 of the total magnetic energy is in the non-axisymmetric part of the field, consisting of multipoles with 10 ≤ ≤ 40.Similarly to the previous cases we have used four normalisation levels for the overall magnetic field strength.
We ran these simulations until the magnetic energy decayed to 0.01 of its initial value.We also present a set of runs where the Hall term is switched off and the field evolves only via Ohmic decay.
Our basic conclusions are the following: 1.In systems where the initial azimuthal field is strong (e t 0.5), and especially when the profiles V pu and V tu are used (Figs. 1 and  2), the amount of energy in the non-axisymmetric part of the field increases exponentially.The growth rate is approximately proportional to the strength of the magnetic field (see S.I. Figure 1), i.e. for the case depicted in Figs. 1 and 2 the non-axisymmetric energy growth timescale is ∼ 1.5kyr.The exponential growth stops if one of the following two conditions is fulfilled: either 0.5 − 0.9 of the total magnetic energy has decayed, or the non-axisymmetric field contains about 50% of the total energy.While initially the m = 0 mode was dominant, once the instability develops energy is transferred and a local maximum appears around 10 m 20 (see S.I. Figure 3), with some further local maxima at higher harmonics, leading to a characteristic wavelength for the instability of about 5km at the equator, with finer structure appearing later giving kilometre-sized features.These features are well above the resolution of the simulation and appear both in the low and high resolution runs.This behaviour is indicative of the density-shear instability in Hall-MHD, where, based on analytical arguments [21,22], the wavelength of the dominant mode is 2πL where L is the scale height of the electron density and the magnetic field, which can be approximated by the thickness of the crust.
2. At the start of each simulation, the structure of the magnetic field changes rapidly, on the Hall timescale, unless the field is initiated in a Hall equilibrium in which case the evolution occurs on the slower Ohmic timescale.This change is accompanied by rapid magnetic field decay, i.e. in the case of QU05-4 (S.I.Table 3) 50% of the energy decays within the first 20kyr.We find that the early decay approximately scales with the magnetic field strength (Tables 1-3 S.I. column t 0.5 ), whereas the later decay has a much weaker dependence on the magnetic field strength.3. Initial conditions with poloidal magnetic fields generate toroidal fields.The energy content of the toroidal field never becomes dominant irrespective of the strength of the initial field, (see S.I. Figure 2).Even if a toroidal component exists in the initial conditions, stronger fields tend to suppress it more efficiently than weaker ones, in favour of the poloidal component.Thus a strong toroidal field is unlikely to be generated by a poloidal field being out of equilibrium.Once the magnetic field has decayed substantially the amount of the energy in the toroidal field may become higher than that of the poloidal.However, at this stage the dominant effect is Ohmic decay rather than Hall drift.This happens after a few Myrs as it is illustrated in the last column of the tables in the S.I. for the simulations with B 0 = 1. 4. The evolution of the axisymmetric part of the magnetic field is in accordance with the results of previous axisymmetric numerical and semi-analytical studies.In particular the toroidal field interacts with the poloidal field winding it up, while magnetic flux is transferred in the meridional direction leading to the formation of zones, Fig. 2. The drift towards the equator is due to the polarity of the toroidal field, an effect that has been discussed in detail in plane-parallel geometry in [35] and the axially symmetric case in [27]. 5. Configurations initially dominated by the poloidal field, tend to remain axially symmetric, especially if the V ps profile is used.Their evolution is in general accordance with previous axially symmetric simulations that have described the interaction of the poloidal and toroidal part of the field.6.Once about 90% of the magnetic energy has decayed the Hall effect saturates, and thereafter the magnetic energy decays on the Ohmic timescale.However, the remaining small-scale features decay at the same rate as the large-scale field, indicating that the Hall effect is still active.These results are indicative of the "Hall attractor" seen in axisymmetric simulations [30].7. Initial conditions corresponding to highly non-axisymmetric structures undergo a weak inverse cascade (see runs Turb in the S.I.).The Hall effect transfers energy across the spectrum both to smaller and larger scales.Despite the great differences in the structure of the magnetic field compared to the axially symmetric case, the evolution follows a qualitatively similar pattern, with the rapid decay at early times, and saturation of the Hall effect later.

Discussion
In all simulations we find that the magnetic field undergoes a major restructuring at the very beginning of its evolution.This process occurs on the Hall timescale, and is thus faster for stronger magnetic fields.We then identified two distinct behaviours.Initial conditions where the poloidal field was dominant remain axially symmetric, while cases with substantial, but not necessarily dominant, azimuthal fields triggered non-axisymmetric instabilities.As magnetic energy is conserved in Hall-MHD, this growth of the non-axisymmetric field occurs at the expense of the axisymmetric components.
In the model with initial dipole poloidal magnetic field of 4×10 14 G and a toroidal component containing the same amount of energy with the poloidal (Figs. 1 and 2), magnetic spots of sizes ∼ 2km appear.These spots contain energies up to 5 × 10 43 erg, while the dipole field at that time is 1 − 2 × 10 14 G, (Fig. 3).The rapid decrease of the dipole component, with a drop by a factor of 2 within 10-20 kyrs, is caused for the simulation shown in Figure 1.
by the latitudinal Hall drift towards the equator.This is more efficient when the stronger fields are in shallow depths, when the potentials V tu and V tu,q are used (runs DU05-4 and QU05-4), where the maximum value of the toroidal field occurs at r = 0.95.On the contrary it takes ten times longer for the same decrease in the dipole component, in the case of the potentials V ts (i.e.run DS05-4), where the maximum of the toroidal field occurs at the base of the crust.Therefore, in the case of shallower fields the dipole field drops quickly, while there is still a large reservoir of energy deeper in the crust.This can be related to AXP 1E2259+586 and SGR 0418+5729 which despite their relative weak dipole field exhibit magnetar activity [8,5].The energy in the magnetic spots is sufficient to power strong individual magnetar bursts or sequences of weaker events [36,37].The energy in the equatorial zone of width ∼ 2km is ∼ 2 × 10 46 erg, representing 20% of the total magnetic energy in the crust at that moment.Such energy is comparable to the energy released by a magnetar giant flare [38,39,40].Note though, that magnetars that have exhibited giant flares have spin-down dipole fields above 5 × 10 14 G, so the overall magnetic energy should be scaled by a factor of ∼ 10 compared to our simulations.These structures take ∼ 10 4 yrs to develop, given the choice of an initial condition where only 10 −4 of the total is in the non-axisymmetric part.Had we chosen initial conditions with more energy in the non-axisymmetric part, the development of the smaller scale structure would have been imminent and in agreement with the ages of the most active magnetars, which are in the range of a few 10 3 yrs.This was confirmed in the extreme case where no energy is in the axisymmetric part (see runs Turb in S.I.Table 3).Finally, we find that the Ohmic dissipation rate in these models provided enough thermal energy to cover the needs of young magnetars, (Fig. 4), even though we have made a conservative choice of the magnetic dipole field.By contrast, in cases where the evolution remains predominantly axisymmetric, instead of magnetic spots the field forms magnetic zones, where the intensity is only a factor of ∼ 3 stronger than the large-scale dipole field.Eventually, once a significant amount of energy has decayed (∼ 90% of the total energy) the magnetic evolution saturates, and any small-scale magnetic features become "frozen in" as in previous studies in Cartesian and axisymmetric simulations [41,30].
The diversity of observational manifestations of neutron stars has called for their Grand Unification [42] under a common theory, with their magnetic field being the key parameter [43,44,45,46].From our 3-D simulation survey we confirm indeed, that the initial structure and magnitude of the magnetic field are critical to the later evolution.We also find that the 3-D evolution is in general qualitatively different from the axisymmetric one.While predominantly poloidal initial magnetic field tend to remain axially symmetric [23], the inclusion of a strong, but not dominant, toroidal field is sufficient to break axial symmetry and has two major impacts, first: magnetic energy is converted faster to heat; second: small scale magnetic features of sizes appropriate for the hotspots observed in magnetars [10] form spontaneously and persist for several 10 4 yrs.This result is supportive of the importance of a second parameter, namely the toroidal magnetic field.However, contrary to the axisymmetric studies which suggested that most of the energy is concealed in the toroidal field in order to explain bursts [18] and the formation of magnetic spots [19], a toroidal field containing the same amount of energy as the poloidal field is sufficient to power magnetars.This leads to a more economical magnetar theory compared to previous works.Moreover, it restores the theoretical consistency between the structure of the magnetic field in magnetars and the fact that the poloidal and toroidal fields inside a star need to be comparable [47,48] a result that has been confirmed numerically in MHD equilibria and dynamo studies [32,16].While the magnetic evolution in the crust is now well understood, it is critical for future studies to address fully the role of the core and its coupling to the crust, as the existing core studies focus on axisymmetric structures [49,28,50].Furthermore, it is important to assess the impact of a magnetic field penetrating the crust-core boundary in the overall stability and Hall drift timescale.Similarly, it is possible that the strong magnetic field may plastically deform the crust, impeding the Hall evolution as the magnetic field lines are no longer frozen into the electron fluid [51].Regarding the thermal evolution, while it was shown here that the Ohmic decay provides sufficient energy to power magnetar X-ray Luminosity, it is important to couple magnetic and thermal evolution [46] in a single 3-D calculation.Such a calculation will consider the possibilities of heat transfer deeper into the crust, suppression of radiation because of the non-radial magnetic field components, or losses to neutrinos.Finally other parameters such as the neutron star mass and crust thickness need to be taken into account in a wider exploration of physical mechanisms.

Conclusions
The observational implications of these results for magnetars are in three basic directions.First, we have shown that the Hall effect has the tendency to produce strong, small-scale magnetic fields.This process does not appeal to an external source of energy, as the Hall effect exactly conserves magnetic energy, but rather to redistribution of the existing magnetic energy.The local magnetic field strength can exceed 10 15 G for a background dipole that is an order of magnitude weaker.Local magnetic fields of such strength are sufficient to reach the breaking stress [52,53] and produce bursting activity.Second, the enhanced Ohmic dissipation within the localised features generates sufficient heat to power their thermal X-ray luminosity.Therefore, the magnetic field will be used more efficiently to produce heat, i i "pnas˙arxiv" -2018/2/12 -18:16 -page 5 -#5 i i i i i i making the overall energy requirements smaller, while the fact that the strong magnetic field is localised in smaller areas, gives a viable explanation for the existence of hotspots.Finally, the evolution of the magnetic field, under the Hall effect, saturates once the energy has decayed to ∼ 10% of its initial value, within a few 10 5 yrs; in this case the neutron star still hosts a strong magnetic field ∼ 5 × 10 13 G, but its slow magnetic evolution reduces the chances of bursts and flares.Therefore the most energetic behaviour occurs while the star is young, particularly for magnetars with very strong fields.Fig. 3.The maximum magnetic field strength (dashed line) at the surface of the star versus the dipole field in three simulations.In the two cases where the nonaxisymmetric instability is triggered (red and blue lines), the strongest field on the surface of the star exceeds the dipole component by an order of magnitude at a very early time, as opposed to the model that remains axially symmetric (green lines).
Fig. 4. The Ohmic heating rate versus time for the same simulations as Figure 3 (red, green and blue), each initially containing 3.6 × 10 47 erg of magnetic energy.Also shown two simulations: one with an = 1 poloidal and an = 2 toroidal field (yellow line) where the latitudinal and azimuthal components decrease monotonically in the crust and the field remains mostly axisymmetric; and a simulation where the structure of the field is identical to that of the blue line, in which the Hall effect has been artificially suppressed (dashed line).The dots represent observations of X-ray luminosities for magnetars in the McGill Magnetar Catalogue [54] plotted against their characteristic ages (or the age of the associated supernova remnant if available), see S.I.Fig. S1.Ratio of the energy in the non-axisymmetric part of the field over the total magnetic energy for six runs.In all runs except for the DS00-3, which corresponds to a poloidal Hall equilibrium supported by a uniformly rotating electron fluid, the amount of energy in the non-axisymmetric part of the field grows exponentially.The exponential growth is considerably faster when V pu − V tu potentials are utilised (dotted lines) compared to V ps − V ts and when the amount of energy in the toroidal field is larger.

Supplementary Information
We use the following naming convention: Simulations whose names start with DS correspond to combinations of V ps and V ts potentials, DU correspond to combinations of V pu and V tu potentials, QS correspond to combinations of V ps and V ts,q potentials, QU correspond to combinations of V pu and V tu,q and those with Turb correspond to highly non-axisymmetric structures.The two numbers following the letters correspond to e t , the ratio of the initial toroidal energy to the total energy, with 00 corresponding to purely poloidal field and 10 to purely toroidal.The number following the hyphen is related to the increments of the normalisation factor B 0 , while simulations ending −0 correspond to runs where the Hall term is switched off and the evolution is purely Ohmic.Note that the simulations QS05-2M and QU05-2M utilise the negative sign for the potential V ts,q and V tu,q respectively.In the tables we provide the following information.The initial conditions are given in the first three columns (1-3); the following three columns (4)(5)(6) show the time it takes for the energy to decay to 0.5 (t 0.5 ), 0.1 (t 0.1 ) and 0.01 (t 0.01 ) of its initial value; the rest of the table contains the values of the fraction of energy in the non-axisymmetric part of the field E n /E tot (columns 7-10), and the ratio of the energy in the axisymmetric toroidal field over the total energy in the axisymmetric part of the field E φ /E ax (columns 11-14), at the instances determined by the decay levels mentioned before.Simulations appearing in Tables 1 and 2 start with setups corresponding to initial conditions of potentials V ps , V ts (DS) and V pu , V tu (DU) where the amount of energy in the non-axisymmetric part of the field is less than 10 −7 of the total magnetic energy (except for the simulations that evolve only under the effect of Ohmic dissipation where the energy in the non-axisymmetric part of the field is 10 −4 of the total energy).Simulations in Table 3 correspond to systems where the energy in the non-axisymmetric part is initially ∼ 10 −4 of the total energy.
In the figures we illustrate the basic results of our simulations.In Fig. S1 we plot the fraction of the energy in the non-axisymmetric part of the magnetic field.We find that (with the exception of the purely poloidal initial condition DS00-3) in all other cases the energy in the non-axisymmetric part grows exponentially.In Fig. S2 we plot the fraction of the energy in the axisymmetric toroidal component for various runs, showing that a poloidal field does not generate a dominant toroidal field, and it actually even suppresses the initial toroidal field.In Figs.(S3-S5) we plot the energy decomposition in and m modes of the spherical harmonics for various runs.

Fig. 1 .
Fig. 1.Magnetic field lines in the simulation with an initially differentially twisted magnetic field consisting of an = 1 (dipole) poloidal and an = 2 toroidal field, starting with only 10 −4 of the magnetic energy in the non-axisymmetric field, (run QU05-4, S.I.Table3).The snapshot is at t = 15kyr overplotted with the magnetic energy density e m = B 2 /(8π) on part of the surface.The surface field is highly

Table 1 .
Fig.S2.Ratio of the energy in the axisymmetric toroidal component of the field over the total energy in the axisymmetric part of the field for eight different models.The DU00-1,2,3,4 runs start with a purely poloidal field, with intensities doubling between consecutive runs.The magnetic field configuration generates a toroidal component, as it is out of Hall equilibrium, however, the toroidal component does not become dominant, with stronger fields suppressing it efficiently and remaining effectively poloidal.Note that once a significant fraction of the field has decayed and the Ohmic term becomes strong, the fraction of energy in the toroidal field rebounds and increases (see for instance the red solid curve).Once the energy is equally split between the toroidal and the poloidal components, in the runs DU05-1,2,3,4 plotted with dotted lines, stronger fields again tend to suppress the toroidal field.Simulations corresponding to V ps and V ts initial conditions.Times are expresses in 10 3 years./Etot(t0)En /E tot (t 0.5 ) E n /E tot (t 0.1 ) E n /E tot (t 0.01 ) E φ /E ax (t 0 ) E φ /E ax (t 0.5 ) E φ /E ax (t 0.1 ) E φ /E ax (t 0.01 ) Power spectrum of the Y m decomposition of the fields for run QU05-4, at times t 0 , t 0.5 , t 0.1 and t 0.01 .Top panel: The spectrum, starting from a state where the energy is equally distributed in a poloidal = 1 field and a toroidal = 2 field.Early evolution excites higher modes with local maxima at characteristic scales, which are indicative of the zonal formations (please refer to Figure3of the main paper).Bottom panel: The m spectrum; initially the energy is dominated by the axisymmetric mode m = 0 and the system is in equilibrium in the azimuthal direction.Non-axisymmetric modes grow rapidly which is indicative of the unstable behaviour.Power spectrum as in FigureS3for run DU10-4.Top panel: The spectrum, starting from a state where all energy is contained in an = 1 toroidal mode.Hall evolution pushes energy into higher 's especially early in the neutron star's life.Bottom panel: The m spectrum, where initially the energy is dominated by the axially symmetric component m = 0, however the development of Hall instability generates features with local maxima at m = 7, 14, 17, 21, 28, 31.Power spectrum as in FigureS3for run Turb-4.Top panel: The spectrum, starting from a state where the energy is distributed between 10 ≤ ≤ 40.Energy is transferred both to smaller and larger scales, even to = 1, with the fraction of energy in this mode growing with time.Bottom panel: The m spectrum, initially the energy is set to be equally distributed between 0 ≤ m ≤ 10 and then drops obeying a power-law.

Table 3 .
Simulations starting with 10 −4 of the total energy being in the non-axisymmetric component of the field, times expressed in 10 3 years./E tot (t 0 ) E n /E tot (t 0.5 ) E n /E tot (t 0.1 ) E n /E tot (t 0.01 ) E φ /E ax (t 0 ) E φ /E ax (t 0.5 ) E φ /E ax (t 0.1 ) E φ /E ax (t 0.01 )