**1. Introduction**

The occurrence of highly energetic particles is a ubiquitous feature in space plasmas (e.g., in the ionosphere, the auroral zone, solar wind, and at the mesosphere, etc.) and in laboratory plasmas [1–12]. The velocity distribution in such plasmas may deviate from the usual thermal Maxwellian distribution, developing a long-tail for high-velocity arguments due to an excess in the fast (superthermal) part of the population; such a behavior is effectively modeled by a (kappa) *κ*-type distribution function [1,4,12–19]. The kappa distribution function was initially postulated by Vasyliunas [1] in an effort to reproduce the observed power-law dependence at high velocities [17,20,21]. Since then, a large number of studies adopted kappa distributions, combining theoretical [18,19,22,23], computational [24], and even experimental [25,26] approaches to study the effect of superthermal particle populations on wave dynamics in different plasma environments.

Particle "trapping", i.e., the fact that a portion of, for example, the electron population remains confined in a finite region—thus generating vortices—in phase space, is an intrinsic characteristic of plasma dynamics, often overlooked in studies based on basic fluid theory. Phase-space structures, known as "electron-holes" are thus formed due to particles trapped in the wave potential. This mechanism, initially predicted via kinetic theory [27–29], was later observed in space and in the laboratory [30–36], and it was also shown to occur spontaneously in computer simulations [37]. Of particular relevance to current study is the fact that Simpson et al. [38] reported the presence of trapped electrons in the Saturnian magnetic field, an environment characterized by the existence of *κ*-distributed electrons with values of *κ* 2–4, as confirmed by Schippers et al. [34]. It is therefore important to consider the effect of particle nonthermality and trapping effect simultaneously to explore the properties of different electrostatic modes.

**Citation:** Sultana, S.; Kourakis, I. Dissipative Ion-Acoustic Solitary Waves in Magnetized *κ*-Distributed Non-Maxwellian Plasmas. *Physics* **2022**, *4*, 68–79. https://doi.org/ 10.3390/physics4010007

Received: 21 October 2021 Accepted: 3 January 2022 Published: 20 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

As regards the theoretical modeling of particle trapping, Schamel's original papers [27,28] showed that trapped particles led to a vortex-like electron distribution, and the kinetic model was shown to be associated with a modified version of the known integrable Korteweg–de Vries (KdV) partial differential equation. The "Schamel equation" [28,39–42] describes the evolution of nonlinear electrostatic waves under the influence of a fractional nonlinearity (in contrast with the standard KdV theory where quadratic nonlinearity is dominant). A number of theoretical studies followed [43–49] in an effort to investigate the properties of nonlinear waves (solitary waves, shocks) in the presence of trapped particles, using first principles.

The combined effect of electron superthermality and phase-space trapping was first considered by Williams et al. [40], who adopted the Schamel equation approach to model and characterize ion-acoustic solitary waves in an unmagnetized electron-ion plasma with *κ*-distributed electron populations subject to trapping. Following that study, the combined effect of electron superthermality and trapping was considered by Sultana and coworkers [41,42] (on ion-acoustic modes in collisionless plasmas) and by Hassan et al. [50], who investigated the nonlinear features of electron-acoustic waves in a magnetized plasma and considered the combined effect of electron trapping and electron superthermality. The study led by Hassan et al. [50] focused on electron-acoustic waves, a mode known to occur exclusively in the simultaneous presence of two distinct electron populations (usually referred to as the 'cold' and the 'hot' electrons), as it relies in fact on the inertia being provided by the cold electron component and the restoring force being provided by their hot counterpart. The associated (electronic) dynamical frequency scales are clearly distinct from the (slower, ionic) scales that are typical of the study presented here.

To our best knowledge, there is no rigorous and systematic study of the nonlinear propagation of the ion-acoustic waves in a magnetized *collisional* plasma in the presence of trapped *κ*-nonthermal electrons. The investigation at hand is therefore an attempt to fill in this gap by presenting a rigorous and systematic study of the characteristics of ion-acoustic waves propagating in a magnetized *κ*-nonthermal plasma [41], taking into account the combined impact of electron trapping and of a suprathermal electron distribution, in account of the inherent plasma collisionality. The main focus here is to investigate the influence of particle trapping on the dynamics of dissipative solitary waves, and also to analyze the effect of the ambient magnetic field and its interplay with wave damping and how these affect the characteristics of obliquely propagating ion-acoustic solitary excitations.

This article is organized as follows. The basic formalism is presented in the following Section 2. A dissipative version of the Schamel equation is derived via a multiscale perturbative approach, and the detail about the nonlinear, dispersion, and dissipative term, is discussed in Section 3. The propagation nature (basic features) of dissipative ion-acoustic waves for different relevant plasma (configuration) parameters is studied numerically in Section 4. Finally, the results obtained are summarized in the concluding Section 5.

#### **2. Basic Plasma-Fluid-Dynamic Formalism**

An electron-ion plasma is considered here being embedded in a uniform magnetic field directed along the *z*-axis, i.e., **B**<sup>0</sup> = *B*0*z*ˆ. Due to their large mass (relative to the electrons), inertial ions are modeled as a cold fluid, i.e., their thermal pressure is neglected for simplicity. At the ionic scale of interest, the electron inertia may be neglected: the electrons are assumed to deviate from thermal equilibrium, and hence, a *κ*-type distribution will be explicitly adopted to model their distribution. For the purpose of this analysis, the combined effect of electron trapping and superthermality is considered, following the steps outlined in Ref. [40].

Charge neutrality at equilibrium imposes: *zini*<sup>0</sup> − *ne*<sup>0</sup> = 0, where *ni*<sup>0</sup> and *ne*<sup>0</sup> denote the unperturbed ion and electron number densities, respectively, while *zi* is the charge state of the ion component (e.g., 1, 2,. . .; the value of *zi* is left arbitrary here, for generality).

We are interested in modeling the dynamics of ion-acoustic excitations, whose phase speed *vph* may well exceed the ion thermal speed, but is far smaller than the electron thermal speed. The following fluid evolution equations are considered:

$$\frac{\partial n\_i}{\partial T} + \nabla . (n\_i \mathbf{u}\_i) = 0 \,, \tag{1}$$

$$\frac{\partial \mathbf{u}\_i}{\partial T} + (\mathbf{u}\_i, \nabla) \mathbf{u}\_i = -\frac{z\_i e}{m\_i} \nabla \Phi$$
 
$$z\_i e B\_0 \ll \delta$$

$$+\frac{z\_i e B\_0}{m\_i} (\mathbf{u}\_i \times \boldsymbol{\varepsilon}) - \nu\_i \mathbf{u}\_i \,\,\,\,\tag{2}$$

$$
\nabla^2 \Phi = 4\pi e (n\_\varepsilon - z\_i n\_i) \,, \tag{3}
$$

where *ni* (*ne*) denotes the number density of the ion (electron) species, **u***<sup>i</sup>* is the ion fluid speed, Φ is the electrostatic wave potential (all these quantities are dynamic functions of space and time), and *e* is electron charge. An ad hoc damping term is introduced in the fluid equation of motion (momentum conservation equation) to account for, e.g., ion-neutral collisions; the collision frequency *ν<sup>i</sup>* was defined to this effect.

The combined effect of electron trapping and deviation from Maxwell-type equilibrium was studied analytically in Ref. [40]; the tedious algebraic procedure need not be reproduced here, but the main steps are summarized. A modified *κ*-distribution function, effectively taking into account the trapped part of the electron population (i.e., for electrons trapped in the wave potential if their energy *Ee* < 0) is given by [40]

$$f(v, \phi) = \frac{\Gamma(\kappa)}{\sqrt{2\pi} (\kappa - 3/2)^{1/2} \, \Gamma(\kappa - 1/2)}$$

$$\times \quad \left[1 + \beta \left(\frac{v^2/2 - \phi}{\kappa - 3/2}\right)\right]^{-\kappa} \text{ for } E\_{\mathfrak{e}} \le 0. \tag{4}$$

Here, *v* is the velocity, *φ* is the electrostatic potential, and *κ* is the superthermality index (measures deviation from the Maxwell–Boltzmann distribution), while *β* (<1) quantifies the efficiency of electron trapping. The known vortex-type distribution for trapped Maxwellian electrons is recovered in the limit *κ* → ∞. The number density of the electrons is obtained by integration as [40]

$$n\_{\varepsilon}(\phi) \quad = \int\_{-\infty}^{-\sqrt{2\phi}} f\_{\varepsilon}^{\chi}(v, \phi) \, dv + \int\_{-\sqrt{2\phi}}^{\sqrt{2\phi}} f(v, \phi) \, dv$$

$$+ \quad \int\_{\sqrt{2\phi}}^{\infty} f\_{\varepsilon}^{\chi}(v, \phi) \, dv,\tag{5}$$

where *f <sup>κ</sup> <sup>e</sup>* (*v*, *φ*) is the *κ*-distribution function for the free electrons; details can be found in Ref. [17].

By normalizing all variables, one obtains the following system of (dimensionless) equations:

$$
\frac{
\partial \mathbf{u}
}{
\partial t
} + \nabla .(\nu \mathbf{u}) = 0 \,, \tag{6}
$$

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \boldsymbol{\nabla})\mathbf{u} = -\mathbf{\tilde{\nabla}}\phi + \Omega\_{\mathbf{c}}(\mathbf{u} \times \boldsymbol{\varepsilon}) - \nu \mathbf{u} \,, \tag{7}$$

$$
\bar{\nabla}^2 \phi \simeq 1 - n + a\_1 \phi + a\_2 \phi^{3/2} + a\_3 \phi^2 \, , \tag{8}
$$

where *n* = *ni*/*ni*0, *u* = [*miu*<sup>2</sup> *<sup>i</sup>* /(*ziTe*)]1/2 with *mi* being the mass of ion, *Te* being the electron temperature (Boltzmann's constant *kB* is omitted where obvious), the electrostatic potential *φ* = *e*Φ/*Te*, *λ<sup>D</sup>* = (*Te*/4*πe*<sup>2</sup>*zini*0)1/2, *t* = *ωp*,*iT* (where *ωp*,*<sup>i</sup>* = (4*πe*2*z*<sup>2</sup> *<sup>i</sup> ni*0/*mi*)1/2 is the ion plasma frequency, and *T* is the inverse of the ion plasma frequency), Ω*<sup>c</sup>* = *ωc*,*i*/*ωp*,*<sup>i</sup>*

(where *ωc*,*<sup>i</sup>* = *zieB*0/*mi*) is the reduced cyclotron frequency, and *ν* = *νi*/*ωp*,*<sup>i</sup>* (note that all frequencies were scaled by the ion plasma frequency for convenience). The information related to electron trapping, for *κ*-distributed electrons, is "hidden" in the coefficients *a*1, *a*2, *a*<sup>3</sup> entering the normalized expression of Poisson Equation (8), which are given by [40]

$$a\_1 = \frac{2\kappa - 1}{2\kappa - 3}, \quad a\_2 = \frac{8\sqrt{2/\pi}(\beta - 1)\kappa}{3(2\kappa - 3)^{3/2}} \frac{\Gamma(\kappa)}{\Gamma(\kappa - 1/2)},$$

$$a\_3 = \frac{4\kappa^2 - 1}{2(2\kappa - 3)^2}.\tag{9}$$

Once a solution for the electrostatic potential *φ* is formally obtained, the trapped electron population's density obtained from Equation (5) can be expressed as [40]

$$m\_{\varepsilon} \simeq 1 + a\_1 \phi + a\_2 \phi^{3/2} + a\_3 \phi^2 \cdot \cdots \cdot,\tag{10}$$

which to be substituted into Poisson Equation (8). The known analogous expression for the trapped electron number density in the case of Maxwellian plasma [28] is readily recovered here, upon considering the limit *κ* → ∞ in the latter relation. On the other hand, the limit of Equation (10) as *β* → 1 leads to the classical expression for *κ*-distributed electrons; see e.g., Ref. [51] and elsewhere. Finally, the Maxwellian limit for free electrons, viz. *a*<sup>1</sup> = 2*a*<sup>3</sup> = 1, is recovered by considering *β* → 1 and *κ* → ∞, reducing the electron density dependence to e*<sup>φ</sup>* <sup>1</sup> <sup>+</sup> *<sup>φ</sup>* <sup>+</sup> *<sup>φ</sup>*2/2, as expected.

The closed system of Equations (6)–(8), describes the evolution of the plasma (fluid) state variables, and forms the basis of the analysis here.

#### **3. A Schamel Equation for Damped Ion-Acoustic Waves (IAWs)**

To model small-amplitude ion-acoustic excitations (dissipative solitary waves) within the model under consideration, one needs to proceed by defining a set of stretched coordinates as

$$\xi = \epsilon^{1/4} \left( l\_x \pounds + l\_y \gplus + l\_z \pounds - v\_p t \right), \quad \tau = \epsilon^{3/4} t,\tag{11}$$

where (1) is a small parameter that measures the strength of the nonlinearity, *vp* is a constant to be determined (in fact, representing the phase speed, scaled by the ion sound speed, *c*<sup>0</sup> = (*ziTe*/*mi*)1/2), and *lx*, *ly* and *lz*, are directional cosines of the wave vector **k** along the *x*, *y* and *z* axes, respectively (for instance, *lz* = (**k** · *z*ˆ)/*k*), hence *l* 2 *<sup>x</sup>* + *l* 2 *<sup>y</sup>* + *l* 2 *<sup>z</sup>* = 1. Let us recall that the position variables *x*, *y* and *z* are all normalized by *λD*, while *τ* is normalized by the ion plasma period *ω*−<sup>1</sup> *<sup>p</sup>*,*<sup>i</sup>* . The above *Ansatz*, which was first introduced in Ref. [28] and then later adopted by various authors (e.g., [40,52]), essentially describes a Galilean transformation into a slowly varying moving frame, wherein the time variation of the structure is even slower in time.

The dependent variables *n*, **u** and *φ* may now be expanded near the equilibrium states as power series of as follows:

$$\begin{cases} n = 1 + \epsilon n\_1 + \epsilon^{3/2} n\_2 + \dots, \\ u\_x = \epsilon^{5/4} u\_{1,x} + \epsilon^{3/2} u\_{2,x} + \dots, \\ u\_y = \epsilon^{5/4} u\_{1,y} + \epsilon^{3/2} u\_{2,y} + \dots, \\ u\_z = \epsilon u\_{1,z} + \epsilon^{3/2} u\_{2,z} + \dots, \\ \phi = \epsilon \phi\_1 + \epsilon^{3/2} \phi\_2 + \dots \end{cases} \tag{12}$$

To close the series expansion of the variables, a weak dissipation [53–55] to be considered due to ion-neutral collisions by assuming that the damping coefficient scales as *ν* = 3/4*ν*0.

Let us now proceed by substituting expansions (11) and (12) into the considered fluid plasma model Equations (6)–(8) and collecting various terms arising in each order in . The phase speed *vp* is obtained as a compatibility constraint upon considering the lowest order contributions in from each of the equations; the resulting expression reads:

$$
v\_p = l\_z / \sqrt{a\_1} \,. \tag{13}$$

This expression for the phase speed *vp* depends on the angle *θ* (via *lz* = cos *θ*) and on the electron superthermality index *κ*. Considering parallel propagation (*lz* = 1), the known expression for the ion sound speed in nonmagnetized plasma [22,40] is recovered as expected. Recovering dimensions for a minute, for physical transparency, Equation (13) leads to

$$V\_p = \frac{\omega}{k} = \lambda\_D \,\omega\_{p,i} \, \frac{l\_z}{\sqrt{a\_1}} = \left(\frac{z\_i T\_\varepsilon}{m\_i}\right)^{1/2} \left(\frac{2\kappa - 3}{2\kappa - 1}\right)^{1/2} l\_z \,\tag{14}$$

where *ω*, *k* and *Vp* here denote the wave (angular) frequency, the wavenumber and the phase speed (in the dimensional form), respectively. The acoustic ("sound") speed is thus recovered for infinitely large *κ*, while a lower value (i.e., predicting slower solitary waves) is predicted for small *κ*, in agreement with earlier theoretical predictions [22] and with space observations [21].

The variation of the ion-acoustic phase speed, *vp*, versus the electron's superthermality index, *κ*, is depicted in Figure 1, suggesting a slower phase speed in a plasma with significant portion of the electrons in the superthermal region (i.e., lower values of *κ*), when compared with the case of thermal (Maxwellian) electron. The phase speed is higher for parallel propagation than for oblique propagation, as shown in Figure 1. As expected, the curve tends to unity, asymptotically (*vp* → 1) for infinite kappa and for parallel propagation, prescribing the acoustic speed (in electron-ion plasma) as the phase speed in the Maxwellian limit.

**Figure 1.** Phase speed, *vp*, versus electron superthermality index, *κ*, and obliqueness angle, *θ* = cos−<sup>1</sup> *lz*, where *lz* is a directional cosine of the wave vector **k** along *z*-axis.

The perpendicular (*x* and *y*) components of the electric field related drift of the ion fluid, in terms of the electric potential *φ*1, can be obtained by separating the *y* and *x*components of the momentum equation, respectively, as

$$
\mu\_{1,\chi} = -\frac{l\_y}{\Omega\_c} \frac{\partial \phi\_1}{\partial \xi} \,, \tag{15}
$$

*∂ξ* . (16)

and *<sup>u</sup>*1,*<sup>y</sup>* <sup>=</sup> *lx*

Following an analogous procedure, the parallel (*z*) component of the ion fluid velocity is obtained as

Ω*<sup>c</sup> ∂φ*<sup>1</sup>

$$
\mu\_{1,z} = \frac{l\_z}{v\_p} \phi\_{1,z} \tag{17}
$$

Finally, in leading order, the perturbation of the ion density *<sup>n</sup>* <sup>1</sup> <sup>+</sup> *n*<sup>1</sup> <sup>+</sup> <sup>O</sup>(3/2) is obtained as

*n*<sup>1</sup> = *lz vp* 2 *φ*<sup>1</sup> . (18)

The next order in (obtained upon separating 3/2 from the momentum equation) leads to the *x* and *y*-components of the second order drift velocity of the ion fluid in the form

$$
\mu\_{2,x} = \frac{l\_x v\_p}{\Omega\_c^2} \frac{\partial^2 \phi\_1}{\partial \xi^2} \,, \tag{19}
$$

$$
\mu\_{2,y} = \frac{l\_y \upsilon\_p}{\Omega\_c^2} \frac{\partial^2 \phi\_1}{\partial \tilde{\xi}^2} \,'\, \tag{20}
$$

Following the same procedure (i.e., separating coefficient of 7/4 from the continuity and the *z*-component of the momentum equations, and then 3/2 from the Poisson equation) and eventually eliminating *n*2, *u*2,*z*, and *φ*2, one is led to a nonlinear partial differential equation (PDE) in the form

$$\frac{\partial \Psi}{\partial \tau} + A\psi^{1/2} \frac{\partial \Psi}{\partial \xi} + B \frac{\partial^3 \Psi}{\partial \xi^3} + \mathcal{C}\psi = 0 \,,\tag{21}$$

where, for brevity, the leading contribution of the electrostatic potential is denoted by *ψ* = *φ*1.

Equation (21), which bears the structure of the original Schamel equation [28], with the addition of the last term (arising due to collisions being taken into account), represents an evolution equation for an electrostatic potential disturbance, *<sup>φ</sup> ψ* <sup>+</sup> <sup>O</sup>(3/2), in a region where the trapped electrons are present. The algebraic scheme implied is obvious: once *ψ* is obtained from Equation (21), the leading contributions for the ion density and for the ion fluid speed (three) components can be obtained from (four) Equations (15)–(18).

The nonlinearity coefficient *A*, which is responsible for wave steepening, is given by

$$A = -\frac{3}{4} \frac{v\_p^3}{l\_z^2} a\_2 = -\frac{3}{4} \frac{a\_2}{a\_1^{3/2}} l\_z \,. \tag{22}$$

The nonlinearity dependence enters via both *θ* and *κ*, as expected: this is seen in Figure 2a.

On the other hand, the dispersion coefficient *B*—which is responsible for wave broadening—is given by

$$B = \frac{\upsilon\_p^3}{2l\_z^2} \left( 1 + \frac{1 - l\_z^2}{\Omega\_c^2} \right). \tag{23}$$

The expression for the coefficient *B* can be simplified upon setting *v*<sup>3</sup> *<sup>p</sup>*/2*l* 2 *<sup>z</sup>* = *lz*/2*a*3/2 <sup>1</sup> , showcasing the dependence of *B* on the propagation angle (via *lz*) and on *κ* (via *a*1), as shown in Figure 2b. The influence of the magnetic field (via Ω*c*) disappears in the case of parallel propagation (*lz* = 1), thus recovering a one-dimensional damped Schamel equation for unmagnetized plasma (this was intuitively expected, since the Larmor force has no component in the direction of the magnetic field, and thus does not affect parallel wave propagation).

Finally, the dissipative term *C* is given by

$$
\mathbb{C} = \frac{\nu\_0}{2},
\tag{24}
$$

as imposed by compatibility requirements (i.e., balancing various terms occurring in the same order in ).

**Figure 2.** Nonlinearity term *A* (**a**) and dispersion term *B* (**b**) versus *κ* and obliqueness angle *θ* for *β* = 0.5 and Ω*<sup>c</sup>* = 0.3, where *β* denotes the efficiency of electron trapping and Ω*<sup>c</sup>* is the ratio of the the reduced cyclotron frequency to the ion plasma frequency.

Interestingly, both *A* and *B* vanish for perpendicular propagation (i.e., for *lz* = cos(*π*/2) = 0), as seen in Figure 2. On the other hand, considering parallel propagation (*lz* = cos(0) = 1) and the Maxwellian limit (infinite *κ*), one finds *vp* = 1 (acoustic speed), while the two coefficients become *A* = (1 − *β*)/ <sup>√</sup>*<sup>π</sup>* and *<sup>B</sup>* <sup>=</sup> 1/2, thus recovering exactly the analytical form of the original Schamel equation [28]. Figure 2 examines the influence of superthermality index *κ* and the obliqueness *θ* on the nonlinear term *A* and the dispersion term *B*. One can see that the (value of the) nonlinearity term increases, while the dispersive term decreases, if one assumes stronger deviation from the Maxwellian equilibrium, i.e., for small value of the *κ* parameter. On the other hand, for fixed *κ*, the nonlinearity term *A* attains its highest value for parallel propagation (*θ* = 0), as shown in Figure 2.

The dispersive term shows slightly more perplex behavior by increasing with growing *θ*, reaching a maximum, and then going to zero—as said above—for *θ* = *π*/2. It is evident in Equation (22) and in Figure 3a that *a*<sup>2</sup> → 0, and, hence, *A* → 0 in the limit *β* → 1; therefore, the nonlinear Equation (21) is not valid in the absence of trapped electrons. On the other hand, the dispersive term decreases with stronger magnetic field, as seen in Figure 3b.

**Figure 3.** (**a**) Nonlinearity term *A* versus *β*. *A* does not depend on the magnetic field. (**b**) Dispersion term *B* versus Ω*c*. *B* does not depend on *β*. *θ* = 10◦ is assumed and *κ* = 2, 3, 5, 10, and ∞ (top to bottom in (**a**) and bottom to top in (**b**)).

#### **4. Parametric Analysis**

In this Section, we are interested in tracing the influence of different plasma configuration parameters, such as the superthermality (spectral) index, *κ*, the electron trapping parameter, *β*, the collisional term, *ν*, the obliqueness angle, *θ*, and the ambient magnetic field (strength), *B*0, on the propagation characteristics of ion-acoustic solitary waves within the model under consideration. To see how these plasma configuration parameters affect the dynamical properties of solitary waves, first, dissipative effect is assumed to be negligible: *ν* → 0. The damped Schamel Equation (21) then reduces to a *κ*-dependent form of the Schamel-type equation [28], which possesses a solitary wave solution in the form [28,40]

$$
\psi(\xi,\tau) = \psi\_0 \text{ sech}^4 \left( \frac{\xi - \mu\_0 \tau}{\delta} \right),
\tag{25}
$$

representing a pulse-shaped excitation with amplitude *<sup>ψ</sup>*<sup>0</sup> = (15*u*0/8*A*)2, width *<sup>δ</sup>* <sup>=</sup> <sup>√</sup>16*B*/*u*<sup>0</sup> and velocity *<sup>u</sup>*<sup>0</sup> in the moving reference frame (note that the actual speed in the laboratory frame is *vp* + *u*0, so the pulse structure is superacoustic; recall that *vp* is essentially the sound speed). The product *ψ*0*δ*<sup>4</sup> = (30*B*/*A*)<sup>2</sup> is constant for a given (fixed) set of plasma parameters, that is in fact independent of *u*0. The electric field *E* (=−∇*ψ*) which is associated with the solitary potential in Equation (25) is of the form:

$$E = -E\_0 \operatorname{sech}^4\left(\frac{\tilde{\xi} - \mu\_0 \tau}{\delta}\right) \tanh\left(\frac{\tilde{\xi} - \mu\_0 \tau}{\delta}\right),\tag{26}$$

where *E*<sup>0</sup> = 225 *u*<sup>2</sup> <sup>0</sup>/(16*A*2*δ*) (this is actually the norm of the vector; the respective components are regulated by the direction cosines *lx*,*y*,*z*; recall that *l* 2 *<sup>x</sup>* + *l* 2 *<sup>y</sup>* + *l* 2 *<sup>z</sup>* = 1). The pulse form for the potential is shown in Figure 4, while the associated bipolar electric field structures are shown in Figure 5b,d.

To trace the dynamical evolution of the solitary wave solution and to elucidate the role of different plasma configuration parameters on the properties of (nonlinear) solitary waves, the nonlinear damped Schamel Equation (21) was solved numerically by using the Wolfram MATHEMATICATM software package, adopting the solitary wave solution (25) as initial condition.

**Figure 4.** (**a**) The ion-acoustic solitary potential pulse *ψ* versus the space coordinate *ξ* and the time *τ*. (**b**) Potential pulse versus *ξ* at different time instants. Here, *β* = 0.5, *κ* = 3, Ω*<sup>c</sup>* = 0.2, *ν* = 0.01, *θ* = 10◦, and *u*<sup>0</sup> = 0.01 are used. See text for details.

The time evolution of dissipative ion-acoustic solitary potential waveforms (pulses) is shown in Figure 4. The analytical solution (25) was adopted at an initial condition, and then Equation (21) was solved numerically for *ν* = 0. As expected, the pulse amplitude decreases in time due to the damping, as illustrated in Figure 4.

The influence of the trapping parameter *β* and also of the superthermality (spectral) index *κ* was investigated numerically; a snapshot at *τ* = 50 (dimensionless units) is shown in Figure 5. Here, Ω*<sup>c</sup>* = 0.2, *ν* = 0.01, *u*<sup>0</sup> = 0.01, *τ* = 50, and *θ* = 10◦ are used. One can see that both the height and the width of the solitary wave are affected by the trapping parameter *β* (Figure 5a) and by the value of *κ* (Figure 5c). As *β* increases, the waves become taller in amplitude, but the width remains unchanged; see Figure 5a. An increase in plasma superthermality (that is, a smaller value of *κ*) results in shorter and narrower solitary waves, as seen in Figure 5c. These results recover the theoretical predictions of Ref. [41] for the collision-free case, i.e., for *ν* = 0.

A similar investigation is shown in Figure 6, where the solitary wave solution (25) was obtained numerically for *κ* = 3, *β* = 0.5, *u*<sup>0</sup> = 0.01, *τ* = 50, and *θ* = 10◦. The role of the magnetic field (strength) *B*<sup>0</sup> (via Ω*c*) is shown in Figure 6a. As *B*<sup>0</sup> increases, the width of the solitary wave decreases, while the amplitude is unaffected.

Finally, in Figure 6b, various values of the collisional parameter *ν* are considered (keeping all other values fixed). As expected, the pulse amplitude decreases with higher dissipation.

**Figure 5.** (**a**) Effect of trapping parameter *β* on electrostatic solitary wave (pulse) profile and (**b**) associated electric field structures for *κ* = 3. (**c**) Effect of superthermality index *κ* on solitary pulse and (**d**) associated electric field for *β* = 0.5. Here, Ω*<sup>c</sup>* = 0.2, *ν* = 0.01, *u*<sup>0</sup> = 0.01, *τ* = 50, and *θ* = 10◦ are used.

**Figure 6.** Effect of (**a**) external magnetic field Ω*<sup>c</sup>* (for *ν* = 0.01) and of (**b**) collisionality parameter *ν* (for Ω*<sup>c</sup>* = 0.1) on electrostatic solitary waves (pulses) for *κ* = 3, *β* = 0.5, *u*<sup>0</sup> = 0.01, *τ* = 50, and *θ* = 10◦.

#### **5. Conclusions**

In this paper, the basic features of damped ion-acoustic solitary waves were investigated in the presence of trapped superthermal electrons described by a *κ*-type (non-Maxwellian) distribution [40]. The effect of ion-neutral collisions was also taken into account, leading to wave damping as expected.

The reductive perturbation approach were adopted to derive a nonlinear Schameltype partial differential equation featuring an additional damping term. The solitary wave solution of the standard (nondissipative) Schamel equation was used to solve the

damped Schamel equation numerically and to analyze the basic features of dissipative ion-acoustic solitary waves. The amplitude of solitary waves was found to decrease, while their width becomes narrower with an increase in superthermality (i.e., for a stronger deviation from Maxwellian equilibrium). The proportion of trapped electrons also affects solitary waves, since their amplitude increases in the presence of a larger proportion of trapped electrons in the plasma; on the other hand, rather counter-intuitively, their width remains the same. The above behavior was also observed via numerical integration of the dissipative Schamel equation.

While the nonlinear term is independent of the external magnetic field, the dispersive term depends strongly on the external magnetic field and in fact decreases for strong magnetic field (strength) values. Therefore, a steeper solitary wave with same maximum amplitude will be expected to occur in the presence of a stronger magnetic field, as confirmed by the numerical simulation here.

The current study focused on the 'simplest' version of a fluid model for magnetized plasma, i.e., assuming a uniform magnetic field and neglecting drift forces. A drift-kinetic approach would require the electrons to create a more complete picture to correctly account for the non-negligible *E* × *B* drift, for example. These aspects, as investigated, e.g., by Jovanoviç et al. [56], and summarized by Eliasson and Shukla [57], lie beyond the scope of the present study (weak ∼ excitations were considered here which, in addition to the absence of an ambient electric field (bias), prescribe a negligible *E* × *B* drift).

The fundamental trapping scenario was considered in this paper, (the so-called *β*trapping effect). It may occur, however, that the filamentation process in the final state of pattern formation in the electron phase space results in multiple electron transfer taking place through the separatrix; in turn, leading to additional trapping scenarios. In that case, the electric wave potential may not be expressed in a closed algebraic form, and new types of nonlinear structures may arise, as recently pointed out by Schamel [58,59]. As argued there, the existing wave theory for phase space holes, based on the linear Landau–van Kampen approach, overlooks these trapping and coherence aspects in pattern formation, and thus fails to account for a plethora of nonlinear phenomena, which are nonetheless predicted by this new approach [58,59]. Covering these aspects may form the focus of future studies.

The results, obtained here, aim to contribute to the understanding of the salient features of nonlinear electrostatic perturbations in non-Maxwellian plasmas, in account of electron trapping in phase space.

**Author Contributions:** Conceptualization, S.S.; formal analysis, S.S.; supervision, I.K.; writing original draft, S.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** Financial support by the project FSU-2021-012/8474000352 *"Modeling of Nonlinear Waves and Shocks in Space and Laboratory Plasmas"* (PI: Ioannis Kourakis), funded by Khalifa University of Science and Technology (Abu Dhabi, United Arab Emirates), is acknowledged, with thanks. Support from ADEK (Abu Dhabi Department of Education and Knowledge, United Arab Emirates), currently ASPIRE, in the form of an AARE-2018 (ADEK Award for Research Excellence 2018) grant (ADEK/HE/157/18—AARE-179) is also gratefully acknowledged.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Both authors had the pleasure and honor to have collaborated with Reinhard Schlickeiser on various projects in the past. This work is therefore dedicated to Reinhard Schlickeiser on the occasion of his 70th birthday. As a renowned plasma physicist and as an efficient supervisor, collaborator, and colleague, Reinhard Schlickeiser was a source of inspiration for both of us.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Felix Spanier 1,\*, Cedric Schreiner 2,3 and Reinhard Schlickeiser 4,5**

	- <sup>4</sup> Institut für Theoretische Physik IV, Ruhr-Universität Bochum, Universitätsstrasse 150, 44801 Bochum, Germany; rsch@tp4.rub.de
	- <sup>5</sup> Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, Leibnizstr. 15, 24118 Kiel, Germany
	- **\*** Correspondence: felix.spanier@uni-heidelberg.de

**Abstract:** Transport of energetic electrons in the heliosphere is governed by resonant interaction with plasma waves, for electrons with sub-GeV kinetic energies specifically with dispersive modes in the whistler regime. In this paper, particle-in-cell simulations of kinetic turbulence with test-particle electrons are performed. The pitch-angle diffusion coefficients of these test particles are analyzed and compared to an analytical model for left-handed and right-handed polarized wave modes.

**Keywords:** cosmic-ray transport; turbulence; plasma waves; particle-in-cell simulations

#### **1. Introduction**

The solar wind and most phases of the interstellar medium are strongly turbulent media with high magnetic Reynolds numbers of 10<sup>14</sup> [1]. Turbulence manifests itself in a spectrum of plasma waves at various length scales and frequencies. The energy distribution as a function of the frequency follows a characteristic power law. The current understanding of the turbulent processes is such that energy is injected at large scales, i.e., at small wave numbers and frequencies, and then cascades to smaller spatial scales.

The energy spectrum can be divided into several regimes, each may span several orders of magnitude in wave number or frequency. At the largest scales, the injection range is found, which then transitions into the inertial range. The inertial range can be described by magnetohydrodynamic (MHD) theory, and turbulence is dominated by the interaction of Alfvén waves. At smaller scales, kinetic effects of the particles come into play.

This high wave number regime is often referred to as the kinetic, dispersive, or dissipation range of the spectrum, since the waves become dispersive and dissipation starts to set in. While the spectrum extends to even smaller scales, damping eventually becomes dominant and leads to an exponential cutoff of the energy spectrum.

Power-law distributions of the fluctuating magnetic energy are expected in the injection, inertial, and dissipation range of the spectrum. However, the spectral indices may differ among the individual regimes. Goldreich and Sridhar [2,3] presented a detailed model of the turbulent energy cascade in the inertial regime. Their model predicts a spectral index of −5/3, which actually seems to be realized in the solar wind [4]. Subsequent models by Galtier et al. [5,6] give rise to a spectral index *k*−<sup>2</sup> ⊥ .

Kinetic turbulence in the dissipation range is an active field of research [7]. In particular, the composition of the wave spectrum is subject to discussion, because a transition from non-dispersive Alfvén waves to dispersive wave modes is expected. Possible candidates for the waves in the dissipation range are the so-called kinetic Alfvén waves (KAWs) and whistler waves [8].

**Citation:** Spanier, F.; Schreiner, C.; Schlickeiser, R. Determining Pitch-Angle Diffusion Coefficients for Electrons in Whistler Turbulence. *Physics* **2022**, *4*, 80–103. https:// doi.org/10.3390/physics4010008

Received: 30 August 2021 Accepted: 7 December 2021 Published: 20 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The impact of kinetic turbulence on the transport of energetic particles is another major topic. The transport of energetic protons is well described by models of Alfvénic turbulence, since the protons mainly interact with these waves at low frequencies. The theoretical framework of quasi-linear theory (QLT) can be used to describe particle transport by a series of resonant interactions with the magnetic fields of Alfvén waves, which leads to scattering of the particles [9–11]. This theory describes changes of the particles' pitch angles (the angle of the velocity vector relative to a background magnetic field), momenta, or positions as diffusion processes and allows to predict diffusion coefficients and other quantities, such as the mean free path, which can be compared to observations.

Dispersive waves are more difficult to handle in (analytical) theory. Nonetheless, QLT can also yield predictions for particle transport in a medium containing dispersive waves [12,13]. The introduction of dispersive waves can even solve some of the problems that are encountered if a purely Alfvénic spectrum of waves is assumed [14]. Still, the model remains an approximation, and computer simulations are used to clarify some of the details that are not included in the analytical theory. Different kinds of simulations are used to study different physical regimes and processes—from the acceleration of particles [15,16] to the transport of energetic particles—considering both non-dispersive [17] and dispersive waves [18,19].

The key problem that has been chosen for the subject of this study is the process of electron transport in dispersive turbulence. The transport of electrons at sub-GeV energies has been of high interest for quite some time [20]. As mentioned above, particle acceleration in Alfvénic turbulence in the inertial range is well understood. However, turbulence on kinetic scales still poses problems for both observations and modeling.

#### **2. Theory**

#### *2.1. Turbulence Theory*

From the observations, it is not entirely clear which types of plasma waves constitute the spectrum of turbulent waves in the dispersive and dissipative regime. KAWs and whistler waves are both possible candidates [8]. While KAWs represent the continuation of the Alfvén mode (Equation (13) below) in the dispersive regime with very large perpendicular wave numbers (perpendicular wavelength in the range of proton gyroradius), whistler waves (Equation (16) below) are right-handed polarized modes at large wave numbers. It is reasonable to assume that non-dispersive Alfvén waves simply transition to KAWs. However, observations reveal that whistler waves are also present in various regions of the heliosphere, such as in the interplanetary medium [21], close to interplanetary shocks [22,23] or planetary bow shocks [24], and also in the Earth's ionosphere and foreshock region [25,26].

Whereas left-handed polarized Alfvén waves are damped by protons and cannot cascade to frequencies above the proton cyclotron frequency, a spectrum of whistler waves may extend to frequencies beyond the ion-cyclotron frequencies. Whistler waves primarily interact with electrons and are also damped by electrons at higher frequencies (close to the electron cyclotron frequency). This is an interesting aspect of kinetic turbulence, since a population of whistler waves can heat the electrons in the solar wind or even accelerate particles in the high-energy tail of the thermal spectrum.

Kinetic turbulence includes the smallest length scales, where the interaction of waves and particles becomes important. Although the wave-particle interactions are not explicitly included in the theory, their effect has to be considered by allowing dispersive waves and damping. This regime is generally more complicated and less understood than MHD turbulence.

The general picture associated with turbulence is as follows. Energy is injected into the system at a large outer scale, small wave numbers *k*. The energy is transported via the interaction of waves to smaller spatial scales (larger wave numbers), and the (magnetic) energy spectrum, *EB*(*k*), follows a power-law distribution. This is the inertial range. The spectrum steepens as the kinetic regime or dissipation range is reached, but energy is still transported to smaller scales. First, only ion effects will start influencing the plasma dynamics, but at even larger wave numbers, the electrons can also interact with the plasma waves. This is where the energy spectrum is cut off. One aspect that has been discussed in greater detail since the seminal paper of Goldreich and Sridhar [3] is the possible anisotropy with respect to the background magnetic field: the turbulent spectrum may behave differently depending on wave numbers parallel (*k*) or perpendicular (*k*⊥) to the background magnetic field.

The special case of whistler turbulence has been discussed in greater detail. The properties of this whistler turbulence have been analyzed using kinetic simulations in two [27,28] and three [29–32] dimensions. These studies suggest a steeper energy spectrum than for Alfvénic turbulence, with a spectral index *s* in the range between −3.7 and −5.5 and a possible break in the energy spectrum [32]. Results by Chang [30] also suggest that the wave vector anisotropy depends on the choice of the plasma beta. A relatively isotropic spectrum is obtained for a plasma beta *β* ∼ 1, whereas *β* < 1 yields an anisotropic cascade which favors the transport of energy to larger *k*⊥. The plasma beta is the ratio of thermal to magnetic energy. The anisotropy additionally depends on the energy deployed to the electromagnetic fields of the turbulent whistler waves [29].

#### *2.2. Subluminal Parallel Waves in Cold Plasmas*

In warm thermal plasmas with low plasma betas, the real part of the dispersion relation agrees well with the cold plasma dispersion relation, so, the latter is used here. In addition, the resonance broadening effects, caused by a finite imaginary part of the dispersion relation in warm plasmas implying a finite weak-damping rate; those effects were considered by Schlickeiser and Achatz [33].

Using the convention of frequencies *ω* > 0, where *ω* is the real part of the generally complex frequency, but *<sup>k</sup>* ∈ [−∞, <sup>∞</sup>] (here, the case *<sup>k</sup>*<sup>⊥</sup> = 0, also known as slab case, is treated), the dispersion relation of right-handed ("R") and left-handed ("L") polarized undamped low-frequency (*ω* ≤ Ω*e*,0, with Ω*e*,0 being the electron gyrofrequency) parallel Alfvén wave in a cold electron-proton background plasma reads [34]:

$$m\_L^2 = 1 - \frac{\omega\_{pi}^2}{\omega(\omega - \Omega\_i)} - \frac{\omega\_{pe}^2}{\omega(\omega + \Omega\_\varepsilon)},\tag{1}$$

$$m\_R^2 = 1 - \frac{\omega\_{pi}^2}{\omega(\omega + \Omega\_i)} - \frac{\omega\_{pe}^2}{\omega(\omega - \Omega\_\ell)}\tag{2}$$

$$\frac{k\_{\parallel}^2 c^2}{\omega\_{L,\text{R}}^2} = 1 - \frac{\omega\_{pi}^2}{\omega(\omega \mp \Omega\_i)} - \frac{\omega\_{pe}^2}{\omega(\omega \pm \Omega\_e)},\tag{3}$$

$$\frac{k\_{\parallel}^2 c^2}{\omega\_{L,R}^2} - 1 = -\frac{c^2 \Omega\_i^2}{v\_A^2} \frac{M+1}{(\omega \mp \Omega\_i)(\omega \mp M\Omega\_i)}\tag{4}$$

with the proton-to-electron mass ratio, *M* = *mp*/*me* = 1836, and the Alfvén speed, *vA* <sup>=</sup> *<sup>β</sup>Ac* <sup>=</sup> 2.18 <sup>×</sup> <sup>10</sup>11*B*[G]*n*−1/2 *<sup>i</sup>* [cm−3] Here, *<sup>ω</sup>pi* is the ion plasma frequency, *<sup>ω</sup>pe* is the electron plasma frequency, Ω*<sup>i</sup>* is the ion gyrofrequency, Ω*<sup>e</sup>* is the electron gyrofrequency, *B* is the magnetic field, *β* is the plasma beta, and *c* is the speed of light. For subluminal wave phase speeds,

$$\left| V\_{\text{phase}} \right| = \left| \frac{\omega\_{L,\mathbb{R}}}{k\_{\parallel}} \right| \ll \mathcal{c} \,\,\,\tag{5}$$

which has to be checked a posteriori, the dispersion relation (4) simplifies to

$$\frac{k\_{\parallel}^2 c^2}{\omega\_{L,R}^2} \simeq -\frac{(M+1)\Omega\_i^2 \omega\_{L,R}^2}{v\_A^2 (\omega \mp \Omega\_i)(\omega \mp M\Omega\_i)}.\tag{6}$$

It is convenient to introduce dimensionless frequencies and parallel wave numbers by defining

$$y\_{L,R} = \frac{\omega\_{L,R}}{\Omega\_i} > 0 \quad \text{and} \quad k = \frac{k\_{\parallel}}{k\_c},\tag{7}$$

respectively, with

$$k\_{\varepsilon} = \frac{\Omega\_i}{v\_A} \,\, \,\, \tag{8}$$

so that the subluminal dispersion relation becomes

$$k^2 = -\frac{(M+1)y\_{L,R}^2}{(y\_{L,R} \mp 1)(y\_{L,R} \pm M)}\tag{9}$$

with the two solutions,

$$k\_{1,2} = -\frac{\sqrt{(M+1)}y\_{L,R}}{\sqrt{(y\_{L,R}\mp 1)(y\_{L,R}\pm M)}}.\tag{10}$$

As *yL*,*<sup>R</sup>* is always positive, the solution *k*<sup>1</sup> > 0 describes forward-moving waves with positive phase speed, whereas the negative solution *k*<sup>2</sup> = −*k*<sup>1</sup> < 0 describes backwardmoving waves.

#### 2.2.1. Left-Handed Modes

Equation (9) indicates that no left-handed polarized solution with *yL* > 1 exists, so, a further simplification of Equation (9) for left-handed polarized waves is possible, assuming that *M* 1:

$$k^2 \simeq \frac{y\_L^2}{1 - y\_L} \tag{11}$$

with the solutions,

$$y\_L(k) = \frac{k^2}{2} \pm \sqrt{\frac{k}{2}}\sqrt{k^2 - 4} \simeq \begin{cases} |k| & \text{for } k \ll 1, \\ 1 - 1/k^2 & \text{for } k \gg 1, \end{cases} \tag{12}$$

corresponding to

$$
\omega\_L \cong \begin{cases}
V\_A |k\_{\parallel}| & \text{for } k\_{\parallel} \ll k\_c, \\
\Omega\_i \left(1 - k\_c^2 / k\_{\parallel}^2\right) & \text{for } k\_{\parallel} \gg k\_c.
\end{cases} \tag{13}
$$

#### 2.2.2. Right-Handed Modes

The right-handed solutions of Equation (9),

$$k^2 = -\frac{(M+1)y\_R^2}{(y\_R+1)(y\_R+M)'} \tag{14}$$

can be approximated under the assumption that *M* 1:

$$y\_R = \frac{(M+1)k^2}{2(1+k^2+M)} \left(1 - \frac{2}{M+1} + \sqrt{1 + \frac{4M}{k^2(M+1)}}\right). \tag{15}$$

Depending on *k*, different regimes can be identified:

$$\omega\_R = \begin{cases} V\_A |k\_{\parallel}| & \text{for} \quad |k\_{\parallel}| < k\_c \\ \Omega\_i + \frac{V\_A^2 k\_{\parallel}^2}{\Omega\_i} & \text{for} \quad k\_c \le |k\_{\parallel}| \le M^{1/2} k\_c \\ \Omega\_c \left(1 - \frac{M k\_c^2}{k\_{\parallel}^2}\right) & \text{for} \quad |k\_{\parallel}| > M^{1/2} k\_c \end{cases} \tag{16}$$

The first range describes the linear dispersion regime, the second range is the whistler regime, and the last range is the electron–cyclotron range. While these approximate solutions provide good estimates to the real solution, there is a major problem: the solutions do not provide continuous coverage. An alternative approximation is

$$y\_{\mathcal{R}}(k) \simeq |k|(1+|k|). \tag{17}$$

In the following, the particle scattering by parallel waves at electron or ion cyclotron frequencies are ignored as soon as these are highly damped in a realistic warm thermal plasma, so that the resonant interaction does not apply; see [33] for a discussion of waveparticle interactions with damped waves. For left-handed and right-handed polarized waves, this restricts the normalized wave numbers to values of *k* ≤ 1 and *k* ≤ *M*, respectively.

#### *2.3. Particle Transport*

Any charged particle of given velocity *<sup>v</sup>*, Lorentz factor *<sup>γ</sup>* = (<sup>1</sup> <sup>−</sup> (*v*2/*c*2))−1/2, pitchangle cosine *<sup>μ</sup>* = *<sup>v</sup>*/*v*, mass *<sup>m</sup>*, charge *qi* = *<sup>e</sup>*|*Zi*|*<sup>Q</sup>* with *<sup>Q</sup>* = sgn(*Zi*), *Zi* being the ion charge number, and relativistic gyrofrequency Ω = *Q*Ω/*γ* with positive Ω = |*q*|*B*0/(*mc*) with Ω being the particle's non-relativistic gyrofrequency, *q* being the electric charge, and *B*<sup>0</sup> being the magnetic background field, interacts with parallel right-handed and lefthanded polarized plasma waves whose wave number, *k*, and real frequency, *ωR*,*L*, obey the resonance condition,

$$
v\mu k\_{\parallel} - \omega\_{\mathbb{R},L}(k\_{\parallel}) \mp \frac{Q\Omega}{\gamma} = 0.\tag{18}$$

Introducing

$$\alpha = \frac{p}{m\_c \varepsilon}, \quad \epsilon = \frac{v\_A}{v} = \beta\_A \frac{(1+\varkappa^2)^{1/2}}{\varkappa}, \tag{19}$$

with the particle momentum *p* and the dimensionless frequency and wave number (7), the resonance condition (18) reads:

$$
\Omega \left( \frac{\mu k}{\epsilon} - y\_{R,L}(k) \mp S\_{\dot{l}} \right) = 0 \tag{20}
$$

with

$$S\_i = \frac{\mathbb{Q}|Z\_i|m\_p}{m\gamma} = \frac{1}{\gamma} \begin{cases} 1 & \text{for protons,} \\ -M & \text{for electrons.} \end{cases} \tag{21}$$

The quasilinear Fokker–Planck coefficient for the pitch-angle cosine, *μ*, is given by

$$\begin{split} D\_{\mu\mu}(\mu) &= \quad \frac{\pi^2 \Omega^2 (1 - \mu^2)}{B\_0^2} \int\_{-\infty}^{\infty} \mathbf{d}k\_{\parallel} \\ &\times \quad \left[ \mathcal{g}\_R(k\_{\parallel}) \, \delta(\upsilon k\_{\parallel} - \omega\_R - \Omega') \left( 1 - \frac{\mu \omega\_R}{k\_{\parallel} \upsilon} \right)^2 \\ &+ \quad \mathcal{g}\_L(k\_{\parallel}) \, \delta(\upsilon k\_{\parallel} - \omega\_L + \Omega') \left( 1 - \frac{\mu \omega\_L}{k\_{\parallel} \upsilon} \right)^2 \right] \end{split} \tag{22}$$

with the magnetic fluctuation spectra of right/left-handed polarized waves *gR*,*L*(*<sup>k</sup>*), where the total magnetic field fluctuations are determined as in [35]:

$$\mathbb{E}(\delta B)^2 = 2\pi \int\_{-\infty}^{\infty} \mathrm{d}k\_{\parallel} [\mathcal{g}\_L(k\_{\parallel}) + \mathcal{g}\_R(k\_{\parallel})].\tag{23}$$

In Equation (22), the frequencies *<sup>ω</sup>R*,*L*(*<sup>k</sup>*) are determined by the solutions of the dispersion relations, discussed in Section 2.2 above. The function *δ*(*x*) is the Dirac delta function.

In terms of the normalized wave number, *<sup>k</sup>* = *kck*, and frequency, *<sup>ω</sup>R*,*<sup>L</sup>* = <sup>Ω</sup>*iyR*,*L*, the Fokker–Planck coefficient (22) reads:

$$\begin{split} D\_{\mu\mu}(\mu) &= \quad \frac{\pi^2 \Omega^2 k\_c (1 - \mu^2)}{\Omega\_i^2 B\_0^2} \int\_{-\infty}^{\infty} \mathbf{d}k\_{\parallel} \\ &\times \quad \left[ \lg\_R(k\_{\parallel}) \,\delta(k\mu/\epsilon - y\_R(k) - S\_i) \left( 1 - \frac{\mu \omega\_R}{k\_{\parallel} v} \right)^2 \\ &+ \quad \lg\_L(k\_{\parallel}) \,\delta(k\mu/\epsilon - y\_L + S\_i) \left( 1 - \frac{\mu \omega\_L}{k\_{\parallel} v} \right)^2 \right]. \end{split} \tag{24}$$

The calculation of the Fokker–Planck coefficient requires a knowledge of the magnetic fluctuation spectra. The correct theoretical description is complicated, as described in Section 2.1 above, but results from numerical calculations can be inferred.

Deriving the Fokker–Planck coefficients in general is quite an involved task but one can derive some limiting cases. It is helpful to account for the relative abundance of forwardmoving and backward-moving waves and the corresponding polarization states. Let us introduce the cross helicities, *HL*,*<sup>R</sup>* ∈ [−1, 1] for left/right-handed polarized waves to express the spectra (23) of backward-propagating ("b") and forward-propagating ("f") waves:

$$\mathcal{g}\_{L,R}^{b} = \frac{1 - H\_{L,R}}{2} \mathcal{g}\_{L,R}(k) \Theta(-k),\tag{25}$$

$$\mathbf{g}\_{L,\mathbb{R}}^f = \frac{1 + H\_{L,\mathbb{R}}}{2} \mathbf{g}\_{L,\mathbb{R}}(k) \, \Theta(k). \tag{26}$$

The step function Θ(±*k*) appears because backward-moving and forward-moving waves only occur at negative and positive wave numbers, respectively. In general, these cross helicities can depend on the wave number, but throughout this article isospectral turbulence is adopted, where *HL*,*<sup>R</sup>* are constants (independent of *k*). The magnetic helicity *σ*(*k*) ∈ [−1, 1] characterizes the relative abundances of left-handed and right-handed polarized waves:

$$g\_L(k) = \frac{1 + \sigma(k)}{2} g\_{\text{tot}}(k),\tag{27}$$

$$g\_R(k) = \frac{1 - \sigma(k)}{2} g\_{\text{tot}}(k) \tag{28}$$

where *g*tot(*k*) is the total wave abundance at a specific wave number. For parallel plasma waves, *σ*(*y* > 1) = −1 as soon as no left-handed polarized waves exist at these normalized frequencies.

Using the helicities introduced, the Fokker–Planck coefficient (24) reads:

$$\begin{split} D\_{\mu\mu}(\mu) &= \quad \frac{\pi^2 \Omega^2 k\_c (1 - \mu^2)}{\Omega\_i^2 B\_0^2} \int\_{-\infty}^{\infty} \mathrm{d}k\_{\parallel} g\_{\text{tot}}(k) \\ &\quad \times \quad \left[ (1 - H\_R)(1 - \sigma(k)) \, \delta(k\mu/\epsilon + y\_R + S\_i) \left( 1 + \frac{\epsilon \mu y\_R(k)}{k} \right)^2 \right. \\ &\quad + \quad \left( 1 - H\_L \right)(1 + \sigma(k)) \, \delta(k\mu/\epsilon + y\_R(k) - S\_i) \left( 1 + \frac{\epsilon \mu y\_L(k)}{k} \right)^2 \\ &\quad + \quad \left( 1 + H\_R \right)(1 + \sigma(k)) \, \delta(-k\mu/\epsilon + y\_R + S\_i) \left( 1 - \frac{\epsilon \mu y\_R(k)}{k} \right)^2 \\ &\quad + \quad \left( 1 + H\_L \right)(1 - \sigma(k)) \, \delta(-k\mu/\epsilon + y\_R(k) - S\_i) \left( 1 - \frac{\epsilon \mu y\_L(k)}{k} \right)^2 \right]. \end{split} \tag{29}$$

#### 2.3.1. Interactions in the Whistler Regime

For frequencies above the ion cyclotron frequency, only right-handed waves exist obeying the dispersive whistler mode dispersion relation. As discussed above, the turbulent spectrum typically has a much softer spectral index than 2 (theoretical values are in the range of 3 to 6 [7,36]) in this case.

We consider the case *HR* = −1 and *σ*(*k*) = −1, only backward-moving right-handed polarized waves, which reduces the Fokker–Planck coefficient (29) to

$$D\_{\mu\mu}(\mu) = \begin{array}{c} \frac{\pi^2 \Omega^2 k\_c (1 - \mu^2)}{\Omega\_i^2 B\_0^2} \int\_{k\_0}^{\infty} \mathrm{d}k\_{\parallel} \mathcal{g}\_{\text{tot}}(k) \\\\ \delta(k\mu/\epsilon + y\_R + \mathcal{S}\_i) \left(1 + \frac{\epsilon \mu y\_R(k)}{k}\right)^2 \end{array} \tag{30}$$

The result for forward-moving waves is similar:

$$D\_{\mu\mu}(\mu) = \frac{\pi^2 \Omega^2 k\_c (1 - \mu^2)}{\Omega\_i^2 B\_0^2} \int\_{k\_0}^{\infty} \mathbf{d}k\_{\parallel} \mathbf{g}\_{\text{tot}}(k) \tag{31}$$

$$\times \quad \left[ \delta(-k\mu/\epsilon + \mathcal{y}\_R + S\_i) \left( 1 - \frac{\epsilon \mu \mathcal{y}\_R(k)}{k} \right)^2 \right].$$

With Equation (17), the resonance condition with positive values of *k* reads:

$$0 = S\_l + k \frac{\mu}{\epsilon} + \begin{cases} |k| & \text{for } k \le 1, \\ k^2 & \text{for } 1 \le k \le M^{1/2} = 43. \end{cases} \tag{32}$$

It can be shown that this equation for protons and electrons has only one solution *kr* > 0. In the Alfvénic wave number range (*k* ≤ 1), this is trivial: *kr* = −*Si*/(1 + (*μ*/)), which can be positive depending on the signs of *Si* and *μ*.

In the whistler wave number range (0 ≤ *k* ≤ 43), the proof is a bit more involved. Here, Equation (32) has two solutions:

$$k\_1 = \frac{1}{2} \left( \sqrt{\frac{\mu^2}{\epsilon^2} - 4S\_{\bar{i}}} - \frac{\mu}{\epsilon} \right),\tag{33}$$

$$k\_2 = -\frac{1}{2} \left( \sqrt{\frac{\mu^2}{\epsilon^2} - 4S\_i} + \frac{\mu}{\epsilon} \right). \tag{34}$$

To obtain real-valued solutions (33) and (34), the condition *<sup>μ</sup>*<sup>2</sup> <sup>≥</sup> <sup>4</sup>*Si*<sup>2</sup> must be fulfilled. Assuming that this is fulfilled, one notes that for non-negative values of *μ* ≥ 0, the solution *k*2(*μ* ≥ 0) < 0 is always negative, leaving only one solution for *kr* = *k*1(*μ* ≥ 0) > 0. Alternatively, for negative values of *μ* = −|*μ*| the solutions (33) and (34) become

$$k\_1(\mu < 0) = \frac{1}{2} \left( \sqrt{\frac{\mu^2}{\epsilon^2} - 4S\_i} - \frac{|\mu|}{\epsilon} \right),\tag{35}$$

$$k\_2(\mu < 0) = \frac{1}{2} \left( \frac{|\mu|}{\epsilon} - \sqrt{\frac{\mu^2}{\epsilon^2} - 4S\_i} \right). \tag{36}$$

To further evaluate the solution, it is necessary to distinguish between positive (for protons) and negative values (for electrons) of *Si*. For electrons the solution *k*2(*μ* < 0) < 0 is again always negative. For protons, both solutions (35) and (36) are positive, but the second one is

$$k\_2(\mu < 0, S\_i > 0) \le S\_i^{1/2} = \frac{1}{2\gamma} < 1\tag{37}$$

always (for protons *Si* = 1/ √2*γ*). This solution is positive but is in contradiction with the above-made assumption of the modes in the whistler regime with *k* > 1. This leaves us with only one solution for the resonant wave number *kr* = *k*1(*μ* < 0, *Si* > 0) in the whistler wave number range.

One then obtains:

$$D\_{\mu\mu}(\mu) \quad = \quad \frac{\pi^2 \Omega^2 k\_c (1 - \mu^2)}{\Omega\_i^2 B\_0^2} \Theta(k\_r - k\_0) \, \Theta(M^{1/2} - k\_r)$$

$$\times \quad \frac{\mathcal{g}\_{\text{tot}}(k\_r)}{\left| \frac{\text{d} \eta\_R}{\text{d} \text{k}} + \frac{\mu}{\text{d}} \right|\_{k = k\_r}} \left( 1 + \frac{\epsilon \mu \mathcal{y}\_R(k\_r)}{k\_r} \right)^2. \tag{38}$$

The case of forward-moving right-handed polarized waves is similar to the backwardmoving ones. The main difference is the resonant wave number,

$$k\_r = \frac{1}{2} \left( \sqrt{\frac{\mu^2}{\epsilon^2} - 4S\_i} + \frac{\mu}{\epsilon} \right). \tag{39}$$

#### 2.3.2. Alfvén and Whistler Contributions

The total Fokker–Planck coefficient can be written as a sum of Alfvén and whistler contributions:

$$D\_{\mu\mu}(\mu) = D\_{\mu\mu}^{\mathcal{A}} + D\_{\mu\mu}^{\mathcal{W}}.\tag{40}$$

For the Alfvén wave Fokker–Planck coefficient, inserting the asymptotic expansions *yR*,*L*(*k* ≤ 1) *k* of Equations (12) and (15) one obtains:

$$\begin{split} D^{\rm A}\_{\mu\mu}(\mu) &= \quad \frac{\pi^2 \Omega^2 k\_c (1 - \mu^2)}{\Omega^2\_l B\_0^2} \int\_{k\_0}^1 \mathrm{dk}\_{\parallel} g\_{\rm tot}(k) \\ &\times \quad \left\{ (1 + \epsilon \mu)^2 \left[ (1 - H\_R)(1 - \sigma) \, \delta(k(1 + \mu/\epsilon) + S\_i) \right. \right. \\ &\left. + \quad (1 - H\_L)(1 + \sigma) \, \delta(k(1 + \mu/\epsilon) - S\_i) \right] \\ &\left. + \quad (1 - \epsilon \mu)^2 \left[ (1 + H\_R)(1 - \sigma) \, \delta(k(1 - \mu/\epsilon) + S\_i) \right. \\ &\left. + \quad (1 + H\_L)(1 + \sigma) \, \delta(k(1 - \mu/\epsilon) - S\_i) \right] \right\}. \end{split} \tag{41}$$

The whistler contribution is calculated above, while for completeness is given here in the same form:

$$\begin{split} D^{\mathsf{W}}\_{\mu\mu}(\mu) &= \ \frac{\pi^2 \Omega^2 k\_c (1 - \mu^2)}{\Omega^2\_i B\_0^2} \int\_1^{M^{1/2}} \mathrm{d}k \, \vert \, \mathsf{g}\_{\mathsf{tot}}(k) \\ &\times \ \left[ (1 + \epsilon \mu)^2 (1 - H\_R)(1 - \sigma) \, \delta(k^2 + k\mu/\epsilon + S\_i) \\ &+ \ (1 - \epsilon \mu)^2 (1 + H\_R)(1 - \sigma) \, \delta(k^2 - k\mu/\epsilon + S\_i) \right]. \end{split} \tag{42}$$

#### 2.3.3. Electrons

The Fokker–Planck coefficient derived in the previous section hold for any particle. In general, the integration is not complicated as the delta function of the resonance condition helps to simplify the calculations. However, a specific turbulent spectrum has to be defined. We refrain from performing the integral here but point out which interaction will take place. There is a clear difference between electrons and protons, and the discussion in this Section is limited to the electron case.

For Alfvén waves, one can distinguish the interaction of electrons with forward/backwardmoving left/right-handed modes. Defining

$$
\mu\_{\mathcal{R}}(\mathbf{x}) = \frac{\mathcal{B}\_A(M - \sqrt{1 + \mathbf{x}^2})}{\mathbf{x}},\tag{43}
$$

$$
\mu\_L(\mathbf{x}) = \frac{\beta\_A(M + \sqrt{1 + \mathbf{x}^2})}{\mathbf{x}},\tag{44}
$$

one can constrain the waves for which the resonant interaction with electrons is possible:


$$
\mu\_0(\mathbf{x}) = \frac{\beta\_A M^{1/2} (\sqrt{1 + \mathbf{x}^2} - 1)}{\mathbf{x}} \tag{45}
$$

is defined. The resonant interaction takes place within the following range:

$$-\mu\_0(\mathbf{x}) \le \mu \le \mu\_R(\mathbf{x}),\tag{46}$$

The resulting total Fokker–Planck scattering coefficient is shown in Figure 1.

**Figure 1.** The pitch-angle dependence of the Fokker–Planck scattering coefficient model calculation for electrons at extreme ends. A power-law spectrum with with spectral index *s* = 3 in the wave number range 0.01 < *k* < *M*1/2, where *M* is the proton-to-electron mass ratio, is assumed. Electrons at <sup>80</sup> keV and 3MeV are considered. The proton gyrofrequency <sup>Ω</sup>*<sup>p</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> 104 Hz, and the Alfvén speed *vA* = 0.001*c* in fractions of the speed of light *c*.

#### **3. Numerical Methods**

#### *3.1. Particle-in-Cell Simulations*

To be able to model dispersive waves of different kinds and to obtain a self-consistent description of electromagnetic fields and charged particles in the plasma, a fully kinetic particle-in-cell (PiC) approach [37] is employed here. In particular, the explicit second-order PiC code ACRONYM [38] is used which is fully relativistic, parallelized and three-dimensional. Although the PiC method might not be the most efficient numerical technique when dealing with proton effects, we still favor this approach because of its versatility. A more detailed discussion of advantages and drawbacks as well as a direct comparison of PiC and MHD approaches with the specific problem of the interaction of protons and left-handed polarized waves can be found in Sections 3 and 6 of Ref. [39]. However, the PiC approach is well-suited for the study of electron scattering as soon as the time and length scales of electron interactions

are closer to the scales of time step lengths and cell sizes in PiC simulations, thus reducing computing time compared to simulations, in which proton interactions are studied.

The details of the PiC code are not discussed here. The simulation technique used here does not differ from standard techniques. Two points, however, to be mentioned: the initialization of turbulence and the tracking of test particles. Turbulence is discussed in Section 3.1.1 just below, because the numerical treatment is inherently connected to the physical processes of turbulence. The treatment of energetic particles is divided into to two parts: the injection of particles and the analysis of the particle data.

#### 3.1.1. Setup of Turbulence Simulations

Here, a simulation setup that was inspired by Gary et al. [27] is used. Waves are initialized at low *k*, *k* < Ω*i*/*VA*. The layout of the initial waves in the wave number space is explained below and is drawn in Figure 2 for the two-dimensional setup. However, in the PiC simulation, the velocity space and the electromagnetic fields are considered threedimensional. As shown below, the two-dimensional simulations give an energy cascade similar to that of the three-dimensional simulations. As highlighted by Gary et al. [29], the three-dimensional case differs mainly in the anisotropy and a break at *k*⊥*c*/Ω*e*. The question whether particle transport is different in two and three dimensions may be understood from the fact that particle motion is still three-dimensional. The theoretical description is assuming gyrotropy anyway, so the different perpendicular directions are therefore averaged out.

In the simulations, the natural mass ratio, *mp*/*me* = 1836, is used. Waves are initialized with equal amplitudes and a random phase angle. The total magnetic energy density of the initial waves is set to 10% of the energy density of the background magnetic field. This can be expressed by *δB*2/*B*<sup>2</sup> <sup>0</sup> = 0.1, where *<sup>δ</sup>B*<sup>2</sup> = <sup>∑</sup>*<sup>j</sup> <sup>δ</sup>B*<sup>2</sup> *<sup>j</sup>* and *j* denotes individual waves.

To analyze the simulations, the spectra of the magnetic energy density, *EB* <sup>=</sup> <sup>|</sup> *B*2( *k*)|/(<sup>8</sup> *<sup>π</sup>*), in wave number space are considered. A two-dimensional energy spectrum, *EB*(*<sup>k</sup>*, *<sup>k</sup>*⊥), can be obtained by Fourier transforming the field data to obtain the perpendicular coordinate *k*⊥. The parallel direction is equivalent to the *z*-direction of the simulation, whereas the perpendicular direction is represented by the *x*-direction in a two-dimensional simulation or by the *x*-*y*-plane in a three-dimensional simulation. A one-dimensional spectrum *EB*(*k*) can be obtained by integrating over the angle *<sup>θ</sup>* in the *<sup>k</sup>*-*k*⊥-plane. Additional onedimensional spectra, *EB*(*<sup>k</sup>*) and *EB*(*k*⊥), are obtained by integrating *EB*(*<sup>k</sup>*, *<sup>k</sup>*⊥) over *<sup>k</sup>*<sup>⊥</sup> and *k*, respectively.

Electron transport is studied in two simulations, S1 and S2. The aim is to resolve several wave numbers in both the undamped and damped regimes of the whistler mode. This should allow us to see differences in the spectral slope or the anisotropy in both regimes. To resolve the relatively large spatial scales of the undamped regime, large simulation boxes are required. Thus, the decision was taken to restrict the investigations to two-dimensional setups.

Simulations S1 and S2 are characterized by the physical and numerical parameters listed in Tables 1 and 2, respectively.

The setups are aimed to simulate decaying turbulence with a set of 42 initially excited whistler waves according to Figure 2.

**Table 1.** Physical parameters for simulations S1 and S2: electron plasma frequency, *ωp*,*e* electron cyclotron frequency, Ω*e*, and thermal speed, *v*th,*<sup>e</sup>* of electrons, sum *δB*<sup>2</sup> of the squares of the magnetic field amplitudes of the individual waves, and plasma beta *β*.



S1 2048 2048 1.0 <sup>×</sup> 105 7.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 4.1 <sup>×</sup> <sup>10</sup>−<sup>2</sup> <sup>256</sup> S2 2048 2048 1.0 <sup>×</sup> 105 3.5 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 2.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> <sup>256</sup>

**Table 2.** Numerical parameters for the two-dimensional simulations S1 and S2: number of cells, *<sup>N</sup>* and *N*⊥, in the directions parallel and perpendicular to the background magnetic field, respectively;


**Figure 2.** Schematic representation of two-dimensional wave number space. The discretized wave vectors are represented by the gray boxes, and the axes mark the directions parallel (*kz*) and perpendicular (*kx*) to the background magnetic field *B*<sup>0</sup> in the case of a two-dimensional simulation. For the simulation of decaying turbulence in two dimensions, a set of 42 initial waves is excited, where each wave occupies one position on the grid. These positions are indicated by the blue boxes in accordance with the setup specified by Gary [27].

#### 3.1.2. Turbulence Spectra

Here, the simulations S1 and S2 with the parameters from Tables 1 and 2 are briefly discussed.

Figure 3 shows the perpendicular spectra *EB*(*k*⊥) of the magnetic field energy from simulations S1 (Figure 3a) and S2 (Figure 3b). At small, perpendicular wave numbers, the magnetic energy distribution follows a power-law with spectral index *<sup>s</sup>*<sup>⊥</sup> = −3.1 and −3.0 for S1 and S2, respectively. After the break, the spectra steepen, and significant differences between both simulations become obvious in the different spectral indices.

The numerical noise level in simulation S2 is about one order of magnitude lower than in S1, which allows an energy cascade to higher wave numbers. This can be explained by the lower plasma temperature in S2, leading to less kinetic energy of the particles and therefore less fluctuations in the electromagnetic fields. The flatter spectrum in S2 (after the break; compared to S1) agrees with results from Chang [30], who reported a more efficient perpendicular energy transport with decreasing plasma beta *β*.

Chang [30] also observed stronger anisotropy in simulations with lower *β*. However, this is not supported by the data from simulations S1 and S2. The parallel spectra *EB*(*<sup>k</sup>*) are depicted in Figure 4. In both cases, the parallel spectra do not contain a break and are steeper than the perpendicular spectra at small wave numbers. For S1, the parallel spectrum reaches the noise level approximately at the position where cyclotron damping is

assumed to set in. Figure 4a,b, however, shows that the parallel spectrum in simulation S2 extends to wave numbers in the damped regime. The slope does not change at the transition into the dissipation range and is flatter than the slope in the perpendicular spectrum at corresponding *k*⊥. Thus, the parallel energy transport is assumed to dominate at large wave numbers. Unfortunately, the turbulent cascade reaches the numerical noise level prior to that.

**Figure 3.** Normalized magnetic field energy distribution *EB*(*k*⊥)/*EB*<sup>0</sup> over the perpendicular wave number, *k*<sup>⊥</sup> normalized to *c*/*ωp*, for simulations S1 (**a**) and S2 (**b**). Here, *EB*<sup>0</sup> is the magnetic field energy of the background, *ωp* is the plasma frequency, and *c* is the speed of light. The data are obtained at four times *t*|Ω*e*| as indicated, where Ω*<sup>e</sup>* is the electron gyrofrquency. At the earliest time steps, the spectra reach their steady states. Later in the simulations, the shapes of the spectra do not change significantly. Power-law fits to the data are indicated by the black lines at times *t* |Ωe| = 1090.1 (**a**) and 547.1 (**b**). The corresponding spectral indices are highlighted by the arrows.

**Figure 4.** Normalized magnetic field energy distribution *EB*(*<sup>k</sup>*)/*EB*<sup>0</sup> over the parallel wave number, *<sup>k</sup>* normalized to *<sup>c</sup>*/*ωp*, for simulations S1 (**a**) and S2 (**b**). The data are obtained at four times *<sup>t</sup>*|Ω*e*<sup>|</sup> as indicated. The spectra reach their steady state at the earliest time steps. Power-law fits to the data are indicated by the black lines at times *t* |Ωe| = 1090.1 (**a**) and 547.1 (**b**). The spectral indices are highlighted by the arrows. The dashed vertical lines indicatethe expected onset of cyclotron damping for purely parallel propagating waves.

The energy distribution *EB*(*<sup>k</sup>*, *<sup>k</sup>*⊥) in two-dimensional wave number space supports the claim that parallel energy transport becomes important in simulation S2, as Figure 5 shows. Figure 5b shows the distribution of magnetic field energy in simulation S2. Although hardly any (quasi-)parallel waves are produced above *<sup>k</sup> <sup>c</sup>*/*ω*<sup>p</sup> ≈ 1 (where cyclotron damping sets in), this critical parallel wave number can be passed at higher *k*⊥. At small wave numbers, however, the perpendicular cascade clearly dominates. In simulation S1, the situation is different, as Figure 5a shows. The perpendicular cascade at small wave numbers is similar to S2, as expected, but at larger *k*⊥, there is hardly any energy transport to higher parallel wave numbers, in agreement with the spectra given in Figures 3 and 4.

**Figure 5.** Two-dimensional magnetic field energy distribution in wave number space for simulations S1 (**a**) and S2 (**b**). Note different scales on the axes of (**a**) vs. (**b**).

#### *3.2. Simulation of Energetic Particles*

In order to study wave-particle scattering, a specific initialization of a test particle population is prepared. The ACRONYM code allows for different particle species (typically protons and electrons, but also positrons or heavier ions) and different particle populations (a background plasma and, for example, additional jet populations, non-thermal particles, etc.).

The simulations S1 and S2, discussed here, employ a thermal background plasma (see Table 1) and an additional population of non-thermal test particles to study the transport of energetic electrons. It is found that the test particles have no noticeable influence on the background plasma, even if the ratio *N*t/*N*bg of numerical particles in the test to the background particle population is of the order of unity.

#### 3.2.1. Initialization and Analysis

Test particles are initialized as a mono-energetic population; i.e., the particles have the same absolute speed, but their direction of motion is chosen randomly. The speed is calculated from the resonance condition for waves in the plasma. Solving Equation (18) for the speed of a particle of species *α* yields:

$$v\_a = \left| \frac{k\_{\parallel} \,\omega \, |\mu\_{\rm res}| \pm |\Omega\_{\rm n}| \, \sqrt{k\_{\parallel}^2 \mu\_{\rm res}^2 + (\Omega\_{\rm n}^2 - \omega^2)/c^2}}{k\_{\parallel}^2 \mu\_{\rm res}^2 + \Omega\_{\rm n}^2/c^2} \right|,\tag{47}$$

where *μ*res is the desired resonant pitch-angle cosine. The sign in the numerator changes depending on the polarization of the wave, its direction of propagation, and the particle species.

The directions of motion of the bulk of the test particles are chosen at random, using the speed calculated from Equation (47), a random polar angle cosine *μ*, and a random azimuth angle *φ*. This yields an isotropic distribution of the velocity vectors in *μ*-*φ* space. It is convenient to choose an isotropic distribution in *μ* = cos *θ* (instead of *θ*), because the analysis of pitch-angle scattering relies on the pitch-angle cosine and not on the pitchangle itself.

A fraction of the test particle population is not initialized as described above, but instead uses a parabolic distribution of polar angle cosines. This is done by assigning

$$
\mu = A \left( R + B \right)^{1/3} - \mathbb{C} \tag{48}
$$

to the particles, where *R* is a random number between zero and one, and *A*, *B*, and *C* are parameters describing the shape of the parabola. The parabolic distribution is required for the analysis of pitch-angle scattering using the method of Ivascenko et al. [40]. Ivascenko et al. suggest the use of a half-parabola, i.e., *A* = 2, *B* = 0, *C* = 1, but other distributions are also possible (see Figure 6). The resulting angular distribution of the entire test particle population is, therefore, not entirely isotropic.

**Figure 6.** A fraction of the test particle population has pitch-angle cosines assigned according to a parabolic distribution. **Left panel:** the assigned pitch-angle cosine *μ* as a function of the random number *R* ∈ [0, 1], which is used to generate the distribution. The curves follow Equation (48) and employ two sets of parameters *A*, *B*, and *C*, as indicated. **Right panel:** the resulting particle distribution *f*(*μ*) as a function of *μ*. The purple lines employ the parameters suggested by Ivascenko et al. [40], whereas the green curves show an improved implementation that is used in the ACRONYM code. Note that the derivative *d f*(*μ*)/*dμ* = 0 over the whole range of pitch-angle cosines in the latter case, whereas it becomes zero at *μ* = −1 when the parameters of Ivascenko et al. [40] are used.

The technique described above to create a population of energetic test particles for the study of wave-particle scattering was designed for a single plasma wave in the simulation [41]. However, it can also be applied to simulations with several plasma waves.

The test particle population is not injected at the start of the simulations S1 and S2 but at a later time *t*inj for the following reason: it is expected that turbulence develops from the initial conditions of the simulation, i.e., from a small set of seed waves that interact and start the turbulent cascade. This process takes time, and it may be desired to wait until a turbulent cascade is established before the transport of energetic test particles can be studied. Therefore, an optional deployment of test particles at later times in the simulation is favored. The particles are created at a a pre-defined time step, and the initialization is carried out as described earlier. Those particles can then be tracked for the rest of the simulation.

To evaluate particle transport in the turbulent plasma, the test particle data can be analyzed after the simulation to obtain the diffusion coefficient *Dμμ*. The diffusion coefficient is calculated from a simplified Fokker–Planck equation, Equation (22), where pitch-angle diffusion is assumed to be the only relevant diffusion process:

$$\frac{\partial f\_{\mu}}{\partial t} - \frac{\partial}{\partial \mu} D\_{\mu \mu} \frac{\partial f\_{\mu}}{\partial \mu} = 0. \tag{49}$$

This equation can be rewritten to yield

$$\frac{\partial f\_{\hbar}(\mu, t)}{\partial t} = \left(\frac{dD\_{\mu\mu}(\mu)}{d\mu}\right) \frac{\partial f\_{\hbar}(\mu, t)}{\partial \mu} + D\_{\mu\mu}(\mu) \frac{\partial^2 f\_{\hbar}(\mu, t)}{\partial \mu^2}.\tag{50}$$

The method described by Ivascenko et al. [40] is based on integrating Equation (49) over *μ*, which yields the pitch-angle current *jμ*:

$$\int\_{-1}^{\mu} \mu \frac{\partial f\_a(\mu, t)}{\partial t} = D\_{\mu \mu}(\mu) \frac{\partial f\_a(\mu, t)}{\partial \mu} = -j\_{\mu}. \tag{51}$$

The diffusion coefficient is then obtained by dividing *jμ* by *∂ fα*/*∂μ*.

#### 3.2.2. Physical Parameters

Using the setups of simulations S1 and S2, the transport of energetic electrons in kinetic turbulence is studied. In the following, the exact parameters for the test particle energy distribution are presented.

In simulations S1 and S2, decaying whistler turbulence is simulated, as was shown in the Section 3.2. As can be seen in the magnetic energy spectra presented in Figures 3 and 4, a steady state in terms of the power-law slope of the spectral energy distribution is established after a given time in each of the two simulations. As soon as this stage of the simulation is reached, a population of energetic test electrons can be injected as described in Section 3.2.1.

The time step for the checkpoint and subsequent restart is chosen to be *t* |Ω*e*| = 726.8 for S1 and *t* |Ω*e*| = 364.7 for S2. For each of these two setups, six test electron configurations are prepared. The simulations are labeled according to the physical setup (S1 or S2) followed by a letter referring to the test particle configuration ("a" through "f"). The parameters of the test particles can be found in Table 3 and describe the test electron speed *ve* and kinetic energy *Ee*.

**Table 3.** Test electron characteristics for the simulations of particle transport: test electron speed *ve* and corresponding kinetic energy *E*kin,e. The individual simulations (letters "a" through "f") are based on the simulations of kinetic turbulence S*j* with *j* = 1 and 2, which are described in Section 3.1.2 (Tables 1 and 2). Note that simulations S*j*e and S*j*f employ the same test electron energies. However, they differ in the way the test electron distribution is initialized (see text).


The test electron energy is increased from simulation S1a (S2a) to S1e (S2e). Simulation S1f (S2f) uses the same particle energy as S1e (S2e), but a different parabolic angular distribution of the particles. Here the particle density *f*(*μ*) increases with increasing *μ*, while in the other simulations it decreases with increasing pitch-angle cosine. This change in the pitch-angle distribution allows to check for systematic errors in the particle data.

#### **4. Results**

#### *Pitch-Angle Diffusion Coefficients*

The test particle simulations are analyzed as described in Section 3.2.1. The energetic electrons are tracked for several electron cyclotron time scales, and the resulting pitch-angle diffusion coefficients *Dμμ* are presented in Figures 7 and 8 for data based on the setup of S1 and S2, respectively. Time is measured as the interval Δ*t* from the time of the injection of the particles to the current time step.

The results of both sets of simulations, one based on S1 and the other based on S2, do not differ qualitatively, as would be expected from the two setups. The only difference between the physical parameters for S1 and S2 is the plasma temperature, which has no direct influence on the test electrons. Although the magnetic energy spectrum differs at high wave numbers (see Figure 3), the distribution of magnetic energy at small wave numbers is very similar. As it is explained below, this low-wave-number regime represents the

dominant influence on particle transport. Thus, the two sets of simulations are discussed simultaneously in what follows.

Although particle data can be obtained for an arbitrary number of time steps, the interval that can be used for the analysis is still limited. The method of Ivascenko et al. [40] critically depends on the particle distribution *f*(*μ*) in pitch-angle space. In order for the method to work, the initial distribution must be slightly disturbed, but the perturbations must not be too strong. This leaves only a brief period of time for the optimal efficiency of the method.

Figures 7a and 8a show the typical behavior of the derived *Dμμ* over time. Shortly after the injection of the test electrons, the perturbations of *f*(*μ*) are small, resulting in a low amplitude of *Dμμ* (purple lines). With increasing time, the amplitude grows and reaches a maximum (green and blue lines). At later times, the amplitude decreases again, as the perturbations become too strong and the method becomes unreliable (orange lines). The other panels of the two figures show the time evolution until the maximum amplitude of *Dμμ* is reached for other test electron energies.

Although the turbulent cascade is assumed to be symmetric about *μ* = 0, the panels of Figures 7 and 8 show an obvious asymmetry in the pitch angle diffusion coefficients derived from the test electron data. The amplitude of *Dμμ* is generally larger for *μ* < 0. While the energy spectrum itself is isotropic in *μ*, one could argue that the polarization of the waves' magnetic fields relative to the direction of the background magnetic field *B*<sup>0</sup> is different (i.e., the plasma physics definition of the polarization).

The magnetic helicity of the plasma waves is one of the possible causes of this anisotropy. Another reason for the asymmetry found in Figures 7 and 8 is that the parabolic distribution *f*(*μ*) of the test particles implies that there are more test electrons at negative pitch-angle cosines (except for simulations S1f and S2f). Therefore, the particle statistics is more reliable for negative *μ*, and the method of Ivascenko et al. [40] produces more accurate diffusion coefficients. While *Dμμ* can also be calculated for *μ* > 0, it is more prone to errors, and statistical fluctuations play a more important role, as Figure 9 indicates. However, small-scale statistical fluctuations can be suppressed by use of a Savitzky-Golay filter, as suggested by Ivascenko et al. [40].

Especially for early time steps, it can be seen that *Dμμ* is found to diverge at *μ* = 1. This is, of course, not a physical effect. At *μ* = 1, the derivative of the initial parabolic distribution *f*(*μ*) becomes (almost) zero. In this case, the method of Ivascenko et al. [40] becomes numerically unstable.

Another numerical effect causes *Dμμ* to become negative. This can be seen in Figure 8d and Figure 8e for early times. Negative solutions are most likely related to statistical fluctuations in the particle distribution, which drown the signal at early times, when the physically motivated perturbations of *f*(*μ*) are still developing.

Besides these flaws, the derived pitch-angle diffusion coefficients appear reasonable. They develop a (more or less) symmetric shape about *μ* = 0, indicating that neither direction is preferred. This is expected from the setup of the turbulence simulations S1 and S2, which employ a symmetric layout of initial waves and therefore should produce turbulent cascades that are symmetric in *μ*. This, however, cannot be proven by the plots of the energy distribution in wave number space, since the information about the direction of propagation of the waves is lost.

An interesting observation is that the pitch-angle diffusion coefficients grow in amplitude with the particle energy increasing from 100 keV to 2 MeV. At the highest test electron energy, however, the amplitude of *Dμμ* is significantly lower than in all other cases. Both Figures 7e and 8e also show that *Dμμ* forms a single peak close to *μ* = 0 in the case of the highest electron energy, whereas all other simulations produce a double peak structure. The reason for these differences is not clear. However, it is assumed that the different behavior of the 10 MeV electrons is related to their scattering characteristics. These high energy particles resonate with all of the initially excited waves in the simulations (see Figure 2), which is not the case in the simulations of less energetic electrons. Since the initial waves

contain the most energy, they are also assumed to significantly influence particle transport, especially if wave–particle resonances may occur.

In fact, the *Dμμ* in Figures 7e and 8e exhibit distinct peaks at early times (purple and green curves). Similar behavior is also found in simulations S1f and S2f, which are not included in Figures 7 and 8. For the example of one time step in simulation S1f, the peak structures in *Dμμ* are related to wave–particle resonances calculated according to Equation (18). The result is shown in Figure 10, where the colored vertical lines mark the expected positions of resonances. It can be seen that the resonances coincide with the positions of the peaks in *Dμμ*. The region around *μ* = 0 is most densely populated by resonances, which might explain the single peak in *Dμμ* at later times as seen in Figures 7e and 8e.

**Figure 7.** Pitch-angle diffusion coefficients *Dμμ* for test electrons with different energies for simulations S1a to S1e (**a**–**e**). The colored lines denote the diffusion coefficients derived from the simulation data at various times. The black lines represent the model predictions derived in Section 2.3.3.

**Figure 8.** Pitch-angle diffusion coefficients *Dμμ* for test electrons with different energies for simulations S2a to S2e (**a**–**e**). The colored lines denote the diffusion coefficients derived from the simulation data at various times. The black lines follow the model predictions derived in Section 2.3.3.

Finally, Figures 7 and 8 also include the model predictions from Equations (41) and (42). Some of the parameters required can be directly obtained from the setup of the simulations: the ratio *δB*2/*B*<sup>2</sup> <sup>0</sup> is listed in Table 1, and the test electron speed speed *ve* and the electron cyclotron frequency Ω*<sup>e</sup>* are found in Table 3. However, the minimum wave number *k*min, the spectral index *s*, and the cross helicity and magnetic helicity are not as trivial to find.

**Figure 9.** Test electron distribution in pitch-angle space: (**a**) the initial distribution *f*(*μ*) at the time of the injection of the test electrons (Δ*t* = 0) and at a later time in simulation S1c; (**b**) the relative deviation Δ*f* / ¯ *f* of the distributions at these two time steps. The deviation is defined as the difference of the two distributions over their mean value. Statistical fluctuations are visible to become more significant for larger *μ* values.

**Figure 10.** Pitch-angle diffusion coefficient *Dμμ* at one point in time as derived from the data of simulation S1f (black line). A noticeable number of peaks in *Dμμ* coincide with the positions of wave–particle resonances predicted by the resonance condition (18), which are marked by the colored, vertical lines. The colors denote the parallel wave numbers *<sup>k</sup>* (in numerical units) from one to four: purple, green, red, and orange. The line style refers to perpendicular wave numbers *k*<sup>⊥</sup> (also in numerical units) from zero to three: solid, dashed, dotted, dot-dashed. For example, the red dashed lines represent the resonance with a wave at (*<sup>k</sup>* <sup>=</sup> <sup>±</sup>3, *<sup>k</sup>*<sup>⊥</sup> <sup>=</sup> <sup>1</sup>). Only resonances of the first order, i.e., N = ±1 in Equation (18), are shown. Note that *<sup>k</sup>*<sup>⊥</sup> does not enter the resonance condition explicitly but is required to calculate the frequency *<sup>ω</sup>*(*<sup>k</sup>*, *<sup>k</sup>*⊥) according to the cold plasma dispersion relation.

For the minimum wave number, the magnetic energy spectra in Figures 3 and 4 have been considered. The second smallest resolved wave number *k*min = 2 Δ*k* has been chosen, where Δ*k* is the grid spacing in wave number space. In a square simulation box, where the

numbers of grid cells *N* and *N*<sup>⊥</sup> in the parallel and perpendicular directions are equal, the grid spacing is given by <sup>Δ</sup>*<sup>k</sup>* = <sup>2</sup> *<sup>π</sup>*/(*<sup>N</sup>* <sup>Δ</sup>*x*) = <sup>2</sup> *<sup>π</sup>*/(*N*<sup>⊥</sup> <sup>Δ</sup>*x*). The minimum wave number *k*min marks the beginning of the downward slope of the energy spectrum. Waves at small wave numbers are assumed to dominate the interaction with the particles due to their high energy content and the steep spectral slope. Therefore, the spectral index *<sup>s</sup>* = |*s*⊥| = 3.1 was chosen, in accordance with the index of the perpendicular spectrum in Figure 3a. This spectral index corresponds to simulation S1, but since the index in S2 is similar, *s* = 3.1 was used in both cases.

Finally, the magnetic helicity *σ* was chosen to be 0. The effect is in fact rather small since electrons mostly resonate with right-handed polarized modes.

From this starting point, the three parameters *k*min, *s*, and *HR* were fitted according to the numerical data from each simulation. The resulting parameters, which are used in the plots in Figures 7 and 8, are listed in Table 4. It can be seen that most simulations can be described with the initial choices for *k*min and *s* explained above.

**Table 4.** Parameters assumed for the model: spectral index *s*, cross-helicity *HR*, and minimum wave number *<sup>k</sup>*min. The latter is given in units of the grid spacing <sup>Δ</sup>*<sup>k</sup>* <sup>=</sup> {4.4, 8.7} × <sup>10</sup>−<sup>2</sup> *<sup>ω</sup>p*/*<sup>c</sup>* in wave number space in simulations S1 and S2, respectively.


In general, the model describes the data surprisingly well. Position and amplitude of the maxima and the inclination of the flanks are in good agreement. The contributions at *μ* = 0 are in disagreement; this is however not unexpected as for quasi-linear theory. Still, a non-zero contribution at medium energies is found, which is different from a nondispersive quasi-linear approach. The agreement of the model and the simulation results also supports the claim that the waves at small wave numbers dominate the interactions with the particles. Otherwise, the spectral index *s* would have to be changed according to the particle energy. The low energy particles, e.g., 100 keV, resonate with plasma waves in the high wave number regime, where the spectrum is steeper. Thus, according to Figures 3 and 4, the effective spectral index *s* should increase for these particles, if the resonant interactions with high-*k* waves were important. However, this seems not to be the case. But as the results in Figures 7 and 8 show, a change of the model equations for *Dμμ* is not necessary.

The model only fails for the simulations of 10 MeV-electrons (S1e, S1f, S2e, S2f). This might already be expected from the considerations discussed above: the high energy electrons are able to resonate with the initially excited waves. These waves contain the most energy and thus dominate the interaction of the particles with the turbulent spectrum. However, the initial waves cannot be considered to be part of the power-law spectrum itself. As Figures 3 and 4 show, the energy distribution forms a plateau at smallest wave numbers, where the initial waves are located. The initial waves are also only few in number, thus not forming a continuous spectrum but a population of distinct, individual waves. Representing a continuous spectrum on a discretized grid is always doomed to fail, but at larger wave numbers, the higher number of individual waves at least creates a rudimentary approximation of a continuum. Thus, the whole model assumption, i.e., a continuous power-law spectrum, is invalid. As Figure 10 shows, the pitch-angle diffusion coefficient derived from the simulation data can be described reasonably well by individual resonances with a number of waves.

Finally, it is worth taking a look at simulations S1f and S2f, which have not been discussed so far. These simulations, which employ the same test electron energies as S1e and S2e, were carried out to test whether the initial particle distribution *f*(*μ*) has an influence on the resulting *Dμμ*. It was already discussed above that the statistical fluctuations tend to become more noticeable at those *μ* where fewer particles are located. Thus, reversing the slope of the initial parabola should shift the dominant influence of statistical fluctuations from positive *μ* to negative.

Figure 11 depicts the pitch-angle diffusion coefficients derived from simulations S2e and S2f in panels a and b, respectively. This example is chosen because physical results for *Dμμ* are only obtained at negative *μ* at early times in S2e. Should this also be the case in S2f, this would mean that some physical process prefers the interaction of waves and energetic particles that propagate opposite to the background magnetic field. However, as Figure 11 shows, this is not the case. The pitch-angle diffusion coefficient derived from S2f appears to be more symmetric about *μ* = 0 at early times.

**Figure 11.** Comparison of the pitch-angle diffusion coefficients *Dμμ* derived from the test electron data of simulations S2e (**a**) and S2f (**b**). The two simulations differ by the slope of the parabolic particle distribution *f*(*μ*) used to initialize the test electrons. The asymmetry of *Dμμ* in S2e at early times is not reproduced by S2f, which suggests a numerical or statistical reason for the asymmetry. At late times, the *Dμμ* become similar, with a single peak near *μ* = 0 in both simulations.

At late times both simulations produce a single peak in *Dμμ*, which is located near *μ* = 0. The peak is slightly shifted to negative *μ* in both S2e and S2f. This may hint at a physical process leading to the peak not being centered exactly around *μ* = 0. Such an asymmetry is sometimes predicted in theoretical models, e.g., Schlickeiser [11]. However, considering the results of simulations of energetic particles and their interaction with individual waves, an asymmetry is not expected here.

Thus, the results of simulations S2e and S2f depicted in Figure 11 do not entirely agree with the expectations. It might be worthwhile to investigate the behavior of the pitchangle diffusion coefficient in more detail in a future project. Changing the initial particle distribution *f*(*μ*) once more (e.g., by altering the parameters *A*, *B*, and *C* in Equation (48)) or reversing the direction of the integration over *μ* in the method of [40] might help to distinguish between a physically motivated asymmetry and numerical artifacts.

#### **5. Conclusions**

In this paper, a set of pitch-angle diffusion coefficients for dispersive whistler waves are derived. Using a particle-in-cell code turbulence in the dispersive regime was simulated. Test particle electrons were injected into the simulated turbulence and their transport parameters were derived.

The conducted turbulence simulations yield power-law spectra of the magnetic field energy in wave number space. The measured spectral indices are in agreement with the findings of Refs. [27,29]. Numerical noise limits the energy spectra at high wave numbers, thus hindering the production of an energy cascade in the dissipation range.

While the theory is limited to parallel waves, simulations were performed in twodimensional wave-number space. The theoretical description of oblique, dispersive waves is not practically doable, while one-dimensional turbulence simulations are not producing an energy cascade. The approximation of a parallel spectrum makes this difference between dimensionalities reasonable.

The simulations of energetic particle transport in kinetic turbulence show that the steep energy spectrum leads to wave–particle interactions primarily in the low wave number regime. While low-energy particles, in principle, resonate with waves in the dispersive or dissipative regime of the turbulent cascade, these interactions are subordinate to interactions with non-resonant waves at lower wave numbers. The reason for this is that the energy content of dispersive waves decreases rapidly with increasing wave number due to the steep power-law spectrum. Thus, the waves at low wave numbers dominate the spectrum as far as particle transport is concerned.

This can be seen when comparing simulation data to the theoretical model. The test electron data from the simulations allows us to derive pitch-angle diffusion coefficients *Dμμ* using the method of [40]. The presented model for *Dμμ* in plasma turbulence with dispersive waves allows for the prediction of pitch-angle diffusion coefficient for Alfvén and whistler turbulence.

Simulation data and model match rather well for low-energy electrons. Contributions at *μ* = 0 are not modeled correctly as is expected for a quasi-linear model. The crosshelicity assumed in model parameters may not necessarily represent the cross-helicity of the plasma, but may be to some degree a numerical artifact. At higher electron energies, particles interact with the small number of excited plasma waves, which are used as a seed population for the generation of kinetic turbulence. The resulting *Dμμ* does not match the prediction for the interaction with the (continuous) turbulent spectrum but can be explained by resonant scattering with several waves at discrete wave numbers.

In general simulations, dispersive whistler turbulence and the corresponding particle transport are possible but are also still too expensive in terms of computing resources.

**Author Contributions:** Conceptualization, F.S., C.S. and R.S.; methodology, F.S. and R.S.; software, C.S.; writing—original draft preparation, F.S., C.S. and R.S.; writing—review and editing, F.S., C.S. and R.S.; visualization, F.S. and C.S.; funding acquisition, F.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Deutsche Forschungsgemeinschaft (DFG) through Grant No. SP 1124/9.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** F.S. would like to thank the Deutsche Forschungsgemeinschaft (DFG) for support through Grant No. SP 1124/9. The authors gratefully acknowledge the data storage service SDS@hd supported by the Ministry of Science, Research and the Arts Baden-Württemberg (MWK) and the German Research Foundation (DFG) through grant INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gausscentre.eu, accessed on 30 November 2021) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de, accessed on 30 November 2021) through grant pr84ti. We acknowledge the use of the *ACRONYM* code and would like to thank the developers (Verein zur Förderung kinetischer Plasmasimulationen e.V.) for their support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Borovsky, J.E.; Funsten, H.O. Role of solar wind turbulence in the coupling of the solar wind to the Earth's magnetosphere. *J. Geophys. Res. Space Phys.* **2003**, *108*, 1246. [CrossRef]

