---
abstract: 'We present three-dimensional force-free electrodynamics simulations of magnetar magnetospheres that demonstrate the instability of certain degenerate, high energy equilibrium solutions of the Grad-Shafranov equation. This result indicates the existence of an unstable branch of twisted magnetospheric solutions and allows to formulate an instability criterion. The rearrangement of magnetic field lines as a consequence of this instability triggers the dissipation of up to 30% of the magnetospheric energy on a thin layer above the magnetar surface. During this process, we predict an increase of the mechanical stresses onto the stellar crust, which can potentially result in a global mechanical failure of a significant fraction of it. We find that the estimated energy release and the emission properties are compatible with the observed giant flare events. The newly identified instability is a candidate for recurrent energy dissipation, which could explain part of the phenomenology observed in magnetars.'
author:
- |
    J. F. Mahlmann$^1$[^1], T. Akgün$^2$, J. A. Pons$^2$, M.A. Aloy$^1$ and P. Cerdá-Durán$^1$[^2]\
    $^1$Departament d’Astronomia i Astrofísica, Universitat de València, 46100, Burjassot, Spain\
    $^2$Departament de Fisica Aplicada, Universitat d’Alacant, 03690, Alicante, Spain
bibliography:
- 'literature.bib'
date: 'Accepted 2019 September 25. Received 2019 September 22; in original form 2019 July 31'
title: Instability of twisted magnetar magnetospheres
---

\[firstpage\]

stars: magnetars – magnetic fields – methods: numerical – stars: neutron – X-rays: bursts

Introduction {#sec:introduction}
============

Soft gamma-ray repeaters (SGRs) are neutron stars with recurrent X-ray activity in the form of short bursts with duration $\sim 0.1$s and luminosities in the range $10^{36}-10^{43}$ergs$^{-1}$. Over the last 40 years, three bursts have been uniquely energetic, the so-called [*giant flares*]{} (GFs) with luminosities of the order of $10^{44}-10^{47}$ergs$^{-1}$ [SGR0525-66, SGR1900+14 and SGR1806-20, see @Cline1980; @Hurley1999; @Hurley2005]. In the three referenced cases, a short initial peak was followed by a softer X-ray tail lasting for $50-400$s. The engine behind these extraordinary events are [*magnetars*]{}, neutron stars with the strongest known magnetic fields [$10^{14}-10^{16}$G; see comprehensive reviews of magnetar observations and physics, e.g. in @Woods2006; @Rea2011; @Turolla2015; @Mereghetti2015; @Kaspi2017].

The precise mechanism producing such energetic events is still unclear. Strong magnetic fields are a gigantic energy reservoir in magnetars, generally of the order $$\begin{aligned}
{\mathcal{E}}_{\rm magnetar} \sim 1.6\times 10^{47} \,\text{erg}\,\left ( \frac{B}{10^{15}\,\text{G}}\right )^{2}\left ( \frac{R_*}{10\,\text{km}}\right )^3,\end{aligned}$$ where we consider a neutron star with radius $R_*$.

The timescale on which the magnetar is evolving, mainly due to Hall drift and Ohmic dissipation in the crust, is of the order of $10^3-10^6$yrs [@Jones1988; @Goldreich1992; @Pons2007; @Pons2009; @Gourgouliatos2016], by itself too slow to explain this phenomenology. Two complementary models have tried to explain these observations. In the [*crustquake*]{} model [@Thompson1996; @Perna2011] the dynamical trigger is the mechanical failure of patches of the magnetar crust due to large stresses built during its magneto-thermal evolution. Numerical simulations of the Hall evolution of the crust [@Vigano2013] show that it is possible to recurrently reach the maximum stress supported by the very same [@Horowitz2009; @Baiko2018]. At this point, the crust likely becomes plastic [@Levin2012], i.e. the crust generates thermo-plastic waves emerging from such a localized trigger, or in other words *yields* [@Beloborodov2014; @Li2016]. The waves propagate into the magnetosphere, probably resulting in rapid dissipation through a turbulent cascade triggered by reconnection on slightly displaced flux surfaces [@Thompson1996; @Thompson2001; @Li2018]. The energy released in those events suffices to explain the observed luminosities, even for GFs [@Thompson1996; @Lander2015]. The burst duration ($\sim 0.1$s) is related to the crossing time of shear waves through the whole crust ($1-100$ ms). A limitation is that, if stressed for long periods of time ($\sim 1\,$yr) as it is the case due to the slow magneto-thermal evolution, the crust may yield at significantly lower breaking stresses [@Chugonov2010]. In that case, it would effectively deform as a plastic flow, and, depending on its (unknown) properties, cease to yield altogether [@Lyutikov2015; @Lander2019]. @Thompson2017 has argued that even in this case the crust could yield.

The [*magnetospheric instability*]{} model requires a strongly twisted magnetosphere that becomes unstable and leads to a rapid reconnection event [@Lyutikov2003]. The existence of long-lived magnetospheric twists is supported by the observation of hard X-ray emission in persistent magnetars [@Beloborodov2013b; @Hascoet2014]. During the magneto-thermal evolution of the crust, the displacement of the magnetic field footprints can generate large twists in the magnetosphere [@Akgun2017; @Akgun2018b]. Above a critical twist, the magnetosphere becomes unstable and undergoes a rapid rearrangement where energy is dissipated by reconnection [@Lyutikov2003; @Gill2010; @Elenbaas2016] in a similar fashion as in the crustquake model. The main challenge of this scenario is the ability of the crust to produce significant twists in the magnetosphere. @Beloborodov2009 estimated that currents supporting magnetospheric twist are bound to dissipate on timescales of years, effectively leading to a progressive untwisting. Therefore, Hall evolution is required to proceed relatively fast in order to allow for significant twists. Plastic viscosity may also be a problem for similar reasons [@Lander2019]. The latter authors have also suggested that the dynamical crust fractures of the crustquake model could be substituted by sustained episodes of accelerated plastic flows which are able to generate large magnetospheric twists on times shorter than the untwisting timescale.

Numerical simulations by @Parfrey2012 [@Parfrey2013] and @Carrasco2019 confirm the instability of the magnetosphere beyond a critical twist, accompanied by the formation of plasmoids. These results are an analogy to the context of eruption processes in the solar corona as found in numerical experiments by @Roumeliotis1994 [@Mikic1994]. The energy dissipated in the reconnection events is sufficient to explain the GF processes [@Parfrey2012]. A caveat to these simulations is that the applied twisting rate is larger than the one expected from the respective magneto-thermal evolution, although it would be fine if the trigger was a rapid plastic deformation.

An alternative approach to the above is the study of stability properties of magnetospheres. A number of authors have constructed equilibrium solutions to the Grad-Shafranov equation (GSE) for neutron star magnetospheres [@Glampedakis2014; @Fujisawa2014; @Pili2015; @Akgun2016; @Kojima2017; @Kojima2018; @Kojima2018b; @Akgun2018a]. @Akgun2017 performed magneto-thermal evolutions coupling the crustal magnetic field at the stellar surface with an exterior equilibrium solution. The results showed that large twists grow in the magnetosphere up to a critical point beyond which no stable equilibrium solutions where found. A more detailed analysis by @Akgun2018a showed that, for sufficiently large twists, the solutions of the GSE are degenerate with several possible configurations of different energies but matching boundary conditions at the surface. This suggests the possibility of an unstable branch of the solutions and, thus, a possible explanation for the occurrence of bursts and GFs. In this work we explore this possibility by performing three-dimensional (3D) numerical simulations of the equilibrium models in @Akgun2018a. We asses their stability properties and their potential as candidates for transient magnetar phenomenology.

This work is organised as follows. In section \[sec:physics\] we review and discuss the physics involved in magnetars relevant to the processes that we want to study. In section \[sec:ff\_equations\] we briefly review the equations of force-free electrodynamics (FFE) implemented for simulations conducted on the infrastructure of the `Einstein Toolkit` (supplemented by appendix \[sec:augmented\_system\]). A detailed description of the derivation of initial models according to @Akgun2018a is given in section \[sec:twisted\_models\]. In section \[sec:simulations\] we present the numerical setup of our simulations as well as the outcome of the conducted 3D force-free simulations of twisted magnetospheres (reviewing details on maintaining the force-free regime in appendix \[sec:ff\_breakdown\]). The observed rapid dissipation of electromagnetic energy through the magnetar crust is interpreted and related to observable quantities, such as luminosity estimates, shear stresses on the stellar crust, and opacity models, in section \[sec:discussion\]. Along this paper we use Gaussian units in CGS, except for section \[sec:ff\_equations\] in which we use Heaviside-Lorentz with geometrised units ($G=c=M_\odot=1$). For convenience we express current densities in A m$^{-2}$ and voltages in V, instead of the corresponding CGS units.

Physics of magnetars {#sec:physics}
====================

The basic structure of the magnetar interior is a (likely) fluid core of $\sim 10$ km radius, amounting for most of the mass of the object, surrounded by a solid crust of about $1$ km size. Outside, there is a tenuous, co-rotating magnetosphere connected to the NS by magnetic field lines (threading the central object) that extend up to the light cylinder at distances larger than $10^5$ km. We start by discussing some basic properties of the different parts of the magnetosphere relevant for the interpretations and models presented later in this work.

Currents supporting the magnetosphere {#sec:current_magnetsophere}
-------------------------------------

For the typical rotation periods of magnetars ($P\sim 1$ - $10$s) the Goldreich-Julian particle density [@GJ1969] for a magnetar magnetosphere has the typical value $$\begin{aligned}
n_{\textsc{gj}} = 7\times 10^{12} \,\text{cm}^{-3}\,\left ( \frac{B_{\rm pole}}{10^{15} \,{\rm G}}\right) \left( \frac{R_{*}}{r}\right )^3 \left ( \frac{10 \,{\rm s}}{P}\right),\end{aligned}$$ where $B_{\rm pole}$ is the magnetic field strength at the magnetar pole, $R_*$ the magnetar radius and $r$ the distance to the center of the star. This limits the magnetospheric current density close to the surface to $J<e\,c\,n_{\textsc{gj}} \approx 3\times 10^8$A/m$^2$, much below the typical values needed to support currents in strongly twisted magnetospheres of magnetars, of the order of $$J\sim \frac{B c}{4 \pi r}\sim 8.2\times 10^{12}\,\text{A}\,\text{m}^{-2}\,\left( \frac{B_{\rm pole}}{10^{15}\,\text{G}}\right ) \left( \frac{R_*}{10\,\text{km}}\right )^{-1}\label{eq:GJcurrent}.$$ In general, magnetospheric currents in magnetars cannot be supported neither by Goldreich-Julian charges nor by charges lifted from the surface. @Beloborodov2007 proposed that the currents are supported by $e^+$-$e^-$ pairs generated in the magnetosphere in an intermittent discharge process that can be sustained for voltages along magnetic field lines of about $10^8 -10^9$V. This voltage can be maintained by self-induction in untwisting magnetospheres [@Beloborodov2009]. This untwisting is driven by the effective resistivity of the magnetosphere; the thermal photons from the magnetar’s surface scatter resonantly off the charges supporting the magnetospheric currents, taking energy away, at the same time that pairs are produced. The untwisting timescale is $\sim 1$yr, and it may explain the spectral evolution of some magnetars [@Beloborodov2009].

Timescales {#sec:timecales}
----------

Changes in magnetars take place during two different timescales. On the one hand, there is a [*secular timescale*]{} of thousands of years during which the magnetar is essentially in equilibrium. On the other hand, there is a [*dynamical timescale*]{} associated to energetic events (burst, flares) that can produce observable variations on timescales as fast as $0.1$s. The latter are likely associated to out-of-equilbrium states.

### Secular timescales {#sec:secular}

The [*secular timescale*]{} is set by the slow magnetothermal evolution of the cooling object. The interior magnetic field evolution is dominated by Hall drift and Ohmic diffusion at the crust [see, e.g. @Vigano2012; @Fujisawa2014 and references therein], which proceeds on typical timescales of $10^3-10^6$ yr. The long-term evolution of the magnetosphere is driven by the changes in the crustal magnetic field, which displaces the footprints of the magnetospheric magnetic field lines. Since this evolution is much slower than the dynamical timescale of the magnetosphere (see below), it can be considered that the magnetosphere evolves through a series of equilibrium states. This evolution creates a twist in the magnetosphere supported by currents - until a critical maximum twist is reached ($\varphi_{\rm crit} \sim 1$rad) beyond which no magnetospheric equilibrium solutions exist [@Akgun2017]. The stability of the magnetosphere close to this critical point is the subject of this paper.

At the same time as the crustal magnetic field evolves, other processes in the magnetosphere can also contribute to the evolution. The untwisting of the magnetosphere on timescales of $\sim1$yr [@Beloborodov2009 and discussion in section \[sec:secular\]], may be a competing action to the twisting process described above.

Although the velocity of the footprints is typically very slow, numerical simulations of the magnetothermal evolution of magnetars including the magnetosphere show that, close to the critical point, it can be as fast as $v_{\varphi}\sim 1$kmyr$^{-1}$ [see @Akgun2017] in the most optimistic scenario. Therefore, close to the critical twist, the magnetosphere twists slowly ($\dot\varphi_{\rm max,crit} \,\lesssim\, 0.1$radyr$^{-1}$), evolving on timescales $\,\gtrsim\, 10$yrs. In the best case scenario, this timescale is comparable to the untwisting timescale ($\sim1$yr) and, hence, parts of the magnetosphere could sustain a significant twist. This timescale is still much longer than the dynamical timescale of the system (see below). Therefore, in our study of the dynamical behavior we can neglect the secular evolution of the field.

### Dynamical timescales {#sec:dynamical}

The [*dynamical timescale*]{} is set by the travel time of waves propagating in the different regions of the magnetar. In the magnetosphere, the mass density can be neglected in view of the dominating magnetic field energy density. Also, the velocity of Alfvén and fast magnetosonic waves is degenerated to the speed of light. Hence, within $\sim 100$km the whole magnetosphere is coupled through timescales smaller than $1$ms, which sets the scale for the dynamical evolution of the magnetosphere. In this region it is possible to neglect the inertia of the fluid in the evolution equations of so-called FFE, which is used in the numerical simulations of this work.

In the outermost parts of the crust, the force-free condition still holds because of low densities. At sufficiently high densities, elastic forces of the solid crust and pressure gradients break this condition. To estimate the transition density one may consider the depth at which waves propagate at a velocity significantly different to the speed of light. Two possible waves can travel in the interior of the magnetised crust, the so-called magnetosonic (ms) waves, related to sound waves, and magneto-elastic (me) waves, a combination of Alfvén and shear waves. The complete eigenvalue structure of relativistic ideal MHD equations in the presence of an elastic solid is not known. To make a simple order of magnitude estimate of the different wave speeds, we use the expression of magneto-elastic torsional waves parallel to the magnetic field derived in [@Gabler2012] as well as the expression for fast magnetosonic waves perpendicular to the field[^3]: $$\begin{aligned}
v_{\rm me} /c = \sqrt{\frac{\mu_s+B^2}{e + B^2}} \qquad \qquad v_{\rm ms} /c= \sqrt{\frac{e c_s^2 + B^2}{e+B^2}},\end{aligned}$$ where $e$ is the energy density and $\mu_s$ the shear modulus. Note that in the limit of low magnetic field ($B^2\ll \mu_s, \, B^2\ll e$) we recover the shear and sound speed, respectively. In the high magnetic field limit ($B^2\gg \mu_s,\, B^2\gg e$) both, $v_{\rm me}$ and $v_{\rm ms}$, coincide with the speed of light. Inside the fluid core ($\mu_s=0$) the magneto-elastic speed becomes the Alfvén speed.

Figure \[fig:speeds\] shows the value of the characteristic speeds in the outer layers of a typical NS model for different magnetic fields in the magnetar range. Indeed, both fast magnetosonic waves and Alfvén waves have a degenerate speed equal to the speed of light in the magnetosphere. Inside the outer crust ($\rho < 4\times 10^{11}$gcm$^{-3}$), all characteristic speeds transition from the speed of light to a significantly lower value, in a region that can still be considered force-free. This transition depends on the magnetic field strength, happening deeper inside the star for larger values of $B_{\rm pole}$. Given these characteristic speeds, any global rearrangement of the magnetosphere can modify the entire structure of the crust (of size $\sim 2 \pi R_*$) on a timescale of $\sim 1$ ms for magnetosonic waves and $\sim 10$ms for magnetoelastic waves.

One last aspect to consider is the ability of magnetospheric waves to transmit energy into the crust. The discussion should be limited to Alfvén waves, which become magnetoelastic waves once they penetrate the crust; the energy carried by fast magnetosonic waves in the magnetosphere can be neglected due to the small density, which renders the compressibility effects of fast-magnetosonic waves unimportant.

Since the characteristic time in the magnetosphere is $\sim 1$ms, the typical frequency of the waves generated during its dynamics is in the kHz range. At this frequency, the crust can be considered as a thin layer because its thickness ($\sim 1$km) is much smaller than the typical wavelength in the magnetosphere ($\lambda \sim 100\,$km). In this case the energy transmission coefficient for waves perpendicular to the surface is approximately [cf. @Link2014; @Li2015] $$\mathcal{T} =\frac{4 v_{\rm me}/c}{(1+v_{\rm me}/c )^2}\approx 0.04 \left (\frac{v_{\rm me}/c}{0.01} \right), 
\label{eq:transmission}$$ for typical physical conditions in the magnetar crust. Given the low transmission coefficients of magnetospheric Alfvén waves hitting the crust as well as the differences on timescales between the crust and the magnetosphere (typically $\sim 10$ times shorter in the later) it is reasonable to consider that most of the crust remains rigid during any dynamical rearrangement of the magnetosphere.

In our magnetar model we will consider two regions: A force-free region ([*exterior*]{}, hereafter) consisting of the magnetosphere and the force-free part of the outer crust as well as the magnetar [*interior*]{} for the remainder of the NS, which we will consider to be fixed during our simulations. The limit between both regions is a spherical surface below the NS surface, where magnetic field lines are anchored, and is located below the transition density between inner and outer crust at a density $\rho<4\times 10^{11}\,$gcm$^{-3}$. For the purpose of describing the simulations we will refer to this transition point simply as *surface*.

![Fast magnetosonic (solid lines) and magnetoelastic (dashed lines) speed in the outer layers of a magnetar, for different magnetic field strengths ranging from 0 to $10^{16}$G. The neutron star model corresponds to the $1.4 M_\odot$ mass APR+DH model of @Gabler2012. The magnetic field is considered to be constant for simplicity. []{data-label="fig:speeds"}](images/figure1.pdf){width="50.00000%"}

Force-free Electrodynamics {#sec:ff_equations}
==========================

In analogy to @Komissarov2004 and @Parfrey2017 we solve Maxwell’s equations in the force-free limit: $$\begin{aligned}
\frac{\partial {\boldsymbol{{\tilde{B}}}}}{\partial t} = -{\boldsymbol{\nabla}}\times{\boldsymbol{{\tilde{E}}}}\qquad\text{and}\qquad\frac{\partial {\boldsymbol{{\tilde{E}}}}}{\partial t} ={\boldsymbol{\nabla}}\times{\boldsymbol{{\tilde{B}}}} - {\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}}\:,
\label{eq:maxwell}\end{aligned}$$ where ${\boldsymbol{{\tilde{E}}}}$, ${\boldsymbol{{\tilde{B}}}}$, and ${\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}}$ are the electric field, the magnetic field, and the so-called force-free current, respectively. We place a tilde to distinguish quantities expressed in our Heaviside-Lorentz *geometrised* (HLG) system of units, while the same symbols without tilde express quantities in the Gaussian *non-geometrised* (GNG) system of units (see Table \[tab:unit\_conversions\]). We explicitly include the charge conservation equation $$\begin{aligned}
 \frac{\partial {\tilde{\rho}_{\rm e}}}{\partial t}+{\boldsymbol{\nabla}}\cdot{\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}} = 0\,,
 \label{eq:chargecons}\end{aligned}$$ where ${\tilde{\rho}_{\rm e}}$ represents the charge density. Furthermore, we use a mixed hyperbolic/parabolic correction by the introduction of additional potentials (further discussed in appendix \[sec:augmented\_system\]) in order to numerically ensure the constraints ${\boldsymbol{\nabla}}\cdot{\boldsymbol{{\tilde{B}}}}=0$ and ${\boldsymbol{\nabla}}\cdot{\boldsymbol{{\tilde{E}}}}={\tilde{\rho}_{\rm e}}$ [@Dedner2002; @Palenzuela2009; @Mignone2010].

In the force-free limit it is necessary to guarantee that there are either no forces acting on the system or, more generally, that the forces of the system balance each other. This is equivalent to a vanishing net Lorentz force on the charges ${\tilde{\rho}_{\rm e}}$ [see, e.g. @Camenzind2007]: $$\begin{aligned}
{\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}}&=0\label{eq:electromagnetic_work}\\
{\tilde{\rho}_{\rm e}}{\boldsymbol{{\tilde{E}}}}+{\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}}\times{\boldsymbol{{\tilde{B}}}}&=0\label{eq_lorentz_force}\end{aligned}$$ From equation (\[eq\_lorentz\_force\]) one readily obtains the degeneracy condition $$\begin{aligned}
{\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{{\tilde{B}}}}&=0.\label{eq:force_free_crossfield}\end{aligned}$$ Additionally, force-free fields are required to be magnetically dominant, the magnetic field being always stronger than the electric one, such that the following condition must hold: $$\begin{aligned}
{\boldsymbol{{\tilde{B}}}}^2-{\boldsymbol{{\tilde{E}}}}^2&\geq 0.\label{eq:force_free_dominance}\end{aligned}$$ Conditions (\[eq:force\_free\_crossfield\]) and (\[eq:force\_free\_dominance\]), as well as the conservation condition $\partial_t\left({\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{{\tilde{B}}}}\right)=0$ can be combined in order to obtain an explicit expression for ${\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}}$ [cf. @Komissarov2011; @Parfrey2017]: $$\begin{aligned}
{\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}}=\left[{\boldsymbol{{\tilde{B}}}}\cdot{\boldsymbol{\nabla}}\times{\boldsymbol{{\tilde{B}}}}-{\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{\nabla}}\times{\boldsymbol{{\tilde{E}}}}\right]\frac{{\boldsymbol{{\tilde{B}}}}}{{\tilde{B}}^2}+{\tilde{\rho}_{\rm e}}\frac{{\boldsymbol{{\tilde{E}}}}\times{\boldsymbol{{\tilde{B}}}}}{{\tilde{B}}^2}\label{eq:ff_current}\end{aligned}$$ Across the literature [e.g. @Komissarov2004; @Alic2012; @Parfrey2017] we find various modifications in the definition of ${\boldsymbol{{\tilde{J}}}}_{\rm \text{\tiny FF}}$ in order to drive the numerical solution of the system of partial differential equations (\[eq:maxwell\]) towards a state which fulfills equation (\[eq:force\_free\_crossfield\]) by introducing a suitable cross-field conductivity. In the numerical setup, we choose to combine the prescription of @Komissarov2004 with the force-free current given above. This strategy effectively minimises the violations of equations(\[eq:force\_free\_crossfield\]) and (\[eq:force\_free\_dominance\]) by exponentially damping the (numerically induced) components of the electric field parallel to ${\boldsymbol{{\tilde{B}}}}$ and suitably adjusting the electric field in magnetospheric current sheets in order to obtain ${\boldsymbol{{\tilde{B}}}}^2-{\boldsymbol{{\tilde{E}}}}^2\rightarrow 0$ at these locations.

Throughout the literature, the magnetic dominance condition (\[eq:force\_free\_dominance\]) condensates to a necessary condition of FFE [e.g. @Uchida1997; @McKinney2006]. For some authors [e.g. @McKinney2006] the breakdown of the magnetic dominance implies the invalidity of the numerical model. Others [e.g. @Uchida1997] claim that some physical processes (e.g. radiation losses) taking place in the regions where condition (\[eq:force\_free\_dominance\]) is breached may restore the magnetic dominance condition. Indeed, [@Uchida1997] explicitly allows for transient phases violating condition (\[eq:force\_free\_dominance\]) - these regions are then interpreted as abandoning the freezing of magnetic flux onto the flux of matter, being necessarily accompanied by dissipation. Following @Uchida1997, the force-free regime continues to be a valid approximation as long as the dissipative effects are only a small fraction of the total energy. The violation of the perpendicularity condition (\[eq:force\_free\_crossfield\]) is an additional source of (Ohmic) dissipation [studied for example in the context of Alfvén waves in force-free electrodynamics by @Li2018]. In practice, this channel of dissipation occurs when ${\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{{\tilde{B}}}}\neq 0$ such that ${\boldsymbol{{\tilde{J}}}}\cdot{\boldsymbol{{\tilde{E}}}}\neq 0$. Currently used force-free codes aim to avoid the transient into this regime by numerically cutting back all violations of condition (\[eq:force\_free\_dominance\]) [e.g. @Palenzuela2010; @Paschalidis2013; @Carrasco2016] or include a suitable Ohm’s law [e.g. @Komissarov2004; @Spitkovsky2006; @Alic2012; @Parfrey2017] in order to minimise these violations during a transient phase. Figure \[fig:energy\_dissipation\_channels\] shows the breakdown of condition (\[eq:force\_free\_dominance\]) during the simulation and hints towards the aforementioned dissipative processes. We refer to appendix \[sec:ff\_breakdown\] as well as, for example, @Lyutikov2003 for further details on the necessary constraint preservation and limitations of the highly magnetised regime (such as the lack of physical reconnection). We will give a thorough review of the procedures employed in our code in a subsequent technical paper.

Twisted magnetar magnetosphere models {#sec:twisted_models}
=====================================

Magnetospheres {#sec:initial_data}
--------------

      Quantity                   Nongeometrised unit                            Conversion factor
  ----------------- ---------------------------------------------- -------------------------------------------
        Mass                              M                                         $M_\odot$
       Length                             L                                    $M_\odot G c^{-2}$
        Time                              T                                    $M_\odot G c^{-3}$
   Electric charge   $\text{L}^{3/2}\text{M}^{1/2}\text{T}^{-1}$         $(4\pi)^{-1/2}M_\odot G^{1/2}$
   Electric field    $\text{L}^{-1/2}\text{M}^{1/2}\text{T}^{-1}$   $(4\pi)^{1/2}M_\odot^{-1} G^{-3/2}c^{4}$
   Magnetic field    $\text{L}^{-1/2}\text{M}^{1/2}\text{T}^{-1}$   $(4\pi)^{1/2}M_\odot^{-1} G^{-3/2}c^{4}$
   Current density   $\text{L}^{-1/2}\text{M}^{1/2}\text{T}^{-2}$   $(4\pi)^{-1/2}M_\odot^{-2} G^{-5/2}c^{7}$
     (EM) Energy        $\text{L}^{2}\text{M}\:\text{T}^{-2}$                   $ M_\odot\:c^{2}$
     (EM) Stress        $\text{L}^{-1}\text{M}\:\text{T}^{-2}$             $ M_\odot^{-2} G^{-3}c^{8}$

  : Conversion table between code output in Heaviside-Lorentz *geometrised* units ($M_\odot=G=c=1$) and *non-geometrised* Gaussian units. In order to transform the respective quantities from code quantities to the *non-geometrised* system, one has to multiply the *geometrised* quantity by its conversion factor expressed in CGS.[]{data-label="tab:unit_conversions"}

![Magnetospheric energy normalised to the vacuum dipole energy (equation \[eq:edipole\]) of the initial equilibrium models, for different values of the parameter ${P_{\rm c}}$ (in units of $P$ at the equator). The solid and dashed lines correspond to a series of models with constant $s$ and $\sigma$. The colored dots correspond to the initial data models used in our simulations.[]{data-label="fig:initial_model"}](images/figure2.pdf){width="45.00000%"}

Due to the long rotational period of observed magnetars pushing the location of the light cylinder to great distances, it is possible to neglect the rotation of the neutron star when building numerical models of magnetospheres in the near zone. The equilibrium structure of a non-rotating axisymmetric force-free magnetosphere is given through the well-known GSE [@Luest1954; @Grad1958; @Shafranov1966]. This approach has been followed in several recent papers [e.g. @Spitkovsky2006; @Beskin2010; @Vigano2011; @Glampedakis2014; @Fujisawa2014; @Pili2015; @Akgun2016; @Akgun2018a; @Kojima2017; @Kojima2018; @Kojima2018b]. In most of these works, the toroidal field is confined within a magnetic surface near the equator, smoothly transitioning to vacuum at large distances. In stationary, non-rotating, axisymmetric magnetosphere models, the toroidal field cannot extend to the poles. Otherwise, the toroidal field would extend all the way to infinity, thus, violating the requirements of finite magnetic energy. Following the notation of @Akgun2016 [@Akgun2018a], we write the axisymmetric magnetic field in terms of its poloidal and toroidal components: $$\begin{aligned}
 {\boldsymbol{B}} = {\boldsymbol{\nabla}} P \times {\boldsymbol{\nabla}}\varphi + T {\boldsymbol{\nabla}}\varphi \ ,\end{aligned}$$ where $\varphi$ is the azimuthal angle in spherical coordinates. Here, $P$ and $T$ are the poloidal and toroidal stream functions. Expressed in the orthonormal spherical basis corresponding to the coordinates $(r,\theta,\varphi)$, the magnetic field can be explicitly computed from the potentials $P$ and $T$ as $$\begin{aligned}
B^r &= \frac{1}{r^2\sin\theta}\partial_\theta P, \label{eq:Br}\\ 
B^\theta &= - \frac{1}{r \sin\theta} \partial_r P, \label{eq:Btheta}\\
B^\varphi &= \frac{T}{r \sin\theta}. \label{eq:Bphi}\end{aligned}$$ For an axially symmetric force-free field, the functions $T$ and $P$ may be expressed in terms of each other and appear as solutions of the force-free GSE: $$\begin{aligned}
\left[\partial_r^2 + \frac{1 - \mu^2}{r^2}\partial_\mu^2\right] P + T \frac{dT}{dP} = 0,\end{aligned}$$ where $\mu=\cos\theta$. $P$ and $T$ are constant on magnetic surfaces or, equivalently, along magnetic field lines. $P$ is related to the magnetic flux passing through the area centered on the axis and delineated by the magnetic surface. Therefore, its value at the poles is zero and increases towards the equator. The function $T$ is related to the current passing through the same area. Its functional dependence on $P$ can be chosen freely (consistently with any continuity and convergence requirements, particularly for the currents), which is equivalent to setting boundary conditions for $T$ at the surface of the star. Here, we invoke the same functional form for $T(P)$ as in @Akgun2016 [@Akgun2018a]. Thus, the toroidal field is confined within some [*critical*]{} magnetic surface ($P = {P_{\rm c}}$), $$\begin{aligned}
T\left(P\right)=\left\{\begin{array}{ccc}
s\times (P-{P_{\rm c}})^\sigma& : & P \geqslant {P_{\rm c}}\\ 
0 & : & \text{else}
\end{array}\right.,
\label{eq:toroidal}\end{aligned}$$ $s$ being a parameter determining the relative strength of the toroidal field with respect to the poloidal field. In order to avoid divergences in the currents we must demand that the power index satisfies $\sigma \geqslant 1$. For a pure dipolar field, the poloidal stream function in the magnetosphere is $$\begin{aligned}
P = \frac{1}{2}B_{\rm pole}\frac{R_*^3}{r}\sin^2\theta,\end{aligned}$$ while the toroidal stream function is $T=0$ everywhere. We will consider the simplest cases where the boundary value of $P$ at the surface of the magnetar coincides with that of a dipolar field, and, therefore, the initial data is symmetric with respect to the equator. For different choices of the functional relation $T(P)$ given by equation (\[eq:toroidal\]) we solve the GSE and obtain a twisted magnetospheric initial model. We would like to note that all equations can be rescaled with $B_{\rm pole}$, hence, the results of our numerical simulations can be normalised to the field strength of interest.

The energy stored in the magnetosphere can be computed as a volume integral $$\begin{aligned}
 {\mathcal{E}}= \frac{1}{8\pi}\int ({\boldsymbol{B}}^2 + {\boldsymbol{E}}^2)\, {\rm d}V .
 \end{aligned}$$ For later reference and in order to normalize the energetic content of our models, we provide the energy stored in the magnetosphere of a pure dipolar magnetic field ($\vec{E}=0$, $B^r=B_{\rm pole}(R_*/r)^3\cos{\theta}$, $B^\theta=(B_{\rm pole}/2)(R_*/r)^3\sin{\theta}$, $B^\varphi=0$): $$\begin{aligned}
\begin{split}
 {\mathcal{E}}_{\rm d} = \frac{1}{12}B_{\rm pole}^2 R_*^3 =8.3\times 10^{46}\,\text{erg}\,\left(\frac{B_{\rm pole}}{10^{15}\,\text{G}}\right)^2\left(\frac{R_*}{10\,\text{km}}\right)^3.
\end{split}
 \label{eq:edipole}
 \end{aligned}$$ Once the surface value of $P$ and the functional relation $T(P)$ are defined, one can solve the GSE iteratively (as it is a non-linear equation), while imposing vacuum boundary conditions at large distances. We use the numerical code described in @Akgun2018a to build our initial data. Using this parametrisation, the boundary condition at the surface of the neutron star for the GSE (values of $P$ and $T$) is fully determined by four parameters $B_{\rm pole}$, $s$, ${P_{\rm c}}$ and $\sigma$. However, the solution of the GSE with this fixed boundary condition is not necessarily unique. @Akgun2018a showed that for sufficiently large magnetospheric twists, there exist degeneracies, i.e. different solutions of the GSE for the same boundary conditions (the same set of four parameters). These solutions differ in their energy, twist and the radial extent of the toroidal currents.

Table \[tab:model\_overview\] shows the parameters used to construct the initial data for our numerical simulations. Each of the series A, B and C of initial models were chosen to have identical parameters but [*different*]{} magnetospheric energies and, hence, represent degenerate magnetospheric models. We would like to point out that the value of ${P_{\rm c}}$ is only equal, within each series, up to the second significant digit, due to numerical reasons. Figure \[fig:initial\_model\] shows the energy of the initial models as a function of the parameter ${P_{\rm c}}$. Models within each spiral curve (constant $s$ and $\sigma$) and with the same value of ${P_{\rm c}}$ have identical boundary conditions but different energies. In the interpretation made by @Akgun2018a, the lower energy state for each series of degenerate models (i.e., A1, B1 and C1) corresponds to stable configurations, while high energy states (i.e., A2, B2, C2 and C3) may be unstable and would evolve towards the stable configuration releasing the respective energy difference. This instability is a possible scenario for the flare activity observed in magnetars.

The lowest energy solutions are the ones that are most similar to the vacuum solutions, with all field lines connected to the surface, while the higher energy solutions are more radially extended, and can contain disconnected field lines.

![image](images/figure3.pdf){width="80.00000%"}

Magnetar interior {#sec:interior}
-----------------

The initial models described above provide solutions only for the magnetosphere. For each possible magnetospheric model one can build infinite solutions to describe the neutron star interior. The magnetospheric (exterior) values of $P$ and $T$ determine the magnetic field ${\boldsymbol{B}}$ at the exterior side of the surface (equations \[eq:Br\] to \[eq:Bphi\]). To match this solution to the interior, one has to ensure the continuity of $B^r$ at the surface. This is valid if $P$ is continuous and, hence, $T$ and $B^\varphi$ are continuous as well. However, $B^\theta$ does not necessarily match continuously to the neutron star interior because current sheets (thin current-carrying layers across which the magnetic field changes either direction or magnitude) in the $\varphi$ direction may occur. Even if all components of ${\boldsymbol{B}}$ are continuous at the surface, the magnetic field structure in the interior depends completely on how currents are internally distributed.

In the astrophysical scenario we are considering, the magnetar reaches the initial state in which we start our numerical simulation after a slow magnetothermal evolution that proceeds in a long timescale compared to the dynamical timescales (cf. section \[sec:timecales\]) of the magnetosphere ($\sim 1$ms) or the crust ($\sim 10$ms). On such long timescales, any current close to the surface of the NS is expected to be dissipated by Ohmic diffusion. Therefore, we consider that initially all fields are continuous across the surface. We build our interior solution by extrapolating the exterior magnetic field towards the stellar interior across a number of grid cells as needed by the reconstruction algorithm used for the magnetospheric evolution in our simulations. Since the neutron star is basically a perfect conductor, the initial charge density and electric field in the interior (and the magnetosphere) are set to zero.

The surface values of $B^r$ and $B^\varphi$ are coincident for degenerate models (e.g. within the series C1, C2 and C3 in Figure \[fig:initial\_model\]) because $P$ and $T$ at the surface are identical. However, since $P$ and $T$ may have a different radial dependence outside of the magnetar, and $B^\theta$ depends on the radial derivative of $P$ (equation \[eq:Btheta\]), it is different for every model of the same series.

Simulations {#sec:simulations}
===========

We have performed numerical simulations of the neutron star magnetosphere using the initial models in Table \[tab:model\_overview\]. For all the simulations we employ our own implementation of a General Relativistic FFE code in the framework of the `Einstein Toolkit`[^4] [@Loeffler2012]. The `Einstein Toolkit` is an open-source software package utilizing the modularity of the `Cactus`[^5] code [@Goodale2002a] which enables the user to specify so-called `thorns` in order to set up customary simulations. There exist other code packages such as `GiRaFFE` [@Etienne2017], which integrate the equations of force-free electrodynamics employing an evolution scheme based on the Poynting flux as a conserved quantity [cf. @McKinney2006; @Paschalidis2013] rather than the electric field and its current sources [as formulated in, e.g. @Komissarov2004; @Parfrey2017]. The `Einstein Toolkit` employs units where $M_\odot=G=c=1$, which sets the respective time and length scales to be $1M_\odot\equiv 4.93\times 10^{-6}\text{ s}\equiv 1477.98\text{ m}$. This unit system is a variation of the so-called system of *geometrised units* [as introduced in appendix F of @Wald2010], with the additional normalisation of the mass to $1M_\odot$ (i.e. our HLG units, as introduced in section \[sec:ff\_equations\]). For easy reference, we provide a set of conversion factors for relevant physical quantities in Table \[tab:unit\_conversions\].

Numerical setup {#sec:numerical_setup}
---------------

All shown simulations are conducted on a 3D box with dimensions $\left[4741.12 M_\odot\times4741.12 M_\odot\times4741.12 M_\odot\right]$ with a grid spacing of $\Delta_{x,y,z}=74.08M_\odot$ on the coarsest grid level. For the chosen magnetar model of radius $R_*=9.26M_\odot$ ($\simeq {13.7\,\text{km}}$km) this corresponds to a $\left[512R_*\times 512R_*\times 512R_*\right]$ box with a grid spacing of $\Delta_{x,y,z}=8R_*$. For the low and high resolution tests we employ seven and eight additional levels of mesh refinement, each increasing the resolution by a factor of two and encompassing the central object, respectively. This means that the finest resolution of our models (close to the magnetar surface) are $\Delta_{x,y,z}^{\rm min}=0.0625\times R_* = 0.5787M_\odot$ and $\Delta_{x,y,z}^{\rm min}=0.03125\times R_*= 0.2894M_\odot$ for the low and high resolution models, or in other words 16 and 32 points per $R_*$, respectively. The initial data is evolved for a period of $t=1185.28M_\odot\simeq 5.84\,$ms, which is chosen to be well below the dynamical timescale of the magnetar crust, which can be considered as a fixed boundary (see section \[sec:timecales\]).

In order to ensure the conservation properties of the algorithm, it is critical to employ [*refluxing*]{} techniques correcting numerical fluxes across different levels of mesh refinement [see, e.g. @Collins2010]. Specifically, we make use of the thorn `Refluxing`[^6] in combination with a cell-centered refinement structure [cf. @Shibata2015]. We highlight the fact that employing the refluxing algorithm makes the numerical code $2-4$ times slower for the benefit of enforcing the conservation properties of the numerical method (specially of the charge). Refluxing also reduces the numerical instabilities, which tend to develop at mesh refinement boundaries.

In conservative schemes, numerical reconstruction algorithms [we employ an MP7 scheme, cf. @Suresh1997] derive inter-cell approximations of the conservative variables by making use of their values at several adjacent grid-points (for MP7, one requires seven points). As a result of the numerical coupling between the magnetosphere and the magnetar crust introduced by the inter-cell reconstruction at the stellar surface, the field dynamics induce a mismatch in the current flowing through the surface and effectively trigger a (numerical) flow of charges leaving or entering the domain. In order to avoid this artifact, we replace the reconstructed values of the radial current ${\tilde{J}}^r_{\rm FFE}$ at interfaces between the stellar interior and exterior by the cell-centered value in the stellar interior. This procedure ensures a conservation of magnetospheric charge.

The (3D) initial data is imported from the (2D) initial models (see se. \[sec:initial\_data\]) by bicubic spline interpolation. Throughout the numerical evolution, all quantities on grid-points inside of the magnetar radius are fixed to their initial values.

Instability onset and magnetospheric energy balance {#sec:instability_relax}
---------------------------------------------------

![Energy balance during the evolution of the high resolution model C2 (Table \[tab:model\_overview\]). *Top:* Comparison of the change in total magnetospheric energy, normalized to the energy of a magnetosphere equiped with a pure dipolar magnetic field, $\Delta {\mathcal{E}}/{\mathcal{E}}_{\rm d}$, as well as the Poynting flux through the magnetar surface. Up to a simulation time of $t~\sim 3.33\,$ms the energy change is dominated by Poynting flux onto the magnetar crust. *Middle:* Maximum violation of the ${\boldsymbol{{\tilde{B}}}}^2-{\boldsymbol{{\tilde{E}}}}^2\ge 0$ condition throughout the numerical grid. *Bottom:* Maximum violation of the ${\boldsymbol{{\tilde{B}}}}\cdot{\boldsymbol{{\tilde{E}}}}= 0$ constraint throughout the numerical grid. At the time of the breakdown of conditions (\[eq:force\_free\_crossfield\]) and (\[eq:force\_free\_dominance\]), the energy change is dominated by secondary (possibly numerical) effects. []{data-label="fig:energy_dissipation_channels"}](images/figure5.pdf){width="45.00000%"}

We have performed simulations with initial models in the low energy branch (A1, B1 and C1) and in the high energy branch (A2, B2, C2, C3). We observe a differentiated behavior in the evolution of the system depending on the class of initial model. For models in the low energy branch we find that the magnetosphere is stable and that the system remains essentially unchanged. The energy of the system remains constant throughout the simulation (see blue lines in Figure \[fig:absolute\_energy\_evolution\]), confirming the stability of these configurations, at least on dynamical timescales. This is specially true in the high resolution models, which exhibit a smaller numerical dissipation. The slightly larger numerical dissipation of the low resolution models explains the small drift in time with respect to the initial energy displayed by the blue dashed lines in Figure \[fig:absolute\_energy\_evolution\]. On the other hand, models in the high energy branch become unstable on a timescale of a few milliseconds and the magnetosphere changes its shape roughly at the same time as the energy of the magnetosphere decreases (see red and green lines in Figure \[fig:absolute\_energy\_evolution\]). This numerical experiment confirms the hypothesis of @Akgun2018a that, for degenerate initial models, only the lowest energy state is stable, and that all corresponding degenerate cases of high energy are unstable. In addition, we note that the lower energy states are closer to a purely dipolar magnetosphere, hence, the minimised circumference of the magnetic surfaces minimise the magnetospheric energy content [cf. @Thompson1996].

For configurations in the unstable branch, the onset of the instability proceeds earlier for lower numerical resolution. This is expected because a coarser grid contains larger numerical discretisation errors acting as a seed for the instability onset. However, the rapid drop in energy during the instability proceeds in a similar fashion for both numerical resolutions, indicating that the instability has a physical origin and is not a numerical artifact. In the case of the high energy initial model C2 we observe a rearrangement of the lobes of magnetic twist towards a dipolar structure (see Figure \[fig:high\_energy\_relaxation\]) prior to a significant drop of magnetospheric energy (by approximately $30\%$ of its initial value). During the phase of full validity of the force-free condition (see equation \[eq:force\_free\_dominance\]) the loss of magnetospheric energy is dominated by an outgoing Poynting flux at the innermost boundary (see Figure \[fig:energy\_dissipation\_channels\]). For our boundary condition it can be interpreted as the formation of a strong current on a thin layer below the surface, where energy can be efficiently dissipated.

Following @Parfrey2013 in the context of twisted magnetar fields and @Li2018 in a study of energy dissipation in collisions of force-free Alfvén waves, the onset of the (topological) relaxation is likely to be linked to Ohmic heating ${\boldsymbol{J}}\cdot{\boldsymbol{E}}\neq 0$, which occurs as a result of (minor) violations of the force-free condition (\[eq:force\_free\_crossfield\]), as can be seen in the bottom panel of Fig. \[fig:energy\_dissipation\_channels\] (note the much smaller scale of that panel compared to the middle one). We give a more detailed review of the treatment of these violations in our code and throughout the literature in appendix \[sec:ff\_breakdown\].

Surface currents and long-time evolution {#sec:surface_currents}
----------------------------------------

Following the initial instability and subsequent rapid rearrangement of the magnetar magnetosphere (section \[sec:instability\_relax\]), thin currents form at the magnetar surface (see Figures \[fig:surface\_currents\_xz\] and \[fig:surface\_currents\_avg\]). These currents are expected to appear as the initial model in the high energy state tries to relax to the lowest energy magnetospheric configuration, while keeping the interior field fixed (see the discussion in section \[sec:interior\]). There are two possible fates for these currents: i) they could propagate inwards, inside the magnetar crust, deforming the magnetic field inside and creating a mechanical stress in the crust, on a timescale of several $10 \text{ ms}$, or ii) they could form a thin surface current dissipating on a timescale shorter that the time it takes to deform the crust. These two possibilities are not mutually exclusive and a combination of both is possible. In none of the cases our simulations can give a conclusive answer because i) we are not evolving the magnetar interior as we are considering only timescales smaller than the dynamical timescale of the crust, ii) the formation of thin surface currents is numerically challenging (would require a computationally prohibitively high resolution near the magnetar surface), and iii) it would eventually violate the FF conditions (\[eq:force\_free\_crossfield\]) and (\[eq:force\_free\_dominance\]), hence invalidating our current numerical approach.

![*xz*-cross-sections of the toroidal current in geometrised units showing the development of strong surface currents during the evolution, in addition to other currents extended on larger magnetospheric volumes. *Top:* Low resolution model C2 (16 points per $R_*$). *Bottom:* High resolution model C2 (32 points per $R_*$). The high resolution evolution shows currents located around the magnetar surface with more detailed structures, emphasizing their interpretation as surface currents. The spatial coincidence of the currents in both resolutions reinforce the argument that the observed currents are of physical origin (in spite of the - relatively small - differences among different resolutions).[]{data-label="fig:surface_currents_xz"}](images/figure6.pdf){width="49.00000%"}

The aforementioned current layers are expected to be regions of strong energy dissipation and the breakdown of the force-free conditions [see, e.g. @Uchida1997; @McKinney2006; @Palenzuela2010; @Parfrey2013]. Figures \[fig:energy\_dissipation\_channels\] and \[fig:surface\_currents\_avg\] link the breakdown of the force-free condition (\[eq:force\_free\_dominance\]) and the occurrence of surface currents with the opening of dissipation channels different to the energy flow through the magnetar surface (see appendix \[sec:ff\_breakdown\] for a short review of the force-free breakdown). We find the violation of condition (\[eq:force\_free\_crossfield\]) to be continuously occurring with peaks at the instance of rapid energy dissipation. Condition (\[eq:force\_free\_dominance\]) starts to fail on longer timescales at the moment of fastest transfer of magnetic energy through the surface. At this time, further dissipation mechanisms (see Figure \[fig:energy\_dissipation\_channels\]) come into play, as is expected throughout the literature [@Uchida1997; @McKinney2006; @Li2018].

It should be noted that the total magnetospheric energy for the models B2, C2, and C3 drops below the energy of their respective low energy equilibrium solutions, and even below the magnetospheric energy of the vacuum dipole (equation \[eq:edipole\]). However, this energy drop is (slightly) smaller for the high resolution simulations, and shows some dependence on the chosen setup of the hyperbolic/parabolic cleaning procedures (see appendix \[sec:augmented\_system\]) at the magnetar surface. The sensitivity of this behavior to the numerical details at the location of the (3D Cartesian) crust may be attributed to the numerical dissipation of the employed code.

![Azimuthal angular averages of the toroidal current (normalised to its initial value at the stellar surface) in the equatorial plane showing the development of surface currents during the evolution of the C2 initial model. We display the current evolution for both low resolution (16 points per $R_*$, denoted by dashed lines), and high resolution (32 points per $R_*$, denoted by solid lines) models. The increase of the toroidal current during the transient of energy dissipation (see Figure \[fig:absolute\_energy\_evolution\]) at the lower resolution (compare the two blue lines) may be attributed to a faster onset of the twist instability for this model.[]{data-label="fig:surface_currents_avg"}](images/figure7.pdf){width="45.00000%"}

  -------------- -------------------------------- ---------------------------------------------- --------------------- -------------------------------- -------------------------------
      Model       $\Delta t_{\rm r} \text{(ms)}$   $\Delta {\mathcal{E}}/{\mathcal{E}}_{\rm d}$   $\Delta_{\rm mx} J$   $\Delta_{\rm mx} T^{r\varphi}$   $\Delta_{\rm mx} T^{r\theta}$
   \*\[0pt\] A1               5.8400                                  0.0033                            0.0159                      0.0012                          0.0010
        A2                    1.4162                                  0.0963                            1.6350                      0.0295                          0.0150
        B1                    5.8400                                  0.0042                            0.0363                      0.0012                          0.0014
        B2                    3.0427                                  0.1002                            0.9805                      0.0358                          0.0232
        C1                    5.8400                                  0.0009                            0.0640                      0.0008                          0.0013
        C2                    2.1604                                **0.2808**                        **3.5400**                    0.0851                          0.0414
        C3                    1.0490                                  0.1962                            3.1720                    **0.1008**                      **0.0811**
  -------------- -------------------------------- ---------------------------------------------- --------------------- -------------------------------- -------------------------------

  : Selection of electromagnetic quantities monitored throughout the (high resolution, 32 points per $R_*$) simulation of the models of Table \[tab:model\_overview\]. The total change in energy $\Delta {\mathcal{E}}$ (displayed as a fraction of the vacuum dipole energy; equation \[eq:edipole\]) corresponds to the maximum drop of electromagnetic energy during the total runtime (see section \[sec:energy\_release\]). The operator $\Delta_{\rm mx}$ acting on any quantity $A(t,{\boldsymbol{x}})$ is defined as $\Delta_{\rm mx}A:=\max_{\{t,|{\boldsymbol{x}}|=R_*\}}{ \{ A(t,{\boldsymbol{x}})-A(0,{\boldsymbol{x}}) \}/\max_{\{|{\boldsymbol{x}}|=R_*\}}A(0,{\boldsymbol{x}}) }$. Hence, $\Delta_{\rm mx} J$ is the maximum increase in current density in the magnetosphere during the relaxation relative to the initial values (see section \[sec:current\_magnetsophere\]). In the right columns, $\Delta_{\rm mx} T^{r\varphi}$ and $\Delta_{\rm mx} T^{r\theta}$ denote the maximum increase of electromagnetic stresses relative to their corresponding initial values (see section \[sec:crust\_stresses\]) on the stellar surface compared to its initial value. We highlight with bold face the maximum values of each of the last four columns.[]{data-label="tab:model_results"}

Discussion {#sec:discussion}
==========

Conclusions {#sec:conclusion}
===========

In this work, we explore the stability properties of force-free equilibrium configurations of magnetar magnetospheres by performing numerical simulations of a selection of the models computed in @Akgun2018a. For the case of degenerate magnetospheres (i.e. the same boundary conditions but different energies) we validate the hypothesis of @Akgun2018a that configurations in the high-energy branches are unstable while those in the lowest energy branch are stable. This confirms the existence of an unstable branch of twisted magnetospheres. It also allows to formulate an instability criterion for the sequences of models computed in @Akgun2018a. Our results are consistent with an interesting scenario where bursts and GFs in magnetars are triggered without involving crustal failures. The twist that is naturally produced in the magnetosphere by the Hall evolution of the crust [@Akgun2017] can lead to unstable configurations that will release up to a 10% of the energy stored in the magnetosphere, sufficient to explain the observations.

[@Akgun2017] have shown that the magneto-thermal evolution of the crust leads naturally to configurations close to the instability threshold. However, the amount of energy released depends on how far away from the stable branch can the evolution drive the configuration. This is essentially a problem of comparing the evolution timescale and the instability timescale. For the models studied in this work the instability timescale is of the order of milliseconds, much shorter that the magneto-thermal evolution timescales of the object (see Sect. \[sec:secular\]). However, close to the critical point, the growth rate of the instability could be significantly smaller (actually, it should be zero at the critical point) which would allow to overshoot the instability threshold. Note that, since the energy reservoir is large ($\sim 10^{46}$ erg), even a very small fraction of energy release could explain many of the phenomenology of magnetars. Alternatively, there could be phenomena leading to fast dynamics in the crust such as sustained episodes of accelerated plastic flows triggered by the magnetic stresses in the crust [@Lander2019].

For the unstable models, we observe the development of almost axisymmetric instabilities on a timescale of a few ms rearranging the magnetic field to a configuration similar to those in the (stable) lower energy branch. The energy of the magnetosphere also decreases towards the value of the stable configuration. Differences with respect to the corresponding stable configuration can be attributed to the influence of the non-preservation of the force-free constraints (\[eq:force\_free\_crossfield\]) and (\[eq:force\_free\_dominance\]). Using (much) larger numerical resolution (beyond the scope of our computational resources) we envision that the violation of the force-free constraints would be significantly reduced and the expected (low-energy) states would be the endpoint of the evolution after a full relaxation of the magnetosphere takes place. The energy decrease is explained, mainly, by a flow of energy towards the surface of the star, where it is dissipated efficiently. A large fraction of this energy is also dissipated in the magnetosphere at locations where the force-free conditions break. This contrasts with the work of @Beloborodov2011, @Parfrey2013 and [@Carrasco2019] in which most of the energy is dissipated by the formation and ejection of plasmoids. The different setup used in these workst (dynamically twisting vs. unstable equilibrium configurations) makes a direct comparison difficult. A possible source for the qualitative discrepancy may be differences in the boundary condition at the surface of the star. While we use a boundary condition that dissipates very efficiently any strong currents formed at the surface, in their work, their use of essentially non-dissipative boundary conditions make the surface perfectly reflective. For the future it would be interesting to compare more closely the differences in the boundary condition and to develop a better physical model for dissipation at the NS surface.

The magnetic field remains nearly axisymmetric throughout the simulation indicating that the instability is mostly an m=0 instability. A complete theoretical analysis of the origin of the instability and its properties is beyond the scope of this paper. However, we anticipate that such analysis has to be carried out on a global scale either by calculating the eigenmodes or by using the so-called energy principle of @Bernstein1958 and is not trivial due to the presence of both poloidal and toroidal components [@Akguen2013 and references therein]. However, we note that, since the poloidal field structure changes somewhat less than the toroidal field, this instability could be compared to the interchange instability discussed by @Tayler1973, where displacing the toroidal field radially decreases the energy (even in the absence of a fluid).

We have made a crude estimation of the observational properties of the energy liberated in the magnetosphere as a result of the instability. The fact that large amounts of energy (in excess of $10^{46}\,$erg) are released on milliseconds timescales results in dynamical luminosities significantly larger than $10^{48}\,$ergs$^{-1}$ (reaching in some models $4\times 10^{49}\,$ergs$^{-1}$). This should trigger the expansion of a pair-photon fireball polluted with baryons unbound from the magnetar crust. The bolometric signature of these fireballs seems incompatible with the observations of the initial spikes observed in GFs. With our simple analytic model, most of the unstable magnetospheres produce over-luminous, too cool and excessively short flashes. However, this problem can be solved if the energy can be liberated on longer timescales, of the order of the observed GF spikes ($\Delta t_{\rm spike}\sim 0.1\,$s). This could be possible in a scenario of slow energy dissipation as the one proposed by [@Li2018], which we plan to explore in the future.

The currents produced during the instability increase significantly the amount of pairs in the magnetosphere, a large fraction of which, of size $\sim 10 R_*$, becomes optically thick. The hot plasma magnetically confined in this region could be responsible for the extended thermal X-ray emission lasting for $50-300$s after GFs.

Our force-free numerical method cannot properly deal with the evolution of extremely thin surface currents. Therefore, the dynamical millisecond timescales computed in our models should be taken as a lower bound for the physical timescales. The magnetic dissipation taking place at these locations can be due to, e.g. Ohmic processes or to non-linear Alfvén wave interactions. Assuming that energy is released on $\sim\Delta t_{\rm spike}$, our estimate of the electromagnetic signature yields photospheric luminosities and temperatures compatible with observational data. Since this is a sound physical assumption, we conclude that observed GFs in SGRs are broadly compatible with the development of instabilities in twisted magnetospheres.

Acknowledgements
================

We thank Amir Levinson for his support in challenging our numerical code prior to the production of the presented results. We also thank Oscar Reula, Federico Carrasco and Carlos Palenzuela for their valuable discussions on the boundary conditions. We acknowledge the support from the grants AYA2015-66899-C2-1-P and PROMETEO-II-2014-069. JM acknowledges a Ph.D. grant of the *Studienstiftung des Deutschen Volkes*. PC acknowledges the Ramon y Cajal funding (RYC-2015-19074) supporting his research. We acknowledge the partial suport of the PHAROS COST Action CA16214 and GWverse COST Action CA16104. The shown numerical simulations have been conducted on *MareNostrum 4* of the *Red Española de Supercomputación* (AECT-2019-1-0004) as well as on the computational infrastructure of the *University of Valencia*. We thank the *EWASS 2019* and *GR22* conferences for the possibility to disseminate the results of this work.

Numerical details {#sec:num_details}
=================

The augmented system {#sec:augmented_system}
--------------------

![Energy evolution of the high energy initial data models A2, B2, and C2 using different damping constants $\kappa_\psi$ (divergence cleaning) in a low resolution study (16 points per $R_*$). While one observes a converging evolution for the lower cleaning potentials $\kappa_\psi=0.03125$ and $\kappa_\psi=0.125$, the energy evolution shows a strong (non-physical) dependence on $\kappa_\psi$ for larger damping constants. This effect is amplified in the high resolution (32 points per $R_*$).[]{data-label="fig:psi_optimization"}](images/figureA1.pdf){width="49.00000%"}

In order to preserve the physical conditions $\text{div}{\boldsymbol{{\tilde{B}}}}=0$ and $\text{div}{\boldsymbol{{\tilde{E}}}}={\tilde{\rho}_{\rm e}}$ we make use of hyperbolic/parabolic cleaning potentials [@Dedner2002; @Palenzuela2009; @Mignone2010]. Specifically, we implement an augmented system of Maxwell’s equations as follows [@Palenzuela2009; @Miranda-Aranguren2018]: $$\begin{aligned}
\partial_t\phi-\partial_i {\tilde{E}}^i&=-\tilde{\rho}_{\rm e}-\kappa_\phi\phi\label{eq:scalar_phi}\\
\partial_t {\tilde{E}}^i-\partial_j\left(\epsilon^{ijk}{\tilde{B}}_k+\delta^{ij}\phi\right)&=-{\tilde{J}}^i_{\rm \textsc{FF}}\label{eq:vector_phi}\\
\partial_t\psi+\partial_i {\tilde{B}}^i&=-\kappa_\psi\psi\label{eq:scalar_psi}\\
\partial_t {\tilde{B}}^i+\partial_j\left(\epsilon^{ijk}{\tilde{E}}_k+\delta^{ij}\psi\right)&=0\label{eq:vector_psi}\end{aligned}$$ Here, $\psi$ (divergence cleaning) and $\phi$ (charge conservation) are the scalar potentials, $\kappa_\phi$ and $\kappa_\psi$ the respective damping constants and $\delta^{ij}$ denotes the Kronecker delta. As for the practical implementation, we follow a Strang splitting approach [as employed, e.g. in @Komissarov2004], effectively solving part of the scalar equations (\[eq:scalar\_phi\]) and (\[eq:scalar\_psi\]) analytically. Prior (before `MoL_Step`) and after (before `MoL_PostStep`) the time integration of the Einstein Toolkit `thorn` `MoL` we evolve in time the equations $$\begin{aligned}
\phi\left(t\right)&=\phi_0\exp\left[-\kappa_\phi t\right],\\
\psi\left(t\right)&=\psi_0\exp\left[-\kappa_\psi t\right],\end{aligned}$$ for a time $t=\Delta t/2$. The coefficients $\kappa_\phi$ and $\kappa_\psi$ have to be chosen by optimisation in accordance with the grid properties.

We find it beneficial to choose a large value for $\kappa_\phi$, effectively dissipating charge conservation errors on very short timescales. As for the divergence cleaning, we conducted a series of tests, optimizing $\kappa_\psi$ to yield stable and converging evolution for all shown resolutions, ultimately resorting to $\kappa_\psi = 0.125$ (see Figure \[fig:psi\_optimization\] for a review of the optimisation process).

It should be noted at this point that @Mignone2010 present a promising scheme of choosing $\kappa_\psi$ according to the grid resolution that has also been used in [@Miranda-Aranguren2018]. In the framework of mesh refinement of the Einstein Toolkit, this would result in a different damping of the cleaning potentials across the refinement levels. We have found that the optimisation of the hyperbolic/parabolic cleaning becomes a very subtle issue and may experience strong numerical effects when increasing the overall resolution. This observation may, however, be an artifact of the fixed boundary of the magnetar surface - which on a Cartesian grid, resembles an accumulation of boxes rather than a perfectly aligned spherical boundary. The exploration of these effects and the transition to a fully spherical version of this force-free `thorn` [as introduced in @Baumgarte2013; @Montero2014] will be a subject of future efforts.

Conservation of force-free constraints {#sec:ff_breakdown}
--------------------------------------

FFE codes are valid in the limit of high electromagnetic energy compared to the rest mass and thermal energy of the respective plasma. The dynamics of force-free fields is described entirely without the plasma four-velocity. However, demanding the existence of a physical, timelike velocity field ${\boldsymbol{u}}$ with $F_{\mu\nu}u^\nu = 0$, as well as the degeneracy condition $F_{\mu\nu}J^\nu=0$ [see @Uchida1997 for a detailed algebraic review] one is left with the aforementioned constraints: $$\begin{aligned}
{\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{{\tilde{B}}}}&=0 \tag{\ref{eq:force_free_crossfield}}\\
{\boldsymbol{{\tilde{B}}}}^2-{\boldsymbol{{\tilde{E}}}}^2&\geq 0 \tag{\ref{eq:force_free_dominance}}\end{aligned}$$ Within the shown simulations we find it beneficial to employ an approach presented in @Komissarov2011 and @Parfrey2017 in order to archive $\partial_t\left({\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{{\tilde{B}}}}\right)=0$ throughout the evolution (by making use of the force-free current as in equation \[eq:ff\_current\]) without the employment of target currents [as discussed in @Parfrey2017]. Additionally, we include a suitable Ohm’s law [@Komissarov2004 section C3] into our Strang splitting approach aiming towards an evolution minimizing the violation of conditions (\[eq:force\_free\_crossfield\]), and (\[eq:force\_free\_dominance\]).

In order to build up a force-free current, @Komissarov2004 introduces a generalised Ohm’s law in the context of FFE: $$\begin{aligned}
\begin{split}
{\boldsymbol{{\tilde{J}}}}=\sigma_\parallel{\boldsymbol{{\tilde{E}}}}_\parallel+\sigma_\perp{\boldsymbol{{\tilde{E}}}}_\perp+{\boldsymbol{\tilde{j}}}_d\label{eq:GeneralizedOhmsLaw},
\end{split}\end{aligned}$$ where the subscripts $\parallel$ and $\perp$ denote the components parallel and perpendicular to the magnetic field, ${\boldsymbol{{\tilde{B}}}}$. A to be specified model for $\sigma$ introduces a suitable resistivity into the force-free system [see also @Lyutikov2003 for further comments on resistive FFE], while ${\boldsymbol{\tilde{j}}}_d$ is the drift current perpendicular to the electric and magnetic fields. In its general form, (\[eq:GeneralizedOhmsLaw\]) plays the central role in ensuring the force-free conditions (\[eq:force\_free\_crossfield\]) and (\[eq:force\_free\_dominance\]). @Komissarov2004 suggests a resistivity model that depends on the time-step of the evolution $\Delta t$ (throughout the presented simulations we employ CFL $=0.2$), where $$\begin{aligned}
\begin{split}
\sigma_\parallel = \frac{d}{\Delta t}\label{eq:sigmapar}.
\end{split}\end{aligned}$$ The cross-field resistivity $\sigma_\perp$ is strongly linked to the violation of condition (\[eq:force\_free\_dominance\]), $$\begin{aligned}
\arraycolsep=1.4pt\def{1.5}{1.5}
\begin{split}
\sigma_\perp=\left\{\begin{array}{ccc}
0 & \qquad :\qquad & {\boldsymbol{B}}^2\geq{\boldsymbol{E}}^2\\
b\displaystyle \frac{\left({\tilde{E}}_\perp-{\tilde{E}}_\perp^*\right)}{{\tilde{E}}_\perp^*} & \qquad :\qquad & {\boldsymbol{{\tilde{B}}}}^2<{\boldsymbol{{\tilde{E}}}}^2
\end{array}\right.\label{eq:sigmaperp},
\end{split}\end{aligned}$$ where ${\tilde{E}}_\perp=\left|{\boldsymbol{{\tilde{E}}}}_\perp\right|$ and $\left({\tilde{E}}_\perp^*\right)^2=\left({\boldsymbol{{\tilde{B}}}}-{\boldsymbol{{\tilde{E}}}}_\parallel\right)^2$ and $b$ is an scalar parameter controlling the magnitude of $\sigma_\perp$. Equations (\[eq:sigmapar\]) and (\[eq:sigmaperp\]) have a pair of analytic solutions: $$\begin{aligned}
\begin{split}
{\boldsymbol{{\tilde{E}}}}_\parallel\left(t\right)=\:&{\boldsymbol{{\tilde{E}}}}_\parallel\left(0\right)\times e^{-\sigma_\parallel t}
\end{split}\label{eq:StrangCurrentPar}\\
\begin{split}
{\boldsymbol{{\tilde{E}}}}_\perp\hspace{-3pt}\left(t\right)=\:&\left[{\tilde{E}}_\perp^*\hspace{-3pt}\left(0\right)+\frac{{\tilde{E}}_\perp^*\hspace{-3pt}\left(0\right)\left[{\tilde{E}}_\perp\hspace{-3pt}\left(0\right)-{\tilde{E}}_\perp^*\hspace{-3pt}\left(0\right)\right]\times e^{-b\sigma_\parallel t}}{{\tilde{E}}_\perp\hspace{-3pt}\left(0\right)-\left[{\tilde{E}}_\perp\hspace{-3pt}\left(0\right)-{\tilde{E}}_\perp^*\hspace{-3pt}\left(0\right)\right]\times e^{-b\sigma_\parallel t}}\right]\\
&\times\frac{{\boldsymbol{{\tilde{E}}}}_\perp\hspace{-3pt}\left(0\right)}{{\tilde{E}}_\perp\hspace{-3pt}\left(0\right)}.
\end{split}\label{eq:StrangCurrentPerp}\end{aligned}$$ During our numerical simulations), we usually choose $d=5.0$, and $b=0.1$, and solve equation (\[eq:StrangCurrentPar\]) prior to equation (\[eq:StrangCurrentPerp\]) in a Strang splitting scheme in direct analogy to the implementation described in section \[sec:augmented\_system\]. This resistivity model ensures the validity of the force-free regime throughout time, in other words, the evolution is driven towards a force-free state $$\begin{aligned}
\begin{split}
{\boldsymbol{{\tilde{E}}}}\cdot{\boldsymbol{{\tilde{B}}}}&\rightarrow0 \\
{\boldsymbol{{\tilde{B}}}}^2-{\boldsymbol{{\tilde{E}}}}^2&\rightarrow 0 \qquad\text{ : }{\boldsymbol{{\tilde{B}}}}^2<{\boldsymbol{{\tilde{E}}}}^2.
\end{split}\end{aligned}$$

Optical depth to resonant cyclotron scattering {#sec:optical_thickness}
==============================================

For the presented modeling of the optical thickness of highly magnetised force-free plasmas around magnetars (see section \[sec:emission\_processes\]), we adapt the techniques describing resonant scattering as presented by @Beloborodov2013 (from now on Be13). In the following, we will give a short review of the underlying equations. In order to derive the optical thickness $\tau$, we integrate equation (Be13/A15), $$\begin{aligned}
\frac{\text{d}\tau}{\text{d} s}=2\pi^2 r_e\frac{c}{\omega}\frac{\xi}{\left|\tilde{\mu}\right|}n_e\left[f_e\left(p_1\right)+f_e\left(p_2\right)\right].
\label{eq:optical_thickness_integrand}\end{aligned}$$ Here, $r_e=e^2/m_e c^2$ denotes the photon wavelength, $\omega$ the frequency of the seed photon (we consider 1keV photons), and $\xi=1$ or $\xi=\tilde{\mu}^2$ depending on the photon polarisation ($\perp$ or $\parallel$, respectively). The relativistic particles require the specification of the quantities $\mu=\cos\vartheta$ and $\tilde{\mu}=\cos\tilde{\vartheta}$, where $\vartheta$ is the angle between the photon path and the magnetic field ${\boldsymbol{B}}$ in the lab frame and $\tilde{\vartheta}$ in the rest frame of the electron. The dimensionless momenta $p_{1,2}$ correspond to the electron (or positron) velocities favored by the resonant scattering model. As both polarisations yield similar results, we only consider the slightly dominant $\perp$ orientation for our model. @Beloborodov2013 estimated that the contribution of non-resonant scattering to the optical depth is negligible and will not be considered in our calculations (see, however, footnote \[foot:one\]).

Following Be13, we employ the so-called waterbag model as a distribution function for electron (or positron) momenta. In analogy to a two-fluid model, the distribution function is characterised by the two parameters (dimensionless momenta) $p_+$ and $p_-$, with the overall shape $$\begin{aligned}
f_e\left(p\right)=\left\{\begin{array}{ccc}
(p_+ - p_-)^{-1}& : & p_-<p<p_+ \\ 
0 & : & \text{else}
\end{array}\right..
\label{eq:waterbag}\end{aligned}$$ Applying the waterbag model (\[eq:waterbag\]) in equation (\[eq:optical\_thickness\_integrand\]) selects the relevant electron (or positron) momenta for the scattering process. The distribution of this normalisation factor throughout the magnetosphere especially depends on the flow direction of charges along ${\boldsymbol{B}}$. As described in section 5.2 of Be13, we adjust their model according to a flow of electrons (or positrons) which turns back to the central object when field lines cross the equator. We apply this to all field lines crossing regions with $B<10^{13}$G (this holds everywhere except in the inner coronal region of strong closed magnetic field lines).

\[lastpage\]

[^1]: jens.mahlmann@uv.es

[^2]: pablo.cerda@uv.es

[^3]: Slow magnetosonic waves are also possible but their velocity is much smaller and not relevant for this work, in fact, for the case of waves perpendicular to the magnetic field their speed is zero.

[^4]: <http://www.einsteintoolkit.org>

[^5]: <http://www.cactuscode.org>

[^6]: Refluxing at mesh refinement interfaces by Erik Schnetter: <https://svn.cct.lsu.edu/repos/numrel/LSUThorns/Refluxing/trunk>
