Magnetar bursts due to Alfv\'{e}n wave nonlinear breakout

The most common form of magnetar activity is short X-ray bursts, with durations from milliseconds to seconds, and luminosities ranging from $10^{36}$ to $10^{43}\ {\rm erg}\,{\rm s}^{-1}$. Recently, an X-ray burst from the galactic magnetar SGR 1935+2154 was detected to be coincident with two fast radio burst (FRB) like events from the same source, providing evidence that FRBs may be linked to magnetar bursts. Using fully 3D force-free electrodynamics simulations, we show that such magnetar bursts may be produced by Alfv\'{e}n waves launched from localized magnetar quakes: a wave packet propagates to the outer magnetosphere, becomes nonlinear, and escapes the magnetosphere, forming an ultra-relativistic ejecta. The ejecta pushes open the magnetospheric field lines, creating current sheets behind it. Magnetic reconnection can happen at these current sheets, leading to plasma energization and X-ray emission. The angular size of the ejecta can be compact, $\lesssim 1$ sr if the quake launching region is small, $\lesssim 0.01$ sr at the stellar surface. We discuss implications for the FRBs and the coincident X-ray burst from SGR 1935+2154.


INTRODUCTION
Magnetars are young, strongly magnetized neutron stars with surface magnetic fields reaching B ∼ 10 14 G, beyond the quantum critical field B c = 4.4 × 10 13 G (Duncan & Thompson (1992), see Kaspi & Beloborodov (2017) for a recent review). Their quiescent X-ray luminosity is usually far larger than the spin down luminosity, so the emission is believed to be powered by dissipation of the strong magnetic field instead of rotation. They often display dramatic variability in the X-ray and soft γ-ray band. This activity includes short (milliseconds to seconds duration) "bursts" with peak X-ray luminosity ranging from 10 36 to 10 43 erg s −1 , much longer (weeks to months) "outbursts", and sometimes "giant flares" with sudden release of > 10 44 erg of energy.
The short bursts are by far the most common type of magnetar activities. Thompson & Duncan (1995) first proposed a picture for magnetar bursts: the internal field evolution could build up stress locally on the neutron star crust; the stress could become strong enough to cause mechanical failure of the crust, which leads to a sudden shift in the magnetospheric footpoints; this sends Alfvén waves into the magnetosphere, and the subsequent dissipation of the Alfvén waves in the magnetosphere could power the X-ray emission. However, the dissipation mechanism and the accompanying radiative processes are not established. The radius of burst emission is also unknown.
Magnetars have also been proposed as sources of mysterious fast radio bursts (FRBs)-the bright, millisecond-long GHz bursts detected from cosmological distances (e.g., Thornton et al. 2013). Recent detection of FRB-like bursts from a galactic magnetar, SGR 1935+2154, provides evidence that magnetars are indeed capable of producing at least some of the FRBs (The Chime/Frb Collaboration et al. 2020;Bochenek et al. 2020). The radio bursts were detected during an active period of the magnetar, and were coincident with an X-ray burst of energy ∼ 10 40 erg (Mereghetti et al. 2020;Ridnaia et al. 2021;Li et al. 2021). This detection provides a good opportunity to understand more about magnetar activity and its relation to FRBs.
Our previous work ) proposed a possible scenario for the simultaneous generation of the X-ray and radio bursts from SGR 1935+2154. Using 2D axisymmetric simulations, we showed that lowamplitude Alfvén waves from a magnetar quake may propagate to the outer magnetosphere, become nonlinear and convert to "plasmoids" (closed magnetic loops) that accelerate away from the star. An Alfvén wave packet with an energy E A ∼ 10 40 erg, as required by the energetics of the X-ray burst from SGR 1935+2154, forms freely expanding ejecta at a radius R ∼ 10 8 cm, where the wave energy exceeds the local magnetospheric energy. The ejecta pushes out the magnetospheric field lines, and a current sheet forms behind it, leading to magnetic reconnection. Such reconnection events must produce X-ray emission. The spectrum of the resulting X-ray burst was calculated by Beloborodov (2021) and showed good agreement with observations. Magnetospheric ejecta play a significant role in FRB models. They were proposed to launch blast waves in the magnetar wind (Beloborodov 2017), which are capable of emitting coherent radio waves by the synchrotron maser process (e.g. Lyubarsky 2014; Beloborodov 2017Beloborodov , 2020Plotnikov & Sironi 2019;Sironi et al. 2021). GHz waves can also be seeded by the process of magnetic reconnection triggered by magnetospheric ejecta (Lyubarsky 2019;Philippov et al. 2019;Yuan et al. 2020;Lyubarsky 2020;Mahlmann et al. 2022); these waves may be released by the ejecta when it expands to a large radius.
In the present paper, we extend the axisymmetric simulations of Yuan et al. (2020) to a full 3D model of a magnetospheric explosion from a magnetar quake. We still use the framework of force-free electrodynamics (FFE) designed for magnetically-dominated systems, such as magnetospheres of neutron stars. FFE is essentially the infinite magnetization limit of magnetohydrodynamics; the plasma is treated as a massless conducting fluid, moving under the electromagnetic stress, while providing the necessary charge and current densities. We describe the problem setup and our numerical method in §2, and present the results in §3. The results are discussed and compared with the previous 2D simulations in §4. Our main conclusions are summarized in §5.

PROBLEM SETUP AND NUMERICAL METHOD
We consider a localized star quake that exerts a twisting Alfvénic perturbation on a small patch offset from the magnetic pole on the neutron star surface. The neutron star is assumed to have a simple dipole magnetic field. We consider the parameter regime as suitable for SGR 1935+2154. The light cylinder of the magnetar is located at r LC = cP/(2π) ∼ 1.6 × 10 10 cm, where P ≈ 3.25 s is the spin period of the magnetar, and we expect the Alfvén wave to propagate in the closed magnetosphere along a flux tube that extends to a radius R ∼ 100r * ∼ 10 8 cm , where r * is the radius of the neutron star. Since R/r LC ∼ 10 −2 , the rotation induced electric field at R is where B 0 is the background magnetic field of the magnetar. On the other hand, an Alfvén wave that becomes nonlinear at R will have a wave electric field δE B 0 (R) E 0 , therefore, we can neglect the rotation of the neutron star to a good approximation in this study.
Furthermore, for the Alfvén wave to reach a radius of R ∼ 100r * , the launching region on the stellar surface should be located at a polar angle θ ∼ 0.1 with respect to the magnetic pole. For our simulation, to cover as much dynamic range as possible, we put our inner boundary at r in = 10r * . The evolution of the Alfvén wave before reaching r in should be purely linear and well described by WKB theory. The wave packet will reach a polar angle θ ∼ 0.32 at r = r in . Our simulations will then self-consistently track the evolution of the Alfvén wave beyond r in . Our detailed setup is as follows.
We introduce the Alfvén wave perturbation by twisting a small, circular region on the r = r in surface, as shown in Figure 1. The circular region is centered at a polar angle θ = 0.4 and azimuth angle φ = 0, with a radius r 1 = 0.2r in . This motion twists back and forth one foot point of a closed flux bundle, and breaks the axial symmetry of the initial dipole configuration. The twist angular velocity with respect to the twisting center has the following profile where δω 0 is the amplitude, r is the distance to the twisting center, T is the duration of the twist, and n determines the number of wave periods. The factor cos 2 (πr /2r 1 ) ensures that the perturbation smoothly goes to zero at the boundary of the circular region, while the factor sin 2 (πt/T ) allows the perturbation to gradually transition to zero at the beginning and the end. We use these smooth profiles to avoid any numerical pathology.
For our simulation domain, we employ a uniform, 3dimensional Cartesian grid, with the neutron star located at the origin. The inner boundary radius r in is typically resolved by 64 grid points (the highest resolution run uses 128 cells per r in length). At the inner boundary, we enforce the perfectly conducting boundary condition. To avoid the stair stepping at r in , we force the fields to known values inside r in with a smoothing kernel (Spitkovsky 2006). The grid covers the region 0 ≤ x ≤ 40, −20 ≤ y ≤ 20, −20 ≤ z ≤ 20 (lengths are in units of r in and times are in units of r in /c, same below). The outer boundaries of the computational grid are covered by an absorbing layer that damps the outgoing waves (e.g., Cerutti et al. 2015;Yuan et al. 2019).

1999; Blandford 2002)
with the constraints E · B = 0 and E < B (we use Heaviside-Lorentz units and set c = 1). A brief summary of the basic algorithms used by Coffee and the results from convergence tests can be found in Appendix A.

RESULTS
Let us first show the results from an example run where the initial Alfvén wave perturbation has δω 0 = 2.0, duration T = 10, and n = 4 periods. For these parameters, the initial maximum relative amplitude of the wave is δB/B ∼ 0.05 at the inner boundary r in , and the total injected energy is 3.4 × 10 −3 µ 2 r −3 in /(4π), where µ is the magnetic dipole moment of the star. We choose this initial amplitude such that the Alfvén wave packet successfully breaks out from the magnetosphere at r ∼ 10r in . We also ran a simulation with half the initial amplitude; however, the wave packet was not able to launch an ejecta in that case. We focus on the first case below.

Nonlinear evolution of the Alfvén wave
The perturbation at the inner boundary launches a torsional Alfvén wave (and a small amount of fast magnetosonic wave, see Appendix B). For each half wavelength, the current structure consists of a core, aligned or anti-aligned with the background magnetic field, surrounded by a return current sheath of finite thickness. The field lines within the flux tube perturbed by the Alfvén wave experience an alternation between clockwise and counterclockwise twisting.
The wave initially propagates along the dipole field lines and its relative amplitude grows with radius as r 3/2 . The wave packet becomes significantly nonlinear at r ∼ 10. It is no longer confined to dipole field lines, but instead moves radially outward. The wave packet pushes the dipole field lines to open up, then the stretched field lines start to reconnect near the equator, allowing the twisted field lines in the wave packet to start detaching from the dipole magnetosphere. The wave packet is therefore launched as an ejecta. We take r = 10 as the ejection radius, R ej . Figure 2 shows a few snapshots of the magnetic field and electric field on the φ = 0 plane, which cuts through the azimuthal center Figure 2. Slices of the electromagnetic field on the y = 0 (φ = 0) plane, at three different time steps. In the top row, color shows the magnetic field component perpendicular to the plane, By; in the bottom row, color shows the electric field component perpendicular to the plane, Ey. In all panels, streamlines show the in-plane magnetic field. Lengths are in units of the inner boundary radius rin and times are in units of rin/c (same below). Note that although we showed By and Ey here, the total electric field E is perpendicular to the total magnetic field B.   . Angular distribution of the electromagnetic energy in the ejecta, measured at t = 25 (left) and t = 35 (right). The perturbation electromagnetic energy density has been integrated along the radial direction between two spheres with radius r = t − T and r = t that enclose the ejecta shell. White dashed lines show the 20%, 40%, 60% and 80% containment regions.
of the wave packet. Figure 3 shows the 3D structure of the perturbation electromagnetic energy density δU and magnetic field lines within the flux bundle perturbed by the Alfvén wave. δU is defined as δU = (δB 2 + δE 2 )/2, where δB = B − B 0 , B 0 is the initial background magnetic field, and δE = E − E 0 = E. From the E y plot in Figure 2, and the tenuous, spherical shell-like structure in Figure 3, it can be seen that a low frequency fast magnetosonic wave is generated at the leading edge of the ejecta; this is the consequence of the nonlinear conversion of Alfvén mode to fast mode when the Alfvén wave propagates along curved background magnetic field lines.
As the ejecta moves out, its thickness ∆r remains the same, but it expands laterally, roughly following spherical expansion from the star, so that the solid angle spanned by the ejecta remains more or less the same. The constant thickness can be understood from the con-servation of magnetic energy and flux. Energy conservation requires B 2 r 2 ∆r = const. Since the field in the ejecta is mostly transverse (in the θ, φ directions), the flux conservation can be written as Br∆r = const. The two conditions then suggest ∆r = const and B ∝ r −1 . We will confirm this scaling relation in the following subsection. In addition, within the ejecta, each half wavelength moves slightly sideways following its twisting direction. The ejecta looks like displaced, stacked pancakes, as shown in Figure 4, where we plot the 3D structure of the B φ component at a particular time step t = 20.
After the ejection, at t = 25 when the wave packet has reached r ≈ 2R ej , we find that about half of the initial Alfvén wave energy resides in the ejecta. For the rest of the energy, a significant fraction is used to push on the background field lines, stretching them out radially. We can roughly estimate how much work is done by the Alfvén wave packet on the background magnetosphere as follows. When the ejecta has moved to a radius r > R ej , between R ej and r, the initial dipole field is stretched into a monopole-like field. This takes energy per unit solid angle For our case, taking r = 2R ej , and the total solid angle within which the field lines open up is Ω total ∼ 2, we find that the work done is This is a fixed amount regardless of the initial Alfvén wave energy. It turns out to be about a quarter of the initial Alfvén wave energy in our case. This energy is stored in the stretched field lines, and some of it is dissipated in the formed current sheet due to magnetic reconnection. Besides this, there is also some energy that follows a portion of the Alfvén wave to go back toward the southern pole of the star before the ejection happens, which starts to bounce back and forth in the inner magnetosphere and gradually gets dissipated. Figure 5 shows the angular distribution of the electromagnetic energy in the ejecta. The 60% containment region has a solid angle Ω ∼ 0.5 steradian and the 80% containment region has a solid angle Ω ∼ 1 steradian.
We note that at the inner boundary r = r in , the Alfvén wave perturbation has an angular size Ω 0 ∼ 0.12. The wave first evolves linearly along the background dipole magnetic field, so the angular size grows with radius: Ω ∼ Ω 0 r/r in . Near the ejection radius R ej ∼ 10r in , we should have Ω ∼ 1.2. After the ejection, the ejecta roughly follows spherical expansion from the star, so the angular distribution remains more or less the same. Our measured angular size is indeed consistent with this picture.

Scaling of quantities in the ejecta
To better understand how the ejecta evolves, we measure the scaling of a few quantities in the ejecta. Firstly, we look at the peak electromagnetic field in the ejecta. A convenient measure to use is the peak perturbation energy density δU . We track the location r m of the maximum δU r 2 at each output time step; we make sure that the maximum is located in the main ejecta instead of the current sheet by choosing the maximum location on a data set smoothed using a Gaussian kernel with a standard deviation of 24∆x, where ∆x = 1/128 is the grid resolution. Figure 6 shows the result. On the left panel, we show the distance r m as a function of time. This shows that the peak point is indeed consistently located on one peak/trough of the Alfvén wave packet. In fact, the velocity of the pattern deviates slightly from a purely radial direction, and the speed is indistinguishable from c. The middle panel shows δU m ≡ δU (r m ) as a function of r m . We can see that δU m is consistent with decreasing as r −2 m besides some additional dissipation. This confirms that the magnetic field in the ejecta scales approximately as r −1 , instead of r −2 as suggested by Lyutikov (2021). Now let us turn to the fluid velocity in the ejecta. We can define a velocity field v = E × B/B 2 , then E = −v × B, and the force-free equations (2)-(4) can be cast into the following form (e.g., Gruzinov 1999) This set of equations is very similar to the usual MHD equations, except that the inertia is provided by B 2 .
Here v is essentially the plasma drift velocity. In forcefree electrodynamics, the fluid velocity itself is not defined and cannot be obtained directly from the fields, as there can be an arbitrary velocity component along B. However, the drift velocity can be a good reference to gain insights into the plasma motion. In what follows, we look into the evolution of this drift velocity v in the ejecta. The corresponding Lorentz factor is γ = 1/ √ 1 − v 2 , and the proper velocity is u = γv. Figure 6 right panel shows the evolution of u and its components at the point of the maximum perturbation electromagnetic energy density. The velocity field is also smoothed using a Gaussian kernel with a standard deviation of 24∆x. It can be seen that after the ejection, namely, after t ∼ 15, the drift proper velocity grows more or less linearly, and at large distances, this drift velocity is mostly radial. Another important point to note is that, although the fluid velocity at the peak/trough of the Alfvén wave packet can be quite large, the Alfvén wave packet is still a smooth wave structure even after the ejection: fluid enters from the front of the wave, then exits from behind. Shock formation is possible at large distances; this requires taking into account the inertia of the fluid and going beyond force-free approximation. Including the fluid inertia and pressure may also change the acceleration history of the ejecta. We leave this to future studies.

Magnetic reconnection and dissipation
When the ejection happens, magnetic field lines in the sheared flux bundle and ahead of it are pushed open by the Alfvén wave packet. The left panel of Figure 7 shows the resulting current distribution in the magnetosphere. Behind each of the half-wavelength pancakes in the ejecta, there are current layers (which look more like current filaments) connecting the pancake with the closed zone. A main current sheet forms near the equatorial plane where opposite open magnetic fluxes from the northern and southern hemispheres meet. This is also seen in the third panel of Figure 2. Plasmoidmediated reconnection happens in the equatorial current sheet (middle and right panel of Figure 7), allowing the open field lines to reconnect and return to the closed initial state. From the right panel of Figure 7, we can clearly see that it is primarily the poloidal field component that is reconnecting at the equatorial current sheet.
Dissipation of the electromagnetic energy can happen at the current sheets. In our simulation, the dissipation is numerical and occurs through three channels: (1) Kreiss-Oliger dissipation that filters out the high frequency noise; (2) when E > B, E is reduced to B; (3) when E · B = 0, the component of E that is parallel to B is cut away. It turns out that most of the dissipation is accounted for by the first channel, the Kreiss-Oliger dissipation. In figure 8, we show where this dissipation is triggered in a snapshot. It can be seen that the dis- Figure 8. A slice on the φ = 0.305 plane at time t = 30, showing the energy density of the electromagnetic field U = (B 2 +E 2 )/2, the fluid proper velocity γβ calculated using the E × B drift (the proper velocity γβ has been smoothed using a Gaussian kernel with a standard deviation of 24∆x, where ∆x = 1/128 is the grid resolution), the current density, and the numerical dissipation rateUKO weighted by the local magnetic energy density. sipation is concentrated along the current sheets. This is indeed consistent with our expectation that current sheets are natural sites for energy dissipation. These are likely sites for efficient X-ray emission.
We use the numerically dissipated energy as a proxy, to provide a picture of how the X-ray light curve behaves. In particular, since the physical dissipation preferably happens at reconnection sites where the magnetic field changes direction significantly, these reconnection sites tend to have weaker magnetic field compared to other types of dissipation sites. So we usė U KO /U B as a proxy for the emissivity, where U B = B 2 /2 is the local magnetic energy density. We assume that the emission is isotropic in the fluid rest frame; we also assume that the fluid moves with the drift velocity v = E × B/B 2 , so the beaming of the received emission is affected by the fluid velocity. We calculate the sky map as a function of the observer angle and the observer time, taking into account the light travel time across the simulation box. This is done using a Monte-Carlo approach: we assign 1-2 particles per grid cell; for each particle, the emissivity is assigned to beU KO /U B in the lab frame, and the beaming direction is randomly drawn from an isotropic distribution in the fluid frame then boosted into the lab frame using the fluid drift velocity. Figure 9 left panel shows a few snapshots of the sky map, at different observer times. It can be seen that the emission first beams around the equator, then expands and moves downward. It turns out that most of this beamed emission comes from the portion of the current sheet within the ejecta, namely, the vertical current sheet that shows up in the bottom panels of Figure  8. This part moves relativistically with the ejecta; its Lorentz factor is already a few at r ∼ 2R ej , as shown in the upper right panel of Figure 8. This results in beamed X-ray emission. We see two peaks offset in φ angle in the sky map; this is because the half cycles in the Alfvén wave with different twisting directions move slightly sideways with respect to each other (Figure 4), and the emission from the current sheet within each of the half cycles also beams differently.
In Figure 9 right panel, we show the light curve at a particular observer angle, corresponding to the small white box in Figure 9 left panel. As a comparison, we also show the total dissipated energy in the simulation box, as a function of the simulation time. Although the overall dissipation happens on a time scale ∼ R ej /c, the observed light curve is much more peaked. This is also due to the relativistic effect: as the ejecta moves relativistically toward the observer, the arrival time of the emission is compressed by a factor (1 − β).

DISCUSSION
Firstly, considering the energetics, if we scale our simulation to realistic parameters of SGR 1935+2154, the stellar magnetic field is B 0 = 4.4 × 10 14 G at the pole, and the ejection radius is R ej = 10 8 cm, then the injected Alfvén wave packet has an energy E A = 1.3×10 40 erg, and the initial relative amplitude of the Alfvén wave is δB/B ∼ 10 −3 at the maximum. As a comparison, the background magnetospheric energy at R ej is roughly E bg ∼ B 2 R 3 ej /(8π) ∼ 7.7 × 10 39 erg, so the Alfvén wave can successfully break out from the magnetosphere. Another run we did with half the perturbation magnitude, thus 1/4 of the energy in the Alfvén wave packet, E A ∼ 3 × 10 39 erg, did not successfully produce an ejecta. E A must be well above E bg for the nonlinear wave packet to overcome confinement by the surrounding background field. This threshold is comparable to that found in our axisymmetric simulations , although it is somewhat lower, because the 3D Alfvén wave packet needs to push open only a portion of the magnetosphere to break out.
A few features are robust across 2D and 3D simulations. Although the initial Alfvén wave perturbations are different in the 2D and 3D models, the ejecta structure on a poloidal plane looks remarkably similar. The ejecta is mainly composed of the current carrying, twisted field of the Alfvén wave packet, plus a fast wave in front of it, generated as the initial Alfvén wave propagates along curved background field lines. After the ejection, in both 2D and 3D, the ejecta retains its radial thickness and solid angle, expands balistically from the star, and becomes a pancake-like structure at large distances. As the ejecta pushes open the magnetospheric field lines, the main current sheet is formed near the equatorial plane behind the ejecta. However, the 3D nature of the initial Alfvén wave perturbation and its subsequent evolution does produce a few new features. First, the angular distribution of the ejecta energy is not axisymmetric in the 3D model. The angular size is ultimately determined by the disturbed region on the stellar surface that launches the Alfvén wave. The wave initially evolves linearly along the background dipole field lines, its angular size growing proportional to r; after the ejection, the angular size becomes frozen. A compact wave launching region can thus produce an ejecta compact in angular size. Secondly, the current distribution shows a more complex structure in the case of the localized 3D wave launching. Besides the equatorial current sheet, there are quite a few current filaments, especially near the lateral boundary of the perturbed magnetosphere.
We also find that at the majority of the reconnecting current sheets, it is the poloidal component of the magnetic field that reconnects, not the transverse field in the Alfvén wave. Although most of the magnetic energy in the Alfvén wave initially resides in the transverse component, the wave packet gives part of its energy to the poloidal component by deforming the background magnetic field. The deformed poloidal magnetic field then reconnect and dissipates the energy.
Our simulations are carried out in the FFE limit, neglecting the plasma inertia and pressure effects. As a result, there are only two characteristic wave modes, the Alfvén mode and the fast mode, both having a group speed of c. Therefore, shocks cannot form in the FFE framework. To understand physically how the ejecta accelerates with radius, and how the shock forms as the ejecta runs into the magnetar wind, one would need to go beyond force-free approximation. This will be investigated using relativistic magnetohydrodynamics simulations in the future.
In addition, our force-free simulations cannot capture the microphysics of the dissipation processes happening at the current sheet and other locations. In reality, the reconnection physics at the current sheet is governed by kinetic plasma processes. Furthermore, due to the relatively strong magnetic field and small length scales, the resulting radiation field is compact and the plasma strongly interacts with the radiation. Photons initially emitted through synchrotron radiation can experience additional inverse Compton scatterings and photon-photon pair production; photons can also be regenerated through pair annihilation. These processes will influence the plasma dynamics and shape the emergent radiation spectrum (Beloborodov 2021). Kinetic plasma simulations including all the relevant radiative processes are needed for a complete description of the reconnection process.

CONCLUSION
We have carried out fully 3D force-free electrodynamics simulations of a localized Alfvén wave packet launched by a magnetar quake into the magnetosphere. We find that if the Alfvén wave packet propagates to a radius R and has a total energy greater than the magnetospheric energy B 2 R 3 /(8π), then the wave can become quite nonlinear and get ejected from the magnetosphere. The ejecta can carry a large portion of the initial Alfvén wave energy. The ejecta preserves its radial thickness during its expansion from the star, so it becomes a pancake-like structure. Its angular size Ω is determined by the initial Alfvén wave perturbation at the stellar surface: Ω ∼ Ω 0 R ej /r * , where Ω 0 is the perturbation solid angle at the stellar surface, and R ej is the ejection radius. The ejecta pushes open the magnetospheric field lines, creating current sheets behind it that connect back to the closed zone. Magnetic reconnection can happen at these current sheets; this will lead to plasma energization and X-ray emission. The energy source of this dissipation is the magnetic energy contained in the stretched poloidal field lines. Some of the current sheets move relativistically with the ejecta; they can produce beamed X-ray emission, and may be responsible for the sharp spikes coincident with the radio bursts from SGR 1935+2154.   uses an algorithm similar to East et al. (2015); Zrake & East (2016): we use fourth-order central finite difference stencils on a uniform Cartesian grid and a five-stage fourth-order low storage Runge-Kutta scheme for time evolution (Carpenter & Kennedy 1994). We use hyperbolic divergence cleaning (Dedner et al. 2002) to damp any violations of ∇·B = 0. To enforce the force-free condition, we explicitly remove any E by setting E → E − (E · B)B/B 2 at every time step, and whenever E > B happens, we reset E → E(B/E). We apply standard sixth order Kreiss-Oliger numerical dissipation to all hyperbolic variables to suppress high frequency noise from truncation error (Kreiss & Oliger 1973): where U represents any hyperbolic variables, KO < 1 is a constant parameter, and we use a second order stencil for the sixth order derivative. The code is parallelized and optimized to run on GPUs as well as CPUs with excellent scaling. We carried out convergence tests for the code Coffee, following procedures discussed by Mahlmann et al. (2021). We show the results from two most important tests below.

A.1. Planar Alfvén wave test
In this test, we set up a 3D Cartesian periodic box with size L × L × L, and a uniform background magnetic field B 0 along the x direction. We initialize a planar Alfvén wave, with wave vector k = (2, 0, 1)2π/L, and relative amplitude ξ = δB/B 0 = 0.1. We let the wave evolve for a long time. Due to numerical diffusion, The wave magnetic field will slowly decay with time according to where δB i is the magnetic field of the ideal wave solution in the absence of any numerical diffusion, and D is the damping rate. In force-free codes, D is related to the numerical resistivity through where k is the wave vector, η is the numerical resistivity (Mahlmann et al. 2021, and references therein). The numerical resistivity can be written in the form where R is a resolution independent numerical coefficient, L and V are the characteristic length and speed of the problem, ∆x is the grid spacing, and r is the measured order of convergence.
In this Alfvén wave test, we only change the grid resolution, namely ∆x, to measure the damping rate D and the order of convergence r. The damping rate is more conveniently measured using the total wave energy δE: We use a simulation grid with a number of N 3 points, where N is the number of cells on each side of the box, which ranges from 16 to 320 in the series of simulations. Figure 10 left panel shows the measured D from these runs. We can see that the order of convergence r is around 5 for our scheme. It turns out that the Kreiss-Oliger dissipation is one of the most important source for the numerical resistivity; it seems to determine the convergence order. The wave damping rate D also directly depends on the prefactor KO of the Kreiss-Oliger dissipation term. Figure 10 left panel shows a direct proportionality between D and KO . We use KO = 0.1 for the global simulations presented in the paper. We can also see that the wave damping rate is less than 10 −2 c/L when there are more than 16 grid points per side of the box, or more than 8 points per wavelength.

A.2. Tearing mode test
In this test, we set up a force-free current sheet similar to Mahlmann et al. (2021). Our simulation box has a length of L=2 along x and y directions, and a length of 3L=6 along z direction. The background magnetic field has the following form B 0x = B 0 tanh(z/a), and we set a = 0.1. The field is initially perturbed by where k = 2π/L is the perturbation wavenumber. We set the perturbation amplitude to be = 10 −4 . The boundary condition is periodic in x and y directions, and has zero derivative in the z direction. The growth rate of the tearing mode can be traced using the B z component, which grows exponentially with time: B z = B z (t = 0)e γt . Figure 10 right panel shows the measured tearing mode growth rate for a series of runs with different resolutions. We find that roughly γ ∝ N −1.54 , where N is the number of grid points within the current sheet thickness scale a.
In resistive MHD description of the tearing mode, the growth rate of a single k tearing mode is given by (e.g., Rembiasz et al. 2017;Mahlmann et al. 2021) where η is the resistivity, and v A is the Alfvén speed. On the other hand, the growth rate of the fastest-growing mode is (Furth et al. 1963) Suppose the growth rate we measured is the maximum growth rate, then we would obtain the relation between the resistivity and the grid resolution as η ∝ N −3.08 . This order of convergence seems to be different from what we found in §A.1. This is because in the tearing mode experiment, there can be locations where the force-free condition is violated, therefore the enforcement of force-free condition is activated and the components of the electric field that violate E < B or E·B = 0 are cut away. This will affect the actual numerical resistivity. We expect the convergence order to follow the order of the time integration in this case. The conclusion is similar to Mahlmann et al. (2021).

B. FAST WAVES LAUNCHED FROM THE INNER BOUNDARY
To understand the wave modes launched from the inner boundary, let us consider the following simplified problem. Consider an infinitely large conductor covering the space z < 0, while the region z > 0 is filled with a force-free plasma. There is a uniform magnetic field making an angle θ 0 with respect to the normal of the conductor. Without loss of generality we assume that the magnetic field lies in the xz plane, B 0 = B 0 (sin θ 0x + cos θ 0ẑ ). A circular region on the surface of the conductor is twisted with a radially dependent angular velocity Ω Ω Ω = Ω(R, t)ẑ, as shown in Figure  11. On the surface of the rotating region, a point with cylindrical coordinates (R, φ, z = 0) has the following velocity where v φ = Ω(R)R. The rotation induced electric field at this point is then = B 0 v φ (− cos θ 0R + sin θ 0 cos φẑ) Figure 11. A simplified setup to understand the wave modes launched from the perturbation on the inner boundary of the global simulation.
The magnitude of E is E = B 0 |v φ | cos 2 θ 0 + sin 2 θ 0 cos 2 φ, and for small θ 0 , |E z /E R | tan θ 0 1. Since the conductor is surrounded by a perfectly conducting plasma, immediately outside the conductor, the electric field should be continuous. To determine the nature of the modes, we carry out a local expansion of the electric field immediately outside the conductor around the point (R, φ, z = 0 + ) into force-free normal modes. Since Alfvén modes have E lying in the k-B 0 plane, while the fast modes have E perpendicular to the k-B 0 plane, we can find the component of the fast mode by projecting the electric field on to the normal of the k-B 0 plane.
We look at the E R component first. It does not have φ dependence, therefore the wave vector only hasR and z components: k = k RR + k zẑ . The unit vector along the normal of the k-B 0 plane can be written as where b 0 = B 0 /B 0 is a unit vector along the background magnetic field, k × b 0 = k R cos θ 0 sin φx + (−k R cos θ 0 cos φ + k z sin θ 0 )ŷ + k R sin θ 0 sin φẑ, |k × b 0 | = k 2 R cos 2 θ 0 − 2k R k z sin θ 0 cos θ 0 cos φ + sin 2 θ 0 (k 2 R sin φ 2 + k 2 z ) 1/2 .
Putting together the above results, we can see that if θ 0 = 0, namely, B 0 is perfectly perpendicular to the conductor surface, the launched wave mode is purely Alfvénic. For small angle θ 0 , the fast mode electric field is a factor of max(sin θ 0 k z /k R , sin 2 θ 0 ) compared to the total electric field. In our boundary condition for the global simulation, typically k z /k R 0.1, and θ 0 0.3, therefore the fast mode electric field amplitude is at most 0.1 of the total electric field, and its energy is at most 1% of the total perturbation. Furthermore, fast modes, unlike Alfvén waves, are not collimated by the field lines and therefore propagate more or less isotropically out and decrease more quickly than Alfvén waves; and in the case of our boundary perturbations, the two sides of the rotating region will create fast wave contributions that will be negatively interfering after the wave propagates far enough. Therefore, the effect of the fast waves launched from the boundary is negligible in our simulations.