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

The solar corona between below 10 solar radii is an important region for early acceleration and transport of solar energetic particles (SEPs) by coronal mass ejection-driven shock waves. There, these waves propagate into a highly variable dynamic medium with steep gradients and rapidly expanding coronal magnetic fields, which modulates the particle acceleration near the shock/wave surfaces, and the way SEPs spread into the heliosphere. We present a study modeling the acceleration of SEPs in global coronal shock events in the corona, as well as their transport to 1 au, based on telescopic observations coupled with dynamic physical models. As part of the project Solar Particle Radiation Environment Analysis and Forecasting—Acceleration and Scattering Transport (SPREAdFAST), we model the interaction of observed off-limb coronal bright fronts (CBF) with the coronal plasma from synoptic magnetohydrodynamic (MHD) simulations. We then simulate the SEP acceleration in analytical diffusive shock acceleration (DSA) model. The simulated fluxes are used as time-dependent inner boundary conditions for modeling the particle transport to 1 au. Resulting flux time series are compared with 1 au observations for validation. We summarize our findings and present implications for nowcasting SEP acceleration and heliospheric connectivity.

1 Introduction

One of the most common manifestations of solar activity are coronal mass ejections (CMEs). They are usually defined by observations in white light (Vourlidas et al., 2003; Zhang and Dere, 2006; Bein et al., 2011), but various aspects of these eruptions are observed in ultraviolet and radio wavelengths (Bastian et al., 2001; Veronig et al., 2010). In extreme ultraviolet (EUV) light, in particular, the early stages of CMEs may be observed well by telescopes such as the Atmospheric Imaging Assembly [(Lemen et al., 2012), AIA] on board the Solar Dynamics Observatory (Pesnell et al., 2012). CMEs may drive shock waves in the corona, if their propagation speeds exceed the local speed of information, usually the fast magnetosonic speed. Such shock waves are readily observed in EUV as so-called EUV waves (Thompson et al., 1998), also known as coronal bright fronts [(Long et al., 2011), CBFs].

CMEs are the largest contributor to the production of solar energetic particles (SEPs), ions and electrons of energies several orders of magnitude above the thermal coronal plasma (Reames, 1999). Flares are the other important source of SEPs in solar eruptions. SEP production in CMEs occurs mostly in the magnetized shock waves they drive, as well as in plasma compressions caused by the CMEs. Historically, the bulk of the SEP acceleration was thought to occur in interplanetary space, inferred from in situ observations of energetic storm particle (ESP) fluxes during the passage of interplanetary shocks by spacecraft. However, over the last fifteen years more advanced observations and numerical models have confirmed that in their early stages (below 5–10 R_Sun) CMEs often drive shocks (Ontiveros and Vourlidas, 2009; Gopalswamy and Yashiro, 2011), and those shocks may accelerate SEPs to energies up to and beyond 100 MeV/n (Battarbee et al., 2013; Kozarev et al., 2013; Schwadron et al., 2014; Kong et al., 2017). Consequently, recent work has been devoted to characterizing the dynamics of CMEs and the shocks they drive in the corona, using ever more advanced observations in white light, EUV, and radio.

This has been explored in studies to obtain estimates of the early-stage SEP production in the corona (Kozarev et al., 2013; Schwadron et al., 2015). In particular, Kozarev et al. (Kozarev et al., 2019) studied 9 separate western CBF events, and used the diffusive shock acceleration (DSA) model of Kozarev and Schwadron (Kozarev and Schwadron, 2016) to simulate the particle acceleration very early on, while the CMEs were still below 1.5 R_Sun. They discovered that the SEP production varies among events, and it also changes over the event’s duration. In addition, acceleration efficiency depends strongly on the varied coronal environment, through which the shock waves propagate. In this work, we extend the work of Kozarev et al. (Kozarev et al., 2019), by modeling the CBF-related shock/compression wave dynamics and particle acceleration out to 10 R_Sun, coupling the results to a global numerical particle transport model, and comparing the model results to in situ observations. This is a significant improvement to our method, and the first extended validated application of Sun-to-Earth physics-based modeling of SEP acceleration and transport, of which we are aware. We outline the methodology, which is the central part of the Solar Particle Radiation Environment Analysis and Forecasting—Acceleration and Scattering Transport (SPREAdFAST) framework for forecasting SEP events, in Section 2. The results of our study are presented in Section 3. We discuss them and present a summary in Section 4.

2 Materials and Methods

2.1 Event Selection

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 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¹.

TABLE 1

TABLE 1. Events in this study.

2.2 The SPREAdFAST Framework

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., 2022; 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.

2.3 CBF Kinematics and Geometric Modeling

To characterize the kinematics of CBFs, we applied to the AIA observations the methodology of the Coronal Analysis of SHocks and Waves (CASHeW) framework (Kozarev et al., 2017), updated and implemented in SPREAdFAST. Our framework estimates the CBF kinematics in a similar way to Long et al. (Long et al., 2021) and Downs et al. (Downs et al., 2021), 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) (Sheeley et al., 1999): 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 (Figure 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 (Figure 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.

FIGURE 1

FIGURE 1. Overview of the kinematics measurements for the May 11, 2011 event. (A) Global SDO/AIA 193 channel base difference image, showing clearly the CBF in mid-eruption. The radial measurement direction is shown as a blue solid line. The two lateral directions are shown as red/green dotted/dashed arcs. (B) The radial kinematic measurements, including the J-map and the CBF average front, peak, and back positions. (C) The northward kinematic measurements, including the J-map and the CBF average front, peak, and back positions. (D) The southward kinematic measurements and J-map.

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. (Kozarev et al., 2019). The smoothed positions were extrapolated to 10 Rs using the analytical model of Byrne et al. (Byrne et al., 2013).

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 Figure 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.

FIGURE 2

FIGURE 2. An example S2M geometric shock surface model, including the points belonging to the spheroid “cap” (in blue) and “zones” (green and red). The footpoints on the Sun are shown in dark red.

2.4 Coronal SEP Acceleration

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 (Kozarev and Schwadron, 2016; Kozarev et al., 2019) 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 |B| and shock angle θ_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 (Kozarev and Schwadron, 2016)) gives, for an initial momentum p₀, both the first distribution function (f₁) and momentum (p₁) values, and and an iterative solution for their values (f_i and p_i at subsequent time steps i, separated by the observational cadence δ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 (Torsti et al., 1995), observed during the 24 h 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.

2.5 Interplanetary SEP Transport

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 (Schwadron et al., 2010), 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_Sun. 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 (Kozarev et al., 2010; Schwadron et al., 2014).

For the EPREM model runs we performed on the 62 events here, we used the same input parameters, namely: mean free path λ₀ = 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⁻³; magnetic field magnitude at 1 au of |B| = 5.0, ×, 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 (Zank et al., 1998; Li et al., 2003; Sokolov et al., 2004; Verkhoglyadova et al., 2009):

${\lambda }_{\Vert }={\lambda }_0{\left(\frac{pc}{1GeV}\right)}^{1/3}{\left(\frac{R}{1AU}\right)}^{2/3},(1)$

giving the parallel mean free path for the simulation.

A 20-point energy grid was used with equal log-spacing between 1.0 and 200.0 MeV, as well as a 4-point pitch-angle grid. A simulation time-step of 0.5 au/c (approximately 4 min) 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 h from the onset of the event at the Sun for all events, since the goal is to model their initial stages.

3 Results

We have modeled 62 events with the SPREAdFAST model framework. To our knowledge, this physics-based Sun-to-Earth study is unprecedented in terms of the sheer number of events modeled from telescopic observations to in situ fluxes directly comparable to observations. All these events have related SEP flux increases at 1 au. Summary plots for each of the modeled events are available at the catalog web page of SPREAdFAST². Below we showcase the results for a single simulated event, and give summarized results for all 62 events.

3.1 Plasma Parameters Along Individual Shock-Crossing Magnetic Field Lines

For all events in the sample, the kinematics in the radial and lateral directions were obtained using the procedure described in Section 2.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.

FIGURE 3

FIGURE 3. Radial and lateral kinematics of the May 11, 2011 SPREAdFAST event. Shown are the radial (in blue) and lateral (in red, mean of the northward and southward directions) time-dependent front positions, as well as calculated speeds and accelerations. Also shown are the mean intensities (calculated as the mean value of the intensity in the J-maps between the SPREAdFAST-determined front and back of the CBFs) at each time step. The CBF thickness in each measurement and direction is also calculated from the automatically determined front and back of the CBF.

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-s 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 [(Mikić et al., 1999), 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³. 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 θ_BN, shock upstream magnetic field amplitude |B|, 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 (Cheung et al., 2015) 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 θ_BN, shock upstream magnetic field amplitude |B|, shock speed V_shock, and the shock density jump ratio, r. Similarly to other authors, (Frassati et al., 2019; Frassati et al., 2020; Long et al., 2021), 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 (Cheung et al., 2015), 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. (Kozarev et al., 2017). 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 (Kozarev et al., 2015; Vanninathan et al., 2015). Consequently, beyond the AIA FoV, where we do not have observational information, we set the value to 1.2.

FIGURE 4

FIGURE 4. Time series of the main parameters used for diffusive shock acceleration modeling, estimated at all shock-crossing field lines for the May 11, 2011 event. Note that the density jump ratio, r, is set to 1.2 (representing a weak shock) beyond the field of view of AIA. The different colors are only used to discern between the separate shock-crossing points, and do not carry any additional meaning in this Figure.

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 Figure 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 θ_BN and |B| panels shows a consistent decline in the values as a function of time (and radial distance).

FIGURE 5

FIGURE 5. A dynamic spectrum of the main parameters used for diffusive shock acceleration modeling, estimated at all shock-crossing field lines for the May 11, 2011 event. This is a summary of the result presented in Figure 4.

Perhaps the most important parameter for diffusive shock acceleration is the angle θ_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 θ_BN angle, with a significant decrease occurring in the first 50 min 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.

FIGURE 6

FIGURE 6. Dynamic spectra showing the distributions of the θ_BN shock-field angle as a function of time for the May 11, 2011 event. The (A) shows the distribution over time for the entire S2M surface. The (B–D) plots show the distributions for the S2M points of the cap, zone 1, and zone 2, shown in Figure 2 with blue, green, and red, respectively. The white dot symbols show the mean at each time step.

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 θ_BN time-dependent distribution is extracted separately for each region, as in Figure 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 θ_BN are seen in the cap/nose region, while in zone 2 the high values dominate throughout the event, with means consistently above 60°.

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 θ_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.

3.2 Diffusive Shock Acceleration of Protons

The suprathermal flux spectra derived from observations along the method described in Section 2.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 (Kozarev and Schwadron, 2016) 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.

FIGURE 7

FIGURE 7. The event input spectra for the DSA model, derived from 1 au quiet-time suprathermal observations.

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 min 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 Figure 2), and the consequent reduction in the number of shock-crossing field lines.

FIGURE 8

FIGURE 8. Differential proton fluxes simulated with the DSA model for the event, aggregated in 7 energy channels between 0.4 and 100 MeV. The sharp drop in the fluxes is due to the sudden change of the geometric model between a spheroid and a partial spheroid around 03:10.

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 ∼1000 lines for input to the DSA model.

FIGURE 9

FIGURE 9. Resulting fluence spectra for all shock-crossing field lines modeled with the SPREAdFAST DSA model for the May 11, 2011 event. The sum total of all fluences is also plotted with a light gray color, and the power law index to its fit is shown (in this case, δ = − 2.2.) We also show the input spectrum fluence with a dotted line below the colored output fluence lines. The plot also features reference power law slopes corresponding to shock density jumps r = 1.2 and r = 4.

3.3 Interplanetary SEP Transport Results and Comparisons to In Situ Observations

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.

FIGURE 10

FIGURE 10. The differential proton flux output of the EPREM model for the 110511 event near the L1 point. EPREM flux time series at the geometric mean energies of the ERNE channels are shown as colored solid lines. The ERNE observations are shown for comparison as dot-dashed lines colored the same as the corresponding EPREM energies. Time on the x-axis is in hours from the start of May 11, 2011.

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 ∼1 h 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.

FIGURE 11

FIGURE 11. Dynamic spectra of the differential proton fluxes of the ERNE observations (A) and the EPREM simulation (B) for the 110511 event. Time on the x-axis is in hours from the start of May 11, 2011.

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.

FIGURE 12

FIGURE 12. A comparison between the modeled (green) and observed (yellow) proton fluences for the event. The pre-event background ERNE spectrum is also shown with black symbols. Power law fit amplitudes and indices for all fluence spectra are shown as well.

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.

FIGURE 13

FIGURE 13. Scatter plots of the relationships between the fitted power index (A) for all 62 events, as well as the onset times (B) for EPREM model and ERNE observations.

On 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 (Laitinen et al., 2018; Laitinen and Dalla, 2019).

FIGURE 14

FIGURE 14. Scatter plots of the onset times at four separate energies for EPREM model and ERNE observations, similar to the right plot in Figure 14.

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.

FIGURE 15

FIGURE 15. A comparison between the modeled (yellow) and observed (blue) proton fluences for the event. Power law fit amplitudes and indices for both fluence spectra are shown as well. The pre-event background ERNE fluence, together with the power law fit is also shown for reference.

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 ∼ 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.

4 Conclusion

We have presented the first of its kind multi-event study of detailed physics-based Sun-to-1 au SEP simulations based on coronal diffusive shock acceleration and interplanetary propagation, of which we are aware. In this study, we investigated 62 separate eruptive events, which included an EUV CBF and sizable enhancements of 1 au proton fluxes. All but one had an associated flare. We applied the SPREAdFAST framework to model the evolution of the plasma immediately upstream of the coronal shock associated with the CBFs between 1 and 10 R_Sun, as well as the resulting acceleration and transport of protons between 1.05 R_Sun and 1 au. The SPREAdFAST framework, intended for forecasting the early-stages of SEP events, has proved very useful in the study of such a large set of events.

The input spectra for the coronal proton acceleration were derived from quiet-time suprathermal spectra averaged over one to three days preceding each event; the resulting averages were scaled back to the Sun assuming a simple inverse square proportionality to heliospheric distance, as in (Kozarev et al., 2019). We used an energy-independent mean free path for the coronal acceleration of the protons in the diffusive shock acceleration model of (Kozarev and Schwadron, 2016). Similarly to previous work (Kozarev et al., 2015; Rouillard et al., 2016; Kong et al., 2017; Kong et al., 2019), we found that the conditions in the solar corona influence significantly the acceleration of protons. Specifically, the large gradients in the plasma parameters between neighboring streamers, quiet-Sun areas, and coronal holes have the effect of continuous, smooth changes in the theta_BN angle along the shock wave surface, as well as in density and density enhancements. (Kong et al., 2019) found that the SEPs accelerated at a global shock wave expanding through a streamer in the low and middle corona will concentrate closer to the streamer apex, the higher the particle energy. However, the lower energy (10 MeV) protons were distributed more or less uniformly along the shock surface in their 2D model. Thus, answering the question “Where are SEPs produced?” in the early stages of CME-driven shock and compressive waves, depends strongly on the overlying coronal structure, and the particle energy.

The output from the DSA model was used as time-dependent input to the interplanetary transport EPREM model, which produced as input proton fluxes at 1 au. These were then compared with in situ observations by the SOHO/ERNE instrument. Overall, the results are very encouraging for the efficiency and accuracy of the SPREAdFAST model chain. However, we find that the fluxes at the highest energies show the most disagreement, due mostly to the slope of the increase and the onset times. Figure 15 shows an example of this for the May 11, 2011 event.

In future work, we will carry out more realistic modeling of the events in order to improve the match to observations. First, we will investigate how introducing time-dependent injection of the source spectra at the inner boundary of the EPREM simulation may change the L1 output, especially for those events in which the observed fluxes decay much faster than the modeled ones, or there are other discrepancies with observations. Although the focused transport is dominant, Dalla et al. (Dalla et al., 2020) show that 3D transport effects in a realistic interplanetary magnetic fields are quite important. In future work we will implement perpendicular transport to account for such effects. We will also include location-dependent output, which does not guarantee constant connectivity between the source and observer. Finally, constraining the geometric shock models with existing (such as SOHO/LASCO) and new observations of the CME evolution in the middle corona will help reduce uncertainty in the results. Comparing near-Sun in situ observations of the quiet-time suprathermal populations by Parker Solar Probe and Solar Orbiter, and comparing them with the 1 au fluxes, will help us improve the input spectra estimation.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://spreadfast.astro.bas.bg/catalog/, https://sdac.virtualsolar.org/cgi/search, https://catalogs.astro.bas.bg.

Author Contributions

KK lead the study and development of the SPREAdFAST framework; MN performed the event modeling runs; RM gathered the event observational information and contributed to the analysis; MD contributed to building the framework and to running the software; PZ contributed to the analysis.

Funding

The authors acknowledge support by ESA's PECS Programme for Bulgaria through contract 4000126128 to the Institute of Astronomy and NAO, Bulgarian Academy of Sciences, as well as support by the National Scientific Programme “VIHREN” in Bulgaria through contract \#KP-06-DV-8/18.12.2019 to the Institute of Astronomy and NAO, Bulgarian Academy of Sciences.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹www.lmsal.com/isolsearch

²http://spreadfast.astro.bas.bg/catalog

³http://www.predsci.com/mhdweb/

References

Bastian, T. S., Pick, M., Kerdraon, A., Maia, D., and Vourlidas, A. (2001). The Coronal Mass Ejection of 1998 April 20: Direct Imaging at Radio Wavelengths. apjl 558, L65–L69. doi:10.1086/323421

CrossRef Full Text | Google Scholar

Battarbee, M., Vainio, R., Laitinen, T., and Hietala, H. (2013). Injection of thermal and Suprathermal Seed Particles into Coronal Shocks of Varying Obliquity. aap 558, A110. doi:10.1051/0004-6361/201321348

CrossRef Full Text | Google Scholar

Bein, B. M., Berkebile-Stoiser, S., Veronig, A. M., Temmer, M., Muhr, N., Kienreich, I., et al. (2011). Impulsive Acceleration of Coronal Mass Ejections. I. Statistics and Coronal Mass Ejection Source Region Characteristics. apj 738, 191. doi:10.1088/0004-637X/738/2/191

CrossRef Full Text | Google Scholar

Byrne, J. P., Long, D. M., Gallagher, P. T., Bloomfield, D. S., Maloney, S. A., McAteer, R. T. J., et al. (2013). Improved Methods for Determining the Kinematics of Coronal Mass Ejections and Coronal Waves. aap 557, A96. doi:10.1051/0004-6361/201321223

CrossRef Full Text | Google Scholar

Cheung, M. C. M., Boerner, P., Schrijver, C. J., Testa, P., Chen, F., Peter, H., et al. (2015). Thermal Diagnostics with the Atmospheric Imaging Assembly on Board the Solar Dynamics Observatory: A Validated Method for Differential Emission Measure Inversions. apj 807, 143. doi:10.1088/0004-637X/807/2/143

CrossRef Full Text | Google Scholar

Dalla, S., de Nolfo, G. A., Bruno, A., Giacalone, J., Laitinen, T., Thomas, S., et al. (2020). 3D Propagation of Relativistic Solar Protons through Interplanetary Space. Astron. Astrophys. 639, A105. doi:10.1051/0004-6361/201937338

CrossRef Full Text | Google Scholar

Downs, C., Warmuth, A., Long, D. M., Bloomfield, D. S., Kwon, R. Y., Veronig, A. M., et al. (2021). Validation of Global EUV Wave MHD Simulations and Observational Techniques. Astrophys. J. 911, 118. doi:10.3847/1538-4357/abea78

CrossRef Full Text | Google Scholar

Frassati, F., Susino, R., Mancuso, S., and Bemporad, A. (2019). Comprehensive Analysis of the Formation of a Shock Wave Associated with a Coronal Mass Ejection. Astrophys. J. 871, 212. doi:10.3847/1538-4357/aaf9af

CrossRef Full Text | Google Scholar

Frassati, F., Mancuso, S., and Bemporad, A. (2020). Estimate of Plasma Temperatures across a CME-Driven Shock from a Comparison between EUV and Radio Data. Solar Phys. 295, 124. doi:10.1007/s11207-020-01686-0

CrossRef Full Text | Google Scholar

Gopalswamy, N., and Yashiro, S. (2011). The Strength and Radial Profile of the Coronal Magnetic Field from the Standoff Distance of a Coronal Mass Ejection-Driven Shock. apjl 736, L17. doi:10.1088/2041-8205/736/1/L17

CrossRef Full Text | Google Scholar

Kong, X., Guo, F., Giacalone, J., Li, H., and Chen, Y. (2017). The Acceleration of High-Energy Protons at Coronal Shocks: The Effect of Large-Scale Streamer-like Magnetic Field Structures. apj 851, 38. doi:10.3847/1538-4357/aa97d7

CrossRef Full Text | Google Scholar

Kong, X., Guo, F., Chen, Y., and Giacalone, J. (2019). The Acceleration of Energetic Particles at Coronal Shocks and Emergence of a Double Power-Law Feature in Particle Energy Spectra. apj 883, 49. doi:10.3847/1538-4357/ab3848

CrossRef Full Text | Google Scholar

Kozarev, K. A., and Schwadron, N. A. (2016). A Data-Driven Analytic Model for Proton Acceleration by Large-Scale Solar Coronal Shocks. apj 831, 120. doi:10.3847/0004-637X/831/2/120

CrossRef Full Text | Google Scholar

Kozarev, K., Schwadron, N. A., Dayeh, M. A., Townsend, L. W., Desai, M. I., and PourArsalan, M. (2010). Modeling the 2003 Halloween Events with EMMREM: Energetic Particles, Radial Gradients, and Coupling to MHD. Space Weather 8. doi:10.1029/2009SW000550 https://agupubs.onlinelibrary.wiley.com/action/showCitFormats?doi=10.1029%2F2009SW000550

CrossRef Full Text | Google Scholar

Kozarev, K. A., Evans, R. M., Schwadron, N. A., Dayeh, M. A., Opher, M., Korreck, K. E., et al. (2013). Global Numerical Modeling of Energetic Proton Acceleration in a Coronal Mass Ejection Traveling through the Solar Corona. apj 778, 43. doi:10.1088/0004-637X/778/1/43

CrossRef Full Text | Google Scholar

Kozarev, K. A., Raymond, J. C., Lobzin, V. V., and Hammer, M. (2015). Properties of a Coronal Shock Wave as a Driver of Early SEP Acceleration. apj 799, 167. doi:10.1088/0004-637X/799/2/167

CrossRef Full Text | Google Scholar

Kozarev, K. A., Davey, A., Kendrick, A., Hammer, M., and Keith, C. (2017). The Coronal Analysis of SHocks and Waves (CASHeW) Framework. J. Space Weather Space Clim. 7, A32. doi:10.1051/swsc/2017028

CrossRef Full Text | Google Scholar

Kozarev, K. A., Dayeh, M. A., and Farahat, A. (2019). Early-stage Solar Energetic Particle Acceleration by Coronal Mass Ejection-Driven Shocks with Realistic Seed Spectra. I. Low Corona. apj 871, 65. doi:10.3847/1538-4357/aaf1ce

CrossRef Full Text | Google Scholar

Laitinen, T., and Dalla, S. (2019). From Sun to Interplanetary Space: What Is the Pathlength of Solar Energetic Particles? Astrophys. J. 887, 222. doi:10.3847/1538-4357/ab54c7

CrossRef Full Text | Google Scholar

Laitinen, T., Effenberger, F., Kopp, A., and Dalla, S. (2018). The Effect of Turbulence Strength on Meandering Field Lines and Solar Energetic Particle Event Extents. J. Space Weather Space Clim. 8, A13. doi:10.1051/swsc/2018001

CrossRef Full Text | Google Scholar

Lemen, J. R., Title, A. M., Akin, D. J., Boerner, P. F., Chou, C., Drake, J. F., et al. (2012). The Atmospheric Imaging Assembly (AIA) on the Solar Dynamics Observatory (SDO). Solphys 275, 17–40. doi:10.1007/s11207-011-9776-8

CrossRef Full Text | Google Scholar

Li, G., Zank, G. P., and Rice, W. K. M. (2003). Energetic Particle Acceleration and Transport at Coronal Mass Ejection-Driven Shocks. J. Geophys. Res. (Space Physics) 108, 1082. doi:10.1029/2002JA009666

CrossRef Full Text | Google Scholar

Long, D. M., DeLuca, E. E., and Gallagher, P. T. (2011). The Wave Properties of Coronal Bright Fronts Observed Using SDO/AIA. apjl 741, L21. doi:10.1088/2041-8205/741/1/L21

CrossRef Full Text | Google Scholar

Long, D. M., Reid, H. A. S., Valori, G., and O’Kane, J. (2021). Localized Acceleration of Energetic Particles by a Weak Shock in the Solar Corona. Astrophys. J. 921, 61. doi:10.3847/1538-4357/ac1cdf

CrossRef Full Text | Google Scholar

Mikić, Z., Linker, J. A., Schnack, D. D., Lionello, R., and Tarditi, A. (1999). Magnetohydrodynamic Modeling of the Global Solar corona. Phys. Plasmas 6, 2217–2224. doi:10.1063/1.873474

CrossRef Full Text | Google Scholar

Ontiveros, V., and Vourlidas, A. (2009). Quantitative Measurements of Coronal Mass Ejection-Driven Shocks from LASCO Observations. apj 693, 267–275. doi:10.1088/0004-637X/693/1/267

CrossRef Full Text | Google Scholar

Pesnell, W. D., Thompson, B. J., and Chamberlin, P. C. (2012). The Solar Dynamics Observatory (SDO). Solphys 275, 3–15. doi:10.1007/s11207-011-9841-3

CrossRef Full Text | Google Scholar

Reames, D. V. (1999). Particle Acceleration at the Sun and in the Heliosphere. Space Sci. Rev. 90, 413–491. doi:10.1023/A:1005105831781

CrossRef Full Text | Google Scholar

Rouillard, A. P., Plotnikov, I., Pinto, R. F., Tirole, M., Lavarra, M., Zucca, P., et al. (2016). Deriving the Properties of Coronal Pressure Fronts in 3D: Application to the 2012 May 17 Ground Level Enhancement. apj 833, 45. doi:10.3847/1538-4357/833/1/45

CrossRef Full Text | Google Scholar

Schwadron, N. A., Townsend, L., Kozarev, K., Dayeh, M. A., Cucinotta, F., Desai, M., et al. (2010). Earth-Moon-Mars Radiation Environment Module Framework. Space Weather 8. doi:10.1029/2009SW000523https://agupubs.onlinelibrary.wiley.com/action/showCitFormats?doi=10.1029%2F2009SW000523

CrossRef Full Text | Google Scholar

Schwadron, N. A., Gorby, M., Török, T., Downs, C., Linker, J., Lionello, R., et al. (2014). Synthesis of 3-D Coronal-Solar Wind Energetic Particle Acceleration Modules. Space Weather 12, 323–328. doi:10.1002/2014SW001086

CrossRef Full Text | Google Scholar

Schwadron, N. A., Lee, M. A., Gorby, M., Lugaz, N., Spence, H. E., Desai, M., et al. (2015). Particle Acceleration at Low Coronal Compression Regions and Shocks. apj 810, 97. doi:10.1088/0004-637X/810/2/97

CrossRef Full Text | Google Scholar

Sheeley, N. R., Walters, J. H., Wang, Y. M., and Howard, R. A. (1999). Continuous Tracking of Coronal Outflows: Two Kinds of Coronal Mass Ejections. jgr 104, 24739–24768. doi:10.1029/1999JA900308

CrossRef Full Text | Google Scholar

Sokolov, I. V., Roussev, , Gombosi, T. I., Lee, M. A., Kóta, J., Forbes, T. G., et al. (2004). A New Field Line Advection Model for Solar Particle Acceleration. apjl 616, L171–L174. doi:10.1086/426812

CrossRef Full Text | Google Scholar

Thompson, B. J., Plunkett, S. P., Gurman, J. B., Newmark, J. S., St Cyr, O. C., and Michels, D. J. (1998). SOHO/EIT Observations of an Earth-Directed Coronal Mass Ejection on May 12, 1997. grl 25, 2465–2468. doi:10.1029/98GL50429

CrossRef Full Text | Google Scholar

Torsti, J., Valtonen, E., Lumme, M., Peltonen, P., Eronen, T., Louhola, M., et al. (1995). Energetic Particle Experiment ERNE. solphys 162, 505–531. doi:10.1007/BF00733438

CrossRef Full Text | Google Scholar

Vanninathan, K., Veronig, A. M., Dissauer, K., Madjarska, M. S., Hannah, I. G., and Kontar, E. P. (2015). Coronal Response to an EUV Wave from DEM Analysis. apj 812, 173. doi:10.1088/0004-637X/812/2/173

CrossRef Full Text | Google Scholar

Verkhoglyadova, O. P., Li, G., Zank, G. P., Hu, Q., and Mewaldt, R. A. (2009). Using the Path Code for Modeling Gradual SEP Events in the Inner Heliosphere. apj 693, 894–900. doi:10.1088/0004-637X/693/1/894

CrossRef Full Text | Google Scholar

Veronig, A. M., Muhr, N., Kienreich, I. W., Temmer, M., and Vršnak, B. (2010). First Observations of a Dome-Shaped Large-Scale Coronal Extreme-Ultraviolet Wave. apjl 716, L57–L62. doi:10.1088/2041-8205/716/1/L57

CrossRef Full Text | Google Scholar

Vourlidas, A., Wu, S. T., Wang, A. H., Subramanian, P., and Howard, R. A. (2003). Direct Detection of a Coronal Mass Ejection-Associated Shock in Large Angle and Spectrometric Coronagraph Experiment White‐Light Images. apj 598, 1392–1402. doi:10.1086/379098

CrossRef Full Text | Google Scholar

Zank, G. P., Matthaeus, W. H., Bieber, J. W., and Moraal, H. (1998). The Radial and Latitudinal Dependence of the Cosmic ray Diffusion Tensor in the Heliosphere. J. Geophys. Res. Space Phys. 103, 2085. doi:10.1029/97JA03013

CrossRef Full Text | Google Scholar

Zhang, J., and Dere, K. P. (2006). A Statistical Study of Main and Residual Accelerations of Coronal Mass Ejections. apj 649, 1100–1109. doi:10.1086/506903

CrossRef Full Text | Google Scholar