Magnetic field evolution in magnetar crusts through three-dimensional simulations

Edited by Neta A. Bahcall, Princeton University, Princeton, NJ, and approved February 26, 2016 (received for review November 12, 2015)
March 28, 2016
113 (15) 3944-3949


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 a field is unclear. Our fully 3D 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.


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 modeling of the evolution of the magnetic field in the crust of a neutron star through 3D 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 nonaxisymmetric, kilometer-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 localized enhancement of Ohmic heating can power the star’s persistent emission. Thus, the observed diversity in magnetar behavior can be explained with mixed poloidal−toroidal fields of comparable energies.
An 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 quantum electrodynamic (QED) magnetic field 4.4 × 1013 G are conventionally called magnetars, and typically exhibit highly energetic behavior, 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 behavior (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 (68). Furthermore, despite the high thermal conductivity of neutron stars’ solid outer crusts (9), observations of their thermal emission indicate that, in some cases, the surface is highly anisothermal (10), with kilometer-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 1013 G. It seems likely, then, that the strong fields in magnetars must result from a combination of differential rotation and dynamo action before 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 in 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 thermal radiation. However, axisymmetric simulations of the crustal magnetic field only generate such features if a toroidal magnetic field exceeding 1016 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 localized patches of strong magnetic field via magnetic instabilities. Such instabilities have been demonstrated in local, plane-parallel models (2022) but not yet in a realistic global model.

Hall Evolution in Neutron Star Crusts

Motivated by this puzzle, we study the fully nonlinear 3D problem of the magnetic field evolution in a magnetar crust through numerical simulations, using a modified version (23) of the parallel rotating dynamo (PARODY) code (24, 25), to determine the necessary conditions for the spontaneous generation of strong localized magnetic field out of a large-scale 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−magnetohydrodynamics (MHD) with subdominant Ohmic dissipation (17). This process is described by the Hall−MHD induction equation,
where B is the magnetic induction, ne is the electron number density, σ is the electric conductivity, c is the speed of light, and e is the elementary charge. The first term in the right-hand-side of Eq. 1 describes the Hall effect, and the second one describes Ohmic dissipation. We assume a neutron star radius R=10 km and a crust thickness of 1 km. We shall express the quantities in terms of the normalized radial distance r=R/R, where R is the distance from the center. The electron number density is taken to be ne=2.5×1034 cm[(1.0463r)/0.0463]34, and the electric conductivity is σ=1.8×1023 s{[(1.0463r)/0.0463]}18/3. Such profiles are good analytical fits of the form σne2/3 of more precise crust models (26) at temperatures 108K.
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, whereas, for the nonlinear 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 (see figures 3 and 4 of ref. 28, and figure 9 of ref. 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 3D 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 against 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 Ym up to max,mmax=6080. We have repeated some of the runs with the strongest magnetic fields, using twice the angular resolution, with max=mmax=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 Eq. 1 in each of the three scenarios for the initial magnetic field in the star’s crust (see Tables S1S3). 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 normalization of the magnetic field intensity. The intensity of the poloidal field on the surface takes values between 1013 G and 5 × 1014 G, except for the purely toroidal initial conditions.
Table S1.
Simulations corresponding to Vps and Vts initial conditions
Times are expressed in 103 y.
Table S2.
Simulations corresponding to Vpu and Vtu initial conditions
Times are expressed in 103 y.
Table S3.
Simulations starting with 104 of the total energy being in the nonaxisymmetric component of the field
Times are expressed in 103 y.
The magnetic field is expressed in terms of two functions Vp(r,θ,ϕ) and Vt(r,θ,ϕ),
We express the fraction of the energy in the toroidal component of the magnetic field et, allowing the following values: 0 (purely poloidal), 0.1, 0.5, 0.9, and 1 (purely toroidal). B0 governs the overall amplitude and takes the values 0.5, 1, 2, and 4 expressed in units of 1014 G; the normalization is chosen so that the initial magnetic field energy for a simulation with amplitude B0 is 2.25×1046B02 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
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 B0 for et=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=15B0et1/2. Thus, for the strongest toroidal field simulated (B0=4, et=1), the maximum azimuthal field is 6 × 1015 G.
The second family of =1 profiles is the following:
These profiles have Bθ and Bϕ components with an extremum close to the center of the crust. The choice of normalization is such that the magnetic energy in the crust is the same as the corresponding Vps and Vts 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=12B0et1/2. Thus, for the strongest toroidal field simulated (B0=4, et=1), the azimuthal field has a maximum value of 4.8 × 1015 G.
We also combined the poloidal dipole fields with =2 toroidal fields. In such cases, we have used the following profiles:
which is used in conjunction with Vps, and
which is used in conjunction with Vpu. In the simulations with the quadrupolar toroidal field, we assumed equipartition in energy between the poloidal and the toroidal field (et=0.5) and varied B0=0.5,1,2,4 for the case of the positive sign, whereas we run simulations with B0=2 only for the case of the negative sign. The maximum value of the azimuthal field assumed, using the Vts,q profile, is 5 × 1015 G and occurs at the base of the crust at latitudes θ=π/4,3π/4; for the case of of Vtu,q, the maximum value assumed is 3.8 × 1015 G and occurs at r=0.95 and θ=π/4,3π/4. In all these axially symmetric initial conditions, we have superimposed a nonaxisymmetric random noise, containing 104108 of the total energy. Axially symmetric initial conditions correspond to equilibria with respect to the ϕ coordinate, so any growth of the nonaxisymmetric part of the field may indicate the presence of a nonaxisymmetric instability.
Finally, we have simulated a highly nonaxisymmetric magnetic field. In this case, we have populated the crust with a fully confined nonaxisymmetric magnetic field, where 0.9 of the total magnetic energy is in the nonaxisymmetric part of the field, consisting of multipoles with 1040. Similarly to the previous cases, we have used four normalization 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:
In systems where the initial azimuthal field is strong (et0.5), and especially when the profiles Vpu and Vtu are used (Figs. 1 and 2), the amount of energy in the nonaxisymmetric part of the field increases exponentially. The growth rate is approximately proportional to the strength of the magnetic field (see Fig. S1), i.e., for the case depicted in Figs. 1 and 2, the nonaxisymmetric energy growth timescale is ∼1.5 ky. The exponential growth stops if one of the following two conditions is fulfilled: either 0.50.9 of the total magnetic energy has decayed or the nonaxisymmetric field contains about 50% of the total energy. Although, initially, the m=0 mode was dominant, once the instability develops, energy is transferred, and a local maximum appears around 10m20 (see Figs. S2 and S3 for a purely toroidal field), with some further local maxima at higher harmonics, leading to a characteristic wavelength for the instability of about 5 km at the equator, with finer structure appearing later giving kilometer-sized features. These features are well above the resolution of the simulation and appear both in the low- and high-resolution runs. This behavior 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.
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 (Table S3), 50% of the energy decays within the first 20 ky. We find that the early decay approximately scales with the magnetic field strength (Tables S1S3, column t0.5), whereas the later decay has a much weaker dependence on the magnetic field strength.
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 Fig. S4). Even if a toroidal component exists in the initial conditions, stronger fields tend to suppress it more efficiently than weaker ones, in favor 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 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 megayears, as is illustrated in the last column of Tables S1S3 for the simulations with B0=1.
The evolution of the axisymmetric part of the magnetic field is in accordance with the results of previous axisymmetric numerical and semianalytical 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 toward the equator is due to the polarity of the toroidal field, an effect that has been discussed in detail in plane-parallel geometry in ref. 35 and the axially symmetric case in ref. 27.
Configurations initially dominated by the poloidal field tend to remain axially symmetric, especially if the Vps 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.
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).
Initial conditions corresponding to highly nonaxisymmetric structures undergo a weak inverse cascade (see runs Turb in the Supporting Information and Table S3). The Hall effect transfers energy across the spectrum to both smaller and larger scales (see Fig. S5). Despite the great differences in the structure of the magnetic field compared with 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.
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 nonaxisymmetric field (run QU05-4, Table S3). The snapshot is at t=15 ky overplotted with the magnetic energy density em=B2/(8π) on part of the surface. The surface field is highly anisotropic, with small regions in which the magnetic energy density exceeds the average surface value by at least an order of magnitude.
Fig. 2.
Contour plot of the radial (Left), latitudinal (Middle), and azimuthal (Right) component of the magnetic field at r=0.995R at t=0 (Top) and t=15 ky (Bottom) for the simulation shown in Fig. 1.
Fig. S1.
Ratio of the energy in the nonaxisymmetric part of the field over the total magnetic energy for six runs. In all runs except for DS00-3, which corresponds to a poloidal Hall equilibrium supported by a uniformly rotating electron fluid, the amount of energy in the nonaxisymmetric part of the field grows exponentially. The exponential growth is considerably faster when VpuVtu potentials are used (dotted lines) compared with VpsVts and when the amount of energy in the toroidal field is larger.
Fig. S2.
Power spectrum of the Ym decomposition of the fields for run QU05-4, at times t0, t0.5, t0.1, and t0.01. (Top) 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 Fig. 3). (Bottom) The m spectrum. Initially, the energy is dominated by the axisymmetric mode m=0 and the system is in equilibrium in the azimuthal direction. Nonaxisymmetric modes grow rapidly, which is indicative of the unstable behavior.
Fig. S3.
Power spectrum as in Fig. S2 but for run DU10-4. (Top) 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) The m spectrum, where, initially, the energy is dominated by the axially symmetric component m=0, but the development of Hall instability generates features with local maxima at m = 7, 14, 17, 21, 28, and 31.
Fig. S4.
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, DU00-2, DU00-3, and DU00-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, DU00-2, DU00-3, and DU00-4 plotted with dotted lines, stronger fields again tend to suppress the toroidal field.
Fig. S5.
Power spectrum as in Fig. S2 but for run Turb-4. (Top) The spectrum, starting from a state where the energy is distributed between 1040. Energy is transferred both to smaller and larger scales, even to =1, with the fraction of energy in this mode growing with time. (Bottom) The m spectrum. Initially, the energy is set to be equally distributed between 0m10 and then drops, obeying a power law.


In all simulations, we found 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 behaviors. Initial conditions where the poloidal field was dominant remain axially symmetric, whereas cases with substantial, but not necessarily dominant, azimuthal fields triggered nonaxisymmetric instabilities. As magnetic energy is conserved in Hall−MHD, this growth of the nonaxisymmetric field occurs at the expense of the axisymmetric components.
In the model with initial dipole poloidal magnetic field of 4 × 1014 G and a toroidal component containing the same amount of energy with the poloidal (Figs. 1 and 2), magnetic spots of sizes ∼2 km appear. These spots contain energies up to 5 × 1043 erg, while the dipole field at that time is 1−2 × 1014 G (Fig. 3). The rapid decrease of the dipole component, with a drop by a factor of 2 within 10–20 ky, is caused by the latitudinal Hall drift toward the equator. This is more efficient when the stronger fields are in shallow depths, when the potentials Vtu and Vtu,q are used (runs DU05-4 and QU05-4, as defined in the Supporting Information), where the maximum value of the toroidal field occurs at r=0.95. On the contrary, it takes 10 times longer for the same decrease in the dipole component, in the case of the potentials Vts (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 (5, 8). 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 ∼2 km is ∼2 × 1046 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 (3840). Note, however, that magnetars that have exhibited giant flares have spin-down dipole fields above 5 × 1014 G, so the overall magnetic energy should be scaled by a factor of 10 compared with our simulations. These structures take ∼104 y to develop, given the choice of an initial condition where only 10−4 of the total is in the nonaxisymmetric part. Had we chosen initial conditions with more energy in the nonaxisymmetric 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 103 y. This was confirmed in the extreme case where no energy is in the axisymmetric part (see runs Turb in Table S3). 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 (30, 41).
Fig. 3.
The maximum magnetic field strength (dashed lines) at the surface of the star versus the dipole field (solid lines) 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 Fig. 3 (red, green, and blue), each initially containing 3.6 × 1047 erg of magnetic energy. Also shown are 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 another 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 Table S4.
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 (4346). From our 3D 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 3D evolution is, in general, qualitatively different from the axisymmetric one. Although predominantly poloidal initial magnetic fields 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 104 y. This result is supportive of the importance of a second parameter, namely, the toroidal magnetic field. However, contrary to the axisymmetric studies that suggested, to explain bursts (18) and the formation of magnetic spots (19), that most of the energy is concealed in the toroidal field, 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 with 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 (16, 32). Although 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 (28, 49, 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, although 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 3D calculation. Such a calculation will consider the possibilities of heat transfer deeper into the crust, suppression of radiation because of the nonradial 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.


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 1015 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 localized features generates sufficient heat to power their thermal X-ray luminosity. Therefore, the magnetic field will be used more efficiently to produce heat, making the overall energy requirements smaller, whereas the fact that the strong magnetic field is localized 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 105 y; in this case, the neutron star still hosts a strong magnetic field, ∼5 × 1013 G, but its slow magnetic evolution reduces the chances of bursts and flares. Therefore, the most energetic behavior occurs while the star is young, particularly for magnetars with very strong fields.

SI Text

We use the following naming convention: Simulations whose names start with DS correspond to combinations of Vps and Vts potentials, DU corresponds to combinations of Vpu and Vtu potentials, QS corresponds to combinations of Vps and Vts,q potentials, QU corresponds to combinations of Vpu and Vtu,q, and those with Turb correspond to highly nonaxisymmetric structures. The two numbers following the letters correspond to et, the ratio of the initial toroidal energy to the total energy, with 00 corresponding to purely poloidal field and 10 corresponding to purely toroidal. The number following the hyphen is related to the increments of the normalization factor B0, and 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 use the negative sign for the potential Vts,q and Vtu,q, respectively.
In Tables S1S3, we provide the following information. The initial conditions are given in the first three columns (1–3); the following three columns (4–6) show the time it takes for the energy to decay to 0.5 (t0.5), 0.1 (t0.1), and 0.01 (t0.01) of its initial value; the rest of the table contains the values of the fraction of energy in the nonaxisymmetric part of the field En/Etot (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ϕ/Eax (columns 11–14), at the instances determined by the decay levels mentioned before. Simulations appearing in Tables S1 and S2 start with setups corresponding to initial conditions of potentials Vps, Vts (DS) and Vpu, Vtu (DU) where the amount of energy in the nonaxisymmetric 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 nonaxisymmetric part of the field is 10−4 of the total energy). Simulations in Table S3 correspond to systems where the energy in the nonaxisymmetric part is initially ∼10−4 of the total energy. In Table S4 we present the magnetar data used in Fig. 4 of the main text.
Table S4.
Magnetar data used in Fig. 4 of the main text
NameLxτc, kySNR age, ky
CXOU J010043.1–7211346.52E+346.8
4U 0142+611.05E+3568
SGR 0418+57299.57E+2936,000
SGR 0501+45168.14E+32154–7
SGR 0526–661.89E+353.4∼4.8
1E 1048.1–59374.94E+344.5
1E 1547.0–54081.31E+330.69
PSR J1622-49504.36E+324<6
SGR 1627–413.62E+332.2
CXOU J164710.2–4552164.55E+32420
1RXS J170849.0–4009104.20E+349
CXOU J171405.7–3810315.59E+340.90.650.3+2.5
SGR J1745-29001.07E+324.3
SGR 1806–201.63E+350.24
XTE J1810-1974.25E+3111
Swift J1822.3–16063.98E+296,300
Swift J1834.9–08468.44E+304.9∼100
1E 1841–0451.84E+354.60.5–1
3XMM J185246.6+0033176.03E+301,300
SGR 1900+148.97E+340.9
1E 2259+5861.73E+3423014(2)
PSR J1845-02581.85E+340.73
The data are taken from the McGill Magnetar Catalogue (54). The error bars used in the plot are based on the distance uncertainties. With respect to age, the associated supernova remnant age is used in favor of the characteristic spin-down age τc, if the former is known.
In Figs. S1S5, we illustrate the basic results of our simulations. In Fig. S1, we plot the fraction of the energy in the nonaxisymmetric 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 nonaxisymmetric part grows exponentially. In Fig. S4, 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. S2, S3, and S5, we plot the energy decomposition in and m modes of the spherical harmonics for various runs.


We thank Andrew Cumming, Vicky Kaspi, Scott Olaussen, Yuri Levin, Anatoly Spitkovsky, and Andrei Beloborodov for discussions and comments. We thank an anonymous referee for stimulating comments that improved the paper. K.N.G. and R.H. were supported by Science and Technology Funding Council (STFC) Grant ST/K000853/1. The numerical simulations were carried out on the STFC-funded Distributed Research Utilizing Advanced Computing (DiRAC) I UKMHD Science Consortia machine, hosted as part of and enabled through the Advanced Research Computing (ARC) High Performance Computing (HPC) resources and support team at the University of Leeds. We acknowledge the use of the McGill Magnetar Catalogue (∼pulsar/magnetar/main.html).

Supporting Information

Supporting Information (PDF)
Supporting Information


AJ Deutsch, The electromagnetic field of an idealized star in rigid rotation in vacuo. Ann Astrophys 18, 1–10 (1955).
C Thompson, RC Duncan, The soft gamma repeaters as very strongly magnetized neutron stars - I. Radiative mechanism for outbursts. Mon Not R Astron Soc 275, 255–300 (1995).
C Thompson, RC Duncan, The soft gamma repeaters as very strongly magnetized neutron stars. II. Quiescent neutrino, X-ray, and Alfven wave emission. Astrophys J 473, 322–342 (1996).
F Haberl, The magnificent seven: magnetic fields and surface temperature distributions. Astrophys Space Sci 308, 181–190 (2007).
N Rea, et al., A low-magnetic-field soft gamma repeater. Science 330, 944–946 (2010).
H An, VM Kaspi, R Archibald, A Cumming, Spectral and timing properties of the magnetar CXOU J164710.2-455216. Astrophys J 763, 82 (2013).
P Scholz, et al., Post-outburst X-ray flux and timing evolution of swift J1822.3-1606. Astrophys J 761, 66 (2012).
FP Gavriil, VM Kaspi, Long-term Rossi X-Ray Timing Explorer monitoring of anomalous X-ray pulsars. Astrophys J 567, 1067–1076 (2002).
AY Potekhin, DA Baiko, P Haensel, DG Yakovlev, Transport properties of degenerate electrons in neutron star envelopes and white dwarf cores. Astron Astrophys 346, 345–353 (1999).
S Guillot, R Perna, N Rea, D Viganò, JA Pons, Modelling of the surface emission of the low magnetic field magnetar SGR 0418+5729. Mon Not R Astron Soc 452, 3357–3368 (2015).
F Bernardini, et al., Emission geometry, radiation pattern and magnetic topology of the magnetar XTE J1810-197 in its quiescent state. Mon Not R Astron Soc 418, 638–647 (2011).
A Tiengo, et al., A variable absorption feature in the X-ray spectrum of a magnetar. Nature 500, 312–314 (2013).
T Güver, E Göǧüş, F Özel, A magnetar strength surface magnetic field for the slowly spinning down SGR 0418+5729. Mon Not R Astron Soc 418, 2773–2778 (2011).
GA Rodríguez Castillo, et al., The outburst decay of the low magnetic field magnetar SWIFT J1822.3-1606: phase-resolved analysis and evidence for a variable cyclotron feature. Mon Not R Astron Soc 456, 4145–4155 (2016).
HC Spruit, Origin of neutron star magnetic fields. AIP Conf Proc 983, 391–398 (2008).
P Mösta, et al., A large-scale dynamo and magnetoturbulence in rapidly rotating core-collapse supernovae. Nature 528, 376–379 (2015).
P Goldreich, A Reisenegger, Magnetic field decay in isolated neutron stars. Astrophys J 395, 250–258 (1992).
JA Pons, R Perna, Magnetars versus high magnetic field pulsars: A theoretical interpretation of the apparent dichotomy. Astrophys J 741, 123 (2011).
U Geppert, D Viganò, Creation of magnetic spots at the neutron star surface. Mon Not R Astron Soc 444, 3198–3208 (2014).
M Rheinhardt, D Konenkov, U Geppert, The occurrence of the Hall instability in crusts of isolated neutron stars. Astron Astrophys 420, 631–645 (2004).
TS Wood, R Hollerbach, M Lyutikov, Density-shear instability in electron magneto-hydrodynamics. Phys Plasmas 21, 052110 (2014).
KN Gourgouliatos, T Kondić, M Lyutikov, R Hollerbach, Magnetar activity via the density-shear instability in Hall-MHD. Mon Not R Astron Soc 453, L93–L97 (2015).
TS Wood, R Hollerbach, Three dimensional simulation of the magnetic stress in a neutron star crust. Phys Rev Lett 114, 191101 (2015).
E Dormy, P Cardin, D Jault, MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field. Earth Planet Sci Lett 160, 15–30 (1998).
J Aubert, J Aurnou, J Wicht, The magnetic structure of convection-driven numerical dynamos. Geophys J Int 172, 945–956 (2008).
A Cumming, P Arras, E Zweibel, Magnetic field evolution in neutron star crusts due to the Hall effect and Ohmic decay. Astrophys J 609, 999–1017 (2004).
R Hollerbach, G Rüdiger, Hall drift in the stratified crusts of neutron stars. Mon Not R Astron Soc 347, 1273–1278 (2004).
KT Henriksson, I Wasserman, Poloidal magnetic fields in superconducting neutron stars. Mon Not R Astron Soc 431, 2986–3002 (2013).
SK Lander, The contrasting magnetic fields of superconducting pulsars and magnetars. Mon Not R Astron Soc 437, 424–436 (2014).
KN Gourgouliatos, A Cumming, Hall attractor in axially symmetric magnetic fields in neutron star crusts. Phys Rev Lett 112, 171101 (2014).
KN Gourgouliatos, A Cumming, Hall effect in neutron star crusts: Evolution, endpoint and dependence on initial conditions. Mon Not R Astron Soc 438, 1618–1629 (2014).
J Braithwaite, HC Spruit, A fossil origin for the magnetic field in A stars and white dwarfs. Nature 431, 819–821 (2004).
C Thompson, N Murray, Transport of magnetic fields in convective, accreting supernova cores. Astrophys J 560, 339–357 (2001).
KN Gourgouliatos, et al., Hall equilibria with toroidal and poloidal fields: Application to neutron stars. Mon Not R Astron Soc 434, 2480–2490 (2013).
SI Vainshtein, SM Chitre, AV Olinto, Rapid dissipation of magnetic fields due to the Hall current. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 61, 4422–4430 (2000).
S Mereghetti, et al., Strong bursts from the anomalous X-ray pulsar 1E 1547.0-5408 observed with the INTEGRAL/SPI anti-coincidence shield. Astrophys J 696, L74–L78 (2009).
AJ van der Horst, et al., SGR J1550-5418 bursts detected with the Fermi Gamma-Ray Burst Monitor during its most prolific activity. Astrophys J 749, 122 (2012).
EP Mazets, SV Golentskii, VN Ilinskii, RL Aptekar, IA Guryan, Observations of a flaring X-ray pulsar in Dorado. Nature 282, 587–589 (1979).
K Hurley, et al., A giant periodic flare from the soft γ-ray repeater SGR1900+14. Nature 397, 41–43 (1999).
DM Palmer, et al., A giant γ-ray flare from the magnetar SGR 1806-20. Nature 434, 1107–1109 (2005).
CJ Wareing, R Hollerbach, Forward and inverse cascades in decaying two-dimensional electron magnetohydrodynamic turbulence. Phys Plasmas 16, 042307 (2009).
VM Kaspi, Grand unification of neutron stars. Proc Natl Acad Sci USA 107, 7147–7152 (2010).
DN Aguilera, JA Pons, JA Miralles, 2D cooling of magnetized neutron stars. Astron Astrophys 486, 255–271 (2008).
JA Pons, U Geppert, Magnetic field dissipation in neutron star crusts: From magnetars to isolated neutron stars. Astron Astrophys 470, 303–315 (2007).
JA Pons, JA Miralles, U Geppert, Magneto-thermal evolution of neutron stars. Astron Astrophys 496, 207–216 (2009).
D Viganò, et al., Unifying the observational diversity of isolated neutron stars via magneto-thermal evolution models. Mon Not R Astron Soc 434, 123–141 (2013).
KH Prendergast, The equilibrium of a self-gravitating incompressible fluid sphere with a magnetic field. I. Astrophys J 123, 498−507 (1956).
E Flowers, MA Ruderman, Evolution of pulsar magnetic fields. Astrophys J 215, 302–310 (1977).
SK Lander, Magnetic fields in superconducting neutron stars. Phys Rev Lett 110, 071101 (2013).
JG Elfritz, JA Pons, N Rea, K Glampedakis, D Vigan, Simulated magnetic field expulsion in neutron star cores. Mon Not R Astron Soc 456, 4461–4474 (2016).
AM Beloborodov, Y Levin, Thermoplastic waves in magnetars. Astrophys J 794, L24 (2014).
C Thompson, RC Duncan, The giant flare of 1998 August 27 from SGR 1900+14. II. Radiative mechanism and physical constraints on the source. Astrophys J 561, 980–1005 (2001).
SK Lander, N Andersson, D Antonopoulou, AL Watts, Magnetically driven crustquakes in neutron stars. Mon Not R Astron Soc 449, 2047–2058 (2015).
SA Olausen, VM Kaspi, The McGill Magnetar Catalog. Astrophys JS 212, 6 (2014).

Information & Authors


Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 113 | No. 15
April 12, 2016
PubMed: 27035962


Submission history

Published online: March 28, 2016
Published in issue: April 12, 2016


  1. neutron stars
  2. magnetars
  3. pulsars
  4. magnetohydrodynamics


We thank Andrew Cumming, Vicky Kaspi, Scott Olaussen, Yuri Levin, Anatoly Spitkovsky, and Andrei Beloborodov for discussions and comments. We thank an anonymous referee for stimulating comments that improved the paper. K.N.G. and R.H. were supported by Science and Technology Funding Council (STFC) Grant ST/K000853/1. The numerical simulations were carried out on the STFC-funded Distributed Research Utilizing Advanced Computing (DiRAC) I UKMHD Science Consortia machine, hosted as part of and enabled through the Advanced Research Computing (ARC) High Performance Computing (HPC) resources and support team at the University of Leeds. We acknowledge the use of the McGill Magnetar Catalogue (∼pulsar/magnetar/main.html).


This article is a PNAS Direct Submission.



Konstantinos N. Gourgouliatos1 [email protected]
Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, United Kingdom;
Toby S. Wood
School of Mathematics and Statistics, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom
Rainer Hollerbach
Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, United Kingdom;


To whom correspondence should be addressed. Email: [email protected].
Author contributions: K.N.G. designed research; K.N.G. performed research; K.N.G., T.S.W., and R.H. analyzed data; and K.N.G., T.S.W., and R.H. wrote the paper.

Competing Interests

The authors declare no conflict of interest.

Metrics & Citations


Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Citation statements



If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by


    View Options

    View options

    PDF format

    Download this article as a PDF file


    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    Magnetic field evolution in magnetar crusts through three-dimensional simulations
    Proceedings of the National Academy of Sciences
    • Vol. 113
    • No. 15
    • pp. 3903-E2208







    Share article link

    Share on social media