Energetic Neutral Atoms from Solar Energetic Particles due to Shocks: Inclusion of Upstream Particles

Many aspects of solar energetic particles are not well understood, including their acceleration mechanism. There has been recent interest in the potential of energetic neutral atoms (ENAs) as remote probes of solar energetic particles (SEPs) and their acceleration. The single accidental observation (in physical units) has been modeled as accelerated by a coronal mass ejection (CME)-driven shock by several authors, all of whom have assumed that the upstream component of the shock can be ignored. In this article, we relax this assumption and model the flux of ENAs at 1 au due to a CME-driven shock with an upstream component. We show the effect of varying parameters of the shock acceleration model, specifically α, the exponent of the power law in momentum of the mean free path, and η, a measure of the relative turbulence level. The main result is that including the upstream component significantly increases the flux at 1 au for typically assumed parameters in the energy range of the STEREO observation. We also derive the form of the ENA transport equation that we used in this study. These results enable a better understanding of potential observations of ENAs due to SEPs.


Introduction
Many aspects of solar energetic particles are not well understood, including their acceleration mechanism.There has been recent interest in the potential of energetic neutral atoms (ENAs) as remote probes of solar energetic particles (SEPs) and their acceleration.The single observation that has been reported, at least in physical units (Mewaldt et al. 2009)-there have been recent exciting results reported using SAMPEX data, but these are only shown in terms of instrument count rates (Mason et al. 2021)-has been modeled as accelerated by a coronal mass ejection (CME)-driven shock (Wang et al. 2014;W14 hereafter), which was recently enhanced by a parameter study using a similar model (Wang et al. 2022), in which they varied the CME speed, direction, and opening angle, as well as extended to include flare acceleration (Li et al. 2023).One of the assumptions common to all these studies is that the upstream component of the shock can be ignored.In this article, we relax this assumption and model the flux of ENAs at 1 au due to a CME-driven shock with an upstream component using a Monte Carlo test particle simulation.
Monte Carlo test particle simulations have been used in multiple publications to study shock acceleration in a variety of environments, including interplanetary shocks (Baring et al. 1997;Summerlin & Baring 2006), the solar wind termination shock (Ellison et al. 1999), supernova remnants (Baring et al. 1999;Baring & Summerlin 2007), and in the regime of highly relativistic shocks that are generally encountered in extragalactic contexts (Ellison et al. 1990b;Ellison & Double 2004;Stecker et al. 2007).

Shock Model Overview
We first consider a shock near the Sun traveling at a constant speed u s radially away from the Sun.As in W14, we presume the center of the shock has a trajectory in the ecliptic plane and is perpendicular to the Earth-Sun line.The shape of the shocked region is conical, with a spherical cap at the position of the shock, and a spherical indent at 1.5 R e (see Figure 1 of Wang et al. 2022).The full opening angle of the cone is 60°, again, as in W14.We take u s to be 1800 km s −1 as in W14.Since inside 40 R e , the magnetic field can be approximated accurately by its radial component only, the shock is assumed to be parallel.
For every test particle simulation (see Figure 1 for a general outline of the process), we inject a cold distribution of protons traveling at 1350 km s −1 relative to the shock as the incident distribution function.This simulates a 1800 km s −1 shock overtaking a 450 km s −1 solar wind.The distribution is chosen to be cold so as to more easily distinguish the incident distribution from the upstream precursors.Including and varying the temperature are easily done as well, but its net effect is on injection efficiency rather than the resulting power-law distribution, which is not the focus of this study and the same reason why only parallel shocks are considered.Both of these effects will alter somewhat the normalization of the resulting power law but not its index.The particles are then scattered according to a prescribed phenomenology, described in more detail below.Particle positions and velocities are tracked, and distribution functions are calculated at given energies and locations once the simulation has been completed.Technically, the code can be made to output a distribution function at any synchronized time step; however, for the present work, we only create a distribution function once the code has run to completion.
The particular implementation used in the present work is based on Summerlin & Baring (2006).It follows closely the approach of Bell (1978), which was extended to the relativistic regime by Peacock (1981).New capabilities of the Monte Carlo test particle code include additional particle removal conditions (see below), as well as the inclusion of statistical errors for the distribution function.

Scattering
The code models particle gyration about bulk magnetic fields in convecting fluid flows, while having their trajectories perturbed by embedded hydromagnetic turbulence.The effects of magnetic turbulence are simulated by specifying a local fluid frame mean free path for particle diffusion, given by where r g = pc/(qB) is the gyroradius of an ion or electron of momentum p = mv, mass m, and charge q in a magnetic field of strength B = |B|, where c is the speed of light in vacuum.Also, r g1 = mu 1⊥ c/(qB) is the gyroradius of an ion with a speed equal to the far upstream flow speed normal to the shock plane, u 1⊥ , where ⊥ denotes the direction normal to this plane.Without loss of generality, the mean free path scale is set proportional to r g1 with a constant of proportionality η, while the momentum dependence of the mean free path is determined by α.See Ellison et al. (1990a), Mason et al. (1983), andGiacalone et al. (1992) for discussions about the microphysical expectations for α.The total magnetic field strength, B, can be approximated inside 40 R e by its radial component only, for typical solar wind speed: B = B 0 (10 R e /r) 2 , where we have used B 0 = 1.83 × 10 −6 T (Zank et al. 2000).
The simplest invocation of scattering is to isotropize the fluid frame momentum of a scattered particle over the surface of a sphere in momentum space (Ellison et al. 1990b).This is referred to as large angle scattering (LAS) and physically corresponds to large magnetic disturbances that completely disrupt the trajectories of particles.To model moderate or even smaller disturbances, each scattering event can be restricted to a much smaller solid angle, i.e., it can be isotropic on a conical sector of a momentum sphere.The angular extent of this spherical sector dq max becomes an additional parameter for diffusion.In this case, multiple scattering events are required to realize a full mean free path.The relationship between dq max and λ was originally developed in Ellison et al. (1990b), but is more succinctly presented in Ellison & Double (2004) via g max where r g is the gyroradius, and N is the number of times per gyroradius the particle is scattered.The limit of small angle scattering (SAS) corresponds to N ? 1, for which the increment δp in momentum for a single scattering event satisfies | | d dq p p max .For nonrelativistic shocks and the scattering mechanisms used here, there is no material difference between SAS and LAS; therefore LAS is used as it is less computationally expensive.
The sample particle trajectory in Figure 2 depicts the salient features of the particle motions described above.In the example in Figure 2, particles gyrate in the local magnetic field until their velocity vector is isotropized, modeling an encounter with a large electromagnetic disturbance.Particles crossing the shock experience a sudden change in the direction of the flow and magnetic field vectors downstream of the shock.The kinks in the flow vector and the magnetic field vector are a result of the Rankine-Hugoniot shock jump conditions.The color coding represents changes in the local fluid frame momentum that occur at each shock crossing.The overall energy gain can be seen by comparing the gyroradius (r g = p ⊥ /qB) of the particle at the beginning of its trajectory (top) to the gyroradius at the time the particle trajectory tracking ceases (bottom right).
Cross-field diffusion emerges naturally from this scattering mechanism since, at every scattering, the particle's momentum vector is shifted in the local fluid frame, with the resulting effect that the gyrocenter of the particle is shifted by a distance of order r g , in the plane orthogonal to the local field.Transport perpendicular to the field is then governed by a kinetic theory description, so that the ratio of the spatial diffusion coefficients parallel (κ || = λv/3) and perpendicular (κ ⊥ ) to the magnetic field is given by κ ⊥ /κ || = 1/(1 + η 2 ) (see Forman et al. 1974 andEllison et al. 1995 for detailed expositions).Hence, η couples directly to the amount of cross-field diffusion and is a measure of the level of turbulence in the system, i.e., is an indicator of 〈δB/B〉.The quasiisotropic diffusion case of η = 1 constitutes the Bohm diffusion limit, presumably corresponding to 〈δB/B〉 ≈ 1.However, as the scattering takes place in the local fluid frame, the energy of the particle in that frame is constant.This is seen in Figure 2 where particles only gain or lose energy relative to the local plasma frame at the shock crossing.

Particle Escape
There are two conditions that can cause a particle to be removed from the simulation.The first condition is the distance from the shock.We distinguish between parallel to and perpendicular to the plane of the shock.In the downstream direction (perpendicular to the shock), particles are removed from the simulation via an analytic formula originally developed by Bell (1978) and later shown to be applicable to relativistic shocks by Peacock (1981)see Equation (3) of Summerlin & Baring (2006).In the directions parallel to the shock, particles are removed after a fixed distance from the shock center.This distance is set to the radial distance of the shock from the Sun, r s , which implies that the angular extent of the shock is ∼57°, consistent with the assumed 60°of W14 and in the ENA calculation below.In the case of spatial location, the removal is meant to physically correspond to escape from the shock.Although escaped particles technically exist in the heliosphere, for the purposes of the simulation and for simplicity, they do not count toward the accelerated proton distribution and thus cannot become ENAs.
The second condition for removal is particle age.The age of each particle is tracked in the simulation, and once its age exceeds a certain value, it is removed from the simulation.In this work, the value is set to the lifetime of the shock.Physically, this condition corresponds to the shock reaching its desired location and thus the particle having no more time to be accelerated.The lifetime is taken to be the travel time from the location of the shock's formation to the current location of the shock.Written as a function of the shock's radial distance from the Sun, r s , this lifetime is τ = (r s − R e )/u s .We use a lower limit of R e to ensure the distributions at 1.5 R e have a nonzero lifetime.

The Accelerated Proton Distribution
The parameter α, the exponent of the power law in momentum of the mean free path; η, the overall turbulence level; and r c , the compression ratio of the shock, were all explored.Values of η used are 10 0 , 10 1 , 10 2 , and 10 3 ; values of α used are 1, 1.5, and 2; and values of r c used are 3, 3.5, and 4, corresponding to power laws in momentum of 4.5, 4.2, and 4, respectively.Every combination of the values for each parameter was simulated, at 10 locations between 1.5 and 40 R e (logarithmically spaced), resulting in a total of 360 simulations.
As an example, in Figure 3 we have plotted the distribution function at the shock, f 0 , for a single parameter set (η = 1, α = 1, r c = 4) as a function of the distance from the Sun of the shock, r s .The log of the values between the simulated distance points has been interpolated linearly as a function of the log of the distance, i.e., a linear interpolation was performed in loglog space.Notice that the high-energy rollover changes greatly with shock distance, mostly decreasing, whereas the lowenergy portions maintain a consistent power law.This is understood by the fact that decreasing the upstream magnetic field increases the gyroradius of the particles, which, according to Equation (1), should increase the mean free path of the particles.Therefore, the particles travel farther and take longer to return to the shock, allowing them to escape the shock or just not be accelerated by the time the simulation ends.This effect is mitigated somewhat by the fact that the convection time, τ, increases with distance, being roughly proportional to r at large r; however, it increases more slowly than B decreases (∝r −2 ).
In Figure 4, we show an example of the upstream portion of the distribution function at a given location (r s = 1.5 R e ) and for a given set of shock and scattering parameters (r c = 4, α = 1, η = 1), plotted as a function of radial distance from the Sun, r, relative to the shock location, r s , shifted by 1 so that the shock is at 0. We then fit these points with the following functional form, sometimes referred to as a selfsimilar distribution, where ¢ = r r r 1 s , l is a dimensionless parameter that controls the falloff distance, q is a unitless parameter that controls how sharp the falloff is (a higher q corresponds to a sharper fall-off; q = 1 is a Gaussian; q = 1/2 is an exponential), and f 0 is taken directly from the simulation as the value of the distribution function at the shock.Therefore, only the parameters l and q were allowed to vary. Figure 4 also shows the best fits of Equation (3) to the simulation results obtained through an automated reduced χ 2 minimization tool (Markwardt 2009).The uncertainties used for calculating the reduced χ 2 were taken from the simulation output.We repeated this process for the rest of the 360 simulations at the 16 energies shown in Figure 4.
The best-fit parameters of the fits to the simulation of Figure 4 are shown in Figure 5, including f 0 , even though it was not allowed to vary in the fitting process.The parameter l increases with increasing energy, meaning that higher energy particles have an easier time moving upstream, as one would expect.In Figure 5, the energy dependence is roughly l ∝ E at moderate energies.In Figure 5, the parameter q is between 0.3 and 0.5, increasing monotonically with energy, though much more slowly than l.In general, these trends remain for all simulation runs, though the energy dependence of l changes, and, in general, q is always less than 1 but greater than 0.1.
For use in Section 2.5, we need to formally include the full spatial and time dependence of the accelerated proton distribution: The distribution function of protons at the shock (open circle symbols), f 0 ("f0"), vs. distance from the Sun of the shock, r s , with interpolation between the points (straight line segments).Colors denote energy.The parameters of the shock model are r c , the shock compression ratio, and α ("alpha") and η ("eta"), the exponent of the momentum power law and normalization of the mean free path, respectively, as in Equation (1).
where Θ(ξ) is the Heaviside step function; r is the heliocentric distance and is equal to + + x y z 2 2 2 ; r s is the shock's location and is equal to r s,0 + u s t; the initial shock location r s,0 is equal to 1.5 R e , as in W14; f d is the downstream distribution and is simply equal to f 0 ; n(r s ) is the solar wind density at the location of the shock and is calculated using Equation (2) in W14; n 0 is the density at which the shock simulations were run (see Section 3); and the cone angle of the shock is defined by y = y r cos .In the above, we have used the same coordinate system as Figure 1 of Wang et al. (2022), where the Sun is at the origin, the observer is at 1 au in the positive x-direction, and the CME's trajectory is in the positive y-direction.

ENA Calculation
The next step, and the ultimate goal of this work, is to compute the energetic neutral atom (ENA) differential flux, J,   4 (symbols) as a function of energy in keV; solid lines are guides to the eye and simply connect the points.Note that f 0 was not allowed to vary in the fit, only l and q; f 0 is the distribution function at the shock and comes directly from the simulation.The dashed line is the power law predicted by diffusive shock acceleration theory.
at Earth (see Appendix A): where m H is the mass of hydrogen.Starting from the general formula, Equation (A2), we have converted to flux units via J = fv 2 /m H , dropped the explicit velocity dependence, and converted the velocity into spherical coordinates: ( ) q q f q f = v v cos , sin cos , sin sin .In the actual results reported below, we have taken the observer's location to be (x, y, z) = (1 au, 0, 0), and set x 0 equal to −50 R e .Theoretically, x 0 should be − ∞; however, in this case, because the source region is confined to within 40 R e of the Sun, the range of integration can be shortened for numerical convenience.We have also interpolated linearly in log-log space the best-fit parameters of the fit to the shock simulation results between the 10 spatial locations of the shock for better resolution, ignoring points where the fit could not converge for whatever reason (this is quite rare).
In Equations ( 5)-( 6), we have used the production and loss functions as specified by Equations (A4) and (A5), respectively.In those equations, the specific form used for P CX and L CX is given by Equations (A9) and (A10), respectively.As discussed in Appendix A, we also use the charge-exchange form for the impact ionization loss terms.For photoionization loss, we use the same values as in W14 (D'Amicis et al. 2007), the form of which is simply 2 , where r is the heliocentric distance.We use the same ratios between species as in W14.Finally, the cross sections used, which differ somewhat from W14, are discussed in Appendix B. An example of the survival probability is shown in Figure 6.Note that because the survival probability contains only heliocentric symmetric functions, e.g., the solar wind density, n(r), the f dependence drops out, i.e., ( x tan 2 2 .In the production function, this is not the case, since the CME trajectory is along the y-direction.
Since an observer at Earth may not have the resolution to detect the incoming direction of particles, and in order to condense the information contained in the four-dimensional quantity J, we take an average over all possible viewing angles: , , , sin d d .7 , 0 2 0 An example of the result of performing this average is shown in Figure 7.One can see that the upstream component is highly peaked in time, as to the downstream component, and arrives slightly before the downstream component.The highly peaked property is qualitatively consistent with the time series data shown in Mewaldt et al. (2009).Note that the upstream component is greater (at its peak) than the downstream component at every energy simulated by more than an order of magnitude at most energies.Note that the plot artificially cuts off at the lower flux end to better highlight the energy range of the STEREO measurements.Finally, in order to compare with data reported in Mewaldt et al. (2009), we compute the angle-averaged fluence, which is simply the time integral of the angle-averaged differential flux, <J > , over the duration of the event.This result is shown in Figure 8.Note the excellent agreement of the model with the data.The most significant result to observe in this plot is that the upstream component is the dominant component by more Figure 6.Survival probability, h, as a function of the location of birth along the x-coordinate in units of R e .Colors represent energy in keV.The parameter θ ("theta") is the polar angle of particle velocity relative to the x-direction.Note that there is no dependence of h on the parameter f. than an order of magnitude in the energy range of the STEREO measurement.We believe this to be due to the fact that as particles diffuse upstream, they get farther from the Sun, where the density is lower and thus the probability of reionization is lower, increasing the chance of the particle escaping and reaching 1 au to be observed.Lastly, similar to W14, oxygen is the dominant source of charge exchange beyond ∼1 MeV.
In the above, the example set of shock parameters shown (η = 1, α = 1, r c = 4) represents one of several that have excellent agreement with the data reported in Mewaldt et al. (2009).Figure 9 shows the fluence that arises when η and α are varied while keeping r c fixed at 4. The curves for r c equal to 3.5 and 3 are similar, except that the overall fluence drops by factors of approximately 2 and 5, respectively.

Discussion
Note that the simulations do not naturally produce an absolute distribution function, only a relative one between energies and spatial locations.That is to say, the density needs to be normalized.To compute the normalization, we compared our spectra with the analytical prediction in W14, Equation (1).7 as a function of energy in keV.Shown are the upstream components (dashed curves) and downstream components (thin solid curves), as well as the contribution from oxygen (red), carbon (green), and hydrogen (blue); the sum of all components and charge exchange sources is also shown (bold black curve).Data (filled circular symbols) are adapted from Mewaldt et al. (2009).The power-law fit from that work is also shown (dotted line).Shock parameters are the same as in Figure 3.
A small dependence on r c was found and applied, and the overall normalization was then adjusted to best match the fluence data at r c = 4, η = 1 = α; the adjustment was to multiply everything by a factor of 5.5, assuming an upstream solar wind density of 1 cm −3 .This step is well justified since there are multiple free parameters in the analytical estimate that could account for such a small change.
The profile of the shock with distance is not extremely realistic, since most simulations suggest an evolution of the shock angle and compression ratio with distance (even inside 40 R e ), as well as variation across the surface of the CME (Li et al. 2021;Lario et al. 2020).We chose to keep the shock properties constant with distance and across the shock surface so as to better compare with W14.This issue should be investigated in a future study.
The only parameter in the above that is different from W14 is the shock compression ratio, r c .We have shown a value of 4 as opposed to the value of 3.5 used in W14.As stated above, 4 was a better fit to the data than 3.5; as noted above, the r c = 3.5 curves were mostly simply lower in overall flux by a factor ∼2.However, if the background solar wind density were different, 3.5 could well have been a good fit instead.The model is indeed fairly sensitive to the solar wind density, and the solar wind density at 1 au is known to vary considerably, up to a factor of 10.From our new understanding of solar wind plasma close to the Sun from Parker Solar Probe, we could anticipate even larger variations in density at the region considered in this study.Equivalently, since the contribution due to charge exchange with oxygen dominates the energy spectrum for the energy range of the measurements, a higher ratio of oxygen to hydrogen densities could account for this factor of ∼2, as well, and there is evidence of a range of ratios, at least in the slow wind (Zhao et al. 2022).Slow wind is of course the more relevant, since it is more common at low latitudes.
In W14, the energy range was cut off at 5 MeV, so the comparison with the full reported data is not shown in that work.The baseline model of Wang et al. (2022) is an extension of W14 to higher energies.In Figure 4 of that work, they show a pure power law from 2 to 20 MeV, presumably close to the E −2.43 power law reported in Mewaldt et al. (2009), and Figure 3 of Li et al. (2023) shows something very similar.Our prediction drops well under this power law by around 10 MeV, but it agrees much better with the data at these energies.We believe the main cause of this discrepancy is the oxygen cross section used, specifically the extrapolation beyond 2 MeV using a power law from the points up to and near 2 MeV.We believe this choice of extrapolation to be directly leading to an overestimation of the fluence beyond ∼2 MeV, which grows with energy.
Yet, even at 1 MeV, near the peak of the oxygen cross section, comparison with W14 would indicate that either W14 are overestimating the downstream component or that we are underestimating it.We suspect that their form of the ENA transport equation is the issue, as it differs from the standard form or any seeming variant.In W14 it is Equation (4); in Wang et al. (2022) it does not appear, though they state that they use the same calculation as in W14; and in Li et al. (2023) it is Equation (A3), though in a different form.In none of these works is this equation derived, nor is any reference cited for it.We suspect that, since it is written as a three-dimensional integral, something could be wrong with how parcels from different areas of the sky are weighted, but this must remain speculation.Regardless, we have provided a rigorous proof of the transport equation, and thus we are quite confident of the estimate of the downstream distribution and of the results overall.

Conclusions
We have built upon recent work that has focused on the energetic neutral atoms (ENAs) from solar energetic particles (SEPs).The main conclusion that should be drawn from this work is that the upstream component of shocks is extremely important for the observation of ENAs from SEPs.The reason for the dominance of this component is that high-energy particles have a much higher chance of diffusion far upstream of the shock, and after becoming neutral there, have a much higher probability of escaping.We have also provided a rigorous derivation of the standard ENA transport equation, as well as an improved set of analytical functions for high-energy charge-exchange cross sections for use by the community.
subject to an external force F(r).For tractability and because neutral particles, out of the relevant forces, are only subject to gravity, which we may regard as negligible for energetic enough particles, we set F(r) to be zero.This assumption is not strictly necessary, and a similar calculation for lower energy neutrals has been carried out with the force due to the Sun's gravity and radiation pressure (Fahr 1979).
If we then choose a Cartesian coordinate system in which r = (x, y, z) and v = (v x , v y , v z ), Equation (A1) has the following semianalytical solution: A2 assuming appropriate boundary conditions, specifically that f → 0 as x → x 0 .As can be seen, Equation (A2) consists of a one-dimensional integral, containing a function, h, the survival probability, which itself contains another nested one-dimensional integral.There is a change of the time variable corresponding to velocity dispersion, as well as a change of variables for the location of the source (or loss) corresponding to a trace-back of the particles according to their velocity.
Given that differential flux, J, is related very simply to the distribution function, , Equation (A2) is already quite close to the form of the equation most often used (Gruntman et al. 2001;Prested et al. 2008).In these terms, r is the location of the "observer," and the "line of sight" is antiparallel to the unit vector defined by the velocity at the observer.Since the coordinate system chosen was completely arbitrary, any translation or rotation of the coordinate system is allowed; thus, the direction of integration can be defined to be along the line of sight, as in the standard form, but it is not necessary.
We can gain further insight into Equation (A2) (and motivate the form of Equation (A1)) by breaking up the production and loss into component terms: PI CX II where P CX is the production function due to charge exchange, as will be discussed below, and L PI , L CX , and L II are the loss rates due to photoionization, charge exchange, and impact ionization, respectively.Here we have kept only the production and loss terms that were used in W14; in general, they will all apply.
The simplest of the three is the photoionization loss rate, L PI , which is simply a velocity-independent ionization rate, sometimes referred to as β(r, t).
The next most complex is charge exchange, which we derive from the Boltzmann form for a binary collision reaction: where dΩ is the differential solid angle in the center-of-mass reference frame representing a relative change of direction of particles involved in a binary collision; s B ij is the differential interaction cross section between species i and species j; and v + and v − are the velocities after the collision of species i and j, respectively, which are in general a function of the incoming velocities, v and ¢ v , as well as the scattering angles, represented by Ω.For a charge-exchange reaction, the collision produces no change in direction.However, what does change are the identities of the particles.In a resonant reaction, e.g., between protons and hydrogen atoms, they literally exchange identities, i.e., i ↔ j; whereas, in the more general case, there is a mapping from incoming to outgoing particles, i.e.,  ¢ i i and  ¢ j j .Therefore, we must modify Equation (A6) slightly to account for the mapping, which we accomplish by splitting up the production and loss terms, while introducing more sums: where σ is equal to σ B integrated over the solid angle dΩ.The next step needed to recover the right-hand side of Equation (A1) is to assume that This assumption is appropriate for the case considered that shock-accelerated protons react with thermalized coronal and solar wind matter, with characteristic energies on the order of hundreds of eV to 1 keV, while in the analysis above, we consider only accelerated ions with at least ∼40 keV.With this assumption, we have where ( ) ¢ r n t , i and n j (r, t) are the number densities of the corresponding populations, and where we have assumed that the density of the parent population is approximately constant, i.e., that charge exchange induces a very small percent change on the parent population, allowing for the integrals over ¢ v to be taken.This form matches that in Equation (A1).
The last type of loss is impact ionization, which in general is a three-body interaction.However, for the purposes of this study, it can be treated as a two-body interaction because we are not interested in the distribution function of the reaction products, only the reactants (actually just the ENA reactant), and there are only two reactants in the impact ionization reactions.So, we can use the same formula derived from the Boltzmann equation for impact ionization.However, some of the other assumptions become questionable in the case of electron impact ionization, especially the proton speed (v) being much greater than the thermal electron speed ( ¢ v ).
A thermal electron in the 1 MK corona will have a characteristic speed of 5 Mm s −1 , which is approximately the same as a 100 keV proton.Thus, even our lower limit of ∼40 keV is apparently not enough to be clear of this issue.One could further raise the lowest energy considered here, but the usefulness of the study starts to diminish.We therefore accept that the lowest energies have some level of error, which will need to be evaluated and perhaps corrected for in a future study, as it is outside the scope of this one.

Appendix B Cross Sections
This section describes the cross sections used in this work.Figure 10 shows as curves the functions used for each reaction, as labeled in the figure, while the circular symbols represent tabularized data taken from the sources.The curves without associated symbols are those unchanged from their source reference, whereas the others have been interpreted or altered from the original in some way.
The unaltered cross sections are for the resonant protonhydrogen charge exchange (Lindsay & Stebbings 2005) and the electron impact ionization (Kim & Rudd 1994) reactions.For the charge-exchange reactions involving carbon and oxygen, we took data from Table 1 of Kuang (1992) and fit to those data a form of the lognormal distribution: where E is the energy of the interaction.
For the proton impact ionization reaction, we use the given Chebyshev fits to the cross section data at energies below 80 keV (Barnett et al. 1990).Since the Chebyshev fits diverge to infinity for higher energies, we fit a lognormal function to the cross-section data at energies above 80 keV using a nonlinear least squares minimization algorithm.
The best-fit parameters for the lognormal fits are in Table 1.Note that some of these cross sections, especially the oxygen, differ significantly from those in W14, Wang et al. (2022), andLi et al. (2023).Specifically, our cross sections drop off faster with energy than those in the mentioned works, leading to results that differ both quantitatively and qualitatively.
Figure 10.All cross sections used in this work.Fits to the oxygen and carbon charge-exchange cross sections of Kuang (1992) are shown as blue and red curves, respectively, while the purple curve is a combination of fits to the data in Barnett et al. (1990).See text in Appendix B for more details.

Figure 1 .
Figure1.The Monte Carlo simulation process is a computationally efficient framework that permits gradual improvement through more sophisticated specifications of plasma properties and scattering processes acting in the solar wind, which come from other simulation techniques and theories (dashed lines).

Figure 2 .
Figure 2. A sample particle trajectory from the Monte Carlo simulation used in this study, with different parameters.See Section 2.2 for details.

Figure 4 .
Figure 4. Distribution function, f, as a result of the shock simulation plotted (open circle symbols) as a function of the distance upstream of the shock, r − r s , normalized to the shock's location, r s .Colors represent energy.Also shown are fits of Equation (3) to the simulation results (curves).Simulation parameters are the same as in Figure 3.

Figure 5 .
Figure5.Best-fit parameters corresponding to the fits in Figure4(symbols) as a function of energy in keV; solid lines are guides to the eye and simply connect the points.Note that f 0 was not allowed to vary in the fit, only l and q; f 0 is the distribution function at the shock and comes directly from the simulation.The dashed line is the power law predicted by diffusive shock acceleration theory.

Figure 7 .
Figure 7. Angle-averaged differential flux as a function of time in hours.Shown are the upstream component (dashed), the downstream component (thin solid curves), and the sum (bold solid curves).Colors represent energy in keV.Shock parameters are the same as in Figure 3.

Figure 8 .
Figure8.Fluence of the angle-averaged flux in Figure7as a function of energy in keV.Shown are the upstream components (dashed curves) and downstream components (thin solid curves), as well as the contribution from oxygen (red), carbon (green), and hydrogen (blue); the sum of all components and charge exchange sources is also shown (bold black curve).Data (filled circular symbols) are adapted fromMewaldt et al. (2009).The power-law fit from that work is also shown (dotted line).Shock parameters are the same as in Figure3.

Figure 9 .
Figure 9. Fluence as a function of energy in keV.Various parameters of the shock model are shown.Data are the same as shown in Figure 8.Note the change in scale of the y-axis.

Table 1
Best-fit Parameters of the Lognormal Fit to the Tabular "Data"