Probing the Puzzle of Behind-the-limb γ-Ray Flares: Data-driven Simulations of Magnetic Connectivity and CME-driven Shock Evolution

To reconstruct the global corona and solar wind environment during the SOL2014-09-01 CME eruption, we used the University of Michigan Alfvén Wave Solar Model (AWSoM; Sokolov et al. 2013; van der Holst et al. 2014) within the Space Weather Modeling Framework (SWMF; Tóth et al. 2012). AWSoM is a data-driven global MHD model with the inner boundary specified by observed magnetic maps and the simulation domain extending from the upper chromosphere to the corona and heliosphere. Physical processes implemented in the model include multispecies thermodynamics, electron heat conduction (both collisional and collisionless formulations), optically thin radiative cooling, and Alfvén-wave turbulence that energizes the solar wind plasma. The Alfvén-wave description is physically self-consistent, including non-Wentzel–Kramers–Brillouin reflection (Heinemann & Olbert 1980; Velli 1993; Hollweg & Isenberg 2007) and physics-based apportioning of turbulence dissipative heating to both electrons and protons. AWSoM has demonstrated its capability of reproducing the solar corona condition with high fidelity (e.g., Oran et al. 2013, 2015; Sokolov et al. 2013; van der Holst et al. 2014; Jin et al. 2016, 2017a).

Based on the steady-state global corona and solar wind solution, we initiate the CME by using an analytical Gibson–Low (GL) flux-rope model (Gibson & Low 1998), which has been successfully used in numerous modeling studies of CMEs (e.g., Manchester et al. 2004a, 2004b, 2014; Lugaz et al. 2005; Jin et al. 2016, 2017a). The GL flux rope is mainly controlled by five parameters: the stretching parameter a determines its shape, the distance r₁ of the flux rope center from the center of the Sun determines its initial position before being stretched, the radius r₀ of the flux-rope torus determines its size, a₁ determines its magnetic field strength, and a helicity parameter determines its positive (dextral) or negative (sinistral) helicity. Analytical profiles of the GL flux rope are obtained by finding a solution to the magnetohydrostatic equation $({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}-{\rm{\nabla }}p-\rho {\boldsymbol{g}}=0$ and the solenoidal condition ${\rm{\nabla }}\cdot {\boldsymbol{B}}=0$. This solution is derived by applying a mathematical stretching transformation $r\to r-a$ to an axisymmetric, spherical ball of twisted magnetic flux with radius r₀ centered in the heliospheric coordinate system at $r={r}_{1}$. The transformed flux rope appears as a tear-drop shape. At the same time, Lorentz forces are introduced, which lead to a density-depleted cavity in the upper portion and a dense core at the lower portion of the flux rope, corresponding to a coronal cavity and a dense prominenence, respectively. This configuration can thus readily reproduce the typical three-part structure of an observed CME (Illing & Hundhausen 1985). The GL flux rope and contained plasma are then superposed onto the steady-state AWSoM solution of the solar corona: i.e., $\rho ={\rho }_{0}+{\rho }_{\mathrm{GL}}$, ${\boldsymbol{B}}={{\boldsymbol{B}}}_{0}+{{\boldsymbol{B}}}_{\mathrm{GL}}$, and $p={p}_{0}+{p}_{\mathrm{GL}}$. The temperature will be updated from the new density ρ and pressure p. The resulting combined background-flux rope system is in a state of force imbalance, due to the insufficient background plasma pressure to counter the magnetic pressure of the flux rope, and thus erupts immediately when the numerical model advances in time.

To specify the inner boundary condition of the magnetic field, we utilize a global magnetic map sampled from an evolving photospheric flux transport model (Schrijver & DeRosa 2003), which assimilates new observations within 60° from disk center obtained by the SDO Helioseismic and Magnetic Imager (Schou et al. 2012). The assimilated magnetogram is updated every 6 hr. The Sept14 flare occurred behind the east limb where no direct observation of the magnetic field is available. This means that the magnetic field around the flare site at the time of the event contains the most aged observation obtained from about a half solar rotation earlier when the region was on the western side of the visible solar disk. Therefore, a large amount of magnetic flux could potentially be missing. From the magnetogram closer in time shown in Figure 1(a), we find that the flare source region AR 12158 is indeed completely missing. To alleviate this problem, we choose the assimilated magnetogram on 2014 September 8 00:04:00 UT (Figure 1(b)), about a week after the event on September 1, when the magnetic field around the source region was first assimilated into the flux transport model. The missing flare source region AR 12158 and another large AR 12157 to the south of it are now properly included. The rest of the old and new magnetic maps are qualitatively very similar. As such, the 2014 September 8 magnetogram is a reasonable representation of the photospheric magnetic field at the time of the Sept14 flare and is thus used to specify the inner boundary condition of our global magnetic field model.

Zoom In Zoom Out Reset image size

To configure a proper GL flux rope for initiating the Sept14 CME, we utilize a newly developed tool called the Eruptive Event Generator using Gibson–Low configuration (EEGGL; Jin et al. 2017b), which is designed to determine the GL flux-rope, including its location, orientation, and five key controlling parameters, using the observed magnetogram and CME speed near the Sun. The left panel of Figure 2(a) shows a zoom-in view of AR 12158 with weighted centers of positive/negative polarities and the polarity inversion line determined by EEGGL. The green asterisk marks the central location to insert the GL flux rope, whose calculated key parameters are also listed. The right panel shows the 3D configuration of the global coronal magnetic field, with the inserted GL flux rope shown in red. The white field lines represent the large-scale helmet streamer structures. The selected field lines from surrounding ARs and open field are marked in green. The GL flux rope erupts due to the force imbalance upon insertion into the AR. The simulation is then evolved forward in time and the MHD equations are solved in conservative forms to guarantee the energy conservation across the CME-driven shock (van der Holst et al. 2010; Manchester et al. 2012; Jin et al. 2013). To better resolve the shock structure, two more levels of refinement along the CME path are performed, which make the cell size ∼0.02 R_⊙ at 2 R_⊙ and ∼0.06 R_⊙ at 5 R_⊙. We run the simulation for 1 hr after the initiation, until the CME reaches ∼10 R_⊙.

Zoom In Zoom Out Reset image size

Since the CME propagation near the Sun is mainly observed by coronagraphs, we generate synthesized white-light images (Thomson-scattered white-light brightness) and compare them with observations. The top panel of Figure 3 shows the observations from SOHO/LASCO C2 and STEREO-B/COR1. The bottom panel shows the synthesized white-light images. The color scale shows the relative total brightness changes with respect to the pre-event level. At the time of the Sept14 event, STEREO-B and SOHO were separated by ∼161°, therefore, observing the Sun from nearly opposite directions. By comparing the observation and simulation from two different view points, we find that the observed CME is adequately simulated in terms of the direction of propagation and the width. Note that the absolute brightness comparison between the observation and simulation requires advanced calibration of the observational data as well as the inclusion of the contribution from the F corona (light scattered by interplanetary dust) in the simulation data, which are beyond the scope of this study. The reader is referred to previous studies for such model validation (e.g., Manchester et al. 2008; Jin et al. 2017a). Also, there is a distinct feature (marked in Figure 3) around the CME leading edge in the LASCO C2 image, which might imply that the corresponding shock front can deviate from a typically circular or dome shape. We speculate that this feature might be related to the complex and dynamically changing background solar wind environment (e.g., affected by previous CMEs or coronal disturbances), which is not captured by the current simulation. We further compare the observed and simulated CME speeds by tracking the HT history of the CME leading edge, as shown in Figure 4. The black dots show measurements from SOHO/LASCO C2/3 (left panel) and STEREO-B COR1/2 (right panel), while the red asterisks show corresponding measurements from the synthesized white-light images. We use the moment when the observed and simulated CMEs are around the same height as a guidance to calibrate the start time of the simulation in terms of the real observation. With this assumption, the first appearance of CME in the LASCO C2 field of view (11:12:05 UT) corresponds to t = 10 minutes in the simulation. In general, the CME HT history is well reproduced in the simulation, with the simulated CME being slightly slower by about 10%–14%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the course of the eruption, the flux rope interacts and reconnects with the magnetic fields of the source AR as well as the global coronal field. As a result, the magnetic field configuration and connectivity can change dramatically, which could significantly influence the transport of the accelerated particles. With this global MHD simulation of the Sept14 event, we now investigate the field connectivity evolution in detail during the first hour of CME evolution.

Figures 5(a)–(d) show the 3D magnetic field configurations at selected times (5, 10, 20, and 30 minutes). Magnetic reconnection between the erupting flux rope (red) and the surrounding field lines (green and white) is evident, especially after the first 10 minutes. The interaction between the flux rope and the large-scale helmet streamers significantly changes the global corona configuration around the CME source region. The helmet streamers are opened up by reconnection or stretched by the CME expansion. Specifically, we further examine the field line connectivity around the Fermi γ-ray emission region at t = 30 minutes (shown in Figure 5(e)). The derived Fermi γ-ray emission centroid and 68% uncertainty circle (adapted from Ackermann et al. 2017) are overlaid on the simulation data. The green field lines are the pre-existing open field connected to the CME-driven shock after t ∼ 6 minutes. The red field lines are the closed field connected to the flaring AR. These field lines were not present before, but started to develop ∼5 minutes after the eruption through magnetic reconnection between the flux-rope magnetic field and the global coronal field.

Zoom In Zoom Out Reset image size

To investigate the details of these two types of field lines, we further mark their photospheric footpoints on the magnetic field map in the top panel of Figure 6. The closed field line regions (red) are relatively compact, compared with the elongated open field line region (green) to the south. The bottom panel of Figure 6 shows the 3D configuration of these two types of field lines. As is evident in this plot, the open field lines change directions abruptly due to the rapid expansion of the CME and CME-driven shock. For the closed field lines, the configuration is more complex with twisted large-scale loops. It appears that some of these field lines result from reconnection between the erupting flux rope and the nearby helmet streamers.

Zoom In Zoom Out Reset image size

After the eruption, the flux rope drives a shock in the corona that propagates freely into the heliosphere. CME-driven shocks are believed to be responsible for acceleration of particles through the diffusive shock acceleration (DSA) mechanism (e.g., Axford et al. 1977) that produces the so-called gradual SEP events (Reames 1999). Due to the nonuniform background environment, the CME-driven shock evolution is highly spatially dependent. For example, a shock that is propagating into the fast solar wind could acquire a higher shock speed therefore leading to a higher stand-off distance from the flux rope driver (Jin et al. 2017a). Also, the shock parameters could vary significantly over the shock front, which can significantly affect the acceleration process (Manchester et al. 2005; Li et al. 2012). Based on the white-light observations from SOHO and STEREO, several methods have been developed to derive the shock parameters directly from the data (e.g., Rouillard et al. 2016; Lario et al. 2017; Kwon & Vourlidas 2018). In this study, with the data-driven MHD simulation of the Sept14 event, we can track the shock location and key parameters (e.g., the compression ratio (CR), shock Alfvén Mach number, shock speed, and shock obliquity angle θ_Bn) during the CME evolution. The shock obliquity angle θ_Bn refers to the angle between the shock normal (see Equation (1)) and the upstream magnetic field. As shown below, such an analysis can provide a more comprehensive picture of the shock as to its configuration and properties over the area linking back to the visible side of the Sun, where LAT γ-rays were detected.

We first determine the shock location at each time step by using the proton temperature gradient criteria (Jin et al. 2013). At each shock location, the shock normal is determined by using the magnetic coplanarity $({{\boldsymbol{B}}}_{{\rm{d}}}-{{\boldsymbol{B}}}_{{\rm{u}}})\cdot {\boldsymbol{n}}=0$ (Lepping & Argentiero 1971; Abraham-Shrauner 1972):

$$\begin{eqnarray}&&{\boldsymbol{n}}=\pm \displaystyle \frac{({{\boldsymbol{B}}}_{{\rm{d}}}\times {{\boldsymbol{B}}}_{{\rm{u}}})\times ({{\boldsymbol{B}}}_{{\rm{d}}}-{{\boldsymbol{B}}}_{u})}{| ({{\boldsymbol{B}}}_{{\rm{d}}}\times {{\boldsymbol{B}}}_{{\rm{u}}})\times ({{\boldsymbol{B}}}_{{\rm{d}}}-{{\boldsymbol{B}}}_{{\rm{u}}})| },\end{eqnarray} \tag{ 1 }$$

where B_d and B_u represent downstream and upstream magnetic field respectively. Note that this method fails for θ_Bn = 0° or 90°, which we found to be very rare in the actual simulations. The ± sign is determined by assuming a forward moving shock in the heliocentric coordinate. We then determine the upstream/downstream plasma parameters, from which the shock parameters (e.g., the CR, shock speed, shock Alfvén Mach number, and shock obliquity angle) can be calculated accordingly.

Figure 7 shows the shock geometry evolution during the first hour of the simulation. The color scales on the shock surface represent the four shock parameters (from top to bottom): the CR, shock speed, shock Alfvén Mach number, and shock θ_Bn. The yellow field lines represent the open field near the Fermi emission region (shown in Figures 5 and 6). Based on the simulation, we found that the CME-driven shock started to intersect the open field lines around t = 6 minutes, when the fastest part of the shock reached ∼3 R_⊙. This finding is consistent with the estimation of the CME located at ∼2.5 R_⊙ at the onset of the Fermi-LAT emission (Ackermann et al. 2017). However, we should note that the part of the shock intersecting the open field is closer to the Sun at ∼1.6 R_⊙. At t = 20 minutes, the CME-driven shock covered the entire open field region around ∼3 R_⊙ linking to the front side of the Sun. Furthermore, the derived θ_Bn suggests that this part of the shock is a quasi-perpendicular shock with a mean θ_Bn ∼ 73°. Another observational fact worth mentioning is the EUV wave observed in this event. The EUV wave from the source region arrived at the open field region (connecting to the CME-driven shock) by 11:20 UT as shown in online movies (http://aia.lmsal.com/AIA_Waves), an extension of the Nitta et al. (2013) study. There is a possibility that this EUV wave may trace the low-corona flank of the shock, as see in several other eruptive events (e.g., Carley et al. 2013), including the recent X8.2 flare on 2017 September 10 (Gopalswamy et al. 2018a; Liu et al. 2018; Morosan et al. 2018), which was the first (and the only one so far) long-duration Fermi flare associated with a ground level enhancement event (Omodei et al. 2018).

Zoom In Zoom Out Reset image size

At t = 30 minutes, the shock surface starts to deviate from its initial spherical shape due to the nonuniform background solar wind condition. In particular, one part of the shock (as marked by a white arrow in Figure 7), with open field lines crossing it, propagates into a fast wind region originating from an on-disk coronal hole and acquires a higher speed. This process may also lead to a "shock–shock" interaction at the boundary between the two shock surfaces that causes an elevated shock CR and Alfvén Mach number at t = 60 minutes (marked with a white circle in the upper right panel of Figure 7). Note that the open field lines connected to this shock interaction region are closer to the Fermi γ-ray emission region. Since the CR, shock Alfvén Mach number, and shock geometry are key parameters for DSA of SEPs, this part of the shock can be favorable for accelerating particles to higher energies (Manchester et al. 2005; Sokolov et al. 2006; Li et al. 2012; Zhao & Li 2014; Hu et al. 2017, 2018).

We obtain the shock parameters averaged over the portion of the shock surface that is connected back to the visible side of the Sun. The temporal evolution of the resulting four key shock parameters as well as the average upstream local plasma density is shown in Figure 8 and described as follows:

• 1.  
  The shock CR (Figure 8(a)) increases rapidly from ∼1.8 at t ∼ 10 minutes to ∼4.6 at t ∼ 20 minutes,¹¹ and then gradually decreases to ∼3.7 at t = 60 minutes. This evolution trend closely follows that of the Fermi/LAT γ-ray flux profile (red curve), with a similar ∼10 minute duration of the rapid rise phase.
• 2.  
  The average local plasma density at the shock front (Figure 8(b)) is another important parameter, which is related to the seed population and the number of particles available for shock acceleration. We also plot an empirical quantity CR · ρ^1/3 (blue curve), a product of the shock CR and the ambient density to a 1/3 power (heuristically selected to match the temporal trend of the Fermi γ-ray flux). The density generally decreases with time as the CME travels away from the Sun, which causes this empirical quantity to decrease after its initial increase during the first  20 minutes. This could potentially explain the simultaneous decrease in the Fermi γ-ray flux, even though the shock CR and Mach number remain high.
• 3.  
  The shock speed (Figure 8(c)) shows a gradual increase from ∼400 km s⁻¹ at t ∼ 10 minutes to ∼1000 km s⁻¹ at t ∼ 35 minutes and then remains roughly constant.
• 4.  
  Likewise, the shock Alfvén Mach number (Figure 8(d)) gradually increases from ∼1 to ∼3 during the t ∼ 10–35 minute interval and then grows even more slowly to ∼4 at t ∼ 60 minutes.
• 5.  
  The shock obliquity angle analysis (Figure 8(e)) shows that the shock is originally a quasi-perpendicular shock before t ∼ 30 minutes (with θ_Bn ∼ 75° at t ∼ 10 minutes) and evolves into a quasi-parallel shock (with θ_Bn ∼ 30° at t = 60 minutes). The most rapid decrease in θ_Bn occurs during t ∼ 22–35 minutes, the onset of which coincides with the peak time (t = 22 minutes) of the CR and Fermi γ-ray flux shown in Figure 8(a). Note that the same trend of shock obliquity angle variation during CME evolution was also found by Manchester et al. (2005).

Zoom In Zoom Out Reset image size

As noted earlier in Section 3.2, the Fermi γ-ray emission centroid is closer to the footpoints of the closed field lines connecting the source AR than those of the open field lines connecting the CME shock. It is possible that the particles accelerated at the flaring site or the areas passed by the coronal shock (Hudson 2018) can be trapped in the closed field through some mechanisms (e.g., Sheeley et al. 2004), reaccelerated, and then transported to the front side of the Sun through the connectivity established by the interaction between the erupting flux rope and global corona field. However, we would like to emphasize that, based on the present simulation result, it is difficult to unambiguously distinguish the potential contributions from the two groups of field lines to the observed LAT emission for four reasons. (i) The northern portion of the open field footpoints is adjacent in space to both the closed field footpoints and the LAT centroid. (ii) The appearances of the closed field (t ∼ 5 minutes) and open field (t ∼ 6 minutes) are very close in time. (iii) The current localization of the LAT centroid is based on an assumption of a point source for the γ-ray emission and can change if the actual source shape deviates from this assumption. (iv) Because of the proximity between the open and closed field lines, cross-field diffusion (e.g., Zhang et al. 2003) can allow shock-accelerated energetic particles to access the closed field lines as well.

However, our simulation result of the Sept14 event together with our initial inspection of other Fermi BTL flares, does reveal some unique and attractive features about the CME-driven shock linked to the observed γ-rays:

• 1.  
  Since the shock CR is one of the key parameters in the DSA mechanism that determines the energetic particle production at the shock (e.g., the particle spectral index), the temporal correlation noted in Item 1 (Section 3.3) above indicates an intimate relation between the γ-ray flux and the shock particle production. This provides clear evidence supporting the mechanism that (at least some of) the γ-ray producing particles are accelerated by the CME-driven shock.
• 2.  
  The open field is connected to a quasi-perpendicular shock early on (Item 4 in Section 3.3), which is generally believed to be an efficient particle accelerator if the upstream coronal or heliospheric magnetic field is sufficiently turbulent (e.g., Giacalone 2005; Tylka et al. 2005). Furthermore, a recent in-situ observation of Saturn's bow shock from the Cassini spacecraft (Masters et al. 2017) shows that energetic electrons were only detected downstream of the quasi-perpendicular shock, which suggests the potential importance of a quasi-perpendicular shock in accelerating particles that could escape the downstream and propagate back to the Sun to produce γ-rays.
• 3.  
  Another piece of supporting evidence for shock acceleration is that in all three identified Fermi BTL events (Sept14, SOL2013-10-11, and SOL2014-01-06; Pesce-Rollins et al. 2015; Ackermann et al. 2017), there are pre-existing open field lines (e.g., in on-disk coronal holes) near the γ-ray emission region, which could be potentially connected to the CME-driven shock.

By using a 3D triangulation technique, Plotnikov et al. (2017) reconstructed the CME-driven shock structure from white-light observations of the Sept14 event. With the density and magnetic field information obtained from a static solar corona constructed with the Magnetohydrodynamic Algorithm outside a Sphere (MAS) model (Lionello et al. 2009), the time-dependent distribution of the shock Mach number and obliquity angle were approximately derived. They found that the Mach number shows a rapid increase to supercritical values after the type-II burst onset and the shock has a quasi-perpendicular geometry during the γ-ray emission, which are in general agreement with our results. However, an important distinction between their and our studies is that, instead of using a static coronal model, we self-consistently simulated the dynamic evolution of the CME and the CME-driven shock. This allows us to track the detailed temporal evolution of the shock and derive the shock CR, which is of critical importance to particle acceleration by shocks. In addition, we found that the shock geometry evolves and changes from quasi-perpendicular to quasi-parallel, instead of remaining quasi-perpendicular all the time.

We also briefly discuss the shock Alfvén Mach number evolution derived from the simulation. Note that when this number is around unity (see Item 3 in Section 3.3), stochastic acceleration of particles by plasma turbulence (e.g., in the downstream of the shock) is more efficient than DSA (Petrosian 2016). Stochastically accelerated particles could also serve as the seed population to be further accelerated by the shock. This could be the case early on, when the shock CR is also relatively low, and could be related to the rapid rise in the detected LAT γ-ray flux. Later on, when the shock Alfvén Mach number and CR are sufficiently large, shock acceleration would be more important and can account for the gradual, long-duration γ-ray emission. Therefore, it is likely that, in addition to DSA, stochastic acceleration could also play a role in the Sept14 event, especially in the early stage. On the other hand, as mentioned in Section 3.1, the simulated CME speed is 10%–14% slower than the observed one. Considering the higher shock speed in the simulation, the Alfvén Mach number will also be higher. The critical Alfvén Mach number is estimated between 1 and 1.7 for the solar wind plasma with $\gamma =\tfrac{5}{3}$ and β = 1.0 (Edmiston & Kennel 1984). Therefore, it is also possible that the shock in the simulation already becomes supercritical in the early stage. In summary, we believe multiple acceleration mechanisms could be important in the beginning but DSA might be the dominant one after 20 minutes.

Finally, we discuss the possibility of particle transport back to the Sun from the CME-driven shock. This may appear difficult because of strong magnetic mirroring due to the high degree of convergence of magnetic field lines toward the Sun (with a mirror ratio of η = B_⊙/B_CME ≫ 1) and thus extremely small loss cones (e.g., a few degress; Klein et al. 2018, see their Section 8.4.4). In an ideal scattering-free environment, this could potentially prevent particles from reaching the photosphere to produce γ-rays. In reality, however, the CME environment in general and the shock downstream region in particular are most likely highly turbulent, rendering a sufficiently short scattering mean free path which can continuously scatter particles into the loss cone, thus precipitating to the Sun. How fast and what fraction of the particles can reach the photosphere depend on the relative values of the duration ΔT of the emission and the escape time from the trap, T_esc. The latter depends on the mirror ratio η, the scattering time τ_sc, and the crossing time from the shock back to the Sun, τ_cross ∼ L/v_p, where v_p ∼ c is the velocity of GeV protons responsible for γ-rays and L = v_CMEΔT is the distance between the CME and the Sun. With some analytical and numerical treatments, Malyshkin & Kulsrud (2001) gave an approximate relationship between the escape time and the three variables η, τ_sc, and τ_cross (see Figure 2, Petrosian 2016):

$$\begin{eqnarray}&&{T}_{\mathrm{esc}}={\tau }_{\mathrm{cross}}\left(2\eta +\displaystyle \frac{{\tau }_{\mathrm{cross}}}{{\tau }_{\mathrm{sc}}}+\mathrm{ln}\eta \displaystyle \frac{{\tau }_{\mathrm{sc}}}{{\tau }_{\mathrm{cross}}}\right).\end{eqnarray} \tag{ 2 }$$

Recent numerical simulations by Effenberger & Petrosian (2018) have confirmed this relation. The upshot of this result is that for an isotropic pitch-angle distribution, ${T}_{\mathrm{esc}}\sim 2\eta {\tau }_{\mathrm{cross}}$ for τ_sc ∼ τ_cross and η ≫ 1, which means that T_esc/ΔT ∼ 2ηv_CME/c ∼ η/100 for a CME speed of v_CME = 1500 km s⁻¹. Therefore, for η ≲ 100, we have T_esc ≲ ΔT; i.e., a large fraction of the downstream GeV protons can reach the photosphere within the emission duration ΔT. We have tracked the history of the magnetic field strength along the field lines connecting the shock to the solar surface in our simulation and found that the median value of the mirror ratio η increases from ∼10 to ∼100 during the first hour (Figure 8(f)). Thus, for the first few hours, which is the duration of most Fermi events, the difficulty of particle transport back to the Sun due to magnetic mirroring can be overcome with a scattering mean free path of the order of a solar radius or a scattering time of τ_sc ∼ 2 s. For protons with gyro-frequency of ${{\rm{\Omega }}}_{p}\sim 15\cdot \left[\tfrac{B}{\mathrm{mG}}\right]$ Hz, and magnetic field of ∼100 mG (the typical average value on the shock surface in the first hour of simulation), this would require a fractional turbulence energy (δB/B)² ∼ 1/(Ω_pτ_sc) ∼ 3 × 10⁻⁴. According to Effenberger & Petrosian (2018), the situation is similar for a pancake pitch-angle distribution. But for a distribution beamed along the field lines, the escape time will be shorter and thus facilitate particle precipitation back to the Sun.