COSMIC-RAY PITCH-ANGLE SCATTERING IN IMBALANCED MHD TURBULENCE SIMULATIONS

The diffusion of a beam of charged cosmic-ray particles with velocity v in a turbulent magnetic field ${\boldsymbol{b}}({\boldsymbol{x}})$ is commonly described in terms of a Fokker–Planck equation for the cosmic-ray phase-space density $f(\mu ,z)$, where $\mu =\mathrm{cos}({v}_{z}/v)$ is the cosine of the pitch-angle and z is the coordinate along the global mean magnetic field ${B}_{0}\hat{z}$ (e. g. Schlickeiser 2002; Shalchi 2009):

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+v\ \mu \ \displaystyle \frac{\partial f}{\partial z}=\displaystyle \frac{\partial }{\partial \mu }({D}_{\mu \mu }\ \displaystyle \frac{\partial f}{\partial \mu }).\end{eqnarray} \tag{ 1 }$$

The form of the quasilinear pitch-angle diffusion coefficient ${D}_{\mu \mu }={\mathrm{lim}}_{t\to \infty }\langle {[\mu (t)-\mu (0)]}^{2}\rangle /2t$ was first derived for cosmic-ray particles in magnetostatic slab turbulence by Jokipii (1966). Expanding the turbulent magnetic field ${\boldsymbol{b}}({\boldsymbol{x}})={B}_{0}\hat{z}\;+\delta {{\boldsymbol{b}}}_{\perp }(z)$ in a perpendicular perturbation $\delta {{\boldsymbol{b}}}_{\perp }(z)\perp \hat{z}$ that only depends on the coordinate along the zeroth-order mean field ${B}_{0}\hat{z}$, he found that

$$\begin{eqnarray}&&{D}_{\mu \mu }=\displaystyle \frac{\pi }{4}\ {{\rm{\Omega }}}^{2-s}(1-{\mu }^{2})\ \displaystyle \frac{{\boldsymbol{\delta }}{{\boldsymbol{B}}}^{2}}{{B}_{0}^{2}}\ {| v\mu | }^{s-1}(s-1){k}_{\mathrm{min}}^{s-1}.\end{eqnarray} \tag{ 2 }$$

Here, the power spectrum of the magnetic field perturbations is assumed as a slab spectrum ${| \delta {{\boldsymbol{b}}}_{\perp }({k}_{z})| }^{2}{{dk}}_{z}\;={\boldsymbol{\delta }}{{\boldsymbol{B}}}^{2}(s-1){k}_{\mathrm{min}}^{s-1}{k}_{z}^{-s}{{dk}}_{z}$ with a spectral index s, a lower cutoff wavenumber at ${k}_{\mathrm{min}}$, and root-mean-square amplitude $\langle \delta {{\boldsymbol{b}}}_{\perp }^{2}{\rangle }^{1/2}=| {\boldsymbol{\delta }}{\boldsymbol{B}}|$.

As the typical velocity of cosmic-ray particles $v\sim c$ exceeds the typical Alfvén velocity in the interstellar medium ${v}_{{\rm{A}}}\sim 10\;\mathrm{km}\;{{\rm{s}}}^{-1}$ by several orders of magnitude, the magnetic turbulence can be approximated as being time-independent. In the limit of very small wave frequencies ${\rm{\Omega }}\ll {{vk}}_{\mathrm{min}}$, we can also presume that the characteristic electric field strength $\delta E$ becomes negligibly small since Faraday's law implies $\delta E\sim {\rm{\Omega }}| {\boldsymbol{\delta }}{\boldsymbol{B}}| /k$. As a consequence of these assumptions, the cross helicity does not enter into Equation (2).

For charged particles with $v\sim {v}_{{\rm{A}}}$, however, the propagation speed of shear-Alfvén waves must be considered. As first shown by Schlickeiser (1989), all powers of $\gamma ={v}_{{\rm{A}}}/v$ must be included in the derivation of ${D}_{\mu \mu }$ to obtain the correct result. Accounting for the time-dependence of the magnetic turbulence and, through Faraday's law, the resulting electric fields self-consistently, one finds for the quasilinear pitch-angle diffusion coefficient in incompressible MHD turbulence:

$$\begin{eqnarray}{D}_{\mu \mu } & = & \displaystyle \frac{\pi }{2}\ {{\rm{\Omega }}}^{2}\ {B}_{0}^{-2}\ (1-{\mu }^{2})\\ & & \times \displaystyle \sum _{\beta =\pm 1}{| v\mu +\beta {v}_{{\rm{A}}}| }^{-1}\ {{\mathcal{E}}}_{B}^{\beta }({k}_{\mathrm{res}}^{(\beta )}).\end{eqnarray} \tag{ 3 }$$

Here $\beta \in \{+1,-1\}$ distinguishes shear-Alfvén waves traveling along ($\beta =-1$) or opposite ($\beta =+1$) the mean-field axis, the only types of waves that we consider in our simplified model of incompressible slab turbulence. We neglect additional contributions to the Fokker–Planck equation proportional to ${\partial }_{p}f$ and ${\partial }_{p}^{2}f$ that result in adiabatic focusing and stochastic heating, focusing only on the pitch-angle evolution of a monoenergetic particle distribution.

The total turbulence energy of the system can be separated into the energies contained in counter- and co-propagating waves traveling opposite to or along the direction of the mean magnetic field, the so-called positive and negative Elsasser energies ${{\mathcal{E}}}^{\pm }=\langle {({\boldsymbol{u}}\pm \delta {\boldsymbol{b}})}^{2}\rangle /4$. In quasilinear theory, an analogous decomposition into two Alfvén-wave populations is used for only the magnetic turbulence energy as well (Dung & Schlickeiser 1990), and we have written the respective components of the magnetic-turbulence power spectrum as ${{\mathcal{E}}}_{B}^{\beta }({k}_{z})$ above. The wave numbers resonant with each wave population are defined as ${k}_{\mathrm{res}}^{(\beta )}={\rm{\Omega }}/(v\mu +\beta {v}_{{\rm{A}}})$. In case the power spectra for both wave populations are not equal, the degree of imbalance between both propagation directions can be captured in the normalized cross helicity

$$\begin{eqnarray}&&{\sigma }_{c}=\displaystyle \frac{{{\mathcal{E}}}^{+}-{{\mathcal{E}}}^{-}}{{{\mathcal{E}}}^{+}+{{\mathcal{E}}}^{-}}.\end{eqnarray} \tag{ 4 }$$

Written in terms of ${\sigma }_{c}$, in the electrodynamic case the pitch-angle diffusion coefficient ${D}_{\mu \mu }$ becomes (Dung & Schlickeiser 1990)

$$\begin{eqnarray}{D}_{\mu \mu } & = & \displaystyle \frac{\pi }{4}(1-{\mu }^{2})\displaystyle \sum _{\beta =\pm }{{\rm{\Omega }}}^{2-{s}^{\beta }}\displaystyle \frac{{({\boldsymbol{\delta }}{\boldsymbol{B}})}^{2}}{{B}_{0}^{2}}({s}^{\beta }-1){k}_{\mathrm{min}}^{{s}^{\beta }-1}\\ & & \times {v}^{{s}^{\beta }-1}{(1+\beta \mu \gamma )}^{2}\ {| \mu +\beta \gamma | }^{{s}^{\beta }-1}\ (1+\beta {\sigma }_{c}).\end{eqnarray} \tag{ 5 }$$

To obtain this equation, we assume that the magnetic turbulence spectrum can be decomposed into contributions by co- and counter-propagating Alfvén waves in the same manner as the Elsasser energies, so that the two spectra

$$\begin{eqnarray}&&{| {\boldsymbol{\delta }}{{\boldsymbol{b}}}^{\pm }({k}_{z})| }^{2}{{dk}}_{z}={({\boldsymbol{\delta }}{{\boldsymbol{B}}}^{\pm })}^{2}\ ({s}^{\pm 1})\ {k}_{\mathrm{min}}^{{s}^{\pm }-1}\ {k}_{z}^{-{s}^{\pm }}\ {{dk}}_{z}\end{eqnarray} \tag{ 6 }$$

describe the magnetic turbulence due to the two Alfvén-wave populations, where ${({\boldsymbol{\delta }}{{\boldsymbol{B}}}^{\pm })}^{2}={\boldsymbol{\delta }}{{\boldsymbol{B}}}^{2}(1\pm {\sigma }_{c})$ is the root-mean-square of the magnetic turbulence in the two wave populations.

In order to distinguish between the pitch-angle scattering caused by the electric-field acceleration and that due to the magnetic field alone, we have performed test-particle runs in which the electric-field component of the Lorentz force was not included. Without electric heating, the pitch-angle diffusion Equation (1) becomes exact again, although ${D}_{\mu \mu }$ must be adapted to this "magnetodynamic" turbulence model. Generalizing the results that Shalchi et al. (2009) published for balanced turbulence to arbitrary cross helicity, we obtain for the pitch-angle diffusion coefficient:

$$\begin{eqnarray}{D}_{\mu \mu } & = & \displaystyle \frac{\pi }{4}\ \displaystyle \sum _{\beta =\pm 1}{{\rm{\Omega }}}^{2-{s}^{\beta }}\displaystyle \frac{{({\boldsymbol{\delta }}{\boldsymbol{B}})}^{2}}{{B}_{0}^{2}}\ ({s}^{\beta }-1)\ {k}_{\mathrm{min}}^{{s}^{\beta }-1}\ {v}^{{s}^{\beta }-1}\\ & & \times (1-{\mu }^{2})\ {| \mu +\beta \gamma | }^{{s}^{\beta }-1}\ (1+\beta {\sigma }_{c}).\end{eqnarray} \tag{ 7 }$$

As one can easily check, Equations (5) and (7) are identical to (2) if the Alfvén velocity becomes negligibly small compared to the particle velocity ($v\gg {v}_{{\rm{A}}}$ or equivalently $\gamma \ll 1$) and both wave populations exhibit the same spectral exponent (${s}^{+}={s}^{-}$). The quasilinear predictions for the pitch-angle diffusion coefficient in all three different cases in zero cross-helicity turbulence and strongly imbalanced turbulence are compared in Figure 1.

Zoom In Zoom Out Reset image size

In order to compare the validity of the quasilinear slab-turbulence results for the pitch-angle diffusion with numerical results, we use the pseudospectral MHD code Turbo (Teaca et al. 2009) to obtain representations of turbulent fields in which we propagate test-particles. Employing Alfvénic units, Turbo solves the equations of resistive incompressible MHD on a cubic grid with periodic boundaries:

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{u}}=-({\boldsymbol{u}}\cdot {\rm{\nabla }}){\boldsymbol{u}}+({\boldsymbol{b}}\cdot {\rm{\nabla }}){\boldsymbol{b}}+\nu \ {{\rm{\nabla }}}^{2}{\boldsymbol{u}}+{{\boldsymbol{f}}}^{u}-{\rm{\nabla }}\tilde{p},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{b}}=-({\boldsymbol{u}}\cdot {\rm{\nabla }}){\boldsymbol{b}}+({\boldsymbol{b}}\cdot {\rm{\nabla }}){\boldsymbol{u}}+\eta \ {{\rm{\nabla }}}^{2}{\boldsymbol{b}}+{{\boldsymbol{f}}}^{b}.\end{eqnarray} \tag{ 9 }$$

Here $\tilde{p}=p/\rho$ is a re-normalized pressure that ensures that the velocity field remains divergence-free, ν and η are kinetic viscosity and magnetic diffusivity, respectively, while ${{\boldsymbol{f}}}_{u}$ and ${{\boldsymbol{f}}}_{b}$ are forcing fields. The magnetic field b includes a constant guiding field ${{\boldsymbol{B}}}_{{\boldsymbol{0}}}$ applied along the z direction.

The cross-helical forcing scheme we use decomposes the Fourier-transformed fields ${\boldsymbol{u}}({\boldsymbol{k}})$ and ${\boldsymbol{b}}({\boldsymbol{k}})$ into eigenmodes of the curl operator. These helical eigenmodes for wave vectors in the ellipsoidal shell $2.5\leqslant {k}_{\perp }^{2}+{10}^{2}{k}_{z}^{2}\leqslant 3.5$ are then forced separately in such a way that the injected amounts of energy $({\mathcal{E}}),$ kinetic helicity (${{\mathcal{H}}}_{\mathrm{kin}}$), magnetic helicity (${{\mathcal{H}}}_{\mathrm{mag}}$), and cross-helicity (${\mathcal{K}}=\langle {\boldsymbol{u}}\cdot {\boldsymbol{\delta }}{\boldsymbol{b}}\rangle ={\sigma }_{c}{\mathcal{E}}$) are constant. Accordingly, we define the perpendicular and parallel correlation lengths as ${\lambda }_{\perp }=2\pi /3$ and ${\lambda }_{\parallel }=10{\lambda }_{\perp }$. The k-space ellipsoid is shortened by a factor of 10 along the k_z axis to allow for anisotropic forcing of the turbulent fields so that the turbulence can attain stable steady-states in the presence of a magnetic guide field. In particular, turbulent energy can be supplied to the system at the rate $\epsilon =\partial {\mathcal{E}}/\partial t$ without affecting any of the helicities, or cross-helicity can be injected at a desired rate $\sigma =(\partial {\mathcal{K}}/\partial t)/\epsilon$ to drive the system toward an imbalanced steady-state with ${\mathcal{K}}\ne 0$.

The fields used for the test-particle runs were obtained on 512³ grids with physical side lengths of $2\pi \;\times \;2\pi \;\times \;20\pi$ and periodic boundary conditions. This anisotropically shaped simulation domain allows the turbulence to reach a steady-state in the presence of a strong magnetic mean-field along the z direction, ${B}_{0}=10\;| {\boldsymbol{\delta }}{\boldsymbol{B}}|$, with $| {\boldsymbol{\delta }}{\boldsymbol{B}}| =\langle {({\boldsymbol{b}}-{{\boldsymbol{B}}}_{{\boldsymbol{0}}})}^{2}{\rangle }^{1/2}$.

Viscosity and magnetic diffusivity are chosen to ensure that Kolmogorov's dissipation length ${{\ell }}_{{\rm{K}}}={({\nu }^{3}/\epsilon )}^{1/4}$ is resolved by the numerical grid in all directions. Since we use isotropic viscosity and diffusivity coefficients, energy dissipation remains an essentially isotropic process, while the physical size of the numerical domain is anisotropic. Kolmogorov's dissipation length therefore corresponds to fewer grid cells in the direction along the mean-field than perpendicular to it, and resolving it with the same accuracy in both directions is impossible. With our choice of dissipation coefficients, Kolmogorov's dissipation length ${{\ell }}_{{\rm{K}}}$ fulfills ${k}_{{{\rm{\Delta }}}_{\perp }}{{\ell }}_{{\rm{K}}}\sim 12$ and ${k}_{{{\rm{\Delta }}}_{\parallel }}{{\ell }}_{{\rm{K}}}\sim 1.2$, where ${k}_{{{\rm{\Delta }}}_{\perp }}=\pi /{\rm{\Delta }}x$ and ${k}_{{{\rm{\Delta }}}_{\parallel }}=\pi /{\rm{\Delta }}z$ are the wavenumbers corresponding to the grid resolution in perpendicular and parallel direction, respectively.

Driving turbulence with an energy-injection rate $\epsilon =0.1$ and cross-helicity injection rate $\sigma =0.0$ leads to a balanced steady-state with ${\sigma }_{c}={\mathcal{K}}/{\mathcal{E}}=0.0$, while changing the cross-helicity injection to $\sigma =0.8$ eventually yields a steady-state cross helicity of ${\sigma }_{c}=0.9$. The latter value is obtained from averaging the scalar product ${\mathcal{K}}=\langle {\boldsymbol{u}}\cdot \delta {\boldsymbol{b}}\rangle$, which saturates at ${\mathcal{K}}=0.9\ {\mathcal{E}}$ when the dissipation of energy and cross helicity matches their injection rates. Similarly, injecting positive or negative cross helicity at a lower rate ($\sigma =\pm 0.5$) results in a steady-state configuration with ${\sigma }_{c}=\pm 0.6$, respectively.

Figure 2 provides comparisons of steady-state spectra for ${\sigma }_{c}=0.0$ and ${\sigma }_{c}=0.9$. For the strongly imbalanced case (solid lines), the Elsasser energy with the same (in this case, positive) sign as the total cross-helicity dominates over the opposite Elsasser energy at low $k\lesssim 20$, while both Elsasser energies are in equipartition in the balanced case. The positive Elsasser energy is generally associated with waves propagating opposite to a strong magnetic mean-field. We deduce that the majority of the propagating shear Alfvén waves in the imbalanced case has a phase velocity directed opposite to the mean field (their energy is respresented by ${{\mathcal{E}}}^{+}$), while the wave distribution in the balanced case exhibits no such anisotropy.

Zoom In Zoom Out Reset image size

The spectra in the direction parallel to the mean field show a similar dominance of the Elsasser energy with the same sign as the cross-helicity (Figure 2). As analytical predictions and recent solar-wind measurements have shown (see Podesta 2009 and references therein), one expects a scaling ${\mathcal{E}}\propto {k}_{z}^{-2}$, and indeed this scaling fits our simulations well in the balanced case. In the cross-helical cases, the spectral indices of the positive and negative Elsasser energies differ slightly (Grappin et al. 1983).

The imbalance of the cross-helical cases is also obvious in plots of the spatial distribution of the cosine between the velocity field and the magnetic field ${\boldsymbol{b}}={{\boldsymbol{B}}}_{{\boldsymbol{0}}}+\delta {\boldsymbol{b}}$. The zero-cross-helicity case exhibits positive and negative values of $\mathrm{cos}\angle ({\boldsymbol{u}},{\boldsymbol{b}})$ in approximately equal proportion (Figure 3(a), (b)). For the strongly cross-helical case with ${\sigma }_{c}=0.9$, the velocity field shows a similarly weak correlation with the total magnetic field (Figure 3(c); however, note the small changes in the slope of the histogram curve at the extremal values $\mathrm{cos}\angle ({\boldsymbol{u}},{\boldsymbol{b}})=\pm 1$ in b), while positive values of $\mathrm{cos}\angle ({\boldsymbol{u}},\delta {\boldsymbol{b}})$ dominate if the alignment between the velocity field and only the fluctuations of the magnetic field is considered (Figure 3(d)). A similar plot of $\mathrm{cos}\angle ({\boldsymbol{u}},\delta {\boldsymbol{b}})$ for the balanced case is not depicted since it varies only negligibly from the total-field alignment already shown in Figure 3(a).

Zoom In Zoom Out Reset image size

Due to the almost identical distribution of the alignment of ${\boldsymbol{u}}$ and $({{\boldsymbol{B}}}_{{\boldsymbol{0}}}+\delta {\boldsymbol{b}})$ in balanced and cross-helical turbulence, both cases exhibit a similar level of isotropy of the motional electric field along the mean-field direction. The distribution of e_z is symmetric around e_z = 0 with an approximately exponential drop-off in both directions (Breech et al. 2003), with the global average $\langle {e}_{z}\rangle =0.000458\langle | {\boldsymbol{e}}| \rangle$ in the balanced case and $\langle {e}_{z}\rangle =-0.000034\langle | {\boldsymbol{e}}| \rangle$ in the strongly cross-helical case (${\sigma }_{c}=0.9$).

In order to investigate the presence of incompressible Alfvén waves, we switched off forcing and dissipation once the balanced and imbalanced field configurations had reached a steady-state. Setting ${{\boldsymbol{f}}}_{u}={{\boldsymbol{f}}}_{b}=\eta =\nu =0$, we ensure that the propagation of shear-Alfvén waves is not disturbed by non-ideal effects, which are not accounted for in the derivation of the quasilinear diffusion coefficients. Comparing the power spectra of the magnetic field Fourier-transformed in space and time (Figure 4), we find evidence of Alfvén waves traveling along and opposite the mean-field direction and therefore confirm that the evolution of the MHD turbulence is shaped mainly by propagating shear-Alfvén waves. The x-shaped dispersion plot of the balanced case is symmetric with respect to the k_z = 0 line, indicating that the waves propagating parallel and antiparallel to the mean field are equally strong. In the strongly cross-helical case (${\sigma }_{c}=0.9$), there is more energy in the modes along $\omega =-{v}_{{\rm{A}}}{k}_{z}$ than in the modes with the oppositely directed phase velocity along $\omega =\;+{v}_{{\rm{A}}}{k}_{z}$, particularly at larger k_z, indicating once more a dominance of Alfvén waves propagating opposite to the mean-field.

Zoom In Zoom Out Reset image size

We continue using zero forcing and dissipation while we compute the forces that test-particles experience in these turbulent steady-state fields. The trajectories of these test-particles are evolved using an implicit fourth-order Runge–Kutta solver with an adaptively determined timestep size. Their velocity at each timestep is calculated using the Lorentz force computed from the MHD magnetic field b and, in some runs, the motional electric field ${{\boldsymbol{e}}}_{\mathrm{mot}}=-{\boldsymbol{u}}\times {\boldsymbol{b}}$, as interpolated from the grid data using third-order spline functions.

At the starting time t = 0 of the test-particle simulation, ensembles of ${N}_{p}=2\times {10}^{4}$ particles are distributed randomly on the grid, with an energy equal to ${E}_{p}=v{(0)}^{2}/2$ and a velocity ${v}_{z}=\mu v(0)$ along the direction of the mean-field magnetic field. The initial pitch-angle cosine μ differs for each ensemble, taking values from −1.0 to 1.0. To demonstrate the interaction with propagating Alfvén waves, we initialize the test-particle ensembles with velocities close to the Alfvén velocity, using both $v(0)={v}_{{\rm{A}}}$ and $v(0)=3{v}_{{\rm{A}}}$. The charge-to-mass ratio of the test-particles is chosen such that their initial gyroradius is ${r}_{g}(0)=v(0)/({{qB}}_{0}/m)=0.12{\lambda }_{\perp }$ for Alfvén-speed particles and ${r}_{g}(0)=0.36{\lambda }_{\perp }$ for $v(0)=3{v}_{{\rm{A}}}$.

Since our goal is to compare analytical pitch-angle diffusion coefficients predicted by quasilinear theory with simulation results, we define μ as the cosine of the angle subtended by the particle velocity and the direction of the global mean magnetic field (the z direction) and thus follow the QLT convention. Using an only locally averaged mean-field direction as reference would yield a pitch-angle diffusion coefficient that would not fully conform to the quasilinear definition. QLT is derived by expanding the equations of motion in a small perturbation of a globally constant magnetic field, hence it is the global mean-field direction that is relevant for a comparison.

In order to compare the pitch-angle scattering of test-particles on short timescales with quasilinear estimates, we define

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mu \mu }(\mu (0))=\displaystyle \frac{1}{2}\langle \displaystyle \frac{{[\mu ({T}_{g})-\mu (0)]}^{2}}{{T}_{g}}\rangle \end{eqnarray} \tag{ 10 }$$

as the scattering rate after one gyroperiod. Here $\mu (0)$ and $\mu ({T}_{g})$ are the pitch-angle cosine at the test-particle injection time and after one gyroperiod, respectively. Whereas ${D}_{\mu \mu }$, as defined above in the context of QLT, corresponds to an infinite-time limit, we decided to measure ${{\rm{\Delta }}}_{\mu \mu }$ instead for several reasons. Although calculating ${{\rm{\Delta }}}_{\mu \mu }$ for timescales significantly longer than T_g could be expected to improve the agreement with quasilinear calculations of ${D}_{\mu \mu }$, we will demonstrate that this is not the case, due to various effects which were not included in the derivation of the quasilinear coefficients above, but which are inevitably present both in our simulations and in realistic MHD turbulence. In the Appendix, we demonstrate that ${{\rm{\Delta }}}_{\mu \mu }$ is bounded from below by ${D}_{\mu \mu }$.

First, as the pitch-angle cosine μ is limited to the interval $[-1,+1]$, its evolution is subject to boundary effects. Measuring the mean-square pitch-angle deviation

$$\begin{eqnarray}&&{\rm{\Delta }}{\mu }^{2}(t)=\langle {[\mu (t)-\mu (0)]}^{2}\rangle \end{eqnarray} \tag{ 11 }$$

for large t will, in general, underestimate the diffusion coefficient ${D}_{\mu \mu }(\mu (0))$ at the initial $\mu (0)$ as some test-particles will have been reflected at the boundaries. Even before this reflection occurs, however, ${D}_{\mu \mu }(\mu )$ varies strongly with μ and is therefore difficult to determine if the pitch-angle distribution spreads too quickly. This is particularly relevant if the turbulence amplitude is too large (for instance $\delta B/{B}_{0}\sim 1$). In that case, even particles with a gyroradius as small as ${r}_{g}(0)=0.36{\lambda }_{\perp }$ isotropize almost completely within a single gyroperiod. Therefore, we only present results for an intermediate turbulence strength of $\delta B/{B}_{0}=0.1$.

Second, our goal is to investigate the role of electric fields, which accelerate or decelerate test-particles on longer timescales. Since the Fokker–Planck Equation (1) does not account for these effects, good agreement can only be expected if the timescales on which ${{\rm{\Delta }}}_{\mu \mu }$ is measured are short enough that the change in kinetic energy is negligible.

Most importantly, we will demonstrate that pitch-angle scattering on longer timescales is severely constrained by the approximate conservation of the magnetic moment,

$$\begin{eqnarray}&&{M}_{{\rm{mag}}}=m\displaystyle \frac{{v}_{\perp }^{2}}{2B},\end{eqnarray} \tag{ 12 }$$

the first adiabatic invariant in slowly varying magnetic fields (Landau & Lifshitz 1982). The standard formulation of quasilinear theory does not account for the approximate conservation of ${M}_{{\rm{mag}}}$ (Goldstein 1980; Jaekel 1998; see however Schlickeiser & Jenko 2010). As an adiabatic invariant, ${M}_{{\rm{mag}}}$ needs only be considered for the evolution of ${\rm{\Delta }}{\mu }^{2}$ on timescales exceeding several gyroperiods, so that ${{\rm{\Delta }}}_{\mu \mu }$, which is computed after only one gyroperiod, is not affected by these constraints and therefore allows for a better comparison with quasilinear predictions. Although the absolute value of ${{\rm{\Delta }}}_{\mu \mu }$ will not necessarily be close to the value of ${D}_{\mu \mu }$ in all cases, the scaling with the pitch-angle cosine μ, which is our primary concern in these investigations of cross-helical turbulence, will be almost identical.

We begin with the most realistic scenario, fully electrodynamic turbulence, in which the evolution of the MHD fields is computed in parallel with the test-particle trajectories. The initial conditions correspond to the field configurations presented in Section 3.1. In Figure 5, we show the evolution of the kinetic energy, the magnetic moment, and the mean-square pitch-angle deviation for two representative ensembles with $v(0)=3{v}_{{\rm{A}}}$ and $\mu (0)=\pm 0.6$ in balanced and strongly imbalanced turbulence.

Zoom In Zoom Out Reset image size

For both ensembles in the balanced-turbulence case, energy and the averaged magnetic moment oscillate with the gyrofrequency, but the values of the minima attained after each full gyroperiod show only negligible variation. The mean-square deviation of ${M}_{{\rm{mag}}}$, as indicated by the shaded regions, converges to about 12% of the initial value as the magnetic moment is conserved to a large degree because of the relatively weak perturbation amplitude $\delta B/{B}_{0}$. This puts a limit to the evolution of ${\rm{\Delta }}{\mu }^{2}(t)$, however, since the kinetic energy, the magnetic moment, and the pitch-angle cosine are approximately connected via

$$\begin{eqnarray}&&E\approx {M}_{{\rm{mag}}}B+E{\mu }^{2}.\end{eqnarray} \tag{ 13 }$$

This relation is valid if we assume that the magnetic turbulence is so weak that the difference between the averaged velocity components $\langle {v}_{\parallel }^{2}\rangle$ along the local magnetic-field direction and $\langle {v}_{z}^{2}\rangle$ along the global mean-field can be neglected, the difference being of order ${(\delta B/{B}_{0})}^{2}$. Consequently, the minima of ${\rm{\Delta }}{\mu }^{2}(t)$ after each gyroperiod increase in value for two gyroperiods before settling on an almost constant value.

In imbalanced turbulence with ${\sigma }_{c}=0.9$, the test-particle ensemble with a positive pitch-angle of $\mu (0)=+0.6$ evolves similarly to the balanced case. However, the counter-propagating particles with $\mu (0)=-0.6$ are much closer to resonance with the dominating Alfén-wave population, which is traveling opposite to the mean-field direction as well. Consequently, these particles undergo deceleration as a result of inverse Landau damping and their kinetic energy exhibits a steady decrease (red lines in Figure 5). While the magnetic moment is approximately conserved again, the variation in energy makes it possible for pitch-angle scattering to continue on longer timescales; hence, the minima of ${\rm{\Delta }}{\mu }^{2}(t)$ after each gyroperiod grow steadily. Obviously, a comparison of the mean-square pitch-angle deviations for co- and counter-propagating ensembles on these longer times would be misleading. Instead, as described above, we compare the values of the first minima of ${\rm{\Delta }}{\mu }^{2}$ (or rather the slopes of the dotted lines through these minima).

These values of ${{\rm{\Delta }}}_{\mu \mu }$ for 21 different initial pitch-angles ($\mu (0)\in \{-1.0,-0.9,\ldots ,+1.0\}$) are then used to determine a best-fit parameter set $\{\delta {B}^{(\mathrm{fit})},{\sigma }^{(\mathrm{fit})},{s}^{+},{s}^{-}\}$ for the quasilinear prediction for ${D}_{\mu \mu }$ given by Equation (5). The remaining parameters (${k}_{\mathrm{min}},{{\rm{\Omega }}}_{0},\gamma$) were fixed at their actual values before the optimization was initialized. Although Equation (5) was derived under the assumption of a highly simplified slab spectrum with a constant power index, hardly what our simulations of realistic MHD exhibit in Figure 2, it will become clear that the agreement of ${{\rm{\Delta }}}_{\mu \mu }$ and the quasilinear predictions of ${D}_{\mu \mu }$ is rather striking. The results from a Levenberg–Marquardt fit using our measurements of ${{\rm{\Delta }}}_{\mu \mu }(\mu )$ are shown in Figure 6.

Zoom In Zoom Out Reset image size

The fit-values of ${\sigma }^{(\mathrm{fit})}$ match the values of σ used for the cross-helicity injection well in all cases. A decrease of the fit quality if one Elsasser energy is very small (that is, for $\sigma =0.8$) is visible but not surprising. The relative turbulence amplitudes $\delta {B}^{(\mathrm{fit})}/{B}_{0}$ are close to the actual value $\delta B/{B}_{0}=0.1$, usually lying slightly below because the particle scattering is mainly due to only the slab component of the turbulence. The spectral indices obtained from the fit procedure are almost always larger than the predictions for the inertial-range spectrum in Kolmogorov's or Kraichnan's theories (5/3 or 3/2, respectively) since, like realistic turbulence, the spectra in our simulations steepen at the transition to the dissipation range (see Figure 2).

The graphs of the quasilinear ${D}_{\mu \mu }$ for these parameter sets perfectly match the scattering rates ${{\rm{\Delta }}}_{\mu \mu }$ determined from the test-particle trajectories, as shown in Figure 7. Whereas the graphs in the balanced case are symmetric with respect to the sign of $\mu (0)$, particles in the imbalanced cases are scattered faster within the first gyroperiod if they travel in the direction of the weaker Alfvén-wave population (if $\mu (0)\gt 0$ for $\sigma \gt 0$ and vice versa).

Zoom In Zoom Out Reset image size

This asymmetry is shown more directly in Figure 8, which contains pitch-angle histograms of particles with $v(0)=3{v}_{{\rm{A}}}$ after ten gyroperiods for test-particles initially propagating exactly parallel or anti-parallel to the mean-field direction. In this case, both the heating that was observed for the negative-μ ensemble in Figure 5(d) and the initial magnetic moment are insignificantly small, and the behavior predicted by quasilinear theory can be observed for much longer. While both ensembles, with $\mu (0)=+1$ and with $\mu (0)=-1$, are scattered equally fast in balanced MHD turbulence, positive cross helicity leads to particles propagating along the mean-field direction being scattered much faster than in the opposite case.

Zoom In Zoom Out Reset image size

In Figure 9, we compare similar histograms of sample pitch-angle distributions after one gyroperiod with solutions that we directly obtained from the Fokker–Planck Equation (1), using $f(\mu ,t=0)=\delta (\mu (0)-\mu )$ as initial condition and assuming a homogeneous density along the z-coordinate. While we showed above that the mean-square displacement of the pitch-angle cosine after one gyroperiod agrees with the scaling of ${D}_{\mu \mu }(\mu )$ predicted by QLT, Figure 9 demonstrates that the shape of the histograms agrees with a diffusive broadening of the pitch-angle distribution within the first gyroperiod. This confirms once more that the quasilinear prediction (5) gives an accurate description of the pitch-angle scattering on short timescales even in evolving imbalanced turbulence with electric fields.

Zoom In Zoom Out Reset image size

There are two possible causes for this asymmetry in imbalanced turbulence: the coherent interaction with propagating Alfvén waves and the spatial structure of the turbulent fields. To identify the more important mechanism, we first ignore the electric-field component of the Lorentz force and propagate the test-particles using only

$$\begin{eqnarray}&&{\dot{{\boldsymbol{v}}}}_{i}=\alpha {{\boldsymbol{v}}}_{i}\times {\boldsymbol{b}}({{\boldsymbol{x}}}_{i},t).\end{eqnarray} \tag{ 14 }$$

In the "magnetodynamic" runs presented in this subsection, the MHD fields were still dynamically evolved in parallel with the test-particles. Since the acceleration is perpendicular to the velocity at each time step, the kinetic energy of each test-particle is perfectly conserved. As shown in Figure 10, the evolution of both ensembles in the balanced case and of the positive-μ ensemble in the imbalanced case is similar to the fully electrodynamic setting presented above.

Zoom In Zoom Out Reset image size

However, the averaged magnetic moment of the ensemble with $\mu (0)=-0.6$ visibly increases on longer timescales. Similarly, the distribution of ${M}_{{\rm{mag}}}$ continues to broaden for at least ten gyroperiods, as indicated by the red shaded region in Figure 10(c). Although the ensemble in Landau resonance with the dominant Alfvén-wave population cannot be decelerated anymore because electric acceleration is not considered, the coherent interaction with the magnetic component of the Alfvén waves is now able to violate the approximate conservation of the magnetic moment. Energy is transfered from the v_z-component into the perpendicular velocity components. Without a constant adiabatic invariant, the pitch-angle cosine can decrease as before, even though the kinetic energy is fixed, and ${\rm{\Delta }}{\mu }^{2}(t)$ keeps increasing at a similar rate as in the fully electrodynamic runs.

We use the first minima of ${\rm{\Delta }}{\mu }^{2}$ for $\mu (0)\;\in \{-1.0,-0.9,\ldots ,+1.0\}$ to determine a best-fit parameter set with the quasilinear ${D}_{\mu \mu }$ in magnetodynamic turbulence for the same variables as before, but with ${D}_{\mu \mu }$ given by Equation (7). Most of the results are similar to those presented in the previous subsection, so we only show the best-fit cross-helicity values ${\sigma }^{(\mathrm{fit})}$ in Figure 11. The results for the Alfvén-speed ensemble ($v(0)={v}_{{\rm{A}}}$) are highly inaccurate, while the ensemble with $v(0)=3{v}_{{\rm{A}}}$ at least exhibits the correct trend: ${\sigma }^{(\mathrm{fit})}$ grows weakly as σ increases.

Zoom In Zoom Out Reset image size

Plotting ${D}_{\mu \mu }$ and ${{\rm{\Delta }}}_{\mu \mu }$ in Figure 12, we see that, without electric fields, the scattering asymmetry is strongly reduced compared to the graphs in Figure 7. Electric fields obviously play a crucial role in establishing this asymmetry in imbalanced turbulence.

Zoom In Zoom Out Reset image size

To determine the relative importance of wave propagation, we have also performed test-particle simulations with single static snapshots of the magnetic-field configurations. In these magnetostatic runs, we calculate the Lorentz force on the particles using only the magnetic field, similar to the previous subsection (Equation (14)) but without evolving the field in time.

Since particles can no longer interact coherently with traveling waves, Landau resonance as observed before is impossible and pitch-angle scattering proceeds at similar rates for positive and negative $\mu (0)$ even in strongly cross-helical turbulence. As shown in Figure 13, the mean magnetic moment stays constant throughout the simulations and ${\rm{\Delta }}{\mu }^{2}(t)$ settles on an approximately constant value after a few gyroperiods.

Zoom In Zoom Out Reset image size

Evolving the trajectories for longer times, which is possible since magnetostatic test-particle simulations are not very demanding compared to fully electrodynamic runs, yields a subdiffusive behavior of the pitch-angle distribution. Thus, we computed ${{\rm{\Delta }}}_{\mu \mu }$ after one gyroperiod again to perform a best-fit operation, varying $\delta B$ and s in Equation (2). As we found statistical noise to be significantly larger in static runs than in the evolving-MHD runs presented above, we show a 10%-error bar in the plots of ${{\rm{\Delta }}}_{\mu \mu }(\mu )$. Repeating each test-particle run in ten different static snapshots of the magnetic field for otherwise identical parameters, we found that the values we obtained for ${{\rm{\Delta }}}_{\mu \mu }$ varied by this amount, whereas the evolving-MHD runs only varied by about 2%.

Taking this statistical noise into account, the graphs of ${{\rm{\Delta }}}_{\mu \mu }(\mu )$ are symmetrical with respect to the sign of μ (see Figure 14). Cross helicity is therefore irrelevant in magnetostatic turbulence, and the spatial structure of the magnetic fields can be ruled out as a contributor to the scattering asymmetry. As a sidenote, we remark that none of our results were affected by the sign of the test-particle charge, proving that magnetic-helicity contributions to the scattering asymmetry, which would be sensitive to the particle charge, are insignificant.

Zoom In Zoom Out Reset image size

Adding the electric-field component back into the Lorentz force while still using only static field snapshots, we found that this asymmetry can be re-established (Figure 15). While the absolute scattering rates ${{\rm{\Delta }}}_{\mu \mu }$ are reduced to about 70% compared to the values for the fully electrodynamic runs in Figure 7, the shapes of the graphs are almost identical to their fully electrodynamic equivalents in each case, especially if one takes the statistical noise into account.

Zoom In Zoom Out Reset image size

Since the electromagnetic fields of the Alfvén waves are fixed in time and space in these runs, it is the spatial structure of the electric fields that must be the predominant soure of the scattering asymmetry. Hence, we conclude that particles in fully evolving MHD turbulence with non-zero cross helicity are scattered differently for positive and negative pitch-angles not only because of coherent wave–particle interactions, but additionally because the spatial structure of electric fields in imbalanced turbulence favors such an asymmetry. Of course, this structure can be viewed as the result of the imbalance between Alfvén-waves propagating parallel and anti-parallel to the magnetic mean-field direction, a property that is conserved when the fields are frozen.