EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 659 (March 2022) A&A, 659 (2022) A187 Full HTML
Free Access
Issue
A&A
Volume 659, March 2022
Article Number A187
Number of page(s) 14
Section The Sun and the Heliosphere
DOI https://doi.org/10.1051/0004-6361/202142698
Published online 25 March 2022

© ESO 2022

1. Introduction

Occasionally, particle detectors on board spacecraft measure strong intensity enhancements when an interplanetary shock driven by a coronal mass ejection (CME) crosses the spacecraft. Such events are referred to as energetic storm particle (ESP) events (Bryant et al. 1962; Gosling et al. 1981) and are believed to be the result of continuous particle acceleration and trapping at the CME shock. A viable acceleration mechanism is diffusive shock acceleration (DSA), which involves a self-generated turbulence process (e.g. Bell 1978; Lee 1983; Vainio & Laitinen 2007; Ng & Reames 2008; Afanasiev et al. 2015). This turbulence may create a foreshock region that traps particles near the shock wave, allowing an efficient DSA process and thereby explaining the sudden particle intensity increases observed near CME-driven shocks.

Not every CME-driven shock produces detectable energetic particle intensity enhancements (e.g. Lario et al. 2003), and it is not yet fully understood why this is the case. A plethora of factors are believed to contribute to the (in)efficiency of a shock as a particle accelerator, such as its geometry and speed, the turbulence conditions near the shock wave, the presence of a suprathermal seed population, and the preconditioning produced by preceding CMEs (see e.g. Desai & Giacalone 2016; Reames 2017; Guo et al. 2021, for recent reviews). Many of these factors depend critically on the properties of the plasma environment through which the CME-driven shock propagates. In addition, the plasma conditions far from the shock wave can strongly influence the transport of the energetic particles that have escaped the shock acceleration site. Hence, to model a specific solar energetic particle (SEP) event, it is desirable to have an adequate description of the background plasma. To some extent, this can be achieved by using a global magnetohydrodynamic (MHD) model. An energetic particle transport model can then be used to study how the modelled CME shock affects the transport and acceleration of particles. Such an approach was, for example, taken by Kozarev et al. (2013), Schwadron et al. (2015), and Young et al. (2021) to model particle acceleration in the low corona using the energetic particle radiation environment module (EPREM) model (Schwadron et al. 2010, 2014) together with the Alfvén wave solar model (AWSoM; van der Holst et al. 2010; Manchester et al. 2012). Similarly, Li et al. (2021) recently used the improved particle acceleration and transport (iPATH; Hu et al. 2017) model together with the AWSoM code to model the SEP event observed on 17 May 2012.

In this work we follow a similar approach as the aforementioned studies, but we focus on modelling the interplanetary acceleration and transport of low-energy protons (50 keV−2 MeV) associated with the passage by 1 au of an interplanetary CME-driven shock. This is done by using the particle transport code named ‘PArticle Radiation Asset Directed at Interplanetary Space Exploration’ (PARADISE; Wijsen 2020), which is coupled to the data-driven inner-heliospheric model EUHFORIA (Pomoell & Poedts 2018). EUHFORIA stands for ‘EUropean Heliospheric FORcasting Information Asset’ and simulates the propagation of the solar wind and CMEs through interplanetary space by solving the ideal MHD equations using boundary conditions derived from a solar magnetogram in a semi-empirical manner. In Wijsen et al. (2021), the EUHFORIA + PARADISE model was used to successfully reproduce an energetic particle event associated with a corotating interaction region (CIR), which was observed in September 2019 by both Parker Solar Probe (PSP) and the Ahead spacecraft of the Solar TErrestrial RElations Observatory (STEREO-A).

In this work we use the EUHFORIA + PARADISE model to simulate the interplanetary acceleration of low-energy protons responsible for the strong ESP event observed by near-Earth spacecraft on 14 July 2012. The CME that generated this ESP event is simulated with the linear force-free spheromak model of EUHFORIA (Scolini et al. 2019; Verbeke et al. 2019) using parameters derived from remote-sensing observations. To reproduce the observed intensity-time profiles at Earth, we include in the PARADISE model the effects produced by the turbulent foreshock and sheath regions formed upstream and downstream of the shock, respectively. Within these regions, the parallel mean free path of the energetic particles decreases towards the shock wave. By comparing the simulation results with in situ data, it is demonstrated that our modelling approach successfully reproduces several features of the observations. This comparison also highlights a number of caveats and limitations that must be kept in mind when using MHD solar wind models together with a test-particle approach to simulate ESP and SEP events. We discuss each of these in detail and suggest how future developments may address them and lead to better models of SEP events with, ultimately, predictive capabilities.

The paper is structured as follows. Section 2 contains a description of the ESP event. In Sect. 3 the EUHFORIA simulation of the solar wind and the CME associated with the ESP event is discussed. Section 4 provides the details of the PARADISE simulation set-up, and in Sect. 5 the results of the PARADISE simulations are discussed. Section 6 provides a summary of the present work and lists future improvements and applications of the EUHFORIA + PARADISE tool to consider.

2. The ESP event on 14 July 2012

On 14 July 2012, an ESP event was observed near Earth that lasted for approximately one day. As shown in Fig. 1a, the ESP onset was characterised by abrupt and simultaneous intensity rises in all ion energy channels of the Low-Energy Magnetic Spectrometer (LEMS120) of the Electron, Proton, and Alpha Monitor (EPAM; Gold et al. 1998) on board the Advanced Composition Explorer (ACE). Most of the ions measured by ACE/EPAM are protons (Marhavilas et al. 2015). Signatures of the ESP event were also detected by the Energetic and Relativistic Nuclei and Electron (ERNE) experiment (Torsti et al. 1995) on board the Solar and Heliospheric Observatory (SoHO) for protons with energies up to ∼30 MeV. The event originated in a halo-CME that erupted from the Sun on 12 July 2012, temporally associated with a X1.4/2B flare at 15:37 UT from Active Region 11520 at S15W01. The CME arrived near Earth in the form of a flux rope with leading and trailing edges observed at ∼06:00 UT on 15 July and at ∼05:00 UT on 17 July, respectively. This flux rope was preceded by a shock observed by ACE at 17:26 UT on 14 July.

thumbnail Fig. 1.

Ion intensity measurements of the ESP event near Earth. (a) 5-min averages of the ion intensity-time profiles measured by ACE/EPAM for different energy channels. Far from the shock, the lowest energy channel of ACE/EPAM/LEMS120 departs from the nominal energy spectrum due to artificial counts produced by electrons bypassing the magnetic deflection system of LEMS120 (Marhavilas et al. 2015). The vertical lines indicate, from left to right, the onset of the ESP event, the shock arrival, and the leading edge of the flux rope. (b) Energy spectrum of the intensities observed by ACE/EPAM/LEMS120 and LEMS30, observed by ERNE/LED, and extracted from the SEPEM RDS v2.0 at the time of the shock arrival.

The peak intensities of the ESP event coincided with the arrival of the CME-driven shock. The energy spectrum of the observed energetic particle intensities at the shock arrival time is shown in Fig. 1b. To derive this energy spectrum, we included the available data from the ACE/EPAM/LEMS30 telescopes (Gold et al. 1998) and the intensity values for E < 15 MeV from the SOHO/ERNE Low Energy Detector (LED) and from the Solar Energetic Particle Environment Modelling Project (SEPEM) reference data set (RDS) v2.0 (sepem.eu; Jiggens et al. 2018). The energy channels used are indicated in the insets of Fig. 1b. As shown in this figure, the spectrum can be fitted by a power law with an exponential rollover. We find that the rollover energy is located around ∼2 MeV and the exponent of the power law at lower energies is equal to −1.1. Moreover, as can be deduced from Fig. 1a, the energy spectrum at low energies (< 1.90 MeV) is relatively flat prior to the passage of the associated interplanetary shock. This flattening is often observed prior to intense ESP events (e.g. Lario et al. 2018) and may result from a balance between the competing processes of DSA at the shock wave and adiabatic cooling in the upstream region (Prinsloo et al. 2019).

The ESP intensities measured by EPAM and ERNE showed a sudden drop on 15 July, around 06:00 UT, which corresponds to the arrival time of the leading edge of the flux rope. During the passage of this structure, the particle intensities remained low in all energy channels.

3. Simulating the CME with EUHFORIA

In order to model the ESP event, it is important to first have a realistic simulation of the large-scale disturbances associated with the propagating CME. Recently, Scolini et al. (2019) used EUHFORIA to model the propagation of the CME through the heliosphere for the event under study. This was done using the spheromak CME model, which approximates the magnetic cloud of a CME by a linear force-free spheroidal magnetic field. The magnetic parameters of the spheromak were estimated by Scolini et al. (2019) using images obtained by the Helioseismic and Magnetic Imager (HMI) instrument and the Atmospheric Imaging Assembly (AIA) on board the Solar Dynamic Observatory (SDO; Lemen et al. 2012). The speed and the shape of the CME were estimated from images taken by the coronagraphs on board the STEREO-A and B spacecraft, using the graduated cylindrical shell (GCS) model (Thernisien et al. 2009; Thernisien 2011). In this work, we use the same solar wind configuration and CME parameters as determined by Scolini et al. (2019). The only difference is that the spheromak used in the present work is assumed to have a density of 1.6 × 10−18 kg m−3 instead of 1.0 × 10−18 kg m−3. This increased CME mass density provides an improved match between the arrival time of the modelled and observed CME. The spheromak simulation parameters are summarised in Table 1.

Table 1.

Input parameters of the spheromak CME model used in the EUHFORIA simulation.

Figure 2 shows the EUHFORIA simulation result at the arrival of the CME at Earth, compared to in situ solar wind and magnetic field measurements provided by the OMNI database1 in geocentric solar ecliptic (GSE) coordinates. It is seen that the simulation captures the CME arrival time and the solar wind speed jump at the shock fairly well. There is also a reasonable qualitative agreement between the observed and modelled magnetic field components inside the CME. In Fig. 2 the arrival time of the simulated flux rope is indicated, which was determined following the methodology of Scolini et al. (2021). At the entry into the CME flux rope, the measured magnetic field exhibited a sudden increase, whereas in the simulation the arrival of the flux rope occurred ∼9 h later, showing a much more gradual increase in the magnetic field magnitude.

thumbnail Fig. 2.

EUHFORIA time series at Earth (green), compared to 1 min OMNI data (blue). From top to bottom: speed, proton number density, magnetic field strength B, and Bx, By, Bz components in GSE coordinates. The vertical lines indicate the onset times of the observed (dotted line) and simulated (dashed line) flux rope.

In order to study CME-driven shocks in EUHFORIA simulations, we developed a shock finder tool (see Appendix A). As explained in Sect. 4, this tool is also used to inject energetic particles continuously at the shock front and to prescribe regions of increased particle scattering in both the foreshock and the sheath region. Since SEPs are believed to propagate mainly along interplanetary magnetic field (IMF) lines, it is important to determine the point on the shock surface with which an observer connects magnetically. This point is called the cobpoint (Connecting with the OBserver point; after Heras et al. 1995). Due to the outward propagation of the CME in space and the nominal IMF topology, the cobpoint typically moves longitudinally in clockwise direction along the shock front as seen from the north ecliptic pole.

Figure 3 shows the evolution of the cobpoint for an observer at Earth’s location in the EUHFORIA simulation. In the simulation, this observer (hereafter Earth) connects with the shock front about ∼30 h after the CME insertion (i.e. on 13 July, around 22:00 UT). However, the observed onset of the SEP event suggests that Earth had likely a direct magnetic field connection to the shock wave shortly after the CME eruption (i.e. on 12 July around 17:00 UT). The discrepancy between the SEP event onset and the late magnetic connection in the EUHFORIA simulation can have two reasons: (i) the shock wave driven by the CME was in reality much wider than what is modelled and/or (ii) the IMF was significantly more radial than the modelled IMF. The latter explanation can be verified by inspecting the radial magnetic field component in Fig. 2 (Bx in GSE coordinates). On average, EUHFORIA underestimated Bx by a factor of ∼4 in the days between the CME eruption and the shock arrival at Earth, whereas the magnitude of the By and Bz components were reasonably well captured by the EUHFORIA simulation. The strongly radial IMF might be due to a preceding high-speed stream (HSS), which crossed Earth around 11 July 2012. Magnetic footpoint motion at the Sun and shearing of the magnetic field in the rarefaction region may produce a so-called sub-Parker IMF (e.g. Fisk 1996; Fisk et al. 1999; Schwadron & McComas 2005; Schwadron et al. 2021), for which the radial magnetic field component is significantly larger than that of a standard Parker spiral magnetic field. Such an IMF configuration was recently utilised by Schwadron et al. (2021) to explain energetic particle observations associated with a CIR. Even though EUHFORIA captures the HSS passing Earth on 11 July (green region in Fig. 3b above the Earth’s location indicated by the black dot), the model does not reproduce the sub-Parker IMF. This is because EUHFORIA does not take into account the footpoint motions and the interchange reconnection that are believed to generate the sub-Parker spiral.

thumbnail Fig. 3.

Simulated shock characteristics at Earth’s cobpoint. Panels a–c: give, respectively, the heliospheric radial distance, latitude, and longitude of Earth’s cobpoint in heliocentric Earth euatorial (HEEQ) coordinates as a function of the time elapsed after the CME injection at 0.1 au. Panels d–f: give, respectively, the shock speed, the shock compression ratio, and the shock angle at Earth’s cobpoint. The dashed vertical lines indicate the shock arrival time at Earth. Panel g: a snapshot of the solar wind radial speed at constant latitude θ = 4.5° close to the time when Earth established magnetic connection with the shock front. Earth is indicated by a black dot; it is magnetically connected to the cobpoint (magenta dot), which resides on the western flank of the CME.

A parameter that strongly influences the particle acceleration efficiency at shock waves is the shock obliquity (e.g. Ellison et al. 1995), which is determined by the angle θBn between the upstream magnetic field and the unit normal to the shock front. Figure 3f shows how θBn at the cobpoint evolves from ∼85° to ∼50° degrees during the interplanetary propagation of the shock prior to its arrival at Earth. Hence, the cobpoint samples shock regions that evolve from a quasi-perpendicular to more oblique shock geometries. We note that a sub-Parker spiral would consistently decrease the shock angle at the cobpoint. For example, if we multiply the radial magnetic field component of EUHFORIA by a factor of 4, a shock angle of ∼30° is obtained when the simulated shock crosses Earth. In the EUHFORIA simulation, such a shock geometry is obtained in the eastern flank of the CME, to which Earth is connected while residing in the sheath. Furthermore, we note that just after the shock crossing, the heliospheric current sheet (HCS) crosses Earth in the simulation, explaining the ∼180° jump in the shock angle in Fig. 3f.

Besides the obliquity, the shock compression ratio rc also plays an important role in DSA. Denoting the upstream and downstream solar wind velocity in a frame co-moving with the shock by, respectively, u1 and u2, we express the shock compression ratio as rc := u1n/u2n, where the subscript n refers to the vector component along the unit shock normal. Figure 3e illustrates that the compression ratio rc is close to 3.9, and remains relatively constant across the simulated shock surface. Steady-state DSA predicts an energy spectrum j ∝ E−Γ, with spectral index Γ = 0.5(rsc + 2)/(rsc − 1) and rsc the scattering-centre compression ratio (Bell 1978; Vainio & Schlickeiser 1999). Under the assumption that the scattering centres are frozen into the solar wind plasma, we have that rsc = rc. Using this approximation, we find that Γ = 1.02, which is slightly harder than the power law component of the observed energy spectrum shown in Fig. 1b.

4. PARADISE set-up

To model the ESP event with PARADISE, we inject 50 keV protons continuously at the CME-driven shock wave. Similar to Prinsloo et al. (2019) and Wijsen et al. (2021), the injected particle intensity jinj is scaled as

(1)

with Vsw denoting the solar wind velocity and r the heliocentric radial coordinate. This is done because −∇⋅Vsw gives a measure of the local compression and therefore of adiabatic heating of the background solar wind. A region where −∇⋅Vsw ≫ Vsw/r may therefore contain a preheated seed population, which can more readily be injected in the DSA mechanism (Prinsloo et al. 2019). The r−2 dependence in Eq. (1) takes into account the expansion of the solar wind and hence a possible dependence for the seed population into the interplanetary medium.

It should be noted, however, that the origin of the particle seed populations of SEP events is poorly understood. This is also tightly coupled to the so-called injection problem of DSA (e.g. Jokipii 1987), which addresses the conundrum that particles need a sufficiently high initial energy to commence DSA. This is necessary to allow the particles to propagate from the downstream to the upstream shock region. The particles’ injection energy is expected to depend strongly on the shock obliquity and the particle diffusion conditions (see e.g. Giacalone & Jokipii 1999; Zank et al. 2006; Neergaard Parker & Zank 2012). In particular, the injection energy may be significantly higher for a quasi-perpendicular shock, unless the particles are subjected to an efficient cross-field diffusion process (e.g. Zank et al. 2006; Neergaard Parker et al. 2014). By injecting 50 keV protons in the PARADISE simulation, we do not address any effect of the energy spectrum of the seed population on the modelled particle acceleration. However, since the injected protons get accelerated in the simulation, we can infer that 50 keV is above the injection energy of the modelled shock under the assumed diffusion conditions.

PARADISE solves the focused transport equation (FTE; e.g. Roelof 1969; Isenberg 1997) including a spatial diffusion process perpendicular to the IMF and a pitch-angle diffusion process. The pitch-angle diffusion process results from the interaction between the energetic particles and Alfvénic slab-like turbulence. PARADISE uses the results of quasi-linear theory (QLT; Jokipii 1966) to prescribe the following pitch-angle diffusion coefficient (see Agueda & Vainio 2013; Wijsen et al. 2019):

(2)

where μ denotes the cosine of the pitch-angle, ϵ = 0.048 is a parameter bridging the resonance gap at μ = 0 (e.g. Klimas & Sandri 1971), and D0 is a scaling factor that depends on the particle’s rigidity. In this work, D0 is determined by prescribing the particles’ parallel mean free path λ∥, which relates to Dμμ through (Hasselmann & Wibberenz 1970)

(3)

with v denoting the particle speed. During ESP events, the level of slab turbulence is typically enhanced in front of the shock wave due to the amplification of Alfvén waves by the energetic particles that are undergoing acceleration (Lee 1983). Similarly, the levels of turbulence are also expected to be enhanced in the sheath as compared to the ambient solar wind (e.g. Masías-Meza et al. 2016). As a result, both the foreshock and the sheath may be characterised by a small λ∥. Since this plays a crucial role in the production of ESP events, it needs to be included in the PARADISE simulations. This is achieved by prescribing a parallel mean free path of the form

(4)

(5)

(6)

(7)

where is the interplanetary parallel mean free path, R is the particle rigidity, d is the distance to the shock, x is the distance to the shock measured along the IMF, Δfeb represents a free escape boundary, and a, b, and β are free parameters which are specified below. The reference rigidity is chosen as Rref = 12.88 MV, corresponding to 88 keV protons. The assumed rigidity dependence is based on QLT and a Kolmogorov turbulence spectrum, that is, we choose β = 1/3 unless specified otherwise. The reference shock angle is chosen as , equal to the modelled shock angle when the CME arrives at Earth. The free-escape boundary is assumed to be Δfeb = 0.17 au. This value was obtained by multiplying the modelled shock speed at Earth (∼700 km s−1) by the duration of the observed ESP onset phase (∼10 h). The chosen value for Δfeb is also close to the presumed width of the sheath. This width can be estimated by taking into account that the arrival times of the shock and the leading edge of the flux rope are ∼12 h apart. Multiplying this by the shock speed gives a width of the sheath of ∼0.2 au.

The linear dependence of f(s) close to the shock is based on Eq. (22) of Vainio et al. (2014). In that work, the authors construct a semi-analytical model for ESP events, by using the analytical model of DSA derived by Bell (1978) as a starting point and re-calibrating the theoretical functional forms using numerical simulations of self-generated waves. Furthermore, we included a dependence on cos θBn in because quasi-parallel shocks are expected to have a more turbulent foreshock. This is because the quasi-parallel geometry allows particles to escape more easily upstream of the shock, where they subsequently induce beam-driven instabilities that excite various plasma waves that are expected to decrease the parallel mean free path of the particles (see e.g. Balogh & Treumann 2013, and reference therein).

When using an MHD model such as EUHFORIA to simulate the evolution of a CME in the inner heliosphere, the shock wave generated in front of the CME will typically be thicker than a real interplanetary shock wave, due to the limited resolution of the numerical grid. As explained in Appendix B, prescribing a small λ∥ across the simulated shock would lead to compressional acceleration instead of first-order Fermi acceleration. The former happens because the diffusion length scale is smaller than the shock width (e.g. Schwadron et al. 2020; Wijsen 2020). To avoid this and hence to allow particles to cross the simulated shock multiple times, we allow the particles to cross the shock wave in a scatter-free manner, that is, the particles propagate across the modelled shock structure while conserving their magnetic moment. Hence, unless a particle is mirrored by the magnetic compression, it will cross the entire shock wave structure and sample the full compression ratio. Since the particles still need a finite time to interact with the modelled shock wave, the acceleration efficiency of the modelled DSA is expected to be slightly reduced as compared to the case where the shock is a sharp discontinuity.

The parameter b in Eq. (7) gives the parallel mean free path for 88 keV protons at the shock (x = 0), and can be estimated from the measured particle intensities. This is done by approximating the observed intensity-time profiles by a function of the form f = αet/t*, with t* being the e-folding timescale. This procedure is illustrated in Fig. 4a for two different time intervals prior to the shock arrival (i.e. T1 = [0 h, 0.5 h] and T2 = [0.5 h, 10.3 h]). We note that the 57 keV energy channel has a high background level, likely affecting the intensity-time profiles in the first two time-intervals. Figure 4a also shows that the intensity time-profiles remain approximately constant after the shock passage, which is in agreement with DSA (e.g. Bell 1978; Drury & Voelk 1981). Assuming constant spatial diffusion coefficients in the two time intervals, steady state DSA can be used to express as (e.g. van Nes et al. 1984; Beeck & Sanderson 1989; Giacalone 2012)

(8)

thumbnail Fig. 4.

Parallel mean free path in the foreshock derived from in situ ion intensity measurements. Panel a: dots give the intensity-time profiles observed by EPAM for different energy channels, shifted downwards to avoid overlap. The straight lines give a time-averaged increase in the time intervals T1 = [0 h, 0.5 h] (blue) and T2 = [0.5 h, 10.3 h] (green), and the values above each line give the slope. Panel b: the dotted lines give the derived from the slopes presented in (a) for T1 (blue background) and T2 (green background). The solid line gives the f(s) (see Eq. (7)) used in the PARADISE simulations.

where Vsh is the shock speed. The corresponding parallel mean free path can then be calculated as

(9)

The values for obtained for the different energy channels are shown as dotted lines in Fig. 4b. To obtain these values, the EUHFORIA simulation was used to prescribe θBn = 49°, Vsh = 683 km s−1, and un1 = −401 km s−1.

The values obtained for should be treated with care as steady-state DSA does not take the time history of the (fore)shock or the effect of magnetic focusing due to the diverging IMF into account. The latter process may especially affect the intensity-time profiles in the time interval T2. For this reason, we determine the parameters a and b solely based on the first interval T1, where the diffusive approximation is assumed to be more valid. We chose the parameters a and b such that

(10)

where t = 0.25 h is the midpoint of T1. This gives b = 1.1 × 10−3 au and a = 2.1 × 10−2. The resulting mean free path in the foreshock for 88 keV protons and θBn = 49° is shown as a solid line in Fig. 4b (see also Fig. B.1). The parameters Δfeb, a, and b are kept constant during the entirety of the simulation. Hence, Eqs. (4)–(7) only contain an implicit time dependence due to the propagation and expansion of the shock wave and due to the dependence of on the shock obliquity. In reality, λ∥ is expected to be time dependent due to the coupling between the accelerated energetic particle distributions and the turbulence (e.g. Lee 1983). Capturing the wave-particle dynamics is, however, outside the scope of the current study, but will be the subject of future improvements to our model.

Our simulations include a weak cross-field diffusion process, which is characterised by a constant perpendicular mean free path λ⊥ = 10−4 au. Hence, we have that λ∥/λ⊥ > 1 everywhere in the simulation, meaning that parallel particle transport will dominate. We leave it for future work to study the effect of different cross-field diffusion conditions.

Finally, we note that the EUHFORIA snapshots are updated every 15 min in PARADISE. Apart from the MHD plasma variables, each snapshot contains the location of the shock wave and λ∥, computed according to Eqs. (4)–(7). A linear interpolation in time is performed to obtain the plasma variables and λ∥ at times between two consecutive snapshots.

5. Results

5.1. Intensities at different longitudes

Figure 5 shows the simulated omnidirectional intensity-time profiles at Earth. The particle distribution is scaled such that the simulated peak intensity matches the observed EPAM peak intensity in the 68.1−115 keV channel. The figure shows that the PARADISE simulation successfully reproduces several features seen in the EPAM data, especially for energies below 1 MeV. For these energies, the simulation reproduces: (1) the onset time of the ESP event, (2) the intensity increase prior to the shock arrival, (3) the energy spectrum at the shock wave, (4) the intensity drop near the onset of the magnetic cloud, and (5) a secondary intensity peak just before the onset of the magnetic cloud. Both in the simulation and in the data, the secondary peak occurs just after the three magnetic field components become negative and subsequently attain a local minimum. This structure is more smeared-out in the EUHFORIA simulation as compared to the data, explaining why the secondary intensity peak in the PARADISE simulation is lower in intensity but spans a longer time interval than in the EPAM data.

thumbnail Fig. 5.

Simulated (solid lines) and observed (dots) omnidirectional intensity-time profiles at Earth. The dots in the background give the ACE/EPAM/LEMS120 data (see also Fig. 1a). The dashed vertical lines indicate the arrival times of the shock (left) and the flux rope (right) in the EUHFORIA simulation.

The simulation captures only the ESP event, and not the preceding SEP component that started on 12 July at around 17:00 UT. This is due to the late (by ∼30 h) magnetic connection of Earth to the simulated shock front (see Sect. 3). Moreover, the peak intensities of the high-energy channels are reached after the shock arrival at Earth. This is because after the shock passage, the cobpoint of Earth moves to the eastern flank of the CME, where the acceleration is more efficient due to a more parallel shock geometry. To illustrate this, we show in Fig. 6 the intensity-time profiles of two virtual observers located at the same radial distance and latitude as Earth but shifted in longitude, that is, one virtual spacecraft located 30° east of Earth, which we refer to as the E30 observer, and another virtual spacecraft located 30° west of Earth, referred to as the W30 observer. These observers are indicated by the green and red symbols in Fig. 6d, which shows the omnidirectional particle intensities together with the solar wind speed in a slice of constant latitude, approximately 6 h before the shock arrival at Earth. The figure illustrates that the particle peak intensities are obtained at the CME shock wave. Moreover, the highest particle intensities are obtained on the eastern CME wing, where the shock is quasi-parallel.

thumbnail Fig. 6.

Simulation results for virtual observers at different longitudes. Panels a and b: give the omnidirectional intensity-time profiles (upper panel) together with solar wind speed (lower panel) for the W30 and E30 observers, respectively. Panel c: gives the peak intensities versus energy for different observers. Panel d: is a snapshot of the omnidirectional particle intensity plotted in grey shades on top of the solar wind speed at constant HEEQ latitude θ = 4.5°.

The W30 observer (green symbol) has an early connection to the shock wave, and hence starts observing energetic particles already on 14 July, around 00:00 UT. This is still more than a day later than the onset of the SEP event observed at Earth. This is because we do not include particle acceleration and transport happening in the corona and we only start injecting 50 keV protons at the time the CME shock is formed in the EUHFORIA domain. Figure 6a shows that during the simulated event onset, the energy spectrum is increasing up to 115 keV, flat from 115 keV to 580 keV, and decreasing for higher energies. This can be attributed to the non-trivial interplay between the acceleration efficiency at the shock, the particle escape efficiency from the foreshock, the adiabatic deceleration in the solar wind, and the time a particle needs to travel from the shock to the spacecraft. All four processes are energy dependent.

In contrast to the simulated intensity profiles at Earth, the particle intensities peak in all energy channels at the time of the shock arrival at the W30 observer. Figure 6c also reveals that the energy spectrum of the W30 observer is slightly harder than at Earth, despite the fact that the shock encountered by the W30 observer is slightly weaker (rc = 3.8) than at Earth (rc = 3.9). This can be attributed to the more parallel shock geometry of the eastern CME wing. The inverse proportionality of on cos θBn (see Eq. (5)) allows particles to be more easily trapped at a quasi-parallel shock, increasing the probability on multiple shock crossings.

Due to its position, the E30 observer (red symbol) connects around 05:00 UT on 14 July to the CME’s western flank, which is characterised by a quasi-perpendicular shock geometry. Along this flank, the acceleration is less efficient partly due to the small foreshock, which is because of the cos θBn dependence in Eq. (5). In addition, the particles residing on the field lines crossing the western flank of the CME only had limited time to interact with the shock wave and hence get accelerated. As for the observer at Earth, the high-energy channels peak after the shock passage, which is because the cobpoint of the E30 observer moves towards the shock nose.

5.2. Intensities at different radial distances

Figure 7 shows the energy spectra measured at the peak intensities for different observers located along the Sun-Earth line, hence illustrating the time-evolution of the particle energy spectrum at a point close to the shock nose. The peak intensities have been scaled by a r2 factor, to cancel out the effect of the solar wind expansion, which is incorporated in our injection distribution through the r−2 dependence in Eq. (1). It is seen that the energy spectrum hardens with radial distance up to ∼0.5 au. After that, the shape of the energy spectrum remains approximately constant and can be characterised by a power-law with an exponential rollover at ∼700 keV. Although not shown here, the energy spectrum shows a similar evolution for observers located along the Sun-W30 line. In this case, the exponential rollover is located around 1 MeV (see also Fig. 6c). In future studies, we plan to compare the modelled radial behaviour of the energy spectrum with observations from PSP (Fox et al. 2016) and Solar Orbiter (SolO; Müller et al. 2020).

thumbnail Fig. 7.

Peak intensities scaled by r2 versus energy for virtual observers located along the Sun-Earth line.

5.3. Intensities for different rigidity dependences

In the simulation results presented above, the parallel mean free path followed the QLT rigidity dependence (i.e. β = 1/3 in Eqs. (4)–(7)). This means that the parallel mean free path increases with the particle energy. However, in the steady-state DSA theory derived by Bell (1978), the parallel mean free path decreases with increasing energy for a shock wave with rc ∼ 4. Beeck & Sanderson (1989) analysed five in situ observed foreshock regions and found both positive and negative rigidity dependences. This result was obtained using a method similar to the one presented in Sect. 4 (see Eq. (8)). Moreover, by calculating the scattering mean free path directly from the observed magnetic field fluctuations, Beeck & Sanderson (1989) obtained comparable rigidity dependences and mean free path values when assuming that the particle scattering resulted from a combination of Alfvén waves and magnetosonic waves. In contrast, when only considering the resonance scattering due to Alfvénic fluctuations, the authors derived smaller mean free paths and positive rigidity dependences. Such a positive dependence was also found by Afanasiev et al. (2015) using the SOLar Particle Acceleration in Coronal Shocks (SOLPACS) model, which takes into account the resonant interactions between energetic particles and Alfvén waves in the foreshock, under the assumption of QLT. The main difference of this modelling to the theory of Bell (1978) (in addition to time dependence) is that the wave-particle resonance condition is taken into account in full, accounting for the effect of the pitch angle that is omitted in the theory. This leads to low-energy protons being scattered also by waves generated by high-energy protons especially further out from the shock.

To better understand the effect of the rigidity dependence on the results presented in this work, we performed simulations with β equal to 1/3, 0, and −1/3. Figure 8 shows the peak intensity energy spectra obtained by a spacecraft located at Earth for these three simulations. It can be seen that decreasing β produces a harder energy spectrum. This is because the high-energy particles become more efficiently trapped near the shock and hence they have a larger probability of multiple shock encounters. The acceleration efficiency of the high-energy particles is thus increased. The β = −1/3 case gives the best correspondence with the EPAM data. However, this contrasts with the positive energy dependence derived from the intensity observations at 1 au in time interval T1, using steady state DSA (see Fig. 4b).

thumbnail Fig. 8.

Peak intensity energy spectra obtained from different λ∥ power-law dependences on the ion rigidity. The orange, green, and red symbols give the results of PARADISE simulations using the indicated rigidity dependence. Blue symbols give the peak intensities measured by EPAM.

6. Summary and conclusions

In this work we have presented the simulation results of the ESP event that was observed near Earth on 14 July 2012. This was done by using the MHD model EUHFORIA together with the energetic particle model PARADISE. EUHFORIA was used to model the background solar wind and the propagation of the CME through interplanetary space. The CME was modelled using the spheromak model of EUHFORIA, for which the magnetic and kinematic parameters had previously been derived by Scolini et al. (2019). The resulting EUHFORIA simulation manages to capture both the shock arrival time and the solar wind speed jump at the shock passage reasonably well, and a qualitative match is obtained between the observed and simulated flux rope passage at Earth. One important discrepancy is that in the simulation Earth connects magnetically to the shock wave at ∼0.5 au, whereas the SEP data suggest that Earth established magnetic connection with the shock when the CME was still in the low corona. This difference is, at least partly, the result of the IMF being a sub-Parker field due to the passage of an HSS on 11 July. Even though this HSS was captured by EUHFORIA, a sub-Parker IMF configuration is not reproduced by the model since EUHFORIA does not take into account effects of differential rotation and magnetic footpoint motion, which are believed to produce sub-Parker spirals. Due to this late connection, the simulation cannot capture the onset phase of the SEP event. This illustrates that one should be careful when using an MHD model such as EUHFORIA or a Parker spiral IMF to determine the magnetic connectivity of an observer.

A good agreement was obtained between the simulated and observed ESP event, especially at low energies (< 700 keV). In our simulations, the energetic protons are accelerated at the CME-driven shock wave, assuming a seed population of 50 keV protons. The assumed seed population (see Eq. (1)) may originate from the solar wind suprathermal tail, especially near the shock wave where the solar wind gets shock-heated, hence producing a larger population of suprathermal protons (Prinsloo et al. 2019). The good match between the observed and modelled intensity-time profiles confirms that the low-energy component of ESP events observed at 1 au may indeed be produced by the interplanetary acceleration of particles at the CME-driven shock wave. The simulations do not reproduce the intensity-time profiles for energies above ∼1 MeV. This energy threshold is close to the rollover energy of the observed energy spectrum. This correspondence suggests that the ESP component above 1 MeV depends on particle acceleration happening in the solar corona. This will generate an extra, more energetic seed population that can be re-accelerated in the interplanetary medium along with the suprathermal ions of the solar wind. In future work, we plan to extend our model to the corona to investigate this possibility further.

In our simulation, the particle intensity-time profiles show a drop when the flux rope arrives. A similar sharp drop is observed in the EPAM data and indicates that the particles have limited access into the flux rope. In addition, the interaction between the flux rope and the IMF produces a secondary intensity peak, which is present both in the simulation and the data. The observed particle intensities do not drop to the pre-event background levels, indicating that some energetic particles reside inside the flux rope. Although not shown here, we tried to reproduce this by increasing the perpendicular mean free path in the simulation to allow the simulated SEPs to penetrate the flux rope more easily through cross-field diffusion. This, however, quickly diffuses the sharp drop in intensity that is seen at the beginning of the magnetic cloud. Hence, we believe that the observed SEPs inside the flux rope are likely not the result of a spatial diffusion process happening everywhere along the interface between the sheath and the flux rope. Instead, particles may enter the flux rope at locations where the sheath magnetic field reconnects with the flux rope magnetic field. This is in line with the findings of Laitinen & Dalla (2021), who showed by using full orbit simulations that energetic particles penetrate flux ropes most easily near magnetic x points.

To reproduce the low-energy component of the ESP event, it was necessary to include a turbulent foreshock and sheath region in PARADISE, which we prescribed by using a parallel mean free path that decreases towards the shock wave. The parallel mean free path inside the foreshock was assumed to scale with 1/cos2θ to take into account the fact that the turbulence upstream of a shock wave is expected to depend on the shock angle, as is the case for Earth’s bow shock (e.g. Eastwood et al. 2005). The width of the foreshock region and the parallel mean free path at the shock surface were estimated from the observed ESP event at Earth and did not contain any explicit time dependence or dependence on the intensity of the modelled energetic particles.

In order to transform a model such as EUHFORIA + PARADISE into a fully physics-based forecasting model for ESP or SEP events, it is important to have a model for the foreshock region and particle injection that does not rely on in situ observations of the ESP event. To achieve this, we plan to couple our model to the SOLPACS code (Afanasiev et al. 2015), which models particle acceleration at oblique shock waves, thereby taking into account the generation of Alfvén waves in the upstream region by the energetic particles themselves. Hence, this model provides a physics-based calculation of the parallel mean free path and its rigidity dependence inside the foreshock. This is imperative as the results presented in this work illustrate the sensitivity of the simulations on the assumed scattering conditions in the foreshock.

Solar energetic particle events are an important component of space weather since they can pose significant hazards for the health of astronauts and the microelectronics on board spacecraft. As a result, a profound understanding of the physical processes behind SEP events is necessary to be able to develop a reliable SEP forecasting tool. In this study, we presented the first attempt at using the EUHFORIA + PARADISE framework to reproduce an observed ESP event. With this work, we contribute to the ongoing worldwide effort (e.g. Li et al. 2021; Young et al. 2021) that aims at the development of a physics-based SEP forecasting tool. In the future, we plan to extend our study by modelling several of the recently observed multi-spacecraft SEP events of solar cycle 25.

2

As seen from a frame co-moving with the shock or compression wave.

Acknowledgments

We acknowledge the use of ACE data from the ACE Science Center and of SOHO/ERNE data from sepem.eu and of the SEPEM Reference Data Set version 2.00, European Space Agency (2016). N.W. acknowledges funding from the Research Foundation – Flanders (FWO – Vlaanderen, fellowship no. 1184319N). A.A. and B.S. acknowledge the support by the Spanish Ministerio de Ciencia e Innovación (MICINN) under grant PID2019-105510GB-C31 and through the “Centre of Excellence María de Maeztu 2020–2023” award to the ICCUB (CEX2019-000918-M). C.S. acknowledges the NASA Living With a Star Jack Eddy Postdoctoral Fellowship Program, administered by UCAR’s Cooperative Programs for the Advancement of Earth System Science (CPAESS) under award no. NNX16AK22G. D.L. acknowledges support from NASA Living With a Star (LWS) programs NNH17ZDA001N-LWS and NNH19ZDA001N-LWS, and the Goddard Space Flight Center Heliophysics Innovation Fund (HIF) program. The work in the University of Turku was performed under the umbrella of Finnish Centre of Excellence in Research of Sustainable Space. This project has received funding from the European Union’s Horizon 2020 research and innovation programs under grant agreement No. 870405 (EUHFORIA 2.0). These results were also obtained in the framework of the ESA project “Heliospheric modelling techniques” (Contract No. 4000133080/20/NL/CRS). Computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Centre), funded by the FWO and the Flemish Government-Department EWI.

References

  1. Afanasiev, A., Vainio, R., & Kocharov, L. 2014, ApJ, 790, 36 [Google Scholar]
  2. Afanasiev, A., Battarbee, M., & Vainio, R. 2015, A&A, 584, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Agueda, N., & Vainio, R. 2013, J. Space Weather Space Clim., 3, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Balogh, A., & Treumann, R. A. 2013, Physics of Collisionless Shocks (New York: Springer), 12 [CrossRef] [Google Scholar]
  5. Battarbee, M., Laitinen, T., & Vainio, R. 2011, A&A, 535, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Beeck, J., & Sanderson, T. R. 1989, J. Geophys. Res., 94, 8769 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bell, A. R. 1978, MNRAS, 182, 147 [Google Scholar]
  8. Bryant, D. A., Cline, T. L., Desai, U. D., & McDonald, F. B. 1962, J. Geophys. Res., 67, 4983 [NASA ADS] [CrossRef] [Google Scholar]
  9. Desai, M., & Giacalone, J. 2016, Liv. Rev. Sol. Phys., 13, 3 [Google Scholar]
  10. Drury, L. O., & Voelk, J. H. 1981, ApJ, 248, 344 [Google Scholar]
  11. Eastwood, J. P., Lucek, E. A., Mazelle, C., et al. 2005, Space Sci. Rev., 118, 41 [NASA ADS] [CrossRef] [Google Scholar]
  12. Ellison, D. C., Baring, M. G., & Jones, F. C. 1995, ApJ, 453, 873 [Google Scholar]
  13. Fisk, L. A. 1996, J. Geophys. Res., 101, 15547 [NASA ADS] [CrossRef] [Google Scholar]
  14. Fisk, L. A., Zurbuchen, T. H., & Schwadron, N. A. 1999, ApJ, 521, 868 [Google Scholar]
  15. Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7 [Google Scholar]
  16. Giacalone, J. 2012, ApJ, 761, 28 [NASA ADS] [CrossRef] [Google Scholar]
  17. Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204 [Google Scholar]
  18. Giacalone, J., Jokipii, J. R., & Kóta, J. 2002, ApJ, 573, 845 [Google Scholar]
  19. Giacalone, J., Jokipii, J. R., & Kóta, J. 2005, Geophys. Monogr. Ser., 156, 41 [Google Scholar]
  20. Gold, R. E., Krimigis, S. M., Hawkins, S. E., et al. 1998, in Electron, Proton, and Alpha Monitor on the Advanced Composition Explorer Spacecraft, eds. C. T. Russell, R. A. Mewaldt, & T. T. Von Rosenvinge (Dordrecht: Springer), 541 [Google Scholar]
  21. Gosling, J. T., Asbridge, J. R., Bame, S. J., et al. 1981, J. Geophys. Res., 86, 547 [NASA ADS] [CrossRef] [Google Scholar]
  22. Guo, F., Giacalone, J., & Zhao, L. 2021, Front. Astron. Space Sci., 8, 27 [Google Scholar]
  23. Hasselmann, K., & Wibberenz, G. 1970, ApJ, 162, 1049 [Google Scholar]
  24. Heras, A. M., Sanahuja, B., Lario, D., et al. 1995, ApJ, 445, 497 [Google Scholar]
  25. Hu, J., Li, G., Ao, X., Zank, G. P., & Verkhoglyadova, O. 2017, J. Geophys. Res.: Space Phys., 122, 10,938 [Google Scholar]
  26. Isenberg, P. A. 1997, J. Geophys. Res., 102, 4719 [Google Scholar]
  27. Jiggens, P., Heynderickx, D., Sandberg, I., et al. 2018, J. Space Weather Space Clim., 8, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Jokipii, J. R. 1966, ApJ, 146, 480 [Google Scholar]
  29. Jokipii, J. R. 1987, ApJ, 313, 842 [Google Scholar]
  30. Jokipii, J. R., & Giacalone, J. 2007, ApJ, 660, 336 [Google Scholar]
  31. Jokipii, J. R., Giacalone, J., & Kota, J. 2003, Int. Cosm. Ray Conf., 6, 3685 [Google Scholar]
  32. Klimas, A. J., & Sandri, G. 1971, ApJ, 169, 41 [Google Scholar]
  33. Kozarev, K. A., Evans, R. M., Schwadron, N. A., et al. 2013, ApJ, 778, 43 [Google Scholar]
  34. Laitinen, T., & Dalla, S. 2021, ApJ, 906, 9 [Google Scholar]
  35. Lario, D., Ho, G. C., Decker, R. B., et al. 2003, in Solar Wind Ten, eds. M. Velli, R. Bruno, F. Malara, & B. Bucci, AIP Conf. Ser., 679, 640 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lario, D., Berger, L., Wilson, L. B., III, et al. 2018, J. Phys. Conf. Ser., 1100, 012014 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lee, M. A. 1983, J. Geophys. Res., 88, 6109 [NASA ADS] [CrossRef] [Google Scholar]
  38. Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17 [Google Scholar]
  39. le Roux, J. A., & Webb, G. M. 2009, ApJ, 693, 534 [Google Scholar]
  40. le Roux, J. A., Webb, G. M., Florinski, V., & Zank, G. P. 2007, ApJ, 662, 350 [Google Scholar]
  41. Lewiner, T., Lopes, H., Vieira, A. W., & Tavares, G. 2003, J. Graphics Tools, 8, 1 [NASA ADS] [CrossRef] [Google Scholar]
  42. Li, G., Jin, M., Ding, Z., et al. 2021, ApJ, 919, 146 [Google Scholar]
  43. Manchester, W. B., IV, van der Holst, B., Tóth, G., & Gombosi, T. I. 2012, ApJ, 756, 81 [NASA ADS] [CrossRef] [Google Scholar]
  44. Marhavilas, P. K., Malandraki, O. E., & Anagnostopoulos, G. C. 2015, Planet. Space Sci., 117, 192 [NASA ADS] [CrossRef] [Google Scholar]
  45. Masías-Meza, J. J., Dasso, S., Démoulin, P., Rodriguez, L., & Janvier, M. 2016, A&A, 592, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Müller, D., Zouganelis, I., St. Cyr, O. C., Gilbert, H. R., & Nieves-Chinchilla, T. 2020, Nat. Astron., 4, 205 [CrossRef] [Google Scholar]
  47. Neergaard Parker, L., & Zank, G. P. 2012, ApJ, 757, 97 [Google Scholar]
  48. Neergaard Parker, L., Zank, G. P., & Hu, Q. 2014, ApJ, 782, 52 [Google Scholar]
  49. Ng, C. K., & Reames, D. V. 2008, ApJ, 686, L123 [Google Scholar]
  50. Pomoell, J., & Poedts, S. 2018, J. Space Weather Space Clim., 8, A35 [Google Scholar]
  51. Pomoell, J., Aran, A., Jacobs, C., et al. 2015, J. Space Weather Space Clim., 5, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Prinsloo, P. L., Strauss, R. D., & Le Roux, J. A. 2019, ApJ, 878, 144 [Google Scholar]
  53. Reames, D. V. 2017, Solar Energetic Particles: A Modern Primer on Understanding Sources, Acceleration and Propagation (Springer-Verlag) [CrossRef] [Google Scholar]
  54. Roelof, E. C. 1969, in Lectures in High-Energy Astrophysics, eds. H. Ögelman, & J. R. Wayland, 111 [Google Scholar]
  55. Schwadron, N. A., & McComas, D. J. 2005, Geophys. Res. Lett., 32, L03112 [Google Scholar]
  56. Schwadron, N. A., Townsend, L., Kozarev, K., et al. 2010, Space Weather, 8, S00E02 [Google Scholar]
  57. Schwadron, N. A., Gorby, M., Török, T., et al. 2014, Space Weather, 12, 323 [NASA ADS] [CrossRef] [Google Scholar]
  58. Schwadron, N. A., Lee, M. A., Gorby, M., et al. 2015, ApJ, 810, 97 [Google Scholar]
  59. Schwadron, N. A., Bale, S., Bonnell, J., et al. 2020, ApJS, 246, 33 [Google Scholar]
  60. Schwadron, N. A., Joyce, C. J., Aly, A., et al. 2021, A&A, 650, A24 [EDP Sciences] [Google Scholar]
  61. Scolini, C., Rodriguez, L., Mierla, M., Pomoell, J., & Poedts, S. 2019, A&A, 626, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Scolini, C., Dasso, S., Rodriguez, L., Zhukov, A. N., & Poedts, S. 2021, A&A, 649, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Thernisien, A. 2011, ApJS, 194, 33 [Google Scholar]
  64. Thernisien, A., Vourlidas, A., & Howard, R. A. 2009, Sol. Phys., 256, 111 [Google Scholar]
  65. Torsti, J., Valtonen, E., Lumme, M., et al. 1995, Sol. Phys., 162, 505 [CrossRef] [Google Scholar]
  66. Vainio, R., & Afanasiev, A. 2018, Particle Acceleration Mechanisms (Cham: Springer International Publishing), 45 [Google Scholar]
  67. Vainio, R., & Laitinen, T. 2007, ApJ, 658, 622 [Google Scholar]
  68. Vainio, R., & Schlickeiser, R. 1999, A&A, 343, 303 [NASA ADS] [Google Scholar]
  69. Vainio, R., Pönni, A., Battarbee, M., et al. 2014, J. Space Weather Space Clim., 4, A08 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. van der Holst, B., Manchester, W. B., IV, Frazin, R. A., et al. 2010, ApJ, 725, 1373 [NASA ADS] [CrossRef] [Google Scholar]
  71. van Nes, P., Roelof, E. C., Reinhard, R., Sanderson, T. R., & Wenzel, K. P. 1984, Adv. Space Res., 4, 315 [Google Scholar]
  72. Verbeke, C., Pomoell, J., & Poedts, S. 2019, A&A, 627, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Wijsen, N. 2020, PhD Thesis, KU Leuven (Belgium) and Univ. Barcelona (Spain), https://ui.adsabs.harvard.edu/abs/2020PhDT.........8W/abstract [Google Scholar]
  74. Wijsen, N., Aran, A., Pomoell, J., & Poedts, S. 2019, A&A, 622, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Wijsen, N., Samara, E., Aran, À., et al. 2021, ApJ, 908, L26 [CrossRef] [Google Scholar]
  76. Young, M. A., Schwadron, N. A., Gorby, M., et al. 2021, ApJ, 909, 160 [Google Scholar]
  77. Zank, G. P., Li, G., Florinski, V., et al. 2006, J. Geophys. Res.: Space Phys., 111, A06108 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: The shock finder algorithm

In order to study CME-driven shocks in EUHFORIA simulations, we updated the shock finder tool that was developed by Pomoell et al. (2015) for a 2D MHD model. This new tool is used in PARADISE to both inject particles continuously at the shock front and prescribe the regions of increased particle scattering conditions in the foreshock and the sheath.

The detection of a CME-driven shock wave is complicated by, for example, the presence of shocks associated with CIRs. The shock-detection algorithm removes shocks driven by high-speed solar wind streams by subtracting the unperturbed background solar wind configuration from the solar wind in which the CME is propagating:

(A.1)

where the subscript ‘sw’ refers to a solar wind simulation that does not contain a CME and q can be any plasma variable. With this definition, a CME corresponds thus to a qrel ≠ 0 region in the simulation domain. Moreover, since the EUHFORIA solar wind is steady state in a frame corotating with the sun, we have that qsw(t, r, θ, ϕ) = qsw(0, r, θ, ϕ − Ωt), with Ω the solar rotation rate. Hence, qsw(t) can be derived from a single snapshot, obtained, for example, just before the CME insertion.

In a next step, we search for a specific point on the shock surface. If we are interested in deriving the entire shock surface, qrel is calculated first along a radial line passing through the centre of the CME insertion location. If instead we are only interested in the shock properties at the cobpoint of a particular observer, qrel is calculated along the IMF line crossing the observer. For q, we typically choose both the solar wind speed V and the entropy S ≡ P/ργ, where P is the thermal pressure, ρ the density, and γ the adiabatic index. Following either a radial line segment or an IMF line from the outer boundary of the domain inwards, the shock point candidate r1 is determined as the first point for which both Vrel > ϵV and Srel > ϵS, where ϵV and ϵS are predefined positive thresholds. Such non-zero thresholds are necessary due to the finite numerical accuracy of grid-based simulations. To ensure that r1 is indeed located at a compression wave, the algorithm checks whether the divergence of the solar wind velocity ∇ ⋅ V is negative. This condition implies the presence of converging solar wind flows, for which shocks wave are extreme examples.

Next, the shock surface S1 is determined around r1 by computing it as the isosurface of qrel = ϵq, with ϵq a positive constant. As before, different choices for q are possible. In this work, we choose qrel = −(∇ ⋅ V)rel, because the resulting isosurface encapsulates the simulated MHD shock wave. The isosurface is computed using the marching cubes algorithm (MCA), which generates the isosurface on a triangle mesh (Lewiner et al. 2003). Figure A.1 illustrates an application of the MCA to determine the shock surface generated by a cone CME in a EUHFORIA test simulation. From the mesh generated by the MCA, it is straightforward to compute the shock-normal at any point on the shock surface. The shock angle θBn can then be determined by interpolating the upstream magnetic field to the shock surface. The shock speed Vsh is calculated by determining the shock surface S2 at a time instance that is typically Δt = t2 − t1 ∼ 30 minutes back in time. The shock speed at a point ri ∈ S1 is then estimated as Δx/Δt, where Δx is the minimal distance from point ri to S2. The shock speed and the solar wind conditions just upstream of ri are then used to solve the Rankine–Hugoniot equations, which give us the downstream solar wind conditions everywhere along the shock surface.

thumbnail Fig. A.1.

Illustration of the shock finder algorithm. Left panel: Snapshot of a EUHFORIA simulation containing a cone CME. Shown is the radial speed in the solar equatorial plane. The red line indicates the shock surface, as determined by the MCA. Right panel: 3D view of the shock surface computed by the MCA. The resolution has been reduced to clearly illustrate the triangle mesh.

Appendix B: Particle acceleration in PARADISE

PARADISE solves the FTE, which contains the necessary physics to describe particle acceleration in converging, accelerating, and shear flows (e.g. le Roux et al. 2007; le Roux & Webb 2009). Moreover, as illustrated by, for example, le Roux et al. (2007), the momentum term of the FTE contains the energy changes due to the particle drifts along the motional electric field. In this work, the FTE was solved without including a diffusion coefficient Dpp, where p denotes the particle’s momentum magnitude. This means that our simulations do not take into account any stochastic acceleration mechanism that may occur downstream of the shock wave (e.g. Afanasiev et al. 2014).

When using a 3D MHD model such as EUHFORIA to simulate the evolution of a CME in the inner heliosphere, the shock wave generated in front of the CME will typically be much wider than a real interplanetary shock wave, due to the limited resolution of the numerical grid. In such MHD simulations, the shock wave is thus similar to a large amplitude compression wave, which can be characterised by the length scale (e.g. Vainio & Afanasiev 2018),

(B.1)

Here, U denotes the velocity of the solar wind plasma measured in a reference frame co-moving with the shock wave, and the subscript ‘n’ refers to the vector component measured along the shock normal. The efficiency of the particle acceleration near compression waves is related to the particle diffusion length, which can be defined as (e.g. Giacalone et al. 2002)

(B.2)

where κnn and λn denote the particles’ spatial diffusion coefficient and the mean free path normal to the shock front, respectively. The latter is related to the particles’ parallel and perpendicular mean free paths through

(B.3)

If LD ≫ LC, the particles can cross the compression wave multiple times and gain energy through first-order Fermi acceleration (e.g. Giacalone et al. 2002; Jokipii et al. 2003; Jokipii & Giacalone 2007). For this acceleration process to happen, the diffusive length scale upstream and downstream of the compression wave cannot be too large either (e.g. compared to the focusing length LF = (∇⋅b)−1), since otherwise particles will escape into the interplanetary medium (Battarbee et al. 2011; Vainio et al. 2014). In the other limiting case, when LD ≪ LC, the particles will be closely tied to the scattering centres embedded in the solar wind and hence they are advected with the plasma flow through the compression wave2. As a result, the particles will gain energy in an adiabatic fashion (e.g. Giacalone et al. 2005; Vainio & Afanasiev 2018; Schwadron et al. 2020). In this work, we propagate the particles scatter free across the shock wave, which is identified as the ∇ ⋅ Vsw < 0 region in front of the CME (see Appendix A). This ensures that LD ≫ LC, since LD = +∞ across the CME-driven shock wave. The foreshock and sheath regions, where a small mean free path is prescribed according to Eqs. (4)–(6), begin only outside the ∇ ⋅ Vsw < 0 region. Figure B.1 shows the parallel mean free path for 88 keV protons together with ∇ ⋅ Vsw. For clarity, we only show the parallel mean free path where au. The simulated shock is visible as the blue arc inside the green shaded region. The time step of the particles inside the modelled shock is calculated as if . This is done to avoid a large time step, which would transport the particles artificially far upstream or downstream of the shock in a single time step. Finally, we refer to chapter 8 of Wijsen (2020) for a discussion of different PARADISE simulations in which the parallel mean free path is kept small across the CME-driven shock wave, hence leading to compressional acceleration.

thumbnail Fig. B.1.

Snapshot of the parallel mean free path of 88 keV protons (green shades) and ∇ ⋅ Vsw (red and blue) at constant HEEQ latitude θ = 4°. The parallel mean free path is only shown where λ∥ < 0.1 au.

All Tables

Table 1.

Input parameters of the spheromak CME model used in the EUHFORIA simulation.

In the text

All Figures

thumbnail Fig. 1.

Ion intensity measurements of the ESP event near Earth. (a) 5-min averages of the ion intensity-time profiles measured by ACE/EPAM for different energy channels. Far from the shock, the lowest energy channel of ACE/EPAM/LEMS120 departs from the nominal energy spectrum due to artificial counts produced by electrons bypassing the magnetic deflection system of LEMS120 (Marhavilas et al. 2015). The vertical lines indicate, from left to right, the onset of the ESP event, the shock arrival, and the leading edge of the flux rope. (b) Energy spectrum of the intensities observed by ACE/EPAM/LEMS120 and LEMS30, observed by ERNE/LED, and extracted from the SEPEM RDS v2.0 at the time of the shock arrival.
In the text
thumbnail Fig. 2.

EUHFORIA time series at Earth (green), compared to 1 min OMNI data (blue). From top to bottom: speed, proton number density, magnetic field strength B, and Bx, By, Bz components in GSE coordinates. The vertical lines indicate the onset times of the observed (dotted line) and simulated (dashed line) flux rope.
In the text
thumbnail Fig. 3.

Simulated shock characteristics at Earth’s cobpoint. Panels a–c: give, respectively, the heliospheric radial distance, latitude, and longitude of Earth’s cobpoint in heliocentric Earth euatorial (HEEQ) coordinates as a function of the time elapsed after the CME injection at 0.1 au. Panels d–f: give, respectively, the shock speed, the shock compression ratio, and the shock angle at Earth’s cobpoint. The dashed vertical lines indicate the shock arrival time at Earth. Panel g: a snapshot of the solar wind radial speed at constant latitude θ = 4.5° close to the time when Earth established magnetic connection with the shock front. Earth is indicated by a black dot; it is magnetically connected to the cobpoint (magenta dot), which resides on the western flank of the CME.
In the text
thumbnail Fig. 4.

Parallel mean free path in the foreshock derived from in situ ion intensity measurements. Panel a: dots give the intensity-time profiles observed by EPAM for different energy channels, shifted downwards to avoid overlap. The straight lines give a time-averaged increase in the time intervals T1 = [0 h, 0.5 h] (blue) and T2 = [0.5 h, 10.3 h] (green), and the values above each line give the slope. Panel b: the dotted lines give the derived from the slopes presented in (a) for T1 (blue background) and T2 (green background). The solid line gives the f(s) (see Eq. (7)) used in the PARADISE simulations.
In the text
thumbnail Fig. 5.

Simulated (solid lines) and observed (dots) omnidirectional intensity-time profiles at Earth. The dots in the background give the ACE/EPAM/LEMS120 data (see also Fig. 1a). The dashed vertical lines indicate the arrival times of the shock (left) and the flux rope (right) in the EUHFORIA simulation.
In the text
thumbnail Fig. 6.

Simulation results for virtual observers at different longitudes. Panels a and b: give the omnidirectional intensity-time profiles (upper panel) together with solar wind speed (lower panel) for the W30 and E30 observers, respectively. Panel c: gives the peak intensities versus energy for different observers. Panel d: is a snapshot of the omnidirectional particle intensity plotted in grey shades on top of the solar wind speed at constant HEEQ latitude θ = 4.5°.
In the text
thumbnail Fig. 7.

Peak intensities scaled by r2 versus energy for virtual observers located along the Sun-Earth line.
In the text
thumbnail Fig. 8.

Peak intensity energy spectra obtained from different λ∥ power-law dependences on the ion rigidity. The orange, green, and red symbols give the results of PARADISE simulations using the indicated rigidity dependence. Blue symbols give the peak intensities measured by EPAM.
In the text
thumbnail Fig. A.1.

Illustration of the shock finder algorithm. Left panel: Snapshot of a EUHFORIA simulation containing a cone CME. Shown is the radial speed in the solar equatorial plane. The red line indicates the shock surface, as determined by the MCA. Right panel: 3D view of the shock surface computed by the MCA. The resolution has been reduced to clearly illustrate the triangle mesh.
In the text
thumbnail Fig. B.1.

Snapshot of the parallel mean free path of 88 keV protons (green shades) and ∇ ⋅ Vsw (red and blue) at constant HEEQ latitude θ = 4°. The parallel mean free path is only shown where λ∥ < 0.1 au.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.