Particle-in-cell simulations of the twisted magnetospheres of magnetars

The magnetospheres of magnetars are believed to be filled with electron-positron plasma generated by electric discharge. We present a first direct numerical experiment showing how the plasma is created in an axisymmetric closed magnetosphere. The $e^\pm$ discharge occurs in response to twisting of the magnetic field lines by a shear deformation of the magnetar surface, which launches electric currents into the magnetosphere. The simulation shows the formation of an electric"gap"with unscreened electric field ($\mathbf{E}\cdot \mathbf{B}\neq 0$) that continually accelerates particles along the magnetic field lines and sustains pair creation. The accelerating voltage is self-regulated to the threshold of the $e^\pm$ discharge. It controls the rate of energy release and the lifetime of the magnetic twist. The simulation follows the global evolution of the twisted magnetosphere over a long time and demonstrates its gradual resistive untwisting. A vacuum cavity forms near the star and expands, gradually erasing magnetospheric electric currents $j$. The active j-bundle shrinks with time and its footprints form shrinking hot spots on the magnetar surface bombarded by the created particles.


INTRODUCTION
Magnetars are neutron stars with ultrastrong magnetic fields (B 10 14 G) that display strong activity fed by dissipation of magnetic energy (see e.g. Mereghetti 2008;Turolla et al. 2015 for reviews). They produce strong outbursts and flares as well as bright persistent emission with a prominent hard X-ray component extending above 100 keV. These activities are associated with strong deformations of the external magnetosphere of the neutron star, resembling the activity of the solar corona (e.g. Thompson & Duncan 1995). The magnetosphere is anchored in the solid crust of the star and its deformation is caused by crustal shear motions driven by ultrastrong internal magnetic stresses.
The speed of the surface motions is poorly known. Recent work suggests that the crust yields to internal stresses through an instability launching a thermoplastic wave (Beloborodov & Levin 2014) or a Hall-mediated avalanche  In both cases the motion is plastic and should occur on a timescale much longer than the Alfvén crossing timescale (10 − 100 ms). It is expected to be fast enough to efficiently twist the external magnetosphere.
The surface shear motion launches Alfvén waves along the magnetic field lines and generates magnetospheric twist ∇ × B = 0 (Thompson et al. 2002;Parfrey et al. 2013, hereafter PBH13). Plasma is required to supply the current j = (c/4π)∇ × B. Beloborodov & Thompson (2007, hereafter BT07) found that plasma must be mainly supplied through e ± discharge in the magnetosphere rather than through extraction of charges from the star. They performed simplified one-dimensional (1D) simulations of the discharge. In the simulations, the magnetosphere was replaced by a fixed, uniform field B(x) connecting anode and cathode -metallic plates at x A and x C . The fixed ∇ × B in this setup turns out equivalent to imposing an electric current through the plates into the computational box. When pair creation was not allowed, the system quickly relaxed to a global "double layer" configuration, with surface charges of the opposite sign induced on the plates. The electric field between them gave a huge voltage Φ e accelerating particles to ultra-high energies. When pair creation process was included in the simulation, the voltage dropped to a much lower value, just sufficient to sustain pair creation, and the current was supported through continual e ± discharge. BT07 concluded that pair creation must be responsible for screening electric fields and regulating the magnetospheric activity of magnetars. The simplified 1D model cannot, however, give a compete picture of the magnetospheric activity, for a few reasons. It does not show how ∇ × B is imparted in the first place, as the 1D model does not support Alfvén waves. The exclusion of this important degree of freedom may also put in question the double layer formation in the absence of pair creation, the necessity of the onset of pair creation, and the self-regulation of the discharge voltage seen in the 1D model. Note also that the electric field in the 1D (slab) geometry does not decrease with distance from the charge, and hence one cannot see a realistic distribution of the accelerating electric field along the magnetospheric field lines. Finally, the 1D model offers no way to follow the gradual resistive "untwisting" of the magnetosphere -its global evolution as a result of ohmic dissipation of the twist energy. The expected evolution must occur on the resistive timescale of months to years (regulated by voltage Φ e ) and can be tested against observations. An axisymmetric electrodynamic model of a resistively untwisting magnetosphere was developed by Beloborodov (2009, hereafter B09). This model assumed that a given fixed voltage Φ e is sustained on currentcarrying field lines, without calculating the discharge that regulates Φ e . A surprising result was the formation of two distinct regions in the untwisting magneto-sphere, with a sharp boundary between them, -a "cavity" (j = 0) and a "j-bundle." In essence, the untwisting process was found to be the growth of the cavity, erasing the currents in the j-bundle. A curious immediate implication was the prediction of shrinking hot spots on magnetars -the footprints of the shrinking j-bundle, where the stellar surface is heated by bombardment of accelerated particles. Shrinking hot spots have been observed in seven objects by now (see data compilation in . All of these objects belong to the class of "transient magnetars" that show a sudden outburst and then gradually decay back to the quiescent state of low luminosity. A key parameter governing the j-bundle evolution is its poorly known voltage Φ e , which depends on how the e ± discharge is self-organized and may be different on different magnetospheric field lines.
The goal of the present paper is to overcome the limitations of the 1D discharge model and perform a first selfconsistent calculation of the e ± discharge in an axisymmetric twisted magnetosphere. The process can be simulated from first-principles using a full kinetic description of the magnetospheric plasma as a large number of charged particles moving in the self-consistent collective electromagnetic field. Such a direct numerical experiment will show how the twist and the electric current are created in the magnetosphere in response to crustal shear, and will follow the ensuing dissipative evolution of the twist.
The self-organization of the e ± discharge should determine where the particles are created and accelerated. Should this occur near the footpoints of the magnetospheric field line or near its apex? Will the acceleration region be steady or move around? Answers to these questions may have important implications for nonthermal emission from magnetars. The voltage drop along the twisted field lines will control the dissipated power which feeds the observed emission. We expect to see how particles are accelerated in the current-carrying magnetic loop and rain down on the stellar surface to create hotspots. Finally, the established discharge voltage will determine the life-time of the magnetic twist and the pattern of its evolution.
A suitable technique for such direct numerical experiments is the particle-in-cell (PIC) method, with pair creation implemented. This method has been successfully applied to the old problem of rotation-powered pulsars (Chen & Beloborodov 2014, hereafter CB14;Philippov et al. 2015;Belyaev 2015;Cerutti et al. 2016) The magnetar problem is different in important ways and in some ways easier to study using a global PIC simulation, as will be described below.
The paper is organized as follows. In Section 2 we describe the theory of twisted magnetospheres in axisymmetric geometry, revisit the double-layer configuration (in the absence of pair creation), describe the mechanism of pair creation and basic electrodynamics of untwisting. This will be useful for understanding the simulation results and also introduces notation used in the paper. Section 3 presents the setup of our numerical experiments. Section 4 describes the results and their implications. Finally in Section 5 we summarize our conclusions and provide an outlook for future studies.
2. SUSTAINING CURRENTS IN THE TWISTED MAGNETOSPHERE Let us consider a dipole magnetic field around the star, and assume that its footpoints on the star are sheared in the azimuthal direction about the magnetic axis. In this case, the implanted twist is axisymmetric and its amplitude ψ is simply given by the azimuthal angle between the two footpoints of the magnetic field line. It is convenient to use spherical coordinates r, θ, φ with the polar axis being the axis of symmetry. The magnetospheric twist implies a toroidal component of the magnetic field B φ = 0, and the twist amplitude on a given magnetic field line is related to B φ by where the integral is taken along the field line, and p, q are the two footpoints where the field line is anchored to the surface. As long as the implanted twist ψ is smaller than unity, the poloidal magnetic field remains close to dipolar, and the deformation can be thought of as simply adding a toroidal component B φ without changing the poloidal dipole component (B09). This induces ∇ × B in the dipolar configuration that was originally curl-free. It must be sustained by an electric current in the magnetosphere, j. The magnetic energy strongly dominates over the plasma energy, and hence the currents must be nearly force-free, j × B = 0, i.e. flowing along the magnetic field lines.
The origin of the plasma that could carry the current is a non-trivial issue. The star can have a gaseous atmosphere, however for the typical surface temperature kT < 1 keV the atmosphere scale-height is tiny (centimeters), because of the strong gravity of the neutron star. The atmosphere does not provide enough plasma to conduct currents at large altitudes r ∼ R , where R = 10 − 13 km is the neutron star radius.
Spinning of the neutron star and its magnetosphere with velocity v rot = Ω × r implies a "co-rotation" electric field E = −v rot × B/c and requires charge density ρ GJ = ∇·E/4π = −Ω·B/2πc (Goldreich & Julian 1969). Magnetars are slow rotators, Ω ∼ 1 Hz, and their ρ GJ is small. The currents demanded by the twisted magnetosphere are typically much stronger than cρ GJ .
The magnetosphere must make a special effort to avoid charge starvation and create sufficiently dense plasma to conduct the current j demanded by the twist. It achieves this by inducing an electric field E (parallel to the magnetic field lines) that can accelerate particles and trigger pair creation. This implies a finite voltage in the magnetospheric electric circuit and a finite rate of ohmic dissipation.

Voltage without pair discharge
In the absence of pair creation, the star is the only available source of magnetospheric plasma. The lack of charges leads to induction of an electric field with a component parallel to the magnetic field, which can pull out charges from the star and accelerate them. Then the electric circuit is expected to relax to a static configuration similar to the relativistic double layer derived by Carlqvist (1982) and observed in the 1D plasma simulations of BT07. It sustains the opposite surface charges at the two footpoints of the magnetic loop where the lifted particles still move slowly, v c, and create a large charge density ρ ∼ j/v.
The high charge density near the footpoints generates E according to the Gauss law, and E accelerates the flow on the plasma timescale ω −1 p = (m e /4πeρ) 1/2 . The flow density ρ is reduced to its minimum where its velocity approaches c. As a result, the characteristic thickness of the surface charge layer is the plasma skin depth λ p = c/ω p evaluated for the plasma density ρ ∼ j/c.
The surface charge Σ ∼ (j/c)λ p generates the selfconsistent electric field that lifts and accelerates particles from the footpoint, where ω p is the plasma frequency defined by In other words, the surface charge near the anode and cathode is organized so that particles extracted from the star are accelerated to v ∼ c over a length comparable to the plasma skin depth. For simplicity, consider a symmetric double layer where the positive and negative charges have the same mass. In the 1D model, the electric field is almost constant between the two surface charges of the double layer, giving a voltage drop, where L is the size of the layer (the distance between the footpoints). Using j ∼ ψB/L, one finds for the typical parameters of a magnetar, which implies a huge voltage Φ e . The estimate in Equation (4) is not valid, however, for a realistic magnetosphere, which is not one-dimensional. The current flows along the curved magnetic field lines and their dipolar geometry significantly changes the distribution of the net voltage sustained between the two footpoints.
The corrected voltage may be estimated as follows. Since λ p is small compared with the thickness of the jbundle, the surface charge remains thin and its structure is not changed from the 1D model. The self-consistent electric field extracting charges from the footpoint is still described by Equation (2). However, with increasing altitude the electric field must be reduced on a scale comparable to the horizontal size of surface charge W (thickness of the j-bundle). The resulting potential drop saturates at Φ e ∼ E W , which gives It is smaller than the 1D estimate by the factor of W/L. For instance a j-bundle of thickness W ∼ 0.1R at the stellar surface and length L ∼ 10R would sustain a voltage ∼ 10 −2 smaller than predicted by the 1D model. This is still a huge voltage and particles that tap the full potential drop will be able to induce pair discharge, making the double layer model inconsistent.
One should also note that E = E · B/B, and hence the voltage, have a pure inductive origin. One should think of E as c −1 ∂A /∂t, the result of the slow decay of the ultrastrong twisted magnetic field (BT07). eΦ e measures the energy gain of charge e completing the electric circuit, and this released energy is extracted from the magnetic twist energy. A potential electric field would be unable to support any significant voltage between the footpoints, as they are connected by an excellent conductor -the crust.
The induction electric field E still satisfies the Gauss law ∇ · E = 4πρ; as long as the untwisting process occurs much slower than the light crossing of the system, one can think of the dissipation as a quasi-steady process. The inductive double layer is similar to a normal electrostatic double layer except that the integral of E along the full closed circuit (including the part closing through the crust, where E = 0) does not vanish and instead equals Φ e . There is no external emf applied to the circuit below the stellar surface; the only emf sustaining the current is the induction emf due to the twist decay in the magnetosphere itself.

Voltage with pair discharge
The mechanism of secondary e ± creation by relativistic particles in the magnetar magnetosphere involves an intermediate step of gamma-ray production. It occurs through resonant Compton scattering of photons flowing from the star by particles accelerated in the magnetosphere. A target photon with energy E t ∼ 1 keV can be resonantly scattered by an electron with Lorentz factor γ if the photon energy measured in the electron rest frame matches ω B , where ω B = eB/m e c. 1 The resonance condition reads where θ X is the angle of the target X-ray with respect to local magnetic field line (the electron moves along the field line). The energy ω B equals m e c 2 for the characteristic magnetic field B Q = m 2 e c 3 /e ≈ 4.4 × 10 13 G and scales linearly with B.
Magnetars supply plenty of keV photons, and the electron Lorentz factor required for resonant scattering at B ∼ B Q is moderate, γ ∼ 10 3 . It is far below the electron Lorentz factors that would be reached in the double layer discussed in the previous section.
After the scattering, the photon energy is boosted by a factor comparable to γ 2 , putting the originally keV photon into the GeV range, E γ ∼ 1 GeV. Such energetic gamma-rays can easily convert to e ± pairs in the strong magnetic field, as soon as the gamma-ray pitch angle with respect to the magnetic field, θ γ , is large enough to satisfy the threshold condition, In the region near the star where B > 10 13 G the conversion occurs practically immediately following resonance scattering (Beloborodov 2013). The efficiency of pair creation implies a quick development of electric discharge until the number of created particles becomes sufficient to screen the accelerating electric field. The process develops in a runaway (exponential) manner and hence the accelerating voltage is unlikely to grow beyond a characteristic value that makes particles capable of resonant scattering. This condition defines a "threshold" for discharge, which corresponds to a characteristic electron Lorentz factor γ thr .

Characteristic timescales and energy scales
The shortest timescale of interest is the plasma scale ω −1 p . It describes the growth rate of the local accelerating electric field in response to charge starvation (BT07). It also determines the thickness of the surface charge c/ω p in the double-layer configuration.
The characteristic dynamic timescale of the electric circuit is the light crossing time or the Alfvén crossing time of the system, where L is the length of the magnetospheric field line. The group speed of Alfvén waves is always directed along the magnetic field lines and its value is close to c in the magnetically dominated corona. The longest timescale in the problem is the lifetime of the magnetic twist. The finite voltage sustaining the magnetospheric current implies a finite ohmic dissipation rate, so the magnetic twist energy E twist must dissipate with time, where I is electric current flowing through the magnetosphere. The voltage Φ e controls the timescale of this evolution, Using the characteristic I < ∼ ψ(c/4π)BR and γ thr ∼ 10 3 one can estimate that t ohm is comparable to one year. This theoretical timescale for untwisting is comparable to the observed decay timescale in transient magnetars following an outburst of activity. Because of the vast separation of timescales, t ohm t A , the ohmic dissipation of the magnetospheric twist can be viewed as a quasi-steady process slowly draining the twist energy. Unsteadiness of the discharge may lead to strong variability in the electric circuit, however it occurs on very short timescales, which would be hard to resolve observationally.
The characteristic scales for energy (or electron Lorentz factor γ) also have an important hierarchy. The highest energy corresponds to γ DL , which would only be achieved in the absence of pair creation. It is given by Equation (6) and can exceed 10 6 . The next characteristic γ is determined by the threshold for e ± discharge γ thr , which is comparable to 10 3 . Both γ DL and γ thr are much greater than unity.

Mechanism of untwisting
An integral form of the Faraday's induction law ∂B/∂t = −c∇ × E leads to a simple equation describing resistive evolution of the axisymmetric twist (Beloborodov 2011),ψ = 2πc ∂Φ e ∂f .
Here f (r, θ) is the poloidal magnetic flux function (constant along a magnetic flux surface), which serves to label the magnetic field lines. For any given point (r, θ), f is defined as the magnetic flux through the circle about the axis of symmetry passing through the point; f = 0 on the axis of symmetry. In particular, for a dipole poloidal field with a dipole moment µ the flux function is given by Note that sin 2 θ/r = const along a dipole field line. It is convenient to use the dimensionless flux function where θ is the polar angle of the magnetic field line footprint on the stellar surface. Equation (13) shows that the twist must decrease where ∂Φ e /∂f < 0 and increase where ∂Φ e /∂f > 0. The fact that Φ e (f max ) = 0 (the field line f max is confined to the star, which we approximate as an ideal conductor) implies ∂Φ e /∂f < 0 at some f < f max . This region with large f , comparable to f max , corresponds to the inner magnetosphere near the equator, with short field lines. B09 showed that this fact leads to immediate formation of a "cavity" with j = 0 in the equatorial region near the star, and the cavity expands on the timescale t ohm , erasing the magnetospheric currents. The currents are "sucked" into the star, so that they close inside the conductor.
From the untwisting equation it is evident that the profile of Φ e (f ) plays the key role for the twist evolution. Voltage regulated by pair discharge is expected to satisfy the condition eΦ e ∼ γ thr m e c 2 . Its variation with f over a region ∆f = f max ∆u gives the characteristic twist evolution timescale, The dimensionless quantities ∆u and ψ are comparable to unity, and the characteristic timescale is set by the ratio µ/Φ thr . Note however that t ohm can strongly differ for different magnetic field lines. In particular, if there is a region with a flat dependence of Φ e (f ), ∂Φ e /∂f = 0, then the local t ohm = ∞ and the twist angle ψ is "frozen", waiting for the cavity expansion to reach the region (B09).
Another interesting implication of Equation (13) is that on some field lines the twist may grow as the magnetosphere untwists. In particular, a decrease of Φ e toward the magnetic axis, ∂Φ e /∂f > 0, leads toψ > 0. This effect will be observed in the simulations below. Together with the cavity expansion, this means that the twist relocates toward the axis with a decreasing energy E twist but possibly with increasing amplitude ψ in some regions before being completely dissipated.
3. SETUP OF THE SIMULATION 3.1. Implanting the twist Our simulation starts with a pure dipole magnetosphere, with a magnetic moment µ and no magnetic twist, B φ = 0. The twist is gradually implanted by shearing the stellar surface with a latitude-dependent angular velocity ω(θ) µ. The profile of ω(θ) determines the profile of the implanted twist; we choose a profile similar to previous magnetohydrodynamic (MHD) and force-free electrodynamic (FFE) simulations of twisted magnetospheres (Mikic & Linker 1994;PBH13), where Θ = (θ − π/2)/∆θ m and ∆θ m = π/4 is a measure of the width of the sheared region. This profile gives a smooth twist that is centered at θ = π/4 and decreases to zero at the equator. The prefactor ω 0 (t) describes the rate of implanting the twist. It is smoothly increased from zero at t = 0 to a chosen maximum value, kept at this value for some time, and then smoothly switched off back to zero. As long as the duration t shear of the surface shear ω = 0 is shorter than the resistive timescale of the magnetosphere, t shear t ohm , ohmic dissipation may be neglected during time t shear . Then the implanted twist profile is given by We choose t shear = 40R /c. Then the shearing stage is sufficiently short compared with the total duration of our simulation t sim = 350R /c but longer than or comparable to the Alfvén crossing time t A of the sheared region, so that twist implanting is a relatively gentle process. The maximum shear angle (near θ = π/4) is ψ max ≈ 1.6 radian in the simulations presented below. After the twist implantation is finished, ω is kept at zero and the boundary condition at the stellar surface becomes simply a perfect static conductor. Magnetars are slow rotators, and their light cylinders R LC > ∼ 10 4 R are well beyond the twisted, dissipative region. The slow spinning of the star is neglected in the present paper, which corresponds to R LC = ∞.
The implanted twist ψ ∼ 1 is moderate and expected to result in moderate inflation of the poloidal magnetic field lines. The main effect of surface shearing is creating a strong B φ in the magnetosphere. Analytical arguments (e.g. Uzdensky 2002) and FFE simulations (PBH13) show that a stronger ψ > ∼ 3 will result in a global instability of the magnetosphere, which we do not intend to study in this paper and defer to future work.

Surface atmospheric layer
We start the simulation with a complete vacuum around the star and create a dense neutral atmospheric layer at the stellar surface by injecting warm electron-ion plasma at R . The atmosphere scale-height h is determined by the particle injection temperature and gravity of the star. We choose a Maxwellian injection velocity with the mean value v 0 ≈ 0.1c and the gravitational acceleration g = g 0 /r 2 with g 0 = 0.5R c 2 . This gives the hydrostatic scale-height This is a much thicker atmospheric layer than the magnetar would have at a surface temperature kT < ∼ 1 keV. However, it is sufficiently thin and still resolved by our numerical grid (see below). The characteristic time it takes to form the atmosphere is short, t atm ∼ h/v 0 = 0.1R /c. Throughout the simulation particles are continually injected and absorbed by the star, sustaining a steady atmosphere at t t atm . The injection rate is chosen high enough to ensure a high density at the base of the atmosphere, The density is exponentially reduced with altitude on the scale h, and steeply drops to a low value below j/ec. Therefore, in the absence of E the hydrostatic plasma is not capable of conducting the electric current j required in the twisted magnetosphere.
Where the atmospheric density n(r) falls below j/ec, electric field E is expected to develop in response to charge starvation and lift particles from the atmosphere. The thin and dense atmospheric layer merely makes plasma available, with no special injection assumptions at the stellar surface. The numerical experiment must show how the system responds to the surface shear described in Section 3.1 and whether the induced E will self-organize to conduct the magnetospheric currents that allow the twist to be implanted.

Creation of e ± pairs
If E accelerates the lifted electrons to high Lorentz factors γ > γ thr , pair creation will be ignited. In this paper, we use the simplest implementation of this process: we choose a fixed value for γ thr and let a new e ± pair be instantaneously created every time an electron (or positron) reaches γ thr . This may be a reasonable approximation for the e ± discharge near the star where B 10 13 G (Beloborodov 2013). However, it becomes poor at larger distances where the magnetic field is weak and resonantly scattered photons have lower energies.
An additional simplification in our implementation is the prescription for the energy of the created pair. We will assume that the pair takes a fixed energy ∆E from the primary particle, and shares it equally, i.e. the new e + and e − each receives ∆E/2 (including the rest mass). Total energy and momentum parallel to B is conserved in the pair creation process.
Thus, we do not track the propagation of any highenergy photons, which is significantly simpler than the discharge model of CB14 developed for pulsars. The simplified version appears adequate for the first axisymmetric PIC model of magnetars. It should be sufficient to demonstrate some basic features of plasma selforganization in response to shearing of the magnetospheric footpoints, followed by ohmic dissipation of the twist. The results may be used as a benchmark for future more advanced simulations. Future simulations will have explicitly implemented resonant scattering process , so that ∆E will be the energy of the resonantly scattered photon, which may convert to e ± with a delay. Both γ thr and ∆E will vary with the local magnetic field, see Beloborodov (2013) for a detailed discussion.

Rescaling of large numbers in the problem
Any PIC simulation must resolve the plasma skin depth λ p = c/ω p , which is a demanding condition on the computational grid, as λ p is a microscopic scale and the ratio R /λ p is huge (comparable to 10 8 in magnetars). Similar to the PIC simulations of rotation-powered pulsars, this issue is resolved by rescaling the parameters of the problem so that λ p remains much smaller than the stellar radius, λ p ∼ 10 −2 R , but becoming sufficiently large to be well resolved. This rescaling has two main implications: (1) Similar to the pulsar problem, the increased λ p implies a reduction of the energy scales (cf. CB14). In particular, the maximum voltage that can be induced in a magnetar magnetosphere is given by γ DL (Equation 6), which now becomes moderate, γ DL ∼ 10 2 . To respect the hierarchy of the energy scales 1 γ thr γ DL , a good choice for the discharge threshold in the numercial experiment is γ thr ∼ 10. Secondary pairs receive the energy ∆E, which must be a fraction of γ thr m e c 2 . We will fix ∆E = 3.5m e c 2 for all simulations presented below.
(2) The rescaling of λ p changes the lifetime of the implanted twist, as seen from the following estimate. The value of λ p = c/ω p is related to the electric current density j by Equation (3), and the characteristic value of j scales with the magnetic dipole moment of the star µ: j ∼ ψ (c/4π)(µ/R 4 ). This gives, Combining this relation with Equation (16) for the resistive evolution timescale, one obtains One can see that the rescaling of λ p to ∼ 10 −2 R reduces the resistive timescale to t ohm ∼ 10 3 (R /c) when γ thr ∼ 10. This is fortunate, as the untwisting evolution can now be observed during a reasonably long simulation. With the realistic λ p /R ∼ 10 −8 and γ thr ∼ 10 3 one would have t ohm ∼ 10 13 R /c. Another large number that should be rescaled in the simulation is the ion-to-electron mass ratio m i /m e ≈ 2 × 10 3 . We use m i /m e = 10. This rescaling is useful for two reasons: (1) The characteristic ion plasma frequency ω p,i = (4πn i e 2 /m i ) 1/2 is not very much smaller than ω p , so that ω p,i < r/c is well satsified, and (2) m i c 2 becomes comparable to γ thr m e c 2 . The latter coincidence is also expected for the real magnetar discharge.
It is also useful to evaluate the surface magnetic field B ∼ µ/R 3 , which can be expressed from Equation (21), and then estimate the characteristic gyro-frequency, where λ p corresponds to the current density supporting a twist ψ ∼ 1. One can see that the particles are very strongly magnetized, ω B ∼ 10 4 c/R , and hence expected to move along the magnetic field lines, similar to real magnetars. The characteristic gyro-frequency is also related to another important parameter of the magnetosphere -the ratio of magnetic and plasma energy densities, For real parameters of magnetars this ratio is q ∼ 10 17 . The characteristic parameters chosen in our simulations give q ∼ 10 3 . This is still very much above unity, so the magnetosphere is nearly force-free as it should be. The parameter q also determines the Lorentz factor of Alfvén waves, γ A ≈ q 1/2 . For a real magnetar, this gives γ A γ ∼ γ thr . This condition is satisfied in our rescaled numerical experiment as long as γ thr 30.
3.5. Evolving the fields and the plasma: Aperture The particle-in-cell (PIC) method provides an efficient technique to simulate plasma from first principles. The electromagnetic fields are evolved on a grid according to Maxwell equations with the source (electric current and charge density) provided by the plasma that is selfconsistently evolved in the electromagnetic field. The plasma is represented directly as a large number of individual particles. The simulation follows the motion of each particle by calculating the applied forces. The motion of the plasma particles creates electric current which is interpolated onto the grid and then used as the source term in the Maxwell equations to update the electromagnetic field. The method well describes the plasma behavior at the microscopic kinetic level as long as the plasma skin depth is well resolved by the grid and the number of particles per grid cell is much larger than one.
Our simulations are performed using the PIC code Aperture. 2 The code was originally developed for the PIC simulations of rotationally powered pulsars (CB14). The code can follow pair creation with or without explicit tracking of high-energy photons. In the present work we use the simplified implementation of pair creation (Section 3.3) and do not use the radiative transfer module. The code is fully relativistic and designed to work on curvilinear grids. This is particularly important for problems with natural spherical geometry, such as the plasma dynamics around a spherical star in a region extending far beyond the stellar radius.
The simulations presented below are done in 2.5D, which means that our grid is 2D (in the poloidal plane) but all vector quantities are fully 3D, and we solve the full Maxwell equations assuming axisymmetry. Particles in the simulation may be thought of as rings with poloidal and toroidal velocity components. We use a spherical r, θ grid with logarithmic spacing in r and uniform spacing in θ. For all of the simulations shown in this paper, the grid size is 384 × 384 and the timestep ∆t = 10 −3 R /c.
The outer boundary of the simulation box is set at r out = 30R and employs a damping condition that lets outgoing electromagnetic waves and particles escape the box, preventing reflection. We did not detect any appreciable reflection of waves from the outer boundary. Note also that most of the active (current carrying) field lines are closed well inside the box and do not cross the outer boundary.
The shear motion of the stellar surface during the twist implantation stage t < t shear = 40R /c is equivalent to imposing a tangential electric field at the boundary. The field corresponding to the surface motion with velocity v in the lab frame is given by E = −v×B/c. It corresponds to zero electric field in the comoving frame of the stellar crust, which is assumed to be an ideal conductor. This gives the following boundary condition at r = R , The initial state is a dipole field and the normal component of the magnetic field at the surface remains unchanged during the simulation.
3.6. Units A set of natural units can be defined as follows. All lengths are measured in units the stellar radius R and time is measured in R /c. The corresponding velocity unit is the speed of light c. We define the dimensionless electromagnetic field and current density as Hereafter we will use tilde to denote dimensionless quantities, e.g.r = r/R ,t = ct/R etc.
4. RESULTS In all simulations presented below the magnetic field strength at the pole of the star isB pole = 4 × 10 4 . It corresponds toω B = 4×10 4 . We focus on the simulation with γ thr = 10, as it gives the best re-scaled model of real magnetars (Section 3.4). Simulations with different γ thr are only discussed in Section 4.3.

Initial relaxation
During the initial stage of the simulationt <t shear = 40 the dipole magnetosphere is twisted by the surface shearing motion described in Section 3.1. The surface motion induces a parallel electric field E , which lifts charges from the atmospheric layer into the magnetosphere and accelerates them. The electron Lorentz factors quickly reach γ thr and e ± discharge is triggered within a single Alfvén time of the twisted field line bundle.
The e ± plasma created by the discharge screens E , and the voltage along the current loop temporarily drops, shutting down the discharge. As the created pairs are lost to the star on the light-crossing time, a charge-starved region with significant E develops again. This first happens near the equatorial plane. As a result, an equatorial gap with strong E emerges and begins to accelerate particles, sustaining the pair creation process. The gap structure and how the e ± discharge is sustained will be described in more detail in Section 4.2.
It is clear from the simulation that a magnetospheric source of pair plasma is established in the twisted magnetosphere on a timescale not much longer than the light crossing time, before the surface shearing ends at t shear . Pair creation becomes the dominant source of plasma; the extraction of particles from the atmospheric layer is only important at the initial stage igniting the e ± discharge. After the pair discharge is activated, only a small fraction of the magnetospheric current is carried by the particles lifted from the surface. In particular, we observed that less than 1% of the current is carried by the ions.
We also observed that the twist implantation at t < t shear is accompanied by excitation of Alfvén waves, which bounce back and forth along the magnetospheric field lines. 3 Similar waves were observed in FFE simulations (PBH13). The waves are damped in the magnetosphere at later times, and the initial relaxation period is followed by the gradual evolution on a much longer timescalet ohm 100. After the surface shearing stopped at t shear , the electric discharge persisted for the rest of the simulation. It continually supported the electric current in the slowly untwisting magnetosphere, and the created particles continually bombarded the star. The duration of the simulationt sim = 350 was approximately 9 times longer than t shear and comparable to the expected resistive timescale t ohm estimated in Section 2. The observed gradual evolution of the magnetospheric twist and currents on the timescale ∼ t ohm will be described in Section 4.4.

The equatorial gap
A key aspect of the discharge self-organization is how and where particles are accelerated. The simulation clearly shows the formation of a quasi-steady "gap" with a strong E concentrated around the equatorial plane ( Figure 1). The gap thickness gap is smaller than radius, and its voltage is near the threshold for e ± discharge, Particles are accelerated in the gap and most of the pair creation events happen around this region. As seen in Figure 1, the gap has a rather sharp boundary; E is screened outside it by the created e ± plasma. The drop of E across the two boundaries of the gap is sustained by the layers of positive and negative charge (±Σ above and below the equatorial plane, respectively), according to Gauss law ∇ · E = 4πρ. The charged layers are self-consistently sustained by the difference in velocities of positive and negative charges passing through them in the self-organized E . In essence, the gap is a double layer. It has been compressed toward the equatorial plane to a minimum thickness gap that is still capable of sustaining particle acceleration to γ thr . Similar to the double layer described in Section 2.1, the charge layers sandwiching the gap have the thickness comparable to the local plasma skin depth λ p (evaluated for charge density ∼ j/c) (Figure 2). The electric field in the gap is E gap ∼ 4π(j/c)λ p and its voltage is The self-regulation of the gap voltage to Φ gap ≈ Φ thr controls the gap thickness gap ∼ γ thr λ p . Unlike normal double layers, particles accelerated in the gap are not brought from outside; instead, the gap feeds itself with particles. The accelerated particles create secondary e ± of lower energies near the gap exit, and some of the secondary particles are reversed by E and accelerated toward the opposite boundary of the gap, where they create new pairs, etc.
The multiplicity of the pair plasma is defined by M = (ρ + − ρ − )/j, where ρ + and ρ − are the charge densities of the positrons and electrons, respectively. One can see in Figure 3 that M in the gap is close to 1, i.e. the gap contains the minimum amount of plasma needed to conduct the electric current. This is consistent with no screening in the gap that allows the strong E to be sustained. Pair multiplicity in other parts of the j-bundle is close to 2, just sufficient to screen E . Apparently, the discharge in the simulation is self-organized to carry the current with the minimum voltage Φ e ≈ Φ gap ≈ Φ thr and the minimum rate of pair creation. Figure 4 shows the average hydrodynamic momenta of electrons and positrons. It is apparent that both species Figure 2. Charge density in the magnetosphere, averaged in the same way as in Figure 1. Note the thin charged layers bounding the equatorial gap across the magnetic field lines. The layers extend into the inner magnetosphere along the inner boundary of the jbundle. The charged structure observed on the field lines extending tor ∼ 9 approximately corresponds to the outer boundary of the j-bundle (see Figure 6). are accelerated across the equatorial gap to the threshold Lorentz factor γ thr = 10. The move with almost speed of light in the opposite directions and make approximately equal contributions to the current density, consistent with M ≈ 1. Outside the gap, M ≈ 2 together with the charge neutrality condition n + ≈ n − implies that the current is carried by one species while the other creates the neutralizing, nearly static, background. This is indeed observed in Figure 4.
The gap voltage is not exactly steady and shows quasi-periodic "breathing" with time. This must assist the gap in reversing some of the secondary particles so that they can cross the gap and accelerate to γ thr , sustaining the pair creation cycle. Most of the accelerated particles escape the gap and get absorbed by the star. Since the magnetosphere was set up to be symmetric about the equatorial plane, the fact that the current is strongly dominated by created pairs implies symmetric bombardment of the two footprints of the j-bundle. Thus, our simulation shows two symmetric hot spots (or rather rings, due to the axial symmetry) in the northern and southern hemispheres of the star.
As discussed in BT07 and Section 2.1, the voltage Φ e in the magnetospheric circuit is purely inductive. The parallel electric field E = −c −1 ∂A/∂t is associated with the slow dissipation of B φ rather than an electrostatic potential. Note also that the dissipation rate E · j = E j is localized in the gap while the untwisting of B φ also occurs outside the gap. The re-distribution of the dissipated B φ along the j-bundle into the screened region with E ≈ 0 occurs through the Alfvén mode, which can propagate without dissipation. The Alfvén timescale t A ∼ r/c is much shorter than the untwisting timescale t ohm , and so the magnetosphere slowly evolves through the sequence of global twist equilibria of a decreasing energy E twist , even though the magnetic energy is converted to heat only near the equator.

Dependence on the threshold voltage
While the simulation with γ thr = 10 is the most adequate re-scaled version of the magnetar magnetosphere (Section 3.4), we also performed simulations with γ thr = 20, 100, and ∞ (no pair creation). All other parameters of the four simulations were identical. Figure 5 shows the evolution of the twist energy E twist in the simulations with the four different values of γ thr . An obvious trend is observed: a higher threshold voltage for discharge, eΦ thr = γ thr m e c 2 , leads to a higher dissipation rate and a shorter lifetime of the magnetic twist. When γ thr 10, the dissipation becomes so strong that it affects the initial stage of the twist implantation at t <t shear = 40, so that a substantial part of the twist amplitude (and the corresponding energy E twist ) is lost before it could be implanted.
The extreme model with γ thr = ∞ gives so strong dissipation that E twist does not reach even 10% of its target value. It is instructive to compare this simulation with the expected dissipation rate in the pair-free configuration described in Section 2.1. From equation (6), we can estimate the voltage drop of the double layer as γ DL =Φ e ∼ √W . The initial width of the j-bundle near the star isW ∼ 1. The target current density reaches  ∼ 3 × 10 4 if the twist is fully implanted. This estimate gives γ DL comparable to 200; the actual voltage in the simulation reaches somewhat higher values. The high voltage develops early during the shearing stage and results in strong dissipation, which does not allowj to approach 3 × 10 4 .
The simulation with γ thr = 100 enables the pair discharge, which buffers the voltage growth in the j-bundle and allows a stronger twist to be implanted. The simulations with γ thr = 20 and, in particular, γ thr = 10, allow almost full implantation of the target twist with small ohmic losses. The subsequent slow resistive evolution is similar in the two models, as both have Φ thr well below the double-layer voltage and sustain a long-lived discharge activity in the j-bundle. As expected, the untwisting timescale t ohm is reduced by a factor of 2 as γ thr is increased from 10 to 20 (see Equation 16).
These results unambiguously demonstrate that the energy dissipation timescale is controlled by the pair creation threshold, confirming the conclusion of BT07. In real magnetars, we expect γ thr γ DL (Section 2). Therefore, the most relevant model is the one with low γ thr = 10, which is still high enough to accelerated particles to ultra-relativistic energies and produce relativistic secondary e ± . Figure 6 shows the resistive evolution of the j-bundle. The untwisting of the magnetic field lines proceeds as anticipated in Section 2.4, through formation of a cavity j = 0 that expands from the inner magnetosphere near the equator (large flux function u). Figure 7 shows the evolution of the poloidal current j p until the end of the simulation att sim = 350. We chose to show j p /B p because this quantity is constant along the magnetic field lines (after averaging over short-timescale fluctuations), as expected in a nearly force-free magnetosphere -currents flow along the magnetic field lines. Therefore, j p /B p is a function of the magnetic field line, which we label by the parameter u = sin 2 θ (see Equation 15). Note the expansion of the region where j p = 0 toward the magnetic axis, from u ≈ 0.75 to u 0.55. Figure 8 shows the evolution of the integrated twist angle ψ defined in Equation (1). The untwisting proceeds from near the equator, where the twist angle decreases over time, but the twist angle is not simply erased, but relocated from the inner magnetosphere to the outer parts, as expected from the untwisting Equation (13).

Expanding cavity
A curious feature is observed to develop on the magnetic field lines with u around 0.22: the twist angle ψ grows and approaches 3.5 toward the end of the simulation. This feature is also seen in the current structure shown in Figures 6 and 7. The strongly twisted, narrow bundle of field lines is inflating with time and eventually opens up, causing a magnetospheric instability (cf. PBH13). Our simulation stopped right at the onset of this development, since we would like to limit our study to the quasi-steady untwisting regime. An important difference from over-twisting studied in PBH13 is that here it is not driven by excessive surface shear. Instead, it results from resistive evolution of the implanted twist while the crust is static.

DISCUSSION
We have performed the first axisymmetric particle-incell simulations of the twisted magnetospheres of magnetars. The simulations demonstrate from first principles that electric e ± discharge is self-organized in the magnetosphere to sustain the electric current j demanded by the magnetospheric twist.
The results of our numerical experiment may be summarized as follows.
1. Shear motion of the stellar surface on a timescale t shear < t ohm successfully implants a magnetic twist . Evolution of the twist magnetic energy E twist . Four simulations are shown with discharge thresholds γ thr = 10, 20, 100, ∞. We use the exact expression for E twist = (B 2 − B 2 0 )/8π dV , where B 0 is the initial dipole field. It takes into account that besides B 2 φ /8π part of the twist energy is stored in the inflated poloidal magnetic field, which becomes important when the twist amplitude ψ exceeds unity.
in the magnetosphere. The twist is supported by continual electric current due to self-organized e ± discharge.
2. Particles are accelerated along the magnetic field lines to Lorentz factors γ ≈ γ thr , just sufficient to ignite pair creation. The voltage sustaining the electric circuit, the dissipation rate, and the lifetime of the twist are all regulated by γ thr .
3. Particle acceleration is localized in a gap near the equatorial plane (Figure 1). The gap has the electric field E ∼ 4π(j/c)λ p and width gap ∼ γ thr λ p , where λ p = (m e c 3 /4πej) 1/2 is the local plasma skin-depth. The plasma density in the gap is close to the minimum value n = j/ec required to conduct the electric current. Continual e ± creation occurs near the two exits from the gap.
4. The magnetospheric current is carried by electrons and positrons created in the magnetosphere rather than electrons and ions extracted from the atmospheric layer on the stellar surface. The created particles rain onto the footprints of the j-bundle, creating two hot spots.
5. Resistive untwisting of the magnetosphere occurs on the timescale t ohm estimated in Equation (16), in agreement with theoretical expectations. The evolution proceeds as predicted in B09: a cavity with j = 0 quickly forms in the inner magnetosphere and gradually expands, erasing the remaining electric currents. 6. A curious feature was observed in the untwisting process: while the twist energy was decreasing as expected from ohmic dissipation, the twist amplitude ψ grew in a narrow bundle of field lines at the outer boundary of the twisted region. This overtwisted bundle inflated so much that it eventually opened up.
Our results confirm that the untwisting magnetospheres naturally create shrinking hot spots (footprints of the shrinking j-bundle), which have been detected in 7 transient magnetars. The evolution timescale inferred from the simulations (Equation 16) is consistent with the decay timescale observed in transient magnetars (months to years).
One unknown in the setup of our numerical experiment is the profile of the surface shear. However, basic features observed in the simulation, in particular voltage regulation through e ± discharge and the cavity expansion, should be generic and independent of the details of the twist profile. It is less clear how generic is the formation of the narrow over-twisted bundle. This could be further explored with simulations of different shear profiles. An important caveat in the simulation setup is the simplified "on the spot" prescription for pair creation, with the created e ± pair taking a significant energy fraction from the primary particle. As briefly discussed in Section 3.3, this prescription is reasonable if the twist is confined to the region of ultrastrong magnetic field near