Quantum-electrodynamic magnetic reconnection as the origin of the FRB 200428-associated X-ray burst from SGR J1935+2154

Magnetars, often under the name soft gamma-ray repeaters (SGRs) or anomalous X-ray pulsars 1,2 , are highly magnetized neutron stars 3,4 that exhibit diverse X-ray activities. Recently, a unique non-thermal X-ray burst 5–8 with cut-oﬀ energy up to 84 keV is detected and thought to be associated with the fast radio burst (FRB) 200428 in the same single explosive event from SGR J1935+2154 9,10 , as their spectra show similar feature of narrow double peaks that are emitted almost simultaneously. However, the physical origin of this FRB 200428-associated X-ray burst is still unknown yet 11 . Here, with the ﬁrst cross-scale numerical simulation in which modeling of particle acceleration by magnetic reconnections is self-consistently coupled with that of photon emission by multiple Compton scatterings, we identify that magnetic reconnection at the quantum-electrodynamic ﬁeld strength inside the magnetar magnetosphere is the much likely driving source of such FRB-associated non-thermal X-ray burst. Both its temporal and spectral features are well reproduced in our simulations by assuming the plasma magnetization parameter σ ∼ 10 2 − 10 3 in consistency with the astronomical observations. The results could preceding local ﬁeld provide a tentative clue to probe the emitting location of such X-ray bursts within the magnetosphere. Similar studies on various X-ray the of our cross-scale simulation method that self-consistently couples both Compton scattering modelings physical roles temporal/spectral behavior

The magnetic field lines disturbed by the Alfvén wave perturbation from the magnetar surface, are projected on a specific plane, showing multiple magnetic reconnection regions at a position of R ∼ 10 2 R * < RLC within the magnetar magnetosphere. The non-thermal energetic particles (electron-positron pairs) are accelerated in these reconnection regions, and then they scatter (i.e., multiple Compton scattering) the soft X-ray photons (bisque spiral arrows) emitted from the photosphere (brown disc) on a large astronomical scale with their frequency continuously upshifted, eventually resulting in the production of hard X-ray photons of the burst. The photons escaping from the scattering region are collected by the right probe plane as the observed spectra. B s = 4.4×10 13 G is the QED "critical" field ("Sauter-Schwinger" field) 20,21 , and the plasma magnetization parameter σ = B 2 /8πnm e c 2 ≫ 1, where n is the plasma density. Therefore, on the one hand, the nonlinear QED effects including photon emission and electron-positron pair creation play significant roles, on the other hand, the plasma is highly magnetized, which both lead to dramatically different dynamics of magnetic reconnection and particle acceleration here. Figure 2 shows the reconnecting field topologies with plasmoid formations in the QED magnetic reconnections for σ = 100 (2a) and σ = 10 4 (2c), respectively. The corresponding simulations with the QED calculation switched off are also carried out for comparisons, shown in 2b and 2d. For low σ = 100, pair creations are comparatively nonsignificant, thus the reconnecting field topologies and the plasmoid structures induced by the tearing instability are similar between the results with and without QED effects (comparing 2a and 2b), similar to previous studies 22, 23 . When σ is increased to 10 4 , the tearing instability is heavily suppressed without the formation of plasmoids if without the QED effects, see 2d. Thus, the particles in the vicinity of the X-point is accelerated rapidly up to the speed of light by the first-order Fermi acceleration along the current sheet, leading to a fast reconnection with reconnection rate R rec ≈ 0.16, see Extend Data Fig. 1b. However, after taking into account the QED effects, due to a large number of pair creations that supplement the plasma current sheet from the vicinity of X-points, the total number density of the reconnecting current sheet with QED effects is much higher (about 17 times higher) than that without QED effects, see 2c and 2d. As a result, with increasing of the thermal pressure, the tearing instability grows up again, leading to the formation of much larger and denser plasmoids than the low σ case, see 2c. These plasmoids moving at a comparatively slow velocity ∼ 0.2c, leading to decreasing of both the outflow velocity from X-points and the corresponding reconnection rate with R rec ≈ 0.12 (see Extend Data Fig. 1b). Furthermore, as a large number of created pairs due to the QED effects are eventually cooled down and trapped in the plasmoids, the current density j z increases, leading to a screening of the reconnecting electric field and reduction of the particle acceleration efficiency, see Extend Data Fig. 1d.
The energy spectra of the accelerated electrons driven by the above reconnections are plotted in Fig. 3a, which exhibit obviously different features. The high-energy segment of these electron spectra could be described by a nonthermal power-law distribution dN/dE ∝ E k with a spectral index k at the energy range of E min < E < E c , where E c is the cut-off energy of the electron spectra. The power-law index k ranges from ≃ −1.33 to −1.30, nearly unaffected by the QED effects for the case of low σ = 100, however, it varies significantly in a range of [−1.39,−1.13] after a b c d Fig 2. Magnetic field topologies with plasmoid formations in reconnection for the low and high σ cases. a and b, The snapshot colormaps of the current density jz normalized by n0qec at ωpet = 700 for the case of low σ = 100, where a and b are those of respectively with and without the QED effect taken into account. c and d, The corresponding results for the case of high σ = 10 4 at ωpet = 400. The magnetic field lines are also shown with black solid lines in a-d.
introducing the QED effects for the case of high σ = 10 4 . The dependence of k on σ from more simulation runs is summarized in Fig. 3b. Without the QED effects, the electron energy spectra become monotonically harder when σ increases, where k gradually increases and approaches the asymptotic value −1.0. By contrast, with the QED effects taken into account, k first increases and then decreases with increasing of σ, where, specifically, for σ 10 3 , the electron energy spectra become softer with the spectral indices smaller than about −1.2. Such non-monotonic dependence of power-law index k on the magnetization parameter σ due to the QED effects determines the following unique X-ray emission features during multiple Compton scatterings. The Compton scattering between non-thermal energetic particles and the soft X-ray photons emitted from the magnetar photosphere is the most promising origin of non-thermal hard X-ray emissions in the magnetar magnetosphere [24][25][26] . Since the mean-free-path of photons, λ ≃ 1/(σ T n GJ ) is much larger than the typical scale of magnetic reconnections (∼ 1000d e ), where σ T is the Thomson cross-section, n GJ is the Goldreich-Julian density and d e is the electron skin depth, the Compton scattering process should occur repeatedly with photon frequency continuously upshifted in a much larger scale space than the local reconnection region. Due to this complex cross-scale simulation challenge, so far, no self-consistent numerical justification has been given for an explanation of the specific X-ray bursts yet. Using the above obtained electron spectra as the initial setup of non-thermal electrons, we further perform numerical simulations for the multiple Compton scatterings between the non-thermal particles and the soft X-ray photons in the magnetar magnetosphere with thick optical depth, where a large-scale multiple Compton scattering calculation module is developed and coupled with the above PIC simulations (details see Methods). Figure 4a shows the time-integrated photon spectra resulting from multiple Compton scatterings with different non-thermal electrons. For electrons from the QED magnetic reconnection with low σ (red line in Fig. 3a), a rather soft scattered photon spectrum (p ∼ −1.95 in the middle energy range) with a lower cut-off energy (< 100 keV) is obtained, see red line in Fig. 4a, where the photon spectrum is described by a power-law distribution of ǫdn/dǫ ∝ ǫ p+1 . It deviates from the observational data of the FRB 200428-associated X-ray burst, in the medium energy (ME) and high energy (HE) range. For the high σ case with QED effects (green line in Fig. 3a), a combination of larger E c = 200 MeV and photon spectral index of p ≃ −1.58 is obtained. Although it matches the observational data well in the ME range, strong QED effects limit numbers of high-energy photons, leading to a deviation in the HE range. This implies that σ cannot be very high, otherwise, the particle spectrum is getting too soft to match the observational data. Among a series of simulation trials, we find that a relatively hard electron spectral index of k = −1.2 together with E c = 200 MeV achieves a good fit to the whole photon spectrum of the observed FRB 200428-associated X-ray burst, where the reduced χ-square is χ 2 /d.o.f = 1.9. This, therefore, indicates that the magnetization parameter of the emitting region of FRB 200428-associated non-thermal X-ray burst cannot be extremely high, which is about at a magnitude of σ ∼ 10 2−3 , see Fig. 3b.
As a comparison, we also calculated the scattered photon spectra by solving the photon Fokker-Planck equation The energy spectra of the accelerated electrons driven by the magnetic reconnections in Fig. 2 and the corresponding power-law indices of them varying with different magnetization parameters σ. a, The electron energy spectra from the reconnections for respectively low and high σ cases in Fig. 2; b, The dependence of the electron spectral indices on the plasma magnetization parameter σ for the cases with (red circle) and without QED effects (black triangle).
numerically (see Methods). As shown in Extended Data Fig. 3, due to the lack of the cooling effect calculation for electrons from synchrotron radiation and Compton scattering, the scattered photon spectra are generally harder at the high energy range than those of PIC simulations. And they also deviate significantly from the observational data at the low energy range due to overestimation of the escaping thermal photons. This further manifests that PIC simulations for the multiple Compton scattering process are crucial to understanding the FRB 200428-associated non-thermal X-ray burst from SGR J1935+2154. Finally, the two hard peaks feature of the FRB 200428-associated X-ray burst is also well reproduced in our magnetic reconnection and multiple Compton scattering scenario, as shown in Fig. 4b-3d. The simulated burst lightcurves consist of two components 5 . One is the broad component with the duration of ∆t b = 0.2 s and exists in all energy ranges, see dashed blue lines in 3b-3d, which comes from the black body emission from the magnetar surface and is temporally described by a Gaussian function. The detailed parameters of the broad component can be seen in Methods. The other component, more notable in the higher energy range, is the two narrow peaks with each duration of ∼ ms, which is also well reproduced through our simulations, see the red lines in 3b-3d. The superposition of the two components can generally match the observed lightcurves, where the corresponding luminosity of our simulation is about L ≃ 1.6 × 10 40 erg s −1 close to the observational data with the value of 10 40 erg s −1 and our theoretical estimation for the X-ray luminosity (see Methods).
In summary, we have demonstrated that the FRB 200428-associated X-ray burst originates from the QED magnetic reconnection in the magnetosphere of SGR J1935+2154, which leads to the production of non-thermal particles and then the emission of hard X-ray photons through multiple Compton scattering between these particles and soft-X-ray photons from magnetar photosphere in a large-scale space. Both the temporal and spectral features of this burst are well reproduced in our simulations, implying that the local plasma magnetization parameter σ is at the order of 10 3 , which makes it distinguished from other preceding bursts. Given that the plasma magnetization parameter σ = B 2 /8πnm e c 2 is critically determined by the local magnetic field strength, our results provide a tentative clue to probe the emitting location of such X-ray bursts within the magnetosphere. Similar studies on various X-ray bursts of magnetars from the perspective of our cross-scale simulation method that self-consistently couples both PIC and multiple Compton scattering modelings are encouraged. More miscellaneous physical processes, like turbulence and absorption, which might play roles in depicting the temporal/spectral behavior of the bursts, remain to be considered in our future work.

METHODS
Our cross-scale numerical simulations self-consistently combine modelings of particle acceleration by magnetic reconnection with that of hard X-ray photon emissions by multiple Compton scattering. First, particle acceleration and magnetic reconnection are modeled by particle-in-cell (PIC) simulations in the quantum-electrodynamic (QED) regime, from which we obtain the non-thermal particle spectra. Then, we apply these non-thermal electrons as the initial setup for the multiple Compton scattering calculations on a large astronomical scale. Besides, we also numerically solve the photon Fokker-Planck equations and estimate the hard X-ray luminosity for comparisons with our simulation results.
PIC simulations for QED magnetic reconnection and particle acceleration. The two-dimensional (2D) PIC simulations start from a strong magnetized force-free current sheet configuration, which can refer to previous publications 22,23 by using the code EPOCH 27 . The simulation box is composed of 1600 × 1600 cells, constituting a domain size of Lx × Ly = 400d e0 × 400d e0 and the time step for each loop is set to be ∆t = 10 −3 √ σω −1 pe , where ωpe = 4πneq 2 e /me is the electron plasma frequency. Periodic and conducting boundary conditions are used in the x-direction and y-direction, respectively. The initial anti-parallel magnetic field configuration in the force-free current sheet is written as and the initial drifting relativistic pair plasma density is written as where j = e stands for electron and j = p for positron, respectively. Γ(r) = 1/( 1 − β 2 j ) = 1 + σ/[κ 2 cosh 2 (y/w)] is the Lorentz factor of the bulk drifting velocity β j along the x and z-direction, where w = κd e0 is the half-width of current sheet and d e0 = c/ωpe is the skin depth of electron. The above initial density n j (r), anti-parallel magnetic field B and the bulk drifting velocity β j (r) satisfy the Ampère law ∇ × B = µ 0 c j=e,p q j n j (r)β j (r). (3) Since the twisted magnetic field may be regarded as a multipolar component near the magnetar surface, we adopt a relatively high magnetic field strength of B 0 ∼ 10 12 G (q ∼ −1) here. According to the blackbody temperature analyzed from the observation 5 , we assume that the pair plasma with the initial temperature Te = Tp = 1.2 keV fills the simulation domain, where 50 quasiparticles per cell are used for both electrons and positrons, respectively. The QED module has been successfully applied in the PIC simulations 28,29 , where both the synchrotron radiation and Breit-Wheeler (BW) pair creation processes 30 are included. In order to explore the QED effects, the plasma magnetization parameter σ is taken from 10 2 to 4 × 10 4 for comparisons. The half-width of the current sheet is fixed as κ = 20.
We present detailed analyses on evolutions of energy conversation, magnetic flux (in the x − y plane), and energy dissipation during the simulation of QED magnetic reconnection in Extend Data Fig. 1. In the low σ cases (σ = 100), about 8% of magnetic energy is converted into plasma for both cases with and without the QED effects (Extend Data Fig. 1a). Meanwhile, a slow reconnection rate Rrec ≃ 0.08 and low dissipation rate D = j · E ≈ 0.08 are shown in the Extend Data Fig. 1b and Extend Data Fig. 1c, which are almost not affected by the QED effects. In the high σ cases (σ = 10 4 ), the energy conversion increases up to 16.13% for the case without QED effects and to 11.47% for that with QED effects, respectively, where a big difference exists between them. The reconnection rate, although, also shows an increasing from low σ to high σ cases, it is clearly suppressed from Rrec ≃ 0.16 down to 0.11 when the QED effects are taken into account (see Extend Data Fig. 1b). Besides, the mean dissipation rate also increases with σ but decreases caused by the QED effects (see Extend Data Fig. 1c). The details of the mean value of jz and Ez evolutions of high σ cases are presented in the inset panel d of Extend Data Fig. 1, showing that plenty of generated pairs supplement the current jz and the reconnection electric field Ez is suppressed at a lower level since the early stage of reconnection ωpet ≃ 140. These results manifest that the reconnection acceleration was strongly suppressed by the QED effects.
Numerical simulations for large-scale multiple Compton scatterings. A Compton scattering calculation module is implemented in our PIC code for simulations of the multiple Compton scattering process, which is similar to that in refs. 31,32 . The non-thermal electron spectra obtained in the above QED magnetic reconnection simulations, as parameterized by power-law index k and cut-off energy Ec, are used as the initial conditions for the non-thermal particle setup in our multiple Compton scattering simulation in the large astronomical spatial and temporal scales. The computational domain covers a region of 200R * × 200R * with 24 × 24 grids, where the boundary is periodic for pair particles and electromagnetic fields but open for photons. In the simulations, the multiple Compton scattering between photons and pairs is calculated by neglecting the plasma state, thus resolving the skin depth or the gyroradius of plasma is unnecessary, which allows us to run such astronomical scale simulations. The region is initially filled with background pair plasmas with a uniform number density of n = 10 5 n GJ , where 10 5 is a typical equilibrium factor 33 , and a temperature of k B T = 1.2 keV with 2048 quasiparticles per cell. The non-thermal pair particles are excited from the background plasma in all grids with the excitation rate following a time profile of Vexc(t) = G(0.0, 0.15, 0.003) + G(0.0, 0.019, 0.003) (see Extend Data Fig.2a), where G is a Gaussian function described by the starting time, the peak time, and the time width at half height. The time interval of 40 ms in Vexc(t) represents the observed separation of ∼ 30 ms in the two hard peaks of the burst from SGR J1935+2154 is represented. The soft X-ray photons are continuously injected from the left boundary and the Compton scattered photons are collected at the right boundary (see Fig.1). Note that resonant cyclotron scattering should be taken into account for radiations in the magnetosphere. However, on the one hand, it strongly depends on the angle between the magnetic field and photon motion, which requires a very fine resolution of grids and a sufficient number of particles. On the other hand, because of the complicated turbulent magnetic field on a small spatial scale, the resonance scattering effect may equivalently contribute to the increase of scatter cross-section for photons of different energy ranges. Thus, for simplicity, we only consider the large-scale magnetic field of B 0 in the simulations, and attribute the resonant scattering effect to modification of the Compton scattering cross-section by multiplying a factor of g ∼ 10 3 .
At first, it takes about 0.1 s for the thermal electrons and background soft X-ray photons to achieve the equilibrium state between them. Then, when the non-thermal particles are excited from plasma, the Compton scatterings between them and the soft X-ray photons occur repeatedly with the photon frequency continuously upshifted from soft (∼ keV) to hard (∼ 100 − 1000 keV), as shown in Extend Data Fig.2, forming the X-ray burst. After the burst, at a time of 0.4 s, the equilibrium between pair particles and photons achieves again.
Numerical solving of photon Fokker-Planck equation. The photon Fokker-Planck equation that describes the time evolution of the photon spectrum scattered by a power-law distributed electrons is written as where ǫ is the photon energy, n ′ = ǫ 2 n, k B is the Boltzmann constant, Q and S are the escaping and source terms respectively. Here, the thermal soft X-ray photons simultaneously interact with the thermal background and power-law components of pair particles. Two characteristic timescales for the thermal Comptonization and the interaction with power-law distributed particles can be written as where N is the pair number density and η = Ec/E min . If we neglect the escaping and source terms, and take δ(t) ≡ 0, Eq. (4) will reduce to the Kampaneets equation. For the case of δ(t) > 0, it means that the specific non-thermal pairs are excited in the systems. Similar to the setup of PIC simulations, we assume δ(t) has the same time profile as the above excitation rate Vexc(t). Furthermore, considering the soft X-ray background and the escape of photons, we take Q = n ′ /τesc and S = ǫ 2 exp(−ǫ/k B T )/τexc, where τesc = τexc = 5t c,pl . Eq. (4) is then solved by Chang-Cooper method 34 using the code Pychangcooper. The numerical solutions of the photon Fokker-Planck equation with the excitation of a power-law distributed electrons are presented in Extend Data Fig. 3. By changing the excited particle power-law index k from −2.2 to −1.2 with different cut-off energy Ec = 20 or 200 MeV, the results show that in the middle energy, the spectral indices are consistent with the observation data. However, due to the lack of cooling effect from synchrotron radiation and Compton scattering, the scattered photons are generally harder at the high energy range than those of PIC simulations. And they also deviate significantly from the observational data at the low energy range due to an overestimation of the escaping thermal photons.
Estimation for the X-ray luminosity. We consider the non-thermal accelerated electrons in QED reconnections have a power-law distribution f (γe) = A 0 γ k e , in which the maximum energy can be estimated as γmax ≈ √ σmec 2 . According to the energy conservation during reconnections, we have where α is the averaged particle injection ratio, Erec ≃ 0.1cB 0 is the reconnection electric field, l is the scale of acceleration region, and V = 4πr 2 ∆r is the volume of the system. Taking the background photons as the blackbody radiation from a surface with temperature of T * ∼ 10 5 K, the Compton scattering luminosity between the accelerated pair particles and photons is where U ph = 4σ sb T 4 * (R/R * ) −2 /c is the background photon energy density, f b is the beaming factor of the emission region, and σ sb is the Stefan-Boltzmann constant. Combining Eqs. (7) and (8), the luminosity can be estimated as which is consistent with the observed value of the FRB 200428-associated X-ray burst from SGR J1935+2154.

DATA AVAILABILITY
All relevant data except the Insight-HXMT data are available from the authors on request.   Fig. 2d (normalized by B0de0). The mean reconnection rate Rrec = |∆Φ x(y) /∆t| is calculated by the magnetic flux change within a certain period of time (normalized by cB0). c, Time evolutions of the dissipation rate D = j · E around the current sheet for different cases. The dissipation rate D is an average value calculated in a rectangle around the current sheet (|∆x| < 100de0 and |∆y| < 30de0) (normalized by nqeΓβzB0c 2 ). d, Time evolutions of the averaged jz and Ez around the current sheet in the high σ = 10 4 cases. To show the reconnection electric field more clearly, jz is normalized by nqeΓβc and Ez is normalized by 0.2B0c, respectively.