On the origin of hard X-ray emissions from the behind-the-limb flare on 2014 September 1

The origin of hard X-rays and gamma-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 Ejections (CME) 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 ten minutes after the estimated on-disk hard X-rays appear to have ceased from the flare site.

) detected more than 30 solar eruptive events with late phase > 100 MeV γ-ray emission during its ten-year mission, among which the 2014 September 1 event was one of largest. The γ-ray emission was generated from 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 > 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. Since the magnetic footpoints of flaring loops associated with the 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 3-D 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 a 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 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. Since 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 powerlaw 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 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 since 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 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 ±1σ statistical uncertainties (Panel b), and > 200 keV electron (producing the hard X-ray emission) power-law spectral indices with ±1σ 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.
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 × 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 Xray 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, 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: firstly to present a new modeling framework to address particle acceleration at a shock propagating through coronal loops, and secondly 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 be isolated 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 Figure 1 b and 1 c with ±1σ statistical uncertainties. We do 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 km/s in quiet regions of the solar atmosphere with typical magnetosonic speeds of 300-400 km/s. 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 elections 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 combination of mechanisms could energise electrons to several hundred keVs or even tens of MeVs? (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 non-thermal emissions?

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 quasiperpendicular 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 October 28, 2003. It is also expected that SDA is responsible for the occurrence of Type II solar radio bursts.  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 participate 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, less shock encounters are required for particles to gain significant energy at quasi-perpendicular shocks. In this first theoretical study we focus on the relative role of SDA and DSA for the acceleration of electrons that propagate along magnetic loops traversed by the expanding shock wave.

Electron transport in the corona
To interpret 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 for 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,b;Jokipii 1982;Vainio et al. 2014;Afanasiev et al. 2018). In order to gain a more complete insight on the possible link between coronal shocks and the occurrence of highenergy 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) and 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 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.

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 3-D global MHD model of the coronal magnetic field. Then the variability of shock and ambient plasma parameters along different magnetic loops is illustrated.

Shock reconstruction and the background corona
The 3-D 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 Dynamic Observatory (SDO) spacecraft 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 higher cadence. The shock geometry is assumed to be an ellipsoidal structure. And the 3-D 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 3-D 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 3-D 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) onboard 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 3-D 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 3-D 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 ac- celeration and transport along such loops has not been studied in great detail, this is a point of focus in this article.

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 θ Bn are then determined at the shock along these 4 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 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 − δt, where δt varies between 0.5 second and 1.5 second. Next, we compute the distance between P and P , which we divide by the time interval δt. The speed of the propagating front along the field line is V sp = V sn / cos θ Bn . The fast-mode speed is then obtained by: where V A is the Alfvén speed and C S is the sound speed.
The fast magneto-sonic number M fm is defined as 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 δt = 0.5 second to produce a sequence of regularly time-spaced ellipsoids.
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 θ Bn exceeding 80 • (Figure 4 c) with very high speeds projected along the field line (Figure 4 b). The typical M fm exceeds 2.7 meaning that the shock is supercritical (Figure 4 d). To compute the gas compression ratio γ g and the magnetic compression ratio γ B , we follow the recipe presented in Livadiotis (2015). Livadiotis (2015) defines the thermal ratio where k B is the Boltzmann constant, T is the temperature and m 0 is the rest mass of the particle (cf Eq. (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,  (2015)). The gas compression ratio γ g is given by the solution of the cubic trinomial, q 1 = 1 + 3 2 d cos 2 θ Bn + 1 + 1 2 d and d indicates the degrees of freedom (cf Eq. (24c) in Livadiotis (2015)). For the adiabatic case here, d = 3.
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.
Since these parameters are spline interpolated at a cadence of δt = 0.5 second, the time difference between neighboring data points is 0.5 second 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 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. And consequently, the fast magnetosonic number M fm , the gas compression ratio γ g and the magnetic compression ratio γ 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 γ-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 to 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 late connection to loops such as FL 2305. In future work we plan to produce a full synthetic time-dependent spectrum of hard X-ray emission integrating the contribution of all field lines.

SDA COMPUTED WITH A MONTE CARLO APPROACH
We have developed a Monte Carlo code which 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, and v i,⊥ . Dividing them by the speed of light c, we ob-tain β i, = v i, /c and β i,⊥ = v i,⊥ /c. Electron velocities parallel and perpendicular to the magnetic field after the reflection are v r, and v r,⊥ . Similarly, β r, = v r, /c and β r,⊥ = v r,⊥ /c. β i, and β i,⊥ are related to β r, and β r,⊥ by Mann et al. (2009)) and where Mann et al. (2009)).
Here β i, , β i,⊥ , β r, and β r,⊥ are all in the upstream plasma frame.
We assume α 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 Otherwise, it will be reflected at the shock. This criterion is equivalent to the reflection conditions (cf Eq. (3) in Mann et al. (2009)) and (cf Eq. (4) in Mann et al. (2009)). After the reflection, the electron gets a boost to its parallel speed, whereas its pitch angle decreases. Eqs. (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).

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 3-D 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. 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 ψ · 1 cm 2 , where ψ is the angle between the magnetic field line and the solar surface at FP2.
We introduce an injection time step ∆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 We assume 117 logarithmic bins in the direction of momentum with a minimum momentum p min = 1.7×10 −19 g·cm/s which corresponds to a minimum energy of E min = 0.01 keV. The maximum momentum p max = 5.4×10 −15 g·cm/s corresponds to a maximum energy E max = 1×10 5 keV. Since 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 κ-distribution of κ = 3.5 with an exponential cutoff at p cut = 1.8×10 −17 g·cm/s that corresponds to the cutoff energy E cut = 100 keV. A value κ = 3.5 is reasonable in the corona (Pierrard et al. 1999). We have also verified that those electrons in the tail of the κ-distribution do not contribute much.
Physical electrons in the momentum bin p i are represented by pseudo-particles assigned to the momentum bin p i . In a small spatial volume A(s)V sp (s)∆t inj , and    for a momentum bin at p i , we inject N i pseudo-particles where N i = N 0 · exp(−p i /p cut ). N 0 = 50000 for simulations that 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 in the momentum bin p i is dN (s) dp dp.
We divide N e (p i ) by N i and obtain the weight W i of the injected pseudo-particle assigned to the momentum bin In this paper, by "inject" we mean that we place N i pseudo-particles at the shock.

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 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 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 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 with N 0 = 50000. Then we perform shock interaction routines: If a pseudo-particle satisfies the reflection conditions Eqs. (21) and (22), it is reflected at the shock. This pseudo-particle is labelled as an accelerated pseudo-particle. If a pseudo-particle fulfills Eq. (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 Eq. (21), it does not interact with the shock. This pseudo-particle is labelled as an escaping background pseudo-particle. Accelerated and escaping background pseudo-particles are collected. And based on the momenta, pitch-angles and weights of these pseudo-particles, we plot the logarithm of 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 ∆p 1 , ∆p 2 , ∆p 3 ... ∆p 117 , we select a subset of 39 bins ∆p 3 , ∆p 6 , ∆p 9 ... ∆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 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 guiding-center's motion along the mag-netic field and pitch-angle focusing: and where µ is the pitch-angle cosine and z = vt with v being the speed of the pseudo-particle. Eqs. (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 labelled 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.
Since 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 to the field line. Thus, the accelerated electron velocity exceeds the shock speed projected to 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 (θ Bn ∼ 89 • ), therefore, the shock speed projected to 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, θ Bn ∼ 84 • , the shock becomes a little less efficient. The red curve almost overlaps the black squares at shock location III, which means the shock does not really accelerate electrons much. This is due to a much smaller θ Bn , which is about 58 • . There is a little more acceleration at shock location VI because θ 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 since 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 electrons 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.

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 up with 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 Sections 7, for each run, our simulation consists of two steps.
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 second, i.e. ∆t inj = 0.5 second, as a means of simulating the number of ambient electrons in the loop as the shock moves along. For each pseudo-particle, Eqs. (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 & 8, pseudo-particles are allowed to interact with the shock multiple times. The length of FL 2305 s 0 ≈ 2.19R . It takes 325 seconds 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 ten-second 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 ten-second 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 the accelerated and escaping background pseudo-particles collected in Step 1 are released at their own final positions in Step 1. They then propagate from the shock toward the solar surface in the absence of shock interaction. We solve Eqs. (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 (θ Bn ∼ 84 • ) near this shock location. In contrast near shock location III, θ Bn ∼ 58 • (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 quasiperpendicular. 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 ten-second 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 is providing the least precipitating electrons. Shock location VI is providing 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.

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 treat the return from downstream and the associated results in Appendix B.

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 results in the so-called DSA. In addition, particle scattering can also 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 shown spectra are for shock locations 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 spectral index of 5 in momentum space. 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.
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 ν 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 smallangle scattering after we obtain the particle position s j and pitch-angle cosine µ j by solving Eqs. (26) and (27) 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 ν = 50 Hz. The dotted gray line in each panel is a power-law with spectral index of 5 in momentum space.
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 ν, similar to the method used in Agueda et al. (2008). For each particle time step ∆t j , the angle θ s between the scattering axis and the direction of propagation of the scattered particle and the azimuthal angle φ s around the scattering axis are randomized according to Kocharov et al. (1998), where θ s = −2ν∆t j log (1 − S) and φ s = 2πS , and S and S are uniformly distributed random numbers in the interval [0,1). The solution to Eq. (27) at time t j is µ j . Thus, the new pitch-angle cosine µ j after scattering is 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 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 λ ∼ (10 8 − 10 9 ) cm for ∼ 30 keV (e.g., the electron speed v = 10 10 cm s −1 ) electrons. The scattering frequency ν = v/λ. So for ∼ 30 keV electrons, the scattering frequency ν ∼ 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 ν = 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 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.

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 account in any calculation of particle acceleration, (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 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 (see Figures 13 and  14). 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 ν 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 seconds 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. Since 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 since 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 cm 2 . According to the bottom panel of Figure 11 in Plotnikov et al. (2017), about 2% of the total visible solar surface is magnetically connected to the shock, which is 2% × 2 × π × R 2 ≈ 6 × 10 20 cm 2 . Assuming that all loops connecting to the visible disk are like FL 2305, there would be about N pr × 6 × 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 × 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 ν = 50 Hz. For ν = 50 Hz, we found that 1 × 10 34 electrons precipitated with energies greater than 200 keV between 11:24:32 UT to 11:24:42 UT. We ran a similar simulation assuming a lower frequency of ν = 30 Hz for which the number of precipitating electrons is reported in Table 2 as well. A value ν = 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 since 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 1c. Our single loop model provides 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 ν = 30 Hz, the spectral index (3.2) is very close to the observed one (3.6).
The analysis made the most of available coronal observations. We have obtained realistic shock reconstructions at 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 ν 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.

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 Xrays 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 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 mag-netic loops first. Then these closed loops reconnect to open field lines, which allows those accelerated energetic electrons to escape to the 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 an 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.

A. TEST OF THE MONTE CARLO SIMULATION
In this section we assume the initial distribution function is an isotropic κ-distribution function of momentum p given by f 0 (p) = 1 (πκp 2 κ ) 3/2 · Γ(κ + 1) Γ(κ − 1/2) · 1 + p 2 κp 2 κ −κ−1 where p κ = m 0 w κ / 1 − w 2 κ /c 2 is the momentum corresponding to the thermal velocity w κ (cf Eq. (17) in Mann et al. (2009)). And w κ = (2κ − 3)k B T /κm 0 with k B being the Boltzmann constant, T being the temperature and m 0 being the rest mass of the particle. Eq. (A1) can be transformed to a velocity distribution function (VDF) as The electron speed is randomized from the κ-distribution. In Figure 12 b we present the simulated initial electron VDF with parameters: temperature T = 1.5 MK and κ = 3. We also assume θ Bn = 89.7 • , V sn = 1000 km/s and γ B = 2. If an initial electron satisfies the reflection conditions Eqs. (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 Eq. (A2) with the same parameters: temperature T = 1.5 MK and κ = 3. In Figure 12 a we present the analytical initial electron VDF. From Eqs. (21) and (22), we can derive the analytical VDF of the accelerated electron is a shifted loss-cone distribution (Mann et al. 2006(Mann et al. , 2009 where Θ is the Heaviside step function and f 0 denotes the initial VDF (cf Eq. (5) in Mann et al. (2009)). Thus, the resulting analytical VDF of the electrons that experience a single reflection is In Figure 12 c we show this analytic VDF of the electrons after one reflection calculated with the same parameters θ Bn = 89.7 • , V sn = 1000 km/s and γ 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.

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 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 Eq. (21), but it is inside the loss cone of the shock, we calculate the return probability P ret . If P 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 ret . 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 spectral index of 5 in momentum space.
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, 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.