A Quasi-linear Diffusion Model for Resonant Wave–Particle Instability in Homogeneous Plasma

To investigate the time evolution of the particle VDF through wave–particle resonances, we study the quasi-linear diffusion equation, given by Stix (1992):

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{\partial {f}_{j}}{\partial t}\right)}_{\mathrm{QLD}}=\mathop{\mathrm{lim}}\limits_{V\to \infty }\displaystyle \sum _{n=-\infty }^{\infty }\displaystyle \int \displaystyle \frac{\pi {q}_{j}^{2}}{{{Vm}}_{j}^{2}}\\ \ \ \times \,\hat{G}[{k}_{\parallel }]\left[\displaystyle \frac{{v}_{\perp }^{2}}{\left|{v}_{\parallel }\right|}\delta \left({k}_{\parallel }-\displaystyle \frac{{\omega }_{k}-n{{\rm{\Omega }}}_{j}}{{v}_{\parallel }}\right){\left|{\psi }_{j}^{n}\right|}^{2}\hat{G}[{k}_{\parallel }]{f}_{j}\right]{d}^{3}{\boldsymbol{k}},\end{array}\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{\psi }_{j}^{n} & \equiv & \displaystyle \frac{1}{\sqrt{2}}\left[{{E}}_{k}^{R}{e}^{i\phi }{J}_{n+1}({\rho }_{j})+{{E}}_{k}^{L}{e}^{-i\phi }{J}_{n-1}({\rho }_{j})\right]\\ & & +\,\displaystyle \frac{{v}_{\parallel }}{{v}_{\perp }}{{E}}_{k}^{z}{J}_{n}({\rho }_{j}),\end{array}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\hat{G}[{k}_{\parallel }]\equiv \left(1-\displaystyle \frac{{k}_{\parallel }{v}_{\parallel }}{{\omega }_{k}}\right)\displaystyle \frac{1}{{v}_{\perp }}\displaystyle \frac{\partial }{\partial {v}_{\perp }}+\displaystyle \frac{{k}_{\parallel }}{{\omega }_{k}}\displaystyle \frac{\partial }{\partial {v}_{\parallel }}.\end{eqnarray} \tag{ 3 }$$

The integer n determines the order of the resonance, where n = 0 corresponds to the Landau resonance and $n\ne 0$ corresponds to cyclotron resonances. In our equations, we label contributions from a given resonance order with a superscript n. The subscript j indicates the particle species. The particle VDF of species j is denoted as ${f}_{j}\equiv {f}_{j}({v}_{\perp },{v}_{\parallel },t)$, which is spatially averaged and gyrotropic, q_j and m_j are the charge and mass of a particle of species j, and ${v}_{\perp }$ and v_∥ are the velocity coordinates perpendicular and parallel with respect to the background magnetic field. We choose the coordinate system in which the proton bulk velocity is zero. We denote the nth-order Bessel function as ${J}_{n}({\rho }_{j})$, where ${\rho }_{j}\equiv {k}_{\perp }{v}_{\perp }/{{\rm{\Omega }}}_{j}$. The cyclotron frequency of species j is defined as ${{\rm{\Omega }}}_{j}\equiv {q}_{j}{B}_{0}/{m}_{j}c$,  B₀ is the background magnetic field, c is the speed of light, ${k}_{\perp }$ and k_∥ are the perpendicular and parallel components of the wavevector  k, and V is the volume in which the wave amplitude is effective so that the wave and particles undergo a significant interaction. We denote Dirac's δ function by δ and the azimuthal angle of wavevector  k by ϕ. The frequency ω is a complex function of  k, and we define ${\omega }_{k}$ as its real part and ${\gamma }_{k}$ as its imaginary part ($\omega ={\omega }_{k}+i{\gamma }_{k}$). Without loss of generality, we set ${\omega }_{k}\gt 0$. Furthermore, we assume that $| {\gamma }_{k}| \ll {\omega }_{k}$, i.e., the assumption of slow growth or damping that is central to quasi-linear theory.

The spatially Fourier-transformed electric field has the form of ${{\boldsymbol{E}}}_{{\boldsymbol{k}}}=\hat{x}{E}_{k}^{x}+\hat{y}{E}_{k}^{y}+\hat{z}{E}_{k}^{z}$ and is defined as (Gurnett & Bhattacharjee 2017)

$$\begin{eqnarray}&&{{\boldsymbol{E}}}_{{\boldsymbol{k}}}=\displaystyle \frac{1}{{\left(2\pi \right)}^{3/2}}\int {{\boldsymbol{E}}}_{{\boldsymbol{r}}}\exp \left[-i{\boldsymbol{k}}\cdot {\boldsymbol{r}}\right]{d}^{3}{\boldsymbol{r}},\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{k}}\cdot {\boldsymbol{r}}={k}_{\perp }x\cos \phi +{k}_{\perp }y\sin \phi +{k}_{\parallel }z$ and  r is the position vector. We take the constant background magnetic field as ${{\boldsymbol{B}}}_{0}=\hat{z}{B}_{0}$ and define the right- and left-circularly polarized components of the electric field as ${{E}}_{k}^{R}\,\equiv ({{E}}_{k}^{x}-i{{E}}_{k}^{y})/\sqrt{2}$ and ${{E}}_{k}^{L}\equiv ({{E}}_{k}^{x}+i{{E}}_{k}^{y})/\sqrt{2}$. The longitudinal component of the electric field is ${{E}}_{k}^{z}$.

Linear instabilities typically create fluctuations across a finite range of wavevectors. The Fourier transformation of such a wave packet in wavevector space corresponds to a wave packet in configuration space. For the sake of simplicity, we model this finite wave packet by assuming that the electric field ${{\boldsymbol{E}}}_{{\boldsymbol{r}}}$ of the unstable and resonant waves has the shape of a gyrotropic Gaussian wave packet,

$$\begin{eqnarray}&&{{\boldsymbol{E}}}_{{\boldsymbol{r}}}={{\boldsymbol{E}}}_{0}\exp \left[-\displaystyle \frac{{\sigma }_{\perp 0}^{2}{x}^{2}+{\sigma }_{\perp 0}^{2}{y}^{2}+{\sigma }_{\parallel 0}^{2}{z}^{2}}{2}\right]\exp \left[i{{\boldsymbol{k}}}_{0}\cdot {\boldsymbol{r}}\right],\end{eqnarray} \tag{ 5 }$$

where ${{\boldsymbol{E}}}_{0}=\hat{x}{E}_{0}^{x}+\hat{y}{E}_{0}^{y}+\hat{z}{E}_{0}^{z}$, ${{\boldsymbol{k}}}_{0}\cdot {\boldsymbol{r}}={k}_{\perp 0}x\cos \phi \,+{k}_{\perp 0}y\sin \phi +{k}_{\parallel 0}z$, and  k₀ is the wavevector of the Gaussian wave packet. We allow for an arbitrary angle θ₀ between  k₀ and ${{\boldsymbol{B}}}_{0}$, which defines the orientation of the wavevector at maximum growth of the wave, and assume that ${k}_{\parallel 0}\ne 0$. The vector ${{\boldsymbol{E}}}_{0}$ represents the peak amplitude of the electric field. The free parameters ${\sigma }_{\perp 0}$ and ${\sigma }_{\parallel 0}$ characterize the width of the Gaussian envelope. Quasi-linear theory requires that ${{\boldsymbol{E}}}_{{\boldsymbol{r}}}$ spatially averages to zero. Therefore, we assume that $| {k}_{\parallel 0}| \gg {\sigma }_{\parallel 0}$ so that the spatial dimension of the Gaussian wave packet is large compared to the parallel wavelength $2\pi /| {k}_{\parallel 0}|$.

The spatial Fourier transformation of Equation (5) according to Equation (4) then leads to

$$\begin{eqnarray}&&{{\boldsymbol{E}}}_{{\boldsymbol{k}}}=\displaystyle \frac{{{\boldsymbol{E}}}_{0}}{{\sigma }_{\parallel 0}{\sigma }_{\perp 0}^{2}}\exp \left[-\displaystyle \frac{{\left({k}_{\parallel }-{k}_{\parallel 0}\right)}^{2}}{2{\sigma }_{\parallel 0}^{2}}-\displaystyle \frac{{\left({k}_{\perp }-{k}_{\perp 0}\right)}^{2}}{2{\sigma }_{\perp 0}^{2}}\right].\end{eqnarray} \tag{ 6 }$$

We identify V with the volume of the Gaussian envelope, $V=1/({\sigma }_{\parallel 0}{\sigma }_{\perp 0}^{2})$. Equation (6) describes the localization of the wave energy density in wavevector space. For the instability analysis through Equation (6), we define the unstable  k spectrum as the finite wavevector range in which ${\gamma }_{k}\gt 0$ and argue that resonant waves exist only in this unstable  k spectrum. We ignore any waves outside this  k spectrum because they are damped.

We define k_∥0 as the value of k_∥ at the center of the unstable  k spectrum. We then obtain

$$\begin{eqnarray}&&{k}_{\perp 0}={k}_{\parallel 0}\tan {\theta }_{0}.\end{eqnarray} \tag{ 7 }$$

In the case of linear plasma instability, we identify ${k}_{\perp 0}$ and k_∥0 with the wavevector components at which the instability has its maximum growth rate as a reasonable approximation. To approximate the wave frequency of the unstable waves at the angle θ₀ of maximum growth, we expand ${\omega }_{k}$ of the unstable and resonant waves around k_∥0 as

$$\begin{eqnarray}&&{\omega }_{k}({k}_{\parallel })\approx {\omega }_{k0}+{v}_{g0}\left({k}_{\parallel }-{k}_{\parallel 0}\right),\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{\left.{v}_{g0}\equiv \displaystyle \frac{\partial {\omega }_{k}}{\partial {k}_{\parallel }}\right|}_{{k}_{\parallel }={k}_{\parallel 0}}.\end{eqnarray} \tag{ 9 }$$

In Equations (8) and (9), ${\omega }_{k0}$ and v_g0 are the wave frequency and parallel group velocity of the unstable and resonant waves, evaluated at ${k}_{\parallel }={k}_{\parallel 0}$. We select the values of ${\sigma }_{\perp 0}$ and ${\sigma }_{\parallel 0}$ as the half widths of the perpendicular and parallel unstable  k spectrum. In the case of linear plasma instability, the numerical values for ${k}_{\perp 0}$, k_∥0, ${\sigma }_{\perp 0}$, ${\sigma }_{\parallel 0}$, ${\omega }_{k0}$ and v_g0 can be found from the solutions of the hot-plasma dispersion relation, which thus closes our set of equations.

By using Equations (6) and (8), we rewrite Equation (1) as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial {f}_{j}}{\partial t}\right)}_{\mathrm{QLD}}=\displaystyle \sum _{n=-\infty }^{\infty }\int \hat{G}[{k}_{\parallel }]\left[{D}_{j}^{n}\hat{G}[{k}_{\parallel }]{f}_{j}\right]{d}^{3}{\boldsymbol{k}},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{D}_{j}^{n} & \equiv & \displaystyle \frac{\pi {q}_{j}^{2}{v}_{\perp }^{2}}{{\sigma }_{\parallel 0}{\sigma }_{\perp 0}^{2}{m}_{j}^{2}}\delta ({k}_{\parallel }-{k}_{\parallel j}^{n})\\ & & \times \,\displaystyle \frac{| {\psi }_{j0}^{n}{| }^{2}}{| {v}_{\parallel }-{v}_{g0}| }\exp \left[-\displaystyle \frac{{\left({k}_{\parallel }-{k}_{\parallel 0}\right)}^{2}}{{\sigma }_{\parallel 0}^{2}}-\displaystyle \frac{{\left({k}_{\perp }-{k}_{\perp 0}\right)}^{2}}{{\sigma }_{\perp 0}^{2}}\right],\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\psi }_{j0}^{n} & \equiv & \displaystyle \frac{1}{\sqrt{2}}\left[{{E}}_{0}^{R}{e}^{i\phi }{J}_{n+1}({\rho }_{j})+{{E}}_{0}^{L}{e}^{-i\phi }{J}_{n-1}({\rho }_{j})\right]\\ & & +\,\displaystyle \frac{{v}_{\parallel }}{{v}_{\perp }}{{E}}_{0}^{z}{J}_{n}({\rho }_{j}),\end{array}\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&{k}_{\parallel j}^{n}\equiv \displaystyle \frac{{\omega }_{k0}-{k}_{\parallel 0}{v}_{g0}-n{{\rm{\Omega }}}_{j}}{{v}_{\parallel }-{v}_{g0}}.\end{eqnarray} \tag{ 13 }$$

We set ${{E}}_{0}^{R}=({{E}}_{0}^{x}-i{{E}}_{0}^{y})/\sqrt{2}$ and ${{E}}_{0}^{L}=({{E}}_{0}^{x}+i{{E}}_{0}^{y})/\sqrt{2}$ as constant, evaluated at ${{\boldsymbol{k}}}_{0}$.

Equation (10) is the quasi-linear diffusion equation describing the action of the dominant wave–particle instability and coexisting damping contributions from other resonances in a Gaussian wave packet. We define the n resonance as the contribution to the summation in Equation (10) with only integer n. We note that any n resonance can contribute to wave instability or to wave damping depending on the resonance's characteristics.

We define stabilization as the process that creates the condition in which ${(\partial {f}_{j}/\partial t)}_{\mathrm{QLD}}\to 0$ for all v_⊥ and v_∥. For our analysis of the stabilization of a VDF through a resonant wave–particle instability, including coexisting damping effects, we use Boltzmann's H theorem, in which the quantity H is defined as

$$\begin{eqnarray}&&H(t)\equiv \int {f}_{j}({\boldsymbol{v}},t)\mathrm{ln}{f}_{j}({\boldsymbol{v}},t){d}^{3}{\boldsymbol{v}}.\end{eqnarray} \tag{ 14 }$$

By using Equation (10), the time derivative of H is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dH}}{{dt}}=\displaystyle \sum _{n=-\infty }^{\infty }\int \int (\mathrm{ln}{f}_{j}+1)\hat{G}[{k}_{\parallel }]\left[{D}_{j}^{n}\hat{G}[{k}_{\parallel }]{f}_{j}\right]{d}^{3}{\boldsymbol{k}}{d}^{3}{\boldsymbol{v}}.\end{eqnarray} \tag{ 15 }$$

The integrand in Equation (15) is equivalent to

$$\begin{eqnarray}&&\begin{array}{l}(\mathrm{ln}{f}_{j}+1)\hat{G}[{k}_{\parallel }]\left[{D}_{j}^{n}\hat{G}[{k}_{\parallel }]{f}_{j}\right]\\ \ \ =\ \hat{G}[{k}_{\parallel }]\left[{D}_{j}^{n}\hat{G}[{k}_{\parallel }]({f}_{j}\mathrm{ln}{f}_{j})\right]-{D}_{j}^{n}{\left[\hat{G}[{k}_{\parallel }]{f}_{j}\right]}^{2}/{f}_{j}.\end{array}\end{eqnarray} \tag{ 16 }$$

Upon substituting Equation (16) into Equation (15), the first term on the right-hand side in Equation (16) disappears after the integration over  v. Then, by resolving the δ function in D_jⁿ through the k_∥ integral, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{dH}}{{dt}}=-\displaystyle \sum _{n=-\infty }^{\infty }{\left(\displaystyle \frac{{dH}}{{dt}}\right)}^{n},\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dH}}{{dt}}\right)}^{n}\equiv \int \left\{{\widetilde{D}}_{j}^{n}{\left[\hat{G}[{k}_{\parallel j}^{n}]{f}_{j}\right]}^{2}/{f}_{j}\right\}{d}^{3}{\boldsymbol{v}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\widetilde{D}}_{j}^{n} & \equiv & {W}_{j}^{n}\displaystyle \frac{\pi {q}_{j}^{2}{v}_{\perp }^{2}}{{\sigma }_{\parallel 0}{\sigma }_{\perp 0}^{2}{m}_{j}^{2}}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{\infty }| {\psi }_{j0}^{n}{| }^{2}\\ & & \times \,\exp \left[-\displaystyle \frac{{\left({k}_{\perp }-{k}_{\perp 0}\right)}^{2}}{{\sigma }_{\perp 0}^{2}}\right]{k}_{\perp }{{dk}}_{\perp }d\phi ,\end{array}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{W}_{j}^{n}\equiv \displaystyle \frac{1}{| {v}_{\parallel }-{v}_{g0}| }\exp \left[-\displaystyle \frac{{k}_{\parallel 0}^{2}}{{\sigma }_{\parallel 0}^{2}}{\left(\displaystyle \frac{{v}_{\parallel }-{v}_{\parallel \mathrm{res}}^{n}}{{v}_{\parallel }-{v}_{g0}}\right)}^{2}\right],\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{v}_{\parallel \mathrm{res}}^{n}\equiv \displaystyle \frac{{\omega }_{k0}-n{{\rm{\Omega }}}_{j}}{{k}_{\parallel 0}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}\begin{array}{rcl}\hat{G}[{k}_{\parallel j}^{n}] & \equiv & \left[\displaystyle \frac{n{{\rm{\Omega }}}_{j}}{{\omega }_{k0}-{k}_{\parallel 0}{v}_{g0}-n{{\rm{\Omega }}}_{j}}\right]\displaystyle \frac{{v}_{\parallel }-{v}_{g0}}{{v}_{\mathrm{ph}}{v}_{\perp }}\displaystyle \frac{\partial }{\partial {v}_{\perp }}\\ & & +\,\displaystyle \frac{1}{{v}_{\mathrm{ph}}}\displaystyle \frac{\partial }{\partial {v}_{\parallel }},\end{array}\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&{v}_{\mathrm{ph}}\equiv \displaystyle \frac{{\omega }_{k}({k}_{\parallel j}^{n})}{{k}_{\parallel j}^{n}}=\displaystyle \frac{({\omega }_{k0}-{k}_{\parallel 0}{v}_{g0}){v}_{\parallel }-n{{\rm{\Omega }}}_{j}{v}_{g0}}{{\omega }_{k0}-{k}_{\parallel 0}{v}_{g0}-n{{\rm{\Omega }}}_{j}}.\end{eqnarray} \tag{ 23 }$$

The function ${\widetilde{D}}_{j}^{n}$ in Equation (19) plays the role of a diffusion coefficient for the n resonance. In ${\widetilde{D}}_{j}^{n}$, the v_∥ function W_jⁿ defined in Equation (20) serves as a window function that determines the region in v_∥ space in which the quasi-linear diffusion through the n resonance is effective. The window function W_jⁿ is maximum at ${v}_{\parallel \mathrm{res}}^{n}$ defined in Equation (21), which is the parallel velocity of the particles that resonate with the waves at ${k}_{\parallel }={k}_{\parallel 0}$ through the n resonance. Our window function W_jⁿ is linked to Dirac's δ function in the limit

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{v}_{\parallel \mathrm{res}}^{n}\to {v}_{g0}}{W}_{j}^{n}\approx \sqrt{\pi }\displaystyle \frac{{\sigma }_{\parallel 0}}{| {k}_{\parallel 0}| }\delta ({v}_{\parallel }-{v}_{g0}),\end{eqnarray} \tag{ 24 }$$

where $| {k}_{\parallel 0}| \gg {\sigma }_{\parallel 0}$. Through this ordering between $| {k}_{\parallel 0}|$ and ${\sigma }_{\parallel 0}$, we assume that W_jⁿ restricts a finite region in v_∥ space and that the W_jⁿ for different resonances do not overlap with each other in v_∥ space.

Only particles distributed within W_jⁿ experience the n resonance and contribute to the quasi-linear diffusion, which is ultimately responsible for the stabilization. Because all terms in Equation (18) are positive semidefinite, all resonances independently stabilize f_j through quasi-linear diffusion in the v_∥ range defined by their respective W_jⁿ, according to Equation (17). Therefore, H decreases and ${dH}/{dt}$ tends toward zero during the quasi-linear diffusion through all resonances while f_j is in the process of stabilization. When f_j reaches a state of full stabilization through all n resonances, the instability has saturated and its growth ends.

The v_∥ function ${k}_{\parallel j}^{n}$ defined in Equation (13) is the resonant parallel wavenumber, fulfilling the condition that ${k}_{\parallel j}^{n}={k}_{\parallel 0}$ at ${v}_{\parallel }={v}_{\parallel \mathrm{res}}^{n}$. It quantifies the k_∥ component of the unstable  k spectrum in the v_∥ range defined by W_jⁿ. Equation (23) defines the phase velocity at ${k}_{\parallel j}^{n}$, which is only constant when ${v}_{g0}={\omega }_{k0}/{k}_{\parallel 0}$, in which case ${v}_{\mathrm{ph}}={v}_{g0}$ for all v_∥. We discuss the diffusion operator $\hat{G}[{k}_{\parallel j}^{n}]$ in Equation (22) in the next section.

According to Equation (18), unless the wave amplitude is zero, the condition for achieving stabilization through the n resonance is

$$\begin{eqnarray}&&\hat{G}[{k}_{\parallel j}^{n}]{{\mathbb{F}}}_{j}^{n}({v}_{\perp },{v}_{\parallel })=0,\end{eqnarray} \tag{ 25 }$$

where ${{\mathbb{F}}}_{j}^{n}({v}_{\perp },{v}_{\parallel })$ represents the stabilized VDF of species j through the n resonance. In Equation (25), $\hat{G}[{k}_{\parallel j}^{n}]$ is a directional derivative along the isocontour of ${{\mathbb{F}}}_{j}^{n}$ evaluated at a given velocity position. Considering the role of W_jⁿ, $\hat{G}[{k}_{\parallel j}^{n}]$ describes only the diffusion of resonant particles within W_jⁿ along the isocontour of ${{\mathbb{F}}}_{j}^{n}$. Consequently, the particles experiencing the n resonance diffuse toward the stable state so that ${({dH}/{dt})}^{n}\to 0$, while the isocontours of ${{\mathbb{F}}}_{j}^{n}$ describe the diffusive velocity-space trajectories for the n resonance.

To find such a trajectory, we express an infinitesimal variation of ${{\mathbb{F}}}_{j}^{n}$ along an isocontour as

$$\begin{eqnarray}&&d{{\mathbb{F}}}_{j}^{n}=\displaystyle \frac{\partial {{\mathbb{F}}}_{j}^{n}}{\partial {v}_{\perp }}{{dv}}_{\perp }+\displaystyle \frac{\partial {{\mathbb{F}}}_{j}^{n}}{\partial {v}_{\parallel }}{{dv}}_{\parallel }=0.\end{eqnarray} \tag{ 26 }$$

Equations (22) and (26) allow us to rewrite Equation (25) as

$$\begin{eqnarray}&&{v}_{\perp }{{dv}}_{\perp }+\left[\displaystyle \frac{n{{\rm{\Omega }}}_{j}}{n{{\rm{\Omega }}}_{j}-{\omega }_{k0}+{k}_{\parallel 0}{v}_{g0}}\right]({v}_{\parallel }-{v}_{g0}){{dv}}_{\parallel }=0.\end{eqnarray} \tag{ 27 }$$

By integrating Equation (27), the diffusive trajectory for the n resonance is then given by

$$\begin{eqnarray}&&{v}_{\perp }^{2}+\left[\displaystyle \frac{n{{\rm{\Omega }}}_{j}}{n{{\rm{\Omega }}}_{j}-{\omega }_{k0}+{k}_{\parallel 0}{v}_{g0}}\right]{\left({v}_{\parallel }-{v}_{g0}\right)}^{2}=\mathrm{const}.\end{eqnarray} \tag{ 28 }$$

Kennel & Engelmann (1966) treat the two limiting cases in which ${v}_{g0}={\omega }_{k0}/{k}_{\parallel 0}$ and ${v}_{g0}=0$. Using their assumptions, our Equation (28) is equivalent to their Equation (4.8) if ${v}_{g0}={\omega }_{k0}/{k}_{\parallel 0}$, and our Equation (28) is equivalent to their Equation (4.11) if ${v}_{g0}=0$. Depending on the dispersion properties of the resonant waves, Equation (28) is either an elliptic or a hyperbolic equation when $n\ne 0$. In the case of electron resonances, it is safe to assume that

$$\begin{eqnarray}&&\displaystyle \frac{n{{\rm{\Omega }}}_{j}}{n{{\rm{\Omega }}}_{j}-{\omega }_{k0}+{k}_{\parallel 0}{v}_{g0}}\geqslant 0\end{eqnarray} \tag{ 29 }$$

in Equation (28) if ${v}_{g0}\lt ({\omega }_{k0}+n| {{\rm{\Omega }}}_{e}| )/{k}_{\parallel 0}$ for all positive n and ${v}_{g0}\gt ({\omega }_{k0}+n| {{\rm{\Omega }}}_{e}| )/{k}_{\parallel 0}$ for all negative n. However, in the case of proton resonances, resonant waves are more likely to violate Equation (29) because ${{\rm{\Omega }}}_{p}\ll | {{\rm{\Omega }}}_{e}|$.

Figure 1 illustrates the diffusive flux of particles experiencing two arbitrary resonances: the n₁ and n₂ resonances for an unstable wave. The dark shaded areas represent isocontours of the VDFs of two particle populations in velocity space. The red and dark-blue solid curves represent the diffusive trajectories according to Equation (28) with $n={n}_{1}$ and $n={n}_{2}$, assuming that the resonant wave fulfills Equation (29). The window functions ${W}_{j}^{{n}_{1}}$ and ${W}_{j}^{{n}_{2}}$ describe the v_∥ ranges in which the n₁ and n₂ resonances are effective. The light-blue dashed semicircles correspond to contours of constant kinetic energy in the proton rest frame, for which

$$\begin{eqnarray}&&{v}_{\perp }^{2}+{v}_{\parallel }^{2}=\mathrm{const}.\end{eqnarray} \tag{ 30 }$$

Zoom In Zoom Out Reset image size

In general, the diffusive flux is always directed from higher to lower phase-space densities during the process of stabilization. At point A, resonant particles in ${W}_{j}^{{n}_{1}}$ diffuse along the red solid curve toward smaller v_∥. Considering the relative alignment between the diffusive flux and the constant-energy contour at point A, the diffusing particles lose kinetic energy. This energy is transferred to the resonant wave, which consequently grows in amplitude. Therefore, this situation corresponds to an instability of the resonant wave. At point B, particles do not diffuse along the red solid curve because this point lies outside ${W}_{j}^{{n}_{1}}$.

At point C, resonant particles in ${W}_{j}^{{n}_{2}}$ diffuse along the dark-blue solid curve toward greater v_∥. Considering the relative alignment between the diffusive flux and the constant-energy contour at point C, the diffusing particles gain kinetic energy. This energy is taken from the resonant wave, which consequently shrinks in amplitude. Therefore, this situation corresponds to the damping of the resonant wave and counteracts the driving of the instability through the n₁ resonance. Because the resonant wave is unstable, the n₁ resonant instability must overcome the counteracting n₂ resonant damping.

According to Equation (18), there are three factors that determine the diffusion rate for the action of an n resonance. The first factor is the particle density f_j within W_jⁿ. The second factor is ${\widetilde{D}}_{j}^{n}$, whose magnitude is determined by the polarization properties of the resonant waves. The third factor is the quantity $\hat{G}[{k}_{\parallel j}^{n}]{f}_{j}/{f}_{j}$, which defines the relative alignment between the isocontours of f_i and the diffusive flux along the diffusion trajectory within W_jⁿ. In Figure 1, the magnitude of $| \hat{G}[{k}_{\parallel j}^{{n}_{1}}]{f}_{j}/{f}_{j}|$ at point A is greater than the magnitude of $| \hat{G}[{k}_{\parallel j}^{{n}_{2}}]{f}_{j}/{f}_{j}|$ at point C.

Because the diffusive flux is directed from higher to lower values of f_j, the quantity $\hat{G}[{k}_{\parallel j}^{n}]{f}_{j}/{f}_{j}$ resolves the ambiguity in the directions of the trajectories for resonant particles. A careful analysis of $\hat{G}[{k}_{\parallel j}^{n}]$ using Equation (29) shows that, if $({k}_{\parallel }/| {k}_{\parallel }| )(\hat{G}[{k}_{\parallel j}^{n}]{f}_{j}/{f}_{j})\gt 0$ at a given resonant velocity, resonant particles diffuse toward a smaller value of v_∥ along the diffusive trajectory, while if $({k}_{\parallel }/| {k}_{\parallel }| )(\hat{G}[{k}_{\parallel j}^{n}]{f}_{j}/{f}_{j})\lt 0$ at a given resonant velocity, resonant particles diffuse toward a greater value of v_∥.

To simulate the VDF evolution and to compare the diffusion rates between resonances quantitatively, a rigorous numerical analysis of Equation (10) is necessary. For this purpose, we develop a mathematical approach based on the Crank–Nicolson scheme and present the mathematical details in Appendix A. Our approach is applicable to all two-dimensional diffusion equations with off-diagonal diffusion terms. Our numerical solution, given by Equation (A28), evolves the VDF under the action of multiple resonances in one time step. We tested our numerical solution by showing that the diffusive flux obeys the predicted diffusion properties discussed in Section 2.3.

To find the characteristics of the unstable oblique FM/W wave, we numerically solve the hot-plasma dispersion relation with the NHDS code (Verscharen & Chandran 2018). We use the same plasma parameters as Verscharen et al. (2019), which are, notwithstanding the wide range of natural variation, representative for the average electron parameters in the solar wind (Wilson et al. 2019). We assume that the initial plasma consists of isotropic Maxwellian protons, core electrons, and strahl electrons. The subscripts p, e, c, and s indicate protons, electrons, electron core, and electron strahl, respectively.

We choose our coordinate system so that the anti-sunward and obliquely propagating FM/W waves have ${k}_{\parallel }\gt 0$. We set ${\beta }_{c}={\beta }_{p}=1$ and ${\beta }_{s}=0.174$, where ${\beta }_{j}\equiv (8\pi {n}_{j}{k}_{B}{T}_{j})/{B}_{0}^{2}$, n_j and T_j are the density and temperature of species j, and k_B is the Boltzmann constant. We set n_p = n_e, ${n}_{c}=0.92{n}_{p}$, ${n}_{s}=0.08{n}_{p}$, T_c = T_p, and ${T}_{s}=2{T}_{p}$. In the proton rest frame, we set ${n}_{c}{U}_{c}+{n}_{s}{U}_{s}=0$. We initialize the core and strahl bulk velocity with ${U}_{c}/{v}_{{Ae}}=-0.22$ and ${U}_{s}/{v}_{{Ae}}=2.52$, where ${v}_{{Ae}}\equiv {B}_{0}/\sqrt{4\pi {n}_{e}{m}_{e}}$ is the electron Alfvén speed. NHDS finds that, under these plasma parameters, ${\gamma }_{k}\gt 0$ at angles between ${\theta }_{0}=51^\circ$ and ${\theta }_{0}=67^\circ$. Our strahl bulk velocity then provides a maximum growth rate of ${\gamma }_{k}/| {{\rm{\Omega }}}_{e}| ={10}^{-3}$ (Verscharen et al. 2019).

Figure 2 shows ${\gamma }_{k}$ and ${\omega }_{k}$ as functions of the k_∥ component of the wavevector  k for three different θ₀. The oblique FM/W instability has its maximum growth rate at ${\theta }_{0}=55^\circ$, while ${\gamma }_{k}\gt 0$ for $0.21\lesssim {k}_{\parallel }{v}_{{Ae}}/| {{\rm{\Omega }}}_{e}| \lesssim 0.28$, which is the parallel unstable  k spectrum. As defined in Section 2.1, we acquire ${k}_{\parallel 0}{v}_{{Ae}}/| {{\rm{\Omega }}}_{e}| \approx 0.245$. This value with Equations (7)–(9) leads to ${k}_{\perp 0}{v}_{{Ae}}/| {{\rm{\Omega }}}_{e}| =0.35$, ${\omega }_{k0}/| {{\rm{\Omega }}}_{e}| \approx 0.07$, and ${v}_{g0}/{v}_{{Ae}}\approx 0.86$. We also acquire ${\sigma }_{\parallel 0}{v}_{{Ae}}/| {{\rm{\Omega }}}_{e}| \approx 0.035$ and ${\sigma }_{\perp 0}{v}_{{Ae}}/| {{\rm{\Omega }}}_{e}| \approx 0.05$ from the unstable  k spectrum.

Zoom In Zoom Out Reset image size

Using the wave and plasma parameters from the previous section, we describe the electron strahl and core diffusion in velocity space. In our analysis, we only consider the n = +1, −1, and 0 resonances, ignoring higher-n resonances due to their negligible contributions.

Upon substituting our wave parameters into Equation (20), we quantify the dimensionless window functions ${W}_{e}^{n}{v}_{{Ae}}$ with n = +1, −1, and 0. In Figure 3, the red, dark-blue and orange lines represent ${W}_{e}^{+1}{v}_{{Ae}}$, ${W}_{e}^{-1}{v}_{{Ae}}$ and ${W}_{e}^{0}{v}_{{Ae}}$, which are maximum at ${v}_{\parallel \mathrm{res}}^{+1}/{v}_{{Ae}}=4.37$, ${v}_{\parallel \mathrm{res}}^{-1}/{v}_{{Ae}}=-3.8$ and ${v}_{\parallel \mathrm{res}}^{0}/{v}_{{Ae}}=0.29$, respectively. We reiterate that the superscripts indicate the n resonance. The black line indicates ${v}_{\parallel }={v}_{g0}$. Each ${W}_{e}^{n}{v}_{{Ae}}$ shows the v_∥ range in which the quasi-linear diffusion through each resonance is effective. We note that the ${W}_{e}^{n}{v}_{{Ae}}$ for the three resonances have different widths in v_∥ space and maximum values due to the different magnitudes of $| {v}_{\parallel \mathrm{res}}^{n}-{v}_{g0}|$ (see Equation (24)). By substituting our wave parameters into Equation (28), the diffusive trajectories for the n = +1, −1, and 0 resonances are given by

$$\begin{eqnarray}&&{\left({v}_{\perp }/{v}_{{Ae}}\right)}^{2}+1.16{\left({v}_{\parallel }/{v}_{{Ae}}-0.86\right)}^{2}=\mathrm{const},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{\left({v}_{\perp }/{v}_{{Ae}}\right)}^{2}+0.88{\left({v}_{\parallel }/{v}_{{Ae}}-0.86\right)}^{2}=\mathrm{const},\end{eqnarray} \tag{ 32 }$$

and

$$\begin{eqnarray}&&{\left({v}_{\perp }/{v}_{{Ae}}\right)}^{2}=\mathrm{const}.\end{eqnarray} \tag{ 33 }$$

Equations (31) and (32) describe ellipses with their axes oriented along the v_⊥ and v_∥ directions. In Equation (33), the perpendicular velocity of resonant particles is constant.

Zoom In Zoom Out Reset image size

Figure 4 illustrates the electron diffusion from these three resonances. We show the v_∥ ranges in which the three resonances are effective according to ${W}_{e}^{+1}{v}_{{Ae}}$, ${W}_{e}^{-1}{v}_{{Ae}}$, and ${W}_{e}^{0}{v}_{{Ae}}$ from Figure 3. The red, dark-blue, and orange solid curves represent the contours given by Equations (31)–(33), respectively. The light-blue dashed semicircles correspond to constant-energy contours in the proton rest frame (see Equation (30)). The black line indicates ${v}_{\parallel }={v}_{g0}$. For the initial strahl and core VDF, we apply the plasma parameters in Section 3.1 to the dimensionless Maxwellian distribution:

$$\begin{eqnarray}&&{f}_{j}^{M}=\displaystyle \frac{{n}_{j}{v}_{{Ae}}^{3}}{{\pi }^{3/2}{n}_{p}{v}_{\mathrm{th},j}^{3}}\exp \left[-\displaystyle \frac{{v}_{\perp }^{2}+{\left({v}_{\parallel }-{U}_{j}\right)}^{2}}{{v}_{\mathrm{th},j}^{2}}\right],\end{eqnarray} \tag{ 34 }$$

where ${v}_{\mathrm{th},j}\equiv \sqrt{2{k}_{{\rm{B}}}{T}_{j}/{m}_{j}}$. The red and blue areas in Figure 4 represent f_s^M and f_c^M, which are normalized by the maximum value of f_c^M and plotted up to a value of 10⁻⁵. In this normalization, Figure 4 does not reflect the relative density ratio between both electron species.

Zoom In Zoom Out Reset image size

Due to the v_∥ profile of ${W}_{e}^{+1}$, the n = +1 resonance has a significant effect on f_s^M. As discussed in Section 2.3, because $({k}_{\parallel }/| {k}_{\parallel }| )(\hat{G}[{k}_{\parallel e}^{+1}]{f}_{s}^{M}/{f}_{s}^{M})\gt 0$, this resonance leads to a diffusion of the resonant strahl electrons in ${W}_{e}^{+1}$ along trajectories represented by the red arrows. According to Equation (30), the phase-space trajectory of particles that diffuse without a change in kinetic energy is described by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{{dv}}_{\perp }}{{{dv}}_{\parallel }}\right)}_{E}=-\displaystyle \frac{{v}_{\parallel }}{{v}_{\perp }}.\end{eqnarray} \tag{ 35 }$$

According to Equation (31), the phase-space trajectory of resonant particles fulfilling the n = +1 resonance, indicated by the superscript +1, is described by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{{dv}}_{\perp }}{{{dv}}_{\parallel }}\right)}^{+1}=-1.16\displaystyle \frac{{v}_{{Ae}}}{{v}_{\perp }}\left(\displaystyle \frac{{v}_{\parallel }}{{v}_{{Ae}}}-0.86\right).\end{eqnarray} \tag{ 36 }$$

Comparing Equations (35) and (36) in ${W}_{e}^{+1}$ shows that $| {({{dv}}_{\perp }/{{dv}}_{\parallel })}^{+1}| \lt | {({{dv}}_{\perp }/{{dv}}_{\parallel })}_{E}|$ for the resonant electrons. Therefore, resolving the ambiguity in the directions of the trajectories, the distance of resonant strahl electrons from the origin of the coordinate system decreases. This decrease in ${v}_{\perp }^{2}+{v}_{\parallel }^{2}$ represents a loss of kinetic energy of the resonant strahl electrons. The n = +1 resonance, therefore, contributes to the driving of the FM/W instability.

Due to the v_∥ profile of ${W}_{e}^{-1}$, the n = −1 resonance has a significant effect on f_c^M. Because $({k}_{\parallel }/| {k}_{\parallel }| )(\hat{G}[{k}_{\parallel e}^{-1}]{f}_{c}^{M}/{f}_{c}^{M})\lt 0$, this resonance leads to a diffusion of the resonant core electrons in ${W}_{e}^{-1}$ along trajectories represented by the dark-blue arrows. According to Equation (32), the phase-space trajectory of resonant particles fulfilling the n = −1 resonance, indicated by the superscript −1, is described by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{{dv}}_{\perp }}{{{dv}}_{\parallel }}\right)}^{-1}=-0.88\displaystyle \frac{{v}_{{Ae}}}{{v}_{\perp }}\left(\displaystyle \frac{{v}_{\parallel }}{{v}_{{Ae}}}-0.86\right).\end{eqnarray} \tag{ 37 }$$

Comparing Equations (35) and (37) in ${W}_{e}^{-1}$ shows that $| {({{dv}}_{\perp }/{{dv}}_{\parallel })}^{-1}| \gt | {({{dv}}_{\perp }/{{dv}}_{\parallel })}_{E}|$ for the resonant electrons. Therefore, resolving the ambiguity in the directions of the trajectories, the distance of resonant core electrons from the origin of the coordinate system increases. This increase in ${v}_{\perp }^{2}+{v}_{\parallel }^{2}$ represents a gain of kinetic energy of the resonant core electrons. The n = −1 resonance, therefore, counteracts the FM/W instability through the n = +1 resonance.

Due to the v_∥ profile of W_e⁰, the n = 0 resonance has a significant effect on electrons in the v_∥ range in which ${f}_{c}^{M}\gt {f}_{s}^{M}$ and $\partial {f}_{c}^{M}/\partial {v}_{\parallel }\lt 0$. Because $({k}_{\parallel }/| {k}_{\parallel }| )(\hat{G}[{k}_{\parallel e}^{0}]{f}_{c}^{M}/{f}_{c}^{M})\lt 0$, the resonant electrons in W_e⁰ diffuse along trajectories represented by the yellow arrows. Because the distance of these electrons from the origin of the coordinate system increases, these resonant electrons diffuse toward greater kinetic energies. This diffusion removes energy from the resonant FM/W waves and thus counteracts the driving of the FM/W instability through the n = +1 resonance.

Figure 4 only illustrates the nature of the quasi-linear diffusion through the n = +1, −1, and 0 resonances in velocity space. It does not give any information regarding the relative strengths of the diffusion rates between the three resonances. Because the FM/W wave is unstable according to linear theory, the n = +1 resonant instability must dominate over any counteracting contributions from the n = −1 and 0 resonances.

We use our numerical procedure from Equation (A28) to simulate the quasi-linear diffusion through the n = +1, −1, and 0 resonances, predicted in Section 3.2. According to the definitions in Appendix A, we select the discretization parameters N_v = 60 ${v}_{\perp \max }/{v}_{{Ae}}={v}_{\parallel \max }/{v}_{{Ae}}=7$, and $| {{\rm{\Omega }}}_{e}| {\rm{\Delta }}t=1$. For the computation of Equation (A28), we use the same parameters of resonant FM/W waves as those presented in Section 3.1 and quantify ${\widetilde{D}}_{e}^{\pm 1}$ and ${\widetilde{D}}_{e}^{0}$ in Equation (19).

In ${\widetilde{D}}_{e}^{n}$ for each resonance, we only consider the J₀ term in ${\psi }_{j0}^{n}$, ignoring higher-order Bessel functions due to their small contributions. Our NHDS solutions show that $| {E}_{0}^{y}| \,\approx 0.39| {E}_{0}^{x}|$ and $| {E}_{0}^{z}| \approx 0.28| {E}_{0}^{x}|$ in the unstable  k spectrum. Then, we set $| {{E}}_{0}^{R}| \approx | {{E}}_{0}^{L}| \approx 0.76| {{E}}_{0}^{x}|$. Faraday's law yields ${{E}}_{0}^{x}\,\approx [{\omega }_{k0}/({k}_{\parallel 0}c)]{{B}}_{0}^{y}$ when ignoring the small contributions from any E₀^z terms. This allows us to express ${{E}}_{0}^{x}$ through ${{B}}_{0}^{y}$ in ${\psi }_{j0}^{n}$, where B₀^y represents the peak amplitude of the wave magnetic-field fluctuations. For simplicity, we assume that ${{B}}_{0}^{y}$ is constant in time during the quasi-linear diffusion. Under these assumptions, we acquire

$$\begin{eqnarray}\begin{array}{rcl}{\widetilde{D}}_{e}^{\pm 1} & \approx & {W}_{e}^{\pm 1}\displaystyle \frac{0.58{\pi }^{2}| {{\rm{\Omega }}}_{e}{| }^{2}{v}_{\perp }^{2}}{{\sigma }_{\parallel 0}{\sigma }_{\perp 0}^{2}}{\left[\displaystyle \frac{{B}_{0}^{y}}{{B}_{0}}\displaystyle \frac{{\omega }_{k0}}{{k}_{\parallel 0}}\right]}^{2}\\ & & \times \,{\displaystyle \int }_{0}^{\infty }{J}_{0}{\left({\rho }_{e}\right)}^{2}\exp \left[-\displaystyle \frac{{\left({k}_{\perp }-{k}_{\perp 0}\right)}^{2}}{{\sigma }_{\perp 0}^{2}}\right]{k}_{\perp }{{dk}}_{\perp }\end{array}\end{eqnarray} \tag{ 38 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{\widetilde{D}}_{e}^{0} & \approx & {W}_{e}^{0}\displaystyle \frac{0.16{\pi }^{2}| {{\rm{\Omega }}}_{e}{| }^{2}{v}_{\parallel }^{2}}{{\sigma }_{\parallel 0}{\sigma }_{\perp 0}^{2}}{\left[\displaystyle \frac{{B}_{0}^{y}}{{B}_{0}}\displaystyle \frac{{\omega }_{k0}}{{k}_{\parallel 0}}\right]}^{2}\\ & & \times \,{\displaystyle \int }_{0}^{\infty }{J}_{0}{\left({\rho }_{e}\right)}^{2}\exp \left[-\displaystyle \frac{{\left({k}_{\perp }-{k}_{\perp 0}\right)}^{2}}{{\sigma }_{\perp 0}^{2}}\right]{k}_{\perp }{{dk}}_{\perp },\end{array}\end{eqnarray} \tag{ 39 }$$

where the relative amplitude ${B}_{0}^{y}/{B}_{0}$ is a free parameter, and we set ${B}_{0}^{y}/{B}_{0}=0.001$. Then, we apply Equations (38) and (39) to Equation (A28).

We initialize our numerical computation with the same f_s^M and f_c^M as defined in Section 3.2. Figure 5(a) represents the normalized ${f}_{e}={f}_{c}^{M}+{f}_{s}^{M}$, plotted up to a value of 10⁻⁵. Figure 5(b) shows f_e evolved through the n = +1, −1, and 0 resonances, resulting from our iterative calculation of Equation (A28). Considering the maximum value of the instability's growth rate, ${\gamma }_{k}/| {{\rm{\Omega }}}_{e}| =4.8\times {10}^{-3}$ in Figure 2, we terminate the evaluation of our numerical computation at $| {{\rm{\Omega }}}_{e}| t=5\times {10}^{2}$, which corresponds to ${\gamma }_{k}t\sim 1$ and thus a reasonable total growth of the unstable FM/W waves.

The strahl electrons at around ${v}_{\parallel }/{v}_{{Ae}}\approx 4.4$ diffuse through the n = +1 resonance, as theoretically predicted in Figure 4. This diffusion increases the pitch angle of the resonant strahl electrons and generates a strong pitch-angle gradient at ${v}_{\parallel }/{v}_{{Ae}}\approx 3.8$. During this process, the ${v}_{\perp }$ of the scattered strahl electrons increases while their v_∥ decreases.

Because the longitudinal component of the electric-field fluctuations is weaker than their transverse components, the diffusion through the n = 0 resonance is only slightly noticeable over the modeled time interval. The diffusion through the n = −1 resonance is not noticeable even though ${\widetilde{D}}_{e}^{-1}$ and ${\widetilde{D}}_{e}^{+1}$ have similar magnitudes. This is because $| \hat{G}[{k}_{\parallel e}^{-1}]{f}_{e}/{f}_{e}|$ in ${W}_{e}^{-1}$ is much smaller than $| \hat{G}[{k}_{\parallel e}^{+1}]{f}_{e}/{f}_{e}|$ in ${W}_{e}^{+1}$, as discussed in Section 2.3, and the number of core electrons in ${W}_{e}^{-1}$ is very small (see Figures 4 and 5).