© ESO, 2013

1. Introduction

Coronal mass ejections (CME) are believed to be an important source of solar energetic particles (SEP). The acceleration mechanism taking place is the diffusive shock acceleration process (Reames 1999). One cornerstone in the theory behind this mechanism is understanding particle scattering of SEP in the background plasma. The general approach used for this purpose is the quasilinear theory (QLT) (Jokipii 1966).

The CME may be represented as a system of a thermal plasma with a turbulent energy spectrum and a nonthermal proton and electron component. These constituents are interacting with each other in a nonlinear way. In order to understand this system, it is necessary to model wave excitation by particle beams, turbulent energy transfer, and the transport of high-energy particles. Models have been developed by, e.g. Vainio & Laitinen (2007) and Ng & Reames (2008). While these models are capable of describing the interaction of particles and waves, the details of the specific processes are implemented in a simple way, basically relying on QLT or its simple extensions as the correct description of particle scattering and ignoring oblique or perpendicular wave modes completely. Lange & Spanier (2012) did numerical simulations of the turbulent cascading process of energy injected via proton beams. The basic concept of this paper is to take first steps towards a consistent description of particle transport within these plasma configurations and investigate how well the QLT describes scattering in such turbulence.

The general problem in the derivation of transport coefficients is a gap between theoretical and numerical approaches. Numerical simulations have used artificial turbulence spectra (Qin et al. 2002), while analytical calculations use any possible spectrum. However, the applicability of these calculations is unclear. This work combines the previous turbulence simulations by Lange & Spanier (2012) with a particle tracking algorithm derived from Spanier & Wisniewski (2011). The results are compared with QLT calculations using an identical input spectrum.

This approach allows a comparison of QLT and numerical simulations. The latter are assumed to be valid for a wide range of turbulent spectra and particle energies. Therefore, a determination of the validity of QLT is possible. Novel visualization and comparison methods are used to probe nonlinear and time-dependent scattering processes. Thus, a more fundamental understanding of the scattering process is achieved.

2. Theory

2.1. Numerical model

The region of the heliosphere we are interested in is within the weak turbulence regime, with magnetic field fluctuations defined as $\mathit{\delta }{B}\mathrm{\equiv }{B}\mathrm{-}{{B}}_{\mathrm{0}}\mathit{,}$(1)where ⟨ δB ⟩  = 0 is the mean value of the fluctuations, which leads to ⟨ B ⟩  = B₀. In that sense, weak turbulence is connected to small fluctuation amplitudes compared to the magnetic background field B₀. In terms of wave turbulence, this means that the linear Alfvén wave timescale is shorter than the intrinsically nonlinear time. It is observed that the solar wind magnetic fluctuations decrease as δB² ∝ r^-3, while the background field decreases as ${\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\propto }{\mathit{r}}^{-4}$(Bavassano et al. 1982; Bruno & Carbone 2005). Consequently the δB/B₀ ratio within the heliosphere is increasing with distance from the Sun (Hollweg et al. 2010).

The magnetic background field B₀ is defined towards the z-direction within our simulations and hence is also denoted as B₀e_z. The parallel direction is, therefore, defined as the z-direction and the x- and y-direction are the perpendicular directions. For symmetry reasons, there will be no further distinction between the two perpendicular directions, and all plots show values averaged over the azimuthal angle in cylindrical coordinates of the x–y plane.

High Reynolds numbers in combination with massive energy injection, as seen in, e.g. the solar wind, are strong indicators of a highly turbulent state. This means that most parts of the heliosphere are dominated by turbulent evolution and consequently energy cascading to smaller spatial scales. In situ measurements of the energy spectrum (Tu et al. 1989) are in agreement with this fact.

To simulate conditions within the turbulent heliospheric plasma with particle transport and scattering, the research group at the University of Würzburg has developed a hybrid simulation code, Gismo. It is divided into two parts, the magnetohydrodynamic (MHD) algorithm Gismo-MHD and the kinetic particle transport simulation tool Gismo-Particles.

Gismo-MHD is an incompressible pseudospectral MHD code that is fully parallelised and capable of efficiently using massive computing clusters. For a more detailed description we refer to Lange & Spanier (2012). A brief overview is presented below. The basis of the simulation software is to solve the following set of incompressible MHD equations: $\begin{array}{ccc}\frac{\mathit{\partial }{u}}{\mathit{\partial t}}& \mathrm{=}& {b}\mathrm{·}\mathrm{\nabla }{b}\mathrm{-}{u}\mathrm{·}\mathrm{\nabla }{u}\mathrm{-}\mathrm{\nabla }\mathit{P}\mathrm{+}{\mathit{\nu }}_{\mathit{v}}{\mathrm{\nabla }}^{\mathrm{2}\mathit{h}}{u}\\ \frac{\mathit{\partial }{b}}{\mathit{\partial t}}& \mathrm{=}& {b}\mathrm{·}\mathrm{\nabla }{u}\mathrm{-}{u}\mathrm{·}\mathrm{\nabla }{b}\mathrm{+}\mathit{\eta }{\mathrm{\nabla }}^{\mathrm{2}\mathit{h}}{b}\\ \mathrm{\nabla }\mathrm{·}{u}& \mathrm{=}& \mathrm{0}\\ \mathrm{\nabla }\mathrm{·}{b}& \mathrm{=}& \end{array}$(2)where ${b}\mathrm{=}{B}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \varrho }}$ is the normalised magnetic field, u is the fluid velocity, and ϱ is the constant mass density. The dissipative terms due to viscous and Ohmic dissipation are denoted by ν_v and η. In the following we introduce the parameter ν as a global diffusivity with η = ν_v ≡ ν. The approach is valid for magnetic Prandtl number of order unity, which is applicable within the regime of Alfvén wave turbulence where an equipartition between magnetic and kinetic energy can be assumed (Bershadskii 2002; Bigot et al. 2008). The artificial increase of the dissipation by the hyperdiffusivity parameter h is often used in spectral methods. The pressure term ∇P fulfils the closure condition for incompressibility (Maron & Goldreich 2001): $\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{P}\mathrm{=}\mathrm{\nabla }{b}\mathrm{:\nabla }{b}\mathrm{-}\mathrm{\nabla }{u}\mathrm{:\nabla }{u}\mathit{.}& & \end{array}$(3)In the incompressible regime of a magnetised plasma, the MHD turbulence consists of only two types of waves, which transport energy along the parallel direction: the so-called pseudo and shear Alfvén waves. The former are the incompressible limit of slow magnetosonic waves and play a minor role within anisotropic turbulence (Maron & Goldreich 2001). The pseudo Alfvén waves polarisation vector is in the plane spanned by the wave vector k and B₀. The shear waves are transversal modes with a polarisation vector perpendicular to the k − B₀ plane. They are circularly polarised for parallel propagating waves. Both species exhibit the dispersion relation ω² = (v_Ak_∥)². We note that the shear mode seems to be dominant because pseudo waves are heavily damped by the Barnes damping process within weakly turbulent regimes (Sridhar & Goldreich 1994). The damping weakens in strong turbulence, but according to Goldreich & Sridhar (1995), the wave generation of pseudo Alfvénic wave modes is only possible via three-wave interactions by two shear wave modes. Barnes damping is important for high-β plasmas. Since this is not fulfilled for the solar corona, the role of pseudo Alfvén waves should not be ignored. However, the solar wind and especially its fast component is strongly dominated by Alfvén waves (Bruno & Carbone 2005; Howes et al. 2012) and therefore incompressible MHD is well applicable.

Because of this model, a suitable description of the system is achieved by the use of Alfvénic waves moving either forwards or backwards. Therefore the Elsässer variables (Elsässer 1950) are introduced: $\begin{array}{ccc}{{w}}^{\mathrm{-}}& \mathrm{=}& {v}\mathrm{+}{b}\mathrm{-}{\mathit{v}}_{\mathit{A}}{{e}}_{\mathrm{\parallel }}\\ {{w}}^{\mathrm{+}}& \mathrm{=}& {v}\mathrm{-}{b}\mathrm{+}{\mathit{v}}_{\mathit{A}}{{e}}_{\mathrm{\parallel }}\mathit{.}\end{array}$(4)Applying this definition to Eqs. (2) and transformation to the Fourier space yields $\begin{array}{ccc}\left({\mathit{\partial }}_{\mathit{t}}\mathrm{-}{\mathit{v}}_{\mathit{A}}{\mathit{k}}_{\mathit{z}}\right)\stackrel{}{{{\mathit{w}}_{\mathit{\alpha }}^{\mathrm{-}}}_{\mathrm{˜}}}& \mathrm{=}& \frac{\mathrm{i}}{\mathrm{2}}\frac{{\mathit{k}}_{\mathit{\alpha }}{\mathit{k}}_{\mathit{\beta }}{\mathit{k}}_{\mathit{\gamma }}}{{\mathit{k}}^{\mathrm{2}}}\left(\begin{array}{c}]\\ {\mathit{w}}_{\mathit{\beta }}^{\mathrm{+}}{\mathit{w}}_{\mathit{\gamma }}^{\mathrm{-}}\end{array}\mathrm{+}\begin{array}{c}]\\ {\mathit{w}}_{\mathit{\beta }}^{\mathrm{-}}{\mathit{w}}_{\mathit{\gamma }}^{\mathrm{+}}\end{array}\right)\\ & & \mathrm{-}\mathrm{i}{\mathit{k}}_{\mathit{\beta }}\begin{array}{c}]\\ {\mathit{w}}_{\mathit{\alpha }}^{\mathrm{-}}{\mathit{w}}_{\mathit{\beta }}^{\mathrm{+}}\end{array}\mathrm{-}\frac{\mathit{\nu }}{\mathrm{2}}{\mathit{k}}^{\mathrm{2}\mathit{h}}\stackrel{}{{{\mathit{w}}_{\mathit{\alpha }}^{\mathrm{-}}}_{\mathrm{˜}}}\\ \left({\mathit{\partial }}_{\mathit{t}}\mathrm{+}{\mathit{v}}_{\mathit{A}}{\mathit{k}}_{\mathit{z}}\right)\stackrel{}{{{\mathit{w}}_{\mathit{\alpha }}^{\mathrm{+}}}_{\mathrm{˜}}}& \mathrm{=}& \frac{\mathrm{i}}{\mathrm{2}}\frac{{\mathit{k}}_{\mathit{\alpha }}{\mathit{k}}_{\mathit{\beta }}{\mathit{k}}_{\mathit{\gamma }}}{{\mathit{k}}^{\mathrm{2}}}\left(\begin{array}{c}]\\ {\mathit{w}}_{\mathit{\beta }}^{\mathrm{+}}{\mathit{w}}_{\mathit{\gamma }}^{\mathrm{-}}\end{array}\mathrm{+}\begin{array}{c}]\\ {\mathit{w}}_{\mathit{\beta }}^{\mathrm{-}}{\mathit{w}}_{\mathit{\gamma }}^{\mathrm{+}}\end{array}\right)\\ & & \end{array}$(5)where the tilde-notation represents quantities in the Fourier space. This is the final set of equations which is iterated by Gismo-MHD.

Obviously, the nonlinearities of Eqs. (5) that describe the turbulent behaviour of the MHD plasma cannot be solved in Fourier space in a straightforward fashion. Hence the main numerical load is the transfomation between real and wavenumber space for each iterative step. For this purpose we used the P3DFFT algorithm, which is an MPI-parallelised fast Fourier transfomation based on FFTW3 (Pekurovsky 2011).

A basic problem of spectral methods that use discrete Fourier transformation is the aliasing effect. Because of discrete sampling in the wavenumber space, high k-values exhibit errors that depend on the structure of the real space fields. Therefore we used zero padding, which is also referred to as Orszag’s 2/3 rule, i.e. 2/3 of the wavenumbers below the Nyquist frequency have to be truncated to achieve maximum anti-aliasing, hence reducing the Fourier space resolution to 1/3 of the original wavenumber range (Orszag 1971). This process is repeated for each step, immediately before calculating the nonlinearities and, accordingly, calculating the right-hand side (RHS) of the MHD equations. Consequently, the change in the antialiasing-range of one MHD step is physically correct, but not the long-term evolution.

The code Gismo is capable of using different foward-in-time schemes, namely, Euler and Runge-Kutta (RK) second as well as fourth order. All the simulations in this paper have been peformed using RK-4.

The Alfvén wave generation mechanism by SEPs is not investigated in detail here. The driving mechanism assumed for our simulations is the streaming instability. The estimation of the wave growth rate is described in Vainio (2003). The streaming instability is caused by energetic proton scattering off the coronal Alfvén waves. During the the scattering process the particle changes its pitch-angle cosine by Δμ, while its wave-frame energy remains constant. Thus the particles’ energy in the plasma frame is changed by v_ApΔμ, where p is the wave-frame particle momentum. Accordingly, the Alfvén wave energy in the plasma frame changes due to energy conservation. Another important instability in the solar corona is the electrostatic instability, which is caused by an electron current and streaming ions. Ion acoustic waves would be generated by this process. However, for the growth rate of these modes a sufficiently high ratio T_e/T_i ≫ 1 of the electron and ion temperatures is crucial. Observations and simulations in the vicinity of three solar radii indicate temperature ratios of the order of unity (Landi 2007; Jin et al. 2012). In this regime the ion acoustic waves are also efficiently suppressed by Landau damping. For further reading about the streaming instability we refer to Gary (1993).

The numerical tool for the kinetic simulation of charged test particles is Gismo-Particles. Like Gismo-MHD it is fully parallelised and capable of using massive computing clusters. Its core is the calculation of the relativistic Lorentz force $\begin{array}{ccc}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\mathit{\gamma }\hspace{0.17em}{v}\mathrm{=}\frac{\mathit{q}}{\mathit{mc}}\left[\hspace{0.17em}\mathit{c}{E}\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{v}\mathrm{\times }{B}\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\hspace{0.17em}\right]\mathit{,}& & \end{array}$(6)which acts on the particles at position x through the electric E and magnetic fields B of the plasma. The particle velocities are denoted by v, c is the speed of light and γ represents the Lorentz factor. The Gismo-Particles tool is capable of simulating test particles using physical masses and charges of electrons and ions. The test particles do not influence the background plasma, i.e. self-consistent backreactions to the plasma are neglected. A charged particle within a magnetic field will perform a gyromotion with the frequency $\begin{array}{ccc}\mathrm{\Omega }\mathrm{=}\frac{\mathit{ZeB}}{\mathit{\gamma mc}}& & \end{array}$(7)and the Larmor radius $\begin{array}{ccc}{\mathit{r}}_{\mathrm{L}}\mathrm{=}\frac{\mathit{v}}{\mathrm{\Omega }}\mathrm{·}& & \end{array}$(8)A useful numerical approach for solving Eq. (6) for gyrating particles is the implicit scheme of the Borispush. The basic idea was given by Boris (1970), where the iterations of the Lorentz force are separated in two partial steps. First, the particles are accelerated by the electric field within a half time step. Second, the gyromotion of the particles, which is caused by the magnetic field, is calculated. After that, the electric fields acts again for another half time step to complete the iteration. This approach leads to a discretisation of the Lorentz force with the following set of equations: $\begin{array}{ccc}{{v}}^{\mathit{t}\mathrm{-}\mathrm{\Delta }\mathit{t}\mathit{/}\mathrm{2}}& \mathrm{=}& {{v}}^{\mathrm{-}}\mathrm{-}\frac{\mathit{q}{E}\mathrm{\Delta }{\mathit{t}}^{\mathrm{num}}}{\mathrm{2}\mathit{m}}\\ \frac{{{v}}^{\mathrm{+}}\mathrm{-}{{v}}^{\mathrm{-}}}{\mathrm{\Delta }\mathit{t}}& \mathrm{=}& \\ {{v}}^{\mathit{t}\mathrm{+}\mathrm{\Delta }\mathit{t}\mathit{/}\mathrm{2}}& \mathrm{=}& {{v}}^{\mathrm{+}}\mathrm{+}\frac{\mathit{q}{E}\mathrm{\Delta }{\mathit{t}}^{\mathrm{num}}}{\mathrm{2}\mathit{m}}\mathit{,}\end{array}$where quantities with Δt^num/2 denote the discrete half-time steps. To solve this set of equations with respect to v⁺, the following steps were used: $\begin{array}{ccc}{{h}}_{\mathrm{1}}& \mathrm{=}& \\ {{v}}^{\mathrm{\prime }}& \mathrm{=}& {{v}}^{\mathrm{-}}\mathrm{+}{{v}}^{\mathrm{-}}\mathrm{\times }{{h}}_{\mathrm{1}}\\ {{h}}_{\mathrm{2}}& \mathrm{=}& \frac{\mathrm{2}{{h}}_{\mathrm{1}}}{\mathrm{1}\mathrm{+}{{h}}_{\mathrm{1}}\mathrm{·}{{h}}_{\mathrm{1}}}\\ {{v}}^{\mathrm{+}}& \mathrm{=}& \end{array}$where the auxiliary vectors h were introduced to calculate v⁺(Birdsall & Langdon 2005). By using these equations, the new velocity and location for each particle are then $\begin{array}{ccc}{{v}}^{\mathit{t}\mathrm{+}\mathrm{\Delta }\mathit{t}\mathit{/}\mathrm{2}}& \mathrm{=}& {{v}}^{\mathrm{+}}\mathrm{+}\frac{\mathit{q}{E}\mathrm{\Delta }\mathit{t}}{\mathrm{2}\mathit{m}}\\ {{x}}^{\mathit{t}\mathrm{+}\mathrm{\Delta }\mathit{t}}& \mathrm{=}& {{x}}^{\mathit{t}}\mathrm{+}\frac{{{v}}^{\mathit{t}\mathrm{+}\mathrm{\Delta }\mathit{t}\mathit{/}\mathrm{2}}}{\mathit{\gamma }\mathrm{\left(}{{v}}^{\mathit{t}\mathrm{+}\mathrm{\Delta }\mathit{t}\mathit{/}\mathrm{2}}\mathrm{\right)}}\mathrm{\Delta }{\mathit{t}}^{\mathrm{num}}\mathit{.}\end{array}$The advantage of the Borispush is the very high numerical stability. The particles are assumed to undergo gyromotions, hence the particle orbits themselves are stable for an arbitrary time discretisation. Even in the limit of Δt^num ≫ Ω^-1, the particle orbit is stable but converges to an adiabatic drift motion. The limitation of this method is the correct resolution of the Larmor radius r_L. If the timestep chosen is too large, this leads to a big deviation from the analytical r_L.

To provide a correct resolution of the gyromotion, Gismo-Particles performs a couple of single-particle gyrations with the used parameters during each start-up phase. The Larmor radius is then measured using a three point-circle approximation. If the deviation to the analytical value is above a given accuracy threshold, the gyrosimulation is performed with a smaller timestep Δt^num. This procedure repeats until the accuracy condition is fulfilled. Specifically, in our simulations an accuracy of the order of |r_L − r_measured| /r_L ≈ 10^-5 was used. This ensures that in each simulation the time discretisation is chosen correctly. Of course, the discretisation of spatial grid must resolve r_L as well. This means that twice the Larmor radius must be much larger than a grid cell and smaller than the whole system length.

Ultrarelativistic particle speeds present limitation to the method of the Borispush. In this case the conservation of energy would be violated, since the ideal ohmic law is not fulfilled anymore. Beyond Lorentz factors of γ ≈ 10³, fictitious forces start to act and this approach is no longer applicable (Vay 2008).

Both parts of Gismo are calculated for each step. After iterating the Elsässer MHD fields w^±, they are transformed into the physical, electric, and magnetic fields which are transferred to Gismo-Particles. Then the Borispush is performed. Each particle responds to its local fields, which are calculated by an averaging method via three-dimensional splines (Spanier & Wisniewski 2011; Wisniewski et al. 2012). Periodic spatial boundary conditions were used; thus the number of particles remained constant in each simulation.

2.2. Statistical description

The general statistical approach of particle transport by the relativistic Vlasov equation, $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{f}}_{\mathit{T}}}{\mathit{\partial t}}\mathrm{+}\frac{{p}}{\mathit{m\gamma }}\mathrm{·}\frac{\mathit{\partial }{\mathit{f}}_{\mathit{T}}}{\mathit{\partial }{x}}\mathrm{+}\frac{\mathit{q}}{\mathit{c}}\left(\mathit{c}{E}\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\mathrm{+}{v}\mathrm{\times }{B}\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\right)\mathrm{·}\frac{\mathit{\partial }{\mathit{f}}_{\mathit{T}}}{\mathit{\partial }{p}}\mathrm{=}{\mathit{S}}_{\mathit{T}}\mathrm{\left(}{x}\mathit{,}{p}\mathit{,t}\mathrm{\right)}\mathit{,}& & \end{array}$(18)describes the development of a particle distribution f_T of species T under the influence of the Lorentz force (Schlickeiser 2002). The momentum is denoted by p and the term on the RHS S_T is a source term. It is useful to transform this equation into guiding centre coordinates, since we are not interested in the actual position of the gyration, but in the averaged position of the particle. The Vlasov equation then reads $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{F}}_{\mathit{T}}}{\mathit{\partial t}}\mathrm{+}\mathit{v\mu }\frac{\mathit{\partial }{\mathit{F}}_{\mathit{T}}}{\mathit{\partial Z}}\mathrm{-}\mathrm{\Omega }\frac{\mathit{\partial }{\mathit{F}}_{\mathit{T}}}{\mathit{\partial \phi }}\mathrm{+}\frac{\mathrm{1}}{{\mathit{p}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{X}}_{\mathit{\sigma }}}\mathrm{(}\mathrm{⟨}{\mathit{p}}^{\mathrm{2}}{\mathit{g}}_{{\mathit{X}}_{\mathit{\sigma }}}\mathit{\delta }{\mathit{f}}_{\mathit{T}}\mathrm{⟩}\mathrm{)}\mathrm{=}{\mathit{S}}_{\mathit{T}}\mathrm{\left(}{\mathit{X}}_{\mathit{\sigma }}\mathit{,t}\mathrm{\right)}\mathit{,}& & \end{array}$(19)where F_T =  ⟨ f_T ⟩ is the expectation value of f_T and X_σ are the guiding centre coordinates with the angle in the perpendicular plane φ and the pitch-angle cosine μ. These coordinates are given by $\begin{array}{ccc}{\mathit{p}}_{\mathit{x}}& \mathrm{=}& \mathit{p}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\cos \mathit{\phi }\\ {\mathit{p}}_{\mathit{y}}& \mathrm{=}& \mathit{p}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\sin \mathit{\phi }\\ {\mathit{p}}_{\mathit{z}}& \mathrm{=}& \mathit{p\mu }\\ \mathit{x}& \mathrm{=}& \mathit{X}\mathrm{-}\frac{\mathit{c}}{\mathit{q}{\mathit{B}}_{\mathrm{0}}}\mathit{p}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\sin \mathit{\phi }\\ \mathit{z}& \mathrm{=}& \end{array}$(20)The generalised forces g_Xσ represent the effects of the electromagnetic fields and are basically the time derivatives of the denoted variables, e.g. g_Xσ = Ẋ_σ. This equation is in general analytically unsolvable.

A common approach to describing particle transport analytically is the QLT, which was first suggested by Jokipii (1966) in the context of energetic charged particle transport in turbulent magnetic fields. Its core is the assumption of unperturbed particle orbits. This implies that the fluctuation amplitudes are small, leading to a quasilinear system. The Vlasov equation Eq. (19) would then simplify to $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{F}}_{\mathit{T}}}{\mathit{\partial t}}\mathrm{+}\mathit{v\mu }\frac{\mathit{\partial }{\mathit{F}}_{\mathit{T}}}{\mathit{\partial Z}}& \mathrm{-}& \mathrm{\Omega }\frac{\mathit{\partial }{\mathit{F}}_{\mathit{T}}}{\mathit{\partial \phi }}\mathrm{=}{\mathit{S}}_{\mathit{T}}\mathrm{\left(}{\mathit{X}}_{\mathit{\sigma }}\mathit{,t}\mathrm{\right)}\\ & \mathrm{+}& \end{array}$(21)where the method of characteristics was applied and the hat symbol represents quantities along the characteristics (Schlickeiser 2002). This equation is known as the Fokker-Planck equation with the Fokker-Planck coefficients ${\mathit{D}}_{{\mathit{X}}_{\mathit{\sigma }}\widehat{{\mathit{X}}_{\mathit{\sigma }}}}$. One of the most interesting variables is the pitch-angle diffusion coefficient D_μμ. It describes the pitch-angle scattering of the particle and is consequently connected to the scattering mean free path, which can be evaluated by the observable angle distribution and Monte Carlo simulations (Agueda et al. 2009). Here, scattering means a resonant wave-particle interaction of nth order which fulfils the condition $\begin{array}{ccc}{\mathit{k}}_{\mathrm{\parallel }}\hspace{0.17em}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{-}\mathit{\omega }\mathrm{+}\mathit{n}\hspace{0.17em}\mathrm{\Omega }\mathrm{=}\mathrm{0}\mathit{,} \mathit{n}\mathrm{\in }& & \end{array}$(22)(Schlickeiser 1989), where ω is the wave frequency and k_∥ its parallel wavenumber. Ω is again the particle gyrofrequency and v_∥ its parallel velocity component. Whether a particle with v_∥ interacts resonantly with a wave with k_∥ is determined by the polarisation (Kennel & Engelmann 1966). Because our MHD model consists of pseudo and shear Alfvén waves only, certain values of n can be connected to the different types of waves. The Cherenkov resonance, n = 0, is generated only by waves with a compressive magnetic field (δB·B₀ ≠ 0), i.e. pseudo Alfvén waves. Purely parallel waves can contribute only to the n =  ± 1 resonance, and resonances with | n |  > 1 occur only for waves with nonvanishing perpendicular wavenumber, i.e. with oblique Alfvén waves. For a mathematical treatment of the resonant interactions, see Appendix B.

For further reading and a detailed insight of the derivation of Eq. (21) we refer to Schlickeiser (2002). A serious problem of the QLT is the limited applicability. The approximation of small fluctuations only holds for weak turbulent systems, where δB/B₀ ≪ 1. This is important for the local field which acts on the invidual particles. For instance, a strong turbulent region would change effectively the direction of the mean magnetic field and consequently the gyromotion of the particles. Hence the assumption of unperturbed orbits would be invalid. Another problem of the QLT is the inadequate description of particles propagating perpendicular to the mean magnetic field, μ ≈ 0. However, regarding Eq. (22) a very small parallel particle velocity will generate a singularity for v_∥ → 0. One aspect of this paper will concentrate on the applicability of the QLT to describe solar energetic particle transport.

2.3. The Fokker-Planck coefficient D_μμ

We will now focus on the pitch-angle scattering coefficient in more detail. The basic idea is to compare the simulation results of Gismo with the analytical approach.

Our simulations provide us with all the information of the trajectory and velocity of each individual particle within a turbulent plasma. Particle scattering will be considered here in terms of the pitch-angle cosine $\begin{array}{ccc}\mathit{\mu }\mathrm{=}\cos \mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\frac{{v}\mathrm{·}{B}}{\mathrm{\left|}{v}\mathrm{\right|}\mathrm{\left|}{B}\mathrm{\right|}}\mathit{,}& & \end{array}$(23)which describes the orientation of the particles’ velocity vector with respect to the local magnetic field. In the approximation of the QLT, there are different ways to calculate D_μμ by using the absolute change of μ during a certain time interval for a single particle. We used the definition $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu }}\mathrm{=}\underset{\mathrm{\Delta }\mathit{t}\mathrm{\to }\mathrm{\infty }}{\lim }\frac{\mathrm{(}\mathrm{\Delta }\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathrm{\Delta }\mathit{t}}\stackrel{\mathit{t}\mathrm{\gg }{\mathit{t}}_{\mathrm{0}}}{\mathrm{\approx }}\frac{\mathrm{(}\mathrm{\Delta }\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathrm{\Delta }\mathit{t}}\mathit{,}& & \end{array}$(24)where the time interval Δt = t − t₀ compared to an initial state at t₀ is assumed to be large, i.e. the time evolution t has to be sufficient to develop resonant interactions. The necessary time development is discussed in the result section. Instead of using the change of the pitch-angle cosine Δμ = μ_e − μ₀, the definition with the change of the pitch-angle Δα = α_e − α₀ itself leads to $\begin{array}{ccc}{\mathit{D}}_{\mathit{\alpha \alpha }}\mathrm{=}\frac{\mathrm{(}\mathrm{\Delta }\mathit{\alpha }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathrm{\Delta }\mathit{t}}\mathrm{=}\frac{{\mathit{D}}_{\mathit{\mu \mu }}}{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\mathit{,}& & \end{array}$(25)which represents the scattering frequency.

For the analytical approach we use the magnetostatic limit for the calculation of D_μμ for pseudo $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,P}}& \mathrm{\approx }& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{n}\ne \mathrm{0}}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\mathit{\pi \delta }\left({\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{+}\mathit{n}\mathrm{\Omega }\right)\\ & & \end{array}$(26)and shear Alfvén waves $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,A}}& \mathrm{=}& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\mathit{\pi \delta }\left({\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{+}\mathit{n}\mathrm{\Omega }\right)\\ & & \end{array}$(27)where P_xx is the xx-component of the magnetic field fluctuation tensor, which represents the turbulence power spectrum. These equations are based on Schlickeiser (2002) and a more detailed derivation is presented in Appendix B. We note that a proper description of the Cherenkov resonance in Eq. (26) is not possible in this formulation. Using n = 0 yields an undefined mathematical expression. These coefficients were semi-analytically solved with the set of equations given in Appendix C. We note that the approach of the magnetostatic QLT (SQLT) assumes the wave frequency to be small compared to the gyro frequency, ω ≪ Ω.

3. Simulation setup

We focus on the environmental conditions in the solar corona at a distance of three solar radii. The magnetic background field in this region is approximately B₀ = 0.174G, which yields, assuming a particle density of 10⁵cm^-3, an Alfvén speed of v_A = 1.2 × 10⁸cms^-1(Vainio et al. 2003). This region is of special interest because particle acceleration by CME-driven shocks is strongest there.

The outer length scale of the simulated system is L_scale = 3.4 × 10⁸cm. This value is estimated using the growth rate from Vainio (2003): $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathit{\pi \gamma }\mathrm{\Omega }}{\mathrm{2}{\mathit{n}}_{\mathrm{p}}{\mathit{v}}_{\mathrm{A}}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}\mathit{p}\mathit{v\mu }\mathrm{\left|}\mathit{k}\mathrm{\right|}\mathit{\delta }\left(\mathrm{\left|}\mathit{k}\mathrm{\right|}\mathrm{-}\frac{\mathrm{\Omega }}{\mathit{v}}\right)\mathit{f,}& & \end{array}$(28)with proton cyclotron frequency Ω, proton speed v, the Lorentz factor γ, the proton number density n_p, as well as μ, the pitch-angle cosine, and f the proton distribution function. As the resonance condition, Eq. (22) is applied, assuming that the particles are streaming along the field.

A remark on the notation: the wavenumber is in general defined as k = (2πj)/L, where j stands for the numerical grid position. For simplicity we used the normalised wavenumbers k′ = kL = (2π·j) throughout.

We generated peaks at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{8}$ and ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$, which can be excited by streaming protons with energies of E_p ≈ 64MeV and E_p ≈ 7MeV respectively. Using the resonance condition, this leads to a length scale of $\begin{array}{ccc}{\mathit{L}}_{\mathrm{scale}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi j}}{\mathit{e}{\mathit{B}}_{\mathrm{0}}}\mathit{\gamma }{\mathit{m}}_{\mathrm{p}}\mathit{cv}\mathrm{\approx }{\mathrm{10}}^{\mathrm{8}}\mathrm{cm}\mathit{.}& & \end{array}$(29)On this scale the magnetic background field is assumed to be homogeneous. The spatial resolution is 256³ gridpoints, resulting in 128³ points in k-space of which | k^′ |  = 2π [0···42] wave modes are active modes that remain unaffected by (anti)aliasing. The resistivity parameter was set to ν = 1 in numerical units and the hyperdiffusivity coefficient to h = 2.

The turbulence was simulated using an anisotropic driving mechanism. Energy is, therefore, constantly injected into the simulation volume at the first 5 numerical wavenumbers in perpendicular (${\mathit{k}}_{\mathrm{\perp }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{\left[}\mathrm{0}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{4}\mathrm{\right]}$) and 15 in parallel direction (${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{\left[}\mathrm{0}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{14}\mathrm{\right]}$). The reference frame is the solar wind frame.

The anisotropy was chosen for two reasons: first, to mimic the preferential direction of the solar wind, where particles streaming radially away from the Sun form the Parker spiral. Consequently, these particles can deposit their energy in a parallel direction on different scales. This is mainly valid in the vicinity of the Sun, which we are interested in. Second, a slab component of solar wind turbulence is also observed at small scales in the parallel direction. To ensure turbulence evolution up to high parallel wavenumbers, the driving range was extended along the parallel axis. This is necessary because the parallel evolution is much weaker than the perpendicular one (Goldreich & Sridhar 1995, 1997). Even though this is primarily a technical aspect to ensure the extent of the spectrum to higher k_∥, it is still in line with observations. An isotropic driver would not yield sufficiently turbulent modes at high k_∥.

The turbulence driving is performed by allocating an amplitude with a phase to the Elsässer fields within the Fourier space. The amplitude follows a power-law of | k | ^-2.5 and is initialised using a normal distribution. The phase was randomly chosen between zero and 2π. These settings are divergence free and Hermitian symmetric. After this initialisation the values were scaled to the desired scenario, which in our case is a δB/B₀ ratio of roughly 10^-2. We note that both species, pseudo and shear Alfvén waves, are excited by this type of turbulence driving, but as presented by Maron & Goldreich (2001) the pseudo wave evolution is strongly suppressed. In this inertial range energy is injected every 0.03s (10 MHD steps), which leads to a saturated turbulence, an equilibrium between dissipation and injection.

A Gaussian-distributed energy peak with purely parallel wavenumber k = ke_∥ was injected within the saturated turbulence. For this purpose, two different wave modes were chosen. To investigate the physics of an SEP-event, a wavenumber of k_∥ = 1.5 × 10^-7cm^-1 is used. This corresponds to a numerical wavenumber of 2π·8, which is still within the driving range of the turbulence. To represent injection at smaller scales, a peak is set at k_∥ = 4.4 × 10^-7cm^-1. With ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$, this wave mode is not in the driving range of the background turbulence. The injection of energy in these modes was done gradually over a specific time interval.

One important aspect of the calculations of D_μμ is the applicability of the QLT. The ratio δB/B₀ is, therefore, the limiting parameter (see Sect. 2.2). Consequently, it is of interest to explore peaks at either position under the influence of different magnetic background fields. For this purpose another turbulence simulation was performed which uses the same initial conditions, except for a higher mean field ${\mathit{B}}_{\mathrm{0}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathit{h}}\mathrm{=}\mathrm{1.74}$ G and for stability reasons a higher dissipation rate ν = 10cm^2hs^-1. In the following the smaller field simulation will be denoted by ${\mathit{B}}_{\mathrm{0}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathit{l}}\mathrm{=}\mathrm{0.174}\mathrm{G}$. We note that the ${\mathit{B}}_{\mathrm{0}}^{\mathit{h}}$ case is an artificial scenario rather than representative of coronal conditions. Peaks in such high magnetic fields represent SEP energies up to 2.65 GeV. In this case the SEPs cannot generate waves by streaming as the relativistic proton intensities are insufficient.

To summarise, these four simulations with excited wave modes at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{8}$ and ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$, both within turbulent fields governed by a strong and a weak B₀, are the starting point for the test particle simulations and the calculations of the scattering coefficient D_μμ.

All of the test particle simulations by Gismo-Particles were performed with 10⁵ protons with an initial uniform distribution in μ and φ and random positions x. This amount has proven to be sufficient in test simulations for good statistics. These initial conditions aim to provide a complete coverage of the phase space in μ for the test particles. Thus, we are not interested in the development of special distribution functions.

A constant absolute value of the momentum was chosen so that particles with a certain parallel velocity component fulfil the resonance condition Eq. (22). Consequently, a resonant value of μ$\begin{array}{ccc}{\mathit{\mu }}_{\mathrm{R}}\mathrm{=}\frac{\mathit{\omega }\mathrm{-}\mathit{n}\hspace{0.17em}\mathrm{\Omega }}{{\mathit{k}}_{\mathrm{\parallel }}\hspace{0.17em}\mathit{v}}\mathrm{=}\frac{\mathit{\omega }\mathrm{-}\mathit{n}\hspace{0.17em}\mathrm{\Omega }}{{\mathit{L}}_{\mathrm{scale}}^{-1}\hspace{0.17em}{\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\hspace{0.17em}\mathit{v}}& & \end{array}$(30)must be within the interval [ − 1;1] for a given particle speed v. Since Eq. (30) is dominated by Ω and hence B₀ for particles propagating significantly faster than the Alfvén speed, the particle speed v has to be different for ${\mathit{B}}_{\mathrm{0}}^{\mathit{l}}$ and ${\mathit{B}}_{\mathrm{0}}^{\mathit{h}}$. The parameter setup and the resonant μ_R for each simulation is summarised in Table 1.

To investigate the particle scattering in different stages of the turbulence evolution, multiple simulations were performed. The most promising match between QLT and the simulation results is expected to be in stages with small δB/B₀. According to this, the particle simulations were performed in the driving phase of the peaks as well as in the decaying phase. A simulation at maximum driven peaks would simply lead to random scattering where no reasonable prediction can be made by QLT.

The simulations with the low magnetic background field ${\mathit{B}}_{\mathrm{0}}^{\mathit{l}}$ were also performed with a higher spatial resolution of 512³ grid cells, e.g. 256³ wave modes. For reasons of clarity the results of these simulations are presented and discussed in the Appendix D.

4. Results and discussion

4.1. A toy model for wave-particle resonance

In order to motivate our particle-scattering simulations, we want to present a simple model. Because of the limited applicability of QLT, it is necessary to understand its range of validity. In detail, the δB/B₀ ratio and the time development are of utmost importance. Different simplified scenarios are therefore presented below.

We consider a circularly polarised Alfvén wave with k = (0,0,k_∥)^T. Its magnetic field is described by $\begin{array}{ccc}{B}\mathrm{=}{{B}}_{\mathrm{0}}\mathrm{+}\mathit{\delta B}\mathrm{(}{{e}}_{\mathit{x}}\cos \mathrm{(}\mathrm{\pm }\mathit{kz}\mathrm{+}\mathit{\psi }\mathrm{)}\mathrm{+}{{e}}_{\mathit{y}}\sin \mathrm{(}\mathrm{\pm }\mathit{kz}\mathrm{+}\mathit{\psi }{\mathrm{)}}^{\mathrm{)}}\mathit{,}& & \end{array}$(31)in the wave frame where B₀ ≫ δB is assumed. The interaction between particles and the wave will then change the particles’ pitch-angles via scatterings. A detailed description of this process as well as a derivation for the time dependency of μ is given in Appendix A. Under the assumption givenby Eq. (31), the expression of the time evolution of μ in Eq. (A.8) will simplify to $\begin{array}{ccc}\mathrm{\Delta }{\mathit{\mu }}^{\mathrm{\pm }}\mathrm{\left(}\mathit{t,\psi }\mathrm{\right)}& \mathrm{=}& \mathrm{\Omega }\sqrt{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}\hspace{0.17em}\frac{\mathit{\delta B}}{{\mathit{B}}_{\mathrm{0}}}\\ & & \end{array}$(32)when the calculation is performed along the unperturbed orbit, μ = constant. The change in μ is then maximal at the phase $\begin{array}{ccc}{\mathit{\psi }}_{\mathit{M}}^{\mathrm{\pm }}\mathrm{=}\arctan \left[\frac{\sin \mathrm{(}\mathrm{\pm }\mathit{kv\mu }\mathrm{-}\mathrm{\Omega }\mathrm{)}\mathrm{\Delta }\mathit{t}}{\mathrm{1}\mathrm{-}\cos \mathrm{(}\mathrm{\pm }\mathit{kv\mu }\mathrm{-}\mathrm{\Omega }\mathrm{)}\mathrm{\Delta }\mathit{t}}\right]& & \end{array}$(33)and its identical solutions at periods ${\mathit{\psi }}_{\mathit{M}}^{\mathrm{\pm }}\mathrm{+}\mathit{n\pi }$.

To study this model, a single linearly polarised Alfvén wave was simulated, propagating undisturbed with v_A = 10⁷cm   s^-1 towards positive z-direction (w⁺) with a purely parallel wave vector k′ = 2π·(0,0,1)^T. It should be noted that the magnetic field given by Eq. (31) does not describe a propagating wave but a static disturbance. Linear polarisation was chosen to cover both resonances at ± μ_R (respectively both signs in Eq. (32)), because interactions with μ < 0 are caused by left-handed circular polarised waves and μ > 0 by right-handed ones. The wave amplitude was set to δB/B₀ = 0.1, where B₀ = 4.34 × 10^-4 G. The outer length scale of the simulation cube was set to 10⁸ cm.

In Gismo-Particles, 5 × 10⁵ protons with random μ, φ and position x were initialised. The absolute value of v was set to 2 × 10⁸cm   s^-1, where particles for n =  −1 and n =  + 1 are resonant at μ_R1 =  −0.28 and μ_R2 =  + 0.38 according to Eq. (22).

In Fig. 1, we show the pitch-angle change of the simulated particles and the QLT-predicted change as a function of the initial pitch-angle for wave amplitude δB/B = 0.1 after two gyration periods. The simulated particles are represented as a scatter plot, whereas the QLT result for different wave modes are represented by curves. As expected, two maxima develop at the predicted positions μ_R1 and μ_R2. The adjoining maxima are caused by particles interacting with the wave nonresonantly, henceforth named nonresonant ballistic interactions. These decrease their amplitudes with increasing time. Also, their positions in μ do not remain constant, like the resonant interactions, because the ballistic maxima become more numerous and move closer to the true resonance. The tilt of the maxima deviates from the prediction of Eq. (32), as clearly shown by comparison with the theoretical curves. Furthermore, the resonances of Eq. (32) are located symmetrically at μ_r =  ± 0.33. This shift can be explained by the choice of the inertial system. The derived equation is within the wave frame. The transformation into the laboratory frame (under the assumption of Galilean transformation) via ${\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\approx }\mathit{\mu }\mathrm{-}\frac{{\mathit{v}}_{\mathit{A}}}{\mathit{v}}$(34)will fit to the observed position in μ. This shift is caused by the magnetostatic approximation, since the difference between Eq. (30) and its magnetostatic version (ω = 0) is exactly the term v_A/v. If this ratio is small, as in most of our simulations, the effect is negligible. However, for the simulations with stronger ${\mathit{B}}_{\mathrm{0}}^{\mathit{h}}$ this poses a problem, which will be discussed in the results section.

The tilt cannot be explained by Eq. (32) due to its dependence on the QLT. It stems from the distortion of the orbit by the wave’s δB and consequently, arises by the presentation of Δμ as a function of the initial pitch-angle, μ₀. Because of the finite δB/B₀-ratio, the assumption of unperturbed orbits, which is the basic approximation in QLT, does not hold anymore. The change of μ should therefore be calculated as an average over the whole trajectory of each particle, which would introduce a nonlinear correction to the theory. A simple solution to correct the tilt is the presentation of Δμ dependent not on the intial μ₀, but on the mean value ⟨ μ ⟩ of the initial and the final value of the pitch-angle cosine. We note that the angle of the tilt does not depend on δB/B₀, but the scattering amplitude does. Consequently, the tilt is more visible and influences further calculations more as the δB/B₀-ratio increases. This visibility is called effective tilt hereafter.

In order to prove the explanations given above, we performed a similar simulation with smaller wave amplitude and longer time development. As shown in Fig. 2 the QLT is valid within this scenario. As the scattering amplitude is significantly smaller, the effect of the tilt resulting from using the initial pitch-angle is negligible and the time evolution leads to a clear resonance as the side maxima become smaller. But even under these convenient conditions the deviation from the QLT is still visible.

The last case of our toy model uses an oblique wave with k′ = 2π·(0,1,1)^T, which was simulated under the same conditions as before. The result in Fig. 3 clearly shows the higher orders of resonance with n =  [−3, −2, −1,1,2]. A prediction of Eq. (32) is not applied, since it would not describe the results of the simulation very well for two reasons. First, the assumption that the prediction by QLT becomes better with t → ∞ would not hold because of the finite δB/B₀ and hence perturbed orbits. Second, the nonvanishing Bessel functions have to be taken into account and would modify Eq. (32). But although the ballistic maxima become significant in number and generate a nonresonant background scattering, the resonances are clearly visible. Another interesting fact is that the scattering caused by resonance at |n|  = 2 exceeds the ones at | n |  = 1. This means that in the case of an oblique wave, higher order resonances can become even more important than the fundamental ones. The time development for the purely parallel wave mode at 50 gyration periods showed expectedly only the |n|  = 1 resonances (not shown in Fig. 3 for reasons of clarity).

4.2. Results of particle scattering in turbulence with amplified modes

We first present selected results from Lange & Spanier (2012), which are the base of our particle simulations. However, we do not discuss the evolution of the peaks in detail, but focus on the stages of the driving and decay of the peaks on which our particle simulations were performed.

In Figs. 4 and 5 we show the situation for the decaying peaks. The timesteps of these figures are the starting points of our particle simulations for the decay stage. We note that the timesteps were chosen to retrieve roughly the same order of magnitude in δB/B₀ at the peak positions, which are 10^-3. All of the peaks show a strong perpendicular evolution at this stage. Furthermore, higher harmonics of ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{8}$ at positions 16, 24, 32, 40 are present. These harmonics, however, do not interact dominantly with the particles as they have already lost most of their energy compared to the maximum driven stage (not shown here,) as indicated by the logarithmic colour scaling.

As shown in the previous section on the toy model, the resonance takes several gyration periods to develop. On the other hand, if the particles are simulated for too long, the interactions with the turbulence will smooth the resonances and lead to unstructured scattering. This is again caused by the finite perturbation of the particle orbits, which are assumed to be negligible for the applicability of QLT. The resonances would not converge to a δ-function as predicted by QLT, but be scattered ballistically at the finite distortions of the fields. We find that a reasonable timescale is between 10 and 30 gyration periods for parameters used in this paper. For weaker turbulence, even results at 50 gyrations might yield a structure, but in most cases the resonances develop clearly in shorter timescales. The pitch-angle scattering coefficient D_αα was calculated for different stages of the evolution. In Fig. 6 we show D_αα for the first simulation setup with ${\mathit{B}}_{\mathrm{0}}^{\mathit{l}}\mathrm{=}\mathrm{0.174}$ G within the 256³ grid, with the peak position ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ at the driven stage. We observe a convergence between 5 and 10 gyration periods. As presented in Table 1, the resonant interaction with the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ mode is located at μ_R = 0.29, whereas Fig. 6 shows two maxima at μ = 0.2 and μ = 0.4, which seem to move away from each other during the time development.

The real resonant structure is revealed by the scatter plot in Fig. 7. Indeed, the resonance is centred at the predicted position, but tilted for the same reasons presented with the toy model. The tilt spreads the particles to a wide Δμ range, which results in splitting of the maximum of the pitch-angle diffusion coefficient into two maxima at both sides of the resonant pitch-angle μ_R = 0.29. This is because the calculation of D_αα in the QLT with Eq. (25) is not dependent on the sign of Δμ. Consequently, the scattering coefficient is mapped due to the square value of Δμ to two different maxima. This also explains the movement of the maxima between five and ten gyration periods in Fig. 6, because Δμ increases with time. The smaller resonance at n =  −1, i.e. μ_R =  −0.27, is caused by the different polarisations of the peaked mode. That means resonances with μ < 0 are caused by left-handed circular polarised parts of the peaked mode and μ > 0 by right-handed ones. Furthermore, in the scatter plot the Cherenkov resonance n = 0 is visible. It is hardly observable in the D_αα plot due to the dominant | n |  = 1 resonances and their tilt.

The test particle simulation with the decay phase of the peaked mode at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ shows resonant interactions beyond the fundamental resonance. For example, the maximum located near μ = 0.6 in Fig. 8 represents the n = 2 interaction. We note that the tilt in the corresponding scatter plot in Fig. 9 causes again the split of the maximum in D_αα.

As already shown by the toy model, these higher orders are generated by oblique Alfvén waves. Those modes clearly exist in the decay stage of the peaked mode. For example, in Fig. 4 the maximum of the turbulence energy spectrum is shifted from the k_∥ axis towards higher perpendicular wave modes. However, within the driven stage of the peak the energy is injected in purely parallel modes, which causes the dominant | n |  = 1 resonances. Furthermore, the resonance pattern in Fig. 9 indicates that the left-handed circular polarised parts of the peaked mode decayed faster, which lowers the scattering frequency for particles with μ < 0. While this observation needs further investigation, it seems to be connected to the turbulence evolution, since the wave-particle toymodel does not show this behaviour. Additionally, the magnetic background field plays a role, which is discussed for the results of the other parameter setup below. Again the Cherenkov resonance n = 0 is visible.

The results for the particle scattering at the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{8}$ driven peak are presented in Fig. 10. During the driven stage of the peak, the test particles interacted resonantly at | n |  = 1 and the predicted pitch-angle μ_R = 0.86. The effect of the tilt on the resonances is smaller compared to the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ peak. Consequently, the maxima in the scattering coefficient are not split.

The situation is different for the decay stage. The energy has spread significantly to oblique modes (see Fig. 4, left-hand frame). This leads to stronger scattering and hence a stronger tilt, which is clearly visible in Fig. 11. Even the sharp lines of the maximum change of μ are observed, because several particles were strongly scattered. In this state, D_αα has no clear structure.

The last set of simulations used an increased magnetic background field B₀. Thus, the δB/B₀ ratio is decreased by an order of magnitude, which is more consistent with the assumptions of the QLT. This claim is supported by Figs. 12 and 13 where each resonance has a reduced effective tilt because of the smaller scattering frequency and represents a dominant structure. We note that the evolution time is roughly twice as long as within the weaker magnetic background field. A strong Cherenkov resonance is observed for the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{8}$ peak. Also both n =  | 1 | resonances are visible. The small increase of the scattering rate at 0.53 could be generated by a resonant interaction with the higher harmonic at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{16}$, which is also a dominant wave mode in this stage (see Fig. 5, left-hand frame).

For the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ peak, the fundamental resonance n = 1 is strongest. A higher order of resonance is just slightly visible at μ_R = 0.70. An indication of the decay of the left-handed modes is the absence of resonant interactions n < 0. These are visible during the driving of the peak, but vanish in the decay stage.

4.3. Comparison between SQLT and particle simulations

In the last section we present results of the SQLT approach and compare them to the particle simulations. For this purpose, the pitch-angle scattering coefficient was calculated by using Eqs. (26) and (27) for each spectrum. In contrast to the particle simulations, the calculations via the SQLT approach produce a clear structure in the scattering coefficient, see Fig. 14. That means that the resonances are not broadened by finite time and the effect of the tilt (see the toy model section) does not influence the values of D_αα.

As presented in Fig. 14, the magnetostatic calculations of the SQLT reproduce the resonant structures for n = 1,2 very well. However, the resonance gap for small pitch-angle cosine values causes a strong drop-off below | μ |  ≲ 0.25. Furthermore, a slab approximation was made by the projection of the total wave energy to the parallel modes, i.e. each wave mode is assumed to be parallel. The slab comparison gives a good approximation of the scattering coefficient, except for the resonant structures caused by the oblique waves (| n |  > 1). The scattering caused by the pseudo waves is comparable to the shear waves at the positions of the peaks. For values of μ =  [0.35;0.55], the pseudo mode has higher scatter rates than the shear mode, whose coefficient D_αα is bigger elsewhere. In the following, we restrict the presentation to the D_αα sum of both wave modes, e.g. the fourth curve in Fig.14.

As expected, the comparison between SQLT and particle simulation is best for the pure turbulent background case, as long as the Cherenkov resonance is small or for higher values of μ. Because of the resonance gap, it cannot be reproduced by the SQLT approach. However, the particle simulation presents a clear and broad interaction for n = 0, as shown in Fig. 15. Especially for the background and driven stage, the curves increase towards μ = 0, which is caused by the Cherenkov resonance. For a clearer interpretation of D_αα it is again helpful to compare the D_αα curve with the corresponding scatter plot. According to the previous section, the D_αα is not only increased at μ = 0 but also broadened by the tilt of the resonant peak. Nevertheless, the SQLT curve becomes comparable to the particle simulation for μ > 0.6. Because the scattering is concentrated at the purely parallel peak during the driven stage, the scattering coefficient of the particle simulation shows the resonance very well and is also comparable to the SQLT calculation. Even more complicated is the situation for the scattering at the peaked modes in the decay stage. The SQLT curve in Fig. 15 shows the resonance at the predicted position μ_R = 0.86, whereas the D_αα of the particle simulation seems to have no structure. This is again caused by the visible effect of the tilt at the resonances (see Fig. 11).

Because the effective tilt within the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ peak simulations is stronger and the resonance pattern for the decay stage is more complex, the particle simulation results for D_αα differ from the SQLT calculations. Nevertheless, the resonances in Fig. 16 are located at the predicted positions.

The results of the simulation with a higher magnetic background field are presented in Figs. 17 and 18. The qualitative comparison to the SQLT is much better with the smaller δB/B₀ ratio. The particle scattering coefficients show clear resonant maxima at the predicted positions. A small shift of the SQLT resonances stems from the magnetostatic assumption, which sets ω in Eq. (30) to zero. This is valid for ω ≪ Ω, which is not fulfilled for this simulation, since ω₈ = 177s^-1, ω₂₄ = 532s^-1 and Ω = 4264s^-1. Thus the resonances are slightly shifted towards smaller μ_R, e.g. the μ_R = 0.37 for the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ peak is μ_R = 0.33 in the magnetostatic limit. Another problem is the resonance gap, which is significantly broader. Consequently, the SQLT curves cannot be compared quantitatively to the particle simulations.

5. Conclusions

To summarise, we firstly presented a simple wave-particle interaction model to achieve a fundamental understanding of the underlying processes and interpretation of numerical simulation results. The effects of tilted and shifted resonances were observed and further compared to QLT predictions. Additionally, the difference between ballistic and resonant interactions could be shown in the context of time development.

Afterwards, we presented three different scenarios:

• 1)

  The first scenario used parameters to resemble conditionswithin a heliospheric plasma at three solar radii.

• 2)

  To study the influence of the resolution, the first simulation setup was used with a finer grid of 512³ cells (Appendix D).

• 3)

  The connection between QLT to a smaller δB/B₀ ratio was investigated by an approach using an artificially high B₀.

Each scenario was simulated with two different amplified wave modes with peak positions at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{8}$ and 24. Furthermore, two different evolution stages of the peak within the turbulence were investigated, which allows interactions with purely parallel peaked modes to be studied as well as oblique modes because of the peak development towards a perpendicular direction. All of the simulations were performed with 10⁵ protons. Their initial speed was chosen to be in resonance with the peaks. In addition, every turbulence spectrum was used for a magnetostatic QLT approach to calculate the pitch-angle scattering coefficient. We presented resonance patterns with different methods: the scatter plot Δμ   (μ₀) and the scattering coefficient D_αα   (μ₀) for both the numerical and semi-analytical quasilinear approach.

Our results show a good agreement between hybrid MHD particle simulations and QLT calculations for the scattering coefficient D_αα for a turbulent broad-band spectrum. The comparison reveals that the SQLT calculations based on the MHD spectrum are missing the Cherenkov resonance (n = 0), which is not covered by the QLT calculations used because quasilinear theory produces a singularity at μ = 0. Furthermore, the resonance gap at small values of | μ | is caused by the power spectrum, which does not have sufficient energy at high wave numbers within the dissipation range of the turbulence. The direct particle simulations, on the other hand, show problems linked to broadened resonances: Because of the finite simulation time, the wave-particle resonance peaks are broadened, which in turn yields a significant scattering background. With increasing simulation time, especially in simulations with higher δB/B₀, the resonances did not become narrow but scattering would start to randomise. Consequently, the assumption of δ-shaped resonances does not apply to real particle scattering. Furthermore, the tilt in Δμ, as shown in the scatter plots (e.g. Fig. 7), causes spreading of the resonant peaks over an interval of | μ | for the particle simulations.

Due to the tilt and the broadening of the resonances, the application of the scattering coefficient resulting from Eq. (25) is questionable. The D_αα results are either difficult to interpret or without any reasonable structure. Consequently, the QLT approach does not compare well to the particle simulation results. We point out that in these cases a very nice tool in numerical simulations is the scatter plot. By evaluating the total change of μ versus its initial state μ₀, even small resonance patterns become visible. This can be used for correct interpretation of the scattering coefficients. Especially in cases with multiple overlapping resonances (e.g. Fig. 9), the scatter plot yields important information. It should be noted, however, that the analysis of particle-scattering data requires more in-depth research.

A nice comparison between QLT and particle simulations is achieved with smaller δB/B₀ ratio of the turbulence and the peaked modes. This is not unexpected, because the assumption of unperturbed orbits is more reasonable. In this case the effective tilt in the scatter plots decreases significantly and D_αα becomes more structured. Unfortunately, the resonance gap in the QLT approach expands in these scenarios, which leads to smaller absolute values of the scattering coefficient. This large gap stems from the used spectrum, where the magnetic background field was increased, leading to less resonant waves.

We conclude that for realistic particle transport, especially around μ ≈ 0, the QLT approach does not yield physical scattering coefficients. This is not unexpected, as it has been discussed in literature already (Shalchi et al. 2004). We are not aware of any self-consistent theory describing the numerical particle results derived from our hybrid simulations. Further investigations are needed to disentangle the different transport processes involved in order to develop a new kind of transport theory. We also note that even our numerical approach has limited validity for high δB/B₀ ratios. We will discuss more sophisticated numerical methods for transport analysis in a forthcoming paper.

Online material

Appendix A: Time dependency of the pitch-angle cosine

The interaction of a charged particle with the fluctuation δB of an Alfvén wave forces the particle to leave its gyro orbit. This procedure is referred to pitch-angle scattering. Consequently the change of parallel momentum is connected to this process via $\mathit{p}\frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{d}{\mathit{p}}_{\mathrm{\parallel }}}{\mathrm{d}\mathit{t}}$. Then the time derivative of p_∥ is given by the Lorentz force $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{p}}_{\mathrm{\parallel }}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{ec}}{\mathrm{2}{\mathit{E}}_{\mathrm{p}}}\sum _{\mathrm{\pm }}\mathrm{(}\mathrm{\mp }\mathrm{i}\mathrm{)}{\mathit{p}}_{\mathrm{\perp }}\mathit{\delta }{\mathit{B}}_{\mathit{l,r}}\mathit{,}& & \end{array}$(A.1)where E_p is the energy of the particle, which is assumed to be constant as the scattering is purely elastic. Since the parallel momentum remains constant during an unperturbed gyration, only the perpendicular momentum changes due to the circular movement p_⊥ = p_x ± ip_y = p_⊥ 0exp [ ∓ i(φ₀ − Ωt)]. This leads to a total change of μ for a particle with intial phase φ₀ to the magnetic field with $\begin{array}{ccc}\frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{ec}}{\mathrm{2}\mathit{p}{\mathit{E}}_{\mathrm{p}}}\sum _{\mathrm{\pm }}\mathrm{(}\mathrm{\mp }\mathrm{i}\mathrm{)}{\mathit{p}}_{\mathrm{\perp }\mathrm{0}}\exp \mathrm{[}\mathrm{\mp }\mathrm{i}\mathrm{(}{\mathit{\phi }}_{\mathrm{0}}\mathrm{-}\mathrm{\Omega }\mathit{t}\mathrm{\left)}\mathrm{\right]}\mathit{\delta }{\mathit{B}}_{\mathit{l,r}}\mathrm{\left(}{r}\mathrm{\right)}\mathit{,}& & \end{array}$(A.2)(Lee & Lerche 1974). A similar approach is given by Schlickeiser (2002) using Eq. (19). In this case the solution along the characteristics of the generalised force term is given by Eq. (12.2.4b) Schlickeiser (2002). This result reduces to Eq. (A.2) by using the magnetostatic approximation, which assumes the electric field fluctuations to be negligible, δE = 0. The coordinates used are still Eqs. 20. The fluctuation δB, i.e. of an Alfvén wave as used in the presented case, is given in Fourier space by $\begin{array}{ccc}\mathit{\delta }{\mathit{B}}_{\mathit{l,r}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{8}{\mathit{\pi }}^{\mathrm{3}}}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}{k}\hspace{0.17em}\mathit{\delta }{\mathit{B}}_{\mathit{l,r}}\mathrm{\left(}{k}\mathrm{\right)}\exp \mathrm{\left(}\mathrm{i}{k}{r}\mathrm{\right)}\mathit{.}& & \end{array}$(A.3)The exponential function can be described by the generating Bessel functions J_n$\begin{array}{ccc}\exp \mathrm{\left(}\mathrm{i}{k}{r}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathit{J}}_{\mathit{n}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\exp \mathrm{i}\mathrm{(}{\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathit{t}\mathrm{+}\mathit{n}\mathrm{(}\mathrm{\Psi }\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}\mathrm{+}\mathrm{\Omega }\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{,}& & \end{array}$(A.4)where the argument of J_n is $\begin{array}{ccc}\mathit{z}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{\perp }}\mathit{v}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}}{\mathrm{\Omega }}& & \end{array}$(A.5)and Ψ = cot^-1(k_x/k_y), ${\mathit{k}}_{\mathrm{\perp }}\mathrm{=}\sqrt{{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}}$. The total time derivative of the pitch-angle cosine, separated in parallel and perpendicular interactions, reads then $\begin{array}{ccc}\mathit{g̅}\mathit{\mu }\mathrm{=}\frac{\mathrm{d}\mathit{\mu }}{\mathrm{d}\mathit{t}}& \mathrm{=}& \frac{\mathit{ec}}{\mathrm{2}\mathit{p}{\mathit{E}}_{\mathrm{p}}}\sum _{\mathrm{\pm }}\mathrm{(}\mathrm{\mp }\mathrm{i}\mathrm{)}{\mathit{p}}_{\mathrm{\perp }\mathrm{0}}\exp \mathrm{[}\mathrm{\mp }\mathrm{i}{\mathit{\phi }}_{\mathrm{0}}\mathrm{]}\\ & & \mathrm{\times }\frac{\mathrm{1}}{\mathrm{8}{\mathit{\pi }}^{\mathrm{3}}}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}{k}\mathit{\delta }{\mathit{B}}_{\mathit{l,r}}\mathrm{\left(}{k}\mathrm{\right)}\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathit{J}}_{\mathit{n}}\mathrm{(}{\mathit{z}}^{\mathrm{)}}\\ & & \mathrm{\times }\exp \mathrm{\left[}\mathrm{i}\mathit{n}\mathrm{\right(}\mathrm{\Psi }\mathrm{-}{\mathit{\phi }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\mathrm{i}\mathit{t}\mathrm{(}{\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{\pm }\mathrm{\Omega }\mathrm{-}\mathit{n}\mathrm{\Omega }\mathrm{\left)}\mathrm{\right]}\mathit{,}\end{array}$(A.6)For Alfvén waves δB is aligned towards k × e_zB₀ and consequently δB_l,r(k) = δB(k)( ∓ i)exp( ± iΨ). The time integration and the identity $\begin{array}{ccc}\frac{\mathrm{2}\mathit{n}}{\mathit{z}}{\mathit{J}}_{\mathit{n}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{J}}_{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{+}{\mathit{J}}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& & \end{array}$(A.7)(Abramowitz & Stegun 1965) then gives the time dependency of the pitch-angle cosine with $\begin{array}{ccc}\mathrm{\Delta }\mathit{\mu }\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \mathit{\mu }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{\mathit{\mu }}_{\mathrm{0}}\mathrm{=}\mathrm{-}\frac{\mathit{e}\hspace{0.17em}\mathit{c}\hspace{0.17em}{\mathit{p}}_{\mathrm{\perp }\mathrm{0}}}{\mathrm{2}\mathit{p}{\mathit{E}}_{\mathrm{p}}\mathrm{8}{\mathit{\pi }}^{\mathrm{3}}}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}{k}\mathit{\delta B}\mathrm{\left(}{k}\mathrm{\right)}\\ & & \mathrm{\times }\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{\left[}{\mathit{J}}_{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right(}\mathit{z}\mathrm{)}\mathrm{+}{\mathit{J}}_{\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{(}\mathit{z}\mathrm{\left)}\mathrm{\right]}\\ & & \end{array}$(A.8)The application of a purely parallel propagating wave, as shown in our toy model, will then simplify this equation because z(k_⊥ = 0) = 0. Thus, all Bessel functions vanish, except J₀(0) = 1. This is the case for n =  ± 1. Furthermore, the k integration reduces by the assumption of a single wave , which leads to the presented form in Sect. 4.1.

Appendix B: Derivation of the pitch-angle diffusion coefficients

In the simulations, the turbulence consists of Alfvén and pseudo Alfvén waves, thus the pitch-angle diffusion coefficient must be calculated separately for the two modes. The two modes are decomposed using the method presented Maron & Goldreich (2001). For the pitch-angle diffusion coefficient for Alfvén waves, Schlickeiser (2002) gives $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,A}}& \mathrm{=}& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\mathrm{ℛ}\left({k}\mathit{,\omega }\right){\left[\mathrm{1}\mathrm{-}\frac{\mathit{\mu \omega }}{{\mathit{k}}_{\mathrm{\parallel }}\mathit{v}}\right]}^{\mathrm{2}}\\ & & \mathrm{\times }{\left[\frac{\mathit{n}\hspace{0.17em}{\mathit{J}}_{\mathit{n}}\left({\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }\right)}{{\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }}\right]}^{\mathrm{2}}{\mathit{P}}_{\mathit{xx,}\mathrm{A}}\mathrm{\left(}{k}\mathrm{\right)}\mathit{,}\end{array}$(B.1)where μ, v and Ω again are the particle’s pitch-angle cosine, speed, and gyrofrequency, ${\mathit{v}}_{\mathrm{\perp }}\mathrm{=}\mathit{v}\hspace{0.17em}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}$, P_xx,A the xx-component of the Alfvén mode turbulence power spectrum tensor, and B₀ the background magnetic field. For the resonance function ℛ(k,ω) we assume no damping, which gives the delta function ℛ(k,ω) = π   δ(k_∥v_∥ − ω + nΩ). Thus, in the magnetostatic limit (ω = 0) we have $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,}\mathrm{A}}& \mathrm{=}& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\mathit{\pi \delta }\left({\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{+}\mathit{n}\mathrm{\Omega }\right)\\ & & \mathrm{\times }{\left[\frac{\mathit{n}\hspace{0.17em}{\mathit{J}}_{\mathit{n}}\left({\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }\right)}{{\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }}\right]}^{\mathrm{2}}{\mathit{P}}_{\mathit{xx,}\mathrm{A}}\mathrm{\left(}{k}\mathrm{\right)}\mathit{.}\end{array}$(B.2)For the magnetosonic waves, Schlickeiser (2002) gives for the fast mode on the cold plasma limit $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,M}}& \mathrm{\approx }& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\mathrm{ℛ}\mathrm{\left(}{k}\mathit{,\omega }\mathrm{\right)}\\ & & \mathrm{\times }{\left[{\mathit{J}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{(}{\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }\right]}^{\mathrm{2}}{\mathit{P}}_{\mathit{xx,}\mathrm{M}}\mathrm{\left(}{k}\mathrm{\right)}\end{array}$(B.3)for high particle velocities, v_A ≪ v. As the polarisation of the fast mode on the cold plasma limit is the same as for the pseudo Alfvén wave, it is straightforward to show that the same form also holds for the pseudo Alfvén waves, with the difference in the dispersion relation. However, by using the magnetostatic limit with ω = 0 and again assuming no damping, we get for the resonance function ℛ(k,ω) = πδ(k_∥v_∥ + Ω), thus giving for the diffusion coefficient $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,}\mathrm{P}}& \mathrm{\approx }& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\mathit{\pi \delta }\left({\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{+}\mathit{n}\mathrm{\Omega }\right)\\ & & \mathrm{\times }{\left[{\mathit{J}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{(}{\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }\mathrm{)}\right]}^{\mathrm{2}}{\mathit{P}}_{\mathit{xx,}\mathrm{P}}\mathrm{\left(}{k}\mathrm{\right)}\mathit{.}\end{array}$(B.4)Unlike for the Alfvén waves, the Cherenkov resonance, n = 0, is nonzero for the pseudo Alfvén waves and has to be considered separately. In this case, the resonance function is $\begin{array}{ccc}{\mathrm{ℛ}}_{\mathit{n}\mathrm{=}\mathrm{0}}\mathrm{=}\frac{\mathit{\pi \delta }\left({\mathit{v}}_{\mathrm{\parallel }}\right)}{{\mathit{k}}_{\mathrm{\parallel }}}\mathrm{·}& & \end{array}$(B.5)This results in $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,P}}^{\mathit{n}\mathrm{=}\mathrm{0}}& \mathrm{=}& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{\delta }\left({\mathit{v}}_{\mathrm{\parallel }}\right)\\ & & \mathrm{\times }\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\frac{\mathit{\pi }}{{\mathit{k}}_{\mathrm{\parallel }}}{\mathrm{[}{\mathit{J}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{(}{\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }{\mathrm{)}}^{\mathrm{]}}}^{\mathrm{2}}{\mathit{P}}_{\mathit{xx,}\mathrm{P}}\mathrm{\left(}{k}\mathrm{\right)}\mathit{,}\end{array}$(B.6)which has a singularity at v_∥ = 0, and equals zero elsewhere. Consequently, this term is not used by our model.

For other terms in the sum over n, we again use the magnetostatic approximation ω = 0, thus giving the resonance condition $\begin{array}{ccc}\mathrm{ℛ}\mathrm{\left(}{k}\mathit{,\omega }\mathrm{\right)}\mathrm{=}\mathit{\pi \delta }\left({\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{+}\mathit{n}\mathrm{\Omega }\right)& & \end{array}$(B.7)and the pitch-angle diffusion coefficient $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,P}}& \mathrm{\approx }& \frac{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{n}\ne \mathrm{0}}\mathrm{\int }\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{k}\hspace{0.17em}\mathit{\pi \delta }\left({\mathit{k}}_{\mathrm{\parallel }}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{+}\mathit{n}\mathrm{\Omega }\right)\\ & & \mathrm{\times }{\left[{\mathit{J}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{(}{\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }\mathrm{)}\right]}^{\mathrm{2}}{\mathit{P}}_{\mathit{xx,}\mathrm{P}}\mathrm{\left(}{k}\mathrm{\right)}\mathit{.}\end{array}$(B.8)

Appendix C: Discretisation of the pitch-angle diffusion coefficients

For the numerical calculation of the pitch-angle diffusion coefficient, we consider the spectrum to be continuous in parallel direction, but discrete in the perpendicular direction, i.e. P_xx(k_x,k_y,k_z) = P_xx(hΔk,iΔk,k_z), with h,   i =  −M,...,M. In this manner, the integrals over k_x and k_y can be written as a sum, while the integral over k_z is evaluated using the delta function. Then, the diffusion coefficient for the Pseudo-Alfvén waves can be written as $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,P}}^{\mathit{n}\ne \mathrm{0}}& \mathrm{=}& \frac{\mathrm{2}\mathit{\pi }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\left|{\mathit{v}}_{\mathrm{\parallel }}\right|}\\ & & \mathrm{\times }\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }\mathit{,n}\ne \mathrm{0}}\hspace{0.17em}\sum _{\mathit{h,i}\mathrm{=}\mathrm{-}\mathit{M}}^{\mathit{M}}\mathrm{\Delta }{\mathit{k}}^{\mathrm{2}}{\left[{\mathit{J}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{(}{\mathit{v}}_{\mathrm{\perp }}{\mathit{k}}_{\mathrm{\perp }}\mathit{/}\mathrm{\Omega }\mathrm{)}\right]}^{\mathrm{2}}{\mathit{P}}_{\mathit{xx,}\mathrm{P}}\mathrm{(}\mathit{h}\mathrm{\Delta }\mathit{k,}\mathrm{i\Delta }\mathit{k,n}\mathrm{\Omega }\mathit{/}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{)}\mathit{.}\end{array}$(C.1)In this equation, (nΩ)/(μv) represents the parallel wavenumber at which δ(k_∥v_∥ + nΩ) is nonzero. Consequently, as turbulence data is available only at discrete wavenumbers, we define (nΩ)/(μv) = lΔk. Thus, with l restricted to integer values, we find the values of (v_∥ = μv) at which the value of D_μμ can be solved. For a particle with v = Ω/(mΔk), with m ≤ l ≤ M, the pitch-angle is given as μ = m/l. Applying this discretisation, we have

$\begin{array}{ccc}& & {\mathit{v}}_{\mathrm{\parallel }}\mathrm{=}\mathit{v}\hspace{0.17em}\mathrm{(}\mathit{m}\mathit{/}\mathit{l}\mathrm{)}\mathit{,}\\ & & {\mathit{v}}_{\mathrm{\perp }}\mathrm{=}\mathit{v}\hspace{0.17em}\sqrt{\mathrm{1}\mathrm{-}{\mathit{m}}^{\mathrm{2}}\mathit{/}{\mathit{l}}^{\mathrm{2}}}\mathit{,}\mathrm{and}\\ & & {\mathit{k}}_{\mathrm{\perp }}\mathrm{=}\sqrt{{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}}\mathrm{=}\mathrm{\Delta }\mathit{k}\hspace{0.17em}\sqrt{{\mathit{h}}^{\mathrm{2}}\mathrm{+}{\mathit{i}}^{\mathrm{2}}}\end{array}$(C.2)

and thus $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,P}}^{\mathit{n}\ne \mathrm{0}}\left(\frac{\mathit{m}}{\mathit{l}}\right)& \mathrm{=}& \mathrm{2}\mathit{\pi }\mathrm{\Omega }\left(\mathrm{1}\mathrm{-}\frac{{\mathit{m}}^{\mathrm{2}}}{{\mathit{l}}^{\mathrm{2}}}\right)\hspace{0.17em}\mathit{l}\\ & & \mathrm{\times }\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathit{M}\mathit{/}\mathit{l}}^{\mathit{M}\mathit{/}\mathit{l,}\hspace{0.17em}\mathit{n}\ne \mathrm{0}}\sum _{\mathit{h,i}\mathrm{=}\mathrm{-}\mathit{M}}^{\mathit{M}}{\left[{\mathit{J}}_{\mathit{n}}^{\mathrm{\prime }}\left(\frac{\mathrm{1}}{\mathit{m}}\sqrt{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{m}}^{\mathrm{2}}}{{\mathit{l}}^{\mathrm{2}}}\right)\left({\mathit{h}}^{\mathrm{2}}\mathrm{+}{\mathit{i}}^{\mathrm{2}}\right)}\right)\right]}^{\mathrm{2}}\\ & & \mathrm{\times }\frac{\mathrm{\Delta }{\mathit{k}}^{\mathrm{3}}{\mathit{P}}_{\mathit{xx,}\mathrm{P}}\left(\mathit{h}\mathrm{\Delta }\mathit{k,}\mathrm{i\Delta }\mathit{k,n}\hspace{0.17em}\mathit{l}\hspace{0.17em}\mathrm{\Delta }\mathit{k}\right)}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{·}\end{array}$(C.3)For the Alfvén waves, using the same discretisation, we get $\begin{array}{ccc}{\mathit{D}}_{\mathit{\mu \mu ,}\mathrm{A}}\left(\frac{\mathit{m}}{\mathit{l}}\right)& \mathrm{=}& \mathrm{2}\mathit{\pi }\mathrm{\Omega }\left(\mathrm{1}\mathrm{-}\frac{{\mathit{m}}^{\mathrm{2}}}{{\mathit{l}}^{\mathrm{2}}}\right)\hspace{0.17em}\mathit{l}\\ & & \mathrm{\times }\sum _{\mathit{n}\mathrm{=}\mathrm{-}\mathit{M}\mathit{/}\mathit{l}}^{\mathit{M}\mathit{/}\mathit{l}}\sum _{\mathit{h,i}\mathrm{=}\mathrm{-}\mathit{M}}^{\mathit{M}}{\left[\frac{\mathit{n}\hspace{0.17em}{\mathit{J}}_{\mathit{n}}\left(\frac{\mathrm{1}}{\mathit{m}}\sqrt{\mathrm{(}\mathrm{1}\mathrm{-}{\frac{{\mathit{m}}^{\mathrm{2}}}{{\mathit{l}}^{\mathrm{2}}}}^{\mathrm{)}}\left({\mathit{h}}^{\mathrm{2}}\mathrm{+}{\mathit{i}}^{\mathrm{2}}\right)}\right)}{\frac{\mathrm{1}}{\mathit{m}}\sqrt{\mathrm{(}\mathrm{1}\mathrm{-}{\frac{{\mathit{m}}^{\mathrm{2}}}{{\mathit{l}}^{\mathrm{2}}}}^{\mathrm{)}}\left({\mathit{h}}^{\mathrm{2}}\mathrm{+}{\mathit{i}}^{\mathrm{2}}\right)}}\right]}^{\mathrm{2}}\\ & & \mathrm{\times }\frac{\mathrm{\Delta }{\mathit{k}}^{\mathrm{3}}{\mathit{P}}_{\mathit{xx,}\mathrm{A}}\left(\mathit{h}\mathrm{\Delta }\mathit{k,i}\mathrm{\Delta }\mathit{k,n}\hspace{0.17em}\mathit{l}\hspace{0.17em}\mathrm{\Delta }\mathit{k}\right)}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{·}\end{array}$(C.4)

Appendix D: 512³ results

In this section we present the results of the investigation of the resolution by using a spatial grid of 512³ cells. Again, the MHD simulations were performed, first for the background turbulence, afterward with the peaked modes at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{8}$ and 24. The results are in conformance with the other setups. In particular, we observed the generation of higher harmonic wave modes again.

These harmonics are also visible for the peak at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ within the higher resolved simulation with a 512³ grid. In this simulation setup the number of active modes is eight times greater, which means the antialiasing edge is shifted by a factor two to k′ = 2π·86. Consequently, the higher harmonics at ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{48}$ and 72 are visible. The evolution of the peaks in both simulations is comparable to the 256³ grid simulations. A dominant energy transport towards high perpendicular wavenumbers is observed.

The test particle simulations led to comparable resonance patterns. As presented in Fig. D.2 the resonant interactions at μ_R = 0,1 and − 0.95 are strongly tilted and a significant amount of particles reach the maximum Δμ, as indicated by the sharp thresholds. The scattering coefficient is in this case again without any structure and hence not shown here. The stronger scattering is primarily caused by the higher wave modes, which are not truncated by the antialiasing anymore and hence contribute to the wave-particle interactions. The test particle simulations were performed for the decay stage of the peaks only because of their huge computational effort.

Additionally, an increase of the scattering rate at the ${\mathit{k}}_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{24}$ peak within the 512³ grid was observed. Consequently, the resonances are not as significant as in the smaller grid. When comparing Fig. 9 and D.3 it is harder to recognize the resonant structures.

As discussed previously, the 512³ gridsize particle simulation within the decay stage led to strong effectively tilted resonances and consequently a very unstructured D_αα curve. This can be observed in both Figs. D.4 and D.5. Thus, the results differ greatly from those given by SQLT.

Acknowledgments

We acknowledge support from Deutsche Forschungsgemeinschaft through grant SP 1124/3. S.L. additionally acknowledges support from European Framework Program 7 Grant Agreement SEPServer − 262773 and the computing time at Höchstleistungs Rechenzentrum Stuttgart (HLRS). T.L. acknowledges support from the UK Science and Technology Facilities Council (STFC) (grant ST/J001341/1)

References

All Tables

All Figures