© ESO, 2012

1. Introduction

The Sun emits energetic particles during coronal mass ejections (CME) and solar flares. The observed energies connected to these particles reach GeVs in extreme cases. At first, solar energetic particle (SEP) events were linked exclusively to solar flares (Forbush 1946), but once CMEs were first observed in the seventies, it became clear that they have a major role in the genesis of SEP events (Kahler et al. 1978). The CME-driven shocks are now believed to be the major acceleration agent during the strongest SEP events (e.g., Reames 1999), and the most plausible acceleration mechanism in operation is diffusive shock acceleration (DSA) (Bell 1978). For further reading into the matter observations and models of CMEs we recommend Forbes et al. (2006), and Chen (2011).

In DSA, particles are accelerated through repeated crossings of the shock compression front, each crossing giving a small boost to the particle energy. The shock crossings are mediated through interactions of particles with background waves. Furthermore, the particles amplify wavemodes. This yields a coupled system of waves and particles as described, e.g., in Vainio & Laitinen (2007). Since DSA is a resonant wave-particle process, it is now interesting to see if energy transfer between wavemodes will affect different particle energies. A detailed treatement of the DSA process itself is given in Drury (1983), and Schlickeiser (2002). To enable acceleration to the GeV energies through DSA, the upstream medium needs to be very turbulent. The ambient solar wind turbulence levels are generally too low to account for this acceleration mechanism to operate beyond the MeV regime (Vainio 2006), but acceleration to the highest energies can still proceed through the strong amplification Alfvén waves in the ambient medium through streaming instabilities driven by the accelerated particles themselves. Analytical (Lee 2005) and numerical models of the self-consistent particle acceleration in coronal plasma have been presented (Vainio & Laitinen 2007; Ng & Reames 2008), showing that the wave generation process is strong enough to account for the turbulence responsible for fast scattering of particles from one side of the shock to the other. The models also show that if particles are accelerated to hundreds of MeVs, the waves will grow to nonlinear amplitudes close to the shock.

It is crucial to understand the processes that govern the turbulent waves because the mechanism of DSA is strongly connected to the wave-particle interactions. Nonlinear Alfvén waves may interact with each other, which may lead to three important effects from the point of view of particle scattering: firstly, wave decay through three-wave interactions may limit the wave amplitudes in the shock environment; secondly, the wavenumber spectrum of the Alfvén waves may be altered so that the waves fall out of resonance with the particles; thirdly, the cross-helicity state of the resonant waves may also change, which affects the scattering center compression ratio of the shock and thus the accelerated particle spectrum (Vainio & Schlickeiser 1998). According to our knowledge, nonlinear wave transport models have not yet been applied to the SEP acceleration problem.

In this paper we will take the first steps towards understanding the nonlinear evolution of waves generated by energetic particle beams in the solar corona. We concentrate on a simplified model, where the beam-generated wave component is represented by a Gaussian peak in parallel wavenumber that follows the interaction of this spectral component with a quasi-isotropic background turbulence driven randomly in an incompressible magnetohydrodynamic (MHD) simulation. The turbulent plasma environment mimics the fast solar wind, where Alfvén waves are observed to be the dominant species (Telloni et al. 2009; Bruno & Carbone 2005), hence incompressibility can be assumed. It is especially interesting to see the energy transport in parallel and perpendicular direction. Taking into account the resonance condition, only energy transport parallel to the background magnetic field will alter the transport of particles with an energy different from the incident particle’s energy. On the other hand, most turbulence theories for incompressible plasmas predict perpendicular energy transport.

2. Numerical model

Despite MHD-turbulence has been studied for roughly 70 years after it was initiated by Hannes Alfvén, it is still a controversial topic. We focus on incompressible turbulence for which some promising progress has been made in recent years.

One commonly observed property is the characteristic energy spectrum E(k) following a power-law with a slope of  −5/3, which is commonly referred to as the Kolmogorov-type spectrum. Kolmogorov (1941) predicted this power-law for hydrodynamic turbulence by using dimensional analysis and scaling behaviour assuming isotropy. The basic system of turbulence evolution has also been explained by Kolmogorov. On large scales, i.e. small wavenumbers, energy is injected into the turbulent fluid. This energy decays by generating smaller structures up to the smallest scales where dissipation becomes dominant. Consequently, dissipation maintains the energy flow from small towards large k. This is also the reason for dividing the spectrum into driving-, inertial-, and dissipation range. Although Kolmogorov model of turbulence was first discussed in connection with neutral fluids, it seems to be valid in the magnetohydrodynamic case as well.

A different approach to Kolmogorov’s theories by Iroshnikov (1964) and Kraichnan (1965) assuming a local mean magnetic field and assuming Alfvén wave packets led to an exponent of  −3/2. The problem of the Iroshnikov-Kraichnan (IK) model is the assumption of isotropy because a background magnetic field will lead to a preferential direction in space caused by wave interaction resonances. The IK model implies resonant three wave interactions within a weak turbulent regime. In magnetised incompressible plasmas these interaction rates are empty (Sridhar & Goldreich 1994).

Goldreich and Sridhar (GS) describe anisotropic turbulence and distinguish between its weak (Sridhar & Goldreich 1994) and strong state (Goldreich & Sridhar 1995). Their assumption of strict separation between three and four wave interactions was controversially discussed and an intermediate state was introduced (Goldreich & Sridhar 1997). These theories describe Alfvénic turbulence evolution towards the perpendicular direction. In the weak turbulence regime four wave interactions are the underlying process, referring to the GS-framework. Due to the resonance condition energy transfer to parallel wavenumbers is not possible. It is still debated wheter four wave interactions are indeed the basic mechanism in weak turbulence is questionable. In recent theories the intermediate turbulence (Goldreich & Sridhar 1997), which is based on three wave interactions, replaces the weak four wave interaction model (Lithwick & Goldreich 2003). However, strong turbulence is dominated by nonresonant three wave interactions, which leads to an anisotropic energy-cascade in the perpendicular direction. Parallel evolution is not caused by cascading. One of the main achievements of the Goldreich and Sridhar theory is that is explains the Kolmogorov-type energy spectrum for an anisotropic regime. This could explain the observed  −5/3 slope in parts of the solar wind (Bruno & Carbone 2005) where Kolmogorov-theory is not applicable.

The region of the heliosphere we are interested in is within the weak turbulence regime, with magnetic field fluctuations defined as (1)with a mean value of  ⟨ dB ⟩  = 0, which leads to  ⟨ B ⟩  ~ B₀. It is observed that the solar wind magnetic fluctuations decrease as dB² ∝ r^-3, while the background field decreases as (Bavassano et al. 1982; Bruno & Carbone 2005). Consequently, the dB/B₀ ratio within the heliosphere ratio is increasing with distance to the Sun (Hollweg et al. 2010).

A remark on the notation, the magnetic background field B₀ is defined towards the z-direction within our simulations and hence is also noted as B₀e_z. The parallel direction is, therefore, defined as the z-directions 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 x and y.

For small perpendicular wave numbers the transport in perpendicular direction is dominating until k _⊥ v _⊥  is of the same order as the Alfvén cascading time k_∥v_A. This means that in addition to the perpendicular cascade, a cascade in the parallel direction will occur as well and energy will be transferred towards smaller parallel spatial scales. Accordingly, the parameter (2)where v_A is the Alfvén velocity, is of the order of unity. This state is called the critical balance and was first introduced by Goldreich & Sridhar (1995). In this state the linear wave period of the Alfvén waves are comparable to the intrisically nonlinear timescale. If ζ ~ 1, the fluctuations become more correlated along the parallel direction, up to l_∥ ~ v_A/(k _⊥ v_⊥) as indicated by Eq. (2). Then the turbulence is clearly within the strong regime. This means that the fluctuations become comparable to B₀ and the nonlinear term is not small anymore (Perez & Boldyrev 2008).

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. In situ measurements of the energy spectrum (Tu et al. 1989) agree with this fact. To simulate conditions within the turbulent heliospheric plasma, the research group at the University of Würzburg has developed a simulation code, Gismo.

Gismo is an incompressible pseudospectral MHD-code that is fully parallelised and capable of efficiently using massive computing clusters. The basis of the simulation software is to solve the following set of incompressible MHD-equations: (3)where is the normalised magnetic field, u is the fluid velocity and ϱ is the constant mass density. The diffusion coefficient related to viscous and Ohmic dissipation is denoted by ν_v and η. A common approach in pseudospectral methods is to amplify the diffusion term by a power of h – resulting in hyperdiffusivity. This artificial enhancement of the dissipation is necessary to reach a saturated state of turbulence within a reasonable timescale. It is a methodic problem, since pseudospectral approaches do not strongly suffer from dissipative numerical effects. The only intrisic energy loss of the system is caused by antialiasing, which we discuss below. Furthermore, the parameter ν is introduced as a global diffusivity with η = ν_v ≡ ν. Hence magnetic resistivity and viscous damping are not distinguishable anymore. This is the case for the magnetic Prandtl number of the order of unity, which is valid within the regime of Alfvén wave turbulence where an equipartion between magnetic and kinetic energy can be assumed (Bershadskii 2002; Bigot et al. 2008). The pressure term ∇P fulfills the closure condition for incompressibility (Maron & Goldreich 2001) (4)These equations are solved in Fourier space by using pseudospectral methods that lead to the componentwise equations (5)where the tilde-notation stands for quantities in Fourier space. The components of the wavevector are written as k_α.

In the incompressible regime of a magnetised plasma the MHD-turbulence consists of only two types of waves, which propagate 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 wavevector 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_∥)². 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). However the damping weakens in strong turbulence, but according to Goldreich & Sridhar (1995), the wave generation of pseudo-Alfvénic wavemodes is only possible via three-wave interactions by two shear wavemodes. Barnes damping is important for high-β plasmas. Since this is not fulfilled for the solar corona, the role of pseudo waves should not be ignored.

Because the model consists only of these two wave types, it is suitable to use a description with Alfvénic waves moving either forwards or backwards. This is achieved by introducing the Elsässer variables (Elsässer 1950) (6)and transforming the Eqs. (5) into a suitable form of (7)Obviously, the nonlinearities of Eqs. (7) that describe the turbulent behaviour of the MHD-plasma cannot be solved in Fourier space. 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 a MPI-parallelised fast Fourier transfomation based on FFTW3 (Pekurovsky 2011).

One 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 Orszags 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 each step, immediately before calculating the nonlinearities and, accordingly, calculating the 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.

Gismo is capable of using different forward-in-time schemes, namely Euler and Runge-Kutta second as well as fourth order. All the simulations in this paper has been performed by using RK-4.

To mimic an SEP-event, understanding the underlying mechanisms is crucial. Before we simulated particle scattering at modified field-fluctuations, we focussed on the evolution of amplified wavemodes within the background turbulence. Since streaming particles ejected by the Sun, e.g., protons from coronal mass ejections, will not sharply amplify a discrete mode but a broader interval, a Gaussian distributed shape was assumed. The streaming instability has the highest growth rate in the background field direction. Consequently, we assumed that only purely parallel-propagating Alfvén waves are modified.

The Alfvén wave generation mechanism by SEPs is not be investigated in detail here. The driving mechanism assumed for our simulations is the streaming instability. The estimatation of the wave growth rate is described in Vainio (2003). The streaming instability is caused by energetic proton scattering off the interplanetary Alfvén waves. During the the scattering process the particle changes its pitch angle cosine by Δμ while its momentum in the wave frame remains constant. Thus the particles’ energy in the plasma frame is changed by v_ApΔμ, consequently also happens to the Alfvén wave energy due to energy conservation. Another important instability in the solar corona is the electrostatic instability. This is caused by an electron current as well as by 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).

These parallel peaked modes are influenced by dissipation, diffusion, and convection. An illustration of this is shown in Fig. 1. Regarding the evolution within the Fourier space, dissipation will damp the wavemode, hence lower the maximum without altering the position or broadness of the peak. The dissipation of wavemodes is caused by the spatial diffusion term in Eqs. (7). Convection in Fourier space would shift the position of the peak. If diffusive transport in Fourier space were the dominant mechanism, it would result in a broader energy distribution. The dynamics of convection and diffusion lie within the nonlinear terms, therefore one cannot distinguish exactly between the responsible terms.

To investigate the effects of spatial diffusion on the peaks, we solved the dissipation equation in wavenumber space (8)and calculated its dissipation coefficient (9)When the spatial diffusion is the only dissipative effect in wavenumber space, this wavenumber dissipation coefficient would be the spatial diffusion coefficient. We emphasise that diffusion in the wavenumber space is a different process and is hence explicitly distinguished from spatial diffusion. In the discussed context above the spatial diffusion is clearly connected to the wavenumber dissipation. We concentrated our investigation on those wavemodes whose peak energy is initiated only. For other modes this approach is not feasible because the spatial diffusion is not necessarily the dominant process on an arbitrary wavemode.

Fig. 1 Possible mechanisms that might influence the evolution of an amplified wavemode embedded in a turbulent plasma with a Kolmogorov-type power-law energy spectrum.

3. Simulation setup

To simulate the turbulent plasma in which the SEP-event is set, we performed the following type of magnetohydrodynamic turbulence.

We used an anisotropically driven turbulence with a driving range in k-space up to the first five numerical wavenumbers in perpendicular () and 15 in parallel direction (). A remark to the notation, the wavenumber is in general defined as k = (2πn)/L where n stands for the numerical grid position. For simplicity we used the normalised wavenumbers k′ = (2π·n) throughout. The anisotropy was chosen for two reasons. First, to mimic the preferential direction of the solar wind, where particles that radially stream 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, in which we are interested. Second, a slab-component of SW turbulence is observed also at small scales in 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. 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 random 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 dB/B₀ ratio of roughly 10^-2. 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-waves evolution is strongly suppressed. In this initial range energy is injected at discrete times, which leads to a saturated turbulence – an equilibrium between dissipation and injection. The spatial resolution is 256³ gridpoints, resulting in 128³ points in k-space of which  | k^′ |  = 2π·42 wavemodes are active modes that remain unaffected by (anti)aliasing. The hyperdiffusivity coefficient was set to h = 2.

The simulations of the turbulent background plasma were performed assuming an outer scale of L_scale = 3.4 × 10⁸   cm. This value was estimated using the growth rate from Vainio (2003)(10)with the proton cyclotron frequency ω_cp, the proton speed v_p, the Lorentz factor γ, the proton number density n_p, μ the pitch angle cosine, and f_p the proton distribution function. Here the resonance condition for the nth order of interaction (11)(cf. Schlickeiser 1989) was used, where ω is the wave frequency, and k_∥ its parallel wavenumber. Ω is the particle’s gyrofrequency and v_∥ its parallel velocity component. Note that Eq. (10) is only valid for purely parallel waves and n =  ± 1. Orders of  | n |  > 1 can only be generated by oblique waves. The perpendicular components of the wave would then modify the scattering process by nonvanishing Bessel functions. This is discussed in detail by Schlickeiser (2002). We used peaks at k = 2π·8 and k = 2π·24, which represent proton energies of E ≈ 64MeV and E ≈ 7MeV respectively. Using the resonance condition, this leads to a length scale of (12)The Alfvén speed was assumed to be v_A = 1.2 × 10⁸cms^-1, which leads along with a particle number density of 10⁵cm^-3 to a background magnetic field B₀ = 0.174G. These values resemble the solar wind environment at three solar radii (Vainio et al. 2003), where particle acceleration by CME-driven shocks is strongest.

The discretisation of the timestep is stable for values up to Δt′ = 1 × 10^-11 in numerical units or Δt′ = 3.4 × 10^-3s for the background turbulence.

If the background plasma simulation reaches the saturated state, a Gaussian-distributed energy peak with purely parallel k = ke _∥  is injected. We chose two different positions of the peak in wavenumber space. To investigate the physics of an SEP-event, a wavenumber of k_∥ = 1.5 × 10^-7cm^-1 was used. This corresponds to a numerical wavenumber of 8, which is still within the driving range of the turbulence. The injection at smaller scales was represented by a peak at k_∥ = 4.4 × 10^-7cm^-1. This value lies at the numerical position 24, which is roughly between the maximum driven wavemode and the anti-aliasing truncation edge (which was at 43). We injected the SEP-energy gradually over a certain time interval to develop a realistic scenario. Multiple situations were explored, by using simulations with peaks at either position, with either large (growth rate Γ₁) or small (growth rate Γ₂) total amplitude of the Gaussian at the final driving step. Because the velocities increase near the peaks, the discretisation of the timestep had to be decreased to values of Δt′ = 5 × 10^-12 in numerical units or Δt′ = 1.7 × 10^-3s to sustain stability.

To allow even more diverse case studies, four different initial conditions were used, as described in Table 1. In each setup a complete evolution of the background turbulence was simulated. The first setup SI uses the standard parameters as described above. The main parameters of interest are the resistivity of the plasma and the background magnetic field. The results in changing these will reveal the effects on the mechanisms described in Fig. 1. An increased resistivity compared to simtype SI is approached in SII. A higher value for ν is expected to make a difference in the spatial diffusion behaviour and would make wavenumber dissipation more dominant to the other transport. As indicated in Fig. 1, this would lead to a significant damping of the peak. The dissipation range of the background turbulence will likely be increased by this parameter as well. The third setup SIII has a magnetic field increased by a factor of 10. This is to examine the influence of a more anisotropic turbulence because the perpendicular evolution should be much stronger according to Goldreich & Sridhar (1995). In general, these values may only be achieved in magnetic clouds, but it gives valuable information on the mechanisms of turbulent transport. The high resistivity is necessary because of stability problems with the accompanying high Alfvén speeds. The last variation of the scenario, which was used in SI, is a strongly decreased magnetic background field B₀, see simtype SIV. The aim of this artificial scenario is to investigate strong turbulence at ζ ≈ 1.

4. Results

Fig. 2 Magnetic energy spectrum of the simulated background turbulence setups. The plot was made by total integration within the Fourier space.

The evolution of the anisotropic background turbulence was simulated up to 30 Alfvén wave crossing times, which corresponds to a simulation of 85 s in physical time. At this point, the turbulence has reached a saturated state and a Kolmogorov-type power-law has evolved over a wide range of wavenumbers (see Fig. 2). According to Sridhar & Goldreich (1994) and Goldreich & Sridhar (1995), the 5/3-spectrum is dominant for perpendicular k _⊥ , whereas the parallel evolution is significantly slower. Our simulations confirm this behaviour very clearly. Note that the total spectra in Fig. 2 deviate from the Kolmogorov shape, especially at high wavenumbers, because the shape of the parallel spectrum is not 5/3. As expected, the spectra are very sensible to ν. A factor of ten increases the dissipation range drastically. This is also due to the hyperdiffusivity we used where higher wavenumbers are damped by higher power of k (dissipation term  ∝ k⁴) (see Sect. 2). The dB/B₀ ratio of the developed turbulence is about 10^-2. The magnetic field fluctuations are defined in Fourier space by (13)The influence of the turbulence strength on the energy transport is an important aspect for the peak simulations. We show this below.

Fig. 3 Simulation setup SI: time evolution of a Gaussian-distributed amplification at numerical wavenumber 24 () at the lower growth rate Γ1. The spectrum is a one-dimensional cut along the parallel wavenumber axis where the peak is located. The peak is clearly shifted towards smaller k∥.

Fig. 4 Simulation setup SI: one-dimensional energy spectrum E(k∥) of the peak at numerical wavenumber 24 () with higher growth rate Γ2. The influence of the diffusion process is significant because the peak is broadening during time evolution. Adjoining maxima develop, e.g. at and 38, highlighted by red circles.

Once the background plasma was simulated, a peak is driven over a time period of ca. 1.7 s. An exemplary time evolution of the peak at normalised wavenumber 24 is shown in Figs. 3 and 4. This is a one-dimensional spectrum of the magnetic field energy in numerical units that shows a cut along the parallel axis. The starting time t_start corresponds to the end of the driving interval, hence the maximum amplification of the wavemode. The timesteps are shown in Table 2. For subsequent use the times t_mid and t_end are introduced. Note that the decay time interval of the excited mode at is significantly longer compared to the peak. This is again because of the hyperdiffusivity, which damps higher modes more strongly because the dissipation term is  ∝ k⁴. The time t_end is defined at the state where the peak has lost nearly its complete energy within the parallel direction. We took this as the final point of the energy diffusion or dissipation.

Both peaks at show broadening of its shape and shifting from the initial value to within 17 s, while peaks at are only slightly shifted to  s within 17 s, but the broadening is clearly visible, as shown in Fig. 5.

Fig. 5 Simulation setup SI: time evolution of a Gaussian-distributed amplification at numerical wavenumber 8 () at the lower growth rate Γ1. The spectrum is a one-dimensional cut along the parallel wavenumber axis where the peak is located. Only a slight shift towards small k∥ is observed. The broadening is clearly visible, especially on the flanks of the Gaussian curve.

The next step is to investigate the development of the amplified wavemodes. If the evolution of the peak were solely governed by spatial diffusion, the dissipation coefficient D of Eq. (9) would stay constant in time. Therefore the change in peak energy was measured. Two intervals were used for the calculation of D, denoted as time-interval τ₁ and τ₂. These intervals are different for both peak positions because of the faster decay at high wavenumbers. At the time intervals are τ₁ = 8.5 s and τ₂ = 17 s, whereas at the intervals τ₁ = 1.7 s and τ₂ = 3.4 s are used. The results of the diffusion coefficients are given in Table 3.

To relate the growth rates Γ_1/2 with the total resulting amplitude of the Gaussians of the simulations SI-SIII, the energy was measured and compared to the background at timestep t_start each. The results are presented in Table 4.

To investigate the direction of the peak evolution in the parallel and perpendicular directions, two-dimensional contour plots were produced. Figure 6 shows the time evolution of a peak at normalised . The two-dimensional spectrum is a contour plot of the power spectral density of the magnetic field that was calculated by cylindrical integration in k-space. The single contours are scaled logarithmically.

During the evolution, higher harmonics of the initial peak arise. To investigate these in greater detail, we measured the energy of the initial peak and its first harmonic. The result is shown in Fig. 7. Note that the generation of these modes starts at higher perpendicular wavenumbers (see Fig. 6, most picture left, first harmonic at k _⊥  ≈ 15) and not at purely parallel k.

In addition to the observed higher harmonics, another basic result of the simulation type SI is the discrete generation of other modes along the k_∥ axis especially at higher amplitudes. As seen in Fig. 4, adjoining maxima develop next to the main Gaussian curve, e.g. at positions and 38. The way the peak develops appears to vary with the amplitude, the peak with Γ₁ generates clearly fewer other modes than the larger peak with Γ₂. Both peaks show a significant change of their original wavenumber position and a clear broadening. The two-dimensional spectra reveal a strong perpendicular development in k-space especially at large  | k | . During the decay the evolution becomes a little more isotropic, but the preferential direction of evolution is clearly perpendicular.

The importance of the amplitudes Γ_1/2 is also visible in Fig. 8. The development of each peak is very different. An interesting result is a more dominant evolution of the peak with the high growth rate Γ₂ that is towards smaller k_∥ directed. In the peak-dominated region the turbulence seems to increase. Therefore the ζ parameter (Eq. (2)) is of interest. Figure 9 shows a map of values of ζ =  [0.01···0.15] , i.e. near the critical balance. The same plot for the lower growth rate would be empty. This also indicates that the high ζ values along the k _⊥ -axis in Fig. 9 stem from interactions with the peak.

The peaks at remain much longer than the higher modes, see Table 2. This is due to the higher k on which the dissipation process depends. Furthermore, the hyperdiffusivity damps higher modes more strongly ( ∝ k⁴).

The simulation of the same peaks in setup SII with higher ν reveals one key feature. The shift towards smaller k_∥ positions occurs on both, and 24 peaks and is stronger compared to SI. The shifting is significant, especially at , which is shown in Fig. 11. Within 2.55 s the original position changes from to 22. The peak amplitude is again important for the evolution. We observed that higher growth rates lead to an isotropic evolution, whereas the lower growth rates show strongly perpendicular development. The effective peak energy is lower and consequently the energy transport towards higher wavenumbers is more restricted. There are also fewer of higher harmonics. This can be observed by direct comparison of the left plot in Fig. 10 with the middle plot in Fig. 6. As expected, the decay of the energy is faster compared to SI (see Table 2).

The simulation setup SIII reveals very strong perpendicular evolution for all peaks. Two examples are given in Fig. 12. The growth rates seem not to have a strong influence on the development. Only more higher harmonics of peak are visible with Γ₂. The time of total energy decay of the peaks is slightly longer than in the simulations SII.

In the last simulation setup SIV the magnetic background field was reduced by two orders of magnitude. The resulting dB/B₀ ratio is of the order of 10 and the turbulence development is highly isotropic. The peaks at and both at growth rate Γ₂ show very interesting features.

In addition to the typical peak evolution and generation of higher harmonics, other structures also arise at high k _⊥ . As shown in the Figs. 13 and 14, both peak positions generate these structures, but at various positions. During the development of the peak at position theses structures arise perdominantly at high locations, e.g. (2π·10,2π·18), (2π·10,2π·38), (2π·18,2π·32) and (2π·24,2π·18), see Fig. 13. The structures evolving with the simulations are instead located at middle k_⊥ but low k_∥, e.g. (2π·5,2π·13), (2π·5,2π·17) and (2π·8,2π·22), see Fig. 14. In contrast to the higher harmonics, these structures are not necessarily integral multiples of the initial peaked mode. A map of the critical balance parameter ζ was calculated for both plots. Around the wavenumber of the peaks and along the k _⊥  axis this parameter is of order 1. Interestingly, during the simulation this parameter increases, especially at the k_⊥ axis. In contrast to the other setups, the higher harmonics are located along the k_∥ axis at low or zero k_⊥. A significant shift along this axis towards smaller parallel wavenumbers is observed, e.g. as shown in Fig. 15.

We discuss possible explanations for these phenomena below.

Fig. 6 Two-dimensional magnetic energy spectra of a peak at normalised wavenumber 8. Red regions are at higher energies compared to the blue ones. The parallel and perpendicular wavenumbers are given as absolute values. The time development is shown for mid-drive (Δt ≈ 0.85s), max-peak (Δt ≈ 1.7s) and the decay 17 s after the driving. The colours of the contours were normalised for comparison between the three plots. The colours indicate the logarithm of the total spectral energy. The simulation setup SI was used.

Fig. 7 Energy dependency of the first harmonic on the initial peak energy. A fit resulted in a quadratic function. The simulation setup SI was used.

5. Discussion

As discussed in Sect. 2, there are three possibilities how an excited wavemode can develop: through diffusion, dissipation and convection.

A general conclusion from our simulations is that the dynamics of these mechanisms are strong at high wavenumbers. Especially for the dissipation process this is not unexpected because of it strongly depends on the wavenumber. The hyperdiffusivity might amplify this effect because of the higher power in k. The Figs. 3 and 4 clearly show a rapid dissipation of energy because the peak loses amplitude very fast. Also Table 2 indicates a decay in short timescales at high wavenumbers. However, a broadening also arises at the flanks of the Gaussian distribution. The broadening is strong for the higher growth rate Γ₂. Within a time interval of 1.7 s after the driving range the FWHM is increased by roughly 70% at Γ₂, whereas the peak with smaller amplitude is broadened by ca. 15%. A possible explanation is the equilibrium between energy- and enstrophy cascade, which causes a similar flow of energy to large and small wavenumbers (Mininni & Pouquet 2009).

The convective mechanism shifts the maximum of the peak. Convection is a slow process compared to dissipation and diffusion within the simulated regime. Nevertheless we were able to observe it within our simulations, e.g. in Fig. 3. The transport is towards smaller wavenumbers, which indicates an enstrophy cascade. This effect is more typical for two-dimensional plasmas. The MHD-development in anisotropic plasmas is mostly effectively two-dimensional. This lead to inverse energy cascades as well as upscaling enstrophy cascades. These mechanisms generate larger vortices and consequently transfer energy to smaller wavenumbers.

Fig. 8 Simulation setup SI: comparison between growth rates Γ1 and Γ2 at peak position at maximum drive time. Although the effective energy input varies by a factor 100 (see Table 4) another transport mechanism becomes dominant. The smaller peak (Γ1) develops dominantly in the perpendicular direction while the evolution of bigger peak (Γ2) is more isotropic and tends towards smaller k∥. The colours indicate the logarithm of the total spectral energy.

Fig. 9 Map of the critical balance parameter ζ for the right-hand plot in Fig. 8. The contours linearly represent values between 0.01 and 0.14 of the integral values of ζ(k∥,k⊥). The peak structure and values near the k⊥-axis are clearly visible.

Fig. 10 Two-dimensional peak evolution in setup SII at maximum drive time. The higher harmonics develop at lower k ⊥  compared to SI. The higher growth rate (Γ2, right panel) shows a strong parallel evolution. The smaller amplitude (Γ1, left panel) develops dominantly towards higher k ⊥ . The colours indicate the logarithm of the total spectral energy.

Fig. 11 Simulation setup SII: time evolution of a Gaussian distributed amplification at at the smaller growth rate Γ1. Again a 1D cut spectrum along k∥ axis is shown. The shift of the peak is stronger compared to SI.

Fig. 12 Evolution of the two peak positions at tmid in setup SIII. Very strong perpendicular development in all SIII simulations is observed. The edge of the parallel driving range at is clearly visible in the right-hand panel.

Fig. 13 Two-dimensional peak evolution in setup SIV. Left: development of the peak at with growth rate Γ2. The figure presents the state 1.28 s after the maximum driven peak (t0). Right panel: corresponding map of the critical balance parameter ζ. Each contour represents integral values above ζ = 0.1. The colour scaling is linear.

Fig. 14 Two-dimensional peak evolution in setup SIV. Left: development of the peak at with growth rate Γ1. The figure presents the state 0.68 s after the maximum driven peak (t0). Right panel: corresponding map of the critical balance parameter ζ. Each contour represents integral values above ζ = 0.1.The colour scaling is linear.

Fig. 15 Comparison of the convective peak shifting between the setups with changing background field B0 SI, SIII and SIV. The strong B0 in SIII causes no observable shifting at all, whereas the weak B0 leads to a significant change of the original peak position. One possible explanation is an enstrophy cascade.

To investigate this behaviour in more detail, we concentrated on the peak with growth rate Γ₁ within the setups of changing magnetic background field SI, SIII and SIV. The comparsion of these setups is shown in Fig. 15. We observe a slight shift of the in setup SI after 17s of the maximum drive from position 8 to 7.7. The simulation SIII with increased B₀ shows no transport of the peak. After the same time interval it is still at position 8. The evolution in setup SIV with a small B₀ is very strong. After 17s the peak has shifted by roughly 1.5 grid positions from to 6.5. The development of an enstrophy cascade explains this behaviour. The strong magnetic field effectively makes the MHD-evolution one-dimensional. As B₀ decreases, the evolution becomes less restricted to the magnetic background field. Consequently, the enstrophy cascade increases. Figure 16 clearly supports this explanation. The enstrophy was calculated by (14)All two-dimensional spectra show a strong evolution in the perpendicular direction. This is consistent with the theories of Sridhar & Goldreich (1994) and Goldreich & Sridhar (1995) within a turbulent plasma, as explained in Sect. 2. The process of this perpendicular evolution is caused by the energy cascade. Mainly Fig. 12 shows strong perpendicular behaviour, whereas the evolution clearly becomes more isotropic in SIV (see Figs. 13 and 14). This is due to the increasing strength of the turbulence for the cascade, which is expressed by the dB/B₀ ratio.

The dissipation coefficients presented in Table 3 do remain constant especially for all simulations with Γ₂ at . This implies that spatial diffusion is the dominant process for these simulations. The most significant change of the dissipation coefficient is at for setup SII, between τ₁ and τ₂. We interpret this as a nonlinear effect where wavemodes are generated, which in turn triggers the cascade. This influences the energy of the initialised Gaussian very strongly. As expected, the dissipation coefficients for setup SII are stronger for the  peaks. The spatial diffusion coefficient is connected with the wavenumber dissipation process via ν. This is because the spatial diffusion process is the only mechanism that leads to energy losses in the k space (see Sect. 2). This is at least valid for wavemodes far below the antialiasing edge, where energy is artificially removed as well. A connection of wavenumber diffusion to ν is observed in Fig. 10 where the Gaussian shape is significantly broadened, which is caused by the diffusion process.

Fig. 16 Time evolution of the enstrophy for the simulation setups SI, SIII and SIV. The strong magnetic background field of SIII prevents the enstrophy from developing while a clear change in SIV is visible.

The observed higher harmonics could resemble a three-wave process k₈ + k₈ → k₁₆. This is supported by Fig. 7. The energy dependency between peak and the next higher harmonic is quadratic. This is also the case for three-wave interactions (Sagdeev & Galeev 1969). As pointed out before, this is forbidden for Alfvén waves by wave interaction processes. Just oppositely directed wavepackages can collide, hence momentum conservation would be violated by the proccess described above. This wave interaction can only take place with oppositely directed waves of the background plasma. This must also be true since a cross-check simulation of the peak without turbulent background did not show these harmonics. Another explanation is given by Galinsky et al. (1997): Alfvén waves interact with themselves, which leads to wave steepening.

Investigating the strong turbulence evolution in simulation SIV shows unexpected structures at high k _⊥  (see Figs. 13 and 14). This effect might be caused by the critical balance of strong turbulence within the Goldreich and Sridhar description. This requires the parameter ζ to be  ~1, which is the case for locations in the vicinity of the peaks and along the perpendicular axis. Mainly the peak seems to amplify the region at and . Values of ζ > 0.1 lead to the development of these structures, but it is not possible to conclude whether the stuctures arise first and then generate higher values of ζ, or vice versa. Nevertheless, the turbulence strength increases within these regions, which agrees with Goldreich & Sridhar (1995) where ζ is assumed to become unity during the turbulence evolution. In addition to setup SIV, SI also shows this behaviour as presented in Figs. 8 and 9. More investigation is needed to clarify the underlying processes.

6. Conclusion

We analysed the evolution of waves generated by proton beams in a turbulent medium. This evolution may play an important role in diffusive shock acceleration in the heliosphere. Our study has revealed that three different processes as sketched in Fig. 1 are taking place.

The most interesting question is wheter wavemodes excited by particles of a certain energy yield wavemodes that interact with other particle energies. The observed shifting of the initial parallel wavenumber position towards smaller k_∥ influences the particle acceleration at these modes. As shown in Eq. (11) and in the corresponding section, this means that higher energetic particles can be accelerated because the wavemode develops towards higher spatial scales (Vainio & Laitinen 2007). The shift towards smaller k_∥ is indeed fairly minor in this simulation. This is because of the injection of energy at only one single wavenumber and the limited simulation time. We also expect an effect through nonlinear amplification: Waves at smaller k_∥ accelerate higher energetic particles, again injecting energy at lower k_∥. On the other hand a strong evolution towards high k_∥ has also been observed in terms of development of higher harmonics of the initialised mode. Consequently, particle acceleration at lower energies is also possible. It should be noted, however, that this is not consistent with isotropic diffusion of energy in wavenumber space as assumed in Vainio & Laitinen (2007). This means that the previous models of the wave-particle system will have to be updated accordingly to account for the strong dependence of energy transport on the direction in wavenumber space. Especially the development of the higher harmonics contradicts an identical forward and backward cascade.

The strong perpendicular evolution of the peak initialised with purely parallel propagating waves causes higher orders ( | n |  > 1, see Eq. (11)) of resonance between solar particles and the amplified mode. This is because Eq. (10) has to be modified in this case (Schlickeiser 2002).

Owing to limited computational power we have not been able to investigate the effect of critical balance in detail. But we note that under certain conditions (k_∥/k _⊥ , amplitude) the evolution may be governed by the critical balance.

Acknowledgments

We express our graditude to Rami Vanio, Timo Laitinen and Markus Battarbee for cooperation and their contributions to this work. We acknowledge support from the Deutsche Forschungsgemeinschaft through grant SP 1124/3. S.L. additionally acknowledges support from the European Framework Program 7 Grant Agreement SEPServer – 262773. We thank the anonymous referee for her/his detailed comments, which improved the paper significantly.

References

All Tables

All Figures

	Fig. 1 Possible mechanisms that might influence the evolution of an amplified wavemode embedded in a turbulent plasma with a Kolmogorov-type power-law energy spectrum.
In the text

	Fig. 2 Magnetic energy spectrum of the simulated background turbulence setups. The plot was made by total integration within the Fourier space.
In the text

	Fig. 3 Simulation setup SI: time evolution of a Gaussian-distributed amplification at numerical wavenumber 24 () at the lower growth rate Γ1. The spectrum is a one-dimensional cut along the parallel wavenumber axis where the peak is located. The peak is clearly shifted towards smaller k∥.
In the text

	Fig. 4 Simulation setup SI: one-dimensional energy spectrum E(k∥) of the peak at numerical wavenumber 24 () with higher growth rate Γ2. The influence of the diffusion process is significant because the peak is broadening during time evolution. Adjoining maxima develop, e.g. at and 38, highlighted by red circles.
In the text

	Fig. 5 Simulation setup SI: time evolution of a Gaussian-distributed amplification at numerical wavenumber 8 () at the lower growth rate Γ1. The spectrum is a one-dimensional cut along the parallel wavenumber axis where the peak is located. Only a slight shift towards small k∥ is observed. The broadening is clearly visible, especially on the flanks of the Gaussian curve.
In the text

Fig. 6 Two-dimensional magnetic energy spectra of a peak at normalised wavenumber 8. Red regions are at higher energies compared to the blue ones. The parallel and perpendicular wavenumbers are given as absolute values. The time development is shown for mid-drive (Δt ≈ 0.85s), max-peak (Δt ≈ 1.7s) and the decay 17 s after the driving. The colours of the contours were normalised for comparison between the three plots. The colours indicate the logarithm of the total spectral energy. The simulation setup SI was used.

In the text

	Fig. 7 Energy dependency of the first harmonic on the initial peak energy. A fit resulted in a quadratic function. The simulation setup SI was used.
In the text

Fig. 8 Simulation setup SI: comparison between growth rates Γ1 and Γ2 at peak position at maximum drive time. Although the effective energy input varies by a factor 100 (see Table 4) another transport mechanism becomes dominant. The smaller peak (Γ1) develops dominantly in the perpendicular direction while the evolution of bigger peak (Γ2) is more isotropic and tends towards smaller k∥. The colours indicate the logarithm of the total spectral energy.

In the text

	Fig. 9 Map of the critical balance parameter ζ for the right-hand plot in Fig. 8. The contours linearly represent values between 0.01 and 0.14 of the integral values of ζ(k∥,k⊥). The peak structure and values near the k⊥-axis are clearly visible.
In the text

	Fig. 10 Two-dimensional peak evolution in setup SII at maximum drive time. The higher harmonics develop at lower k ⊥  compared to SI. The higher growth rate (Γ2, right panel) shows a strong parallel evolution. The smaller amplitude (Γ1, left panel) develops dominantly towards higher k ⊥ . The colours indicate the logarithm of the total spectral energy.
In the text

	Fig. 11 Simulation setup SII: time evolution of a Gaussian distributed amplification at at the smaller growth rate Γ1. Again a 1D cut spectrum along k∥ axis is shown. The shift of the peak is stronger compared to SI.
In the text

	Fig. 12 Evolution of the two peak positions at tmid in setup SIII. Very strong perpendicular development in all SIII simulations is observed. The edge of the parallel driving range at is clearly visible in the right-hand panel.
In the text

	Fig. 13 Two-dimensional peak evolution in setup SIV. Left: development of the peak at with growth rate Γ2. The figure presents the state 1.28 s after the maximum driven peak (t0). Right panel: corresponding map of the critical balance parameter ζ. Each contour represents integral values above ζ = 0.1. The colour scaling is linear.
In the text

	Fig. 14 Two-dimensional peak evolution in setup SIV. Left: development of the peak at with growth rate Γ1. The figure presents the state 0.68 s after the maximum driven peak (t0). Right panel: corresponding map of the critical balance parameter ζ. Each contour represents integral values above ζ = 0.1.The colour scaling is linear.
In the text

	Fig. 15 Comparison of the convective peak shifting between the setups with changing background field B0 SI, SIII and SIV. The strong B0 in SIII causes no observable shifting at all, whereas the weak B0 leads to a significant change of the original peak position. One possible explanation is an enstrophy cascade.
In the text

	Fig. 16 Time evolution of the enstrophy for the simulation setups SI, SIII and SIV. The strong magnetic background field of SIII prevents the enstrophy from developing while a clear change in SIV is visible.
In the text