The Acceleration and Confinement of Energetic Electrons by a Termination Shock in a Magnetic Trap: An Explanation for Nonthermal Loop-top Sources during Solar Flares

Nonthermal loop-top sources in solar flares are the most prominent observational signatures that suggest energy release and particle acceleration in the solar corona. Although several scenarios for particle acceleration have been proposed, the origin of the loop-top sources remains unclear. Here we present a model that combines a large-scale magnetohydrodynamic simulation of a two-ribbon flare with a particle acceleration and transport model for investigating electron acceleration by a fast-mode termination shock (TS) at the loop top. Our model provides spatially resolved electron distribution that evolves in response to the dynamic flare geometry. We find a concave-downward magnetic structure located below the flare TS, induced by the fast reconnection downflows. It acts as a magnetic trap to confine the electrons at the loop top for an extended period of time. The electrons are energized significantly as they cross the shock front, and eventually build up a power-law energy spectrum extending to hundreds of kiloelectron volts. We suggest that this particle acceleration and transport scenario driven by a flare TS is a viable interpretation for the observed nonthermal loop-top sources.


Introduction
Solar flares are the most powerful energy release phenomena and important sites for particle acceleration in the solar system (Benz 2017). It is believed that the stored magnetic energy explosively releases via magnetic reconnection (Shibata & Magara 2011). Observations have shown that a significant fraction of the released magnetic energy goes into the accelerated nonthermal particles (e.g., Emslie et al. 2012;Aschwanden et al. 2017). However, it remains controversial how a large number of particles (>10 36 electrons) are impulsively accelerated to high energies within tens of seconds to minutes (Miller et al. 1997). Proposed acceleration mechanisms include acceleration in the reconnection layers (Drake et al. 2006;Oka et al. 2010;Li et al. 2018aLi et al. , 2018b, stochastic acceleration by plasma turbulence (Miller et al. 1996;Petrosian & Liu 2004;Lazarian et al. 2012;Pongkitiwanichakul & Chandran 2014), and shock acceleration in the reconnection outflows (Mann et al. 2009;Guo & Giacalone 2012;Li et al. 2013a).
The most remarkable evidence for flare particle acceleration is a hard X-ray (HXR) emission source above the top of soft X-ray loops first reported by Masuda et al. (1994). Similar events with high-energy emissions have been observed since then (e.g., Melnikov et al. 2002;Krucker et al. 2010;Liu et al. 2013;Krucker & Battaglia 2014;Oka et al. 2015;Gary et al. 2018). The loop-top or above-the-loop-top (referred to as "loop-top" hereafter) HXR and microwave sources indicate energy release and particle acceleration in the corona with energetic electrons typically taking a power-law energy distribution. In addition, the confinement of energetic electrons at the loop top is another issue for the formation of coronal HXR sources. Simões & Kontar (2013) found that the number of energetic electrons required to explain observations is 2-8 times higher at the loop top than at the footpoints, indicating electron trapping at the loop top. A successful model for the particle acceleration and transport in solar flares must be able to explain these important observations for the loop-top sources.
A fast-mode termination shock (TS) has been proposed to explain electron energization at the loop-top region (Masuda et al. 1994;Tsuneta & Naito 1998). The formation of a TS has also been predicted in numerical magnetohydrodynamic (MHD) simulations when high-speed reconnection outflows impinge on magnetic loops (e.g., Forbes 1986;Magara et al. 1996;Yokoyama & Shibata 1998, 2001Takasao et al. 2015;Shen et al. 2018). However, the existence of a TS and its role in particle acceleration have remained uncertain due to the lack of evidence in observations. Recently, using the Karl G. Jansky Very Large Array, Chen et al. (2015) revealed the presence of a TS in an eruptive flare based on radio spectroscopic imaging and provided observational evidence for its role in accelerating electrons to at least tens of kiloelectron volts. They showed that the TS, which manifests as stochastic radio spikes, was located at the front of reconnection downflows and was slightly above the loop-top HXR source.
Despite promising development in observations and theories of loop-top sources, to our best knowledge no study has successfully modeled the spatial distribution of accelerated electrons for explaining high-energy emissions at the loop-top region. As discussed above, this may require both acceleration and confinement of electrons at the loop-top. Some models such as magnetic mirroring, turbulent pitch-angle scattering, and formation of thermal fronts have been suggested to explain the confinement of energetic electrons within the flare loops (e.g., Li et al. 2013b;Simões & Kontar 2013;Kontar et al. 2014;Musset et al. 2018;Sun et al. 2019). However, it remains unclear how energetic electrons are confined at the loop-top region. In previous models, the TS is often considered as a planar standing shock (e.g., Tsuneta & Naito 1998;Mann et al. 2009). Recent MHD simulations have revealed that the TS can be very dynamic and have complex structures (e.g., Takasao et al. 2015;Takasao & Shibata 2016;Takahashi et al. 2017;Shen et al. 2018). The magnetic structures in the TS region may affect the acceleration and transport of particles near the looptop region. Especially, a concave-downward magnetic structure is found below the TS, which may trap electrons at the loop top, as it is more difficult for particles to travel transverse to the magnetic field than along it Kong et al. 2015Kong et al. , 2016. Although such a magnetic configuration has been shown in earlier MHD simulations (e.g., Magara et al. 1996), its role in confining energetic electrons has not been investigated heretofore.
In this study, we numerically model the acceleration of energetic electrons at the flare TS based on MHD simulations of a two-ribbon solar flare and emphasize the importance of a concave-downward magnetic trap structure in the TS region as a confinement mechanism for the formation of loop-top HXR sources. We find that the accelerated electrons are concentrated in the loop-top region due to acceleration at the TS and confinement by the magnetic trap structure. As far as we know, this is the first model that reproduces the necessary electron acceleration and spatial distribution for the loop-top sources. Section 2 introduces our MHD simulations and particle acceleration and transport model. In Section 3, we present the simulation results, and examine the electron acceleration and confinement. In Section 4, we discuss the conclusions and implications of this work.

MHD Simulations of the TS
We model the reconnection-driven TS in a classic tworibbon flare configuration by solving the two-dimensional resistive MHD equations using the Athena MHD code (Stone et al. 2008). The initial setup is a vertical Harris-type current sheet along the y-direction in mechanical and thermal equilibrium. We employ an initial magnetic field perturbation to speed up the reconnection onset. We use a line-tied boundary condition at the bottom and open conditions at other boundaries so the magnetic field and plasma can move into or out of the simulation domain freely. The heat conduction is not included and the specific heat ratio is 5/3. We use a uniform resistivity that corresponds to a Lundquist number S=10 5 . The simulation box is x=[−1, 1] and y=[0, 2], and 2048×2048 Cartesian grids are uniformly spaced. The simulations are normalized by the length L 0 =75 Mm, the velocity V 0 =810 km s −1 , and the time t 0 =92 s. The parameters and setup in the MHD simulation are the same as those of Case 1 in Shen et al. (2018).

Modeling Electron Acceleration by Solving the Parker Transport Equation
In this study we focus on a period (96.5-97.5t 0 ) where the TS is nearly steady and symmetric (Shen et al. 2018). The acceleration and transport process during the dynamical evolution phase of the shock will be discussed in a future publication. We select the region x=[−0.25, 0.25] and y=[0.2, 0.7] as the simulation domain for particle acceleration.
Electron acceleration and transport are modeled by numerically solving the Parker's transport equation (Parker 1965), is the particle distribution function dependent on position x i , momentum p, and time t; κ ij is the spatial diffusion coefficient tensor, U i is the bulk plasma velocity, and Q is the source. The equation is solved by particle-based stochastic differential equations (Giacalone & Neugebauer 2008;Kong et al. 2017Kong et al. , 2019Li et al. 2018b), and the fluid velocity and magnetic field necessary for solving the equations are from MHD simulations. The temporal cadence of MHD frames is 0.01 t 0 and no interpolation is applied in time, meaning that we assume a steady TS between adjacent MHD frames. We use a bilinear interpolation in space to deduce the physical quantities at the particle position. For the source term Q, we inject a large number of pseudo-particles with the same initial energy of 5 keV in the upstream of the shock. The spatial diffusion coefficient describes particle transport in the magnetic field. The diffusion coefficient tensor is, where κ P and κ ⊥ are the parallel and perpendicular diffusion coefficients, and B i is the average magnetic field vector. The antisymmetric diffusion coefficient related to particle drifts is neglected because the gradient and curvature drifts are in the out-of-plane direction. κ P can be calculated from the quasilinear theory (Jokipii 1971). We assume the magnetic turbulence is well developed and has a Kolmogorov power spectrum P∝k −5/3 , then the resulting diffusion coefficient κ P ∝p 4/3 when the particle gyroradius is much smaller than the turbulence correlation length. We use the following expression for κ P (Giacalone & Jokipii 1999), where v is the particle speed, L c is the turbulence correlation length, σ 2 is the normalized wave variance of turbulence, and Ω is the particle gyrofrequency. The normalization of the diffusion coefficient is κ 0 =L 0 V 0 =6×10 17 cm 2 s −1 . We compression ratio X, and shock angle θ Bn at the TS front from 95 t 0 to 100 t 0 . The shaded period is selected to perform particle simulation. assume the turbulence correlation length L c =1 Mm, the average magnetic field B=100 G, and the normalized wave variance of turbulence σ 2 = d á ñ B 2 /B 0 2 =0.6. Then, κ P0 =3×10 15 cm 2 s −1 for the electron initial energy E 0 =5 keV, corresponding to 0.005 κ 0 . Here we take κ ⊥ /κ P =0.1 similar to results of test-particle simulations in synthetic turbulence (Giacalone & Jokipii 1999). Other kinetic processes that may increase the pitch-angle scattering rate is not included in the simple model. As shown in Shen et al. (2018), the TS can be resolved in several cells in MHD simulations, which means that the shock width is on the order of one cell, ∼0.001 L 0 . Therefore the characteristic diffusion length at the shock L d =κ nn /V sh , where κ nn is the diffusion coefficient in the shock normal direction, is roughly the shock width for a quasi-perpendicular TS. With these parameters, the electron energy distribution resembles a power-law shape, close to that predicted by the diffusive shock acceleration (DSA) theory. Figures 1(a) and (b) show the distributions of y-component of plasma flow velocity V y and the divergence of plasma velocity ∇·V at time 96.5 t 0 , respectively. A TS forms when the fast reconnection outflow collides with closed loops and can be well manifested by negative ∇·V regions located at the front of reconnection outflow, as marked by the arrow in Figure 1(b). The TS is very dynamic, with the morphology and physical quantities varying in space and time (Shen et al. 2018). Figures 1(c)-(e) show the temporal variations of the maximum and average values of fast Mach number M F , density compression ratio X, and shock angle θ Bn at the TS front from 95 t 0 to 100 t 0 , respectively. Consistent with the observational results based on the split-band feature of the TS-associated stochastic spike bursts , the maximum (average) M F ranges from 1.5 to 3.2 (1.4 to 2.4), the maximum (average) X ranges from 1.6 to 3.7 (1.4 to 2.5), and the maximum (average) θ Bn ranges from 66°to 90°(12°to 82°).

Simulation Results
Here we focus on a relatively steady period between 96.5 t 0 and 97.5 t 0 and utilize the plasma velocity and magnetic field to perform particle acceleration modeling. During this period, the average density compression ratio X is ∼2 and shock angle θ Bn is 60°-70°.
In the particle acceleration simulation, a total of 2.4×10 6 pseudo-particles are injected at a constant rate upstream of the TS. The particles have the same initial energy of 5 keV. To improve statistics at high energies, we have implemented a particle-splitting technique so a pseudo-particle will be splitted into more particles at higher energy (Kong et al. 2017(Kong et al. , 2019. Particle acceleration at quasi-perpendicular shocks are usually considered as the so-called shock drift acceleration (Wu 1984;Mann et al. 2009). In a diffusion approximation, the shock drift acceleration can be included in the DSA (Jokipii 1987). For X≈2, the DSA predicts a power-law distribution in momentum ( ) 6 . Therefore the differential distribution in energy has the form = µ d - with the spectral index δ=(X+2)/[2(X−1)]≈2 in the nonrelativistic limit. Figure 2(a) shows the temporal variations of energy spectra of accelerated electrons integrated over the whole simulation domain. The low-energy spectra below 50 keV are approximately power-law distributions with a spectral index δ=2.5, close to the prediction from the 1D steady-state DSA solution. With the thin-target emission model, the photon power-law index is g d = + » 1 3.5 thin , which is generally harder than that deduced from loop-top HXR sources (Oka et al. 2018). As in the DSA the power-law index depends strongly on the compression ratio, a TS with a smaller compression ratio (e.g., X∼1.5) may produce electron spectra consistent with observations (δ≈3.5, then γ thin ≈4.5). The compression ratio can be affected by various parameters such as the plasma β, guide field, and thermal conduction in MHD simulations (e.g., Seaton & Forbes 2009;Takasao & Shibata 2016;Shen et al. 2018). As noted in Oka et al. (2018), in some flares the looptop HXR emission can extend to the γ-ray range but shows hard spectra close to our simulation results in the X-ray range below 100 keV (e.g., Pesce-Rollins et al. 2015). Figure 2(b) shows the number of electrons accelerated to different energies as a function of time. Electrons can be accelerated to 100 keV in a few seconds. In DSA the acceleration time depends mainly on the particle diffusion coefficient. The smaller the diffusion coefficient, the higher energy can be obtained for a given time. Figure 3(a) shows the distribution of ∇·V at 97.5 t 0 , when the TS is at y∼0.6 L 0 . Magnetic field lines in red illustrate the magnetic trap structure. Figures 3(b)-(d) show the spatial distributions of accelerated electrons at 97.5 t 0 for three different energy ranges, 5-10 keV, 20-30 keV, and 50-100 keV, respectively. In all energy ranges, a loop-top source, where energetic electrons are concentrated, appears below the TS and above closed loops, coinciding with the configuration of the magnetic trap structure and consistent with the observations in Chen et al. (2015). It also shows that the higher the electron energy, the smaller the source size is and the closer the source is located to the TS (in higher altitude). The loop-top sources for 20-30 keV and 50-100 keV are located ∼0.1 L 0 above closed loops, corresponding to ∼7 Mm or 10″. Those results are consistent with observations of loop-top HXR sources (Liu et al. 2008(Liu et al. , 2013Krucker et al. 2010;Oka et al. 2015). Note that in observations the dependence of source size with energy for coronal HXR sources within the loop remains controversial . We also found that the looptop source does not qualitatively change when we vary κ ⊥ /κ P . However, it does quantitatively control the efficiency of the confinement and the flux contrast between the loop-top source and the rest of the simulation region. More detailed modelobservation comparisons are required to test this model and its parameter dependence.
To further illustrate the effect of the magnetic trap structure, Figure 4 shows the trajectory of a representative pseudoparticle. The particle moves roughly along the large-scale field lines and travels back and forth as seen in the x position. It is accelerated mainly at the TS and has three rapid acceleration phases as highlighted in orange, magenta, and green, respectively. After the second acceleration phase, the electron energy has reached 70 keV and is trapped in the concave-downward magnetic structure. Later it moves upward and is guided back to the TS and receives a further acceleration to 120 keV. Finally the electron escapes from the magnetic trap. The electron can encounter the TS multiple times via largescale field lines, because the downstream flow can be very slow in the shock frame and the electron can cross field lines due to perpendicular diffusion. This suggests that the concave-downward magnetic trap structure can not only confine electrons but also help particle acceleration.
Overall, the simulation results show a concentrated nonthermal electron source at the loop top similar to what is expected from HXR and microwave observations. Figure 5 shows our picture for explaining the loop-top source by including a concave-downward magnetic trap structure in the loop-top region based on our numerical simulations. As in the Because of this process, the newly reconnected field lines and field lines at the loop-top region form a magnetic trap structure that confines energetic electrons. Interestingly, the magnetic trap is not completely closed, and electrons are still allowed to escape from the loop-top region along some field lines connecting to the footpoints.

Conclusions and Discussion
In this paper, we present an MHD-particle model for studying particle acceleration and transport by a TS in a tworibbon solar flare. For a set of parameters reasonable for solar flares, we obtain fast electron acceleration into a power-law distribution to up to several hundred kiloelectron volts. The accelerated electrons are confined in the loop-top region due to a concave-downward magnetic structure, which we refer to as a magnetic trap, forming a concentrated source as expected for producing high-energy emission in the corona. For the first time, we are able to reproduce energetic electron energy and spatial distributions that are necessary for explaining the looptop nonthermal sources.
This model can readily be coupled with radiation models for calculating HXR and microwave emissions and add more realistic effects such as Coulomb collision. We note that the particle mean free path in our simulation λ P =3κ P /v≈20 km for 5 keV electrons, which is smaller than that deduced in some flares (e.g., Kontar et al. 2014;Musset et al. 2018). However, the turbulence properties and electron mean free path at the shock region remain unknown. To achieve efficient acceleration to hundreds of kiloelectron volts, our model requires the diffusion coefficient to be relatively small at the shock. Using a larger diffusion coefficient will reduce the shock acceleration rate. It is worthwhile to note that other nondiffusive acceleration processes may happen that are not included in the DSA (Tsuneta & Naito 1998;Jokipii & Giacalone 2007;Guo & Giacalone 2010, so the acceleration results obtained here may only be a lower limit, or can be achieved by a lower turbulence amplitude assumed in the simulation. The simulations presented in this paper are for a quasi-steady and nearly symmetric phase of the TS evolution. If the TS has an important contribution to particle acceleration in solar flares, we expect that the dynamical evolution of the TS is important in modulating the acceleration of particles and the associated high-energy emissions. One piece of such observational evidence was provided by Chen et al. (2015), who showed that a temporary disruption of the TS front, revealed by radio dynamic spectroscopic imaging, coincided with a reduction of the nonthermal radio and HXR emission. Another outstanding example is the quasi-periodic pulsations (QPPs), which may be associated with modulations of the particle acceleration processes, among others (see, e.g., reviews by Melnikov 2009 andMcLaughlin et al. 2018). One of the physical scenarios, as proposed by Takasao & Shibata (2016), attributes the QPPs to the spontaneous quasi-periodic oscillation of multiple dynamic shocks (including the fast-mode TS) in the loop-top region. In the future, we will study how the evolution of the TS influences particle acceleration and distribution in the loop-top region spatially and temporally. We also note that our current MHD model does not include the effect of heat conduction. The introduction of heat conduction, as shown in previous studies (Seaton & Forbes 2009;Takasao & Shibata 2016), can affect the formation, geometry, and dynamics of the termination shock(s) and, in turn, the associated electron acceleration. It is hence of particular interest to investigate the effects of heat conduction on electron acceleration and transport by the flare TSs in the future.
To conclude, this work may have a strong implication to high-energy solar flare studies, as future development of this technique may enable detailed comparison between particle