A Multi-Event Study of Early-Stage SEP Acceleration by CME-Driven Shocks - Sun to 1 AU

We selected the events on the basis of pre-identified proton events in 17-22 MeV observed by the SOHO/ERNE instrument in the period 2010–2017, obtaining 216 events. Proton events without identified flares and CMEs and without EUV waves (prior the SEP event) were first dropped from the analysis leading to 177 events. Cases without EUV waves or no EUV data, even though they may have identified flare/CME, were also dropped from any further analysis, as the requirement of the SPREAdFAST model is the presence of an EUV wave, leading to 105 events. Several uncertain EUV waves were dropped due to their relevance to a different solar eruption, leading to 99 events. Of those, a final selection of 62 events had measurable near-limb or off-limb CBFs, which can be analyzed with our framework. All events are given in Table LABEL:table:1. Since there are multiple aspects and phenomena related to the events, we provide as reference the date; start, end times, and class of the associated flare; and the source location on the disk in helioprojective Cartesian coordinates. These were obtained from the Heliophysics Events Knowledge Base ¹¹1www.lmsal.com/isolsearch.

To model the particle fluxes at 1 au and compare them with observations, we utilize the physics-based global modeling framework Solar Particle Radiation Environment Analysis and Forecasting–Acceleration and Scattering Transport (SPREAdFAST; Kozarev et al. 2021, in preparation). The framework combines detailed EUV observations of CBF events with modeling of the interacting coronal plasma and the resulting SEP production and interplanetary transport. We briefly outline the various components of the SPREAdFAST framework used here.

To characterize the kinematics of CBFs, we applied to the AIA observations the methodology of the Coronal Analysis of SHocks and Waves (CASHeW) framework [20], updated and implemented in SPREAdFAST. Our framework estimates the CBF kinematics in a similar way to Long et al. [21] and Downs et al. [22], by following the leading edge of the front on consecutive images. In addition, we calculate the kinematics of the peak and back edge of the CBFs over time, which allows us to estimate their time-dependent mean intensity and thickness. Figure 1 shows an example of its application to the May 11, 2011 CBF event. In the top left panel (A) is a mid-event AIA 193-channel base difference image, clearly showing the CBF. Overlaid are the radial and lateral lines of measurement. The kinematics are determined using time-height maps (also known as J-maps)[23]: images created by stacking horizontally in time columns of pixels in a desired direction from a solar image. The shape of a track on these images depends on the direction and speed of the feature measured. The system identifies the radial and lateral wave front positions over time in the J-maps generated with the CASHeW code for each event. In the radial direction (Fig. 1, panel B), they are measured along the line passing through the solar center and the CBF nose (predominant direction of motion). In the lateral direction (Fig. 1, panels C and D), parallel to the solar limb, the wave front signatures are measured in two directions away from the radial direction. In some cases, lateral wave signatures are only visible in one direction, or are missing altogether.

Summary plots of the J-maps with overlaid estimated positions and errors for all events are available at the SPREAdFAST catalog. The lateral measurements in both directions are averaged to form a single lateral kinematics time series for each event. The CBF width and mean intensity in both directions are also recorded. Based on the kinematic measurements, a Savitzky-Golay fit is performed to obtain the kinematics in the AIA FOV, as in Kozarev et al. [18]. The smoothed positions were extrapolated to 10 Rs using the analytical model of Byrne et al. [24].

Based on the radial and lateral measured front positions over time, a three-dimensional, time-dependent geometric spheroid model - Synthetic Shock Model (S2M) - was developed for all compressive fronts, consisting of a large number of points ($>$1000), used for the estimation of the dynamic shock upstream coronal parameters. Figure 2 demonstrates 9 time steps of the S2M for an eruption near the northwest limb of the Sun. The spheroid remains centered on the eruption source throughout the event, while its aspect ratio varies based on the radial and lateral CBF position measurements and extrapolations. Once the S2M wraps around the Sun, all points behind the plane passing through the eruption source point, and perpendicular to the radial direction in that point, are removed from consideration. The blue, green, and red points of the S2M in Fig. 2 separate a cap and two side zones, used to investigate the differences in the plasma environment through which the CBF passes in the radial and lateral directions. We have defined the ‘nose’ of the shock model as the collection of all model points on the spheroidal cap subtending a half-angle of Pi/7 from the semi-major axis of the spheroid at each time step. The points on the remaining part of the shock surface is divided into two flanks, or zones, with the dividing plane always parallel to the Sun-Earth line. The plasma parameters at the points on these three surfaces can be examined separately, as demonstrated in Figure 2.

After the plasma parameters along the individual shock-crossing field lines have been established, they are fed in a time-dependent manner into the coronal DSA model [19, 18] for calculation of the proton acceleration between the low corona and $10~{}R_{Sun}$. The model was specifically developed to take as input remote solar observations and data-driven model output from the CASHeW framework. The model solves for the coronal charged particle acceleration by large-scale CME-driven shocks. Our model uses time-dependent estimates of shock speed $V_{shock}$, density jump ratio $r$, magnetic field strength $\lvert B\rvert$ and shock angle $\theta_{BN}$, for multiple shock-crossing field lines. The model calculates the minimum shock injection momenta for the particles. It takes as input a particle distribution function, and provides time-dependent distribution function spectra or fluxes, as output. The solution obtained (Eqns. 8-11 in [19]) gives, for an initial momentum $p_{0}$, both the first distribution function ($f_{1}$) and momentum ($p_{1}$) values, and and an iterative solution for their values ($f_{i}$ and $p_{i}$ at subsequent time steps $i$, separated by the observational cadence $\delta t$ of the instrument (SDO/AIA, in this case). The model is run for each individual shock-crossing field line, based on observed and calculated parameters at a single shock-crossing point along it. Flux spectra at each time step are then computed. It has been validated and used for the analysis of a number of SEP events.

We use as input to the model observations-based suprathermal proton spectra derived from 1 au fluxes with the SOHO/ERNE instrument [25], observed during the 24 hours of quiet time preceding each SEP event. We fit power laws to each suprathermal spectrum in the energy range $0.056-3.0~{}MeV$, and scale them to $1.05~{}R_{Sun}$, assuming a simple inverse square dependence on radial distance (implying flux conservation). We currently do not consider adiabatic cooling or other particle transport effects; these are indeed important, but a full study of these effects on suprathermals must be performed in a future study, in order to determine a general trend to use for forecasting. This produces time-independent power law input spectra for the DSA model. The suprathermal spectrum calculated for $1.05~{}R_{Sun}$ is injected at all shock positions and distances - in the current implementation it is not modified to account for the changing shock locations.

The final step of the modeling chain is the transport of the accelerated SEPs to 1 au, and subsequent comparison with particle observations with the ERNE instrument. This is achieved by taking the resulting averaged fluxes for the entire event (here an example with the event on May 11, 2011 is shown) as input to the modified version of the EPREM model [26], and transporting them through a Parker-type static interplanetary medium. The particle injection from the DSA model into ERPEM is continued over the duration of the coronal shock event. The model includes the effects of pitch-angle scattering, adiabatic focusing and cooling, convection and streaming, and stochastic acceleration. The solver requires inner boundary conditions, with no initial conditions imposed. It features a dynamical simulation grid, in which the computational nodes are carried away from the Sun (frozen-in) with the solar wind - thus the connected grid nodes (streamlines) naturally assume the shape of a three-dimensional interplanetary magnetic field, along which energetic particles propagate. By default, EPREM incorporates an interplanetary magnetic field model with a radial field component falling off as the inverse square of radial distance, azimuthal component falling off as the inverse of radial distance, and a constant latitudinal component - the so-called Parker spiral model. The spatial grid is housed in a data structure based on nested cubes, whose surfaces are regularly subdivided into square arrays of square cells. These cells represent the structure of the grid, within which computational nodes propagate. The inner boundary surface rotates with the solar rotation rate, and is expelled outwards at the solar wind speed. At each time step, a new shell of cells is created at the inner boundary of the grid, and starts its propagation outwards. The inner boundary for the EPREM simulation is at 1.05$R_{\odot}$. The outer boundary varies for the individual field lines due to the varying dynamic conditions, but it always exceeds 1 AU. The model has been extensively validated and used for SEP studies [27, 15].

For the EPREM model runs we performed on the 62 events here, we used the same input parameters, namely: mean free path $\lambda_{0}=0.1$ au at 1 au and 1 GV magnetic rigidity; constant solar wind $V_{sw}=500$ km/s; proton number density at 1 au $n=5.0$ cm^-3; magnetic field magnitude at 1 au of $\lvert B\rvert=5.0\times 10^{-5}$ G. The mean free path is additionally scaled with proton rigidity and radial distance from the Sun to reflect the magnetic turbulence spectrum and its radial dependence [28, 29, 30, 31]:

giving the parallel mean free path for the simulation.

A 20-point energy grid was used with equal log-spacing between 1.0 MeV and 200.0 MeV, as well as a 4-point pitch-angle grid. A simulation time-step of 0.5 au/c (approximately 4 minutes) was used with 30 sub-steps, allowing for the proper calculation of the SEP propagation among nodes. All effects of diffusive transport were included in these baseline simulations (adiabatic cooling/heating, adiabatic focusing, pitch-angle scattering, convection with the solar wind, streaming). The effects of perpendicular diffusion and particle drifts will be included in subsequent work. The simulations were stopped at 9.6 hours from the onset of the event at the Sun for all events, since the goal is to model their initial stages.

For all events in the sample, the kinematics in the radial and lateral directions were obtained using the procedure described in Sect. 3.3. As an example, Figure 3 shows a summary plot of the kinematics characterization for a single event, which occurred on 11 May, 2011 (denoted as 110501_01). The top left panel shows the radial (blue) and lateral (red) positions above the solar surface, while the other panels show the derived speeds, accelerations, mean wave intensities, and wave thickness (both calculated by using the SPREAdFAST-determined CBF front and back positions in the J-maps). The middle-row right plot shows the extrapolations of the major and minor axes in orange and blue, while the aspect ratio evolution is shown in black. Uncertainties are shown as light blue and light red shadings.

Based on the extrapolated radial and lateral kinematics results obtained for each event, and the inferred major and minor axes of the spheroid representing the shock wave, a three-dimensional geometric model is created, which describes the shock surface at 24-second intervals from the onset of the CBF to the time when its nose (leading direction’s front position) reaches 10 $R_{Sun}$. This model shock surface is then propagated through the solar corona up to 10 $R_{Sun}$, represented by the domain of the results of a Magnetohydrodynamic Algorithm outside a Sphere [32, MAS] synoptic coronal magnetohydrodynamic (MHD) model for the Carrington rotation of each event. The MAS coronal model results are freely available on Predictive Science Inc’s MHDweb service³³3http://www.predsci.com/mhdweb/. The shock surface thus samples at discrete points the relevant parameters for coronal shock acceleration of SEPs, determined by the magnetic field lines crossing it consecutively at each point.

We sample on the order of 1000 shock-crossing field lines for each modeled event, obtaining time-dependent estimates for the shock-field angle $\theta_{BN}$, shock upstream magnetic field amplitude $\lvert B\rvert$, shock speed $V_{shock}$, plasma density $n$, and Alfven speed $V_{A}$. The shock density jump, $r$, is determined by estimating for the emission measure and the average density at the shock crossing points at each time step, and dividing by a pre-event average density model. The emission measure and the average densities within the field of view of the AIA instrument are estimated by applying a differential emission measure (DEM) model [33] to each group of pixels surrounding the projected shock-field crossing points for every timestep.

Figure 4 shows parameter evolution along the individual shock-crossing field lines, for the May 11, 2011 event, between the solar surface and 10 $R_{Sun}$. The panels show, in different colors, the shock-field angle $\theta_{BN}$, shock upstream magnetic field amplitude $\lvert B\rvert$, shock speed $V_{shock}$, and the shock density jump ratio, $r$. Similarly to other authors, [34, 35, 21], we employ DEM analysis to estimate the shock strength from AIA observations. We find that, while we use a different DEM model than these authors [33], the general results are similar - we find weak coronal shocks. The advantage of our method is that we estimate the density change at all points of the shock model, not only one or several regions. This method has previously been employed in Kozarev et al. [20]. We find that the density jump based on the DEM modeling within the AIA field of view is generally quite small, predominantly below 1.2, which is consistent with many previous studies [36, 37]. Consequently, beyond the AIA FoV, where we do not have observational information, we set the value to 1.2.

A better cumulative view of the evolution of these parameters for this event is shown in Figure 5. It shows the same four parameters from Fig. 4 as dynamic spectra over all shock-crossing field lines, with the colors scaled to the minimum/maximum percentages of shock-crossing points in the corresponding parameter value bins for each timestep. This representation shows better the instantaneous distribution of the DSA-relevant plasma parameters at the shock front. For example, the $\theta_{BN}$ and $\lvert B\rvert$ panels shows a consistent decline in the values as a function of time (and radial distance).

Perhaps the most important parameter for diffusive shock acceleration is the angle $\theta_{BN}$. A clearer picture of its evolution in the model corona for this particular event emerges when the values at each timestep are aggregated into histograms, and shown as a dynamic spectrum, as in Figure 6. The top left plot shows the distribution as a function of time for the entire spheroid surface. The color table represents the percentage of the total occurrences in each bin at each timestep. The overall evolution for this event is from high to low $\theta_{BN}$ angle, with a significant decrease occurring in the first 50 minutes of the event. The detailed dynamic spectra of the shock plasma parameters, as well as all results shown for individual events below, are available on the SPREAdFAST catalog web page.

A more nuanced picture may be observed if the full surface of the S2M spheroid is logically divided into points belonging to a cap (nose) and two zones (flanks) to the north and south - and the $\theta_{BN}$ time-dependent distribution is extracted separately for each region, as in Fig. 2. It shows the distribution of points in separate snapshots of the geometric model, with the cap in blue color, zone 1 in green, and zone 2 in red. These correspond to the top right and bottom plots in Figure 6. The lowest values of $\theta_{BN}$ are seen in the cap/nose region, while in zone 2 the high values dominate throughout the event, with means consistently above $60^{\circ}$.

We note that it is currently not possible for our model framework to account for the changes occurring in the corona between two closely separated in time events (i.e. two events originating from the same active region in two consecutive days), because it relies on a synoptic time-independent MHD model of the corona. The difference between such historical events with respect to modeling is mainly in the kinematics and shape of the compressive front, which however is expressed in differences in the shock-crossed field lines, shock-field angle $\theta_{BN}$, etc. This naturally has an effect on the modeled fluxes. The observed fluxes from such events also differ. In future work we will strive to incorporate daily updated coronal plasma models, in order to relax the assumption of no short-term change in the coronal conditions.

The suprathermal flux spectra derived from observations along the method described in Section 3.4 for all events are plotted in Figure 7. It shows a variety of power law indices, dominated by softest spectra. The corresponding time-constant flux spectrum for each event is extended to the range 0.056-200 MeV, and is fed as initial condition into the DSA model of Kozarev and Schwadron [19] throughout its spatial domain (1.05-10 $R_{Sun}$). The model’s energy grid contains 30 energy steps in this range. For this study, we set a constant upstream scattering mean free path in the corona mfp=0.0055 $R_{Sun}$. The DSA model produces time-dependent distribution function histories for all shock-crossing field lines, which are then aggregated to calculate the average flux time series for each event as a function of energy.

The coronal DSA modeling results for the event are shown on Figure 8. The plot shows flux time series at seven energy ranges between 0.4 and 100 MeV, obtained by summing the resulting fluxes over all shock-crossing points. The simulation covered over three hours of the event, until the shock reached 10 $R_{Sun}$. The time series show clearly the time-varying (mostly increasing) fluxes, with the first 90 minutes being most dynamic. The sudden drop in flux around 3:10 is due to the shift from the fully-connected spheroid to the reduced surface shock model (i.e. from step 5 to 6 in Fig. 2), and the consequent reduction in the number of shock-crossing field lines.

The integrated spectra of all individual modeled shock-crossing field lines for event are shown in Figure 9. The fluences were obtained by integrating over the time period the individual line fluxes. The energy extent is 0.056-200 MeV, and the overall combined slope is -2.2. The sum of all fluences is also plotted in light gray color. This figure illustrates the variety of accelerated proton fluences obtained due to the different coronal conditions in the simulation. In the interplanetary simulation, the average of these fluences is used. The plots show significant acceleration to even the highest energies for most lines. Realistic fluences in acceleration under time-changing shock plasma parameters do not have to be power laws, that is why some fluence lines have more complex shapes. We note that the S2M model provides upwards of 10 000 shock-crossing field lines; for manageability reasons the framework randomly selects a subset of $\sim$1000 lines for input to the DSA model.

Figure 10 shows the output of the EPREM model for the event. Proton differential fluxes as a function of time in hours are shown for the geometric mean energies of all ERNE channels between 3 and 115 MeV. The ERNE observations are also shown for direct comparison with dot-dashed lines, and corresponding colors for each energy channel. In order to improve the calculations of onset times and fluences, we add to the energetic particle model output the average pre-event flux from the ERNE observations. It is calculated by averaging the first five time steps of the ERNE fluxes in each energy channel.

Overall, for this event the model yields good agreement for the flux levels at most energies, compared with the observations. At the low-energy channels the agreement is best, including the onset of the particle event. Above 15 MeV, there is a discrepancy in the time profile, with the observed proton fluxes rising $\sim$1 hr before the simulation. This can also be seen in a more detailed comparison in Figure 11. It shows dynamic spectra of the ERNE and EPREM flux intensities. The fluxes in each energy channel have been normalized to the maximum ERNE flux values in it. The pre-event ERNE flux background is also included, for direct comparison to observations. Such discrepancies will be addressed in future work, and likely corrected with improvements to the SPREAdFAST model chain. However, it is a testament to the strength of our method that such good overall agreement has been achieved with an initial set of input parameters, given the extremely complex nature of the problem and the many unknown parameters.

In Figure 12 we show a comparison between the event-integrated fluence spectra from the EPREM model and ERNE observations for the same event. The agreement is somewhat close, in terms of low-energy amplitude, and in terms of the fitted power law. In terms of radiation hardness, the EPREM spectrum is slightly flatter, which is desirable when making forecasts. In other words, in general it is better to forecast a harder spectrum than the observed, rather than a softer spectrum. This comparison between modeled and observed spectra is one of the metrics we discuss next as a way to evaluate the model performance over the full set of 62 events.

We have further explored the relations between the modeled and observed SEP parameters at 1 au for the full set of 62 studied events. Figure 13 shows scatter plots of the relationships between the fitted power indices of the proton fluences from the EPREM model and the ERNE observations (left panel), and the onset hours for the two sets (right panel), defined as the time of 10% increase above the average pre-event flux at the energy channel centered on 28 MeV. This energy is high enough that it is above the geometric mean of the ERNE energy range, but also low enough that flux enhancements are observed at this energy for most events. On the left panel, although the power law indices are centered around -2.5 for both EPREM and ERNE, there is considerable scatter, with a number of outliers, which we will investigate in future work. It is likely due to the input spectral shape and the scattering mean free path.

In the right panel, we find a quite strong linear relation between the onset times at 28 MeV, which gives us confidence in the modeling framework’s performance. We further demonstrate the onset relationships for four higher energies in Figure 14. In future work, we will focus on narrowing the deviation of onset times from the observations by modifications in the input SEP spectra and the transport parameters. We note also that the EPREM model currently uses a nominal Parker spiral model for the interplanetary magnetic field; a realistic interplanetary model would likely change the length of the large-scale magnetic field lines along which the protons travel, due to meandering effects [38, 39].

We have also computed and show in Figure 15: histograms of the difference in hours between the modeled and observed flux onset at 28 MeV; the modeled (EPREM_powindex) and observed (ERNE_powindex) index of the power law fits to the resulting fluence spectra; and the Mean Squared Logarithmic Error (MSLE) for the fluence spectra and the event onsets. One of the most important parameters is the onset time - it is desired that the modeling framework predicts the proton fluxes arriving at 1 au before rather than after the actual onset, if it is to be used for forecasting. In this first study, we find that for about two-thirds of the events the modeled onset is before or close to the observed time (positive hour values in the top-left panel). The rest of the events will be further investigated in future work.

Comparing the EPREM_powindex and ERNE_powindex, we see an overall agreement of the peaks and extent of the distributions, although the simulations must be improved still. The last two panels show that the spectral and onset errors are sharply peaked near 0, meaning that for most events the resulting spectra are very similar, and the onsets are very close. This is a very positive result, demonstrating the capability of the framework. The individual flux and fluence comparisons are available at the SPREAdFAST catalog web page for all 62 events.

We find that in 2/3 of the events there is good agreement in the fluxes between model and observations. However, discrepancies were found in a number of events. This points to a likely existence of additional seed population beyond the quiet-time suprathermal spectra, possibly from flare-accelerated protons. We further find significant agreement with the observed proton flux profiles in the early stages for energies below 12 MeV. However, the model often produces higher fluxes than observations at the energies above $\sim 30$ MeV, suggesting an energy dependence of the scattering mean free path, which must be investigated.

For most of the events, the ERNE fluxes continue rising for the full duration of the interplanetary simulation after onset, indicating continued injection of particles near the Sun, in line with the model setup. For some events, the flux increases are short-lived, resembling impulsive injections convolved with the transport effects. For other events still, there is a relatively quick increase and drop, followed by more prolonged increase. This may be due to flux dropout effects, and will be further investigated in future work.