© ESO, 2011

1. Introduction

Particle acceleration by coronal and interplanetary shocks driven by coronal mass ejections (CMEs) is widely accepted as the primary source of strong solar energetic particle (SEP) events. Particles scatter off plasma waves, cross the shock front repeatedly and gain energy on each crossing. In large events, turbulence sufficient for extended trapping can be generated by the streaming of the accelerated particles themselves, as plasma waves in the upstream are amplified by scattering particles. This diffusive shock acceleration mechanism, as presented by, e.g., Bell (1978), has recently been studied quantitatively as a model of the acceleration of SEPs by, e.g., Lee (2005), Vainio & Laitinen (2007, 2008) and Ng & Reames (2008). Other noteworthy studies of shock-accelerated SEP events include, e.g., Ng et al. (1999), Tylka et al. (2005), Tylka & Lee (2006) and Sandroos & Vainio (2007, 2009a,b).

Although the most abundant ion species in the solar wind, proton, is likely to dominate wave generation in shocks, minor ions are scattered off the same turbulence and accelerated as well. If the turbulence is proton-generated, the maximum rigidity obtained by the protons determines a low-wavenumber cutoff in the spectrum of plasma waves. It has been suggested that this would prevent any ions from being accelerated beyond the same rigidity in a quasi-parallel shock wave (e.g., Zank et al. 2007). This would then lead to a dependence of the maximum (non-relativistic) energy (per nucleon) of the ions on the charge to mass ratio of the form (Q/A)².

In this paper, we present the first simulations of particle acceleration in self-generated waves with minor ions included as particles contributing to the wave generation process. Instead of using a full spectrum of ion species in the solar wind, we limit ourselves to selected interesting populations, namely ³He²⁺, ⁴He²⁺, ¹⁶O⁶⁺ and ⁵⁶Fe¹⁴⁺. By simulating a propagating coronal shock for an extended period of time, we investigate the spectra and maximum attained energy of each particle population, as well as gauge the effect each population has on wave generation at different wavenumbers.

2. Model

In our simulation model, we numerically solve the particle transport equation through propagating representative particles using the guiding centre approximation. We approximate the quasi-linear theory by employing a pitch-angle independent resonance condition ${\mathit{f}}_{\mathrm{res}}\mathrm{=}{\mathit{f}}_{\mathrm{ci}}\frac{{\mathit{u}}_{\mathrm{sw}}\mathrm{+}{\mathit{v}}_{\mathrm{A}}}{\mathit{\gamma v}}$(1)where u_sw is the solar wind speed, v_A is the Alfvén speed, v is the particle speed, γ is the Lorentz factor, f_ci = (1/2π)q_iB/m_ic is the ion cyclotron frequency, B is the magnetic flux density and q_i = Qe and m_i = Am_p are the charge and mass of the ion in question, and e and m_p are the charge and mass of a proton. We scatter representative particles isotropically off plasma waves with scattering frequency $\mathit{\nu }\mathrm{=}{\mathit{\pi }}^{\mathrm{2}}{\mathit{f}}_{\mathrm{ci}}\frac{{\mathit{f}}_{\mathrm{res}}\mathit{P}\mathrm{\left(}{\mathit{f}}_{\mathrm{res}}\mathrm{\right)}}{{\mathit{B}}^{\mathrm{2}}}$(2)where P(f) is the wave power at frequency f. The initial wave power is scaled to give an ambient 100 keV proton mean free path of λ₀ = 1   R_⊙ at r₀ = 1.5   R_⊙, where R_⊙ is the solar radius.

Once particles have been swept up by the shock, they are propagated using a Monte Carlo method, focused along the mean magnetic field due to adiabatic invariance and traced in a superradially expanding flux tube with the solar wind speed inferred from mass conservation. For additional details of our simulation model, we refer the reader to Vainio & Laitinen (2007), Vainio & Laitinen (2008) and Battarbee et al. (2010).

2.1. Particle-shock interactions

As particles encounter the propagating parallel step-profile shock, they scatter off the compressed plasma and may return upstream with a momentum boost. The plasma compression ratio is solved from the Rankine-Hugoniot jump conditions. The analytical return probability from a shock encounter for an isotropic particle population is given as ${\mathit{P}}_{\mathrm{ret}}\mathrm{=}{\left(\frac{{\mathit{v}}_{\mathrm{w}}\mathrm{-}{\mathit{u}}_{\mathrm{2}}}{{\mathit{v}}_{\mathrm{w}}\mathrm{+}{\mathit{u}}_{\mathrm{2}}}\right)}^{\mathrm{2}}\mathit{,}$(3)where v_w is the speed of the particle in the downstream plasma frame and u₂ the downstream plasma speed in the shock frame. However, the downstream-transferred population is no longer isotropic, unless v_w ≫ u₂. Thus, following Vainio et al. (2000), we propagate and scatter particles in the downstream plasma frame up to a distance of 2λ behind the shock, where λ is the particle mean free path. At this distance, the particle has encountered enough scatterings to warrant the assumption of isotropy and its representative weight is multiplied by the isotropic return probability P_ret, sending it back towards the upstream from the distance of 2λ with a randomized shock-bound pitch-angle. Particles propagate and experience small-angle scatterings in the downstream until they either return to the shock front or their cumulative return probability drops below 0.1%. Simulation time is not advanced while the particles are in the downstream region. This corresponds to a situation where the downstream scattering is so intense that the mean residence time downstream is negligible compared to the time between subsequent shock encounters. Using this assumption, we do not propagate the shape of turbulence into the downstream.

2.2. Injected populations

In our simulation, particles swept up by a shock propagating through the solar corona are accelerated and traced up along a grid extending to a distance of 300   R_⊙. In addition to a proton population n_p(r) based on the solar wind density model of Cranmer & van Ballegooijen (2005), we inject fully ionized helium (³He²⁺ and ⁴He²⁺) and partially ionized heavier elements (⁵⁶Fe¹⁴⁺ and ¹⁶O⁶⁺) according to estimated solar wind abundance values. This results in H⁺-relative abundances for ⁴He²⁺, ¹⁶O⁶⁺, ⁵⁶Fe¹⁴⁺ and ³He²⁺ of 4.0 × 10^-2, 8.0 × 10^-4, 1.0 × 10^-4 and 1.6 × 10^-5, respectively.

As a large portion of the solar wind consists of thermal particles, and our strong step-like shock, with the assumed extremely intense downstream turbulence (see Sect. 2.1), injects an unrealistically large proportion of thermal particles, we model the solar wind as consisting of two kinetic populations. The majority, 99% of particles, represents a thermal core and is passed directly downstream. A minority, 1% of particles, is considered a partially suprathermal halo and follows a κ-distribution (Prested et al. 2008). This population, which has a continuous spectrum, is encountered by the shock in our simulation. The parameter κ receives values of 6...2 going from 1.5   R_⊙ to 3.0   R_⊙ respectively. It should be noted, though, that the composition of the seed population, and its dependence on distance from the sun, is not known and these parameters are arbitrary.

The average thermal speed ${\mathit{w}}_{\mathrm{0}}\mathrm{=}\sqrt{\mathrm{2}\mathit{T}{\mathit{k}}_{\mathrm{B}}\frac{\mathit{\kappa }\mathrm{-}\mathrm{1.5}}{\mathit{\kappa }{\mathit{m}}_{\mathrm{p}}}}$(4)is based on the radial temperature profile given in Cranmer & van Ballegooijen (2005), where T is the proton temperature and k_B is the Boltzmann constant. All minor ion populations are initialised using the same distribution function $\mathit{f}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}\frac{\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\Gamma }\mathrm{(}\mathit{\kappa }\mathrm{+}\mathrm{1}\mathrm{)}}{{\mathit{w}}_{\mathrm{0}}^{\mathrm{3}}{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{\kappa }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\Gamma }\mathrm{(}\mathit{\kappa }\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}}{\left[\mathrm{1}\mathrm{+}\frac{{\mathit{v}}^{\mathrm{2}}}{\mathit{\kappa }{\mathit{w}}_{\mathrm{0}}^{\mathrm{2}}}\right]}^{\mathrm{-}\mathit{\kappa }\mathrm{-}\mathrm{1}}\mathit{,}$(5)where Γ denotes the gamma function.

2.3. Wave evolution

The requirement for particle trapping in front of a shock is strong enough turbulence. This can be attained by wave amplification via the scattering of energetic particles upstream of the shock. A large shock-normal velocity results in large amounts of kinetic energy deposited into accelerated particles and thus larger amounts of energy deposited into upstream waves. Also, as particles reach higher velocities, they become resonant with waves of lower frequencies.

Fig. 1 Left: evolution of wave power spectra. Right: wave amplification factor Γw,i integrated over 640 seconds of simulation, at 2 × 10-5   R⊙ and 6 × 10-4   R⊙ from the shock, where Vs = 1500 km s-1.

In steady-state upstream solar wind, the evolution equation for a normalized wave power spectrum $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{P}\end{array}\mathrm{=}\mathrm{(}{\mathit{V}}^{\mathrm{2}}\mathit{/}\mathit{B}{\mathit{v}}_{\mathrm{A}}\mathrm{)}\mathit{P}$ can be written as $\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ \mathit{P}\end{array}}{\mathit{\partial t}}\mathrm{+}\mathit{V}\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ \mathit{P}\end{array}}{\mathit{\partial r}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{w}}\begin{array}{}\mathrm{˜}\\ \mathit{P}\end{array}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial f}}\left({\mathit{D}}_{\mathit{ff}}\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ \mathit{P}\end{array}}{\mathit{\partial f}}\right)\mathit{.}$(6)Here, P(r,f,t) is the Alfvén wave power spectrum as a function of radial distance (r), frequency (f) and time (t), V = u_sw + v_A is the group speed of the Alfvén waves and Γ_w is the wave growth rate. ${\mathit{D}}_{\mathit{ff}}\mathrm{=}\mathrm{(}\mathit{V}\mathit{/}{\mathit{r}}_{\mathrm{\oplus }}\mathrm{)}{\mathit{f}}^{\mathrm{8}\mathit{/}\mathrm{3}}{\mathit{f}}_{\mathrm{b}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}$ is an ad-hoc diffusion coefficient. This coefficient is chosen so that an unenhanced spectrum tends towards the form of Kolmogorov turbulence P ∝ f ^− 5/3 at r_⊕ = 1   AU above the breakpoint frequency f_b = 1   mHz, as suggested by observations (e.g., Horbury et al. 1996). Turbulence and diffusion magnitudes are normalized to result in an 100 kev proton having mean free paths of 1   R_⊙ at r = 1.5   R_⊙ and 54   R_⊙ at r = 1   AU.

Our simulation coordinates are attached to the propagating coronal shock, allowing for good numerical accuracy in the near-shock region. Wave amplification is calculated over time intervals of 1.6   ms at grid cell boundaries. We simulate the effect of diffusion using a Crank-Nicholson method, and advect wave power along the moving grid with a Lax-Wendroff scheme utilizing a Van Leer flux limiter.

For our simulations, we have extended the wave growth approximation of Vainio (2003) to accommodate for different particle masses and charges. Energy deposited into parallel-propagating Alfvén waves per particle scattering can be written as ΔE_w =  −v_ApΔμ, where p = γmv = γAm_pv is the particle momentum and Δμ is the change in particle pitch-angle. Integrating all particles of a given ion population i at a set position gives the rate of change of the wave energy density as $\frac{\mathrm{d}{\mathit{U}}_{\mathrm{w}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}\mathit{v} {\mathit{v}}_{\mathrm{A}}\mathit{\gamma A}{\mathit{m}}_{\mathrm{p}}\mathit{v}\frac{\mathrm{⟨}\mathrm{\Delta }\mathit{\mu }\mathrm{⟩}}{\mathrm{\Delta }\mathit{t}} \mathit{F}\mathrm{\left(}{r}\mathit{,}{v}\mathit{,t}\mathrm{\right)}\mathit{,}$(7)where F(r,v,t) is the particle distribution function.

Fig. 2 Particle populations for H+ and 56Fe14 along with the critical contour where the focusing velocity V/L exceeds the shock velocity. Results are shown for the simulations where Vs = 1500 km s-1, with colour contours at one magnitude intervals. A radial 3-cell boxcar smooth function has been applied.

The pitch-angle diffusion coeffient, according to quasilinear theory, is ${\mathit{D}}_{\mathit{\mu \mu }}\mathrm{=}\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}{\mathit{\omega }}_{\mathrm{ci}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\frac{\mathrm{\left|}{\mathit{k}}_{\mathrm{r}}\mathrm{\right|}\mathit{W}\mathrm{\left(}{\mathit{k}}_{\mathrm{r}}\mathrm{\right)}}{{\mathit{B}}^{\mathrm{2}}}$(8)where ω_ci = Qω_cp/Aγ = QeB/Am_pcγ is the angular ion cyclotron frequency, k_r is the resonant wavenumber and W(r,k,t)dk is the energy density of waves propagating parallel to the mean magnetic field with wavenumber in the range from k to k + dk. This can be written as ${\mathit{D}}_{\mathit{\mu \mu }}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\frac{\mathit{Q}{\mathit{\omega }}_{\mathrm{cp}}}{\mathit{A\gamma }{\mathit{m}}_{\mathrm{p}}{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{v}}_{\mathrm{A}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\mathrm{\left|}{\mathit{k}}_{\mathrm{r}}\mathrm{\right|}\mathit{W}\mathrm{\left(}{\mathit{k}}_{\mathrm{r}}\mathrm{\right)}$(9)which, because  ⟨ Δμ ⟩ /Δt = ∂D_μμ/∂μ, yields the wave growth rate $\begin{array}{ccc}{\mathrm{\Gamma }}_{\mathrm{w}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{W}}\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\mathit{t}}\\ & \mathrm{=}& \frac{\mathit{\pi }}{\mathrm{2}}\frac{\mathit{Q} {\mathit{\omega }}_{\mathrm{cp}}}{{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}\mathit{r}\mathrm{\right)} {\mathit{v}}_{\mathrm{A}}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}\mathit{v} \mathit{v} \mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}\mathrm{\left|}\mathit{k}\mathrm{\right|}\mathit{\delta }\left(\mathit{k}\mathrm{+}\frac{{\mathit{\omega }}_{\mathrm{ci}}}{\mathit{v\mu }}\right)\frac{\mathit{\partial F}}{\mathit{\partial \mu }}\mathrm{·}\end{array}$(10)As in Vainio (2003), we neglect the μ-dependence of the resonance condition and replace δ(k + ω_ci/vμ) with (1/2)δ(|k| − ω_ci/v). Now, using partial integration in μ, we get ${\mathrm{\Gamma }}_{\mathrm{w}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\frac{\mathit{Q} {\mathit{\omega }}_{\mathrm{cp}}}{{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}\mathit{r}\mathrm{\right)} {\mathit{v}}_{\mathrm{A}}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}\mathit{v} \mathit{v} \mathit{\mu }\mathrm{\left|}\mathit{k}\mathrm{\right|}\mathit{\delta }\left(\mathrm{\left|}\mathit{k}\mathrm{\right|}\mathrm{-}\frac{{\mathit{\omega }}_{\mathrm{ci}}}{\mathit{v}}\right) \mathit{F}$(11)which can be further represented as ${\mathrm{\Gamma }}_{\mathrm{w}\mathit{,}\mathrm{i}}\mathrm{\left(}{r}\mathit{,}{\mathit{v}}_{\mathrm{r}}\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}{\mathit{Q}}_{\mathrm{i}}{\mathit{\omega }}_{\mathrm{cp}}\frac{{\mathit{v}}_{\mathrm{r}}{\mathit{S}}_{\mathit{v}}\mathrm{\left(}{r}\mathit{,}{\mathit{v}}_{\mathrm{r}}\mathit{,t}\mathrm{\right)}}{{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{v}}_{\mathrm{A}}}$(12)where ${\mathit{S}}_{\mathit{v}}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathit{v}}^{\mathrm{2}}{{}^{\mathrm{\int }}}_{-1}^{\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{\mu } \mathit{v\mu F}\mathrm{\left(}{r}\mathit{,}{\mathit{v}}_{\mathrm{r}}\mathit{,\mu ,t}\mathrm{\right)}$ is the particle streaming per unit velocity in the frame of the Alfvén waves evaluated at the resonant particle speed v_r = ω_ci/|k|. We monitor particle streaming by counting whenever a particle crosses a grid cell boundary of the shock-attached tracking grid. To correct for the difference of particle streaming between the frame travelling with the shock and the frame of the Alfvén waves, we use an additional weighing factor of (v_wμ_w)/(v_sμ_s), where v_w, μ_w, v_s, and μ_s are the particle speed and pitch-angle in the Alfvén wave frame and shock-attached frame, respectively.

3. Results

We simulate three parallel, constant velocity shocks in the corona starting from 1.5   R_⊙. The shock-normal velocities V_s are 1250 km s^-1, 1500 km s^-1 and 1750 km s^-1. We use zero cross-helicity for the turbulence in the downstream, while in the upstream waves propagate away from the Sun in the plasma frame. Turbulence is tracked on a logarithmic grid reaching out to 300 R_⊙ in front of the shock.

3.1. Wave turbulence and its generation

The left panel of Fig. 1 shows the normalized wave power spectra (multiplied with the frequency f) in front of the shock after 160 and 640 s of simulation. In the right panel we show the wave amplification, ${\mathit{\eta }}_{\mathrm{i}}\mathrm{\left(}\mathit{f,z}\mathrm{\right)}\mathrm{=}{}^{\mathrm{\int }}\mathrm{\Gamma }_{\mathrm{w}\mathit{,}\mathrm{i}}\mathrm{(}\mathit{f,}{\mathit{r}}_{\mathrm{shock}}\mathrm{+}\mathit{z,t}\mathrm{)} \mathrm{d}\mathit{t}$, for each particle population i, integrated over the whole simulation time. The value z is the distance from the shock front and r_shock is the position of the shock. We display wave amplification for the V_s = 1500 km s^-1 run in both the measurement cell nearest to the shock (z = 2 × 10^-5   R_⊙) and further out (z = 6 × 10^-4   R_⊙). Although minor ions are much less abundant than protons, they have higher charges and are resonant with lower frequencies, allowing them to generate a significant amount of turbulence close to the shock. Such dominance of minor ions has been reported by Lee (1982) for Helium ions at low frequencies in relation to Earth’s bow shock. In our simulations, close to the shock, heavier ions surpass protons in wave generation at a narrow frequency range below 7 Hz. Between 7 Hz and approximately 70 Hz, ${}^{\mathrm{4}}\mathrm{He}^{\mathrm{2}\mathrm{+}}$ displays wave generation equal to  ~ 50% of that of H ⁺ .

For the V_s = 1250 km s^-1 shock, ${}^{\mathrm{4}}\mathrm{He}^{\mathrm{2}\mathrm{+}}$-powered wave generation is equal to H ⁺ -powered generation in the region below 100 Hz, with ${}^{\mathrm{4}}\mathrm{He}^{\mathrm{2}\mathrm{+}}$ dominating below 30 Hz. As the shock-normal velocity increases, the dominance of protons on turbulence amplification increases, with the ${}^{\mathrm{4}}\mathrm{He}^{\mathrm{2}\mathrm{+}}$ takeover moving to 3 Hz. In all cases, however, protons are responsible for the bulk of turbulence amplification, as abundant low-energy protons generate a great deal of turbulence above 100 Hz.

3.2. Accelerated particle populations

In Fig. 2 we demonstrate the evolution of two particle populations, H⁺ and ⁵⁶Fe¹⁴⁺. Particles within the expanding flux tube are efficiently accelerated up to a maximum energy at which the turbulence can no longer trap the particles, and instead the focusing effect of the diverging magnetic field allows them to escape.

The focused diffusion model of particle transport has been examined in detail by Kocharov et al. (1996). Here, we start with Parker’s equation, which in the fixed frame, upstream of our coronal shock reads $\frac{\mathit{\partial }{\mathit{f}}_{\mathrm{0}}}{\mathit{\partial t}}\mathrm{+}\mathit{V}\frac{\mathit{\partial }{\mathit{f}}_{\mathrm{0}}}{\mathit{\partial r}}\mathrm{-}\frac{\mathit{p}}{\mathrm{3}}\frac{\mathrm{1}}{\mathrm{𝒜}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}\mathrm{𝒜}\mathit{V}\mathrm{\right)}\frac{\mathit{\partial }{\mathit{f}}_{\mathrm{0}}}{\mathit{\partial p}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{𝒜}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left(\mathrm{𝒜}\mathit{D}\frac{\mathit{\partial }{\mathit{f}}_{\mathrm{0}}}{\mathit{\partial r}}\right)\mathit{,}$(13)where f₀ is the isotropic part of the distribution function, D = (1/3)λv is the spatial diffusion coefficient, λ = v/ν is the particle mean free path and is the flux-tube cross-sectional area related to the focusing length L by . Parker’s equation can be expressed, using the linear density , as a Fokker-Planck equation $\frac{\mathit{\partial n}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[\left(\mathit{V}\mathrm{+}\frac{\mathit{D}}{\mathit{L}}\right)\mathit{n}\right]\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial p}}\left[\frac{\mathit{p}}{\mathrm{3}}\left(\frac{\mathit{\partial V}}{\mathit{\partial r}}\mathrm{+}\frac{\mathit{V}}{\mathit{L}}\right) \mathit{n}\right]\mathrm{=}\frac{\mathit{\partial }}{\mathit{\partial r}}\left(\mathit{D}\frac{\mathit{\partial n}}{\mathit{\partial r}}\right)\mathit{,}$(14)which shows that the effect of focusing in the particle motion is two-fold: it contributes to the advection velocity by ṙ_foc = D/L and to the adiabatic energy changes by ṗ_foc =  −(p/3)   V/L. The addition to the advection velocity at large distances from the Sun is large, since there the waves have not yet grown to make D small. It is clear that particles will, on average, move away from the shock in the upstream region in areas where V + D/L > V_s. This facilitates the escape of particles from the shock to the upstream and the distance where V + D/L = V_s can be regarded as the boundary of the turbulent trapping region ahead of the shock. This boundary, displayed in Fig. 2 as a dashed line, outlines an escaping population further away from the shock. In addition, the energy at which the boundary intersects the shock surface is representative of the maximum energy that the particles can be accelerated to at a given time.

Fig. 3 Particle spectra after 80 s (left panel) and 640 s (right panel) of simulation, where Vs = 1500 km s-1. A power-law and an exponentional cut-off has been fitted to each population, ignoring the enhancement at the lowest energies.

In the latter stages of the simulation, we find decreased wave generation due to lower particle densities and thus less swept-up particles. This, along with wave diffusion, results in the turbulent trapping boundary at the shock moving to lower energies, which causes high energy particles to escape instead of experiencing further acceleration.

3.2.1. Spectral indices

The accelerated particle populations were integrated over the whole upstream. We found the spectral index α for the power-law part of the particle spectrum by fitting a line to a chosen section of the log-log representation of data points. Observing the one magnitude contours in Fig. 2, we see that the spatial distributions of iron and protons differ for both the escaping population and particle populations within the turbulent trapping boundary. For iron, this results in much harder particle spectra than what the steady-state model of Bell (1978) suggests. Spectra along with more complete parameter fits are exemplified in Fig. 3.

Fig. 4 Temporal evolution of spectral index α for H+ and 56Fe14+ populations, as fitted to the power-law section of the particle spectra.

As the simulation continues the spectrum for heavy ions softens. Figure 4 displays the evolution of the spectral indices α for H⁺ and ⁵⁶Fe¹⁴⁺. Protons, being the dominant particle population, do not exhibit a softening of the spectral index, whereas the effect is exceedingly prominent in the case of iron accelerated by a V_s = 1250   km   s^-1 shock.

3.2.2. Attained maximum energies

When attempting to gauge the maximum energy attained by a particle population, we attempted to fit a power law with an exponential cut-off to the simulated spectrum. First, we found the spectral index α for the power-law part of the particle spectrum, as in Sect. 3.2.1. We then used this as the basis for fitting an exponential cut-off energy E_c. The form used is ${\mathit{y}}_{\mathrm{i}}\mathrm{=}\mathit{C}{\mathit{E}}_{\mathrm{i}}^{\mathit{\alpha }}{\mathrm{e}}^{\mathrm{-}{\left(\frac{{\mathit{E}}_{\mathrm{i}}}{{\mathit{E}}_{\mathrm{c}}}\right)}^{\mathit{ϵ}}}$where C is a fitted constant and ϵ is chosen to fit the sharpness of the cut-off. In our work, we used values of ϵ = 4...2.5, with the value decreasing over simulation time. Figure 5 displays how the cut-off energy follows a ratio of mass-to-charge to the power of 1.5–1.6, where the exponent is significantly smaller than the theoretical estimation of 2.

Fig. 5 Ratios of particle cut-off energy to the charge/mass number as a log-log-plot.

4. Discussion and conclusions

Having studied the acceleration of multiple particle populations through self-generated turbulence with three different coronal shocks, we find that during early phases of the acceleration process, very hard, even flat spectra can be seen for high-mass ions. At all but the highest frequencies, the effect of minor ions on wave generation is non-negligible, especially in the region directly in front of the shock. As the shock-normal velocity increases, the deduced spectra become harder and the maximum energy attained increases. It is also seen that the maximum energy dependence (Q/A)^β does not exhibit β = 2, as suggested by Zank et al. (2007). Rather, the behaviour of cut-off energies is between β = 1.5 and β = 1.6.

At high energies, accelerated particles stream away from the shock due to focusing, as scattering particles supply insufficient wave amplification power to trap high-energy particles to the shock. This causes the turbulent trapping boundary to approach and intersect the shock at the maximum ion energy. As the shock-normal velocity increases, the particle spectra become harder and the energy at which the turbulent trapping boundary intersects the shock increases.

To gauge the effect of focusing versus trapping on particle energy, we note that out of V + D/L = V_s only the diffusion coefficient D ∝ v²/   (f_cif_resP(f_res)) depends on the particle species. Thus, at the maximum energy edge of the trapping boundary, D must be same for all particles. If at low frequencies the wave spectrum increases from ambient to amplified levels with a power law P ∝ f^b, as shown in Fig. 1, we find $\mathit{D}\mathrm{\propto }{\left(\frac{\mathit{Q}}{\mathit{A}}\right)}^{\mathrm{-}\mathit{b}\mathrm{-}\mathrm{2}} {\mathit{v}}^{\mathit{b}\mathrm{+}\mathrm{3}}\mathit{.}$(15)This gives the (non-relativistic) cutoff energy per nucleon a dependency of (Q/A)^{2(b + 2)/(b + 3)}. For a sharp cutoff, i.e. a purely rigidity-limited case, this results in a (Q/A)²-dependence for the cutoff energy, while a smoother transition results in a significantly weaker dependence. Due to the weak dependence of β on b, and the dynamic evolution of turbulence, care should be taken when inferring the turbulence spectrum shape from ion cutoff energies.

To examine particle acceleration dynamics, we can calculate the time τ_R required to accelerate a particle from an injection rigidity R₀ to a given rigidity R = p/q_i ∝ Av/Q, where R ≫ R₀. Assuming zero residence time in the downstream, we find ${\mathit{\tau }}_{\mathit{R}}\mathrm{=}\mathrm{𝒞}{\mathrm{\int }}_{{\mathit{R}}_{\mathrm{0}}}^{\mathit{R}}\frac{\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}}{{\mathit{R}}^{\mathrm{\prime }}}\mathit{D}\mathrm{\left(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathit{Q}}{\mathit{A}}\mathrm{𝒞}\left(\mathit{G}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{-}\mathit{G}\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}\right)\mathit{,}$(16)where is a constant and G(R) is a function of rigidity based on the shape of the wave power spectrum. If P(f) ∝ f^b where b is constant, the acceleration time to a given rigidity R is directly proportional to Q/A. Another item of interest is the time required to accelerate a heavy ion from injection speed v₀ to the maximum speed v_ion = (Q/A)^{(b + 2)/(b + 3)}v_p, where v_p is the proton speed at the turbulent trapping boundary. This gives $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{trap}}& \mathrm{=}& \mathrm{𝒞}{\mathrm{\int }}_{{\mathit{v}}_{\mathrm{0}}}^{{\mathit{v}}_{\mathrm{ion}}}{\left(\frac{\mathit{Q}}{\mathit{A}}\right)}^{\mathrm{-}\mathit{b}\mathrm{-}\mathrm{2}}\mathrm{(}{\mathit{v}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathit{b}\mathrm{+}\mathrm{2}}\mathrm{d}{\mathit{v}}^{\mathrm{\prime }}\\ & \mathrm{\propto }& \end{array}$(17)which, using previous assumptions for b, yields an acceleration time independent of the charge-to-mass ratio. Thus, it is clear that minor ions are accelerated to the maximum rigidity of protons much faster than the protons themselves, after which the

ions slowly continue to gain energy until they reach the turbulent trapping boundary. This results in minor ions gaining harder, even flat spectra, especially in early phases of the simulation. In latter phases of the simulation, the value of b increases at low frequencies, which leads to an increase in minor ion acceleration time and thus softer minor ion spectra.

At later stages of the simulations, wave amplification rates decay in line with the decay of injection efficiency. Thus, accelerated particles can stream away from the shock at lower energies, and further acceleration to higher energies ceases. Particles reaching the turbulent trapping boundary propagate in space, forming a plateau which is not consistent with Bell’s steady-state result.

In conclusion, diffusive shock acceleration of protons and minor ions cannot realistically be represented by a steady-state approximation, but instead requires numerical simulations to reveal the full dynamics of the acceleration process and the various particle populations.

Acknowledgments

The authors would like to thank the IT Centre for Science Ltd (CSC) for computational services and the Academy of Finland (AF) for financial support of projects 122041 and 121650. T.L. acknowledges support from the UK Science and Technology Facilities Council (standard grant ST/H002944/1).

References

All Figures

	Fig. 1 Left: evolution of wave power spectra. Right: wave amplification factor Γw,i integrated over 640 seconds of simulation, at 2 × 10-5   R⊙ and 6 × 10-4   R⊙ from the shock, where Vs = 1500 km s-1.
In the text

	Fig. 2 Particle populations for H+ and 56Fe14 along with the critical contour where the focusing velocity V/L exceeds the shock velocity. Results are shown for the simulations where Vs = 1500 km s-1, with colour contours at one magnitude intervals. A radial 3-cell boxcar smooth function has been applied.
In the text

	Fig. 3 Particle spectra after 80 s (left panel) and 640 s (right panel) of simulation, where Vs = 1500 km s-1. A power-law and an exponentional cut-off has been fitted to each population, ignoring the enhancement at the lowest energies.
In the text

	Fig. 4 Temporal evolution of spectral index α for H+ and 56Fe14+ populations, as fitted to the power-law section of the particle spectra.
In the text

	Fig. 5 Ratios of particle cut-off energy to the charge/mass number as a log-log-plot.
In the text