5.1.1 Derivation

For the remainder of § 5.1, we assume that the whistlers pitch-angle scatter electrons in a frame moving at the whistler phase speed, $v_{{w}}$. This assumption is consistent with the quasi-linear operator for slow electromagnetic modes (see § 5.2 of Kennel & Engelmann Reference Kennel and Engelmann1966) and is the same as chosen by Drake et al. (Reference Drake, Pfrommer, Reynolds, Ruszkowski, Swisdak, Einarsson, Thomas, Hassam and Roberg-Clark2021). What follows is a brief outline of the full calculation; see Appendix A for more details. The model collision operator is given by

(5.1)\begin{equation} C_{\rm CE}[f]=\frac{\partial }{\partial \xi^{\prime}}\left[\frac{1-\xi^{\prime 2}}{2}\nu_{\rm CE}(v^{\prime},\xi^{\prime})\frac{\partial f}{\partial \xi^{\prime}}\right], \end{equation}

where the prime denotes quantities evaluated in the reference frame of the whistlers; note that $C_{\rm CE}[f]$ is independent of the gyroangle, $\phi$. We take $v_{{w}}$, $\beta _e^{-1}$, the electron gyroradius and the effective electron mean free path to all be small quantities of the same order $\epsilon$, viz.

(5.2)\begin{equation} \epsilon\sim\frac{v_{{w}}}{v_{\text{th}e}} \sim \frac{1}{\beta_e}\sim \frac{\rho_e}{L_T}\sim\frac{\lambda_{\text{mfp},e}}{L_T}\ll 1. \end{equation}

Expanding $C_{\rm CE}[f]$ in $\epsilon$, we find that the electron kinetic equation evaluated to leading order in $\epsilon$ is

(5.3)\begin{equation} 0 = \varOmega_e\frac{\partial f_{e0}}{\partial \phi}+ \frac{\partial }{\partial \xi}\left[\frac{1-\xi^{2}}{2}\nu_{\rm CE}(w,\xi)\frac{\partial f_{e0}}{\partial \xi}\right], \end{equation}

where $\boldsymbol {w}=\boldsymbol {v}-\boldsymbol {u}_e$ is the velocity of the electrons in the frame of any bulk flow $\boldsymbol {u}_e$. This constraint is satisfied trivially by an isotropic $f_{e0}=f_{e0}(v)$. To next order, after gyroaveraging, we obtain the correction equation for parallel transport,

(5.4)\begin{equation} w_\parallel\nabla_{{\parallel}} f_{e0}+\frac{\nabla_{{\parallel}} p_e}{m_en_e} \frac{w_\parallel}{w}\frac{{\rm d} f_{e0}}{{\rm d} w} =\frac{\partial }{\partial \xi}\left[\frac{1-\xi^2}{2}\nu_{\rm CE}(w,\xi)\left(\frac{\partial \langle\, f_{e1}\rangle_{\phi}}{\partial \xi} + v_{{w}}\frac{{\rm d} f_{e0}}{{\rm d} w}\right) \right]. \end{equation}

The term proportional to $v_{{w}}$ arises as a consequence of the ordering of the whistler phase speed and ensures that the scattering occurs in the wave frame. The correction equation (5.4) can be solved for $\nu _{\rm CE}$ by using $w_\parallel =w\xi$ and integrating both sides with respect to $\xi$. Choosing the integration constant to keep ${\partial \langle\, f_{e1}\rangle _\phi /\partial \xi }$ finite, we find that

(5.5)\begin{equation} \nu_{\rm CE}(w,\xi)= \left. -\left(w\nabla_\parallel f_{e0}+\frac{\nabla_{{\parallel}}p_e}{m_en_e}\frac{{\rm d} f_{e0}}{{\rm d} w}\right) \middle/ \left(\frac{\partial \langle\, f_{e1}\rangle_{\phi}}{\partial \xi} + v_w\frac{{\rm d} f_{e0}}{{\rm d} w}\right) \right. . \end{equation}

5.1.2 Calculating the distribution function

In order to evaluate (5.5), we require expressions for $f_{e0}$ and $\langle\, f_{e1}\rangle _\phi$. We calculate these quantities from our simulations by binning individual particle data on a finite phase-space grid. To do so, we break up the central $60\,\%$ of the simulation domain into $4\,\%$ intervals and integrate over the $y$ direction; this ensures that we sample a distribution function with nearly uniform temperature in each interval. In velocity space we bin using a $60\times 60$ grid in $(v,\xi )$ coordinates, from $v=0$ to $v=3v_{\text {th}e}$ (where $v_{{\rm th}e}$ is the electron thermal speed local to that interval) and from $\xi =-1$ to $\xi =1$. Because the magnitude of the fluid velocity $\boldsymbol {u}_e$ in our simulations is at most a few percent of the whistler phase speed $v_{{w}}$, we do not distinguish between $v$ and $w$ in what follows; in real astrophysical contexts with bulk flows, $\boldsymbol {u}_e$ may not be small compared with $v_{{w}}$. A consequence of our choice to ignore the $\phi$ coordinate is that our measured distributions only contain parallel transport, naturally precluding any transport perpendicular to the background field.

The algorithm is as follows. For every electron at a certain time step, we calculate the bin indices that correspond to its location in velocity-space coordinates ($v,\xi$) and add $1$ to that bin. After normalizing, this procedure yields $f_e(v,\xi )$ for each $4\,\%$ interval of the simulation domain. We then determine $f_{e0}$ by requiring that all fluid quantities ($n_e$, $\boldsymbol {u}_e$, $T_e$) are contained in $f_{e0}$, so that the transport is determined solely by $f_{e1}$ (Krommes Reference Krommes2018), viz.

(5.6)\begin{equation} \left(n_e,\boldsymbol{u}_e,p_e\right)\equiv\int \,{\rm d}^3\boldsymbol{v}\left(1,\frac{\boldsymbol{v}}{n_e},\frac{m_e}{3}w^2\right)f_e=\int \,{\rm d}^3\boldsymbol{v}\left(1,\frac{\boldsymbol{v}}{n_e},\frac{m_e}{3}w^2\right)f_{e0}, \end{equation}

where $p_e$ is the scalar electron pressure associated with $T_e$. As a consequence of this procedure, $f_{e0}$ is naturally isotropic, as required by the Chapman–Enskog expansion at zeroth order in $\epsilon$ (5.3). Finally, we calculate the first-order deviation in $f_e$ from the definition

(5.7)\begin{equation} f_{e1}=f_e - f_{e0} , \end{equation}

so that

(5.8)\begin{equation} \int \,{\rm d}^3\boldsymbol{v}\, (1,\boldsymbol{v},v^2)f_{e1}=0 , \end{equation}

i.e. $f_{e1}$ is fully kinetic.

The only subtlety in the above procedure is the choice of $f_{e0}$, since (5.3) demands only that $f_{e0}$ is isotropic. In what follows, we choose $f_{e0}$ to be an isotropic Maxwellian distribution: $f_{e0}(v,\xi )=f_{{\rm M}e}(v)$. Alternatively, we could have chosen $f_{e0}=\langle\, f_e\rangle _{\phi,\xi }$, which puts fewer restrictions on $f_{e0}$. Thankfully, the exact choice of $f_{e0}$ does not substantially change the following analysis or the resulting computed collision frequency. To determine $\nu _{\rm CE}(v,\xi )$, we first normalize each $4\,\%$ domain segment to a single ${v_{\text {th}e}}$ and then average both $f_{e0}$ and $f_{e1}$ across all segments. Given $f_{e0}(v)$ and $\nabla _\parallel \ln T_e$ from each simulation, we then need only to determine $v_{{w}}$ and $\partial \langle\, f_{e1}\rangle _{\phi }/\partial \xi$. The former is straightforwardly computed for each simulation using the method explained in § 4.3. Unfortunately, the latter presents some difficulties, which are addressed in the next subsection.

5.1.3 Numerical results

In principle, the Chapman--Enskog method can resolve the full velocity-space dependence of the collision frequency (5.5), at least within the model operator. However, in practice we find that spurious zeros in ${\partial \langle\, f_{e1}\rangle _\phi /\partial \xi }$ caused by statistical noise in the binned distribution function as well as shallow gradients in $f_{e1}$ result in a noisy calculation of $\nu _{\rm CE}(v,\xi )$, which can have unphysical (i.e. negative) values in certain regions of phase space. See figure 9 for plots of $\langle\, f_{e1}\rangle _\phi /f_{e0}$ and $(\partial \langle\, f_{e1}\rangle _\phi /\partial \xi + v_{{w}} \,{\rm d} f_{e0}/{\rm d} w)/f_{e0}$ for run b40x4 at reduced resolutions of $10\times 10$ and $20\times 20$, which exhibit noise and negative values for $v/v_{\text {th}e}\gtrsim 2$. Increasing the grid resolution only increases and worsens the artifacts already present in figure 9. It is possible that even lower resolutions might produce favourable numbers; however, much fewer than $10\times 10$ grid cells runs the risk of not resolving the structure of $\nu _{\rm CE}$ in both $v$ and $\xi$. To circumvent these issues, we average $\partial f_{e1}/\partial \xi$ over the pitch angle to find a $v$-dependent collision frequency:

(5.9)\begin{equation} \nu_{\rm CE}(v)= \left. -\left(v\nabla_\parallel f_{e0}+\frac{\nabla_{{\parallel}}p_e}{m_en_e}\frac{{\rm d} f_{e0}}{{\rm d} v}\right) \right/ \left(\left\langle\frac{\partial f_{e1}}{\partial \xi}\right\rangle_{\!\xi,\phi} + v_w\frac{{\rm d} f_{e0}}{{\rm d} v}\right) . \end{equation}

A particularly unfortunate consequence of having to average over $\xi$ is that any resonant structure in $\nu _{\rm CE}(v,\xi )$ (e.g. due to Landau-cyclotron resonances occurring along lines of constant $v_\parallel$) is lost. A distinct benefit, however, is that we are able to resolve precisely the velocity dependence of the scattering frequency.

Figure 9.

Two-dimensional plots of $\langle\, f_{e1}\rangle _\phi /f_{e0}$ (a,c) and $(\partial \langle\, f_{e1}\rangle /\partial \xi + v_{{w}} \,{\rm d} f_{e0}/{\rm d} w)/f_{e0}$ (b,d) for run b40x4 at grid resolution $10\times 10$ (a,b) and $20\times 20$ (c,d). The noise and negative values for $v/v_{\text {th}e}\gtrsim 2$ only get worse at increasing resolution and imply a $\nu _{\rm CE}(v,\xi )$ that is noisy and contains unphysical negative values. Dashed lines correspond to constant $v_{\parallel }$.

In figure 10 we plot (5.9) as a function of speed $v$ for each of the $\boldsymbol {\nabla } p_0=0$ runs. The discontinuities and spikes in the range $1\lesssim v/v_{\text {th}e} \lesssim 2$ are the result of zeros in the denominator not exactly cancelling with the zero in the numerator at $v/v_{\text {th}e}=\sqrt {5/2}\simeq 1.6$. Within this model, subthermal electrons are scattered at a rate that is weakly dependent upon speed, with only a slight increase as $v$ becomes very subthermal. Superthermal electrons, on the other hand, have a pitch-angle-scattering rate that increases steeply with $v$, approximately as $(v/v_{\text {th}})^{3}$ (black dashed line) though with some significant spread. (Note that this is rather different than in Coulomb-collisional plasmas, for which the scattering rate decreases with $v$ at superthermal speeds.) While the model operator is unable to resolve the pitch-angle dependence of the scattering operator, it does serve as an important data point for comparison with the other methods used in §§ 5.2 and 5.3.

Figure 10.

Chapman–Enskog whistler scattering frequency $\nu _{\rm CE}$ (5.9) as a function of speed $v/v_{\text {th}e}$. A Gaussian filter with standard deviation of one cell ($v_{\text {th}e}/20$ and $1/30$ in $v$ and $\xi$, respectively) is applied for smoothing. The distribution function used to calculate $\nu _{\rm CE}$ was taken at the end of each of the runs (see table 1). To guide the eye, we include a black dashed line ${\propto }(v/v_{\text {th}e})^3$.

5.2.1 Definition

We adopt the electromagnetic quasi-linear operator given in Stix (Reference Stix1992) (also given by (B1)–(B3)) under two simplifying assumptions: (i) that the turbulence is only two dimensional (as in our simulations), and (ii) that the phase speed of the scattering whistlers relative to the electron thermal velocity is small and order $\epsilon$, i.e.

(5.10)\begin{equation} \epsilon \sim \frac{v_{{w}}}{v_{\text{th}e}} \sim \frac{1}{\beta_e}\sim \frac{\rho_e}{L}\sim\frac{\lambda_{\text{mfp},e}}{L}\ll 1. \end{equation}

The first assumption allows us to write the original quasi-linear diffusion operator (B1), which is expressed in terms of the electric field, in terms of the magnetic and parallel electric fields (see (B4)). The second assumption allows us to expand the electron kinetic equation in $\epsilon \ll 1$ (cf. § 5.1.1). At zeroth order, we find that the quasi-linear operator is simply a pitch-angle-scattering operator associated with the magnetic-field fluctuations. At this order, there is no Landau damping from the parallel electric field, but there is transit-time damping as the guiding centres of Landau-resonant particles surf the mirror force associated with low-frequency fluctuations in the magnetic-field strength (Stix Reference Stix1992; Barnes Reference Barnes1966). At first order in $\epsilon$, we find a correction equation in which the distortion in the distribution function caused by the parallel temperature gradient is balanced by two terms: pitch-angle scattering of $f_{e1}$ and a term proportional to $v_{{w}}$. Except for the nature of the collision frequency, this correction equation is identical to that associated with a collision operator that pitch-angle scatters in a frame co-moving with whistler waves at $v_{{w}}$ (cf. (5.4)). We refer the reader to Appendix B for a more detailed derivation.

Keeping only relevant terms, the simplified operator is

(5.11)\begin{equation} C_{\rm QL}[f]=\frac{\partial }{\partial \xi} \left[\frac{1-\xi^2}{2}\nu_{\rm QL}(v,\xi) \left(\frac{\partial f}{\partial \xi} + v_{{w}}\frac{\partial f}{\partial v}\right)\right] , \end{equation}

in which the quasi-linear collision frequency is given by

(5.12a)\begin{equation} \nu_{\rm QL}(v,\xi)\simeq 2{\rm \pi} \frac{\varOmega_e}{N_xN_y} \sum_n \sum_{k_x}\sum_{k_y}\delta\left(\frac{\omega(k_{{\parallel}}\rho_e)}{\varOmega_e}-k_{{\parallel}}\rho_e\frac{v_{{\parallel}}}{v_{\text{th}e}}+n \right) \left|\frac{\varPsi_{n,k}}{B_0}\right|^2 . \end{equation}

The magnetic field enters this expression via

(5.12b)\begin{equation} \varPsi_{n,k} \simeq B_k^- {\rm J}_{n-1}(z)+ B_k^+ {\rm J}_{n+1}(z) , \end{equation}

with the argument of the Bessel functions ${\rm J}_{n}$ being $z=k_yv_{\perp } /\varOmega _e$, and the complex quantities

(5.12c)\begin{equation} B_k^\pm{=}\frac{B_{kz}\pm {\rm i} B_{ky}}{2} \end{equation}

are the Fourier-transformed magnetic fields corresponding to right-handed ($+$) and left-handed ($-$) polarizations. For parallel waves, $k_y=k_z=0$ and the cyclotron resonances correspond to purely right-handed modes for $n=-1$ and left-handed modes for $n=1$. For $n=0$, $\varPsi _{0,k}={\rm i} B_{ky}{\rm J}_{1}(z)= -{\rm i} (k_x/k_y)B_{kx} {\rm J}_{1}(z)$; at long wavelengths such that $z\ll 1$ and ${\rm J}_{1}(z)\approx z/2$, we have $2\varOmega _e|\varPsi _{0,k}/B_0|^2 \approx (v^2_\perp /2\varOmega _e) k_x^2 |B_{kx}/B_0|^2$, which corresponds to transit-time damping.

5.2.2 Resonance condition

The argument of the delta function in (5.12a), commonly called the resonance condition, defines the parallel wavenumber that resonantly interacts with an electron at a given parallel velocity. Using the whistler dispersion relation (4.5), the resonance condition is

(5.13)\begin{equation} {-}{\rm sign}(n)\frac{0.23}{\beta_e}(k_{{\parallel}}\rho_e)(k\rho_e)-k_{{\parallel}}\rho_e\frac{v_{{\parallel}}}{v_{\text{th}e}}+n=0. \end{equation}

The factor of ${\rm sign}(n)$ in front of the whistler dispersion relation ensures that (5.13) has the correct solution for electrons with $v_{\parallel }/v_{\text {th}e}>0$ in resonance with oblique whistler waves; we take ${\rm sign}(0)=-1$. The general solution to (5.13) for quasi-parallel ($k_{\perp }/k_{\parallel }\ll 1$) waves is

(5.14)\begin{align} k_{{\parallel}}\rho_e={-}{\rm sign}(n)\frac{1}{2}\frac{\beta_e}{0.23}\frac{v_{{\parallel}}}{v_{\text{th}e}}\left[1\pm \sqrt{1+4\left(\frac{\beta_e}{0.23}\frac{v_{{\parallel}}}{v_{\text{th}e}}\right)^{{-}2}\left(\frac{\beta_e}{0.23}|n|+(k_{{\perp}}\rho_e)^2\right)}\right]. \end{align}

For $n=0$, exactly parallel whistlers have their resonance at

(5.15)\begin{equation} k_{{\parallel}}\rho_e=\frac{\beta_e}{0.23}\frac{v_{{\parallel}}}{v_{\text{th}e}}\quad (n=0), \end{equation}

excluding $k_{\parallel }\rho _e=0$. For $n\ne 0$, we expand (5.14) in $\beta _e\gg 1$ to find the resonance conditions

(5.16)\begin{equation} k_{{\parallel}}\rho_e\approx n\left(\frac{v_{{\parallel}}}{v_{\text{th}e}}\right)^{{-}1}\quad (n\ne 0) , \end{equation}

as well as

(5.17)\begin{equation} k_{{\parallel}}\rho_e={-}{\rm sign}(n)\frac{\beta_e}{0.23}\frac{v_{{\parallel}}}{v_{\text{th}e}} + n\frac{v_{\text{th}e}}{v_{{\parallel}}}\quad (n\ne 0). \end{equation}

We omit the solution (5.17) for the reason that, in the $\beta _e\gg 1$ limit, it yields a negative $k_{\parallel }$ for $|v_{\parallel }|/v_{\text {th}e}\sim 1$ electrons and a proper choice of $n$; such wave vectors are stable to the HWI. Keeping only solutions (5.15) and (5.16) and performing some simple manipulations of the delta function, we have, for purely parallel waves,

(5.18)\begin{equation} \delta\left(\frac{\omega(k_{{\parallel}}\rho_e)}{\varOmega_e}-k_{{\parallel}}\rho_e\frac{v_{{\parallel}}}{v_{\text{th}e}}+n\right)=\left|\frac{v_{\text{th}e}}{v_{{\parallel}}}\right| \begin{cases} \delta\left(k_{{\parallel}}\rho_e-\dfrac{\beta_e}{.23}\dfrac{v_{{\parallel}}}{v_{\text{th}e}}\right), & n=0,\\ \delta\left(k_{{\parallel}}\rho_e-n\dfrac{v_{\text{th}e}}{v_{{\parallel}}}\right), & n\ne 0. \end{cases} \end{equation}

For the numerical results provided in § 5.2.3, we solve the full oblique resonance condition (5.13) numerically.

5.2.3 Numerical results

In figure 11 we plot $\nu _{\rm QL}(v,\xi )L_T/(\rho _{e,0}\beta _{e,0}\varOmega _{e,0})$ from (5.12) for simulations b40 and b40x4. One property of the operator that is readily apparent is the asymmetry about $\xi =0$. The scattering frequency for $\xi <0$ is an order of magnitude higher than that for $\xi >0$, and is consistent with the fraction of energy in the left-hand-polarized component of an elliptically polarized oblique whistler wave (cf. (2.11)). As discussed in § 2.2, whistler waves appear to be left-hand polarized to any electron with $v_{\parallel }>v_{{w}}$, so no gyroresonance can occur. Resonance can only be facilitated by the left-hand component of oblique whistlers, which will appear right-handed in the frame of these electrons. This asymmetry is prominent and a significant refinement over the results in the previous section.

Figure 11.

Two-dimensional plots of the quasi-linear pitch-angle collision frequency $\nu _{\rm QL}(v,\xi )$ (5.12) for all $\beta _{e,0}=40$ runs normalized to $\beta _{e,0}\rho _{e,0}/L_T\varOmega _{e,0}$. Dashed lines correspond to contours of constant $v_{\parallel }$.

It is also clear from figure 11 that the collision frequency is approximately constant along lines of constant $v_{\parallel }$ (dashed lines), a functional dependence that is expected given the form of (5.12). Suppose the magnetic energy scales with the parallel wavenumber as ${\delta B^2/B_0^2\propto k_{\parallel }^{-\alpha }}$; ignore any $k_y$ dependence for simplicity. This scaling implies a resonant scattering frequency $\nu _{\rm QL}\propto |v_{\parallel }|^{\alpha -1}$. In figure 12 we plot $\nu _{\rm QL}(v,\xi )$ as a function of $v_{\parallel }$ for electrons with $\xi <0$ (a) and $\xi >0$ (b). Regardless of the sign of $\xi$, the scattering frequency grows as a power law in $|v_{\parallel }|$. For electrons with $v_{\parallel } < 0$, $\nu _{\rm QL} \sim (v_{\parallel }/v_{\text {th}e,0})^{2.5}$, implying a field spectrum ${\sim }k_{\parallel }^{-3.5}$ that is compatible with the results shown in figure 6. There is a break in this power law at $|v_{\parallel }|/v_{\text {th}e}\sim 2$, which corresponds to the break in magnetic energy spectrum around $k_x\rho _{e,0}\sim 0.6$ for the $n=\pm 1$ principal resonances. There is another break near $|v_{\parallel }|/v_{\text {th}e}\sim 0.3$ that is due to contributions from both the $(n=\pm 1)$ cyclotron resonances and the $(n=0)$ resonance. The contribution from the principal cyclotron resonances is the smallest of the two and corresponds to the transition of the spectrum from being dominated by the $k_x\rho _{e,0}>1$ cascade to being dominated by noise with a spectral index of $0$. The dominant contribution is from the $n=0$ resonance and is the result of transit-time damping. The consequence is scattering at very small parallel velocities $v_{\parallel }/v_{\text {th}e}\sim v_{{w}}/v_{\text {th}e}\le 0.1$ at a rate generally much smaller than that at thermal velocities; however, the effect is most pronounced for positive parallel velocities.

Figure 12.

Quasi-linear collision frequency calculated for each of the runs as a function of $v_{\parallel }$. Electrons with $\xi <0$ are plotted on (a) and those with $\xi >0$ are on (b); the former are scattered at a rate ${\sim }10$ times faster than the latter due to the relative amounts of energy in the right- and left-hand-polarized components of the wave spectrum. We include lines with power laws $(v_{\parallel }/v_{\text {th}e,0})^{2.5}$ and $(v_{\parallel }/v_{\text {th}e,0})^{2}$ on (a) and (b), respectively.

In general we find rough agreement with $\nu _{\rm QL}/\varOmega _e\propto \beta _{e,0}\rho _{e,0}/L_T$, matching the predicted scaling (2.8) and the observed scaling using our Chapman–Enskog method in § 5.1. In figure 13 we plot $\nu _{\rm QL}$ as a function of pitch angle (a) and pitch angle averaged as a function of $v/v_{\text {th}e,0}$ (b). The latter shows a strong resemblance to the pitch-angle-averaged Chapman–Enskog results in figure 10: the scattering frequency shows a clear power law in $v/v_{\text {th}e,0}$ for both superthermal and subthermal velocities. The increase in scattering for subthermal velocities is due to transit-time damping and is significantly more pronounced when plotting $\nu _{\rm QL}$ as a function of $v$ versus as a function of $v_{\parallel }$ because the increased scattering at small parallel velocities applies for electrons with any $v_{\perp }$ and, thus, any $v$. The scattering frequency for superthermal electrons grows as ${\sim }(v/v_{\text {th}e,0})^{1.2-2.9}$, which is compatible with the power laws observed in the Chapman–Enskog case. It is possible that the superthermal power-law indices have a $\beta _e$ or $L_T$ dependence; however, comparing figures 10 and 13 – and indeed looking ahead to figure 23 – it is not clear that there is a consistent dependence between all three of our methods, and it seems best not to over-interpret these indices.

Figure 13.

Quasi-linear scattering frequency averaged over the range $v/v_{\text {th}e}=[2,3]$ as a function of pitch angle $\xi$ (a) and averaged over pitch angle as a function of $v/v_{\text {th}e}$ (b), the latter of which exhibits power laws in $v$ with indices from $1.2$ to $2.9$ (dotted lines).

5.3.1 Definition

The Fokker–Planck equation describes the time rate of change of a distribution function $f(t,\boldsymbol {v})$ undergoing drag, modelled by the vector $\boldsymbol {A}(\boldsymbol {v})$, and diffusion, modelled by the tensor ${{\boldsymbol{\mathsf{B}}}}(\boldsymbol {v})$. These coefficients may be calculated by tracking the Lagrangian changes in velocity,

(5.19)\begin{equation} \Delta \boldsymbol{v}(t;\Delta t) = \boldsymbol{v}(t+\Delta t,\boldsymbol{x}(t+\Delta t)) - \boldsymbol{v}(t,\boldsymbol{x}(t)) , \end{equation}

of a large number of particles having the trajectories $\boldsymbol {x}(t)$ and computing various ‘jump moments’ given by statistical averages over all of these particles within a time interval $t\in [t_{s},t_{e}]$, in which the system is assumed to be in steady state. The specific interval for each run is given in table 1. The Fokker–Planck drag vector and diffusion tensor are given by

(5.20a)\begin{equation} \boldsymbol{A}(\boldsymbol{v}) \doteq \lim_{\Delta t \rightarrow \textrm{`0'}} \frac{\langle\Delta\boldsymbol{v}\rangle}{\Delta t} \quad{\rm and}\quad{{\boldsymbol{\mathsf{B}}}}(\boldsymbol{v}) \doteq \lim_{\Delta t\rightarrow \textrm{`0'}} \frac{\langle\Delta\boldsymbol{v}\Delta\boldsymbol{v}\rangle}{\Delta t} , \end{equation}

respectively, where

(5.21)\begin{equation} \langle\cdots\rangle \doteq \frac{1}{t_{e}-t_{s}} \int^{t_{e}}_{t_{s}} \,{\rm d} t \, \int \,{\rm d}{\boldsymbol{v}} \, \left(\cdots\right) \, f(t,\boldsymbol{v}) . \end{equation}

The jump interval $\Delta t$ must be taken to be much smaller than the characteristic ‘collisional’ time scale over which $f$ evolves and yet not so small that it is comparable to the autocorrelation time of the (assumed random) kicks experienced by the particles as they interact with the electromagnetic fluctuations, i.e.

(5.22)\begin{equation} \tau_{\text{ac}} \ll \Delta t\ll \nu^{{-}1}, \end{equation}

where $\tau _{\rm ac}$ is the autocorrelation time and $\nu$ is the effective collision frequency (thus, the notation $\Delta t\rightarrow \textrm {`0'}$ rather than $\Delta t\rightarrow 0$ in (5.20a,b)). Assuming that the jumps are small, so that the consequent changes in the distribution function can be approximated by a Taylor expansion in $\Delta \boldsymbol {v}$, one can show that the Fokker–Planck operator is given by

(5.23)\begin{equation} C_{\rm FP}[f] ={-}\frac{\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot} \left[\boldsymbol{A}(\boldsymbol{v})f(t,\boldsymbol{v})\right] + \frac{1}{2}\frac{\partial^2}{\partial \boldsymbol{v}\partial\boldsymbol{v}}\,\boldsymbol{:}\,\left[ {{\boldsymbol{\mathsf{B}}}}(\boldsymbol{v})f(t,\boldsymbol{v})\right]. \end{equation}

Provided that $\boldsymbol {A}$ and ${{\boldsymbol{\mathsf{B}}}}$ can be reliably computed, (5.23) is the model collision operator for our Fokker–Planck method.

In what follows, we work in spherical velocity coordinates, $(v,\xi,\phi )$. We further assume that the interactions between the electrons and the fluctuations do not push the distribution function sufficiently far from gyrotropy (i.e. $\partial f/\partial \phi \approx 0$) to warrant retaining jump moments in $\phi$. To transform (5.23) into this coordinate system, we follow Rosenbluth, MacDonald & Judd (Reference Rosenbluth, MacDonald and Judd1957) in writing (5.23) as

(5.24) \begin{align} C_{\rm FP}[f] & ={-}\left({\mathsf{A}}^\mu f\right)_{;\mu} + \frac{1}{2}\left( {\mathsf{B}}^{\mu\nu} f \right)_{;\mu\nu} \nonumber\\ & ={-}\frac{1}{\sqrt{g}} \frac{\partial }{\partial q^\mu} \left(\sqrt{g} {\mathsf{A}}^\mu f \right) + \frac{1}{2\sqrt{g}}\frac{\partial^2}{\partial q^\mu q^\nu} \left(\sqrt{g} {\mathsf{B}}^{\mu\nu} f\right) + \frac{1}{2\sqrt{g}}\frac{\partial }{\partial q^\nu} \left(\sqrt{g} \textsf{$\mathit{\Gamma}$}_{\lambda\mu}^{\,\nu} {\mathsf{B}}^{\mu\lambda} f \right), \end{align}

where the semi-colon denotes the covariant derivative, $g\doteq \det ({\mathsf{g}}_{\mu \nu })$ is the determinant of the metric tensor ${\mathsf{g}}_{\mu \nu }$, $\textsf{$\mathit{\Gamma}$}_{\lambda\mu}^{\,\nu}$ are the associated Christoffel symbols and $q^\mu$ ranges over the coordinates $(v,\xi,\phi )$. With the metric given by

(5.25) \begin{equation} {\rm d} s^2 \doteq {\mathsf{g}}_{\mu\nu} \,{\rm d} q^\mu \,{\rm d} q^\nu = {\rm d} v^2 + \frac{v^2}{1-\xi^2} \,{\rm d} \xi^2 + v^2(1-\xi^2)\,{\rm d}\phi^2 , \end{equation}

it is then a straightforward calculation to show that (5.23) becomes

(5.26) \begin{align} C_{\rm FP}[f] & = \frac{1}{v^2} \frac{\partial }{\partial v} \left( - v^2 {\mathsf{A}}^v f + \frac{1}{2} \frac{\partial }{\partial v} v^2 {\mathsf{B}}^{vv} f \right) + \frac{\partial }{\partial \xi} \left( - {\mathsf{A}}^\xi f + \frac{1}{2} \frac{\partial }{\partial \xi} {\mathsf{B}}^{\xi\xi} f \right) \nonumber\\ & \quad + \frac{1}{v^3} \frac{\partial^2}{\partial v\partial\xi} \left( v^3 {\mathsf{B}}^{v\xi} f \right) + \frac{1}{2v^2} \left( \xi \frac{\partial }{\partial \xi} - v \frac{\partial }{\partial v} \right) \left( \frac{ v^2 }{1-\xi^2} {\mathsf{B}}^{\xi\xi} f \right), \end{align}

where the jump moments are given by

(5.27) \begin{equation} \left. \begin{gathered} {\mathsf{A}}^{v} \doteq \lim_{\Delta t \rightarrow \textrm{`0'}} \frac{\langle\Delta{v}\rangle}{\Delta t}, \quad {\mathsf{A}}^{\xi} \doteq \lim_{\Delta t \rightarrow \textrm{`0'}} \frac{\langle\Delta{\xi}\rangle}{\Delta t},\\ {\mathsf{B}}^{vv} \doteq \lim_{\Delta t\rightarrow \textrm{`0'}} \frac{\langle(\Delta{v})^2\rangle}{\Delta t} , \quad {\mathsf{B}}^{\xi\xi} \doteq \lim_{\Delta t\rightarrow \textrm{`0'}} \frac{\langle(\Delta{\xi})^2\rangle}{\Delta t}, \quad {\mathsf{B}}^{v\xi} \doteq \lim_{\Delta t\rightarrow \textrm{`0'}} \frac{\langle\Delta{v}\Delta\xi\rangle}{\Delta t} . \end{gathered} \right\} \end{equation}

Note that (5.26) can be simplified further if we demand that the operator annihilates a particular form of distribution function. For example, if we demand that $C_\textrm {FP}[f] = 0$ for an isotropic Maxwellian distribution having temperature $T=(1/2) m v^2_\textrm {th}$ and bulk velocity $v_{{w}}\hat {\boldsymbol {b}}$, then (5.26) becomes

(5.28) \begin{equation} C_{\rm FP}[f] = \frac{1}{2} \frac{\partial }{\partial \xi} {\mathsf{B}}^{\xi\xi}\left(\frac{\partial f}{\partial \xi}+v_{{w}} \frac{\partial f}{\partial v} \right) + \frac{1}{v^2}\frac{\partial }{\partial v} \frac{v^2}{v_{\text{th}}^2} {\mathsf{B}}^{vv} \left[ (v-v_{{w}}\xi) f + \frac{v_{\text{th}}^2}{2} \frac{\partial f}{\partial v} \right] . \end{equation}

Each of the three terms in this equation has a clear physical interpretation. The first term corresponds to perpendicular diffusion at fixed energy, i.e. pitch-angle scattering, occurring in a frame moving at speed $v_{{w}}$ with a velocity-dependent rate given by

(5.29) \begin{equation} \nu_{\rm FP}(v,\xi)=\frac{{\mathsf{B}}_{\xi\xi}}{1-\xi^2}. \end{equation}

This term then recovers the form of both the Chapman–Enskog (5.1) and quasi-linear (5.11) pitch-angle-scattering operators, viz.

(5.30)\begin{equation} C_{\rm FP}[f]=\frac{\partial }{\partial \xi} \left[ \frac{1-\xi^2}{2} \nu_{\rm FP}(v,\xi) \left(\frac{\partial f}{\partial \xi}+v_{{w}} \frac{\partial f}{\partial v} \right) \right]. \end{equation}

The second term in (5.28) proportional to $(v-v_{{w}}\xi )$ corresponds to drag, and the third term proportional to $\partial f/\partial v$ to energy diffusion. The physical consequence of these final two terms are that the distribution function relaxes to a Maxwellian with thermal speed $v_\textrm {th}$ and bulk velocity $v_{{w}}\hat {\boldsymbol {b}}$ at a rate given by ${\mathsf{B}}^{vv}/v_{\textrm {th}}^2$. Notwithstanding the relatively pleasant form of (5.28), we caution that the HWI is not guaranteed to push the distribution function towards an isotropic, drifting Maxwellian; indeed, a more natural expectation is that the instability pushes the system towards marginal stability, with a non-zero heat flux and its associated asymmetry in $f(v,\xi )$. In this case, (5.30), with the jump moments (5.27) calculated from the particle trajectories, is arguably a more appropriate operator than (5.28).

We calculate the expectation values $\langle \Delta v\rangle$, $\langle (\Delta v)^2\rangle$, $\langle \Delta \xi \rangle$, $\langle (\Delta \xi )^2\rangle$ and $\langle \Delta v\Delta \xi \rangle$ from the simulated particle data for a wide variety of $\Delta t$. Without a priori knowledge of the autocorrelation and collisional time scales, we select values of $\Delta t$ that are logarithmically spaced between $\varOmega _{e,0}\Delta t\in [ 0.1,500]$ and use the method described in the next subsection (§ 5.3.2) to choose the $\Delta t$ most consistent with the limit $\Delta t\rightarrow \textrm {`0'}$. As with our calculation of the distribution function, we focus on the central 4 % of the domain, which restricts our sample of particles to an area of nearly uniform temperature that is far removed from the domain boundaries. After computing indices describing the velocity-space location of each particle on a $60\times 60$ grid in $(v,\xi )$, we add the individual Lagrangian jump moment to the grid at the indices that correspond to the particle's initial phase-space location. This ensures that electrons that jump out of the 4 % domain under consideration in any given interval $\Delta t$ are included in the jump moment calculations inside the domain. Doing this for all particles in the sample, we generate an ensemble average. We repeat this process at different times throughout the interval $t\in [t_{s},t_{e}]$ (again, reported in table 1) to give a time-averaged operator, which drastically reduces the amount of noise in the computed jump moments. The accuracy of the time-averaged operators relies on the system being in steady state – i.e. in saturation – over the measurement interval. While this is not exactly satisfied for all runs, they are all very near saturation: the heat flux varies at most by $10\,\%$ across any given interval. Given the expected sampling noise, particularly out to $v/v_{\textrm {th}e,0}=3$, we take this amount of error to be acceptable.

5.3.2 Time scales and hierarchies

The limit $\Delta t\rightarrow \textrm {`0'}$ in the definitions of the jump moments (5.20a,b) ensures that the jumps are sufficiently small that the Taylor expansion that relates them to the rate of change of the distribution function converges to the true value. We estimate the autocorrelation time of the HWI fluctuations as the quasi-linear autocorrelation time $\tau _{\textrm {ac}}^{\textrm {lin}}$, the time it takes a resonant particle to interact with a wave packet (Krommes Reference Krommes2002). Given a packet of waves with central wavenumber $\bar {k}$ and width $\Delta k$,

(5.31)\begin{equation} \tau_{\text{ac}}^{\text{lin}}\sim \left(|v_{\rm p}(\bar{k}) - v_{\rm g}(\bar{k})|\Delta k\right)^{{-}1} , \end{equation}

where $v_\textrm {p}$ is the wave phase velocity and $v_\textrm {g}$ is the group velocity of the packet. From the spectra in our simulations (figure 6) we see that whistlers are energetically peaked around $k\rho _{e,0}\sim 1$ and that most of the energy is contained within a few $k\rho _{e,0}$ of the peak. Taking $k \rho _e\sim 1$ as our characteristic wavenumber and $\Delta k/k\sim 1$, we find that

(5.32)\begin{equation} \tau_{\text{ac}}^{\text{lin}}\varOmega_e\sim \left(n+1/\beta_e\right)^{{-}1}. \end{equation}

For our simulations in the limit of high $\beta _e$, the cyclotron $(n\ne 0)$ resonances will therefore have an autocorrelation time of the order of an inverse cyclotron frequency:

(5.33)\begin{equation} \tau_{\text{ac}}^{\text{lin}}(n\ne 0)\varOmega_e\sim 1. \end{equation}

However, the Landau $(n=0)$ resonance can have an arbitrarily long autocorrelation time

(5.34)\begin{equation} \tau_{\text{ac}}^{\text{lin}}(n=0)\varOmega_e\sim\beta_e \end{equation}

in the limit of high $\beta _e$. It is conceivable that in our simulations, there does not exist a $\Delta t$ for electrons with $v_{\parallel }\sim v_{{w}}$ that satisfies the ordering (5.22). We therefore rely on a model process, called an Ornstein–Uhlenbeck process, to inform our choice of an appropriate $\Delta t$ with which to calculate our drag and diffusion values.

Consider a 1-D Fokker–Planck equation in $v$ with drag to a mean velocity $\bar {v}$ at a rate $\nu$ and a diffusion coefficient $D$:

(5.35)\begin{equation} \frac{\partial f(t,v)}{\partial t}=\frac{\partial }{\partial v}\nu (v-\bar{v})f(t,v)+\frac{D}{2}\frac{\partial^2 }{\partial v^2}f(t,v). \end{equation}

This equation represents a particular statistical process called an Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein Reference Uhlenbeck and Ornstein1930) and will serve as a prototypical model for understanding our numerical results. The Green's function for the differential equation (5.35) has a well-known solution (Risken & Caugheyz Reference Risken and Caugheyz1991), which can be written in terms of the jump moments as a Gaussian distribution:

(5.36)\begin{equation} G(\Delta v,\Delta t)=\frac{1}{\sqrt{2{\rm \pi}{\rm Var}(\Delta v)}}\exp\left[-\frac{(\Delta v-\langle \Delta v\rangle)^2}{2\hspace{2pt}{\rm Var}(\Delta v)}\right], \end{equation}

with the first and second jump moments given by

(5.37a)\begin{gather} \langle \Delta v\rangle = \left(\bar{v}-v\right) \left(1-{\rm e}^{-\nu\Delta t}\right), \end{gather}

(5.37b)\begin{gather}{\rm Var}(\Delta v)=\langle \Delta v^2\rangle - \langle \Delta v\rangle^2 = \frac{D}{2\nu}\left(1-{\rm e}^{{-}2\nu\Delta t}\right), \end{gather}

for stationary solutions.

The expressions for the jump moments (5.37), taken in the limit $\Delta t\rightarrow 0$, recover the drag and diffusion coefficients in (5.35). In the opposite limit of $\nu \Delta t\gg 1$, $\langle \Delta v\rangle \simeq \bar {v}-v$ and $\textrm {Var} (\Delta v)\simeq D/\nu \sim v_{\textrm {th}}^2/2$. On this long time scale, particles drift at the mean velocity $\bar {v}$ and equilibrate to a mean temperature corresponding to $v_{\textrm {th}}$. If one were to use such a $\Delta t$ when calculating the coefficients, all kinetic information would be lost and the measured coefficients would be incorrect. We must therefore be cautious in choosing an appropriate $\Delta t$ when calculating jump moments. For cases where the jump moments are constant for a range of $\Delta t\gg \tau _{\textrm {ac}}$, the proper choice of $\Delta t$ is unambiguously in that range and the coefficients can be simply read off from the jump moment at a proper $\Delta t$. This is the case for drag and diffusion in speed, as we show first in § 5.3.3. If no such interval is present, for instance, because $\nu ^{-1}\sim \tau _{\textrm {ac}}$, the choice of an appropriate $\Delta t$ is much more complicated. One must choose the $\Delta t\ge \tau _{\textrm {ac}}$ that corresponds to the range in which the jump moments are most constant with $\Delta t$, much before the jump moments begin to scale as $\Delta t^{-1}$ and represent thermalized, non-kinetic physics. If the scattering rate is near the autocorrelation time, often the best choice is to take $\Delta t = \tau _{\textrm {ac}}$. Critically, if the scattering frequency varies in velocity space, one must repeat this process for every point in phase space one wishes to compute coefficients, or risk reporting a drag or diffusion coefficient that corresponds to thermalized physics in a portion of phase space. This is the case for drag and diffusion in pitch angle in many of our runs; we show in detail the process and results in § 5.3.4.

5.3.3 Fokker–Planck jump moments: velocity

The drag and diffusion coefficients in velocity are remarkably simple to recover from the jump moments. In figure 14 we show 2-D plots of drag and diffusion coefficients for run b40x4 at $\Delta t\varOmega _{e,0}=10$. In figure 15 we plot the values of these coefficients for all runs as a function of $\Delta t\varOmega _{e,0}$ at the point in phase space denoted by the ‘X’ in figure 14, approximately $v/v_{\textrm {th}e,0}=2.5$ and $\xi =0.45$. Across all runs, both the drag and diffusion coefficients are nearly constant as a function of $\Delta t$ for $\Delta t\varOmega _e< 25$. Runs at low $\beta _{e,0}$ and small $L_T/\rho _{e,0}$ (runs b10, b25 and b40) show a non-Markovian increase in the jump moments for $25 \le \Delta t\varOmega _{e,0}\le 500$.

Figure 14.

First (a) and second (b) velocity jump moments with $\Delta t\varOmega _{e,0}=10$ from simulation b40x4. The cross marks the point in phase space where the jump moments are plotted as functions of $\Delta t$ in figure 15 and where the probability density functions (p.d.f.s) of the jump moments are plotted for $\Delta t\varOmega _{e,0}=10$ in figure 16.

Figure 15.

Fokker–Planck drag (a) and diffusion (b) coefficients as a function of $\Delta t$ for all $\boldsymbol {\nabla } p_0=0$ simulations at the location in phase space marked by the cross in figure 14. While the moments are roughly constant for $\Delta t\varOmega _{e,0}<25$, as expected for an Ornstein–Uhlenbeck process, there is clearly some non-Markovian behaviour for larger $\Delta t$.

In the following, we assume the underlying process is Ornstein–Uhlenbeck and take as an appropriate jump interval $\Delta t\varOmega _{e,0}=10$. This assumption is well justified. The jump interval lies comfortably in the interval $\tau _{\textrm {ac}}\varOmega _{e,0}\sim 1\ll \Delta t\varOmega _{e,0}=10\ll (\nu /\varOmega _{e,0})^{-1}\sim 10^3$, and both coefficients are approximately constant with $\Delta t$ in the vicinity of $\Delta t\varOmega _{e,0}=10$ for all points in $(v,\xi )$ space. In figure 16, for all runs, we plot with solid lines histograms of individual electron jumps in speed for the jump interval $\Delta t\varOmega _{e,0}=10$ at the same location in velocity space as in figure 15. Using dashed lines we overlay Gaussian distributions with the same mean and variance as each histogram. The correspondence between the data and the Gaussian distributions further suggest that the observed process is Ornstein–Uhlenbeck.

Figure 16.

Probability densities for jump in velocity $\Delta v/v_{\textrm {th}e}$ for $\Delta t/\varOmega _e=10$ at the location in phase space denoted by the cross in figure 14. Gaussian p.d.f.s constructed from the moments of the densities are plotted in dashed lines, showing the densities are in fact Gaussian.

In figure 17 we plot normalized velocity drag ${\mathsf{A}}^v$ and diffusion ${\mathsf{B}}^{vv}$ coefficients as functions of both $v/v_{\textrm {th}e,0}$ and $\xi$. The drag coefficient is linear in $v$, with a zero near $v=v_{\textrm {th}e,0}$. The diffusion coefficient peaks at $v/v_{\textrm {th}e,0}\sim 2$, decreasing at higher and lower speeds. The dependence on $\xi$ of these coefficients, however, is more complicated. For runs b40 and b100 – i.e. those with the highest $\langle \delta B^2\rangle /B_0^2$ – drag and diffusion are greatest at $\xi =0$ and decrease towards $\xi =\pm 1$. For runs with lower fluctuation amplitudes, the dependence inverts and appears to converge towards a single shape as $L_T/\rho _{e,0}$ increases. This feature is evident in figure 14(b), where the velocity diffusion coefficient for run b40x4 peaks near $\xi =\pm 1$ for $v/v_{\textrm {th}e,0}\in [1,2]$.

Figure 17.

Drag (a,b) and diffusion (c,d) coefficients in speed as functions of $v/v_{\textrm {th}e,0}$ (a,c) and $\xi$ (b,d) for all $\boldsymbol {\nabla } p_0=0$ runs. All coefficients are approximately constant except for drag, which linearly decreases with increasing speed.

Assuming that the speed jump moments are well approximated by their pitch-angle averages, we find that the resulting ${\mathsf{A}}^v(v)$ and ${\mathsf{B}}^{vv}(v)$ can be accurately modelled by

(5.38a)\begin{equation} {\mathsf{A}}^v(v)={-}\nu_v(\beta_e,L_T)(v-v_{\text{th}e})\quad\text{and}\quad {\mathsf{B}}^{vv}(v) = D_v(\beta_e,L_T), \end{equation}

with the constant $\nu _v$ and $D_v$ plotted in figure 18 as functions of $\beta _{e,0}$ and $L_T/\rho _{e,0}$. These rates appear to be independent of $L_T/\rho _{e,0}$ but exhibit power-law dependencies on $\beta _{e,0}$ (shown by the dashed line fits). In §§ 5.1.3, 5.2.3 and 5.3.4, we show that these rates are much smaller than those associated with pitch-angle scattering. We therefore neglect ${\mathsf{A}}^v$ and ${\mathsf{B}}^{vv}$ in what follows.

Figure 18.

Drag rate $\nu _v$ (green) and diffusion coefficient $D_v$ (red) versus $\beta _{e,0}$ (a) and $L_T/\rho _{e,0}$ (b). Both $\nu _v$ and $D_v$ are nearly constant with $L_T/\rho _{e,0}$, but scale as $\beta _{e,0}^{0.53}$ and $\beta _{e,0}^{0.77}$, respectively. These values are subdominant to drag and diffusion in pitch angle (see § 5.3.4).

5.3.4 Fokker–Planck jump moments: pitch angle

Obtaining the Fokker–Planck jump moments in pitch angle, ${\mathsf{A}}^\xi$ and ${\mathsf{B}}^{\xi \xi }$, turns out to be considerably more complicated than for ${\mathsf{A}}^v$ and ${\mathsf{B}}^{vv}$, the principal reason being insufficient lack of scale separation. Analogously to figure 15, we plot in figure 19 the values of the drag and diffusion coefficients in pitch angle as functions of $\Delta t\varOmega _{e,0}$ at $(v,\xi )\simeq (2.5v_{\textrm {th}e,0},0.45)$ for all of our analysed runs. There is a stark difference between runs with higher $\langle \delta B^2\rangle /{B_0^2}$, namely b40 and b100, and those where it is lower, like b10 and b40x8. Runs with large fluctuation amplitude have jump moments ${\sim }0.1\varOmega _{e,0}$ that are comparable to the inverse autocorrelation time $\tau _{\textrm {ac}}^{-1}\sim \varOmega _{e,0}$. It is impossible to be certain whether the coefficients are constant in $\Delta t$ within such a small range. Likewise, the probability distributions of individual particle pitch-angle jumps $\Delta \xi$ are far from Gaussian; these are shown in figure 20 for the jump interval $\Delta t\varOmega _{e,0}=10$ (solid lines) and compared with Gaussian distributions having identical means and variances (dashed lines). At this jump interval, all of the distributions stray from Gaussian, with runs having larger fluctuation amplitudes exhibiting almost completely flat distributions. The distributions do, however, appear to converge towards Gaussian in the limit of low fluctuation amplitude (see b10, b40x4 and b40x8).

Figure 19.

Pitch-angle drag (a) and diffusion (b) as a function of $\Delta t$ for all $\boldsymbol {\nabla } p_0=0$ simulations at the location in phase space denoted by the cross in figure 14. These moments exhibit clear non-Markovian behaviour and in general do not conform to an Ornstein–Uhlenbeck process.

Figure 20.

Probability densities for jumps in pitch angle $\Delta \xi$ for $\Delta t\varOmega _{e,0}=10$ calculated at $(v,\xi )=(2.5,0.45)$. Gaussian p.d.f.s constructed from the moments of the densities are plotted in dashed lines; the p.d.f.s differ significantly from Gaussian due to strong scattering.

Further complicating matters, the drag and diffusion coefficients depend on $v$ and $\xi$, so in most all cases there is no single values of $\Delta t\varOmega _{e,0}$ appropriate for use in calculating these coefficients (unlike for the speed jump moments). As explained in § 5.3.2, we pick an appropriate $\Delta t$ at a subset of locations in phase space and then use this value to calculate ${\mathsf{A}}^\xi$ and ${\mathsf{B}}^{\xi \xi }$ for each $(v,\xi )$ combination. The results of this procedure for runs b40 and b40x4 are shown in figure 21. For run b40, the collisional time scale is close to the autocorrelation time across all velocity space. Run b40x4 does have some separation of time scales, particularly close to $\xi =1$ where the scattering is measured to be weak. The normalized diffusive scattering frequency $\nu _\textrm {FP}L_T/\rho _{e,0}\beta _{e,0}\varOmega _{e,0}$ resulting from the appropriate jump times are shown in figure 22. For $v/v_{\textrm {th}e,0}> 1$, there is a clear difference in symmetry about $\xi =0$ between these two runs. This asymmetry shows up in all runs with $L_T/\rho _{e,0}>250$, as shown in figure 23(a). The Fokker–Planck results at large $L_T/\rho _{e,0}$ as a function of pitch angle are structurally similar to the quasi-linear results (cf. figure 13). The ratio $\nu _\textrm {FP}(\xi =-0.83)/\nu _\textrm {FP}(\xi =0.83)\approx 6$ for run b40x8 is comparable to the quasi-linear result $\nu _\textrm {QL}(\xi =-1)/\nu _\textrm {QL}(\xi =1)\sim 10$; however, the scattering rate around $\xi =0$ is much higher for the Fokker–Planck case. Figure 23(b) shows the pitch-angle average of $\nu _\textrm {FP}$ as a function of $v$. The velocity dependence is a fast and slow power law, similar to § 5.1.3 and § 5.2.3 (cf. figures 10 and 13, respectively), which can be attributed to transit-time damping and cyclotron resonances, respectively, per the analysis in § 5.2.3. Here, as with the other models, the scattering amplitude at small velocities increases with $L_T/\rho _{e,0}$, implying that the relative contribution from transit-time damping also increases with $L_T/\rho _{e,0}$. The velocity dependence of the pitch-angle-average scattering frequency for superthermal particles here is ${\sim }(v/v_{\textrm {th}e,0})^{1.5}$, which is generally shallower than other models (see figures 10 and 13).

Figure 21.

Appropriate jump intervals $\Delta t\varOmega _{e,0}$ for a subset of $(v,\xi )$ points for all $\beta _{e,0}=40$ runs. Due to a lack of scale separation, all jump intervals for run b40 coincide with the quasi-linear autocorrelation time $\tau _{\textrm {ac}}^{\textrm {lin}}\varOmega _{e,0}\sim \Delta t\varOmega _{e,0}=1$.

Figure 22.

Normalized Fokker–Planck diffusive scattering frequency for all $\beta _{e,0}=40$ runs as a function of $v$ and $\xi$, calculated from appropriate jump times $\Delta t$, as shown in figure 21. Dashed lines correspond to contours of constant $v_{\parallel }$.

Figure 23.

Normalized $\nu _\textrm {FP}$ for all runs as a function of $\xi$ at $v/v_{\textrm {th}e,0}=2.5$ (a) and pitch angle averaged as a function of $v/v_{\textrm {th}e,0}$ (b). Runs b40–b40x8 show a transition from nearly symmetric scattering in $\xi$ to electrons with $v_{\parallel }<0$ scattering significantly faster than those with $\xi >0$. All runs show a double power law similar to § 5.1.3, with the normalized scattering frequency of the slow power law increasing with larger $L_T/\rho _{e,0}$.

5.3.5 Summary of Fokker–Planck results

In § 5.3 we calculated an effective Fokker–Planck collision operator for the HWI by constructing drag and diffusion coefficients as functions of velocity space for both velocity and pitch angle. The drag and diffusion coefficients were constructed from jump moments of tracked particles in our simulations, which we showed how to construct rigorously using jump intervals informed self-consistently by the particle statistics. The statistics for velocity were found to be consistent with Ornstein–Uhlenbeck statistics while the statistics for pitch angle were less conclusively Ornstein–Uhlenbeck due to the limited scale separation of our runs. Through a judicious choice of jump interval, we obtained the drag and diffusion coefficients and demonstrated that pitch-angle diffusion clearly dominated drag and diffusion in velocity, as is predicted by the quasi-linear model in § 5.2. The pitch-angle dependence of the pitch-angle-scattering frequency obtained by the Fokker–Planck method, however, is much flatter than the quasi-linear result (§ 5.2.3) and will prove to be a distinct advantage of the former over the latter in § 6.

6.4.1 Evidence for trapped electrons

One avenue to reduce the large contribution of low-pitch-angle electrons to the heat flux is to treat them as trapped particles that do not participate in diffusive transport. In figure 27 we show various quantities versus time for three selected electrons (blue, orange and green) from run b40x4 that exhibit periods of trapping (denoted by a background of the corresponding colour). We select time periods representative of trapping by eye, looking for large-amplitude, quasi-periodic oscillations in $\xi$. This method is by no means perfect, yet we can be confident in stating the selected electrons are largely trapped, rather than scattered: the tracks of velocity, which we have shown to be consistent with an Ornstein–Uhlenbeck process in § 5.3.3, are clearly more reminiscent of a random walk than the tracks of $\xi (t)$ and $v_{\parallel }(t)$, which exhibit periods of coherent, quasi-periodic oscillations about $0$. By sampling and averaging the period of the largest-amplitude oscillations, we estimate that $\omega _{b}/\varOmega _{e,0}\sim 0.2$. This value is approximately half what results from a naïve calculation using (6.7) after assuming that electrons with a thermal perpendicular velocity ($v_{\perp }/v_{\textrm {th}e}\sim 1$) are trapped by $k_{\parallel }\rho _e\sim 0.5$ fluctuations with root-mean-square (r.m.s.) amplitude $\langle \delta B^2\rangle /B_0^2$. Our simple expression for the bounce frequency of an electron in a single mode, (6.7), can therefore be reasonably applied to our turbulent simulations. It is more difficult to say the same of our expression for the trapped-passing boundary (6.4). For run b40x4, $\xi _\textrm {crit}\lesssim 0.5$; our tracked electrons in figure 27 show trapped oscillations largely inside the boundary, but the electrons do occasionally veer outside it. That said, what really matters is how the boundary scales with fluctuation amplitude – of which we can be confident – and the conclusions of § 6.4.2 stand regardless of its precise value.

Figure 27.

Parallel position (a), velocity (b), pitch angle (c) and parallel velocity (d) for three electrons from run b40x4 that exhibit periods of extended wave trapping, denoted by the background of the corresponding colour.

6.4.2 Heat-flux contribution from passing electrons only

As an improvement to our model, we allow only passing particles to contribute to the diffusive part of the heat flux. The trapping/passing boundary occurs at a critical pitch angle $\xi _\textrm {crit}$ given by

(6.4)\begin{equation} \xi_{\rm crit}\simeq \sqrt{\frac{\delta B/B_0}{1+\delta B/B_0}}, \end{equation}

where we have taken $\delta B/B_0=\sqrt {2}\sqrt {0.1\langle \delta B^2\rangle /B_0^2}$, i.e. the r.m.s. of the box-averaged fluctuation energy that is sampled by electrons at the $n=0$ resonance. For trapped electrons to have no contribution to the heat flux, their impact on the perturbed electron distribution function $f_{e1}$ must be neglected. We accomplish this by enforcing $f_{e1}=\textrm {const.}$ for $|\xi |\le \xi _\textrm {crit}$, requiring $f_{e1}$ to be continuous across the trapped-passing boundary $\xi =\pm \xi _\textrm {crit}$ and choose $\textrm {const}$ to satisfy the solvability conditions (5.8). We construct $f_{e1}$ subject to these conditions by simply altering the collision frequency to be infinite for trapped electrons, i.e.

(6.5)\begin{equation} \nu_{\xi_{\rm crit}}(v,\xi)=\begin{cases} \nu_{\rm model}(v,\xi), & (|\xi|> \xi_{\rm crit}),\\ \infty, & (|\xi|\le \xi_{\rm crit}). \end{cases} \end{equation}

Our approach here is similar to that of Felice & Kulsrud (Reference Felice and Kulsrud2001), but much more simple. In that work, the authors find that adiabatic mirroring regularizes the singularity in $f_1$, effectively setting the distribution to a constant for trapped particles. Their treatment included a self-consistent solution to match the distribution function between the quasi-linear and adiabatic regions; our inner and outer solutions match automatically by construction of (6.5).

We evaluate the heat flux implied by (6.5) according to (6.3) and plot the result in figure 26 using green markers. Note that the heat flux implied by this method is at or below the results for $q_{v_\textrm {w}}$, indicating that the diffusive contribution to the heat flux is zero or negative. The reason the diffusive heat flux is negative comes down to how setting $f_{e1}=\textrm {const}$ inside $|\xi |<\xi _\textrm {crit}$ affects the integrand of (6.3). Without any modification to $f_{e1}$ calculated using (6.1c), the integrand is dominated by values where $1/\nu$ is the largest; these regions are positive, hence, the positive heat flux. When we set $f_{e1}=\textrm {const}$ for $|\xi |<\xi _\textrm {crit}$, the regions of largest $1/\nu$ no longer contribute to the integrand; the positive contribution is therefore suppressed relative to the negative one and the overall diffusive heat flux becomes negative. Perhaps a more self-consistent treatment like Felice & Kulsrud (Reference Felice and Kulsrud2001) is required when applying this method, as it clearly does not produce physical results for all physically motivated collision frequencies.

6.4.3 Refined resonance condition at $|v_{\parallel }/v_{\textrm {th}e}|\ll 1$

With the cyclotron resonance condition (5.16) requiring $k_{\parallel }\rho _e\rightarrow \infty$ as $|v_{\parallel }/v_{\textrm {th}e}|\rightarrow 0$, it ought to be checked whether or not the dispersion relation used to derive this result is in fact still appropriate as $k_{\parallel }\rho _e\rightarrow \infty$. As it turns out, (5.13) breaks down for ${|v_{\parallel }/v_{\textrm {th}e}|\lesssim \beta _e^{-1/2}}$. The breakdown occurs for two reasons: firstly, the term under the square root in (5.14) is no longer small; secondly, and most importantly, $k_{\parallel }\rho _e\rightarrow \beta _e^{1/2}$, in which case $\omega (k_{\parallel }\rho _e)\rightarrow \varOmega _e$ under (2.3). It is clear from (2.5) that $\omega \rightarrow \varOmega _e$ as $|v_{\parallel }/v_{\textrm {th}e}|\rightarrow 0$, even without any assumptions on the size of $\beta _e$. However, the dispersion relation (2.3) used in (5.13) is only valid for $k_{\parallel }\rho _e\ll \beta _e^{1/2}$. One could use the full parallel cold plasma dispersion relation (Stix Reference Stix1992),

(6.6)\begin{equation} \frac{\omega(k_{{\parallel}}\rho_e)}{\varOmega_e}=\frac{(k_{{\parallel}}\rho_e)^2/\beta_e}{1+(k_{{\parallel}}\rho_e)^2/\beta_e}, \end{equation}

which would lead to a different resonance condition for $|v_{\parallel }/v_{\textrm {th}e}|\le \beta _e^{-1/2}$. However, even (6.6) is liable to be wrong in detail due to thermal corrections. Thus, a detailed understanding of the dispersion relation is necessary to understand the resonant wave vector at small parallel velocities.

6.4.4 Resonance broadening

One additional possibility for whistler waves to scatter electrons self-consistently across the gap is resonance broadening. The delta function that features in the quasi-linear operator is customarily replaced with a peaked function having finite width (Berk et al. Reference Berk, Breizman, Fitzpatrick and Wong1995). Resonance broadening can result from particle trapping or deviation from linear orbits resulting from statistical variation in the integrated trajectories when calculating the quasi-linear operator. In the former case, the resonance broadening function is often taken to be a table-top distribution of width $\Delta \varOmega =4\omega _{b}$ (Karimabadi, Krauss-Varban & Terasawa Reference Karimabadi, Krauss-Varban and Terasawa1992), where (Roberg-Clark et al. Reference Roberg-Clark, Drake, Reynolds and Swisdak2016; Cai, Wu & Tao Reference Cai, Wu and Tao2020)

(6.7)\begin{equation} \frac{\omega_{b}(k\rho_e)}{\varOmega_e}=\sqrt{k\rho_e\frac{v_{{\perp}}}{v_{\text{th}e}} \frac{\delta B(k\rho_e)}{B_0}} \end{equation}

is the bounce frequency for a deeply trapped particle under the pendulum approximation with perpendicular velocity $v_\perp$ bouncing in a magnetic mirror of amplitude $\delta B$ and wavneumber $k$, viz.

(6.8)\begin{equation} R_{\rm K92}(x/\varOmega_e,\omega_{b}/\varOmega_e)= \begin{cases} \dfrac{\varOmega_e}{4\omega_{b}}, & |x|\le 2\omega_{b},\\ 0, & |x|>2\omega_{b}. \end{cases} \end{equation}

Here, $x=\omega (k_{\parallel })-k_{\parallel } v_{\parallel }+n\varOmega _e$ is the exact linear resonance condition. A recent numerical study of resonance broadening for a single wave mode (Meng et al. Reference Meng, Gorelenkov, Duarte, Berk, White and Wang2018) has verified that the broadening width $\Delta \varOmega =4\omega _{b}$ is valid at small wave amplitudes.

In the case of resonance broadening by statistical variation in the linear particle orbits, one obtains a resonance function (Dupree Reference Dupree1966)

(6.9)\begin{equation} R_{D}= \frac{1}{\rm \pi}\int_0^\infty\,{\rm d}\tau\, \exp\left({{\rm i}(\omega(k_{{\parallel}})-k_{{\parallel}}\langle v_{{\parallel}}\rangle -n\varOmega_e)\tau-\frac{1}{3}k_{{\parallel}}^2\nu(v_{{\parallel}},v_{{\perp}})v_{\text{th}e}^2\tau^3}\right). \end{equation}

Because the scattering frequency appears as a parameter in (6.9), one is left with an implicit expression for $\nu$ when solving an expression like (6.1c) with the delta function replaced with (6.9). In either the case of (6.8) or (6.9), because the scattering frequency or fluctuation energy controls the width of the broadened resonance, the resulting scattering frequency does not in general scale with the energy in the fluctuations. The heat flux implied by a resonance-broadened operator calculated from the magnetic energy spectrum observed in our simulations, therefore, does not equal the HWI threshold heat flux, except perhaps in special circumstances, casting doubt on whether current theories of resonance broadening can explain our results. This is not to say, however, that there could not be a physically relevant regime in which resonance broadening might control HWI saturation and in which the scattering rate scales non-diffusively with the magnetic energy – we just have yet to see it.

Closely related to, and intertwined with the concept of resonance broadening, is the concept of resonance overlap. For a discrete wave mode of infinitesimal amplitude, resonant diffusion occurs along quasi-linear contours in velocity space within an infinitesimal neighbourhood of the resonant parallel velocity. As the wave amplitude is increased, particles within a finite neighbourhood of the resonance become trapped. Resonances are broadened by the trapping according to, e.g. (6.8) with (6.7), leading to diffusion with finite extent along the quasi-linear contours (Karimabadi et al. Reference Karimabadi, Krauss-Varban and Terasawa1992). Particles can therefore only diffuse across discrete wave modes when the wave amplitudes are large enough that their resonances overlap (Chirikov Reference Chirikov1960). Earlier work on the HWI by Roberg-Clark et al. (Reference Roberg-Clark, Drake, Reynolds and Swisdak2016) concluded that the overlap of Landau ($n=0$) and cyclotron ($n=\pm 1$) resonances was necessary to explain the saturation of the HWI in two dimensions. Overlap was found to occur for wave amplitudes $\delta B^2/B_0^2\gtrsim 0.09$; all of our simulations saturate in this regime. Whether or not the HWI saturates in the way we have observed in this work at wave amplitudes below this overlap threshold remains to be seen. Our analysis, however, suggests a more strict condition for HWI saturation than Roberg-Clark et al. (Reference Roberg-Clark, Drake, Reynolds and Swisdak2016), namely, that resonances not only need to overlap, but the scattering frequency where they overlap must also scale like the marginal stability criterion (2.6). We present a model that satisfies this condition in § 6.4.5.

6.4.5 Semi-empirical scattering frequency

While we are so far unable to affix an analytical theory to our quasi-linear model that satisfactorily scatters electrons across the $|v_{\parallel }/v_{\textrm {th}e}|\ll 1$ gap, our results clearly indicate that they are. Not only do all of our simulations saturate at the threshold heat flux (§ 4.2), but our Fokker–Planck analysis (§ 5.3.4) directly reports a finite pitch-angle-scattering rate across the gap that is also consistent with the threshold heat flux (§ 6.3). Motivated by these observations and eschewing any detailed theoretical arguments, we propose a semi-empirical model that is identical to (6.1) except that the scattering frequency (6.1c) has a floor that scales with $\beta _e\rho _e/L_T$, i.e.

(6.10)\begin{equation} \frac{\nu_{{\rm SE}}(v,\xi)}{\varOmega_e}=\max \left[\frac{\nu_{\rm model}(v,\xi)}{\varOmega_e}, \frac{1}{2}\frac{\beta_e\rho_e}{L_T}\right]. \end{equation}

The heat fluxes implied by (6.10) are shown by the grey dots in figure 26 and closely agree with the heat fluxes measured in our simulations. While we cannot guarantee the asymptotic relevance of this operator, it is at least consistent out to the scale separations we were able to simulate. One can also argue that the scattering rate must be this in order to saturate self-consistently under the diffusive scaling predicted by theory and observed in this paper. The fact that our simulations have extended out to larger scale separations numerically than previously published ones and yet the diffusive scaling still holds lends us a certain degree of confidence that the scaling might continue to astrophysical scale separations. We discuss the possibility that the asymptotic HWI does not saturate in this manner in § 6.5.