Nonaxisymmetric Hall instability: A key to understanding magnetars

It is generally accepted that the non-linear, dynamical evolution of magnetic fields in the interior of neutron stars plays a key role in the explanation of the observed phenomenology. Understanding the transfer of energy between toroidal and poloidal components, or between different scales, is of particular relevance. In this letter, we present the first 3D simulations of the Hall instability in a neutron star crust, confirming its existence for typical magnetar conditions. We confront our results to estimates obtained by a linear perturbation analysis, which discards any interpretation as numerical instabilities and confirms its physical origin. Interestingly, the Hall instability creates locally strong magnetic structures that occasionally can make the crust yield to the magnetic stresses and generates coronal loops, similarly as solar coronal loops find their way out through the photosphere. This supports the viability of the mechanism, which has been proposed to explain magnetar outbursts.

It is generally accepted that the non-linear, dynamical evolution of magnetic fields in the interior of neutron stars plays a key role in the explanation of the observed phenomenology. Understanding the transfer of energy between toroidal and poloidal components, or between different scales, is of particular relevance. In this letter, we present the first 3D simulations of the Hall instability in a neutron star crust, confirming its existence for typical magnetar conditions. We confront our results to estimates obtained by a linear perturbation analysis, which discards any interpretation as numerical instabilities and confirms its physical origin. Interestingly, the Hall instability creates locally strong magnetic structures that occasionally can make the crust yield to the magnetic stresses and generates coronal loops, similarly as solar coronal loops find their way out through the photosphere. This supports the viability of the mechanism, which has been proposed to explain magnetar outbursts.
Among the many observational faces of neutron stars, magnetars -highly magnetized, slowly rotating, isolated neutron stars that sporadically show violent transient activity -have received much attention in the last decade. Understanding their puzzling behaviour, arguably caused by their magnetic activity, is one of the most active research areas of neutron star physics. In the magnetar scenario, one usually appeals to the creation or existence of small scale magnetic structures to justify the observed phenomenology [1][2][3][4]. For instance, observational evidence [5] favors small structures emerging from the surface, similar to the sunspots in the solar corona [6]. The details about how exactly such small-scale magnetic structures are created are under debate, but they must have their origin in the dynamics of the interior, in particular of the neutron star crust.
The evolution of the magnetic field in a neutron star (NS) crust is governed by the combined action of Ohmic dissipation and the Hall drift [7,8]. In the presence of a strong field, as in magnetars, the non-linear Hall drift dominates and has been proposed to be the main responsible for the generation of small scales through the so-called Hall instability [9] (hereafter RG02). The occurrence of the Hall instability in a real neutron star has been somewhat controversial, in part because of the lack of numerical simulations able to reproduce the strongfield regime under realistic conditions. The first 2D simulations of the magnetic field evolution in NS crusts [10][11][12][13][14] did not find strong evidence of such instability, but they were restricted to axisymmetric models and moderate magnetization, for numerical reasons. A first 3D non-linear study in a periodic box [15] concluded that any instabilities are overwhelmed by a turbulent Hall cascade. Conversely, [16] confirmed the occurrence of the instability showing that, because the unstable modes have a relatively long wavelength, it is suppressed on a cubic domain with periodic boundary conditions, but it arises on a thin slab where one of the spatial lengths is longer than others. And this is precisely the geometry of a neutron star crust, a thin spherical shell (≈ 1 km) of radius R ≈ 12 km.
We must note the distinction between the resistive Hall instability, which is essentially a tearing mode [17,18], and ideal instabilities which operate under infinite conductivity, for example, the density-shear instability requiring a density gradient [19,20], or the fast collisionless reconnection observed in the whistler frequency range [21]. Indeed, the growth times of the Hall instability modes become increasingly large with vanishing resistivity. Another relevant issue is that, in spherical geometry, the instability may be suppressed for the axisymmetric modes. Thus, 2D simulations could not help resolving the controversy, and only recently [22][23][24], 3D simulations have been possible.
In this letter, we present the first 3D simulations of the Hall instability for a model with a strong toroidal field in a NS crust, confirming its occurrence. We show how small magnetic structures are naturally created and drift toward the star's poles. Although similar structures had been observed before in previous simulations [22,23,25], their actual origin and the true nature of a possible instability had not been settled. We further present a linear stability analysis in a spherical shell that gives similar results to the non-linear simulations, thus reinforcing our conclusion that the Hall instability is actually at the origin of the observed magnetar activity.
To a very good approximation, the crust can be considered a one component plasma, where only electrons can flow through the ionic lattice. Since ions are fixed, there is no mass motion, and the induction equation reduces to: where c is the speed of light, η = c 2 /4πσ is the magnetic diffusivity, σ is the electrical conductivity, e is the elementary charge and n e is the electron number density. We consider a background toroidal field of the form and evolve it by numerically solving eq. (1) inside a spherical shell, using a version of the PARODY 3-D MHD code [26,27] suitably adapted to NSs [22,23]. We impose vacuum boundary conditions at the exterior of the star and superconductor boundary conditions at the base of the crust, not allowing the magnetic field to penetrate into the core. We consider a crust of uniform electron number density n e = 2.5 × 10 34 cm −3 , electric conductivity σ = 1.8 × 10 23 s −1 , and we express the magnetic field in units of B 14 = B/10 14 G. We have decided to study simplified models with constant density and conductivity for two reasons. First, this allows us to distinguish between the families of resistive and ideal instabilities that may operate in the electron-MHD regime. Namely, a realistic crust with stratified density and composition is subject to both instabilities, whereas our set up suppresses the ideal density-shear instability [20], but permits the resistive one [9] to grow. Second, this choice also allows us a more direct comparison to previous works. To explore how our results depend on the thickness of the layer where the magnetic field is confined, we have considered two cases: a realistic one where the inner radius of the crust R in = 0.9R N S (hereafter model A) and another with R in = 0.5R N S (hereafter model B). For structures with a length-scale L, measured in km, the Hall timescale is t H = 4πen e L 2 /cB = 16L 2 /B 14 kyr, and the Ohmic dissipation timescale t D = 4πσL 2 /c 2 = 0.8L 2 Myr. The ratio t D /t H (equivalent to the magnetic Reynolds number) is 50B 14 .
In general, purely azimuthal magnetic fields, as our background field (2), cannot be in Hall equilibrium; however, their evolution maintains their azimuthal structure and does not generate any radial or meridional compo- nent, a result that has been confirmed by axisymmetric simulations. This is also the case in 3-D simulations provided that non-axisymmetric perturbations are not included in the initial conditions. Therefore, a strong indication of the operation of the resistive Hall instability is the growth of non-axisymmetric structures, even if the background field drifts at the same time in the meridional and radial direction.
To test this hypothesis, we have performed 3-D simulations of the initial field given by eq. (2) with whereB is a normalisation parameter, chosen so that the maximum value attained by the magnetic field inside the crust is 2 × 10 15 G. We have used this particular profile due to its similarity to the profile studied in RG02. Furthermore, we also include non-axisymmetric perturbations, with a flat spectrum, exciting the azimuthal l = m modes from 1 to 40. Starting from these initial conditions, we find that there is a continuous growth of the non-axisymmetric modes with the appearance of zones with alternating inwards and outwards magnetic field. While this is happening, the magnetic field drifts towards the northern hemisphere, as in the axisymmetric simulations.
In Fig. 1 we show an illustration of a snapshot during the evolution of model A. Besides the global drift toward the north pole, we find that many small-scale structures arise. Locally, the magnetic field intensity in some regions is about one order of magnitude higher than the average value. This figure corresponds to model A, but the qualitative behaviour is similar in both models. However, there is a clear quantitative distinction. In model A, the number of structures in the azimuthal direction is 18, whereas in model B the number of zones is 8 as we can see in Fig. 2, where we compare their internal structure (at a radius of r = 0.95R N S and r = 0.75R N S for models A and B respectively), after the saturation of the instabilities. Besides, the growth time is consistently slower for model B. This pattern is visible in spectral space, as the peak of the cumulative distribution is m ∼ 20 for the thin crust (Fig. 3) and m ∼ 10 for the thicker one.
To confirm that the results of non-linear simulations actually correspond to the so-called Hall-instability, we have performed a linear perturbation analysis of a background toroidal field of the same form as eq. (2), in a spherical shell. We decompose the perturbation b into its poloidal and toroidal components b = b p + b t , which in turn can be expressed in terms of two scalar functions S, T as follows: where e r is the unit vector in the radial direction. Following the standard procedure, but in spherical coordinates, we expand the scalar functions in spherical harmonics: where we use τ to denote the time variable to avoid confusion with the toroidal radial function t(r). This notation allows us to directly compare with RG02. Here, p is a complex number to be determined. The eigenvalue p with the largest positive real part represents the fastest growing mode (with growth time 1/Re(p)). With help of an algebraic manipulator, we obtain the following equations for the perturbations: where the primes denote radial derivatives and C m l±1 and D m l±1 are a shorthand for the coupling terms with the l±1 coefficients (always with the same m index). The C m l±1 term only involves combinations of f s l±1,m and f s l±1,m , while the D m l±1 only involves the product f t l±1,m and its first radial derivative.
For reference, in Cartesian coordinates and planar symmetry, RG02 expanded the perturbation in plane waves and obtained The similitude between equations is clearly visible by identifying k 2 → l(l + 1)/r 2 , k x f → mf /r 2 . We note that our definition of f in the spherical case differs from the planar case (RG02) in a factor r (f has units of magnetic field times length). In principle, one should solve the fully coupled system of equations for a large number of l, m's. However, as a first approximation, and in part motivated by the similitude of the equations with the cartesian case, we neglect coupling terms. We also note that, in RG02, the fastgrowing mode always has k y = 0. This cancels the term proportional to sf , the counterpart to our coupling C m l±1 . Our purpose is not to perform a complete and detailed linear analysis, but rather to asses on the interpretation of our 3D non-linear simulations.
In Fig. 3 we compare the results of this truncated linear analysis (right) to the non-linear simulations (left). The right panels show, in a red color scale, the inverse growth time of the unstable models for each l and m. Interestingly, the linear analysis (even neglecting couplings between modes) agrees very well with the simulations. For each multipole, the fastest growing modes always correspond to l = m. With some visual effort translating the two scenarios, one can compare our results in a spherical shell with Fig 1 in RG02 in a rectangular slab. Note that the correct analogy should compare l ≈ |k| = (k 2 x +k 2 y ) 1/2 and m ≈ k x . Therefore, the fastest growing modes with k y = 0 in RG02 should correspond to |k| = k x . This is l = m, as we obtain here. If we look for example at model A, the analytical estimate of the growth time is 1/Re(p max ) ≈ 0.6 kyrs, in excellent agreement with the observed growth in the left panel during the linear phase, t = 0, 0.5, and 1 kyrs. After a few growth times, the background has begun to change, the non-linear evolution sets in, and the full spectrum is filled by the Hall cascade (see curves at t = 1.75, 2 and 3 kyr), but we still find a significant excess power around the m = 20 region. We must note that the linear analysis is a first order approximation (because we omit couplings between neighbour modes), and one should not expect an exact identification of a single fast growing mode in the nonlinear results. Moreover, there are a few modes (in the m = 10 − 20 range) with very similar growth times.
The linear growth phase in Fig. 3a corresponds to the first 3 snapshots (t=0,0.5,1 kyrs), and the fact that there is a range of modes (not exactly peaking at m=18) is not due (yet) to a fast evolution of the background, but to the fact that the linear analysis is approximate (we truncate couplings with neighbour modes to make it simple). So, this numbers can be taken as a good indication, but we do not claim that the fast growing mode is exactly m=18, with t=1.75 kyrs. We can simply conclude that the most unstable modes are in the range m=10-25, with typical growth rates of 1-2 kyrs. Similar considerations apply to model B. We have also obtained a few eigenfunctions and checked that they are qualitatively similar to the eigenfunctions of RG02.
Thus, we confirm that the instability observed in the 3D simulations affects approximately the same wavelength range and grows on the same timescale of the linear analysis estimates. We should stress that the typical wavelength of the most unstable mode is closely correlated to the thickness of the shell where the toroidal field is confined. In a realistic crust, we expect structures with m ≈ 20 and a typical size of 2πR/20 ≈ 3 − 4 km (and proportionally smaller for toroidal rings shifted to higher latitudes). We have considered a purely toroidal field for simplicity of the analysis, but we have obtained very similar results in 3D simulations adding an initial poloidal component and a stratified density [25], concluding that the instability operating here is the Hall instability (RG02), rather than the ideal one.
A critical ingredient for the instabilities is the choice of the appropriate boundary conditions. Throughout this work, we have considered a non-permeating boundary condition at the inner crust and a vacuum potential solution in the outer region. A more realistic approach would consider the role of a magnetic field threading through the core, as at magnetar field strengths the assumption of a Type-I superconductor may not hold. In the exterior, a current filled-magnetosphere relaxing to a forcefree equilibrium or even dynamically evolving may be more suitable.
As our main purpose here is to investigate the development of the instability, we have chosen highly axisymmetric initial conditions, where the energy in the nonaxisymmetric component is six orders of magnitude less than the axisymmetric part. Because of that, we see the formation of multiple zones. In a realistic configuration, the initial conditions may not be that symmetric, therefore instead of the excitation of a higher multipolar structure, a less ordered magnetic field configuration may develop.
The main implication of our result is that a sufficiently strong toroidal field, as most magnetized NS models assume, is subject to this non-axisymmetric instability, and will break into small structures (typically [10][11][12][13][14][15][16][17][18][19][20] in the azimuthal direction. Such structures can occasionally make the crust yield to the magnetic stresses [2,28,29], leading to the formation of magnetic loops similar to the solar coronal loops. This mechanism is believed to be at the origin of magnetar outbursts. We also note that the loops created with this mechanism have magnetic field strengths typically one order of magnitude larger than the dipolar large scale field, in line with the observations [5]. Besides, our findings also have implications for the quiescent emission. It is has been proposed [30] that the high temperatures of magnetars are due to the dissipation of currents in a shallow layer when magnetospheric currents return to close the circuit inside the star (see similar arguments in [31,32]). The creation of small, force-free magnetic spots with the right size (a few km 2 ) is consistent with the typical sizes of the hot emitting spots of magnetars in quiescence [33,34]. Further works studying the coupled 3D mag-netic and thermal evolution of magnetars are needed to understand when, and how often, one of this spots results in a coronal-like flare and locally high temperatures.