Modelling Solar Energetic Neutral Atoms from Solar Flares and CME-driven Shocks

Solar flares and coronal mass ejections (CMEs) are two of the most energetic processes in the solar system. Efficient particle acceleration can occur in both solar flares and at CME-driven shocks. Energetic protons accelerated at either CME-driven shocks or solar flares can precipitate down to the Sun’s surface or propagate into the interplanetary medium along open interplanetary magnetic field (IMF) lines. During their propagation, they can interact with ions and thermal neutral atoms in the solar atmosphere via charge exchange, and produce energetic neutral hydrogen atoms. Once produced, energetic neutral hydrogen atoms (hereafter referred as ENAs) do not feel solar magnetic field and propagate along straight lines. They are subject to loss processes wherein they lose the electron and become an energetic proton again. Because the density of the solar wind drops quickly with the heliocentric distance, and because the loss rate of ENAs is proportional to the solar wind density, ENAs reaching 20 Rs suffer no further loss. Since the IMF does not affect the propagation of ENA hydrogen, these ENAs therefore provide a powerful avenue in probing the acceleration processes and plasma properties of the underlying acceleration site.

Because the production cross section is small, the flux of ENAs at a distance of 1 au from the Sun can be extremely small. To date, only a few observational clues of ENAs accompanying SEP events were reported (Mewaldt et al., 2009; Mason et al., 2021). Mewaldt et al. (2009) reported $1.6$ to $5$ MeV energetic neutral atoms (ENAs) from STEREO-A/B observations. They inferred a power-law spectrum of $dJ/dE\sim E^{-2.46}$ accompanying an X9-class solar flare and suggested that these ENAs are produced via charge exchange of SEP protons with O⁶⁺ ions. Following (Mewaldt et al., 2009), Wang et al. (2014) performed a simulation and showed that sufficient counts of ENAs are expected for typical gradual SEP events where particles are accelerated at CME-driven shocks. This stimulated interests in observational effort. More recently Mason et al. (2021) examined $18$ SEP events with SAMPEX, and found indirect, but compelling evidence of solar ENAs near the geomagnetic equator at low altitudes where the geomagnetic field filters out all charged SEPs. This new insight also shed light on three previously reported puzzling $\sim$ MeV ion intensity increases that were also observed near the equatorial regions about $\sim 3$ hrs after the occurrence of the corresponding X-ray flares (Greenspan et al., 1999). The discovery of ENAs by STEREO, and confirmation from SAMPEX, shows that solar ENAs can be expected to accompany many large SEP events.

If ENAs can be detected in SEP events, one of the pressing questions would be where do they originate. Are they accelerated at a rather confined reconnection site at flares or at a broader shock front driven by CMEs? To answer such a question, we examine ENA productions in two different scenarios: CME-driven shock and large post flare loops in this work. A schematic of the two acceleration sites are shown in Figure 1. Note that in many large SEP events, CMEs and flares often occur together. However, the spatial extension of the flare is much smaller than the CME. Ions can be efficiently accelerated at both the flare site and the CME-driven shock front. In the case of CME-driven shocks, protons and ions are accelerated at the shock front via the first order Fermi acceleration mechanism. Once accelerated, they can escape upstream propagating along IMF, or trapped downstream for an extended period of time. They may precipitate down to the solar surface, causing, for example, long duration gamma-ray events (Share et al., 2018; Jin et al., 2018). In the case of flares, particles can be accelerated at the reconnection exhausts and in solar flare loops (Petrosian, 2012; Ryan, 2000) by e.g. the second order Fermi acceleration mechanism. Continued magnetic reconnection can lead to a rising of the post-flare loops (West & Seaton, 2015). Accelerated particles may be trapped in post-flare loops for very long period of time, serving as an alternative candidate for the long-duration gamma ray events (Ryan, 2000; de Nolfo et al., 2019).

We note that in large SEP events it is possible that flare accelerated particles can be re-accelerated at the accompanying CME-driven shocks (Li & Zank, 2005; Petrosian, 2012). Simulations by Li & Zank (2005) showed that depending on if the observer is magnetically connected to the flare and/or the shock surface, the characteristics of the ion time profiles differ and may show two peaks as reported in (Cane et al., 2003). However, because the presence of solar wind MHD turbulence can affect the propagation of charged ions, the interpretation of 1-au ion observations is often complicated. ENAs, with ballistic propagation, do not follow IMF and are not affected by the solar wind MHD turbulence. Therefore, ENA observations with enough angular resolution can clearly distinguish ENAs from flare sites and those from much broader CME-driven shocks.

We examine solar ENA production from CME-dirven shocks in section 2 and from solar flare loops in section 3. Production and loss processes of ENAs are discussed in Appendix A.

In this section, we consider the observation of ENA particles generated at a propagating CME-driven shock. The first ENA simulation was done by Wang et al. (2014) who simulated a CME-driven shock from a side-on orientation and suggested that the observed flux in (Mewaldt et al., 2009) is consistent with ENA production at a CME-driven shock. More recently, following the work of (Wang et al., 2014), Wang et al. (2022) examined a variety cases with different CME speeds, open angles, and CME propagation directions. They also examined the effect of solar wind density variation near the Sun on the production of ENAs. These authors found similar results as Wang et al. (2014).

Here we reexamine the case considered in (Wang et al., 2014) and include another two cases with different CME propagation directions to obtain an estimate of ENA flux range at 1 au. Our treatment is similar to our previous work (Wang et al., 2014) but with a few differences. As in (Wang et al., 2014), we assume protons are accelerated at the shock and then distributed uniformly downstream of the shock. This is based on the DSA mechanism and has been adopted in our previous large SEP event simulations (Li et al., 2003, 2005, 2012b, 2021). Since the turbulence downstream of the shock is a lot stronger than that upstream of the shock (see e.g. (Lee, 1983; Zank et al., 2000; Li et al., 2003)), accelerated particles can be kept downstream of the shock for a long period of time. In (Wang et al., 2014), we assumed there is no leakage of accelerated particles from downstream of the shock. This was mostly for simplicity since accelerated particles can precipitate back to the sun along open field lines. Indeed, Jin et al. (2018) has explored the possibility that the long duration gamma ray events are due to shock acceleration protons. In such a scenario, accelerated protons downstream of the shock can steadily precipitate to the solar surface. Therefore, in this work, we include a decay of the accelerated protons downstream of the shock. As an estimate of the decay time, we refer to Li et al. (2012a), who, from a statistical study of twin-CME events, suggested that a decay time scale of the turbulence in large SEP events is around $9$-$13$ hours. We use a decay time $\tau=10$ hours in this work. We also set our inner boundary at $r=1.02R_{s}$, which differs from that used in (Wang et al., 2014), $1.5R_{s}$. We further improve the treatment of ENA propagation from downstream of the shock to the observer. In (Wang et al., 2014), downstream medium was divided into shells and ENAs produced in individual shells are assumed to propagate to the observer all from the shell center. This is refined in our current work. We now divide the downstream region of the shock into multiple parcels, as shown in the left panel of Figure 2. ENAs are produced and followed in individual parcels. Since ENAs in different parcels propagate to the observer along different paths, our current treatment will lead to a more accurate survival probability computation. Finally, a correction factor $cos(\theta)$ to the flux expression, equation (4) in (Wang et al., 2014) is included, see equation (A3).

Figure 2 shows the configuration of the ENA production process for the CME shock case. The left panel depicts the base case: the observer locates at 1 au along the $X$ axis and the CME is propagating to the right along the $+Y$ direction, i.e., $\phi=90^{\circ}$ where $\phi$ is the angle between the sun-observer line and the CME propagation direction. The plasma downstream of the shock is divided into multiple parcels. ENA production is followed in these parcels. ENAs produced in these parcels can propagate along straight lines to the observer. These trajectories differ for different parcels, and lead to different survival probabilities. Right panel of Figure 2 shows two other cases with different CME propagation directions. In case II, the CME propagates toward the observer with $\phi=0^{\circ}$. In case III, the CME propagates $45^{\circ}$ off from the $+X$ and $+Y$ directions, i.e. $\phi=45^{\circ}$. For our simulation, the shock has a constant speed of $V_{sh}=1500$ km/s and a constant compression ratio of $s=3.5$. The open angle of the shock is $60^{\circ}$ and the shock is followed up to $30R_{s}$. As in the flare ENA case, we use the Leblanc model (Leblanc et al., 1998) to compute the solar wind density.

Figure 3 plots the time profiles and the fluence of ENAs for the three cases shown in Figure 2. The upper left, upper right, and lower left panels show the time profiles for the base case, case II and case III, respectively. For all three cases, ENAs of $11$ energies are considered. The three time profiles are similar. Consider the base case (upper left panel). The x-axis is the time after shock initiation, in unit of $10$ minutes; and the y-axis is the ENA flux at the observer, in unit of #/(cm${}^{2}\cdot$ sec $\cdot$ keV). The observer first see the $20$ MeV ENAs arriving $\sim 40$ minutes after the shock initiation. The flux can reach $7*10^{-3}$ cm${}^{2}\cdot$ sec $\cdot$ keV. It then decreases, reflecting the fact that the density of energetic protons decreases with time as the shock propagates out. As the ENA energy become smaller, their first arrival times become later and the flux increases with decreasing energy, till $E=0.75$ MeV. Below $E=0.75$ MeV the flux shows a more plateau feature and drops slightly. This behaviour is due to the energy dependence of the charge exchange cross sections that are responsible for the ENA production. See Figure 6 in Appendix A. Comparing to the base case, cases II and III are comparable and show larger fluxes than the base case. This is easily understood from Figure 2 because the ENAs produced in these two cases travel shorter distances and through less dense solar atmosphere to the observer and consequently have larger survival probabilities. The lower right panel of Figure 3 plots the fluence for these three cases. Note the relatively plateau-like behavior below $E=0.75$ MeV, which is the consequence of the energy dependence of the relevant charge exchange cross section. Above $1$ MeV, the ENA fluence spectrum shown here is comparable to that inferred in Mewaldt et al. (2009). The general shape of the CME shock ENA fluence is similar to the parent energetic ion spectrum which is a power law. This is in stark contrast to the flare ENA case (see next section) where the ENA fluence does not resemble the parent energetic ion spectra.

We examine ENA production by solar flares in this section. Both electrons and ions are efficiently accelerated at solar flares, and the accelerated electrons and ions lead to the emission of hard X-rays and gamma rays. It is generally accepted that acceleration may occur at reconnection current sheets and/or by turbulence in the flare loops. Observations of hard X-ray and gamma rays suggest that the accelerated electron and ion spectra can be approximated by a power law. Power law like spectra are supported by earlier theoretical works by (Miller & Roberts, 1995; Petrosian, 2012), where ions are accelerated in flare loops by MHD turbulence via second order Fermi acceleration. More recent PIC simulations of ions in flare reconnection site also found a power law spectrum (Zhang et al., 2021).

Once accelerated, ions precipitate down to the solar surface along post-flare loops. The density of a post-flare loop can be constrained by free–free continuum emission for hot loops. In a recent work, Jejčič et al. (2018) reported an electron density as high as $10^{13}$ cm^-3, $10$ to $100$ higher than that at typical flare loops. At a density of $\sim 10^{11}$ cm^-3, ENAs can be easily produced in these loops. Once produced, ENAs are not constrained in the loops and can propagate in all directions. However, if the loops are low, the solar atmosphere density in the surrounding environment can be too dense to allow these ENAs to escape from the Sun. Therefore to observe flare ENAs, the flare loops must be high. The height of flare loops can be estimated from the looptop hard X-ray observations. A recent study of looptop hard X-ray source of solar flares (Effenberger et al., 2017) showed that the height of a typical flare ranges from $10$ to $50$ Mm. If ions are accelerated at and below this height at flares, no ENAs can survive as they propagate out. However, using the Sun Watcher with Active Pixels (SWAP) EUV imaging solar telescope, West & Seaton (2015) examined an M2.2 flare which occurred on 2014 October 14 and found that the post-flare loops were long-lasting, and reached a height of over $400$ Mm ($\geq 0.5\ {{R}_{\odot}}$) $\sim 48$ hours after the eruption. West & Seaton (2015) argued that the giant arches in this event are similar to ordinary post-flare loops and are the results of a long-lasting ($48$ hours) magnetic reconnection occurred along a large-scale current sheet (Forbes & Lin, 2000). This continuous magnetic reconnection provides the energy source to heat the loop and can accelerate particles. Besides magnetic reconnection, turbulence inside the loop can also lead to stochastic acceleration of ions (Ryan, 2000). We note that the magnetic reconnection at the current sheet and the enhanced turbulence inside the large postflare loop may be intimately related. In a recent work by Cheng et al. (2018), the authors examined the 2017 09 10 flare and showed that a Kolmogorov-like turbulence spectrum can develop in the current sheet above the flare loops. Presence of such a turbulence implies that particles can be accelerated in the turbulent current sheet in a similar way as in flare loops through a second order Fermi acceleration process (Miller & Roberts, 1995; Petrosian, 2012). The spatial extension of the current sheet is similar to the flare loops that are beneath it (Cheng et al., 2018; French et al., 2019), but the density of the current sheet is, however, smaller than the density in the postflare loops. Indeed French et al. (2019) concludes that the density in the current sheet is $\sim<10^{10}$/cm³. This is $100$ times smaller than the density inferred in the lower flare loop as reported in (Jejčič et al., 2018), and is $10$ times smaller than what we assume for the post flare loops, $10^{11}$/cm³. The ENA production in these current sheets is therefore much smaller than in postflare loops. So we only consider ENA producion in postflare loops in this work. However, we remark that these turbulent current sheets can be potential sites of ENA production, and if future ENA probes have high enough sensitivities, it is possible to obtain direct observations of these current sheets through ENA observations. A cartoon showing these post flare loops are shown in Figure 4.

Continuous acceleration, as suggested by Ryan (2000), has been identified as a possible scenario for the long duration gamma ray events de Nolfo et al. (2019). Long duration gamma ray events are not uncommon. Recently Share et al. (2018) examined $\sim 30$ long duration gamma ray events and found that the energy spectral indices of $>300$ MeV proton producing gamma rays range from $2.5$ to $6.5$, similar to typical flare events. In a recent study, de Nolfo et al. (2019) compared the gamma-ray-producing proton numbers with the in-situ SEP proton numbers in long duration gamma ray flares and found a poor correlation. Their study supports the continuous acceleration in the post-flare loop scenario, as suggested by Ryan (2000). We point out that the event reported in (West & Seaton, 2015), despite having large post-flare loops, was not a long duration gamma ray event. This is possible if particles are not accelerated to high enough energies ($\sim 100$ MeV/nuc) to produce gamma rays.

We now examine ENAs from post-flare loops. We model the post-flare loops as semi-circle tubes. We assume that the loop has a height (radius) of $h(t)$, which increases with time. We assume the starting height of the post flare loop is $0.04R_{s}$ ($\sim 28$Mm), and a rising rate of $V_{r}=3$ km/s (West & Seaton, 2015). This gives a height of $H=0.22$, $0.41$, and $0.60$ $R_{s}$ when $t=$ 12, 24, and 36 hours, respectively. The cross section of the tube can be assumed to be a circle with a radius as $a$. One can take $a$ to be $\sim 700$ km, which is comparable to the half width for a typical flare ribbon. However, as we will see below, the ENA production depends on the total number of accelerated protons and does not depend on the choice of $a$ and the number of loops we consider.

We also assume a constant proton density inside the flare loop. By way of example, we assume a loop density of $10^{11}$ cm^-3. This is smaller than that obtained in (Jejčič et al., 2018), but larger than the density at the solar surface, which is $\sim 10^{9-10}$ cm^-3. As a simplification, we assume the acceleration process (Ryan, 2000) is time independent and the production rate of energetic protons, $\alpha$, is a constant during the rising phase of the post-flare loop. We denote the duration of the rising phase to be $T$, and the total number of accelerated particle $N_{0}=\alpha T$. Once accelerated these particles can precipitate to the solar surface. We model this as a loss process with an energy-independent decay time $\tau$. The total number of accelerated particles $N(t)$ in the loop is given by,

where $\theta(t)$ is the Heaviside function. The solution of equation (1) is,

In equation (2), $N_{0}$ can be constrained from the following consideration. In the long duration gamma ray events examined by (Share et al., 2018), the authors inferred that accelerated particles at high energies ($>$300 MeV) in the loops is about $0.01$ to $0.5$ of that of the accompanying SEP events, presumably accelerated at the CME-driven shocks. Assuming this ratio is energy independent, then one can estimate the range of $N_{0}$ from the CME-driven shock case. Alternatively, one can estimate $N_{0}$ from an energy budget point of view. In a study of the CME/Flare Energy Budget for two (Emslie et al., 2005), and subsequently for $38$ large SEP events, Emslie et al. (2012) found that the energy budget for $>1$ MeV flare ions can reach $\epsilon\sim 4*10^{31}$-$10^{32}$ erg, which can be comparable and even larger than those observed in-situ. In this work, we estimate $N_{0}$ by assuming the total energy for the accelerated particles ($>1$ MeV) is $\epsilon=10^{31}$ erg. With a source spectrum of the accelerated protons given by,

where $E_{0}$ is the injection energy, $E_{b}$ is the break energy, $\gamma_{1}=2.5$ is the spectral index at energies below $E_{b}$ and $\gamma_{2}=5.5$ is the spectral index at energies above $E_{b}$. This gives,

For a choice of $E_{0}=0.02$ MeV, $E_{b}=30$ MeV, $\gamma_{1}=2.5$ and $\gamma_{2}=5.5$, we find $N_{0}=9*10^{38}$. Equation (3), together with equations (2) and  (4) describe the energetic proton source, as a function of time, for the ENAs inside the post-flare loop.

One can now compute the production of ENAs and obtain the time profiles and fluence of ENAs as observed at 1 au. We consider three cases with $T=12$, $24$, and $36$ hrs, corresponding to a final loop height of $H=0.22$, $0.41$, and $0.60$ R_s, respectively. In all cases $\tau=3$ hr. We further assume that the flare locates at $\phi=0$ degree, i.e. in a face-on situation. For other viewing angles, the results are qualitatively similar.

Figure 5 plots the time profiles and fluence of the flare ENAs. The upper left, upper right and lower left panels are time profiles for the three choices of the final flare loop heights. Seven energies are considered. These are $1.0$, $1.5$, $2$, $5$, $10$, $15$, and $20$ MeVs. As can be seen from these panels, high energy ENAs arrive earlier due to a short propagation time from the Sun to 1 au. In all three panels, the peak of the time profiles occur shortly after the loops reach the maximum height. The energy dependence of the peak intensity (and the fluence, see the lower right panel) strongly depends on the loop height. If the loop height is $0.22R_{s}$ (upper left panel), the peak intensity of $T=2$ MeV ENAs is $5$ orders of magnitude smaller than that of the $T=15$ MeV ENAs. Furthermore, there is no ENAs with $T<2$ MeV. In comparison, when the loop height is $0.41$ or $0.6$ R_s, the peak intensity of $T=2$ MeV ENAs is similar to that of the $T=15$ MeV ENAs. This energy dependence can be also seen from the fluence plot shown in the lower right panel of Figure 5. When the loop height is $0.60R_{s}$, the fluence has a maximum $\sim 800$/cm² at $T=2$ MeV, and at $T=20$ MeV, the fluence is about $10$. When the loop height is $0.22R_{s}$, however, the fluence of 2 MeV ENAs drop by a factor of $4*10^{8}$ to $\sim 2*10^{-6}$/cm². In comparison, the flunece of $20$ MeV ENA drops only by a factor of $50$, to $0.2$/cm². This big difference of ENA fluence at 1 au for different flare loop height is due to efficient loss of ENA close to the Sun. Although plenty of ENAs are produced in the flare loop, they can not escape the high density solar atmosphere if the flare loop is not high enough. Note that during the eruption phase of solar flares, the height of flare loops, as seen from X-ray imaging, is a lot smaller than $0.22R_{s}$ (Effenberger et al., 2017), therefore we expect no ENAs during the eruption phase of solar flares. However, large post flare loops, as those reported in (West & Seaton, 2015), can reach $0.5R_{s}$. Our calculations show that there will be clear ENA signals from such a flare. We do note that the absolute amplitude and the shape of the ENA fluence depend on the solar atmosphere density model as well as the relevant charge exchange cross sections (see Appendix A). Nevertheless, because the ENA fluence, and in particular, its energy dependence, sensitively depend on the flare loop height, one can use the ENA fluence as a probe of the flare loop height. We point out that these large post flare loops may not be common. Consequently, flare ENAs may not be common either. Note that both the time profiles and the fluence for flare ENAs shown in Figure 5 are vastly different from their counterparts in shock accelerated ENAs shown in Figure 3. This suggests that one can use ENA observations to discern if the parent energetic ions are accelerated at CME-driven shock or at solar flares.

Understanding the underlying particle acceleration process in large SEP events has been one of the central problems in heliophysics research. With only in-situ observations of energetic ions, questions such as the relative roles of magnetic reconnection in flares vs shock acceleration at CME shocks, and how to discern the effects of acceleration from that of transport, can be very hard to answer. In part, this is because our basic understanding of the near-Sun conditions and the physical processes involved in the production of SEP events is hampered by our inability to make direct measurements near the acceleration sites and to remove the effects of transport. ENA observations can significantly advance our understanding of SEP acceleration at its source because ENAs do not interact with IMF and is not affected by the transport effect.

In this paper, we examine the production of ENAs at CME-driven shock fronts and in solar flares. We compute the time profiles and fluence of ENAs for these two scenarios. Our calculations suggest that in large SEP events where ions are efficiently accelerated at CME-driven shocks, ENAs are copiously produced behind the shock. At 1 au the flux of these ENAs are at a level that can be readily measured by a dedicated ENA detector. ENAs can also be produced in flares where large scale and high postflare loops exist. The time profiles and fluence of ENAs for these two scenarios differ considerably. This offers us an opportunity to constrain the underlying particle acceleration process via ENA observations. Our work also forms a theoretical basis for interpreting future ENA observations.