Dynamical modulation of solar flare electron acceleration due to plasmoid-shock interactions in the looptop region

A fast-mode shock can form in the front of reconnection outflows and has been suggested as a promising site for particle acceleration in solar flares. Recent development of magnetic reconnection has shown that numerous plasmoids can be produced in a large-scale current layer. Here we investigate the dynamical modulation of electron acceleration in the looptop region when plasmoids intermittently arrive at the shock by combining magnetohydrodynamics simulations with a particle kinetic model. As plasmoids interact with the shock, the looptop region exhibits various compressible structures that modulate the production of energetic electrons. The energetic electron population varies rapidly in both time and space. The number of 5$-$10 keV electrons correlates well with the area with compression, while that of $>$50 keV electrons shows good correlation with strong compression area but only moderate correlation with shock parameters. We further examine the impacts of the first plasmoid, which marks the transition from a quasi-steady shock front to a distorted and dynamical shock. The number of energetic electrons is reduced by $\sim 20\%$ at 15$-$25 keV and nearly 40\% for 25$-$50 keV, while the number of 5$-$10 keV electrons increases. In addition, the electron energy spectrum above 10 keV evolves softer with time. We also find double or even multiple distinct sources can develop in the looptop region when the plasmoids move across the shock. Our simulations have strong implications to the interpretation of nonthermal looptop sources, as well as the commonly observed fast temporal variations in flare emissions, including the quasi-periodic pulsations.


INTRODUCTION
In solar flares, magnetic reconnection is believed to play a crucial role in the explosive release of magnetic energy in the corona (Shibata & Magara 2011;Benz 2017). Observations have shown that an enormous number of particles are accelerated to high energies (up to tens of MeV for electrons) within tens of seconds to minutes and contain a significant fraction (as high as 10%−50%) of the released flare energy (Lin & Hudson 1976;Emslie et al. 2012;Aschwanden et al. 2017). Various mechanisms may contribute to particle energization in flares, including acceleration by parallel electric field, contracting/merging magnetic islands, or large-scale compression in the reconnection layer (e.g., Drake et al. 2006;Oka et al. 2010;Zhou et al. 2015;Li et al. 2018a,b), stochastic acceleration by magnetic turbulence or plasma waves (e.g., Miller et al. 1996;Petrosian & Liu 2004;Fu et al. 2020), acceleration by the shrinkage of reconnected magnetic field lines (e.g., Somov & Kosugi 1997;Karlický & Bárta 2006), and acceleration by a fast-mode shock, often referred to as the flare termination shock (TS), driven by the reconnection outflows (e.g., Tsuneta & Naito 1998;Mann et al. 2009;Warmuth et al. 2009;Guo & Giacalone 2012;Li et al. 2013;Kong et al. 2013;Nishizuka & Shibata 2013;Chen et al. 2015;Kong et al. 2019). However, it remains controversial which process plays a dominant role and can explain various observational signatures of particle energization in solar flares (Miller et al. 1997;Aschwanden 2002;Zharkova et al. 2011).
Hard X-ray (HXR) and radio observations provide primary diagnostics of the acceleration and transport of energetic electrons in solar flares. Nonthermal looptop sources suggest that particle acceleration takes place above the top of flare loops and the TS is one of the promising candidates as the acceleration mechanism (e.g., Masuda et al. 1994;Melnikov et al. 2002;Krucker et al. 2010;Liu et al. 2008Liu et al. , 2013Krucker & Battaglia 2014;Oka et al. 2015;Gary et al. 2018). Although the TS has long been predicted in magnetohydrodynamic (MHD) simulations when fast reconnection outflows impinge upon the top of newly reconnected magnetic loops (e.g., Forbes 1986; Magara et al. 1996;Yokoyama & Shibata 1998), there is rarely solid observational evidence. One of the reasons is that the spatial size of the TS is expected to be very limited (Chen et al. 2019). In several studies, slow-drift radio emissions similar to Type II radio bursts (associated with shocks driven by coronal mass ejections) were interpreted as the radio signature of the TSs (Aurass et al. 2002;Mann et al. 2009;Warmuth et al. 2009). Recently, using the Karl G. Jansky Very Large Array, Chen et al.
(2015) presented a high-cadence radio spectroscopic imaging observation of stochastic spike bursts at decimetric wavelengths during a solar flare, which had a morphology and dynamic evolution most likely representative of a flare TS as suggested by numerical simulations. The observation suggested that the TS is located at the ending points of plasma downflows and slightly above a coronal HXR source. Chen et al. (2019) further performed detailed analysis of the split-band feature in the same radio spike burst event, and found that the high-frequency band is located slightly below the lowfrequency band which, in turn, supported the interpretation that the two bands were emitted in the downstream and upstream side of the TS, respectively. In addition to efficient particle acceleration, in some flare events, confinement of electrons in the looptop region is also required to account for the nonthermal looptop emissions (e.g., Simões & Kontar 2013). Several mechanisms have been suggested, which include magnetic mirroring and turbulent pitch-angle scattering (e.g., Simões & Kontar 2013;Kontar et al. 2014;Musset et al. 2018;Ruan et al. 2020). In theoretical models of particle acceleration in solar flares, the TS is often considered as a planar standing shock (e.g., Tsuneta & Naito 1998;Mann et al. 2009;Nishizuka & Shibata 2013).
However, recent MHD simulations have shown that the TS evolves dynamically and exhibits complex structures due to the impacts of jets/plasmoids (e.g., Takasao et al. 2015;Takasao & Shibata 2016;Takahashi et al. 2017;Shen et al. 2018;Cai et al. 2019;Zhao & Keppens 2020). Particularly, a concave-downward magnetic structure has been shown to be present below the TS in the looptop region (Takasao et al. 2015;Shen et al. 2018;Kong et al. 2019). Such a magnetic configuration is favorable for trapping electrons because it is more difficult for particles to travel transverse to the magnetic field than along it (Guo et al. 2010;Kong et al. 2015Kong et al. , 2016. In order to investigate the role of the flare TSs in electron acceleration in more detail, however, there has been a lack of studies that couple a kinetic energetic-particle model with a realistic MHD simulation of solar flare region. Recently, Kong et al. (2019) presented a model by numerically solving the Parker transport equation (Parker 1965) with the plasma velocity and magnetic field obtained from MHD simulations. They found that the electron spectrum in the low energy resembles a power-law, which agreed well with the prediction of diffusive shock acceleration (DSA) theory.
The accelerated electrons are concentrated in the looptop region due to the acceleration at the TS and confinement by the magnetic trap structure, in agreement of HXR and microwave observations. Therefore, the model in Kong et al. (2019) can have strong implication to the acceleration and confinement of electrons in the nonthermal looptop sources.
During solar flares, fast temporal variations of flare emissions, including the so-called quasi-periodic pulsations (QPPs), are commonly observed in multiple wavelengths including radio and HXRs (Nakariakov & Melnikov 2009;McLaughlin et al. 2018). In addition, radio and HXR imaging occasionally show multiple sources at or above the flare looptops (e.g., Tomczak 2001;Petrosian et al. 2002;Sui & Holman 2003;Liu et al. 2008Liu et al. , 2013Gary et al. 2018;Chen et al. 2020;Yu et al. 2020), suggesting more complicated dynamics during electron acceleration and transport. Kong et al. (2019) explored electron acceleration and distribution in a quasi-steady and nearly symmetric phase of the TS. Recent theoretical work and numerical simulations have shown that a large-scale current sheet can break into numerous plasmoids (e.g., Shibata & Tanuma 2001;Loureiro et al. 2007;Bhattacharjee et al. 2009). As shown in Shen et al. (2018), plasmoids are produced intermittently in the reconnection current sheet and interact dynamically with the TS. As a result, key properties of the TS, such as the compression ratio, Mach number and shock oblique angle, can significantly vary with time.
In this study, we investigate the dynamical modulation of electron acceleration and transport in the looptop region due to plasmoid-shock interactions. The paper is structured as follows. We briefly introduce our numerical methods in Section 2 and present detailed analysis of the simulation results in Section 3. We summarize and discuss the implications of this work in Section 4. Based on the plasma velocity and magnetic field from the MHD simulation, we model the acceleration and transport of electrons by numerically solving the Parker transport equation (Parker 1965).
It is achieved by integrating stochastic differential equations corresponding to the Fokker−Planck form of the transport equation using a large number of pseudo-particles (e.g., Zhang 1999;Giacalone & Neugebauer 2008;Guo et al. 2010;Kong et al. 2017;Li et al. 2018b). The transport of charged particles in the magnetic field is described by the spatial diffusion coefficient. The diffusion coefficient parallel to the magnetic field is calculated from the qusai-linear theory (Jokipii 1971;Giacalone & Jokipii 1999), with an energy dependence κ = κ 0 (E/E 0 ) 2/3 . Following Kong et al. (2019), we take κ 0 = 0.005 κ 0 , where κ 0 = L 0 V 0 . We also consider perpendicular diffusion similar to results of test-particle simulations in synthetic turbulence (Giacalone & Jokipii 1999) and set κ ⊥ /κ = 0.1. In MHD simulations the TS can be resolved in several cells, which means that the shock width is on the order of one grid cell, ∆x ∼0.001 L 0 . Therefore, in nearly all regions the characteristic diffusion length at the lowest energies is larger than the grid cell, κ nn /V 0 > ∆x, where κ nn is the diffusion coefficient in the shock normal direction. We set the time step ∆t = 10 −5 t 0 at the injection energy to ensure that pseudo-particles can "see" the shock transition. Kong et al. (2019) has shown that with these parameters the electron spectrum in the low-energy range resembles a power-law, close to that predicted by the DSA theory. Because the shock shape is very complex and evolves dynamically, here we inject particles uniformly in the particle simulation domain. The initial energy of electrons is fixed to be E 0 = 0.5 keV (corresponding to the typical energy of ∼6 MK plasma in the flaring region) and a total of 10 8 pseudo-particles are injected at a constant rate.

SIMULATION RESULTS
In the MHD simulation, plasmoids form intermittently at the primary X-point at y ∼1.5 L 0 and move upward/downward surfing in the reconnection outflows. We examine three downward-moving plasmoids as they crash into the looptop region during 97.4-98.5 t 0 and interact with the TS. Figures   1(a)- (b) show the maps of magnetic field strength B and plasma velocity divergence ∇ · V at 97.4 t 0 .
The magnetic field is stronger in the plasmoids than that of the ambient, likely due to compression.
The TS locates at y ∼0.6 L 0 , as well manifested by negative ∇ · V regions. Shen et al. (2018) showed that the TS is a sharp transition layer, across which the flow speed decreases and the Mach number quickly drops from ∼2 to less than 1. Figure 1(c) shows the time-distance plot of ∇ · V across the TS. We find that the height of the TS varies dynamically. The average shock compression ratio X ranges from 1.7 to 2.5, and shock angle θ Bn ranges from 30 • to 75 • (see We expect that the dynamical plasmoid-shock interactions will affect both the shock properties and the acceleration of electrons. The production of energetic electrons due to DSA depends on a number of factors such as the compression ratio, the shock angle, and magnetic field configuration in the loop-top region. Here we attempt to explore how the variation of number of energetic electrons depends on the properties of the TS. From the Parker equation we can find that the rate of particle acceleration is related to the magnitude of plasma compression, dp/dt ∼ −p∇ · V , where p is the particle momentum and V We find that the number of energetic electrons exhibit rapid variations at short (<0.1 t 0 or <9 s) time scales. However, the variations during each plasmoid-shock interaction are not exactly the same. This is because various factors, such as the size of plasmoids, the TS properties, and the magnetic field configuration, can affect the acceleration and transport of electrons. In particular, the first plasmoid, denoted "P1" , is relatively small in size, which merges into the looptop region immediately after crossing the TS front. The size of the plasmoid is about 0.02 L 0 = 1.5 Mm, or about 2 on the Sun. Chen et al. (2015) showed that a temporary disruption of the TS by a fast plasma downflow with about the same size coincides with the reduction of the HXR and radio flux.
In this study, we mainly focus on examining the first plasmoid-shock interaction, which is comparable to the situation in Chen et al. (2015).
The first plasmoid-shock interaction occurs between 97.5−97.85 t 0 . As shown in Figures 2(b)-(e), at low energy of 5−10 keV, the electron number first decreases slightly and then increases gradually.
At higher energies, the numbers decrease continually. The amount of decrease is more significant at higher energies. Eventually, the number of 15−25 keV electrons is reduced by ∼ 20% and nearly 40% for 25−50 keV. Both the trends and magnitudes are generally consistent with the evolution of X-ray flux observed by RHESSI and Fermi/GBM in Chen et al. (2015). Figure 3(a) shows the energy spectra of accelerated electrons at three times as marked by black arrows in Figure 2(a). The energy spectra below 10 keV are close to a power-law with a spectral index of δ = 2.5. The spectral index is consistent with the DSA prediction: if we take X ≈ 2, a typical value for the flare TS as predicted in the MHD simulation (Shen et al. 2018) and inferred from the observed split-band feature (Chen et al.

2019)
, the DSA predicts a power-law distribution with the spectral index δ = (X + 2)/[2(X − 1)] ≈ 2 in the non-relativistic limit. We note, however, in RHESSI observations, the X-ray emission from the low-energy part of the nonthermal electron spectrum (at 10-20 keV) is usually "buried" under that from flare-heated thermal plasma (Holman et al. 2011). The energy spectra above 20 keV deviate from a simple power-law and get softer with time. Figure 3(b) shows a more detailed view of the spectra at >10 keV. At the energy range of several tens of keV, the degree of softening is ∼1. It is similar to that reported in Chen et al. (2015), but we note that in their event the nonthermal "tail" is only observed (i.e., above the background) up to ∼25 keV, and they fitted the RHESSI X-ray spectra in the 10-25 keV range with a single power-law.
In Figure 4, the first column plots the maps of ∇ · V, and the other two columns plot the spatial We now analyze the effects of the other two plasmoids, P2 and P3. In Figure 5, the first two rows show the distributions of accelerated electrons during the impact of plasmoid P2. In the course of P2 moving across the TS, a very compact source appears both at low-and high-energies. When P2 reaches the downstream region, an intense source is distributed around the TS, similar to the case at 97.56 t 0 in Figure 4. As shown in Figure 2(d) and (e), during 97.9−98.0 t 0 , the TS has a nearly perpendicular shock geometry and a high compression ratio, both of which are favorable for efficient particle acceleration. This explains why the electron numbers increase when P2 moves across the TS. We can also find another smaller source at y ∼0.55 L 0 . The distance between the two sources is ∼0.06 L 0 ∼5 Mm. Later, during 98.0−98.1 t 0 , the horizontal shock is deformed. As shown in Figure   2, both the shock compression ratio and shock angle decrease. As a result, the growth rate of electron number becomes slower, and the number of >50 keV electrons even is reduced. In the bottom row, it shows that multiple sources appear after the collision of plasmoid P3. The upper source is located near the TS. The lower source is located at ∼0.6 L 0 , which coincides with a compression region existing for nearly 0.1 t 0 , as shown in Figure 1(c).
In this study, we perform a numerical simulation of electron acceleration and transport in the flare looptop region. We explore the impacts of plasmoids that intermittently crash into the TS. As a result, the TS evolves dynamically and can be distorted and restored. We find that the energetic electron population shows a rapid variation in both time and space, with distinct differences at low and high energies. The number of low-energy (5-10 keV) electrons correlates strongly with the total area of compression regions, while the high-energy (>50 keV) electrons show good correlation with strong compression cells but moderate correlations with the shock compression ratio and shock angle.
We have focused on studying the impacts of the first plasmoid-shock interaction, when the shock front first experiences a transition from a quasi-steady state to a highly dynamic period. In solar flares, QPPs are ubiquitously observed in almost all wavelengths including HXR and radio emissions (Nakariakov & Melnikov 2009;McLaughlin et al. 2018). The typical periods range from a fraction of a second to serval minutes (e.g., Tan et al. 2010;Yuan et al. 2019;Li et al. 2020), therefore possibly related to processes in the MHD regime. In our simulations, in response to the arrival of the three plasmoids, the total electron number at different energies shows a quasi-periodic temporal variation with a time interval of tens of seconds. Although the limited number of plasmoid interactions in our simulation renders it difficult to directly associate the temporal evolution of the energetic electrons (and the nonthermal emissions) to the so-called QPPs, our results suggest that the spatial and energy distribution of the electrons in the looptop region can indeed exhibit rapid temporal evolution in response to the arrival of the plasmoids at the TS front. If the formation of the plasmoids exhibits a quasi-periodic behavior, as was demonstrated in a number of numerical studies (e.g., Shen et al. 2011), our results may provide a viable explanation for the formation of the QPPs.
We note that this scenario was previously suggested by Takasao & Shibata (2016), although they did not have particle simulations to explicitly show the temporal and spatial evolution of energetic electrons in the vicinity of the flare TS. The periodicity in our simulation is highly correlated with the production rate of plasmoids in the flare current sheet. Both in linear theory (Loureiro et al. 2007) and 2D Petrosian et al. 2002;Sui & Holman 2003;Liu et al. 2008Liu et al. , 2013Yu et al. 2020). Particularly, the double coronal X-ray sources often display energy-dependent feature, i.e., the centroid locations becoming closer at higher energies. The distance between the upper and lower sources range from a few arcseconds to tens of arcseconds and the sources can be visible for several minutes and even longer. The double X-ray sources are attributed to particle acceleration taking place both in the upward and downward reconnection outflow regions (Liu et al. 2013). MHD simulations have shown that TSs can be generated both at the flare looptop and the bottom of the flux rope (e.g., Takahashi et al. 2017;Zhao & Keppens 2020). In our work, we have only focused on the lower TS in the looptop region. The distance between the upper and lower sources is a few Mm and only lasts for a few seconds. Therefore, it cannot explain the double coronal X-ray sources as shown in Liu et al. (2013), but predicts that multiple sources can occur sporadically in the looptop region due to the moudulations in electron acceleration and transport when jets/plasmoids interact with the TS. We note that a recent study by Yu et al. (2020) does show such a double X-ray (and microwave) source in the looptop region during the post-impulsive phase of the 2017 September 10 X8.2 flare.
There has been lacking studies that couple a kinetic energetic-particle model with a realistic MHD simulation of solar flare region. We present a model by numerically solving the Parker transport equation with the plasma velocity and magnetic field from MHD simulations. Our simulations can reproduce the energy spectrum and spatial distributions of energetic electrons necessary for explaining the nonthermal looptop sources. In our future work, we will combine our MHD-particle model with radiation models to investigate the dynamical evolution in HXR and microwave emissions during solar flares.   The time in each row is the same and marked by black arrows in Figure 2(a).