# On the Origin of Hard X-Ray Emissions from the Behind-the-limb Flare on 2014 September 1

Yihong Wu, Alexis P. Rouillard, Athanasios Kouloumvakos, Rami Vainio, Alexandr N. Afanasiev, Illya Plotnikov, Ronald J. Murphy, Gottfried J. Mann, Alexander Warmuth

► To cite this version:

Yihong Wu, Alexis P. Rouillard, Athanasios Kouloumvakos, Rami Vainio, Alexandr N. Afanasiev, et al.. On the Origin of Hard X-Ray Emissions from the Behind-the-limb Flare on 2014 September 1. The Astrophysical Journal, 2021, 909, 10.3847/1538-4357/abdc20. insu-03672451

HAL Id: insu-03672451

https://insu.hal.science/insu-03672451

Submitted on 19 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# On the Origin of Hard X-Ray Emissions from the Behind-the-limb Flare on 2014 September 1

Yihong Wu<sup>1,2,3</sup>①, Alexis P. Rouillard<sup>1,2</sup>①, Athanasios Kouloumvakos<sup>1,2</sup>, Rami Vainio<sup>4</sup>①, Alexandr N. Afanasiev<sup>4</sup>①, Illya Plotnikov<sup>1</sup>①, Ronald J. Murphy<sup>5</sup>, Gottfried J. Mann<sup>3</sup>, and Alexander Warmuth<sup>3,4</sup>①

Received 2020 July 17; revised 2021 January 11; accepted 2021 January 13; published 2021 March 15

## Abstract

The origin of hard X-rays and γ-rays emitted from the solar atmosphere during occulted solar flares is still debated. The hard X-ray emissions could come from flaring loop tops rising above the limb or coronal mass ejection shock waves, two by-products of energetic solar storms. For the shock scenario to work, accelerated particles must be released on magnetic field lines rooted on the visible disk and precipitate. We present a new Monte Carlo code that computes particle acceleration at shocks propagating along large coronal magnetic loops. A first implementation of the model is carried out for the 2014 September 1 event, and the modeled electron spectra are compared with those inferred from Fermi Gamma-ray Burst Monitor (GBM) measurements. When particle diffusion processes are invoked, our model can reproduce the hard electron spectra measured by GBM nearly 10 minutes after the estimated on-disk hard X-rays appear to have ceased from the flare site.

Unified Astronomy Thesaurus concepts: Solar coronal mass ejection shocks (1997); Solar x-ray emission (1536); Solar particle emission (1517); Solar electromagnetic emission (1490)

## 1. Introduction

### 1.1. The 2014 September 1 Event

NASA's Fermi-Large Area Telescope (LAT; (Atwood et al. 2009) detected more than 30 solar eruptive events with late-phase >100 MeV γ-ray emission during its 10 year mission, among which the 2014 September 1 event was one of the largest. The γ-ray emission was generated from the decay of pions, which were likely created by interactions of high-energy protons that impacted chromospheric material on the visible disk. The NaI and BGO detectors on Fermi's second instrument, the Gamma-ray Burst Monitor (GBM; Meegan et al. 2009), also observed hard X-ray emission in the MeV range during this event, which originated from thin- or thick-target electron bremsstrahlung (Share et al. 2018).

Several studies have focused on this solar event launched 36° beyond the east limb of the Sun (e.g., Ackermann et al. 2017; Plotnikov et al. 2017; Grechnev et al. 2018; Jin et al. 2018; Petrosian 2018; Share et al. 2018; Gopalswamy et al. 2020). While the source of the >100 MeV γ-ray emission is observed on the visible disk, the origin of the hard X-ray emission is not clear. Two strong candidates are the flare loop tops via electron thin-target emission and thick-target emission resulting from electrons impacting the visible disk or perhaps a combination of both. Because the magnetic footpoints of flaring loops associated with the coronal mass ejection (CME) eruption were occulted, it was argued by Plotnikov et al. (2017) that flare-accelerated electrons were unlikely to be responsible for the measured hard X-ray emission. By means of detailed 3D coronal modeling, Plotnikov et al. (2017) showed that the onset of the hard X-ray and γ-ray emissions occurred when the CME-driven coronal shock became magnetically connected to the visible disk. This suggested that electrons and protons accelerated at the coronal shock had the means to propagate

toward the visible disk and impact the chromosphere to produce high-energy radiation. For this event, the shock is a strong candidate for particle acceleration. The moving shock-wave front could be an important electron accelerator far from the flare site, as illustrated in the studies of Krucker et al. (1999) and Rouillard et al. (2012). Krucker et al. (1999) has presented a good correlation between electron beams detected by Wind and extreme ultraviolet (EUV) waves observed by the Solar and Heliospheric Observatory (SoHO) spacecraft. Rouillard et al. (2012) argued that 0.67–3.08 MeV electrons measured near Earth during the 2011 March 21 event could have been accelerated by the CME shock. The flaring site was on the far side of the Sun as viewed from Earth (W132°).

Ackermann et al. (2017) and Grechnev et al. (2018) showed that the time history of the hard X-rays measured by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI), Fermi/GBM, and Konus-Wind matched closely that of the microwaves measured by the Radio Solar Telescope Network. Because microwaves are likely produced via synchrotron emission of relativistic electrons in the flaring loops, it was argued that the hard X-rays must also be produced in this region. Petrosian (2018) presented a detailed physics-based modeling of the event concluding that thin-target processes occurring at the loop tops could explain the hard X-ray spectra.

In contrast, Ackermann et al. (2017) found that the electrons producing the >1 GHz emission had a power-law index of about 3.1; this is also similar to the value found by Share et al. (2018) for electrons producing the bremsstrahlung in thick-target bremsstrahlung. They suggest, at face value, that the electrons emitting at microwave come from the same accelerator as those emitting via thick-target bremsstrahlung. Share et al. (2018) also argued that, due to the absence of rotational modulation in RHESSI detector 9, >20 keV hard X-ray

<footnote n="1">¹ Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse III (UPS), France; yhwyihongwu@gmail.com</footnote>

<footnote n="2">² Centre National de la Recherche Scientifique, UMR 5277, Toulouse, France</footnote>

<footnote n="3">³ Leibniz-Institut für Astrophysik Potsdam (AIP), Potsdam, Germany</footnote>

<footnote n="4">⁴ University of Turku, Turku, Finland</footnote>

<footnote n="5">⁵ Space Science Division, Naval Research Laboratory, Washington DC 20375, USA</footnote>
emission must have been distributed over a much broader region than the RHESSI source imaged in the 6–12 keV range shown in Ackermann et al. (2017). In addition, RHESSI provided only a limited coverage because it was in the South Atlantic Anomaly (SAA) from 10:55:00 to 11:11:00 UT. Furthermore, the NaI detectors on GBM measured enhanced 10–25 keV emission about 2 minutes after the hard X-ray onset observed by the Solar Assembly for X-rays (SAX) on MESSENGER (Schlemm et al. 2007). Share et al. (2018) suggested this might be due to the emergence of the flare’s coronal hard X-ray source $2 \times 10^5$ km above the east limb. Higher-energy emission detected into the MeV range started as the 10–25 keV X-ray flux from the flare decreased. It is possible that this higher-energy emission came from the visible disk. We assume that the hard X-ray source extends over tens of heliographic degrees on the visible disk, far beyond the LAT centroid estimated for this event.

We conclude that the origin of the hard X-rays is unclear and that further studies must quantify the acceleration and propagation of energetic electrons in both the flaring regions and at the shock waves driven by the emerging CME. In the present study, we search for mechanisms that could sustain electron acceleration and hard X-ray production over tens of heliographic degrees on the visible disk via thick-target bremsstrahlung in the chromosphere. For that to happen, a magnetic connection must be established between the visible disk and the candidate accelerator. We will model electron acceleration to high energies during the lateral motion of the shock as it propagated in the corona toward the visible disk. The aims of this paper are twofold: first to present a new modeling framework to address particle acceleration at a shock propagating through coronal loops, and second to combine this new model with a shock-wave reconstruction carried out for a CME event that erupted on the far side of the Sun and during which a delayed and extended hard X-ray emission was measured by a near-Earth spacecraft.

Several studies have looked at the link between the timing of the hard X-ray emission, the propagation of the EUV waves, and the onset of the type II burst emission (e.g., Klassen et al. 1999; Vršnak et al. 2006). The formation of a quasi-perpendicular shock detected as the type II burst onset is sometimes accompanied by hard X-ray enhancements, e.g., the 1993 September 27 and the 1994 July 7 events (Klassen et al. 1999), as well as the 2003 November 3 event (Vršnak et al. 2006). For these flares visible from Earth, the strong hard X-ray emissions are likely dominated by the bright sources at the chromospheric footpoints and the flaring loop tops. Any contributions to the hard X-ray emission from electrons accelerated by a concomitant CME shock would be extremely hard to isolate from the strong flare component. Hard X-ray flares with occulted chromospheric sources (i.e., not visible from Earth), such as the event studied in the present paper, are better suited to investigate other electron acceleration mechanisms that may contribute to hard X-ray emissions particularly during the late phase of a flare–CME event.

Figure 1(a) provides the time history of the arbitrarily scaled 100–300 keV hard X-ray count rate measured by Fermi GBM. The spectra from Fermi GBM were fitted using a thick-target electron–proton model to obtain the average number of >200 keV electrons per second and the power-law spectral index in each 1 minute interval, which are plotted in Figures 1(b) and (c) with $\pm 1\sigma$ statistical uncertainties. We do

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Time series of the arbitrarily scaled 100–300 keV hard X-ray count rate observed by Fermi GBM (panel (a)), the average number of >200 keV electrons (producing the hard X-ray emission) per second in each 1 minute interval with $\pm 1\sigma$ statistical uncertainties (panel (b)), and >200 keV electron (producing the hard X-ray emission) power-law spectral indices with $\pm 1\sigma$ statistical uncertainties (panel (c)). Times when the shock was intersecting FL 2305 are shown as yellow shaded areas. The shock reaches locations I–VI at 11:19:19, 11:19:31, 11:20:45, 11:22:21, 11:24:29, and 11:24:42 UT, respectively. The times that the shock crosses locations I and VI on FL 2305 are indicated by vertical solid lines. The times that the shock crosses locations II, III, IV, and V on FL 2305 are indicated by vertical dashed lines.</figcaption></figure>

not show data after 11:24:42 UT due to significant uncertainties in the estimate of the electron number and the spectral index.

Collisionless shocks with fast magnetosonic Mach numbers ($M_{fm}$) greater than 2.7 are supercritical (Mann 1995; Schwartz 1998; Marcowith et al. 2016). Supercritical shocks cannot be stabilized by ohmic dissipation alone, and they are able to accelerate particles efficiently to high energies. During the onset of the Fermi event, the area of the shock magnetically connected to the visible disk was supercritical and quasi-perpendicular. The quasi-perpendicular geometry resulted from the fact that magnetic connection occurred on the flank of the CME in a predominantly vertical magnetic field. The flank of the shock was moving at speeds in excess of $1000 \text{ km s}^{-1}$ in quiet regions of the solar atmosphere with typical magnetosonic speeds of 300–400 km s$^{-1}$. In these regions, a supercritical shock could form (Plotnikov et al. 2017). These properties provided favorable conditions for the acceleration of electrons to high energy via shock-drift acceleration (SDA; Holman & Pesses 1983; Wu 1984). In this paper, we go one step further than Plotnikov et al. (2017) and Kouloumvakos et al. (2019). We employ a test-particle simulation to compute the acceleration process of electrons under these conditions.

Low in the corona, the CME shock propagated through both open and closed magnetic field lines. During the event of interest here, detailed modeling showed that the coronal shock was predominantly connected magnetically to the chromosphere by magnetic loops during the first hour of the eruption (Plotnikov et al. 2017). Very few studies have investigated particle acceleration by shock waves along magnetic loops. Recent theoretical studies have investigated the effects of field-intensity gradient and field line curvature on particle acceleration at shocks (Sandroos & Vainio 2006). Sandroos & Vainio (2009) have computed the impact on the shock acceleration process of the curvature and expansion rate of open magnetic fields that form helmet streamers. Other studies (Kong et al. 2017)
have evaluated the efficiency of particle acceleration along streamer loops, where particles are mainly trapped and accelerated around the loop top. In this paper, we study particle acceleration along loops that connect the expanding shock to the visible disk. In particular, we study whether such geometries would allow particles accelerated at the shock to be mirrored back from the upstream footpoint on the visible disk (FP2 in Figure 6) by strong magnetic mirroring in the lower corona and chromosphere. We only consider precipitation at the footpoint in the shock's upstream in this paper. In a future study, we will study precipitation at the footpoint situated in the shock's downstream region (FP1 in Figure 6), as pointed out by Gopalswamy et al. (2020).

In order to explain the origin of the hard X-ray emissions measured by Fermi GBM on 2014 September 1, we address a number of fundamental questions: (1) how efficiently are electrons accelerated by the mechanism of SDA during the event? (2) If SDA is insufficient, then what other mechanisms or combinations of mechanisms could energize electrons to several hundred keV or even tens of MeV? (3) If we assume that electrons are indeed accelerated to high energies in the corona, what fraction of these particles can reach the solar surface to produce nonthermal emissions?

## 1.2. Electron Acceleration at Coronal Shocks

Fast-mode shocks can be treated as moving magnetic mirrors because they are associated with a magnetic field compression. In the de Hoffmann-Teller (HT) frame, the plasma flow is along the magnetic field, and thus, the motional electric field is removed (de Hoffmann & Teller 1950). In this frame, when a particle encounters a shock with a pitch angle such that it is outside of the loss cone of the shock, it is reflected rather than transmitted. The acceleration of reflected particles is due to reflection off a moving magnetic mirror, i.e., Fermi acceleration (Ball & Melrose 2001). Comparison of observed events with models of particle acceleration at interplanetary shocks verifies that when the shock is quasi-perpendicular, SDA plays a key role in particle acceleration when the level of magnetic fluctuations is low and scattering is unimportant. Past studies have argued that SDA is responsible for the shock spikes, a strong enhancement in flux perpendicular to the magnetic field direction, measured predominantly at quasi-perpendicular interplanetary shocks (Erdos & Balogh 1994). Mann et al. (2009) applied this mechanism to shock waves produced during solar flares. They computed the production rate and the power of accelerated electrons in the flare region that compared favorably to RHESSI observations of the solar event on 2003 October 28. It is also expected that SDA is responsible for the occurrence of type II solar radio bursts. Schmidt & Cairns (2012a; 2012b) computed the properties of accelerated electron beams that can drive Langmuir waves and produce radio emissions.

By means of fully kinetic particle-in-cell (PIC) simulations, Guo et al. (2014) traced the evolution of the injected electrons and confirmed that the properties of the reflected electrons match the SDA predictions. This demonstrates that a fraction of electrons indeed participates in SDA.

Magnetic fluctuations can come from photospheric convection or nanoflaring. When magnetic fluctuations are introduced, particles can be scattered back to the shock multiple times, gaining more energy with each SDA cycle, which is known as the diffusive-shock acceleration (DSA) process (Decker 1988).

We refer to SDA as simply the drift process, which is the principle acceleration mechanism in this paper. The drift process remains evident in DSA.

Under flare or shock conditions, at the onset of particle acceleration, several cycles of SDA can provide seed particles for the process of DSA. Share et al. (2018) suggested that those seed particles are accelerated by the shock to produce the late-phase >100 MeV γ-ray emission and contribute to the solar energetic particles (SEPs). Each SDA cycle can lead to a significant energy gain of electrons at quasi-perpendicular shocks due to the high projected shock speeds. And overall, fewer shock encounters are required for particles to gain significant energy at quasi-perpendicular shocks. In this first theoretical study, we focus on the relative roles of SDA and DSA for the acceleration of electrons that propagate along magnetic loops traversed by the expanding shock wave.

## 1.3. Electron Transport in the Corona

To interpret the Fermi observations, we must determine under what conditions particles are able to overcome the strong mirroring from FP2 in Figure 6 as they propagate toward the stronger magnetic fields that prevail near the solar surface. Indeed, the mirror force will strongly limit the number of energetic electrons and protons that can hit the solar surface to produce hard X-rays and γ-rays, respectively.

Hua et al. (1989) and later Murphy et al. (2007) discussed how turbulence in flare closed loops affect the ability of protons to reach the solar atmosphere in the presence of strong magnetic mirroring. The scattering of particles in pitch angle off magnetic irregularities that may exist in coronal loops could play an important role in counteracting the effect of the magnetic mirror by shifting particles into the loss cone of FP2 in Figure 6. These particles would then be able to precipitate. As already discussed, pitch-angle scattering of particles near the shock can also increase the number of shock crossings via the DSA process (e.g., Bell 1978a, 1978b; Jokipii 1982; Vainio et al. 2014; Afanasiev et al. 2018). In order to gain a more complete insight into the possible link between coronal shocks and the occurrence of high-energy radiation, it is therefore important to model the acceleration and transport of particles in terms of the pitch-angle scattering efficiency.

This paper then generally proceeds through a number of increasingly complex simulations. In Section 2, we first describe our shock-fitting technique and the method for obtaining shock parameters along magnetic loops. In Section 3, we present our new Monte Carlo code that computes particle acceleration for a single interaction with the shock (one SDA cycle). We validate the computation by comparing our results with analytical SDA calculations (Appendix A). In Section 4, we set up the new model by using the shock properties in Section 2. We compute the SDA induced by realistic shock and coronal conditions at different times and study the effect of mirroring of particles propagating toward the visible disk (Section 5). We then allow for multiple shock encounters along realistic magnetic loops (Section 6) and allow electrons to return from downstream (Section 7 and Appendix B). The next part of the paper investigates DSA (Section 8). We compare our one-loop calculation with Fermi GBM observations in Section 9.
## 2. Modeling the Corona and the Shock

In this section, we reconstruct the shock front and determine shock parameters. We also present magnetic loops connecting the shock front to the solar disk visible from Earth by employing a 3D global MHD model of the coronal magnetic field. Then, the variability of shock and ambient plasma parameters along different magnetic loops is illustrated.

### 2.1. Shock Reconstruction and the Background Corona

The 3D evolution of the expanding shock wave is determined by fitting geometrical shapes to the multipoint observations provided by the EUV and white-light observations taken by the NASA Solar Terrestrial Relation Observatory (STEREO-A), the Solar Dynamics Observatory (SDO), and SoHO. The technique is presented in Rouillard et al. (2016) and exploited in Plotnikov et al. (2017) and improved recently in Kouloumvakos et al. (2019) in order to obtain reconstructions at a higher cadence. The shock geometry is assumed to be an ellipsoidal structure. The 3D kinematic properties of the shock are discussed below.

To obtain the properties of the shock such as the Mach number, compression ratio, and its geometry, we must know the conditions upstream of the shock and thus determine the properties of the background corona into which the shock is propagating. We follow a similar procedure detailed in our past studies (Rouillard et al. 2016; Plotnikov et al. 2017; Kouloumvakos et al. 2019) by using the coronal magnetic fields, plasma density, and temperature obtained from a 3D global MHD model developed and run by Predictive Sciences Inc. The Magneto-Hydrodynamic Around a Sphere Thermodynamic (MAST) model uses magnetograms as boundary conditions at the Sun and solves the 3D MHD equations in spherical coordinates from the upper chromosphere to the upper corona with an outer boundary set at 30 solar radii. The magnetogram used here is obtained by the heliospheric imager (HMI) on board SDO on 2014 September 1 as input at the lower boundary (Lionello et al. 2009). In order to simulate coronal properties accurately, past studies (Withbroe 1988; Lionello et al. 2001) have shown that energy exchanges that occur in the low corona must be solved for explicitly. This includes heating, thermal conduction along the magnetic field, and radiative losses. The MAST model does that and also includes the effects of Alfvén waves on heating and pressure.

The 3D speed of the shock and the conditions upstream can be used to derive Mach numbers at the shock, and by solving the Rankine-Hugoniot relation (see e.g., Priest 1982; Salas 2007), the gas and magnetic compression ratios can also be obtained. From the 3D modeling of the shock and the tracing of magnetic field lines, we can infer how the shock connects magnetically to the visible disk and thus interpret how particles accelerated at the shock may stream down to the lower corona. This was discussed in some detail for that event by Plotnikov et al. (2017). We trace a large number of these magnetic field lines in Figure 2. It was found in our previous study that the onset of the high-energy emissions in hard X-rays and γ-rays measured by Fermi occurred when the shock connected to the visible disk. Moreover, Plotnikov et al. (2017) showed the shock was initially connected to magnetic loops. As particle acceleration and transport along such loops have not been studied in great detail, this is a point of focus in this article.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** All field lines that start from the visible disk and intersect the shock at some point. The east limb of the Sun as visible from Earth is indicated by a light green arc. The cyan arc shows the central meridian of the Sun.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Four field lines of particular interest. The ellipsoids show the shock locations at 11:06, 11:10, 11:14, 11:18, 11:22, and 11:26 UT. The east limb of the Sun as visible from Earth is indicated by a light green arc. The cyan arc shows the central meridian of the Sun.</figcaption></figure>

### 2.2. Shock Parameters along Selected Magnetic Loops

We model particle acceleration along a single magnetic loop to demonstrate the important physical processes at play in the corona that accelerate electrons to relativistic energies. We begin by illustrating the large variability of coronal loop properties (e.g., length and orientation) by selecting four representative field lines that intersect different regions of the shock flank. These lines are shown in Figure 3 together with the expanding shock. The plasma temperature *T*, magnetic field magnitude *B*, plasma density *n*, and the value of the magnetic field inclination to the shock normal $\theta_{Bn}$ are then determined at the shock along these four field lines as a function of time as the shock expands.

The finite resolution of both the shock-reconstruction technique and the background coronal model leads to fluctuations in the shock parameters derived along specific field lines. Such fluctuations can create spurious effects that
<figure><img src="image_6.png" /><figcaption>Figure 4. Shock parameters along the four selected field lines: FL 1487 (magenta), FL 1690 (red), FL 2305 (black), and FL 3268 (blue). Shock-normal speed $V_{sn}$ (panel (a)), projected shock speed $V_{sp}$ (panel (b)), the value of the magnetic field inclination to the shock-normal $\theta_{Bn}$ (panel (c)), fast-mode Mach number $M_{fm}$ (panel (d)), gas compression ratio $\gamma_g$ (panel (e)), and magnetic compression ratio $\gamma_B$ (panel (f)).</figcaption></figure>

must be avoided in order to model accurately the acceleration and transport of particles. For instance, in Afanasiev et al. (2018), analytical expressions are used to fit the shock parameters derived by Rouillard et al. (2016). Here we take a slightly different approach to derive the shock parameters.

The initial step to derive shock speed is similar to Rouillard et al. (2016). We consider a set of 200 × 200 grid points distributed over the shock ellipsoid. To compute the shock-normal speed $V_{sn}$, we first find the grid point $P$ on the ellipsoid at time $t$ that is closest to the shock-field line intersection at time $t$, whose index is $(i, j)$. We then find the grid point $P'$ with index $(i, j)$ on the ellipsoid at time $t - \delta t$, where $\delta t$ varies between 0.5 and 1.5 s. Next, we compute the distance between $P$ and $P'$, which we divide by the time interval $\delta t$. The speed of the propagating front along the field line is $V_{sp} = V_{sn}/\cos\theta_{Bn}$. The fast-mode speed is then obtained by

$$V_{fm} = \sqrt{\frac{1}{2}[V_A^2 + C_S^2 + \sqrt{(V_A^2 + C_S^2)^2 - 4V_A^2 C_S^2 \cos^2 \theta_{Bn}}]}, \quad (1)$$

where $V_A$ is the Alfvén speed and $C_S$ is the sound speed. The fast magnetosonic number $M_{fm}$ is defined as

$$M_{fm} = \frac{V_{sn} - V_w \cdot n}{V_{fm}}, \quad (2)$$

where $n$ is the local shock-normal vector and $V_w$ is the solar wind bulk velocity. Once the parameters of the evolving shock are determined along the four selected field lines, these parameters are spline interpolated at a cadence of $\delta t = 0.5$ s to produce a sequence of regularly time-spaced ellipsoids.

<figure><img src="image_7.png" /><figcaption>Figure 5. Ambient plasma parameters along the four selected field lines: FL 1487 (magenta), FL 1690 (red), FL 2305 (black), and FL 3268 (blue). Plasma temperature (panel (a)), magnetic field magnitude (panel (b)), and plasma density (panel (c)).</figcaption></figure>

Shock parameters and ambient plasma parameters along the four field lines are shown in Figures 4 and 5. As shown in Figure 5(a), the temperature of the magnetic loops drops to chromospheric values at the footpoints and increases radially outward to coronal values. Figure 4(a) confirms that for most of its propagation through the loop, the shock-normal speed exceeds 1000 km s⁻¹. During part of its propagation along the loops, the shock exhibits a quasi-perpendicular geometry with values of $\theta_{Bn}$ exceeding 80° (Figure 4(c)), with very high
speeds projected along the field line (Figure 4(b)). The typical $M_{\text{fm}}$ exceeds 2.7, meaning that the shock is supercritical (Figure 4(d)). To compute the gas compression ratio $\gamma_g$ and the magnetic compression ratio $\gamma_B$, we follow the recipe presented in Livadiotis (2015). Livadiotis (2015) defines the thermal ratio

$$ \zeta_{\text{TV}} \equiv \frac{2k_B T}{m_0(V_{\text{sn}} - V_{\text{w}} \cdot \mathbf{n})^2}, \quad (3) $$

where $k_B$ is the Boltzmann constant, $T$ is the temperature, and $m_0$ is the rest mass of the particle (see Equation (4) in Livadiotis 2015). Livadiotis (2015) also defines the ratio of the perpendicular component of the flow kinetic energy over that of the magnetic energy,

$$ \zeta_{\text{VB}\perp} \equiv \frac{(V_{\text{sn}} - V_{\text{w}} \cdot \mathbf{n})^2}{V_A^2 \cos^2 \theta_{\text{Bn}}}, \quad (4) $$

$$ \zeta_{\text{VB}} \equiv \frac{(V_{\text{sn}} - V_{\text{w}} \cdot \mathbf{n})^2}{V_A^2} \quad (5) $$

(see Equation 23(b) in Livadiotis 2015). The gas compression ratio $\gamma_g$ is given by the solution of the cubic trinomial,

$$ \Pi_3(\gamma_g) \equiv \sum_{k=0}^{3} q_k \gamma_g^k = 0 \quad (6) $$

(see Equation 24(a) in Livadiotis 2015), where

$$ q_0 = -(1 + d)\zeta_{\text{VB}}^2, \quad (7) $$

$$ q_1 = \left\{ \left(1 + \frac{3}{2}d\right) \cos^2 \theta_{\text{Bn}} + \left(1 + \frac{1}{2}d\right) + \left[1 + \left(1 + \frac{1}{2}d\right)\zeta_{\text{TV}}\right]\zeta_{\text{VB}} \right\}\zeta_{\text{VB}}, \quad (8) $$

$$ q_2 = -\left\{ \left[ (d+1) + \left(1 + \frac{1}{2}d\right)(1+2\zeta_{\text{TV}})\zeta_{\text{VB}} \right] \times \cos^2 \theta_{\text{Bn}} + \left(1 - \frac{1}{2}d\right)\zeta_{\text{VB}} \right\}, \quad (9) $$

$$ q_3 = \left[ 1 + \left(1 + \frac{1}{2}d\right) \zeta_{\text{TV}} \cos^2 \theta_{\text{Bn}} \right] \cos^2 \theta_{\text{Bn}}, \quad (10) $$

and $d$ indicates the degrees of freedom (see Equation 24(c) in Livadiotis 2015). For the adiabatic case here, $d=3$.

According to Press et al. (1992), the three solutions of $\Pi_3(\gamma_g)$ are given by

$$ \gamma_g = -\frac{1}{3q_3} \cdot \left( q_2 + hY + h^* \frac{\Delta_0}{Y} \right) \quad (11) $$

(see Equation 25(a) in Livadiotis 2015), where

$$ Y \equiv \sqrt[3]{\frac{\Delta_1 + \sqrt{\Delta_1^2 - 4\Delta_0^3}}{2}}, \quad (12) $$

$$ \Delta_0 \equiv q_2^2 - 3q_1q_3, \quad (13) $$

$$ \Delta_1 \equiv 2q_2^3 - 9q_1q_2q_3 + 27q_0q_3^2 \quad (14) $$

(see Equation 25(b) in Livadiotis 2015), and $h$ represents any of the (complex) cubic roots of unity

$$ h = 1, (-1 + i\sqrt{3})/2, (-1 - i\sqrt{3})/2, \quad (15) $$

$(h^*)$ is the conjugate of $h$; see Equation 25(c) in Livadiotis 2015). We limit solutions to fast-mode shocks. The magnetic compression ratio $\gamma_B$ is derived from Equation 22(c) in Livadiotis (2015),

$$ \gamma_B = \cos \theta_{\text{Bn}} \cdot \sqrt{1 + \gamma_g^2 (\xi + 1)^2 \cdot \tan^2 \theta_{\text{Bn}}}, \quad (16) $$

where

$$ \xi = \frac{1 - \gamma_g^{-1}}{\gamma_g^{-1} \zeta_{\text{VB}\perp} - 1} \quad (17) $$

(see Equation 23(a) in Livadiotis 2015).

Both the gas (Figure 4(e)) and magnetic (Figure 4(f)) compression ratios often exceed 2. Such a fast supercritical and quasi-perpendicular shock provides favorable conditions for SDA to occur along these loops.

Because these parameters are spline interpolated at a cadence of $\delta t = 0.5$ s, the time difference between neighboring data points is 0.5 s in Figures 4 and 5. At the footpoints of loops, the magnetic fields are almost radial, which results in large shock angles and large $V_{\text{sp}}$. Thus, at the footpoints of loops, the difference in heliocentric distances between neighboring data points is large. Therefore, in Figure 5, the plasma temperature shows an abrupt decrease at the footpoints, while the plasma density shows an abrupt increase at the footpoints. Consequently, the fast magnetosonic number $M_{\text{fm}}$, the gas compression ratio $\gamma_g$, and the magnetic compression ratio $\gamma_B$ increase abruptly at the footpoints.

From the four field lines shown in Figure 3, the extended loops FL 1487 and FL 2305 were well suited to study the bouncing of particles between converging magnetic mirrors imposed by the moving shock front and the strong magnetic fields at the footpoint of the loop, an effect our new model was designed to address. FL 1487 was actually connected at one footpoint very close to the $\gamma$-ray source and we have started simulating protons and electrons on FL 1487. This ongoing work will be presented in a future paper. The present paper is focused on the hard X-ray source, and FL 2305 was particularly interesting for the hard X-ray problem. The time interval over which the shock intersects this magnetic field line is shaded in yellow in Figure 1. The shock propagated through FL 2305 after the hard X-rays derived from soft X-ray measurements of the on-disk flare viewed from MESSENGER had ceased. There appears to be no softening of the GBM hard X-ray spectra during that phase, and our hypothesis was that hard X-rays could be further produced by the shock during a late connection to loops such as FL 2305. In future work, we plan to produce a full synthetic time-dependent spectrum of the hard X-ray emission integrating the contribution of all field lines.

### 3. SDA Computed with a Monte Carlo Approach

We have developed a Monte Carlo code that is capable of simulating SDA. We compute the change in pitch angle and speed of an electron reflected at the moving shock in the HT frame, in which the convective electric field disappears. Electron velocities parallel and perpendicular to the magnetic field before the reflection are $v_{i,\parallel}$ and $v_{i,\perp}$. Dividing them by the speed of light $c$, we obtain $\beta_{i,\parallel} = v_{i,\parallel}/c$ and $\beta_{i,\perp} = v_{i,\perp}/c$. Electron velocities parallel and perpendicular to the magnetic field after the reflection are $v_{r,\parallel}$ and $v_{r,\perp}$. Similarly, $\beta_{r,\parallel} = v_{r,\parallel}/c$
and $\beta_{r,\perp} = v_{r,\perp}/c$. $\beta_{i,\parallel}$ and $\beta_{i,\perp}$ are related to $\beta_{r,\parallel}$ and $\beta_{r,\perp}$ by

$$ \beta_{i,\parallel} = \frac{2\beta_s - \beta_{r,\parallel}(1 + \beta_s^2)}{1 - 2\beta_{r,\parallel}\beta_s + \beta_s^2} \quad (18) $$

(see Equation (1) in Mann et al. 2009), and,

$$ \beta_{i,\perp} = \frac{(1 - \beta_s^2)}{1 - 2\beta_{r,\parallel}\beta_s + \beta_s^2} \cdot \beta_{r,\perp}, \quad (19) $$

where $\beta_s = V_{sp}/c$ (see Equation (2) in Mann et al. 2009). Here $\beta_{i,\parallel}$, $\beta_{i,\perp}$, $\beta_{r,\parallel}$, and $\beta_{r,\perp}$ are all in the upstream plasma frame.

We assume $\alpha$ is the electron pitch angle in the HT frame, and $B_1$ and $B_2$ are the upstream and downstream magnetic field strengths. The electron will be transmitted to downstream if

$$ \sin^2 \alpha \le \sin^2 \alpha_c = B_1/B_2 = 1/\gamma_B. \quad (20) $$

Otherwise, it will be reflected at the shock. This criterion is equivalent to the reflection conditions

$$ \beta_{i,\parallel} \le \beta_s \quad (21) $$

(see Equation (3) in Mann et al. 2009) and

$$ \beta_{i,\perp} \ge \frac{\tan \alpha_c}{\sqrt{1 - \beta_s^2}} \cdot (\beta_s - \beta_{i,\parallel}) \quad (22) $$

(see Equation (4) in Mann et al. 2009). After the reflection, the electron gets a boost to its parallel speed, whereas its pitch angle decreases. Equations (21) and (22) provide the basic conditions for electron-shock interactions, which we implement in the Monte Carlo code.

In Appendix A, we test our Monte Carlo simulation of a single interaction of electrons with the shock by comparing our simulated distribution functions with those derived analytically (Mann et al. 2009).

## 4. Numerical Setup for a Magnetic Flux Tube

We now present the results of several numerical experiments illustrating the high variability and dynamic nature of particle acceleration that unfolds when shock waves propagate along a single loop. As already stated, we focus here on FL 2305. We run our Monte Carlo code along the magnetic loop by employing the plasma properties derived from the 3D MHD MAST model. The shock properties, the injected particle population, and the magnetic field lines along which the particles are propagating are all based on MAST. As illustrated in Figures 4 and 5, we can define the time evolution of the shock and the properties of the plasma immediately upstream of the shock from the MAST model.

A single magnetic flux tube is assumed to be surrounding FL 2305. Electrons are considered confined within this magnetic flux tube. The cross-sectional area of the magnetic flux tube A (s) scales as $B(s)^{-1}$, where B(s) is the magnetic flux density. Both A(s) and B(s) are functions of the position along FL 2305 s. The length of FL 2305 $s_0 \approx 2.19R_\odot$, where $R_\odot$ is the solar radius. The area of the solar surface in the shock's upstream covered by this magnetic flux tube is assumed to be 1 cm$^2$. The choice of 1 cm$^2$ is made to make the magnetic flux tube relatively thin. It has no impact on the total number of precipitating electrons over the entire visible disk. The cross-sectional area of the magnetic flux tube at FP2, A(s$_0$), is then formulated as $A(s_0) = \sin \psi \cdot 1 \text{ cm}^2$, where $\psi$ is the angle between the magnetic field line and the solar surface at FP2.

We introduce an injection time step $\Delta t_{inj}$. The number of physical electrons N(s) in a small spatial volume A(s)V$_{sp}$(s)Δt$_{inj}$ ahead of the shock is computed using plasma density n(s) from MAST:

$$ N(s) = n(s)A(s)V_{sp}(s)\Delta t_{inj}. \quad (23) $$

We assume 117 logarithmic bins in the direction of momentum with a minimum momentum $p_{min} = 1.7 \times 10^{-19} \text{ g cm s}^{-1}$, which corresponds to a minimum energy of $E_{min} = 0.01 \text{ keV}$. The maximum momentum $p_{max} = 5.4 \times 10^{-15} \text{ g cm s}^{-1}$ corresponds to a maximum energy $E_{max} = 1 \times 10^5 \text{ keV}$. Because the extremely high-energy electrons would not be sustained in the loops without a strong acceleration mechanism that can overcome the continuous precipitation, in all the following simulations in this paper, we assume that the initial electrons follow an e-folding spectrum built from a $\kappa$-distribution of $\kappa = 3.5$ with an exponential cutoff at $p_{cut} = 1.8 \times 10^{-17} \text{ g cm s}^{-1}$ that corresponds to the cutoff energy $E_{cut} = 100 \text{ keV}$. A value $\kappa = 3.5$ is reasonable in the corona (Pierrard et al. 1999). We have also verified that those electrons in the tail of the $\kappa$-distribution do not contribute much.

Physical electrons with momentum p between the grid points $p_{i-1/2}$ and $p_{i+1/2}$ are represented by pseudo-particles assigned to the interval $[p_{i-1/2}, p_{i+1/2}]$. In a small spatial volume A(s)V$_{sp}$(s)Δt$_{inj}$, and between the grid points $p_{i-1/2}$ and $p_{i+1/2}$, we inject N$_i$ pseudo-particles where $N_i = N_0 \cdot \exp(-p_i/p_{cut})$. $N_0 = 50,000$ for simulations where only a single interaction with the shock is allowed, while $N_0 = 2000$ for simulations that allow for multiple shock interactions per electron. The number of physical electrons between the grid points $p_{i-1/2}$ and $p_{i+1/2}$ is

$$ N_c(p_i) = \int_{p_{i-1/2}}^{p_{i+1/2}} \frac{dN(s)}{dp} dp. \quad (24) $$

We divide $N_c(p_i)$ by $N_i$ and obtain the weight $W_i$ of the injected pseudo-particle assigned to the interval $[p_{i-1/2}, p_{i+1/2}]$,

$$ W_i = \frac{N_c(p_i)}{N_i}. \quad (25) $$

In this paper, by “inject” we mean that we place $N_i$ pseudo-particles at the shock.

## 5. Single Interaction

We begin our numerical study by illustrating the efficiency of particle acceleration at different shock locations along a single magnetic loop (FL 2305) and do not model yet that a particle may encounter the shock multiple times. Thus, we compute the particle acceleration assuming that each pseudo-particle can interact at most once with the shock (one SDA cycle), and the shock is fixed at location I, II, III, ... or VI as shown in Figure 6. We have calculated the electron spectra for these six runs. The shock reaches locations I–VI at 11:19:19, 11:19:31, 11:20:45, 11:22:21, 11:24:29, and 11:24:42 UT, respectively. In Figure 1, we plot the times when FL 2305 intersects the shock as yellow shaded areas and denote the times at which the shock crosses locations I, II, III, ... VI on FL 2305 with vertical lines. Table 1 gives the shock properties at all these six shock locations. Our simulation here consists of two steps.

**Step 1:** The shock is fixed at location I, II, III, ... or VI shown in Figure 6. We inject pseudo-particles at each shock location
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Shock locations I (purple), II (red), III (green), IV (orange), V (magenta), and VI (brown) on the loop FL 2305 where the shock was fixed for a single interaction or forced to stop for multiple interactions (panel (a)). The east limb of the Sun as visible from Earth is indicated by a light green arc. The cyan arc shows the central meridian of the Sun. Magnetic field strength *B* vs. distance *s* along FL 2305 (panel (b)). Footpoint 1 (FP1, i.e., the downstream footpoint) and footpoint 2 (FP2, i.e., the upstream footpoint) are indicated by a black star and a gray star, respectively. Both FP1 and FP2 are on the visible disk.</figcaption></figure>

<table><caption>**Table 1**

Shock Properties at Locations I, II, III, IV, V, and VI</caption><thead><tr><th>Shock Location</th><th>Shock Normal Speed (km s<sup>-1</sup>)</th><th>Projected Shock Speed (km s<sup>-1</sup>)</th><th>Shock Normal Angle (°)</th><th>Fast-mode Mach Number</th><th>Gas Compression Ratio</th><th>Magnetic Compression Ratio</th></tr></thead><tbody><tr><td>I</td><td>1.02 × 10<sup>3</sup></td><td>5.48 × 10<sup>4</sup></td><td>88.9</td><td>3.77</td><td>3.44</td><td>3.44</td></tr><tr><td>II</td><td>1.10 × 10<sup>3</sup></td><td>9.83 × 10<sup>3</sup></td><td>83.6</td><td>4.08</td><td>3.56</td><td>3.54</td></tr><tr><td>III</td><td>1.25 × 10<sup>3</sup></td><td>2.38 × 10<sup>3</sup></td><td>58.4</td><td>6.05</td><td>3.88</td><td>3.34</td></tr><tr><td>IV</td><td>1.16 × 10<sup>3</sup></td><td>3.32 × 10<sup>3</sup></td><td>69.5</td><td>4.47</td><td>3.65</td><td>3.46</td></tr><tr><td>V</td><td>1.16 × 10<sup>3</sup></td><td>4.14 × 10<sup>3</sup></td><td>73.7</td><td>1.87</td><td>2.30</td><td>2.27</td></tr><tr><td>VI</td><td>1.16 × 10<sup>3</sup></td><td>4.81 × 10<sup>3</sup></td><td>76.0</td><td>2.48</td><td>2.77</td><td>2.73</td></tr></tbody></table>

with $N_0 = 50,000$. Then we perform shock interaction routines: if a pseudo-particle satisfies the reflection conditions Equations (21) and (22), it is reflected at the shock. This pseudo-particle is labeled as an accelerated pseudo-particle. If a pseudo-particle fulfills Equation (21), but it is inside the loss cone of the shock, it is absorbed, i.e., removed from the simulation. If a pseudo-particle does not satisfy Equation (21), it does not interact with the shock. This pseudo-particle is labeled as an escaping background pseudo-particle. Accelerated and escaping background pseudo-particles are collected. Based on the momenta, pitch angles, and weights of these pseudo-particles, we plot the logarithm of the physical electron number per momentum as a function of the logarithm of momentum using a base 10 logarithmic scale on the vertical and horizontal axes. Then we only place $10^{-2}$, $10^{-1}$, ...$10^5$ keV along the horizontal axis. The horizontal axis still represents the logarithm of momentum. Thus, $10^{-2}$, $10^{-1}$, ...$10^5$ keV are not evenly spaced along the horizontal axis. This is how we compute electron spectra in all the following figures as well. From 117 logarithmic bins in the direction of momentum $\Delta p_1$, $\Delta p_2$, $\Delta p_3$, ... $\Delta p_{117}$, we select a subset of 39 bins $\Delta p_3$, $\Delta p_6$, $\Delta p_9$, ... $\Delta p_{117}$. For all the spectra, we only plot the values for these 39 bins. The red curves in Figure 7 show the accelerated electron momentum spectra at each shock location. The initial

spectra are shown as black squares for comparison. Among all the initial electrons, those with a parallel speed that is higher than the projected speed of the shock along the magnetic field do not interact with the shock. We call them escaping background electrons.

Step 2: The shock is still fixed at location I, II, III, ... or VI. The simulation is performed in a fixed spatial simulation box in the upstream region of the shock. The box is bounded by the shock front and the visible disk at FP2 in Figure 6. All the accelerated and escaping background pseudo-particles collected in Step 1 are released at the shock. They then propagate from the shock toward the solar surface in the absence of shock interaction. The propagation of the pseudo-particle can be described by the guiding-center's motion along the magnetic field and pitch-angle focusing,

$$ \frac{ds}{dz} = \mu \qquad (26) $$

and

$$ \frac{d\mu}{dz} = -\frac{(1 - \mu^2)}{2B(s)} \cdot \frac{\partial B(s)}{\partial s}, \qquad (27) $$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Simulated electron spectra for which only a single interaction with the shock is allowed. The spectra of accelerated electrons are plotted as red curves. Black squares show the initial spectra. The spectra of precipitating electrons are presented as green squares. The dotted gray line in each panel is a power law with a spectral index of 5 in momentum space.</figcaption></figure>

where $\mu$ is the pitch-angle cosine and $z = vt$ with $v$ being the speed of the pseudo-particle. Equations (26) and (27) are solved with an Euler scheme. Pseudo-particles that cross the boundary at FP2 in Figure 6 are absorbed, i.e., precipitated. These pseudo-particles are labeled as precipitating pseudo-particles. Based on their momenta, pitch angles, and weights, we compute the spectra of precipitating electrons and plot them as green squares in Figure 7.

Because only a single shock encounter is allowed, the acceleration is overall limited, and the spectra of accelerated electrons (red curves) are very close to the initial spectra (black squares) in Figure 7. The accelerated electron velocity component parallel to the magnetic field has to be greater than the shock speed projected onto the field line. Thus, the accelerated electron velocity exceeds the shock speed projected onto the field line, which varies between different locations. Thus, the low-energy cutoff of the red curve in Figure 7 varies between the different spectra. At location I, the shock is nearly perpendicular ($\theta_{Bn} \sim 89^{\circ}$); therefore, the shock speed projected onto the magnetic field line is much higher (see Table 1). Thus, the shock at location I accelerates electrons to much higher energies than the shock at any other locations, with a peak occurring near 26 keV. At location II, $\theta_{Bn} \sim 84^{\circ}$, the shock becomes a little less efficient. The red curve almost overlaps the black squares at shock location III, which means the shock does

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Simulated spectra of precipitating electrons for which only a single interaction with the shock is allowed: shock location I (solid purple curve), shock location II (dashed red curve), shock location III (green dots), shock location IV (orange pentagons), shock location V (magenta crosses), and shock location VI (brown diamonds).</figcaption></figure>

not really accelerate electrons much. This is due to a much smaller $\theta_{Bn}$, which is about 58°. There is a little more acceleration at shock location VI because $\theta_{Bn}$ increases from shock location III to shock location VI. Figure 8 shows the varying precipitation for different injected locations. The significant particle acceleration computed near shock location I occurs soon after the shock first intersects the loop near one footpoint on the visible disk (FP1 in Figure 6). Precipitation of the electrons at FP2 in Figure 6 is also very high because the two footpoints share similar magnetic field strengths and therefore particles experience little mirroring from FP2 during their propagation. One can see that the closer to the loop top the acceleration occurs, the fewer the electrons that can precipitate due to the large magnetic field ratio between FP2 and the top of the magnetic loop that imposes a stronger mirror force on the propagating particles.

## 6. Multiple Interactions

Some electrons interacting once with the shock can be mirrored back during their propagation to FP2 and be naturally brought back to the shock. This could force a second shock interaction and further acceleration of electrons, which was one of the mechanisms that lead to multiple interactions with the shock studied by Sandroos & Vainio (2006). Some electrons interacting once with the shock can also be caught by the shock when the shock turns more and more perpendicular. These electrons receive a boost to their parallel speed even without the mirror force, which was described as the HT resonance effect in Sandroos & Vainio (2006) as well. In this section, we allow electrons to interact with the shock multiple times to illustrate the dramatic increase in energy experienced by these electrons.

We now compute particle acceleration and transport by allowing multiple shock interactions per electron. We have calculated the electron spectra for three runs. Run 1: the shock travels from the first shock-field line intersection to location II. Run 2: the shock travels from the first shock-field line intersection to location III. Run 3: the shock travels from the first shock-field line intersection to location VI. In this section as well as in Section 7, for each run, our simulation consists of two steps.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** In the same format as Figure 7, we show particle energy spectra obtained upstream of the shock for simulations that allow for multiple shock interactions per particle but do not allow particles to be scattered downstream of the shock to return to the shock. The spectra shown are for shock location II (a), III (b), or VI (c) only because they exhibit the most interesting features. As in Figure 7, they include accelerated electrons (red curves), the escaping background electrons (blue squares), and the electrons that precipitate (green squares). The dotted gray line in each panel is a power law with a spectral index of 5 in momentum space.</figcaption></figure>

Step 1: The simulation is performed in a shrinking spatial simulation box in the upstream region of the shock. The box is bounded by the moving shock front propagating along the magnetic loop and the visible disk at FP2 in Figure 6. We inject $N_0 = 2000$ new pseudo-particles every 0.5 s, i.e., $\Delta t_{\text{inj}} = 0.5$ s, as a means of simulating the number of ambient electrons in the loop as the shock moves along. For each pseudo-particle, Equations (26) and (27) are solved with a fourth-order Runge-Kutta scheme (Runge 1895; Kutta 1901). The position of a pseudo-particle is compared with the shock position. If a pseudo-particle encounters the propagating shock, the time step is reduced, and shock interaction routines described in Section 5 are performed. In this section as well as Sections 7 and 8, pseudo-particles are allowed to interact with the shock multiple times. The length of FL 2305 $s_0 \approx 2.19R_⊙$. It takes 325 s for the shock to travel from the first shock-field line intersection to location VI. In Run 1, Run 2, or Run 3, we collect those pseudo-particles that cross location II, III, or VI over a 10 s time interval just before the shock reaches location II, III, or VI, respectively. Those pseudo-particles can be divided into accelerated pseudo-particles and escaping background pseudo-particles. Accelerated pseudo-particles are those that have interacted with the shock at least once. Escaping background pseudo-particles are those that have never interacted with the shock. The 10 s interval is considered just for getting better statistics. We compute the spectra of accelerated electrons and the spectra of escaping background electrons in a similar manner to Section 5. They are shown as red curves and blue squares in Figure 9, respectively. We end the simulation in Step 1 when the shock hits location II, III, or VI in Run 1, Run 2, or Run 3, respectively.

Step 2: In Run 1, Run 2, or Run 3, the shock has stopped at location II, III, or VI, respectively. The simulation is performed between the shock front and the visible disk at FP2 in Figure 6. All of the accelerated and escaping background pseudo-particles collected in Step 1 are released at their own final

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Figure 8, but for multiple interactions with the shock. We assume the shock only stops at locations II, III, or VI rather than all the six locations (I–VI). We show the spectra of precipitating electrons, calculated from multiple interactions with the shock without particle return from downstream.</figcaption></figure>

positions in Step 1. They then propagate from the shock toward the solar surface in the absence of shock interaction. We solve Equations (26) and (27) and compute the spectra of precipitating electrons in a similar manner to Section 5. The spectra of precipitating electrons are plotted as green squares in Figure 9.

Comparing the escaping background electron spectrum (blue squares) with the accelerated spectrum (red curve) at the three shock locations, we see that significant acceleration occurs up to 43 keV by the time the shock reaches location II and hardly any particles are accelerated by the time the shock reaches location III. In contrast, a huge boost in acceleration has occurred by the time the shock reaches location VI near FP2 in Figure 6.

Comparing Figures 9 and 7, we see that the inclusion in our simulations of multiple interactions produces a wider enhancement in the spectrum (up to 43 keV) than for the case for single interaction near shock location II. This is because electrons are efficiently accelerated by a quasi-perpendicular shock ($\theta_{\text{Bn}} \sim 84^\circ$) near this shock location. In contrast, near shock location III, $\theta_{\text{Bn}} \sim 58^\circ$ (see Table 1) is significantly smaller than that at shock location II and the electrons do not gain much energy per SDA cycle. The strongest increase in the hardening of the electron spectrum is observed near shock location VI, shown in Figure 9(c), where energies exceeding 10 MeV are produced. As the shock propagates toward FP2, many more electrons were both injected and mirrored multiple times. As the shock propagates toward location VI, it also becomes increasingly quasi-perpendicular. Thus, the projected speed of the shock increases as a function of time, and the shock catches the electrons it had accelerated a moment ago, even without the mirror force (Sandroos & Vainio 2006). Particles could gain energy throughout the propagation of the shock. In order to get better statistics, we consider a 10 s time interval before the shock reaches a given location to calculate the spectra.

Figure 10 presents a comparison of the number of precipitating electrons corresponding to shock locations II, III, and VI taken from Figure 9. As in the case for single interaction (Figure 8), shock location III provides the fewest precipitating electrons. Shock location VI provides the most precipitating electrons because for each interaction with the
shock, the pitch angle of electrons decreases. Particles enter the loss cone and precipitate at FP2 in Figure 6.

The trapping of electrons between the moving shock and FP2 forces significant particle energy gain for as long as the particles are trapped between the two magnetic barriers i.e., until they enter the loss cone of either FP2 or the shock itself. HT resonance is also efficient in our magnetic field configurations.

## 7. Particle Return From Downstream

Electrons transmitted downstream of the shock can scatter and return to the shock front. This is because the region behind the shock is highly turbulent and particles can easily bounce back to the upstream side of the shock. To simulate this return process from the downstream region, we follow the recipe presented in Battarbee et al. (2013). Allowing for return from the downstream region of the shock did not change the simulation results significantly, and for clarity purposes and completeness, we present our approach to treating the return from downstream and the associated results in Appendix B.

## 8. Diffusive-shock Acceleration

Magnetic fluctuations such as plasma waves and turbulence may exist all along coronal loops. These fluctuations may cause some level of particle scattering that will also force some particles to diffuse back to the shock for further acceleration, thus resulting in the so-called DSA. In addition, particle scattering can also force some particles to precipitate at FP2. These effects of diffusion are investigated in this section.

The level of plasma turbulence in the corona is unknown at both the large fluid scales and the small electron scales. Insights have been gained by studying the level of electron scattering in the vicinity of flaring active regions in Kontar et al. (2014) by modeling radiation in these regions. Recent RHESSI observations of the flares of 2002 July 23, 2003 November 2, 2011 February 24, and 2011 September 24 have put constraints on the mean free path associated with the particle scattering. However, no quantitative conclusions about the strength of pitch-angle scattering have been made for the 2014 September 1 event. The level of turbulence in the quieter solar corona is even less well understood. The frequency of pitch-angle scattering $\nu$ is therefore unknown.

In this section, in Run 1 and Run 2, we only perform Step 1. In Run 3, we perform both Step 1 and Step 2. Here, Step 1 is implemented in a similar manner to Step 1 in Section 7, except we now perform an isotropic small-angle scattering after we obtain the particle position $s_j$ and pitch-angle cosine $\mu_j$ by solving Equations (26) and (27) at time $t_j$ with a fourth-order Runge-Kutta scheme. To model the DSA process, pseudo-particles are scattered through isotropic small-angle scattering at frequency $\nu$, similar to the method used in Agueda et al. (2008). For each particle time step $\Delta t_j$, the angle $\theta_s$ between the scattering axis and the direction of propagation of the scattered particle and the azimuthal angle $\phi_s$ around the scattering axis are randomized according to Kocharov et al. (1998), where $\theta_s = \sqrt{-2\nu\Delta t_j \log(1 - S)}$ and $\phi_s = 2\pi S'$, and $S$ and $S'$ are uniformly distributed random numbers in the interval [0,1). The solution to Equation (27) at time $t_j$ is $\mu_j$. Thus, the new pitch-angle cosine $\mu'_j$ after scattering is

$$ \mu'_{j} = \mu_{j} \cos \theta_{s} + \sqrt{1 - \mu_{j}^{2}} \sin \theta_{s} \cos \phi_{s}, \quad (28) $$

given by Ellison et al. (1990).

Pitch-angle scattering is performed in the entire upstream region of the shock, which can scatter electrons back to the shock. Electrons transmitted downstream of the shock can scatter and return to the shock front according to the probability of return $P_{\text{ret}}$. This results in the DSA process.

In Run 3, after Step 1, we also perform Step 2 in the absence of shock interaction and scattering. Step 2 is implemented between location VI and FP2 in Run 3. The segment of FL 2305 between location VI and FP2 is very short. We ignore scattering on this very short segment in Step 2. We have to admit that this two-step process is not entirely physical but it serves as a first approximation to consider DSA. Kontar et al. (2014) found that, in coronal loops, the mean free path $\lambda \sim (10^8-10^9)$ cm for $\sim 30$ keV (e.g., the electron speed $\nu = 10^{10}$ cm s$^{-1}$) electrons. The scattering frequency $\nu = \nu/\lambda$. So, for $\sim 30$ keV electrons, the scattering frequency $\nu \sim 10-100$ Hz. We have tried several different frequencies in the 10–100 Hz range. In Figure 11, we plot the spectra of accelerated and escaping background electrons for $\nu = 50$ Hz as red curves and blue squares. We also plot the spectrum of precipitating electrons in Run 3. One can see that the inclusion of particle diffusion has given a significant boost in the acceleration of particles and precipitation rate. The spectra for all three shock locations are harder and now extend up to 100 MeV. Electrons are reflected upstream via SDA first, gaining the energy necessary to be injected into DSA. Then, these reflected electrons are scattered back toward the shock for multiple cycles of acceleration. Thus, pitch-angle scattering is keeping particles in the vicinity of the shock, as expected in DSA. Pitch-angle scattering is also forcing a higher rate of precipitation by shifting particles into the mirroring loss cone of FP2 in Figure 6 sooner.

Strong diffusion is more likely to occur in the downstream region than in the upstream region. This strong turbulence combined with the stronger coronal magnetic field strength compressed in the downstream region will likely lead to an increased precipitation rate at the footpoint in the shock’s downstream (FP1 in Figure 6). We will discuss this topic in future work.

## 9. Discussion

The main goals of this study were to illustrate the evolution of electron acceleration at expanding shock waves and to determine whether electrons can be efficiently accelerated to relativistic energies and produce the hard X-rays measured by Fermi on 2014 September 1 shown in Figure 1 (Plotnikov et al. 2017). The study employed realistic shock properties constrained by multi-point imaging (Rouillard et al. 2016; Plotnikov et al. 2017), numerical models of the solar atmosphere, and a new numerical model of particle acceleration.

As discussed in Plotnikov et al. (2017), coronal shock waves propagate during the first hour of expansion through a complex system of magnetic loops. Figure 2 presented different loops magnetically connecting the shock wave of 2014 September 1 to the visible solar disk. To illustrate the variability of shock and ambient plasma parameters along different field lines, four different loops that magnetically connected the shock to the visible disk during the hard X-ray event were considered (Figure 3) and plotted in Figures 4 and 5. These figures illustrate (1) the considerable time variability in shock properties sampled by a single magnetic loop that must be taken into
<figure><img src="image_13.png" /><figcaption>**Figure 11.** In the same format as Figure 9, we show the spectra of accelerated (red curves) and escaping background (blue squares) electrons calculated from DSA. For shock location VI (c), we also show the spectra of precipitating electrons (green squares). The spectra are computed for a particle-scattering frequency $v = 50$ Hz. The dotted gray line in each panel is a power law with a spectral index of 5 in momentum space.</figcaption></figure>

account in any calculation of particle acceleration, and (2) that shock properties can differ greatly along different coronal loops. In order to illustrate some important mechanisms at play, the present study focused on one of these loops only that was crossed by the shock late in the hard X-ray event but just after the peak in γ-ray emission (Plotnikov et al. 2017).

We developed and tested a new model of particle acceleration at collisionless shocks that can take this variability into account. The model was tested first against the well-established SDA theory (discussed in Appendix A). Combining the particle acceleration code with the shock modeling we determined for a single magnetic field line (FL 2305), the important processes at play during electron acceleration to high energy.

The combination of the loop topology and the changing location of the shock-loop connection means that the shock geometry can evolve very rapidly along a single field line. SDA was efficient when the shock first intersected the loop (see Figure 7). As the shock reached the loop top, it temporarily became more quasi-parallel where few electrons could be accelerated by a single SDA cycle. An additional factor is related to the pool of particles available for acceleration. When the shock intersects the loop in the dense plasma prevalent in the low corona, the number of electrons available for acceleration is necessarily higher than near the more tenuous loop top.

We proceeded by allowing for multiple encounters with the shock, which clearly increased the number of accelerated and precipitating electrons (demonstrated in Figure 9). The effect of the mirror force and the resulting precipitation of electrons (see Figure 10) were similar to the previous case. We then allowed electrons transmitted through the shock to return from downstream (Appendix B). The number of accelerated electrons did not increase significantly, but slightly more electrons were able to precipitate. Finally, we moved on to the acceleration of

electrons through DSA. As demonstrated in Figure 11, particle scattering can boost the acceleration of electrons by forcing additional interactions with the shock and the precipitation rate by shifting particles into the loss cone of FP2 in Figure 6.

As mentioned above, particle-scattering frequencies $v$ are not known in the corona and only estimates have been inferred from past modeling of solar flares. We allow the shock to propagate through the entire loop to location VI (which is very close to FP2) and sum accelerated and escaping background electrons during the last 10 s just before the shock wave reaches location VI at 11:24:42 UT. Thus, accelerated and escaping background electrons that reach shock location VI from 11:24:32 UT to 11:24:42 UT are summed. We then propagate them to the chromosphere and count the number of precipitating electrons $N_{pr}$.

To evaluate whether the simulated number of precipitating electrons could even roughly explain the hard X-rays measured by Fermi GBM during the 2014 September 1 event, we compare the number of simulated electrons that precipitate along FL 2305 with the number of electrons inferred from the hard X-ray spectra (Plotnikov et al. 2017). The hard X-rays measured by GBM are integrated over the entire visible disk. For a meaningful comparison, we must therefore consider how the shock connects magnetically to the visible disk via many different realistic flux tubes. In this study, we assume that all field lines connecting the shock to the visible disk are like FL 2305. Because the shock properties along this particular loop were highly favorable for particle acceleration, our estimate of the number of energetic electrons precipitating at the surface of the Sun is likely to be an upper limit as most other loops will be less efficient.

The area of the solar surface in the shock's upstream covered by the magnetic flux tube surrounding FL 2305 is assumed to be $1 \text{ cm}^2$. According to the bottom panel of Figure 11 in Plotnikov et al. (2017), about 2% of the total visible solar
<table><caption>**Table 2**
A Comparison of Observationally Derived and Simulated Numbers of >200 keV Electrons and Spectral Indices between 11:24:32 UT and 11:24:42 UT</caption><thead><tr><td>Acceleration Process</td><td>Number of Electrons</td><td>Spectral Index</td></tr></thead><tbody><tr><td>SDA (without particle return from downstream)</td><td>6.6 × 10<sup>33</sup></td><td>7.8</td></tr><tr><td>SDA (with particle return from downstream)</td><td>2.9 × 10<sup>34</sup></td><td>7.7</td></tr><tr><td>SDA+DSA (ν = 30 Hz)</td><td>2.6 × 10<sup>34</sup></td><td>3.2</td></tr><tr><td>SDA+DSA (ν = 50 Hz)</td><td>1.1 × 10<sup>34</sup></td><td>1.4</td></tr><tr><td>Inferred from GBM</td><td>1.1 × 10<sup>32</sup></td><td>3.6</td></tr></tbody></table>

Note. This corresponds to a shock passage through location VI along FL 2305.

surface is magnetically connected to the shock, which is $2\% \times 2 \times \pi \times R_{\odot}^2 \approx 6 \times 10^{20} \text{cm}^2$. Assuming that all loops connecting to the visible disk are like FL 2305, there would be about $N_{pr} \times 6 \times 10^{20}$ precipitating electrons from 11:24:32 UT to 11:24:42 UT.

We fit Fermi GBM spectra using a thick-target electron-proton bremsstrahlung model and obtain $1.13 \times 10^{32}$ electrons above 200 keV from 11:24:32 UT to 11:24:42 UT. We report the simulated numbers of precipitating electrons and compare them to the number of electrons inferred from GBM in Table 2.

Simulations that take multiple interactions by SDA into account provide more than enough particles to explain the number of electrons inferred from GBM. For DSA simulations, shown in Figure 11, we assumed a particle-scattering frequency of $\nu = 50$ Hz. For $\nu = 50$ Hz, we found that $1 \times 10^{34}$ electrons precipitated with energies greater than 200 keV from 11:24:32 UT to 11:24:42 UT. We ran a similar simulation assuming a lower frequency of $\nu = 30$ Hz for which the number of precipitating electrons is reported in Table 2 as well. A value of $\nu = 30$ Hz is also reasonable in the loop (Kontar et al. 2014).

We remind the reader that the comparison between observed and modeled precipitating electrons is derived from only one loop with very favorable conditions, and this comparison does not consider the variability induced by different loops. In a future study, the numbers may change when we consider all loops, but this first simple comparison is very promising.

Estimates of the number of electrons from our simulations are interesting but prone to considerable uncertainty because we have not accounted yet for the variability of all field lines. A unique feature of GBM hard X-ray data is the relatively constant and hard power-law spectrum throughout the event shown in Figure 1(c). Our single-loop model provides a good ground for comparison with such a spectrum. The right-hand column of Table 2 presents the spectral indices derived from the different simulated spectra. The spectral indices derived from SDA simulations are greater than 7 even when multiple encounters are allowed and return from downstream is taken into account. Simulations that fold in diffusion (DSA) on the other hand are able to produce very hard spectra. In particular, for a particle-scattering frequency of $\nu = 30$ Hz, the spectral index (3.2) is very close to the observed one (3.6).

The analysis made the most of the available coronal observations. We have obtained realistic shock reconstructions at a high cadence that address to some extent the variability of acceleration processes along a magnetic loop. We have used realistic ambient coronal conditions derived from established benchmarked MHD simulations. There are nevertheless a number of limitations in the approaches taken in this work. The pitch-angle scattering frequency $\nu$ is unknown and will remain for some time an open question, and as already mentioned, we have performed our simulations on only one magnetic loop. A

future study must extend the analysis to calculations of particle acceleration and transport along multiple magnetic field lines. We have also ignored cross-field diffusion and electron transport between adjacent loops that might be possible in response to either a dynamic reconfiguration of coronal loops or a significant level of particle scattering.

## 10. Conclusions

Our model roughly reproduces the very hard electron spectrum and the number of electrons inferred from Fermi GBM. Thus, we have obtained more evidence for the CME-driven shock being an important accelerator of energetic electrons in the corona to produce hard X-rays far away from the flare. In this study, we have not addressed acceleration on open magnetic field lines that are connected to interplanetary space and potentially to probes situated in the inner heliosphere. Some electrons may be accelerated by the same shock along open magnetic fields and produce SEPs measured in situ. For the particular event studied here, Plotnikov et al. (2017) reported the occurrence of a strong SEP event with 2.8–4.0 MeV electrons measured at STEREO-B. Accelerating electrons to 2.8–4.0 MeV along open magnetic field lines will require a single open field line intersecting with the shock front multiple times (Bale et al. 2019; Kasper et al. 2019). We have found such cases in our shock modeling for this event. This mechanism will also be quantified in a future study. Alternatively, magnetic fluctuations may trap electrons in the vicinity of the shock, which leads to efficient DSA acceleration. In the presence of magnetic disturbances, electrons can reach 2.8–4.0 MeV along open magnetic field lines. Or potentially, electrons are accelerated along closed magnetic loops first. Then, these closed loops reconnect to open field lines, which allows those accelerated energetic electrons to escape to interplanetary space. Another possibility is that electrons first gain energy along closed loops, then reach open field lines via cross-field diffusion jump. This will be studied further in a future study.

In coming years, Parker Solar Probe and Solar Orbiter will no doubt measure many energetic particle events closer to the Sun in times of elevated electromagnetic radiation (which also will be measured with the STIX instrument on Solar Orbiter). This study is excellent preparation for future studies that will compare the populations of energetic particles released from the solar corona with those that impact the solar surface.

We acknowledge Gerald H. Share (Department of Astronomy, University of Maryland, College Park, MD 20742, USA) for providing the number and the power-law spectral index of electrons inferred from the Fermi/GBM observations. We also thank the referee for several helpful comments. The IRAP team acknowledges support from the French space agency (Centre
National des Etudes Spatiales, CNES; https://cnes.fr/fr). Y.W. and A.K. acknowledge financial support from the ANR project SLOW_SOURCE (ANR-18-ERC1-0006-01), COROSHOCK (ANR-17-CE31-0006-01), and FP7 HELCATS project, https://www.helcats-fp7.eu/, under the FP7 EU contract number 606692. The work of A.P.R. was funded by the ERC SLOW_SOURCE project (SLOW_SOURCE—DLV-819189). The work at the University of Turku was conducted in the framework of the Finnish Centre of Excellence in Research of Sustainable Space, funded by the Academy of Finland (grant 312357). R.V. and A.A. also acknowledge the financial support of the Academy of Finland of the project 309939. The work of A.W. was supported by DLR grant No. 50 QL 1701. Support by ISSI and ISSI-BJ through the international teams 469 is also acknowledged.

## Appendix A
### Test of the Monte Carlo Simulation

In this section, we assume the initial distribution function is an isotropic $\kappa$-distribution function of momentum $p$ given by

$$f_0(p) = \frac{1}{(\pi \kappa p_{\kappa}^2)^{3/2}} \cdot \frac{\Gamma(\kappa + 1)}{\Gamma(\kappa - 1/2)} \cdot \left[ 1 + \frac{p^2}{\kappa p_{\kappa}^2} \right]^{-\kappa - 1}, \quad (A1)$$

where $p_{\kappa} = m_0 w_{\kappa} / \sqrt{1 - w_{\kappa}^2/c^2}$ is the momentum corresponding to the thermal velocity $w_{\kappa}$ (see Equation (17) in Mann et al. 2009), and $w_{\kappa} = \sqrt{(2\kappa - 3)k_B T / \kappa m_0}$ with $k_B$ being the

Boltzmann constant, $T$ being the temperature, and $m_0$ being the rest mass of the particle. Equation (A1) can be transformed to a velocity distribution function (VDF) as

$$f_0(\beta_{i,\parallel}, \beta_{i,\perp}) = \left( \frac{\pi \kappa w_{\kappa}^2}{c^2 - w_{\kappa}^2} \right)^{-3/2} \cdot \frac{\Gamma(\kappa + 1)}{\Gamma(\kappa - 1/2)} \times \left\{ 1 + \frac{(c^2 - w_{\kappa}^2)(\beta_{i,\parallel}^2 + \beta_{i,\perp}^2)}{\kappa w_{\kappa}^2 [1 - (\beta_{i,\parallel}^2 + \beta_{i,\perp}^2)]} \right\}^{-\kappa - 1} \quad (A2)$$

The electron speed is randomized from the $\kappa$ distribution. In Figure 12(b), we present the simulated initial electron VDF with parameters: temperature $T = 1.5$ MK and $\kappa = 3$. We also assume $\theta_{Bn} = 89^{\circ}7$, $V_{sn} = 1000$ km s$^{-1}$, and $\gamma_B = 2$. If an initial electron satisfies the reflection conditions Equations (21) and (22), we reflect this electron. In Figure 12(d), we show the simulated VDF of these electrons that experience a single reflection.

In order to validate our Monte Carlo code, we solve analytically the VDF of the electrons that experience a single reflection. We consider the same initial electron VDF Equation (A2) with the same parameters: temperature $T = 1.5$ MK and $\kappa = 3$. In Figure 12(a), we present the analytical initial electron VDF. From Equations (21) and (22), we can derive that the analytical VDF of

<figure><img src="image_14.png" /><figcaption>Figure 12. Electron initial VDFs: analytical $f_0$ (panel (a)) and simulated initial VDF (panel (b)). Electron VDFs after one reflection: analytical $f_{acc}$ (panel (c)) and simulated VDF (panel (d)).</figcaption></figure>
the accelerated electron is a shifted loss-cone distribution (Mann et al. 2006, 2009),

$$f_{\text{acc}} = \Theta(\beta_s - \beta_{i,\parallel}) \times \Theta\left(\beta_{i,\perp} - \frac{\tan \alpha_c}{\sqrt{1 - \beta_s^2}} \cdot (\beta_s - \beta_{i,\parallel})\right) \times f_0(\beta_{i,\parallel}, \beta_{i,\perp}), \quad (A3)$$

where $\Theta$ is the Heaviside step function and $f_0$ denotes the initial VDF (see Equation (5) in Mann et al. 2009). Thus, the resulting analytical VDF of the electrons that experience a single reflection is

$$f_{\text{acc}} = \left( \frac{\pi \kappa w_{\kappa}^2 c^2}{c^2 - w_{\kappa}^2} \right)^{-3/2} \cdot \frac{\Gamma(\kappa + 1)}{\Gamma(\kappa - 1/2)} \\ \times \Theta\left( \beta_s - \frac{2\beta_s - \beta_{r,\parallel}(1 + \beta_s^2)}{1 - 2\beta_{r,\parallel}\beta_s + \beta_s^2} \right) \\ \times \Theta\left( \frac{(1 - \beta_s^2)\beta_{r,\perp}}{1 - 2\beta_{r,\parallel}\beta_s + \beta_s^2} - \frac{\tan \alpha_c}{\sqrt{1 - \beta_s^2}} \cdot \left[ \beta_s - \frac{2\beta_s - \beta_{r,\parallel}(1 + \beta_s^2)}{1 - 2\beta_{r,\parallel}\beta_s + \beta_s^2} \right] \right) \\ \times \left\{ 1 + \frac{c^2 - w_{\kappa}^2}{\kappa w_{\kappa}^2 (1 - \beta_s^2)^2 (1 - \beta_{r,\perp}^2 - \beta_{r,\parallel}^2)} \right. \\ \times [(1 - 2\beta_{r,\parallel}\beta_s + \beta_s^2)^2 - (1 - \beta_s^2)^2 (1 - \beta_{r,\perp}^2 - \beta_{r,\parallel}^2)] \}^{-\kappa-1}. \quad (A4)$$

In Figure 12(c) we show this analytic VDF of the electrons after one reflection calculated with the same parameters $\theta_{Bn} = 89.7^\circ$, $V_{sn} = 1000 \text{ km s}^{-1}$ and $\gamma_B = 2$. We compare the analytical and

simulated VDFs of the initial electrons and the electrons that experience a single reflection in Figure 12. The VDFs in the left column and the right column are identical, as expected.

## Appendix B
### Allowing for Return from the Downstream Region

The downstream plasma is assumed to be very turbulent so that the fate of the particle is basically determined instantly and based on local fluid properties. For isotropic downstream populations, the probability of return to the shock front was given in Jones & Ellison (1991) as

$$P_{\text{ret}} = \begin{cases} \left( \frac{v' - u_2}{v' + u_2} \right)^2, & v' > u_2 \\ 0, & v' \le u_2, \end{cases} \quad (B1)$$

where $v'$ is the particle velocity in the downstream rest frame and $u_2$ is the plasma speed in the downstream region in the HT frame. Successfully returning electrons re-enter the upstream region and may experience energy gains.

A similar simulation as presented in Section 6 is repeated here, except we now modify shock interaction routines to allow pseudo-particles transmitted through the shock to return from downstream. If a pseudo-particle fulfills Equation (21), but it is inside the loss cone of the shock, we calculate the return probability $P_{\text{ret}}$. If $P_{\text{ret}} > 0$, we return this pseudo-particle back to the upstream side of the shock and adjust its weight by multiplying its original weight with $P_{\text{ret}}$.

In Figure 13, we plot the spectra of accelerated, escaping background and precipitating electrons with the same scheme as before. It should be noted that returning from downstream is now taken into account.

When particles are allowed to return from downstream, slightly more particles precipitate from all shock locations compared with simulations when particles are not allowed to return (see Figure 14). This is due to high-energy particles that,

<figure><img src="image_15.png" /><figcaption>Figure 13. Same as Figure 9, but calculated from multiple interactions with the shock with particle return from downstream. The dotted gray line in each panel is a power law with a spectral index of 5 in momentum space.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** In the same format as Figure 10, we show the spectra of precipitating electrons, calculated from multiple interactions with the shock with particle return from downstream, in dashed red curve, green dots, and brown diamonds. For comparison, we show the spectra of precipitating electrons for multiple interactions with the shock without particle return from downstream in gray.</figcaption></figure>

in the previous set of simulations, were crossing the shock and removed from the simulation. Return from downstream being a diffusive process, these high-energy particles are sent back to the upstream region with the appropriate modifications to speed and pitch angle. If their new pitch angles are small enough, they can enter the loss cone and precipitate at FP2 in Figure 6. Figure 14 shows that the particles released from all three shock locations can now precipitate at high energies. These high-energy particles were in the previous simulations but were lost behind the downstream region and are now returned to the shock.

## ORCID iDs

Yihong Wu https://orcid.org/0000-0003-1660-0513

Alexis P. Rouillard https://orcid.org/0000-0003-4039-5767

Rami Vainio https://orcid.org/0000-0002-3298-2067

Alexandr N. Afanasiev https://orcid.org/0000-0001-9325-6758

Illya Plotnikov https://orcid.org/0000-0002-0074-4048

Alexander Warmuth https://orcid.org/0000-0003-1439-3610

## References

Ackermann, M., Allafort, A., Baldini, L., et al. 2017, ApJ, 835, 219

Afanasiev, A., Vainio, R., Rouillard, A. P., et al. 2018, A&A, 614, 4

Agueda, N., Vainio, R., Lario, D., et al. 2008, ApJ, 675, 1601

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Bale, S. D., Badman, S. T., Bonnell, J. W., et al. 2019, Natur, 576, 237

Ball, L., & Melrose, D. B. 2001, PASA, 18, 361

Battarbee, M., Vainio, R., Laitinen, T., et al. 2013, A&A, 558, 110

Bell, A. R. 1978a, MNRAS, 182, 147

Bell, A. R. 1978b, MNRAS, 182, 443

de Hoffmann, F., & Teller, E. 1950, PhRv, 80, 692

Decker, R. B. 1988, SSRv, 48, 195

Ellison, D. C., Jones, F. C., & Reynolds, S. P. 1990, ApJ, 360, 702

Erdos, G., & Balogh, A. 1994, ApJS, 90, 553

Gopalswamy, N., Mäkelä, P., Yashiro, S., et al. 2020, SoPh, 295, 18

Grechnev, V. V., Kiselev, V. I., Kashapova, L. K., et al. 2018, SoPh, 293, 133

Guo, X., Sironi, L., & Narayan, R. 2014, ApJ, 794, 153

Holman, G. D., & Pesses, M. E. 1983, ApJ, 267, 837

Hua, X.-M., Ramaty, R., & Lingenfelter, R. E. 1989, ApJ, 341, 561

Jin, M., Petrosian, V., Liu, W., et al. 2018, ApJ, 867, 122

Jokipii, J. R. 1982, ApJ, 255, 716

Jones, F. C., & Ellison, D. C. 1991, SSRv, 58, 259

Kasper, J. C., Bale, S. D., Belcher, J. W., et al. 2019, Natur, 576, 228

Klassen, A., Aurass, H., Klein, K. L., et al. 1999, A&A, 343, 287

Kocharov, L., Vainio, R., Kovaltsov, G. A., et al. 1998, SoPh, 182, 195

Kong, X., Guo, F., Giacalone, J., et al. 2017, ApJ, 851, 38

Kontar, E. P., Bian, N. H., Emslie, A. G., et al. 2014, ApJ, 780, 176

Kouloumvakos, A., Rouillard, A. P., Wu, Y., et al. 2019, ApJ, 876, 80

Krucker, S., Larson, D., Lin, R., et al. 1999, ApJ, 519, 864

Kutta, W. 1901, Z. Math. Phys., 46, 435

Lionello, R., Linker, J. A., & Mikić, Z. 2001, ApJ, 546, 542

Lionello, R., Linker, J. A., & Mikić, Z. 2009, ApJ, 690, 902

Livadiotis, G. 2015, ApJ, 809, 111

Mann, G. 1995, in Coronal Magnetic Energy Release, ed. A. O. Benz & A. Kruger (Berlin: Springer), 183

Mann, G., Aurass, H., & Warmuth, A. 2006, A&A, 454, 969

Mann, G., Warmuth, A., & Aurass, H. 2009, A&A, 494, 669

Marcowith, A., Bret, A., Bykov, A., et al. 2016, RPPh, 79, 046901

Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791

Murphy, R. J., Kozlovsky, B., Share, G. H., et al. 2007, ApJS, 168, 167

Petrosian, V. 2018, ApJ, 865, 99

Pierrard, V., Maksimovic, M., & Lemaire, J. 1999, JGR, 104, 17021

Plotnikov, I., Rouillard, A. P., & Share, G. H. 2017, A&A, 608, 43

Press, W. H., Teukolsky, S. A., Vetterling, W. T., et al. 1992, Numerical Recipes in FORTRAN. The Art of Scientific Computing (Cambridge: Cambridge Univ. Press)

Priest, E. 1982, Solar Magneto-hydrodynamics (Dordrecht: Reidel)

Rouillard, A. P., Plotnikov, I., Pinto, R. F., et al. 2016, ApJ, 833, 45

Rouillard, A. P., Sheeley, N. R., Tylka, A., et al. 2012, ApJ, 752, 44

Runge, C. 1895, MatAn, 46, 167

Salas, M. D. 2007, ShWav, 16, 477

Sandroos, A., & Vainio, R. 2006, A&A, 455, 685

Sandroos, A., & Vainio, R. 2009, A&A, 507, L21

Schlemm, C. E., Starr, R. D., Ho, G. C., et al. 2007, SSRv, 131, 393

Schmidt, J. M., & Cairns, I. H. 2012a, JGRA, 117, 4106

Schmidt, J. M., & Cairns, I. H. 2012b, JGRA, 117, 11104

Schwartz, S. J. 1998, ESA/ISSI, 1, 249

Share, G. H., Murphy, R. J., White, S. M., et al. 2018, ApJ, 869, 182

Vainio, R., Pönni, A., Battarbee, M., et al. 2014, JSWSC, 4, 8

Vršnak, B., Warmuth, A., Temmer, M., et al. 2006, A&A, 448, 739

Withbroe, G. 1988, ApJ, 325, 442

Wu, C. S. 1984, JGR, 89, 8857