Triggering magnetar outbursts in 3D force-free simulations

TL;DR

This study uses 3D force-free general relativity simulations to explore magnetar outbursts, revealing how magnetic twist leads to instabilities, reconnection, and periodic energy releases, with implications for magnetar surface heating.

Contribution

First 3D force-free simulations of magnetar magnetosphere dynamics, demonstrating the effects of magnetic helicity injection and the resulting instabilities and energy releases.

Findings

Twisting beyond a critical point triggers magnetic instabilities.

Recurrent energy releases occur with each helicity injection.

Surface temperature increases with injected twist, especially in tropical regions.

Abstract

In this letter, we present the first 3D force-free general relativity simulations of the magnetosphere dynamics related to the magnetar outburst/flare phenomenology. Starting from an initial dipole configuration, we adiabatically increase the helicity by twisting the footprints of a spot on the stellar surface and follow the succession of quasi-equilibrium states until a critical twist is reached. Twisting beyond that point triggers instabilities that results in the rapid expansion of magnetic field lines, followed by reconnection, as observed in previous axi-symmetric simulations. If the injection of magnetic helicity goes on, the process is recurrent, periodically releasing a similar amount of energy, of the order of a few % of the total magnetic energy. From our current distribution, we estimate the local temperature assuming that dissipation occurs mainly in the highly resistive outermost layer of the neutron star. We find that the temperature smoothly increases with injected twist, being larger for spots located in the tropical regions than in polar regions, and rather independent of their sizes. After the injection of helicity ceases, the magnetosphere relaxes to a new stable state, in which the persistent currents maintain the footprints area slightly hotter than before the onset of the instability.

Figures

Equations

$$\Omega (\theta ,t)=\frac{\omega (t)}{1+\exp \left[\kappa (\theta -{\theta }_s)\right]}$$

$$T_{\text{bb}}\approx 0.18\mathrm{keV}{\left[\frac{\Delta r}{1\mathrm{m}}\right]}^{\frac{1}{4}}{\left[\frac{10^{17}{\mathrm{s}}^{-1}}{{\sigma }_e}\right]}^{\frac{1}{4}}{\left[\frac{J}{10^{18}G/s}\right]}^{\frac{1}{2}}.$$

Full text

††pagerange: Triggering magnetar outbursts in 3D force-free simulations–LABEL:lastpage

1 Introduction

Magnetars Olausen & Kaspi (2014); Kaspi & Beloborodov (2017); Mereghetti et al. (2015), arguably the most magnetized stable objects in the universe, are relatively slowly rotating neutron stars with spin periods of $P\sim 1-10$ s and estimated magnetic fields of $B\sim 10^{14}-10^{15}$ G. The number of the known magnetars has been steadily increasing Olausen & Kaspi (2014)111See the McGill catalog: http://www.physics.mcgill.ca/$\sim$ pulsar/magnetar/main.html since the mid 90’s where the main aspects of the magnetar theory were proposed Duncan & Thompson (1992); Thompson & Duncan (1993, 1995, 1996). Its most relevant distinguishing feature, compared to common radio-pulsars, is that the observed electromagnetic luminosity is well above the loss of rotational energy, indicating that the magnetic field is the energy source, instead of rotation. Magnetars are also characterized by energetic transient phenomena, detected in X-rays and $\gamma-$rays, ranging from short bursts lasting fractions of second, to long outbursts during which their luminosity suddenly increases, slowly returning to a low-state on long timescales (months or even years) Rea & Esposito (2011); Coti Zelati et al. (2018), and sporadic giant flares (observed in three cases). All these events are supposed to be linked to reconnection events of the magnetospheric magnetic field, possibly triggered by the star’s interior evolution. Thus, understanding the cause and the dynamics of magnetospheric activity becomes a key issue in magnetar theory.

In the rarefied magnetospheric environment, the electromagnetic forces dominate over particle inertia and allows to assume the force-free approximation, which has been widely used to study global properties of pulsars Gruzinov (1999); Contopoulos et al. (1999); McKinney (2006); Timokhin (2006); Spitkovsky (2006). This approximation has also been adopted to model magnetar outbursts, studying the dynamical response of the magnetosphere under the injection of magnetic helicity from a localized spot on the surface. If the injection timescale (determined by the internal field evolution) is much longer than the rapid response of the magnetosphere (milliseconds), the evolution of the twisted magnetic configurations can be modeled as an adiabatic sequence of quasi-stationary states, at least as long as the configurations are stable. Once a magnetic line bundle is twisted above a critical value, the system becomes unstable and eventually undergoes a large scale magnetic reconnection and plasmoid emission. The dynamics of these events has been extensively studied by analytical and numerical models considering axi-symmetric configurations Beloborodov & Thompson (2007); Beloborodov (2009); Parfrey et al. (2013); Pili et al. (2015); Chen & Beloborodov (2016).

In this work we aim at extending previous models in two different directions. First, we include General Relativity (GR) effects by considering the Schwarzschild metric, as a good approximation to describe the spacetime surrounding a relatively slow rotating star. Second, and more important, by performing for the first time 3D simulations without any symmetry restriction.

2 Numerical Setup

We use the version of the GR force-free electrodynamics formalism derived in Carrasco & Reula (2016), which is well-posed and involves the full force-free current density. More concretely, we implement the evolution system given by Eqs. (8)-(10) in Carrasco et al. (2018). The numerical scheme to solve these equations is based on the multi-block approach Lehner et al. (2005); Carpenter et al. (1994); Carpenter et al. (1999); Nordström & Carpenter (2001); Calabrese et al. (2004). We employ difference operators which are eight-order accurate on the interior and fourth-order at the boundaries, with adapted Kreiss-Oliger dissipation operators Diener et al. (2007). A classical fourth order Runge-Kutta algorithm is used for time integration. Our numerical grid extends from an interior sphere at radius $r=R$, that represents the star radius, to an exterior spherical surface at $r\sim 50R$. We adopt a resolution $N_{\theta}\times N_{\phi}\times N_{r}$ with $N_{\phi}=321$, $N_{\theta}=161$, while $N_{r}$ is adjusted to satisfy $\Delta r\lesssim r\Delta\theta$ everywhere in the domain.

The stellar surface is assumed to behave as an idealized perfect conductor. Thus, at the surface, the normal component of the magnetic field is fixed and the electric field is prescribed from the ideal MHD condition consistent with the shear perturbations at the stellar crust. At the outer boundary, on the other hand, we set non-reflective conditions to allow all perturbations to propagate away. The numerical implementation of such boundary conditions and how to handle current sheets have been detailed in Carrasco et al. (2018); Carrasco & Reula (2017).

We consider an initial magnetic dipolar field configuration in a Schwarzschild external region $r>R$, of mass $M$ Wasserman & Shapiro (1983). We adopt geometrized units in which $c=G=1$, and Lorentz-Heaviside units for the electromagnetic field. A typical neutron star of compactness $M/R=0.2$ is considered. In our coordinates, the magnetic moment forms an angle $\chi$ with the $z$-axis, which represents the center of the circular surface region of the sheared footprints. This misalignment angle $\chi$ is therefore the angle between the NS’s magnetic dipole and the center of the twisted region. The profile describing the perturbation is given by:

[TABLE]

where $\theta_{s}$ is the angular extension of the spot, and we set $\kappa=30$. The value of $\omega(t)$ smoothly increases from zero, at $t=0$, to its final value $\omega_{o}=2\times 10^{-3}\left[1/M\right]$, at $t=400M$.

We quantify the accumulated injected twist as, $\psi_{\rm inj}(t)=\int_{0}^{t}\omega(t)dt$. The magnetic energy normalized to the initial potential state, $W/W_{0}$, serves as a measure of the additional magnetic energy stored in the magnetosphere. We also monitor the currents during the dynamics, which will allow later to estimate the temperature of the surface.

3 Results

We first focus on a representative case, with a misalignment angle $\chi=0.15\pi$ and angular cap size $\theta_{s}=0.15\pi$, for which the injection is smoothly turned off after some time. The shearing of the footprints launches Alvén waves along the magnetic filed lines, which are reflected back and forth within the closed region of the magnetosphere. The complicated pattern of wave-fronts eventually reaches a quasi-equilibrium state, provided that the shearing timescale is much longer than the wave-crossing times of the twisted magnetic lines. As a consequence, the following description dynamics are valid for a large range of $\omega_{o}$.

The qualitative behavior of the magnetospheric response can be described by several stages, illustrated in Fig. 1, where some representative magnetic field lines and the radial component of the current $J^{r}$ at the stellar surface are shown. After the first short initial transient described above, the injected helicity gradually modifies the magnetic configuration in an adiabatic succession of non-potential equilibrium states (top panel). The evolution of the stored magnetic energy, $W(\psi_{\rm inj})$, and the associated currents, is independent of the injection rate $\omega_{o}$ (provided it is not too fast). Beyond a critical twist, $\psi_{\rm crt}\sim 3$–$4$, the system becomes unstable and departs from the quasi-equilibrium configurations. In this stage (middle panel), some magnetic field lines quickly expand (approaching an open state) and produce a tangential discontinuity in the magnetic field near the equatorial plane 222 Such critical behavior, in which the field experiences an explosive expansion due to the loss of ideal magnetic equilibrium, has been studied in more detail (in axial symmetry) at e.g. Mikic & Linker (1994); Uzdensky (2002); Parfrey et al. (2013).. Then, the plasma resistivity (physical in nature, numerical in our case) leads to magnetic reconnection in these current sheets and the subsequent expulsion of plasmoids. Finally, after the reconnection, the system relaxes into a quasi-equilibrium configuration which differs from those reached during the first quasi-adiabatic stage (bottom panel). Part of the twist and the associated currents remain stable in a rearranged configuration, with a larger magnetic energy than the initial vacuum dipole.

These sudden reconnection events, release and dissipate significant fractions of the energy stored in the magnetosphere. Figure 2 illustrates the magnetic field dynamics in a sequence of one such events: the quasi-equilibrium configurations followed by the unstable expansion of the magnetic field, the formation of a current sheet in the equatorial plane, and the subsequent ejection of plasmoids with different sizes. The main feature introduced by the misalignment being that reconnection and emissions are no longer axially-symmetric as for the case $\chi=0$. There is now a preferred direction, defined by the intersection of the $x-z$ plane with the dipole’s equator, along which these plasmoids propagate at nearly the speed of light. We turn to a detailed analysis of the energetics of these events, potentially available to power outbursts and flares. The accumulated energy prior to reconnection, due to the presence of twisted magnetic fields, constitutes about $2.7\%$ of the initial energy of the dipole, which can be estimated as $W_{0}\sim B_{15}^{2}R_{10}^{3}\times 10^{47}$ erg, being $B_{15}\equiv\frac{B}{10^{15}{\rm G}}$ and $R_{10}\equiv\frac{R}{10\rm{km}}$ the normalized values (assumed to be 1 hereafter) of the dipolar magnetic field at the polar surface and of the star’s radius. From this available energy reservoir, approximately $20\%$ is expelled from the system, while nearly $30\%$ is dissipated 333We compute this quantity by tracking the energy associated with the electric field removed from the system at the current sheet, in order to maintain the condition $B^{2}-E^{2}>0$, as in other force-free studies. The numerical dissipation is instead negligible. during the process.

The plasmoids are presumably heated at their formation in the current sheet and propagate both outwards and inwards. Those “plasma fireballs” trapped by the NS magnetic field were suggested to explain the few minutes tail after the main flash in giant flares, which is modulated by the spin period of the NS. On the other hand, part of the energy dissipated at the current sheet might be radiated right away at the explosive reconnection, producing the short $0.1-0.2\,s$ flash in luminosity (e.g. Gourgouliatos & Esposito, 2018). The remaining $50\%$ of the energy stored during the injection remains in the magnetosphere, such that the new equilibrium attained after the event is not a potential solution, but a different twisted magnetosphere configuration.

Following Akgün et al. (2018), we assume that, to close the circuit, part of the magnetospheric currents must be dissipated as they go through the very thin and resistive outer layer of the star ($\Delta r\sim$ few m), with electrical conductivity444Calculated with the public routines available at: http://www.ioffe.ru/astro/conduct/ $\sigma_{e}$, (see Potekhin et al. (2015) and references therein). We can then obtain the map of the estimated surface temperature by equating the local Joule dissipation rate to the blackbody emission from the heated surface: $\sigma T_{\text{bb}}^{4}=\Delta rJ^{2}/\sigma_{e}$, where $\sigma$ is the Stefan-Boltzmann constant, which gives:

[TABLE]

On Fig. 3 we show the effective temperature, defined by the following average over the magnetic spots surface $T_{\rm eff}^{4}=<T_{bb}^{4}>$, where $T_{bb}$ is the local surface temperature estimated as described above, assuming the same values as in Akgün et al. (2018): $\Delta r=1$ m, and $\sigma_{e}=7\times 10^{16}$ s*-1*. During the quasi-adiabatic stage, the effective temperature grows with the injected twist (as $T_{\rm eff}\sim\sqrt{\psi_{\rm inj}}$ in the perturbative regime, $\psi_{\rm inj}\lesssim 1$), independently of the particular twist rate employed. As it is shown in the figure, if the injection is stopped before reaching the critical value, the system remains stable with a constant distribution of the currents. Otherwise, the instability is triggered and the temperature quickly increases during the expansion phase. During the reconnection event, and soon after, the fireball internal energy and other physical mechanisms (like accelerated particles) will likely provide additional contributions to the X-ray spectra, which can be temporarily dominant. However, when the system reaches a stable state (after a few crossing times in our case), the final temperature settles down into an approximately constant value which is always slightly larger than those attained prior to the onset of the instability. Even though this behavior is not intuitive, our interpretation is that the expansion phase rapidly amplify the magnetospheric currents. While part of these currents are lost during the reconnection, a significant fraction still remain in the final relaxed system, resulting in the observed increase of temperature.

The most extreme case, when a constant twist is injected continuously even after reaching the critical value, shows a similar behavior, with an asymptotic temperature that grows very slowly in the timescale of our simulations. The plot also illustrates the role of GR effects by comparing the solutions at two different stellar compactness. It can be seen that when curvature is neglected ($M=0$), the effective temperature reduces by apparently $\sim 20\%$. Such increase is due to the combination of the Schwarzschild metric factors entering in the definition of the fields and their derivatives. However, we note that this is the local temperature in the star. For an observer at infinity the temperature must be red-shifted in the relativistic case, and therefore the observed temperatures would be similar in both cases.

We also explored different misalignment angles for a fixed spot size of $\theta_{s}=0.15\pi$, and different spot sizes for a given angle $\chi=0.15\pi$. The results are summarized on Fig. 4, where the total magnetic energy in the magnetosphere is displayed as a function of $\psi_{\rm inj}$. In this sequence of models, we have applied a continuous shear, without turning it off at any moment. Thus, driven by this continuous twist injection, reconnection events occur repeatedly. Since each event dissipates or expels only part of the additional magnetic energy, there is an slight increase of the stored energy with time.

These results provide an estimate of how much energy can be stored and later released depending on the two main parameters of the problem, namely: $\chi$ and $\theta_{s}$. The available energy stored prior to reconnection increases with both the misalignment angle and with the size of the shearing spots, and is of the order on only a few percent of $W_{0}$. The total twist at which the first reconnection takes place, $\psi_{\rm rec}$, varies between 4 and 11 radians. However, some cautionary comments are necessary. First, $\psi_{\rm rec}$, depends on numerical resolution and –fundamentally– on the twisting rate $\omega_{o}$. As observed in Parfrey et al. (2013), large values of $\omega_{o}$ produces a dynamical stabilizing effect which tends to delay the onset of reconnection beyond the critical twist value $\psi_{\rm crt}$. Thus, in general one has that $\psi_{\rm rec}\geq\psi_{\rm crt}$. We find, nevertheless, that such dynamical stabilization is very prominent in the aligned cases, where the symmetry favors the effect, but it is not that important in the misaligned cases studied here.

The maximum temperatures attained near the critical twist value, $T_{\rm crit}$, also increase with the misalignment $\chi$, as shown in Fig. 5: the position and geometry of the twisted lines is crucial in determining the heat released at the surface. On the other hand, we find that changing the magnetic spot size only affects the global energetics (larger volumes involved), but not the estimated effective temperature.

4 Discussion

In the framework of the standard magnetar model, we have studied the magnetospheric response to the smooth injection of currents from the star interior, supposedly originated by the slower internal magnetic field evolution. Our results from GR3D simulations show that, starting with a potential model, the magnetosphere gradually stores helicity and energy until a critical twist value is reached, after which a reorganization occurs through magnetic reconnection and a new, quasi-steady, magnetospheric configuration is reached. The misaligned cases exhibit important differences respect to the axi-symmetric scenarios, in terms of the dependence of the additional magnetic energy stored for a given twist and the temperature of the heated surface on the position of the twisting footprints.

Connecting our results with observations pre- and post-outburst (for the about ten cases where pre-outburst data are available, Coti Zelati et al. (2018)) we can draw the following conclusions. First, we find that stable, non-potential configurations can persist after an event, and the estimated temperatures of the magnetic spots are in line with the high observed values of $0.5$–$1.0$ keV, after the outburst occurs. Second, outburst observations show a sudden flux enhancement by up to a few orders of magnitude, accompanied by an increase of the inferred blackbody temperature by a factor of a few (even considering the many caveats related with spectral fitting uncertainties). We have no precise information about the outburst rise timescale, because no X-ray dedicated observation has ever been lucky enough to catch a magnetar while undergoing an outburst. However, we can place an upper limit at about 1.5 days, which is the closest X-ray observation in quiescent state immediately preceding an outburst Coti Zelati et al. (2018). Since the injection of twist is coupled to a continuous increase of temperature, the fact that the temperature suddenly increases on fast timescales ($<$ one day) indicates that the last stages of the helicity injection must also be very rapid. This suggests that the outburst trigger is not purely driven by the slow, smooth crustal magnetic evolution (timescales of kyrs at least), because that alone should cause a smaller pre- and post-outburst temperature variation (no more than $\sim 20\%$ in our simulations). Instead, the final pre-outburst stage would be compatible with shorter timescale phenomena of sudden magnetic helicity injection (likely triggered, in turn, by the long-term crustal evolution), such as thermo-plastic waves Beloborodov & Levin (2014) or crustal fractures Thompson & Duncan (1993); Perna & Pons (2011). Moreover, for the handful of sources with good-quality data spanning many years, the so-called quiescent state shows no trend of significant spectral/flux evolution over years, indicating that stable twisted magnetospheres supported by slow helicity injections maintaining a small twist below the critical value may be very long lived.

A work in preparation is going to further describe the numerical details and explore a wider range of cases, in terms of spot size, location and injection rate. It will be important in the future to also consider different initial base configurations and injection distribution, including for instance higher multipoles, or allowing the possible interaction between different twisted bundles. A connection with the internal evolution of isolated neutron stars will allow to better assess the shape and distribution of very slow twist injection, and the resulting magnetosphere. Last, observationally the few magnetars showing repeated outbursts can helpful in partially constraining the location and size of the twisted footprints.

Acknowledgments

We acknowledge support from the Spanish Ministry of Economy, Industry and Competitiveness grants AYA2016-80289-P and AYA2017-82089-ERC (AEI/FEDER, UE), and AYA2015-66899-C2-2-P. CP also acknowledges support from the Spanish Ministry of Education and Science through a Ramon y Cajal grant. This work used computational resources from Pirayu Cluster (supported by the Agencia Santafesina de Ciencia, Tecnología e Innovación, Gobierno de la Provincia de Santa Fe, Proyecto AC-00010-18), Centro de Computación de Alto Desempeño (CCAD) and Centro de Cómputos de Alto Rendimiento (CeCAR). All part of the Sistema Nacional de Computación de Alto Desempeño, MinCyT-Argentina. The work has been done within the PHAROS COST action CA16214.