SEP ACCELERATION IN CME DRIVEN SHOCKS USING A HYBRID CODE

Coronal mass ejections (CMEs) are large ejections of solar material that periodically erupt from the Sun (Gopalswamy 2003; Forbes et al. 2006; Gopalswamy et al. 2008). As CMEs propagate out into interplanetary space, they can produce transient bursts of extremely energetic particles referred to as solar energetic particle (SEP) events (Sheeley Jr. et al. 1983; Kahler 2001; Gopalswamy 2003; Kahler 2004).

To be identified as an SEP, the flux of particles (protons, electrons with trace higher Z ions) with energies above 10 MeV, must be greater than 10 particle flux units (pfu = particles cm⁻² s⁻¹ str⁻¹; Gopalswamy 2003).

The energy spectra of the SEP populations vary considerably (Lin 1974; Hollebeke et al. 1975; Kallenrode et al. 1992) and show a dependence on the associated parameters of the originating CME (Park et al. 2012). Often the observed particle energies reach several hundred megaelectronvolts (Reames 1999), and even gigaelectronvolt energies (Ryan et al. 2000; Kim et al. 2013). SEP events can last from a period of several hours up to several days (Reames 1999). The combination of high flux and high penetrating particles means that SEP events intersecting the Earth and man-made technology in space present a significant "space weather" risk of damage and disruption to vulnerable systems (Feynman & Gabriel 2000) and to human tissue of astronauts (Wu et al. 2011). The SEPs from CME shock events tend to be the more extended in duration, or "gradual events," and have the harder energetic particle spectrum (Kahler 2001, 2004; Cliver et al. 2003) and therefore the most interest for space weather mitigation.

The characteristics of high-flux and high-energy spectra suggest a very effective acceleration mechanism associated with CME shock. While acknowledging that SEP-type events maybe associated with other phenomena (Tylka & Lee 2006; Rouillard et al. 2012), here, we consider the acceleration mechanism of CMEs propagating faster than ∼800 km s⁻¹ (Gopalswamy et al. 2008). At these propagation speeds, if the local plasma density n and magnetic field strength B encountered by the CME are such that the wave front is traveling super-Alfvénically, then it will create an interplanetary shock (Gopalswamy 2003; Park et al. 2012). Correlations between CME parameters of linear speed, angular extent, and relative location on the Sun have shown that the greatest predictor of SEP event occurrence and particle flux goes with increasing CME speed, 30% for 800 km s⁻¹ to 100% for CME speeds of 1800 km s⁻¹ (Gopalswamy et al. 2008; Hwang et al. 2010; Park et al. 2012). The presence of preceding CMEs has also been found to decide the peak solar energetic particle flux (Gopalswamy 2003; Gopalswamy et al. 2004), further indicating that the important parameter is the local Alfvén speed, which is being reduced ahead of the second CME front resulting in harder SEP spectra.

Current particle acceleration mechanisms from collisionless shocks (Sagdeev 1966) include shock drift acceleration and diffusive shock acceleration (DSA). The shock drift acceleration mechanism, dominant for perpendicular shocks, was originally studied by Dorman & Freidman (1959) and Schatzmann (1963); more recent reviews (Decker 1983; Toptygin 1983) estimate that the maximum energy gains attainable are ∼5 × the initial particle energies and depend on the magnetic field compression ratio due to the shock. The DSA mechanism (or first-order Fermi acceleration; Bell 1978; Blandford & Ostriker 1978) is thought to be responsible for the highest-energy particles observed at quasi-parallel shocks, thus being the preferred mechanism for cosmic-ray acceleration, and it is also being used to explain some features of particle spectra from SEP events. In DSA, particles crossing the shock front are accelerated by successive reflections downstream and upstream due to turbulence, potentially reaching very high energies. Fundamental theory on shock acceleration can be found in Toptygin (1983; Stone & Tsurutani 1985; Völk 1987; Jones & Ellison 1991).

Although turbulence exists in the solar wind for particle reflection to occur, its level is not sufficient to explain the production of megaelectronvolt and gigaelectronvolt particles during the time CMEs and interplanetary shocks take to reach the Earth (Sagdeev & Kennel 1991). Instead, turbulence is produced at the shock by waves arising at the shock front and propagating upstream (McKenzie & Völk 1982; Gordon et al. 1999; Zank et al. 2000; Ng et al. 2003; Li et al. 2003).

In this paper, we use a kinetic ion/fluid electron numerical simulation approach commonly known as a hybrid code (Dawson 1983; Fonseca et al. 2002) to study the propagation of a quasi-parallel CME shock in the solar wind environment. The code was originally developed to study the interaction of artificial plasmas released in the solar wind (Bingham et al. 1991; Gargate et al. 2008) and is now a massively parallel three-dimensional (3D) hybrid particle code, dHybrid (Gargaté et al. 2007). The code has been successfully used to investigate cosmic-ray acceleration at collisionless shocks (Gargaté & Spitkovsky 2012). The hybrid model uses massless fluid electrons and kinetic ions. The parallel implementation of this model allows the study of large regions of space (e.g., hundreds of ion gyro radius) over extended periods of time (e.g., tens of ion gyro periods), ideal for space plasma studies. Here, we consider a CME driving a fast magnetosonic shock, with shock parameters known to correlate well with SEP events (Park et al. 2012). In our simulations, a CME structure propagates at speeds of up to 1000 km s⁻¹ interacting with the slower solar wind. The interactions cause the formation of a large scale quasi-parallel shock structure due to the flowing CME. The acceleration mechanisms of high-energy particles are studied in this scenario. In the early acceleration phase, our results show that particles crossing the shock front accelerate perpendicularly to the shock front while maintaining their parallel velocity, supporting a surfatron-like acceleration model. The importance of this acceleration model as a means of providing a seed particle population for further acceleration is studied.

We explore the scenario of SEP acceleration and wave formation at CME-driven quasi-parallel shocks using a hybrid model; the shock evolution can be followed on the ion timescale, the ion acceleration at the shock front is correctly modeled, and the smaller electron timescales can be neglected by using an ideal fluid model for this species.

In comparison with MHD simulations, which do not capture kinetic effects and follow the evolution of a CME on a global scale, and over a time period relevant for the propagation of a CME in interplanetary space, hybrid simulations are localized in space, modeling a small part of the CME shock front and running over a time period relevant for the ion dynamics.

Results from dHybrid show the self-consistent formation of Alfvèn waves upstream of the shock, with turbulence building up due to wave breaking and strong particle acceleration. Energy gains of up to 110 times the maximum possible energy gain in one shock crossing are measured.

For the most accelerated particles, the observed energy gain is approximately quadratic in time during the simulation time frame, consistent with surfatron acceleration (Katsouleas & Dawson 1983; Üçer & Shapiro 2001; Lee et al. 1996), while for another less energetic set of particles, the energy scales with t^1/2 are consistent with DSA. The observed energy gain would allow for a typical solar wind proton to reach an energy of hundreds of megaelectronvolts in minutes. A thorough discussion about the observed acceleration mechanisms will be presented.

This paper is organized as follows. In the next section, we present the numerical model in detail, describe the simulation setup, and present the plasma parameters assumed. In the Results section, we investigate the wave formation, the wave–particle interaction mechanisms, and particle acceleration. We also include a simple single-particle theoretical model that clarifies how particles are accelerated in the upstream Alfvén waves, consistent with the observed simulation results. Finally, in the last section, we present the conclusions.

2.1. The Hybrid Model

Hybrid models, with kinetic ions and fluid electrons, are commonly used in many problems in plasma physics (Lipatov 2002). While MHD simulations are used to model CMEs globally, we use hybrid simulations to study shock properties locally, providing a new perspective over the problem of particle acceleration in gradual events.

The hybrid set of equations is derived neglecting the displacement current in Ampére's law, considering quasi-neutrality and calculating moments of the Vlasov equation for the electrons in order to obtain the generalized Ohms Law. In our implementation of the hybrid model in the massively parallel 3D code dHybrid (Gargaté et al. 2007), the effects of electron mass, resistivity, and electron pressure are not considered; thus, the electric field is simply given by $\boldsymbol {E} = -\boldsymbol {V_e} \times \boldsymbol {B}$, which can also be expressed as

$$\begin{equation} \boldsymbol {E} = -\boldsymbol {V_i} \times \boldsymbol {B} + \frac{1}{n} \left(\nabla \times \boldsymbol {B} \right) \times \boldsymbol {B}, \end{equation} \tag{ 1 }$$

where we have used $\boldsymbol {V_e} = - \boldsymbol {J}/(|e|n) \times \boldsymbol {V_i}$. $\boldsymbol {V_i} = ({1}/{n}) \int f_i \boldsymbol {v}\, d\boldsymbol {v}$ is the ion fluid velocity, and n is the electron/ion density. Normalized simulation units are used: time is normalized to $\omega _{ci0}^{-1}$, space is normalized to c/ω_pi0, charge is normalized to the proton charge |e|, and mass is normalized to the proton mass, where ω_ci0 is the ion cyclotron frequency, ω_pi0 is the ion plasma frequency, and c is the speed of light in vacuum. The magnetic field is advanced in time through Faraday's law $({\partial B}/{\partial t}) = - \nabla \times \boldsymbol {E}$, with $\boldsymbol {E}$ calculated from Equation (1).

In dHybrid, the ions are represented by finite-sized particles moving in a 3D simulation box and are treated as kinetic particles, with their velocity updated via the Lorentz force equation. The fields and fluid quantities, such as the density n and ion fluid velocity $\boldsymbol {V_i}$, are interpolated from the particles using quadratic splines (Decyk et al. 1996) and defined on a 3D regular grid. These quantities are then interpolated back to push the ions using quadratic splines, in a self-consistent manner. Equations are solved explicitly, based on a Boris pusher scheme to advance the particles (Boris et al. 1970) in the hybrid approach and on a two step LaxWendroff scheme to advance the magnetic field (Birdsall & Langdon 1985; Hockney & Eastwood 1994). Both schemes are second-order accurate in space and time and are space and time centered.

The present version of dHybrid uses the MPI framework as the foundation of the communication methods between processes and the Osiris visualization package (Fonseca et al. 2002) as the basis for all diagnostics. The three-dimensional simulation space is divided across processes; 1D, 2D, and 3D domain decompositions are possible, and dynamic load balancing is enabled, optimizing parallel efficiency by ensuring that the computational load is similar across processors. The code can simulate an arbitrary number of particle species with arbitrary charge to mass ratios, arbitrary initial thermal velocity and drift velocity distributions, as well as arbitrary spatial configurations. Periodic boundary conditions, open boundary conditions, and configurable particle injectors are used for the particles, and periodic boundary conditions are used for the fields.

Particle tracking techniques are also used in dHybrid and are of particular relevance for the problem at hand, allowing the study of the particle acceleration mechanisms in great detail.

Typically, a simulation is run twice: the first time all usual diagnostics can be analyzed (e.g., electric field, magnetic field, fluid phase spaces), and a special kind of diagnostics, the raw diagnostics, are produced. In these raw diagnostics, a sample of raw simulation particles are stored at given intervals, including the positions, velocities, and charge. A specific set of these particles is then chosen according to specified criteria (e.g., the hundred most energetic particles, a random sample of particles). The list of particles is then supplied as input for the second run, and all the positions, velocities, electric field, and magnetic field at the particle positions are stored for every iteration.

2.2. Simulation Setup

3.1. Wave Generation Upstream of the Shock

Figure 1 shows the charge density of the solar wind superimposed with the magnetic field lines for three distinct moments in time. The shock front is defined at x = 10 by the density jump between the solar wind upstream and the CME plasma downstream, as well as by the jump in the direction of the magnetic field. The solar wind plasma reflected at the shock is strongly modulated in the upstream region, and the magnetic field intensity does not suffer dramatic changes (δB/B ≪ 1), although the field direction varies slightly. The plasma density perturbations upstream of the shock become more turbulent with time, which is an indication of wave breaking that produces turbulence in the non-linear regime.

Zoom In Zoom Out Reset image size

By looking at the shock behavior, it is apparent that the solar wind plasma reflected at the shock front and propagating upstream drives electromagnetic waves in the upstream region. Figure 2 shows the time evolution of the v_xx phase space for the solar wind plasma, with the superimposed electric field intensity lineout along the shock direction, and shows that a wave is formed by the interaction of the two counter-streaming ion populations. The same oscillations are observed in the magnetic field (Figure 3), indicating the presence of an electromagnetic wave. Also, from Figure 3, it is seen that δE/E ≫ δB/B and that oscillations occur in the y and z components of both the electric field and the magnetic field, while there is a smaller amplitude oscillation of the x component of the electric field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Measuring the wavelength of the wave yields λ = 3 c/ω_pi0 so that k = 2.09 ω_pi0/c, and measuring the propagation velocity of the wave front yields v = 6.2v_A0, which is consistent with an Alfvén wave with frequency ω = 12.99 ω_ci0 in the simulation reference frame. This wave is actually supported by the reflected solar wind plasma, and in the reference frame moving with this plasma, the wave actually propagates in the −x direction with the Alfvén velocity of 3.5397 v_A0, yielding a frequency ω = 7.41 ω_ci0. The wave is then a rotating elliptically polarized Alfvén wave propagating along x with main components in the y and z directions and with a smaller (δE_x ≪ δE_y and δE_x ≪ δE_z) electrostatic component directed along the propagation axis x.

Wave formation due to counter-streaming super-Alfvénic ion populations is a known effect (McKenzie & Völk 1982; Gordon et al. 1999; Zank et al. 2000; Ng et al. 2003; Li et al. 2003). Different modes can be excited from MHD modes (Lee 1983) to kinetically driven Alfvén waves and to purely growing instabilities, relevant for cosmic-ray acceleration mechanisms (Lucek & Bell 2002; Bell 2004). Our results show an elliptically polarized Alfvén wave and also include an electrostatic component along the x direction. This component is due to the quasi-parallel magnetic field configuration that increases the complexity of the configuration, in comparison with the parallel magnetic fields usually assumed in the theoretical models. The wavelengths and growth rates are compatible with the instabilities described by Lucek & Bell (2002) and Bell (2004). For this instability, small wavelengths grow with time until a maximum wavelength is reached, beyond which the instability saturates. The quasi-linear MHD theory of the instability predicts a growth rate of $\gamma _{{\rm max}} = \zeta \, v^2_{sh}/(2\,v_A\,r_{{\rm Lsp1}})$ for the fastest growing wave number $k_{{\rm max}} = \gamma _{{\rm max}} v_A^{-1}$, with $\zeta = B_0\,j\,r_{{\rm Lsp1}}\,\rho ^{-1}\,v_{sh}^{-2}$ in the non-relativistic regime. The unstable wave vector range is $1< kr_{{\rm Lsp1}} < \zeta v_{sh}^2v_A^{-2}$. The instability works for parallel and quasi-parallel shocks, as in our case, and when one of the species is unmagnetized and the other species is magnetized. Here, v_sh is the relative velocity between the two plasma species, ρ is the mass density for the background (magnetized) species, j is the current density of the unmagnetized species, r_Lsp1 is the Larmor radius for the unmagnetized species, and v_A is the Alfvén velocity.

In our case, the ions reflected at the shock front get unmagnetized due to scattering, while the ions that are streaming toward the shock front are magnetized (they are streaming along a quasi-parallel magnetic field with a relatively low temperature).

The ion Larmor radius of r_Lsp1 ∼ 0.73 c/ω_pi0 can be measured directly in the simulation, but the density ratio n_sp1/n_sp0 of the two counter-streaming ion populations, controlling the parameter ζ through the current j and mass density ρ, varies strongly during the simulation.

This is not accounted for in the theoretical model, which assumes a constant current driving the instability, an isotropic, or power-law, particle distribution for the unmagnetized species, and propagation parallel to the magnetic field (Bell 2004). Since the current driving the instability is not constant, the propagation is quasi-parallel, and the particle distribution, at the spatial lengths considered, is not isotropic (see Figures 2 and 4). Only an order of magnitude estimation can be done for the theoretical values of γ_max and k_max. For n_sp1/n_sp0 ∼ 0.08, as in the early stages of the simulation, a growth time of $\gamma _{{\rm max}} \sim 7\,\omega _{ci0}^{-1}$ and a wave number of k_max ∼ 2 ω_pi0/c can be estimated.

Zoom In Zoom Out Reset image size

A significant number of particles reflected at the shock front are seen to interact with the previously formed waves and accelerate. From the inspection of Figures 2 and 4, showing the v_xx, v_yy, and v_zz phase spaces, it is seen that the energy gain is mostly in the y and z directions, that is, in the directions perpendicular to the shock propagation.

Another interesting observation is that a part of the particle population that is streaming in the +x direction is being reflected back to the shock at x ∼ 15 c/ω_pi0 for later times, visible in Figure 2(c), where some of the particles in the upper branch of the phase-space plot have negative velocities at that point. Also, the particles with the greatest perpendicular velocities up to $t = 2.53\, \omega _{ci0}^{-1}$ (c) are moving away from the shock front.

At later times, $t \sim 4 \, \omega _{ci0}^{-1}$, in the simulation, there is indication of a new group of particles gaining energy in the region of the shock front around x ∼ 10 c/ω_pi0.

Looking directly at the kinetics of the most energetic particles provides insights on the physical processes that dominate particle acceleration. The particle tracks for the top 80 most energetic particles in the simulation reveal two distinct groups of particles. Figure 5 shows the five most energetic particles from each of these two groups. Particles from group 1 are accelerated very early in time, and the acceleration occurs in two distinct phases: a strong energy gain around $1.3 < t < 2.0 \omega _{ci0}^{-1}$ and a weaker increase in energy from $t \sim 2 \omega _{ci0}^{-1}$ onward. These are the particles that move away from the shock front, reaching energies up to 22.5 times their initial energy (dashed lines in Figure 5), which are in a zone dominated by the Alfvén wave and causing the energetic particle population seen in Figures 4(b) and (c).

Zoom In Zoom Out Reset image size

Particles from group 2 (Figure 5) start gaining energy only later, around $t \sim 2\, \omega _{ci0}^{-1}$. This kind of behavior for interplanetary shocks, with two distinct phases, is predicted by Lee (1983), who presents a model for turbulence enhancement due to counter-streaming ion populations, generation of waves in the upstream media, and DSA acceleration of particles due to this turbulence.

At $t \sim 2\omega _{ci0}^{-1}$, the fields still preserve a wave-like structure in regions away from the shock, as can be seen in the electric field lineout, 17 < x < 25c/ω_pi0, in Figure 2(c). Spatial regions near the shock front, however, start to exhibit turbulence, mostly visible in the density plot of Figure 1(c), upstream of the shock front. This turbulence intensifies with time and is the reason why particles in group 2 are trapped in the shock front and start to shock drift, gaining energies up to 50 times their initial energy in $t \sim 2\omega _{ci0}^{-1} \sim 7.5$s. The total energy gain for these particles is ∼110 times the maximum energy that a particle could gain by simply crossing the shock front once.

With an energy increase in time that is approximately quadratic (Figure 5) it would be possible for a typical 1 keV proton to reach an energy of ∼200 MeV in around 10 minutes time, if the energy gain could be sustained for that period of time.

Figure 6 shows the velocity components of the most accelerated particle in group 1 (dashed) and the most accelerated particle in group 2 (full line). The fundamental difference between the two particles is that particle 1 (from group 1) has a positive v_x velocity and traverses the most efficient acceleration zone, situated in front of the shock, very quickly, gaining most of its energy in a time interval t ∼ 0.7. The velocity increase in this period is, however, approximately linear (1.3 < t < 2; Figure 6). Particle 2 (from group 2) exhibits the same behavior: the mean velocity increases linearly in the time interval 2 < t < 4, although the acceleration is more efficient, due to the initial v_x ∼ 0 velocity. The result is that particle 2 drifts along the shock front in the upstream region, while particle 1 follows the wave propagating away from the shock.

Zoom In Zoom Out Reset image size

The velocity profile increase seen in Figure 6 in the acceleration phases is consistent with the picture of a particle trapped in a circularly polarized Alfvén wave. In the simplest form, we consider a zero-order magnetic field parallel to the shock normal, B_x, and an Alfvén wave with amplitude A₀ and components in $\boldsymbol {e_y}$ and $\boldsymbol {e_z}$, with an electric field component $\boldsymbol {E} = A_0 \,[\, \cos \,(kx - \omega t)\,\boldsymbol {e_y} - \sin \,(k x - \omega t)\,\boldsymbol {e_z}\,]$ and a magnetic field component $\boldsymbol {B} = A_0\, [\, \sin \,(k x - \omega t)\,\boldsymbol {e_y} + \cos \,(k x - \omega t)\,\boldsymbol {e_z}\,]$. An ion can be trapped in this wave if, in zeroth-order, v_x = ω/k + ω_cx with ω_cx = q B/m; solving the single particle motion using the Lorentz force equation and using the above trapping condition yields v_y(t) = K₁cos (ω_cx t) t and v_z(t) = −K₁sin  (ω_cx t) t with K₁ constant.

The above picture recovers the behavior seen in Figure 6 for the perpendicular velocity components v_y and v_z and does not explain any energy gain along the magnetic field in the x direction. If we refer again to Figure 3, we can see that the simple assumptions made are oversimplistic, and instead, an electrostatic wave component would have to be considered, along with different electric field and magnetic field amplitudes, and also finite values for the static magnetic field components. The model describes the main qualitative features of the acceleration well, while the quantitative behavior in the much more realistic simulation scenario is more complex.

We have presented 2D hybrid simulation results of a quasi-parallel shock, with realistic shock parameters. In the shock reference frame used in the simulation, the solar wind plasma flows along the quasi-parallel magnetic field, hits the coronal mass ejection plasma, and is scattered. The upstream population of scattered ions induces the formation of an electromagnetic Alfvén wave. In a completely self-consistent picture, the Alfvén waves create turbulence upstream and accelerate a significant population of ions.

The results presented are qualitatively different from those provided by the usual MHD simulation techniques: the shock propagation is followed on a different timescale, relevant for the ion dynamics, and shorter than the typical MHD simulation timescale.

Also, the simulation is localized in space in comparison with MHD simulations than can model the global behavior of a CME. The detailed spatial and temporal resolution attained results in a much more complete physical picture, in which there are electromagnetic waves propagating due to the counter-streaming ion populations in the upstream region and in which particles are accelerated in two distinct phases.

The shock propagates for a time interval $T = 4\,\omega _{ci0}^{-1} \sim 15\, {\rm s}$, and the most accelerated particles start gaining energy at $t = 2\, \omega _{ci0}^{-1}$ from an initial thermal distribution. The energy gain is approximately quadratic in time up to 50 times the initial energy, meaning that if part of these ions could be trapped for periods of time of ∼10 minutes, the energy gain could lead to particles with ∼200 MeV, consistent with observations.

The crucial question of whether the energy gain is sustainable for long periods of time requires further investigation. Shock-drift theory dictates that, for a particle in a perpendicular shock, the energy gain is proportional to the magnetic field compression ratio which, depending on the shock strength, means an energy gain of up to five times the initial energy of a particle. In this case, due to the wave structure present at the shock front, the observed energy gain is much greater, going up to 50 times the initial energy of a particle that is shock drifting. Particles will always drift away from the shock front and that means that for further acceleration a diffusive shock mechanism, driving the particles back to the shock, is necessary.

Our simulations provide evidence of particles being reflected back to the shock (see Figure 5), suggesting that DSA actually occurs in the configuration considered and indicating that the hybrid simulation model is capable of correctly modeling the mechanism. While a better characterization of the DSA mechanism can be done because particle kinetics can be directly observed, a complete understanding of the mechanism involves modeling the shock propagation for times longer than those presented in this paper.

This work was supported by the European Research Council (ERC-2010-AdG grant 267841) and by Fundação para a Ciência e a Tecnologia (FCT/Portugal) under grant SFRH/BD/17750/2004. The simulations presented in this paper were produced using the IST Cluster (IST/Portugal). The authors also thank the U.K. STFC Center for Fundamental Physics for its support.