SIMULATIONS OF ENERGETIC PARTICLES INTERACTING WITH DYNAMICAL MAGNETIC TURBULENCE

A fundamental problem in space science is the motion of energetic particles such as solar energetic particles or cosmic rays through the solar system. Similar transport phenomena occur in fusion plasmas and different astrophysical scenarios such as the interstellar space (see, e.g., Schlickeiser 2002; Spatschek 2008; Shalchi 2009 for reviews). Various analytical theories for describing the transport were developed in the past. The first and also simplest approach is provided by quasilinear theory (see Jokipii 1966). In his pioneering work, Palmer (1982) compared quasilinear diffusion coefficients with observational data. The data he considered can be represented by a band or box known as a Palmer consensus range (see, e.g., Figure 2 of the current paper). Palmer (1982) concluded that there is no agreement, and thus our understanding of the particle motion was incomplete.

Bieber et al. (1994) extended the standard quasilinear approach by modifying the used turbulence model. They suggested using a so-called two-component model rather than the pure slab model employed before. Furthermore, they modified the used spectrum by incorporating dissipation effects. Most importantly, however, they replaced the magnetostatic model with a dynamical turbulence model. The quasilinear calculations performed by Bieber et al. (1994) agreed well with the aforementioned Palmer consensus range.

Before Bieber et al. (1994) published their pioneering paper, it was already clear that quasilinear theory itself can be problematic. Jones et al. (1973), Völk (1973), as well as Goldstein (1976) pointed out that quasilinear theory cannot describe pitch-angle scattering at large pitch-angles correctly. Later it was shown that there are more issues with this theory. Shalchi et al. (2004) have argued that the quasilinear parallel mean free path does not agree with test-particle simulations as soon as the slab model is replaced by a two-component model. Shalchi (2015) has shown that quasilinear theory only works for perpendicular diffusion if the Kubo number is small and if pitch-angle scattering is suppressed. The role of the Kubo number in classifying the transport regimes in space plasmas was also considered in the work of Zimbardo et al. (2004, 2012).

Since quasilinear theory is questionable, one has to revisit the problem of particle transport in the solar system. In Tautz & Shalchi (2013) detailed simulations have been performed by employing a turbulence model that is very similar to the one used before by Bieber et al. (1994). A major difference, however, is that Tautz & Shalchi (2013) employed a plasma wave propagation model rather than the dynamical turbulence approach used by Bieber et al. (1994).

The purpose of the this article is to perform for the first time test-particle simulations in full dynamical turbulence. We use exactly the model suggested by Bieber et al. (1994), meaning that we use the same analytical form for the spectrum, turbulence geometry, and dynamical correlation function. Of course our numerical approach is not based on the questionable quasilinear approach used before. We compare our numerically obtained parallel mean free paths ${\lambda }_{\parallel }$ with the Palmer consensus range representing the different solar wind observations.

From our test-particle simulations we can also obtain the perpendicular mean free path ${\lambda }_{\perp }$ and the ratio of the two mean free paths ${\lambda }_{\perp }/{\lambda }_{\parallel }$. Our findings are compared with the corresponding Palmer consensus range, Jovian electrons (see Chenette et al. 1977), and Ulysses measurements of Galactic protons (see Burger et al. 2000). This type of comparison was performed before in the context of analytical theories for perpendicular diffusion (see, e.g., Bieber et al. 2004 and Shalchi et al. 2006) where good agreement was found. In the current paper, however, we compare computer simulations directly with the aforementioned solar wind observations.

The remainder of the paper is organized as follows. In Section 2 we explain the physics of turbulence by focusing on dynamical effects. The methodology that is used to perform particle transport simulations in dynamical turbulence is explained in Section 3. In Sections 4 and 5 we show our findings for the two considered dynamical turbulence models and we compare them with different solar wind observations. In Section 6 we conclude and summarize.

Magnetic turbulence, in general, is described by the two-point-two-time correlation tensor whose components are given by

$$\begin{eqnarray}&&{R}_{{mn}}({\boldsymbol{x}},{{\boldsymbol{x}}}_{0},t,{t}_{0})=\langle \delta {B}_{m}({\boldsymbol{x}},t)\delta {B}_{n}({{\boldsymbol{x}}}_{0},{t}_{0})\rangle ,\end{eqnarray} \tag{ 1 }$$

where we have used the ensemble average operator $\langle ...\rangle$. The vectors ${\boldsymbol{x}}$ and ${{\boldsymbol{x}}}_{0}$ denote arbitrary positions in space. The times t and t₀ correspond to current and initial times, respectively. It is usually convenient to describe magnetic turbulence in a Fourier or wavenumber space. Thus we can use the tensor components

$$\begin{eqnarray}&&{P}_{{mn}}({\boldsymbol{k}},t)=\langle \delta {B}_{m}({\boldsymbol{k}},t)\delta {B}_{n}^{*}({\boldsymbol{k}},{t}_{0})\rangle .\end{eqnarray} \tag{ 2 }$$

The latter components are linked to those of Equation (1) via a Fourier transformation. A standard assumption in turbulence and transport theory is that all tensor components obey the same temporal behavior and thus,

$$\begin{eqnarray}&&{P}_{{mn}}({\boldsymbol{k}},t)={P}_{{mn}}({\boldsymbol{k}}){\rm{\Gamma }}({\boldsymbol{k}},t),\end{eqnarray} \tag{ 3 }$$

where we have used the static tensor components ${P}_{{mn}}({\boldsymbol{k}})$ and the so-called dynamical correlation function$\;{\rm{\Gamma }}({\boldsymbol{k}},t)$. Magnetic turbulence is an essentially nonlinear phenomenon, so a simple description in terms of linear waves is not appropriate. In other words, the usual dispersion relation ω = ω (k) cannot be used, thus a number of empirical models are employed in test-particle simulations and analytical theories (see Bieber et al. 1994 and Shalchi et al. 2006). Some of them are summarized in Table 1. In analytical treatments of the transport, we directly use such models. As described in Section 3, this is not the case in test-particle simulations where a Fourier transformation has to be employed for the dynamical correlation function. We define

$$\begin{eqnarray}&&\chi ({\boldsymbol{k}},\omega ):= \displaystyle \frac{1}{\pi }{\mathfrak{R}}{\displaystyle \int }_{0}^{\infty }{dt}\ {\rm{\Gamma }}({\boldsymbol{k}},t){e}^{-i\omega t}.\end{eqnarray} \tag{ 4 }$$

We can solve this time integral easily for the different models listed in Table 1. The corresponding analytical forms for the function $\chi ({\boldsymbol{k}},\omega )$ are also listed in the aforementioned table. Using $\chi ({\boldsymbol{k}},\omega )$ instead of ${\rm{\Gamma }}({\boldsymbol{k}},t)$ means that we describe the turbulence in a four-dimensional Fourier space with the coordinates ${\boldsymbol{k}}$ and ω.

Table 1.  Different Models for Dynamical Magnetic Turbulence

Dynamical Turbulence Model	Correlation Function ${\rm{\Gamma }}({\boldsymbol{k}},t)$	Fourier Representation $\chi ({\boldsymbol{k}},\omega )$
Magnetostatic Turbulence	1	δ (ω)
Undamped Plasma Waves	${e}^{i{\omega }_{{\rm{p}}}t}$	$\delta ({\omega }_{{\rm{p}}}-\omega )$
Damping Model of Dynamical Turbulence	${e}^{-\gamma t}$	$\displaystyle \frac{1}{\pi }\displaystyle \frac{\gamma }{{\gamma }^{2}+{\omega }^{2}}$
Random Sweeping Model	${e}^{-{\gamma }^{2}{t}^{2}}$	$\displaystyle \frac{1}{{(2\pi )}^{3/2}\gamma }{e}^{-{\omega }^{2}/{(2\gamma )}^{2}}$
Nonlinear Anisotropic Dynamical Turbulence	${e}^{i{\omega }_{{\rm{p}}}t-\gamma t}$	$\displaystyle \frac{1}{\pi }\displaystyle \frac{\gamma }{{\gamma }^{2}+{(\omega -{\omega }_{{\rm{p}}})}^{2}}$

Note. ${\rm{\Gamma }}({\boldsymbol{k}},t)$ denotes the dynamical correlation function and $\chi ({\boldsymbol{k}},\omega )$ is its Fourier transform. The parameter ω_p corresponds to the plasma wave frequency and γ is a damping parameter. The latter parameter is explained in the main part of the text.

Download table as:  ASCIITypeset image

For the damping model of dynamical turbulence, which is one of the two models used by Bieber et al. (1994), the dynamical correlation function is an exponential (see Table 1 of the current paper). The inverse correlation time γ entering the exponential is, in this case, given by $\gamma =\alpha {v}_{{\rm{A}}}| {\boldsymbol{k}}|$, where we have used the wave vector ${\boldsymbol{k}}$ and the Alfvén speed v_A. The parameter α is a numerical factor, which, according to Bieber et al. (1994), is a number between zero (static turbulence) and one (strong dynamical turbulence). The same γ is used in the random sweeping model where the dynamical correlation function is a Gaussian (see again Table 1). In the general case, however, the parameter γ can be more complicated (see, e.g., Shalchi et al. 2006).

Of course, one has to specify the static tensor components ${P}_{{mn}}({\boldsymbol{k}})$ as well. In the current paper we employ a so-called slab/2D model, which is also known as two-component model due to the fact that there are two different wave modes, namely slab and two-dimensional (2D) modes. This type of turbulence model is frequently used to approximate solar wind turbulence and is supported by observations (see, e.g., Matthaeus et al. 1990; Osman & Horbury 2009a, 2009b; Turner et al. 2012), turbulence simulations (see, e.g., Oughton et al. 1994; Matthaeus et al. 1996; Shaikh & Zank 2007), and analytical treatments of turbulence (see, e.g., Zank & Matthaeus 1993). More details concerning the model used can be found in those papers, or, in the context of test-particle simulations, in Hussein et al. (2015).

In the following we briefly describe the most important features of the static correlation tensor we used. In slab/2D turbulence, the components of the static correlation tensor are written as ${P}_{{mn}}={P}_{{mn}}^{\mathrm{slab}}+{P}_{{mn}}^{2{\rm{D}}}$, where we have used the components of the slab tensor defined as

$$\begin{eqnarray}&&{P}_{{lm}}^{\mathrm{slab}}({\boldsymbol{k}})={g}^{\mathrm{slab}}({k}_{\parallel })\displaystyle \frac{\delta ({k}_{\perp })}{{k}_{\perp }}{\delta }_{{lm}},\end{eqnarray} \tag{ 5 }$$

and the components of the 2D tensor defined as

$$\begin{eqnarray}&&{P}_{{lm}}^{2{\rm{D}}}({\boldsymbol{k}})={g}^{2{\rm{D}}}({k}_{\perp })\displaystyle \frac{\delta ({k}_{\parallel })}{{k}_{\perp }}\left({\delta }_{{lm}}-\displaystyle \frac{{k}_{l}{k}_{m}}{{k}_{\perp }^{2}}\right).\end{eqnarray} \tag{ 6 }$$

In Equations (5) and (6) we have used the two spectra ${g}^{\mathrm{slab}}({k}_{\parallel })$ and ${g}^{2{\rm{D}}}({k}_{\perp })$, respectively. Furthermore, l, m = x, y and P_lm = 0 if one of two indexes are equal to z. Obviously the 2D tensor (6) is designed so that δ B_z = 0 which is in agreement with the model proposed originally in Bieber et al. (1994, 1996). The assumption δB_z = 0 is motivated by the fact that in the solar wind the power in parallel fluctuations is usually small in the inertial range (see Belcher & Davis 1971).

For the slab spectrum we use the model proposed by Bieber et al. (1994)

$$\begin{eqnarray}{g}^{\mathrm{slab}}({k}_{\parallel }) & = & \displaystyle \frac{C(s)}{2\pi }{l}_{\mathrm{slab}}\delta {B}_{\mathrm{slab}}^{2}\\ & & \times \;\left\{\begin{array}{ccc}{(1+{k}_{\parallel }^{2}{l}_{\mathrm{slab}}^{2})}^{-s/2} & {\rm{for}} & {k}_{\parallel }\leqslant {k}_{{\rm{d}}}\\ {(1+{k}_{{\rm{d}}}^{2}{l}_{\mathrm{slab}}^{2})}^{-s/2}{({k}_{{\rm{d}}}/{k}_{\parallel })}^{p} & {\rm{for}} & {k}_{\parallel }\geqslant {k}_{{\rm{d}}}\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where we have used the slab bendover scale l_slab, the dissipation wavenumber k_d, the inertial range spectral index s, and the dissipation-range spectral index p. Furthermore, we have employed the normalization function

$$\begin{eqnarray}&&C(s)=\displaystyle \frac{{\rm{\Gamma }}\left(\displaystyle \frac{s}{2}\right)}{2\sqrt{\pi }{\rm{\Gamma }}\left(\displaystyle \frac{s-1}{2}\right)},\end{eqnarray} \tag{ 8 }$$

with the Gamma function Γ (z).

For the 2D spectrum we employ an extension of the model proposed by Bieber et al. (1994). By combining the spectrum used in the latter paper with the ideas discussed in Matthaeus et al. (2007) and Shalchi & Weinhorst (2009), we suggest using

$$\begin{eqnarray}{g}^{2{\rm{D}}}({k}_{\perp }) & = & \displaystyle \frac{2D(s,q)}{\pi }{l}_{2{\rm{D}}}\delta {B}_{2{\rm{D}}}^{2}\\ & & \times \;\left\{\begin{array}{ccc}\displaystyle \frac{{({k}_{\perp }{l}_{2{\rm{D}}})}^{q}}{{(1+{k}_{\perp }^{2}{l}_{2{\rm{D}}}^{2})}^{(s+q)/2}} & {\rm{for}} & {k}_{\perp }\leqslant {k}_{{\rm{d}}}\\ \displaystyle \frac{{({k}_{{\rm{d}}}{l}_{2{\rm{D}}})}^{q}}{{(1+{k}_{{\rm{d}}}^{2}{l}_{2{\rm{D}}}^{2})}^{(s+q)/2}}{\left(\displaystyle \frac{{k}_{{\rm{d}}}}{{k}_{\perp }}\right)}^{p} & {\rm{for}} & {k}_{\perp }\geqslant {k}_{{\rm{d}}}.\end{array}\right.\end{eqnarray} \tag{ 9 }$$

The only parameter that is new compared to the slab spectrum is the energy range spectral index q controlling the spectral shape at large turbulence scales. Furthermore, we have used the extended normalization function

$$\begin{eqnarray}&&D(s,q)=\displaystyle \frac{{\rm{\Gamma }}\left(\displaystyle \frac{s+q}{2}\right)}{2{\rm{\Gamma }}\left(\displaystyle \frac{s-1}{2}\right){\rm{\Gamma }}\left(\displaystyle \frac{q+1}{2}\right)}.\end{eqnarray} \tag{ 10 }$$

In Table 2 we list the values we have used for the different turbulence and particle parameters. Now our used turbulence model is fully defined and in the next section we describe how this model is incorporated into test-particle codes.

Table 2.  The Parameter Values Used for Our Test-particle Simulations

Parameter	Symbol	Value
2D energy range spectral index	q	2
Inertial range spectral index	s	5/3
Dissipation range spectral index	p	3
Alfvén speed	vA	33.5 km s−1
Slab bendover scale	lslab	0.030 AU
2D bendover scale	${l}_{2{\rm{D}}}$	0.030 AU
Dissipation wavenumber	kd	3 × 105 1/AU
Mean field	B0	4.12 nT
Slab fraction	$\delta {B}_{\mathrm{slab}}^{2}$	$0.2\;\delta {B}^{2}$
2D fraction	$\delta {B}_{2{\rm{D}}}^{2}$	$0.8\;\delta {B}^{2}$

Note. The values should be appropriate in the interplanetary space at 1 AU heliocentric distance.

Download table as:  ASCIITypeset image

In order to simulate the propagation of energetic particles through a plasma, we have to perform different steps. The first one is the creation of the turbulent magnetic fields based on a specific turbulence model (e.g., the slab model or two-component turbulence). Thereafter we need to solve numerically the Newton–Lorentz equation for an ensemble of particles to obtain their orbits. The last step is the computation of diffusion parameters by averaging the obtained particle orbits in an appropriate way.

In order to calculate the turbulent magnetic field at the position of the charged particle ${\boldsymbol{x}}$, one can use the Fourier representation

$$\begin{eqnarray}&&\delta {\boldsymbol{B}}({\boldsymbol{x}},t)=\int {d}^{3}k\ \delta {\boldsymbol{B}}({\boldsymbol{k}},t){e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}.\end{eqnarray} \tag{ 11 }$$

In numerical treatments of the transport, the wavenumber integrals in Equation (11) have to be replaced by sums. In the current paper we only simulate turbulence models with reduced dimensionality and a superposition of such models. Therefore, the wavenumber integral can be replaced by a single sum. Since we are dealing with dynamical turbulence, the dynamical correlation function is also replaced by the corresponding Fourier representation as described in Section 2. Therefore, compared to static simulations a second sum occurs in the formula for the magnetic field creation. By following previous work (see, e.g., Michałek & Ostrowski 1996; Giacalone & Jokipii 1999; Tautz 2010; Hussein et al. 2015), in combination with our dynamical turbulence approach, we can write the turbulent magnetic field as

$$\begin{eqnarray}&&\delta B({\boldsymbol{x}},t)=\sqrt{2}\;\delta B\displaystyle \sum _{m=1}^{M}\;\displaystyle \sum _{n=1}^{N}A({k}_{m},{\omega }_{n}){{\boldsymbol{\xi }}}_{m}{e}^{i[{{\boldsymbol{k}}}_{m}\cdot {\boldsymbol{x}}+{\omega }_{n}t+{\beta }_{{mn}}]}\end{eqnarray} \tag{ 12 }$$

with the random phase β_mn. For the slab modes and the 2D modes we have used the same polarization vector ${\boldsymbol{\xi }}$, namely

$$\begin{eqnarray}&&{{\boldsymbol{\xi }}}_{m}=(-\mathrm{sin}{\phi }_{m},\mathrm{cos}{\phi }_{m},0),\end{eqnarray} \tag{ 13 }$$

where ϕ_m is a random angle.

All quantities used above are normalized with respect to the slab bendover scale l_slab. This means, for instance, that k_m in the code represents the physical quantity ${k}_{\parallel }{l}_{\mathrm{slab}}$ or ${k}_{\perp }{l}_{\mathrm{slab}}$ and z stands for $z/{l}_{\mathrm{slab}}$. The frequency ω is normalized with respect to the unperturbed gyrofrequency of the particle Ω, hence ${\omega }_{n}=\omega /{\rm{\Omega }}$.

In Equation (12) we have also used the wavevector ${{\boldsymbol{k}}}_{m}={k}_{m}{\hat{k}}_{m}$ with the random wave unit vector

$$\begin{eqnarray}&&{\hat{k}}_{m}=\left(\begin{array}{c}\sqrt{1-{\eta }_{m}^{2}}\mathrm{cos}{\phi }_{m}\\ \sqrt{1-{\eta }_{m}^{2}}\mathrm{sin}{\phi }_{m}\\ {\eta }_{m}\end{array}\right).\end{eqnarray} \tag{ 14 }$$

As above, the angle ϕ_m is a random number. What the value of η_m is depends on the simulated turbulence model. For slab modes we have η_m = 1 and for 2D modes η_m = 0. In Equation (12) we have also used the amplitude function

$$\begin{eqnarray}&&{A}^{2}({\omega }_{n},{k}_{m})=\displaystyle \frac{G({k}_{m},{\omega }_{n}){\rm{\Delta }}{k}_{m}{\rm{\Delta }}{\omega }_{n}}{{\displaystyle \sum }_{\mu =1}^{M}{\displaystyle \sum }_{\nu =1}^{N}G({k}_{\mu },{\omega }_{\nu }){\rm{\Delta }}{k}_{\mu }{\rm{\Delta }}{\omega }_{\nu }}.\end{eqnarray} \tag{ 15 }$$

The function G(k_μ, ω_ν) represents the spacetime spectrum for which we use

$$\begin{eqnarray}&&G({k}_{m},{\omega }_{n})=G({k}_{m})\chi ({k}_{m},{\omega }_{n}).\end{eqnarray} \tag{ 16 }$$

The function χ (k_m, ω_n) can be obtained from Table 1. In the following we employ the damping model of dynamical turbulence as well as the random sweeping model (see Section 2 for more details). The function G(k_m) used here is the usual spectrum, as used in simulations of magnetostatic turbulence (see, e.g., Hussein et al. 2015). In the current paper we use Equation (7) for the slab modes and Equation (9) for the 2D modes.

In the model used for $\chi ({k}_{m},{\omega }_{n})$, one finds the Alfvén speed v_A. The latter parameter can be normalized with respect to the particle speed v so that

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{{\rm{A}}}}{v}=\displaystyle \frac{{v}_{{\rm{A}}}}{c}\displaystyle \frac{\sqrt{{R}_{0}^{2}+{R}^{2}}}{R},\end{eqnarray} \tag{ 17 }$$

with the speed of light c and

$$\begin{eqnarray}{R}_{0}=\displaystyle \frac{1}{{l}_{\mathrm{slab}}{B}_{0}}\left\{\begin{array}{ccc}0.511{\rm{MV}}\quad & {\rm{for}}\quad & {\rm{electrons}}\\ 938{\rm{MV}}\quad & {\rm{for}}\quad & {\rm{protons}}.\end{array}\right.\end{eqnarray} \tag{ 18 }$$

Here we have used the dimensionless rigidity defined via R = R_L/l_slab, where R_L = v/Ω is the unperturbed Larmor radius and l_slab is the slab bendover scale used above. For l_slab = 0.03 AU and B₀ = 4.12 nT this gives R₀ = 9.2 × 10⁻⁵ for electrons and R₀ = 0.169 for protons. All other parameter values used for our simulations are summarized in Table 2.

In order to perform the simulations with high accuracy, one has to take into account several conditions. In our simulations of particle transport in dynamical turbulence we have to deal with the same issues as in the simulations of static turbulence (see again Hussein et al. 2015 for more details). However, there are a few additional concerns when it comes to dynamical turbulence. For instance, we have to deal with the minimal value of the parameter ${\omega }_{n}$ (see below). Furthermore, we need a certain number of grid points in space and time. For most of our runs we have set N = M = 256 in Equation (12). For lower rigidities we had to set N = M = 64 to avoid computation times that were too long. In all our runs we have computed running diffusion coefficients for times up to at least Ωt = 10⁴ (here Ω denotes again the gyrofrequency) to ensure that we are in the stable regime.

While we have performed the simulations, we have realized that the value of the minimum frequency, ω_min, has a strong influence on both parallel and perpendicular mean free paths. This influence was noticed for both protons and electrons but was much stronger for electrons. To address this issue, we performed a series of test simulations for electrons with $R\in \ [{10}^{-2},5\times {10}^{-2},{10}^{-3}]$, with values of ω_min ranging from 10⁻²⁰ to 10⁻². We define the parameter ζ as the ratio of the outcome ${\lambda }_{\parallel }$ from simulations with a specific value of ω_min to the final value from the same simulations with ω_min = 10⁻²⁰. The results are shown in Figure 1. Small values of ω_min are required in order to obtain the correct results from the simulations. This issue is a consequence of the form of the resonance function used here (see Table 1).

Zoom In Zoom Out Reset image size

Following the ideas presented in Tautz (2010) we also compute the errors of the different mean free paths. The author noted that using the standard deviation as a mean of estimating the error is inappropriate, as the mean square displacement calculated in the Monte Carlo code is the variance of the distribution function for the diffusion equation itself. In addition, test particles interacting with turbulent magnetic fields scatter in a random manner, leading to a huge variance in their square deviation. Hence one has to come up with a method that takes into account the averaging processes used over the number of turbulence manifestations, N_T, for each of which a fixed number of test particles were simulated in space and time, resulting in a diffusion coefficient. The mean error is then defined to be the deviation of the different mean free paths for a specific turbulence, λ_n, from the final averaged mean free path, λ_f. Mathematically that reads as

$$\begin{eqnarray}{\sigma }_{\lambda }^{2} & = & \displaystyle \frac{1}{{N}_{T}-1}\\ & & \times \left\{\displaystyle \sum _{{N}_{T}}{({\lambda }_{n}-{\lambda }_{f})}^{2}-\displaystyle \frac{1}{{N}_{T}}{\left[\displaystyle \sum _{{N}_{T}}({\lambda }_{n}-{\lambda }_{f})\right]}^{2}\right\}.\end{eqnarray} \tag{ 19 }$$

Using Equation (19), both the error in the parallel and perpendicular mean free paths were calculated as ${\rm{\Delta }}{\lambda }_{\parallel }$ and ${\rm{\Delta }}{\lambda }_{\perp }$ respectively. To calculate the error in the ratio of the two mean free paths, ${\lambda }_{\perp }/{\lambda }_{\parallel }$, we use the rule of error combination

$$\begin{eqnarray}&&{\rm{\Delta }}\left(\displaystyle \frac{{\lambda }_{\perp }}{{\lambda }_{\parallel }}\right)=\left(\displaystyle \frac{{\rm{\Delta }}{\lambda }_{\perp }}{{\lambda }_{\perp }}+\displaystyle \frac{{\rm{\Delta }}{\lambda }_{\parallel }}{{\lambda }_{\parallel }}\right)\displaystyle \frac{{\lambda }_{\perp }}{{\lambda }_{\parallel }}.\end{eqnarray} \tag{ 20 }$$

In all plots shown in Sections 4 and 5 we have shown the error bars based on the discussion presented here.

We started our investigations by employing the damping model of dynamical turbulence. We have performed simulations for three sets of parameters. First we have used pure slab turbulence and set δB/B₀ = 1.0. In the second run we replaced the slab model with a 20%/80% slab/2D composite model but we kept the magnetic field ratio the same. The assumed contribution of the slab and 2D models is in agreement with the work of Bieber et al. (1994, 1996). In the last set of our simulations we still used the two-component model but the magnetic field ratio was reduced to δB/B₀ = 0.5, as suggested by Ruffolo et al. (2012).

4.1. Slab Turbulence with δB/B₀ = 1.0

In the first set of simulations we have employed the slab model of magnetic turbulence in combination with the damping model of dynamical turbulence. We have computed the parallel and perpendicular diffusion coefficients versus time for different rigidities.

For parallel transport we clearly obtained a diffusive behavior for all considered cases. The obtained parallel mean free path versus rigidity is visualized in Figure 2. The parallel mean free paths we have obtained numerically agree with the quasilinear results derived in Bieber et al. (1994). However, as in this previous work, the parallel mean free paths are clearly below the Palmer consensus range. Therefore, we agree with the statement of Bieber et al. (1994) that a pure slab setup cannot provide an appropriate approximation for solar wind turbulence. Therefore, the slab model has to be replaced with a more realistic model such as the two-component model. This is done in Sections 4.2 and 4.3.

Zoom In Zoom Out Reset image size

For perpendicular transport we observed subdiffusive transport after the initial ballistic regime. The Unified Nonlinear Transport (UNLT) theory developed in Shalchi (2010) can be used to compute the perpendicular diffusion coefficient analytically for dynamical turbulence (see Shalchi 2014). The latter theory predicts that in the limit $t\to \infty$ perpendicular diffusion is recovered but the corresponding perpendicular mean free path is very small. The findings of the current paper agree with such analytical predictions in the sense that all our numerical results are larger than the analytical values.

4.2. Slab/2D Turbulence with δB/B₀ = 1.0

Since theoretical results obtained by using either quasilinear calculations (see Bieber et al. 1994) or numerical simulations (see previous paragraph) do not agree with observations as long as the slab model is used, we now employ the 20%/80% composite model proposed by Bieber et al. (1994, 1996).

The parallel mean free path versus magnetic rigidity is shown in Figure 3. Although the parallel mean free path is now increased compared to the pure slab result, we are still far below the Palmer consensus range. Whereas quasilinear estimates provide a parallel mean free path that is a factor of five larger compared to the slab result, this effect is much smaller if simulations are used. The physical explanation for that different behavior was presented in Shalchi et al. (2006). For the two-component model used here, resonance broadening due to perpendicular diffusion becomes important. This type of nonlinear effect makes the parallel mean free path smaller compared to the quasilinear result. Thus, we conclude that the idea of replacing the slab model with a two-component model helps to increase the parallel diffusion coefficient, but this effect is not strong enough to achieve agreement with solar wind observations.

Zoom In Zoom Out Reset image size

For the perpendicular diffusion coefficient we found diffusive behavior for all considered cases. The perpendicular mean free path versus magnetic rigidity is visualized in Figure 4 and the ratio of the two diffusion parameters ${\lambda }_{\perp }/{\lambda }_{\parallel }$ versus rigidity is shown in Figure 5. The latter ratio is almost constant for small rigidities and then decreases rapidly if the rigidity is increased. According to Figure 4 the perpendicular mean free path eventually becomes constant in the high rigidity limit. This behavior of the perpendicular diffusion coefficient is exactly what is predicted by the UNLT theory, as shown in detail in Shalchi (2015).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Slab/2D Turbulence with δB/B₀ = 0.5

As argued in the previous section, the parallel mean free paths obtained from our simulations are still too small compared to the observations. Those results, however, were obtained for a magnetic field ratio of δB/B₀ = 1.0. Ruffolo et al. (2012) suggested that the latter ratio is δB/B₀ = 0.5. Furthermore, Tautz & Shalchi (2013) successfully used this magnetic field ratio in their simulations of propagating plasma waves. Thus, we employ the same value here and compute parallel and perpendicular diffusion parameters again.

The parallel mean free path versus magnetic rigidity is shown in Figure 6. We can now clearly see that the electron parallel mean free path goes through the Palmer consensus range. Thus we conclude that we can indeed reproduce the solar wind observations by using our test-particle code. However, besides using the two-component turbulence model, we also have to reduce the magnetic field ratio to δB/B₀ = 0.5.

Zoom In Zoom Out Reset image size

The perpendicular mean free path and the ratio of the two diffusion coefficients are shown in Figures 7 and 8, respectively. Qualitatively, the perpendicular diffusion coefficient obeys the same behavior as above, only the magnitudes are different. We compare our numerical findings with observations of Jovian electrons (Chenette et al. 1977), Ulysses measurements of Galactic protons (Burger et al. 2000), and the Palmer (1982) value. Good agreement between the simulations and observations is found.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the current paragraph we replace the damping model of dynamical turbulence used above with a random sweeping model where the dynamical correlation function is a Gaussian function (see Table 1). Motivated by the results obtained in Section 4, we only run the simulations for two-component turbulence and a magnetic field ratio of δB/B₀ = 0.5.

The parallel mean free path versus magnetic rigidity is shown in Figure 9. Again, the electron parallel mean free path goes through the Palmer consensus range and we can reproduce the solar wind observations by using our test-particle code. By comparing Figures 6 and 9, it can be seen that the results are very similar and the differences between the damping model of dynamical turbulence and the random sweeping model are minor.

Zoom In Zoom Out Reset image size

The perpendicular mean free path and the ratio of the two diffusion coefficients are shown in Figures 10 and 11, respectively. By comparing those results with our findings obtained for the damping model of dynamical turbulence (see Figures 7 and 8), we conclude that there is not a strong influence from the chosen dynamical turbulence model.

Zoom In Zoom Out Reset image size

We have revisited the problem of particle transport in interplanetary space. We have combined test-particle simulations with the model for solar wind turbulence proposed by Bieber et al. (1994). The latter model contains the following key features:

• 1.  
  The spectrum has energy, inertial, and dissipation range and allows us to modify the spectral shape by changing the corresponding indexes q, s, and p;
• 2.  
  The used magnetic correlation tensor contains two contributions or wave modes, namely slab and 2D modes; and
• 3.  
  The turbulence is assumed to be dynamical and we have employed two different models for the dynamical correlations function, namely the damping model of dynamical turbulence and the random sweeping model.

Whereas the used turbulence model is exactly the one described and used in Bieber et al. (1994), our approach to compute particle diffusion parameters is very different. Bieber et al. performed quasilinear calculations, whereas our work is entirely based on computer simulations. Similar simulations have been used before by considering static turbulence (see, e.g., Giacalone & Jokipii 1999; Qin et al. 2002a, 2002b) or undamped propagating plasma waves (see, e.g., Michałek & Ostrowski 1996; Tautz & Shalchi 2013). For the first time we have performed test-particle simulations for dynamical turbulence by employing a more-dimensional Fourier transformation approach, as described in Section 3.

In Figures 2–8 we have shown our numerical findings that were obtained by using the damping model of dynamical turbulence, and compared them with the Palmer (1982) consensus range representing the different solar wind observations. As already concluded by Bieber et al. (1994), we are not able to reproduce the consensus range by simply employing a slab model. The use of a two-component turbulence model is required but even then we are not able to find agreement with the observations. This finding is a major difference between previous work and our simulations. The idea of Bieber et al. (1994) was that by replacing the slab model with a two-component model, we could increase the parallel mean free path drastically. This idea, however, works only in the framework of quasilinear theory. The theory, however, is incorrect due to nonlinear effects, as already pointed out in Shalchi et al. (2004). The simulations are quasi-exact apart from numerical inaccuracies. As shown in the current paper, just replacing the slab model with a slab/2D model alone does not lead to agreement with the observations. Despite the aforementioned issue, we were still able to reproduce the Palmer consensus range by assuming that the turbulent magnetic field is a bit weaker compared to what was assumed before. If we assume that δB/B₀ = 0.5, as suggested by Ruffolo et al. (2012), instead of the standard assumption δB/B₀ = 1.0, agreement between simulations and observations can indeed be found.

We have repeated our simulations by using the random sweeping model proposed by Bieber et al. (1994). We have again employed  the slab/2D model and assumed δB/B₀ = 0.5 for the magnetic field ratio. The corresponding results for the parallel mean free path, perpendicular mean free path, and the ratio of the two latter parameters are visualized in Figures 9–11. Our findings are similar to those obtained for the damping model of dynamical turbulence. Again the parallel mean free paths agree very well with the Palmer consensus range.

Zoom In Zoom Out Reset image size

For the perpendicular diffusion coefficient we found very similar results for the damping model of dynamical turbulence and the random sweeping model. The same applies for the ratio ${\lambda }_{\perp }/{\lambda }_{\parallel }$. We have compared our numerical results with different solar wind observations, namely observations of Jovian electrons (Chenette et al. 1977), Ulysses measurements of Galactic protons (Burger et al. 2000), and the Palmer (1982) value. For two-component turbulence and a magnetic field ratio of δB/B₀ = 0.5 good agreement is found for the damping model of dynamical turbulence as well as the random sweeping model.

Reames (1999) argued that the parallel mean free path can be of the order of 1–2 AU for low rigidity particles at a 1 AU heliocentric distance. The largest value for the parallel mean free path in this regime was obtained for the random sweeping model (see Figure 9) where we found ${\lambda }_{\parallel }\approx 0.25$AU. Obviously our numerical results are about a factor of four below the value of Reames (1999). One could try to further reduce the magnetic field ratio δB/B₀ or to use a more realistic dynamical turbulence model. This will be subject of future work.

We would also like to point out that there is some similarity between our findings and the results found previously by Tautz & Shalchi (2013). The latter authors, however, used undamped propagating plasma waves, which corresponds to a different turbulence setup. We conclude that propagating plasma waves as well as dynamical turbulence allow for a reproduction of the Palmer consensus range.

We also conclude that the Bieber et al. (1994) turbulence model can indeed be seen as a valid model for particle transport. However, this requires setting δB/B₀ = 0.5, a value which is a factor of two smaller compared to the value originally used. Finally, we can state that with this magnetic field ratio we can indeed reproduce solar wind observations by using our test-particle code.

A.S. acknowledges support by the Natural Sciences and Engineering Research Council (NSERC) of Canada. Most simulations shown in this article were obtained by using the national computational facility provided by WestGrid. We are also grateful to S. Safi-Harb for providing her CFI-funded computational facilities for code tests and for some of the simulation runs presented here.