Collisionless conduction in a high-beta plasma: a collision operator for whistler turbulence

Abstract

The regulation of electron heat transport in high-$\beta$, weakly collisional, magnetized plasma is investigated. A temperature gradient oriented along a mean magnetic field can induce a kinetic heat-flux-driven whistler instability (HWI), which back-reacts on the transport by scattering electrons and impeding their flow. Previous analytical and numerical studies have shown that the heat flux for the saturated HWI scales as $\beta _e^{-1}$. These numerical studies, however, had limited scale separation and consequently large fluctuation amplitudes, which calls into question their relevance at astrophysical scales. To this end, we perform a series of particle-in-cell simulations of the HWI across a range of $\beta _e$ and temperature-gradient length scales under two different physical set-ups. The saturated heat flux in all of our simulations follows the expected $\beta _e^{-1}$ scaling, supporting the robustness of the result. We also use our simulation results to develop and implement several methods to construct an effective collision operator for whistler turbulence. The results point to an issue with the standard quasi-linear explanation of HWI saturation, which is analogous to the well-known $90^{\circ }$ scattering problem in the cosmic-ray community. Despite this limitation, the methods developed here can serve as a blueprint for future work seeking to characterize the effective collisionality caused by kinetic instabilities.

1 Introduction

1.1 Physical motivation

1.1.1 The intracluster medium of galaxy clusters

Galaxy clusters are the largest virialized structures in the Universe, spanning more than a megaparsec in diameter. Clusters are bound by dark matter, which forms a deep gravitational well. Dark matter accounts for the majority of cluster mass, ${\approx }84\,\%$. The remaining mass is attributed to baryonic matter in the hundreds to thousands of galaxies scattered across the cluster (${\approx }3\,\%$) or the relatively hot, diffuse plasma between them (${\approx }13\,\%$) known as the intracluster medium (ICM) (see Peterson & Fabian 2006 for a review). Studies of temperature gradients in various clusters (e.g. Binney & Cowie 1981; Ettori & Fabian 2000; Vikhlinin, Markevitch & Murray 2001; Markevitch et al. 2003; Zakamska & Narayan 2003) have suggested that heat conduction must be suppressed by a multiplicative factor of ${\sim }0.3$ to ${\sim }10^{-2}$ relative to the Spitzer value (Fabian, Voigt & Morris 2002), which, because of its strong temperature dependence ${\propto }T^{7/2}$ (Spitzer 1962), would otherwise isothermalize the ICM rapidly or fail to prevent runaway cooling (e.g. Bregman & David 1988; Kim & Narayan 2003; Conroy & Ostriker 2008).

1.1.2 The solar wind

The modern understanding of the solar wind has its origins in seminal work by Parker (1958): the corona, which is heated to temperatures upwards of $10^6\,{\rm K}$ (e.g. by Alfvén-wave turbulence (Chandran & Hollweg 2009; Chen 2022; Squire et al. 2022) or perhaps interchange reconnection (Raouafi et al. 2023)), provides much of the free energy and mass required to accelerate the plasma to velocities exceeding ${\approx }400\,{\rm km}\,{\rm s}^{-1}$ (see Verscharen, Klein & Maruca 2019, for a review). The observed radial electron temperature profile inside a few au does not match predictions from a purely adiabatic wind (Richardson & Smith 2003); heat fluxes are therefore thought to play an important role in the thermodynamic evolution of the solar wind. A statistical study of Wind spacecraft measurements taken from the $\beta \gtrsim 1$, weakly collisional slow wind (Bale et al. 2013) revealed that the scaling of the electron heat flux with the Coulomb-collisional mean free path $\lambda _{\text {mfp},e}$ depends on the local temperature-gradient length scale $L_T$: for $L_T\gtrsim \lambda _{\text {mfp},e}$, the measured heat flux matched collisional predictions (Spitzer 1962), whereas for $L_T/\lambda _{\text {mfp},e} \lesssim 3$, the heat flux was found to be constant with $L_T$. Similar results were also found in data from the Parker Solar Probe (Halekas et al. 2021).

1.1.3 Low-luminosity black-hole accretion flows

Accretion flows onto supermassive black holes are often significantly underluminous when compared with predictions from classical thin-disc theory. For example, if plasma were accreted at the Bondi (1952) rate onto the ${\approx }4\times 10^6\,{\rm M}_\odot$ black hole at the Galactic centre, Sgr A$^\ast$, via a geometrically thin, optically thick accretion disc (Shakura & Sunyaev 1973), then the luminosity would be ${\sim }10^5$ times larger than that observed (see reviews by Quataert 2003; Yuan & Narayan 2014). Such radiatively inefficient accretion flows (RIAFs) can be explained by a combination of substantially sub-Bondi accretion, because much of the inflowing plasma is gravitationally unbound and lost to a magnetically driven wind (e.g. Blandford & Begelman 1999; Hawley & Balbus 2002) and low radiative efficiency, because the liberated gravitational potential energy is stored as thermal energy primarily in the poorly radiating ion population (e.g. Narayan & Yi 1994). The result is a geometrically thick accretion flow, one in which the thermodynamics of the putative collisionless, high-$\beta$ plasma should play an important role, whether by instigating hydrodynamic (Narayan, Igumenshchev & Abramowicz 2000; Quataert & Gruzinov 2000) or magnetothermal (Balbus 2001) convective transport, or by being modified directly by conductive transport (Tanaka & Menou 2006; Johnson & Quataert 2007; Ressler et al. 2015).

1.2 General transport considerations and history

Given the potential importance of conductive and convective heat transport in the ICM, solar wind and RIAFs, it is not surprising that it remains an active area of research within each of the associated communities. Obtaining definitive physical models for this transport, however, is made difficult by the extreme scale separations characterizing these plasmas and by the complexities such multiscale physics brings. Namely, all of these plasmas are weakly collisional, with the ratio of Coulomb mean free path and macroscopic temperature-gradient length scales ranging from ${\sim }10^{-2}$ to ${\sim }1$ in the ICM, in the solar wind and around Sgr A$^\ast$ near the Bondi radius (increasing to values ${\gg }1$ towards the event horizon). They are also highly magnetized, with the ratio of electron gyroradius $\rho _e$ to $L_T$ being (very roughly) ${\sim }10^{-15}$ in the ICM, ${\sim }10^{-8}$ in the solar wind and ${\sim }10^{-12}$ at the Bondi radius of Sgr A$^\ast$. As a result, the motions of charged particles are tightly bound to magnetic-field lines and the transport of both momentum and heat perpendicular to the field is highly suppressed relative to the parallel transport. Thus, both the angle of the magnetic field with respect to the local temperature gradient (e.g. in insulated cluster cold fronts; Vikhlinin et al. 2001; Markevitch & Vikhlinin 2007) and the extent to which the field lines are tangled by turbulence (e.g. Chandran & Cowley 1998; Narayan & Medvedev 2001) have a significant effect on heat transport. If the temperature gradient is aligned or anti-aligned with a confining gravitational field, even energetically weak magnetic fields can drive buoyancy instabilities like the magnetothermal instability (Balbus 2000, 2001) or heat-flux buoyancy instability (Quataert 2008), provided that the conductive heat transport between fluid elements is rapid and restricted along field lines (see Kunz (2011) and Xu & Kunz (2016) for detailed treatments of these instabilities for weakly collisional and collisionless plasmas, respectively).

Kinetic instabilities can also strongly affect transport in weakly collisional, magnetized plasma. The transport of heat or momentum implies a distortion of the distribution function (e.g. Braginskii 1965), which can be a source of free energy for kinetic instabilities (e.g. Bott, Cowley & Schekochihin 2024). If the plasma beta $\beta =8{\rm \pi} p/B^2$ (the ratio of thermal pressure $p$ and magnetic pressure $B^2/8{\rm \pi}$) is large, small departures from local thermodynamic equilibrium, and thus, small amounts of transport implied by the departure, are enough to grow Larmor-scale distortions in the energetically weak magnetic field. Such distortions can scatter or trap particles, pushing the particles’ velocity distribution function back towards isotropy and thereby limiting transport.

One instability of note in this context is the heat-flux-driven whistler instability (HWI), which was first investigated in the case of a weakly collisional plasma by Levinson & Eichler (1992). Those authors derived a growth rate for the case of a wave vector oriented parallel to the local magnetic field; however, they considered the case where the HWI is saturated by wave-mode coupling and found that the saturated heat flux was independent of the quasi-linear scattering rate. A few years later, Pistinner & Eichler (1998) found that oblique whistler waves, which are elliptically (rather than circularly) polarized, could diffusively scatter heat-flux-carrying electrons via cyclotron resonance and limit the heat flux $\propto \beta _e^{-1}$. Following these analytical results has been a relatively recent flurry of numerical works, namely particle-in-cell (PIC) simulations of the HWI. The first of these were one-dimensional (1-D) simulations (Roberg-Clark et al. 2016), which found little reduction of the heat flux by parallel whistlers, confirming the need for oblique waves. It was not until two-dimensional (2-D) PIC simulations were performed that the predicted ${\sim }\beta _e^{-1}$ scaling was shown empirically. These runs were performed concurrently by two independent groups: Komarov et al. (2018) (hereafter K18) and Roberg-Clark et al. (2018) (hereafter RC18).

The literature for heat-flux instabilities in the solar wind evolved largely independently from the more astrophysics-focused work discussed in the preceding paragraph. This is likely due to the emphasis on high-$\beta _e$ plasma for the latter and the multiple electron populations encountered in the former. The ‘whistler heat-flux instability’, or WHFI as it is called in the solar wind literature, was first investigated in Gary et al. (1975) for the case of two bi-Maxwellian electron populations representing the core and the halo. The instability was presented as a viable mechanism for heat-flux reduction and was supported by a semi-empirical model by Gary et al. (1994), who used data from the Ulysses spacecraft to show that the parallel heat flux is suppressed by a factor of $\beta _{\parallel, {c}}^{-0.9}$, where $\beta _{\parallel, {c}}$ is the parallel beta of the electron core population. This scaling was later supported by analytic calculations in the limit of high $\beta _{\parallel, {c}}$ (Gary & Li 2000). A statistical study by Scime et al. (1994) also reported that the WHFI best explained the heat flux measured by Ulysses. A later study using data from the Artemis mission (Tong et al. 2019) revealed correlations between the properties of whistler waves seen in the magnetic-field data and the macroscopic plasma parameters, in a way that is consistent with their excitation by the HWI.

1.3 Purpose and organization

Previous numerical work on the HWI suffers from a number of limitations. First, the equilibria adopted in these works – isobaric in the case of K18 and collisionless in the case of RC18 – are unlikely to be found in actual galaxy clusters. Secondly, the scale separations used in the numerical simulations are limited, with the ratio of temperature-gradient length scale to electron gyroradius being no more than 256. As a result, the saturated fluctuation amplitude of the whistlers, predicted to satisfy $\delta B/B_0 \sim (\beta _e \rho _e/L_T)^{1/2}$, approached the strength of the background field. On these topics, the aim of this paper is to assess the extent to which prior results on the whistler-mediated heat flux by K18 and RC18 are robust with respect to the physical set-up and scale separation. For numerical results of the HWI to be extrapolated reliably to astrophysical scale separations and more realistic astrophysical environments, it ought to be checked that the underlying physics of the saturated instability converges at large scale separation and that the HWI is robust to different physical set-ups. We therefore perform electromagnetic PIC simulations similar to those in RC18 and K18, but with an additional, physically motivated equilibrium (i.e. in which thermal stratification is associated with hydrostatic equilibrium in a gravitational field) and a focus on studying the largest scale separations available to us numerically.

We then use the results of our simulations to obtain a numerical effective collision operator for the saturated HWI using three distinct methods: the first leverages a Chapman–Enskog expansion to calculate a pitch-angle-scattering frequency from the electron distribution function, the second uses a quasi-linear operator to obtain a pitch-angle-scattering rate from the magnetic spectrum, and the third utilizes a Fokker–Planck method to obtain a pitch-angle-scattering operator from tracked particle data. The results of each of these methods will be discussed and synthesized in the context of building a model effective collision operator for the HWI. Our model operator has inherent limitations associated with it, which we show to be a manifestation of a longstanding issue with quasi-linear analyses. Despite our inability to construct a fully self-consistent model operator, we hope the present work will advance the methods available and that our experience will provide a useful road map so that future studies of the saturation of kinetic instabilities may be more successful.

The paper is organized as follows. In § 2 we provide a brief quantitative introduction to the HWI. We then detail our numerical methods and simulation diagnostics in § 3. Section 4 contains simulation results that largely confirm past work – namely, the dependence of the saturated heat flux and wave amplitude on $\beta _e$ and the temperature-gradient length scale. For those already familiar with the instability, the heat flux and wave amplitude dependencies are given by (4.4) and are consistent with an effective collision frequency given by (2.6). All of these relationships hold regardless of the equilibrium state and agree with previous work. Section 5 contains the bulk of our original analysis. Here, we detail three methods for obtaining model collision operators from our simulations, namely: leveraging a Chapman–Enskog expansion (§ 5.1), a quasi-linear method (§ 5.2) and a Fokker–Planck method (§ 5.3). In § 6 we present a physically motivated model collision operator (§ 6.2), discuss its shortcomings in explaining the heat flux observed in our simulations (§ 6.3) and end the section by exploring various ways to deal with those shortcomings (§��6.4). Finally, we present our conclusions in § 7.

2 Background

Given the existing body of work on the HWI, we present here only a brief overview of the instability and refer the reader to the works cited in the introduction (§ 1) for more detail. In § 2.1 we provide a simplified description of the HWI, focusing mostly on its qualitative features. It should be noted that the quantitative details used to support this description differ slightly from those obtained from the fully self-consistent, collisionless simulations presented in subsequent sections. A reader already familiar with the HWI can skip to § 3 for our numerical methods and set-up, § 4 for a summary of our simulation results or § 5 for our model collision operator results.

2.1 The HWI

The HWI is a destabilization of the whistler branch of the plasma dispersion relation in the presence of a heat flux. To see this, consider a plasma threaded by a magnetic field $\boldsymbol {B}$. If the plasma is magnetized, meaning that the electron gyroradius $\rho _e$ is much smaller than some characteristic macroscale length scale $L$ in the plasma, then the transport of heat and momentum will be highly anisotropic, with transport along (‘$\parallel$’) the magnetic field occurring much more rapidly than transport across (‘$\perp$’) the field. For example, if a macroscopic electron temperature gradient $\boldsymbol {\nabla } T_e$ has a component oriented along the magnetic-field direction $\hat {\boldsymbol {b}}=\boldsymbol {B}/|\boldsymbol {B}|$ (i.e. $\nabla _\parallel T_e \ne 0$), then electrons will flow to produce a heat flux $\boldsymbol {q}_e$ that is oriented predominantly along the field, viz. $\boldsymbol {q}_e\simeq q_{\parallel,e}\hat {\boldsymbol {b}}$. If, furthermore, the collisional mean free path $\lambda _{\text {mfp}}$ is much smaller than the temperature-gradient length scale along the field, $L_T=(-\nabla _\parallel \ln T_e)^{-1}$, then the distortion in the velocity distribution function of the electrons associated with the heat flux will be small, ${\sim }\lambda _{\text {mfp}}/L_T$, relative to the equilibrium (Maxwellian) distribution. Assuming that electrons are pitch-angle scattered at some velocity-independent rate $\nu$, these considerations suggest a steady-state electron distribution function given by

(2.1)\begin{equation} f_e(v,\xi)=\frac{n_e}{{\rm \pi}^{3/2}v_{\text{th}e}^3}\,{\rm e}^{{-}v^2/v_{\text{th}e}^2}\left[1-\frac{v_{\text{th}e}}{\nu L_T}\, \xi\,\frac{v}{v_{\text{th}e}} \left(\frac{v^2}{v_{\text{th}e}^2}-\frac{5}{2}\right)\right], \end{equation}

where $\xi \doteq v_{\parallel }/v$ is the cosine of the pitch angle, $n_e$ is the local electron number density and $v_{\text {th}e}$ is the local electron thermal speed (e.g. Braginskii 1965). For a derivation of a similar expression, see Appendix A.

Parallel-propagating whistler waves in a plasma with an electron distribution function given by (2.1) are unstable, with a growth rate dependent upon the product of the mean free path normalized by the temperature-gradient length scale and the electron plasma beta parameter $\beta _e=8{\rm \pi} p_e/B^2$ – the ratio of scalar electron thermal pressure $p_e$ to magnetic pressure $B^2/8{\rm \pi}$. An analytic expression for the growth rate of whistlers with parallel wavenumber $k_{\parallel }=\boldsymbol {k}\boldsymbol {\cdot } \hat {\boldsymbol {b}}>0$ may be obtained in the asymptotic limit $v_{\text {th}e}\beta _e/\nu L_T\ll 1$ and $k_{\parallel }\rho _e/\beta _e^{1/2}\ll 1$ (see § 3.3.1 of Bott et al. 2024); it is

(2.2)\begin{equation} \frac{\gamma}{\varOmega_e}\simeq\frac{\sqrt{\rm \pi}}{k_{{\parallel}}^2\rho_e^2}\left(\frac{v_{\text{th}e}}{2\nu L_T}-\frac{k_{{\parallel}}^3\rho_e^3}{\beta_e}\right)\exp\left(-\frac{1}{k_{{\parallel}}^2\rho_e^2}\right), \end{equation}

where $\varOmega _e=v_{\text {th}e}/\rho _e$ is the (positive) electron gyrofrequency. (The qualitative features of (2.2) do not change outside the limit $v_{\text {th}e}\beta _e/\nu L_T\ll 1$; the same is not true for the $k_{\parallel }\rho _e/\beta _e^{1/2}\ll 1$ limit.) Note that the growth rate is only positive when $L_T>0$, i.e. for whistlers moving down the temperature gradient. The corresponding dispersion relation for the (real) wave frequency $\omega$ is

(2.3)\begin{equation} \frac{\omega}{\varOmega_e} \simeq \frac{kk_{{\parallel}} \rho_e^2}{\beta_e}. \end{equation}

Coincidentally, this is the same as the cold plasma dispersion relation for whistler waves, $\omega /\varOmega _e \simeq kk_\parallel d^2_e$, where $d_e$ is the electron skin depth.

The growth rate (2.2) and dispersion relation (2.3) are plotted schematically in figure 1. The growth rate has a maximum near $k_{\parallel }\rho _e\sim 1$ (while this statement does not strictly follow from (2.2) because the expression was derived in a certain asymptotic limit, the statement is rigorous in the correct limit; see Bott et al. (2024) for specific details). For this wave vector, the whistler phase velocity

(2.4)\begin{equation} v_{{w}}\sim\frac{\omega}{k_{{\parallel}}}\sim \frac{v_{\text{th}e}}{\beta_e} \end{equation}

is very subthermal. Resonant wave–particle interactions occur when the frequency of the wave vanishes in the frame co-moving with a particle having parallel velocity $v_{\parallel }$, i.e.

(2.5)\begin{equation} \omega-k_{{\parallel}} v_{{\parallel}}+n\varOmega_e=0, \end{equation}

where $n$ is an integer. Whistler waves in the HWI, therefore, are only Landau $(n=0)$ resonant with electrons that are buried deep in the core of the distribution and do not contribute to the instability. Instead, the principal cyclotron $(n=\pm 1)$ resonances, which couple the waves with thermal $(v_{\parallel }\sim v_{\text {th}e})$ electrons in the unstable region of the distribution function (2.1), mediate the instability. This is represented schematically in figure 1, where the cyclotron resonances corresponding to thermal electrons intersect the whistler dispersion relation at the $k_{\parallel }\rho _e$ corresponding to maximum growth rate. The decline in growth rate at small $k_{\parallel }\rho _e$ can be identified with electrons that have superthermal parallel velocities, of which there are exponentially few (thus, the multiplicative factor of $\exp (-1/k_{\parallel }^2\rho _e^2$)). The growth rate has a zero at $k_{\parallel }\rho _e=(\beta _ev_{\text {th}e}/2\nu L_T)^{1/3}$; any waves with larger $k_{\parallel }$ are cyclotron damped on subthermal electrons in the core of the distribution.

Figure 1. The HWI growth rate (red line), whistler dispersion relation (solid black line) and cyclotron resonances (dashed lines) plotted as functions of $k_{\parallel }\rho _e$. For thermal electrons, only the cyclotron resonances intersect the whistler dispersion relation at values of $k_{\parallel }\rho _e$ corresponding to the maximum growth rate.

2.2 The saturation of the HWI

As the whistlers grow to finite amplitude, resonant wave–particle interactions satisfying (2.5) begin to diffuse electrons, which pushes the distribution function toward isotropy. This in turn reduces the whistler growth rate. The system eventually reaches saturation, a steady state in which the waves have large enough amplitudes to scatter electrons and thereby hold the distribution function at a marginally stable state in which the growth rate is zero. Denoting the rate at which the nonlinear whistlers diffuse electrons as $\nu _{\rm eff}$, an ‘effective collisionality’, and substituting $\nu _{\rm eff}$ for $\nu$ in (2.2) with $k_{\parallel }\rho _e\sim 1$, we find that the HWI should saturate when

(2.6)\begin{equation} \nu_{\rm eff} \sim \frac{\beta_ev_{\text{th}e}}{L_T} . \end{equation}

The saturated heat flux in the saturated state should then scale as

(2.7)\begin{equation} q^{\rm sat}_{{\parallel}, e }\sim m_en_e\frac{v_{\text{th}e}^4}{\nu_{\rm eff} L_T}\sim m_en_e\frac{v_{\text{th}e}^3}{\beta_e}. \end{equation}

If whistlers do indeed scatter electrons diffusively, the effective scattering frequency can be related to the whistler fluctuation energy via

(2.8)\begin{equation} \nu_{\rm eff}\sim \varOmega_e\frac{\delta B^2}{B_0^2} . \end{equation}

Matching this expression for $\nu _{\rm eff}$ with (2.6) implies that

(2.9)\begin{equation} \frac{\delta B^2}{B_0^2}\sim \frac{\beta_e}{L_T/\rho_e} \end{equation}

in saturation.

The scaling (2.9) was predicted by K18 and obtained empirically using bespoke kinetic PIC simulations. Those authors then argued that the whistlers pitch-angle scatter electrons in the frame of the background flow. Roberg-Clark et al. (2018) argued instead that whistlers advect thermal energy at their phase velocity, so that

(2.10)\begin{equation} q_{{\parallel},e}\sim n_eT_e\langle v_{{\parallel}}\rangle\sim n_eT_ev_{{w}} \sim m_en_e\frac{v_{\text{th}e}^3}{\beta_e}. \end{equation}

Both arguments reproduce the scaling (2.7) for the saturated heat flux. This was recognized by Drake et al. (2021), who concluded that whistlers should act as pitch-angle scatterers in a frame moving at the wave phase velocity. The argument for electromagnetic waves scattering electrons in this way goes back to Kennel & Engelmann (1966), who derived a quasi-linear wave–particle diffusion operator and found that, for particles with $v_{\parallel } \gg \omega /k_{\parallel }$, the gyroresonant scattering is dominated by scattering in pitch angle due to the parallel electric field vanishing in the frame of the particle; only particles with $v_{\parallel } \sim \omega /k_{\parallel }$ have strong velocity diffusion.

Another aspect agreed upon by K18, RC18 and Drake et al. (2021) is the need for oblique – and, thus, elliptically polarized – whistler waves to achieve the heat-flux scaling $q_{\parallel e}\sim \beta _e^{-1}$. This requirement, initially discovered by Pistinner & Eichler (1998), is due to both the right-handedness of purely parallel whistler waves and the fact that the HWI only drives unstable whistlers travelling down the temperature gradient. Any electron travelling down the temperature gradient faster than the whistler phase speed will therefore never be in cyclotron resonance with the wave, which is Doppler shifted to a left-hand polarization in that electron's frame. In contrast, the elliptical polarization of oblique waves can be decomposed into left- and right-handed components, which can be in cyclotron resonance with electrons travelling in either direction along the temperature gradient. What seems to have gone unappreciated, however, is that the ratio of energy in left- to right-handed waves is significantly below unity for all but the most oblique wave vectors. Using the cold plasma whistler dispersion (Stix 1992), the ratio of energies in the left-handed ($-$) and right-handed ($+$) components is

(2.11)\begin{equation} \frac{\mathcal{E}^-}{\mathcal{E}^+}\simeq \left(\frac{\cos\theta -1}{\cos\theta+1}\right)^2 . \end{equation}

Even for oblique waves with $\theta = 45^\circ$, $\mathcal {E}^-/\mathcal {E}^+\simeq 0.03$. (Finite-temperature effects can increase this ratio.) Given (2.8), it is possible an asymmetry will arise in the scattering rate with respect to pitch angle; as far as we are aware, this effect has not been commented on in the literature thus far.

The ion species can also affect the saturation of the HWI. At high $\beta _e$, ions may be resonant with oblique whistler waves. For example, when $\beta _e=(m_i/m_e)^{1/2}$, the ion thermal speed is concurrent with the cold plasma whistler phase speed, viz. $v_{\text {th}i}=(m_i/m_e)^{-1/2}v_{\text {th}e}=v_{\text {th}e}/\beta _e\sim v_{{w}}$, and the waves can Landau damp on the ions. For $\beta _e=m_i/m_e$, the whistler wave frequency is of the order of the ion gyrofrequency, viz. $\omega \sim \varOmega _e/\beta _e=\varOmega _i$, and the ions can cyclotron damp the waves. Such interactions would heat the ion species at the expense of the strength of the turbulent fluctuations and, therefore, might increase the saturated electron heat flux in the HWI.

2.3 Constructing collision operators

In this work we explore three methods to construct a collision operator from our kinetic simulations describing how the electrons interact with the HWI fluctuations: Chapman–Enskog, quasi-linear and Fokker–Planck. In § 5.1 we leverage measurements of the distribution function in the simulations to solve the Chapman–Enskog problem for a postulated pitch-angle-scattering operator and obtain the associated scattering frequency $\nu _{\rm CE}$ as a function of $v$. This method is very powerful in that it solves for the effective collision frequency that explains the observed heat flux for a given scattering model, thereby reproducing the observed heat flux by definition. The downside is that one must assume a scattering model, so the approach is not very flexible. The quasi-linear operator (§ 5.2) computes a collision frequency from the observed spectrum of waves, assuming resonant wave–particle scattering. As we will show, the assumption of purely resonant interactions does not hold because of the large fluctuation amplitudes in our runs (though may be better justified at increasingly large scale separations). Finally, in § 5.3 we construct a Fokker–Planck operator by explicitly computing the short-time velocity and pitch-angle jump moments from ${\sim }500\,{\rm k}$ tracked electrons. This operator is by far the most general, but can be easily misinterpreted.

3 Numerical methods and diagnostics

3.1 Initial equilibrium states

We initialize each run with a thermally stratified ion–electron plasma and a constant background magnetic field $\boldsymbol {B}_0=B_0\hat {\boldsymbol {x}}$. We parameterize the temperature profile of each species as a function of $x$ by a factor $N>1$ using

(3.1)\begin{equation} T(x) = T_{h}\left(1-\frac{x}{NL_x}\right) , \end{equation}

where $T_{h}$ is a ‘hot’ temperature at the left wall, $x=0$. The ‘cold’ temperature at the right wall, $x=L_x$, is then given by

(3.2)\begin{equation} T_{c}=T_{h}\left( \frac{N-1}{N}\right). \end{equation}

The temperature-gradient length scale $L_T\doteq -(\nabla _{\parallel }\ln T)^{-1}$ is therefore $L_T = NL_x$.

We take the scalar pressure of species $\alpha$, $p_{\alpha }$, to be related to its density $n_{\alpha }$ and temperature $T_{\alpha }$ through an ideal equation of state, $p_{\alpha }=n_{\alpha }T_{\alpha }$. We then construct two different equilibria depending on whether a gravitational field $\boldsymbol {g}=-g\hat {\boldsymbol {x}}$ is present. For $g=0$, hydrostatic equilibrium is simply ${\rm d} p_\alpha /{\rm d}\kern0.07em x = 0$. Using (3.1) for $T_\alpha (x)$, we find that

(3.3)\begin{equation} n(x) = n_0 \left(1-\frac{x}{NL_x}\right)^{{-}1}, \quad{\rm when}\ g=0, \end{equation}

which is independent of species. This is the initial condition used in K18. While this set-up is simple, it is difficult to imagine a physically relevant scenario in which the gradients of temperature and density have equal and opposite signs, ${\rm d}\ln T_\alpha /{\rm d}\kern0.07em x = -{\rm d}\ln n_\alpha /{\rm d}\kern0.07em x$. We therefore consider a complementary set-up in which gravity balances a non-zero pressure gradient and the gradients of temperature and density have the same sign. For $g\ne 0$, hydrostatic equilibrium is then given by

(3.4)\begin{equation} \frac{{\rm d} p_{\alpha}}{{\rm d}\kern0.07em x} - q_\alpha n_\alpha E_0 + m_\alpha n_\alpha g = 0 , \end{equation}

an equation that must hold for each species individually. We take $T_i=T_e\doteq T$ in the equilibrium state, in which case $E_0$ is the equilibrium electric field that is required to maintain quasi-neutrality with a gravitational force that is greater on the ions than on the electrons. To determine this electric field, one multiplies (3.4) by $q_\alpha /m_\alpha$, sums over particle species and uses $\sum _\alpha q_\alpha n_\alpha = 0$ to find that

(3.5)\begin{equation} E_{0}= \left.\sum_{\alpha}\frac{q_{\alpha}}{m_{\alpha}}\frac{{\rm d}p_{\alpha}}{{\rm d}\kern0.07em x}\right/ \sum_{\alpha}\frac{q_{\alpha}^2n_{\alpha}}{m_{\alpha}} ={-} \frac{1}{en}\left(\frac{m_i-m_e}{m_i+m_e}\right) \frac{{\rm d} p}{{\rm d}\kern0.07em x} , \end{equation}

where in the second equality we have taken $p_i = p_e \doteq p$, $q_i = -q_e = e$ and $n_i=n_e\doteq n$. Directly summing (3.4) over species yields the standard hydrostatic equilibrium,

(3.6)\begin{equation} \sum_\alpha \frac{{\rm d} p_\alpha}{{\rm d}\kern0.07em x} = 2\frac{{\rm d} p}{{\rm d}\kern0.07em x} ={-}\varrho g, \end{equation}

where $\varrho \doteq (m_i+m_e)n$. Using (3.1) for the temperature profile, (3.6) can be straightforwardly integrated to find that

(3.7)\begin{equation} n(x) = n_0 \left( 1 - \frac{x}{NL_x}\right)^{N-1} \quad{\rm when}\ g = \frac{2T_{h}}{(m_i+m_e)L_x} . \end{equation}

In this case, both the temperature and the density decrease with ‘height’ ($x$).

Before proceeding to discuss our numerical approach for evolving these equilibria, it is worth addressing an issue with the $g\ne 0$ case, which may have caught the attention of erudite students of astrophysical fluid dynamics. The situation in which a thermally conducting magnetic field is aligned with a temperature gradient that points in the direction of gravity is known to be linearly overstable to $g$-mode perturbations when the profiles satisfy the inequalities $1 < (1-1/\gamma )({\rm d}\ln p/{\rm d}\ln T) < 2$, where $\gamma$ is the effective adiabatic index of the plasma (Balbus & Reynolds 2008).¹ Using the profiles (3.1) and (3.7), this criterion requires that $N$ satisfies $1 < (1-1/\gamma )N < 2$. Without specifying $\gamma$, which can differ from the standard monatomic $5/3$ when the plasma in question is collisionless and magnetized, we can nevertheless state that this overstability plays no role in our simulations, because its maximum growth rate is a fraction (typically ${\sim }0.1$) of the characteristic buoyancy frequency, ${\sim }(g/L_x)^{1/2} \simeq \varOmega _e (m_e/m_i)^{1/2}(\rho _e/L_x)$ (Kunz 2011, § 4.3). With our use of $m_i/m_e=1600$ and $L_x/\rho _e\ge 125$ (see § 3.2), none of our simulations are run for long enough to realize such slow growth. In fact, with the whistler instability rapidly producing an effective collisionality $\nu _{\rm eff} \sim \beta _e (v_{\text {th}e}/L_T)$, one can show that the fastest-growing $g$ mode has a parallel wavelength larger than the scale height (and, therefore, the box length).

3.2 Simulation approach and choice of free parameters

The initial conditions described in § 3.1 are advanced forward in time by the Vlasov–Maxwell set of equations using the PIC method implemented in TRISTAN-MP (Spitkovsky et al. 2019). All runs are 2.5D: the spatial simulation domain is an elongated, 2-D grid of size $L_x \times L_y$ with $L_x\gg L_y$, while particle velocities are fully three dimensional. We choose the ‘hot’ temperature at the left wall to be commensurate with an electron thermal velocity $v_{\text {th}e,0}\doteq v_{\text {th}e}(x=0)=\sqrt {2T_{h}/m_e}=c/5$, where $c$ is the speed of light. We take the same temperature ratio, $T_{h}/T_{c}=2$, across all of our runs. Thus, $N=2$ in (3.1)–(3.3) and (3.7). The heat flux driven through the box in the absence of collisions is

(3.8)\begin{align} q_{{\parallel},e 0} & \doteq m_e n_e(x=0) v_{\text{th}e}^{3}(x=0) - m_e n_e(x=L_x)v_{\text{th}e}^3(x=L_x)\nonumber\\ & = \begin{cases}(1-2^{{-}1/2}) m_en_{e0}v_{\text{th}e,0}^{3} & \text{for } \boldsymbol{\nabla} p_0=0, \\ \dfrac{1}{2}m_en_{e0}v_{\text{th}e,0}^{3} & \text{for } \boldsymbol{\nabla} p_0=\rho_0\boldsymbol{g}, \end{cases} \end{align}

which we use to normalize the computed parallel heat flux in all of our runs. We set the density at the left wall to be

(3.9)\begin{equation} n_{e0}\doteq {\rm ppc_0}=\begin{cases} 175\,{\rm ppc} & \text{for }\boldsymbol{\nabla} p_0=0, \\ 350\,{\rm ppc} & \text{for }\boldsymbol{\nabla} p_0=\rho_0\boldsymbol{g}, \end{cases} \end{equation}

where ${\rm ppc}$ is particles per cell; in both cases the electron density at the centre of the domain is ${\simeq } 262$. Electrons are first sampled in the $x$–$y$ plane according to the appropriate density distribution, whereupon their temperature is calculated following (3.1). An ion is deposited at the same spatial coordinates as each electron so that the initial condition is exactly charge neutral. Following Zenitani (2015), the particle speed is sampled from a Maxwell–Jüttner distribution with that temperature and projected onto a uniform sphere to obtain an isotropic velocity distribution. We set the electron inertial length $d_e(x=0)=7.5\Delta x$, where $\Delta x$ is the grid spacing, across all runs. The electron Debye length, $\lambda _{{D},e}(x=0)=1.5\Delta x$, is large enough with respect to the grid to minimize spurious particle heating.

The strength of the background field $B_0$ is controlled by the initial electron plasma beta parameter $\beta _{e,0}$ measured at the left boundary of the domain, $x=0$:

(3.10)\begin{equation} \beta_{e,0}=\frac{8{\rm \pi} p_{e,0}}{B_0^2}. \end{equation}

Here $p_{e,0}=p_e(x=0)=n_0T_{h}$. The size of the domain is given in units of the electron gyroradius measured at $x=0$, $\rho _{e,0}=v_{\text {th}e,0}/\varOmega _{e,0}$, with $L_y/\rho _{e,0} = 25$. We vary $L_x/\rho _{e,0}$ to change the temperature-gradient length scale $L_T$. This ensures that the expected effective whistler mean free path,

(3.11)\begin{equation} \lambda_{\text{mfp},\text{eff}}\sim \frac{L_T} {\beta_{e,0}}\sim \frac{2}{\beta_{e,0}}L_x, \end{equation}

is much smaller than $L_x$ for $\beta _{e,0}\gg 1$.

We conducted a number of simulations varying $\beta _{e,0}$, $L_T/\rho _{e,0}$ and the equilibrium state; these runs are summarized in table 1. To determine how the heat flux scales with $\beta _{e,0}$, we conduct simulations at $\beta _{e,0}=\{10,25,40,50,100\}$ using the K18 equilibrium and $\beta _{e,0}=\{10,25,40,50\}$ for the gravitational equilibrium, keeping $L_T=250\rho _{e,0}$ fixed throughout. Similarly, we fix $\beta _{e,0}=40$ while varying $L_T/\rho _{e,0}=\{250,500,1000,2000\}$ and $L_T/\rho _{e,0}=\{250,500\}$ for the K18 and gravitational equilibria, respectively, to quantify how the heat flux scales with the temperature-gradient length scale. All but one of the aforementioned runs were conducted using $m_i/m_e=1600$; there is one run with $m_i/m_e=100$ and $\beta _{e,0}=100$, designed to test the effect of the ion gyroresonance on the instability.

Table 1. For all runs performed, from left to right: run name, equilibrium initial condition, initial electron plasma beta parameter, proton-to-electron mass ratio $m_i/m_e$, length of the box along the temperature gradient $L_x$, temperature-gradient length scale $L_T$ and run time $T_{\rm run}$. For $\boldsymbol {\nabla } p_0=0$ runs used to compute an effective collision operator, we include the whistler phase speed $v_{{w}}$ measured from spectrograms and the time interval over which our Fokker–Planck coefficients are calculated, $t_{s}-t_{e}$.

3.3 Particle boundary conditions and temperature enforcement

In the $y$ direction perpendicular to $\boldsymbol {B}_0$, the boundary conditions for the particles are periodic. In order to maintain the temperature gradient across the domain and thereby drive the system to a steady state, the boundaries in the $x$ direction fix the temperatures at $x=0$ and $x=L_x$ to $T_{h}$ and $T_{c}$, respectively. This is achieved by placing particle barriers on the edges of the simulation domain and re-sampling any particle that crosses them from a half-shell distribution with the appropriate temperature. The algorithm is as follows. If, after a single particle push a particle crosses a barrier, it is pushed back to the point of collision with the wall and the current generated by the push back is added to the current density on the grid. The particle speed is then sampled from a Maxwell–Jüttner distribution in the same manner as the simulation domain was initialized in (§ 3.2); however, instead of projecting the speed onto a uniform sphere, we project onto a half-shell distribution that pushes the particle back into the active domain. The returning particle therefore acts as though it has been thermalized by a heat bath at the appropriate temperature. When this particle boundary condition is enforced, any currents deposited by a particle being rewound to the domain boundary or subsequently re-injected are summed into the total current density for that time step. That total current density is then used to calculate the fields in the next time step, ensuring that the algorithm is physically consistent and conserves charge.

3.4 Field boundary conditions and wave absorption

The electromagnetic fields, following the particles, are taken to be periodic across the $y$ boundaries. In the $x$ direction we employ the masking method for absorbing boundary conditions detailed in Umeda, Omura & Matsumoto (2001). The electric and magnetic fields are multiplied by a masking function $M(x,L_D,r)$, which is a function of $x$, the length of the damping region $L_D$ and the masking parameter $r$ $(0< r\le 1)$ given by

(3.12)\begin{equation} M(x,L_D,r) = \begin{cases} 1 - \left(r\,\dfrac{x-L_D}{L_D}\right)^2, & x\le L_D, \\ 1, & L_D < x < L_x-L_D ,\\ 1 - \left(r\,\dfrac{x-L_x+L_D}{L_D}\right)^2, & x\ge L_x-L_D. \end{cases} \end{equation}

This mask is multiplied to the right-hand side of the discretized Ampère's and Faraday's laws. This has the effect of strongly suppressing the fluctuations near the $x$-domain boundaries. We set the length of the damping region $L_D=\lambda _{{w}}/{\rm \pi} \sim 2\rho _{e,0}$, assuming a characteristic wavenumber for the whistlers satisfying $k_{\parallel }\rho _e\sim 1$, and the masking parameter $r=0.25$. The ratio of the incident wave intensity to the reflected intensity using this scheme is proportional to the amount of time the wave spends in the buffer.

4 Results: I. Whistler regulation of the heat flux

To help organize the presentation, our results are divided into two separate sections. The first (this section) focuses on how the imposed temperature gradient drives a heat flux unstable to the whistler instability, which subsequently scatters electrons and thereby regulates that heat flux to a value that is marginally unstable. This section generally serves to confirm previous work – in particular, the dependence of saturated heat flux and whistler amplitude on $\beta _{e,0}$ and $L_T$ (§ 4.2). We therefore keep this section brief, referring the reader to both K18 and RC18 for more detail. We do, however, comment on the prolonged secular phase of the instability we observe (§ 4.1), directly measure the whistler dispersion relation (§ 4.3) and note a significant reduction in heat flux by ion-cyclotron-resonant damping (§ 4.4) – all of which are novel contributions. The second section, § 5, documents our efforts to derive, using a combination of analytical arguments and simulation data, an effective collision operator describing how electrons interact with the whistler fluctuations. The latter is the principal contribution of this work.

4.1 Evolution of the instability

Each of our runs follow a qualitatively similar evolution, going through four distinct phases: (i) establishment of an unstable heat flux, (ii) exponential growth of whistler waves, (iii) sustained power-law (‘secular’) growth, and (iv) complete saturation of the HWI. Qualitatively, these four stages are manifest in figure 2, which shows 2-D plots of the magnetic-field energy in run b40 ($\beta _{e,0}=40$ and $L_T=250\rho _{e,0}$) during the latter three stages in the top, middle and bottom panels, respectively. In the exponential phase (ii), figure 2(a) shows oblique, $k\rho _{e,0}\sim 1$, whistlers throughout the domain. The grey bars at the edges of the domain correspond to the masked regions described in § 3.4. In the secular phase (middle panel), we see a shift towards smaller wavenumbers and stronger field fluctuations, a trend that continues until saturation (third panel).

Figure 2. Perturbed magnetic energy in run b40 at the beginning of the exponential phase (a), end of exponential phase (b) and in the saturated state (c).

More quantitatively, in figure 3 we plot the box-averaged magnetic-field fluctuations and heat flux as functions of time, with panel (a) corresponding to a scan in $\beta _{e,0}$ at fixed $L_T=250\rho _{e,0}$ and panel (b) corresponding to a scan in $L_T$ at fixed $\beta _{e,0}=40$. For the first ${\sim }100\unicode{x2013}500\varOmega _{e,0}^{-1}$, depending on $L_T/\rho _{e,0}$, the whistler instability is inactive. This is a consequence of our initial conditions: we initialize the domain with a temperature gradient, but with no heat flux. It therefore takes some fraction of the sound crossing time $T_c=L_x/v_{\text {th}e,0}$ for the heat flux – and the associated anisotropy in the distribution function – to grow large enough for the whistlers to overcome cyclotron damping.

Figure 3. Time evolution of $\delta B^2$ and $q_{\parallel }$ for both the $\boldsymbol {\nabla } p_0=0$ (a i,iii,b i,iii) and $\boldsymbol {\nabla } p_0=\rho _0\boldsymbol {g}$ (a ii,iv,b ii,iv) equilibria (a) as a function of the electron plasma beta $\beta _{e,0}$ for $L_T=250\rho _{e,0}$ and (b) as a function of the temperature-gradient length scale $L_T$ for $\beta _e=40$.

Once the heat flux is sufficiently strong, the instability enters an exponential phase, which is predicted by the whistler growth rate (2.2). At $k_{\parallel }\rho _e= 1$,

(4.1)\begin{equation} \frac{\gamma_{{w}}}{\varOmega_e} = \frac{C}{\beta_e}\left(\frac{v_{\text{th}e}\beta_e}{2\nu L_T}-1\right) , \end{equation}

where $C=\sqrt {{\rm \pi} }{\rm e}^{-1}$. The time rate of change of the whistler fluctuation energy at this scale is therefore

(4.2)\begin{equation} \frac{\partial }{\partial t}\delta B^2=2\gamma_{{w}} \delta B^2=\frac{2C}{\beta_e} \left(\frac{\beta_e v_{\text{th}e}}{2L_T \nu}-1\right)\varOmega_e\,\delta B^2. \end{equation}

In the analysis of Bott et al. (2024), $\nu$ is the scattering rate associated with a weakly collisional background and that results in the whistler-unstable steady-state distribution function (2.1). We argue that we, too, have a (very) weakly collisional background in the form of finite-particle-number noise, which scatters electrons at a rate $\nu _{\rm PIC}$. We additionally have effective scattering from whistler fluctuations at a rate $\nu _{\rm eff}$ (see (2.8)), so we take $\nu =\nu _{\rm PIC} + \nu _{\rm eff}$. Exponential growth (ii) occurs in the limit $\nu _{\rm eff} \ll \nu _{\rm PIC}\ll \beta _ev_{\text {th}e}/L_T$, in which the $-1$ term in (4.2) can be ignored and the scattering rate is approximately independent of the magnetic field. This rapid growth is sustained only for a short time – approximately $10-100\varOmega _e^{-1}$, depending on $\beta _{e,0}$ and $L_T$ – until the effective scattering rate by the magnetic-field fluctuations, $\nu _{\rm eff}$, becomes comparable to $\nu _{\rm PIC}$ and the whistler fluctuations begin to feed back on themselves. In the limit where $\nu \simeq \nu _{\rm eff}\sim \varOmega _e(\delta B/B_0)^2$, the instability grows linearly in time following

(4.3)\begin{equation} \frac{\partial }{\partial t}\delta B^2\sim \frac{Cv_{\text{th}e}}{L_T} B_0^2 \quad\Longrightarrow\quad \frac{\delta B^2(t)}{B^2_0} - \frac{\delta B^2(t_{\rm sec})}{B^2_0} \sim \frac{Cv_{\text{th}e}}{L_T} (t-t_{\rm sec}), \end{equation}

where $t_{\rm sec}$ marks the start of the secular phase. For all runs at $L_T/\rho _{e,0}=250$, $t_{\rm sec}\sim 100\varOmega _{e,0}^{-1}$; this value increases up to ${\sim}500\varOmega _{e,0}^{-1}$ at $L_T/\rho _{e,0}=2000$. The secular phase lasts ${\sim }10^{3-4}\varOmega _{e,0}^{-1}$, depending on the simulation, until the magnetic energy saturates at a value that is dependent upon $\beta _e$ and $L_T$. This phase (iii) was not reported in earlier numerical work (K18; RC18).

4.2 Saturated heat-flux scaling with $\beta _e$ and $L_T$

Once the fields are saturated, they are able to maintain a temperature gradient throughout the box, which we show in figure 4. For the $\boldsymbol {\nabla } p_0=0$ set-up, the saturated temperature gradient in the centre of the domain is close to the initialized temperature gradient. However, for the $\boldsymbol {\nabla } p_0=\rho _0\boldsymbol {g}$ set-up, the gradient is ${\sim }60\,\%$ less, likely due to the stronger equilibrium heat flux. In figure 5 we plot the saturated box-averaged heat flux (c,d) and magnetic-field perturbation energy (a,b) versus $\beta _{e,0}$ (a,c) and $L_T/\rho _{e,0}$ (b,d). These results are consistent with the scalings

(4.4a)\begin{gather} \frac{\langle q_{{\parallel},e}\rangle}{q_{{\parallel},e0}} \sim \frac{1}{\beta_e}, \end{gather}

(4.4b)\begin{gather}\frac{\langle\delta B^2\rangle}{B_0^2} \sim\frac{\beta_e}{L_T/\rho_e}, \end{gather}

and are independent of the equilibrium set-up. Our result (4.4) not only confirms the work of both K18 and RC18, but extends the validity of the scaling to larger scale separations and to an equilibrium set by gravity. The empirical scalings (4.4a) and (4.4b), along with the analysis in § 2.2, strongly suggest that the relevant dimensionless parameter in the saturation of the HWI is $\beta _{e,0}\rho _{e,0}/L_T$. For the remainder of the paper, therefore, we normalize the magnetic energy $\langle \delta B^2\rangle /B_0^2$ and all effective collision frequencies using $(\beta _{e,0}\rho _{e,0}/L_T)\varOmega _{e0}$.

Figure 4. Average temperature as a function of distance along the temperature gradient. The initial temperature profile is shown by the black dotted line. In saturation, the $\boldsymbol {\nabla } p_0=0$ runs (save b10) support a temperature gradient close to the initial gradient near the centre of the simulation domain. The runs with a gravitationally supported pressure gradient, however, only support a temperature gradient that is ${\approx }60\,\%$ of its initial value.

Figure 5. Normalized magnetic fluctuation amplitude (a,b) and heat flux (c,d) versus $\beta _{e,0}$ (a,c) and $L_T\rho _{e,0}$ (b,d) for the K18 set-up (blue) and gravity set-up (orange). The heat flux and fluctuation amplitude show good agreement with (4.4a) and (4.4b), respectively.

4.3 Magnetic-field spectra and spectrograms

We compute the time-averaged magnetic energy spectrum for each run over the final $10^3\varOmega _{e,0}$. Results are plotted in figure 6. The spectra are qualitatively similar between runs and are characterized by a sharp break around $k\rho _{e,0}\sim 0.6\unicode{x2013}0.7$. At smaller $k$, the spectral index is approximately zero, while the index (or indices) at larger $k$ are more difficult to interpret due to the noise in the spectrum. For $k\rho _{e,0}>1$, the spectra can be characterized by an index of $-4$, which was also reported in Zhdankin et al. (2017) for PIC simulations of driven turbulence in magnetized, collisionless, relativistic pair plasma. The two spectra that differ quantitatively are runs b10 and b40x8, which have a pronounced peak at $k\rho _{e,0}\sim 0.6$.

Figure 6. One-dimensional spectra of $\delta B^2$ for all runs versus $k_x$ (a) and $k_y$ (b), computed by integrating along the free dimension. For either direction, the spectrum is consistent with a power-law index of $-4$ for $k\rho _e\ge 1$.

We determine the whistler dispersion relation in our simulations by calculating the spectrogram of the combination $B_y+{\rm i} B_z$, which reveals information about the polarization of the waves: positive wave frequencies $(\omega > 0)$ correspond to left-hand-polarized waves, while $\omega < 0$ corresponds to right-hand-polarized waves. We show results for runs b100 and b40x4 in figure 7; spectrograms from the other runs were qualitatively similar. Data was taken over the final ${\sim }4000\varOmega ^{-1}_e$ of each run to obtain sufficient resolution at low frequencies. We also applied Gaussian windows in $x$ and $t$ with standard deviations $1/4$ and $1/6$ of the total window widths, respectively, to ensure non-periodic edge effects were minimized. The black dashed lines in figure 7 are best fits for the wave dispersion relation multiplied by a free parameter $A$, i.e.

(4.5)\begin{equation} \frac{\omega}{\varOmega_e}\simeq A\frac{(k_{{\parallel}} \rho_e)(k\rho_e)}{\beta_e}, \end{equation}

which implies a phase speed

(4.6)\begin{equation} v_{{w}} \simeq A\frac{v_{\text{th}e}}{\beta_e}k_{{\parallel}}\rho_e. \end{equation}

The best fit value of $A$, expressed as a multiplicative constant of $\beta _{e,0}v_{{w}}/v_{\text {th}e,0}$ with ${k_{\parallel }\rho _{e,0}=1}$, is reported for each $\boldsymbol {\nabla } p_0=0$ run in table 1. In general, $A$ increases with $\beta _{e,0}$: from ${\simeq }0.20$ at $\beta _{e,0}=10$ to ${\simeq }0.28$ at $\beta _{e,0}=100$. Our analysis indicates that the whistler phase speed is slower than that of the cold plasma dispersion prediction at $k_{\parallel }\rho _e=1$ (cf. (2.3)). For the remainder of the work, we treat $v_{{w}}$ as a constant, given by the $k_{\parallel }\rho _{e,0}=1$ values reported in table 1. Given the complexity of the analysis in the following sections, we view this as an expedient assumption and leave the effects of a more general treatment to future work.

Figure 7. Spectrograms of $B_y+{\rm i}B_z$ for runs b100 (a,c) and b40x4 (b,d) for $k_y\rho _e=0$ (a,b) and $k_y\rho _e=1$ (c,d). Energy is concentrated in right-hand polarization ($\omega <0$) for parallel modes, but energy does go into left-hand polarization for oblique modes, as expected with whistler waves. Black dots represent the frequency $\omega$ with the highest Fourier amplitude at each $k_{\parallel }$; black dashed lines are the best fits to these points. The best fit lines differ from the cold plasma dispersion by a factor of $0.32\unicode{x2013}0.17$, with an average value of $0.23$.

4.4 Effect of ion-Landau and ion-cyclotron resonances

Runs b40 and b40g, with $\beta _{e,0}=40$ and $m_i/m_e=1600$, allow us to assess any impact from the ion-Landau resonance. Run b100m, with $\beta _{e,0}=100$ and $m_i/m_e=100$, tests the effect of the ion-cyclotron resonance. As explained at the end of § 2.2, ions may be Landau resonant with whistler waves when the phase velocity of the former is of the order of the thermal speed of the latter, i.e. when $\beta _e\sim (m_i/m_e)^{1/2}$. The ions can also be cyclotron resonant with whistlers if the wave frequency equals the ion-cyclotron frequency – this occurs when $\beta _e\sim m_i/m_e$. Re-examining figure 5(a,c) shows that the runs at $\beta _{e,0}=40$ do not stray appreciably from the $q\sim \beta _e^{-1}$ scaling and so we conclude that effects due to the ion-Landau resonance are negligible. Komarov et al. (2018) report a $20\,\%$ increase in saturated heat flux and $15\,\%$ decrease in magnetic-field fluctuation energy. Their conclusion comes from comparing a $\beta _e=15$, $m_i/m_e=225$ run to one at the same $\beta _e$, but at $m_i/m_e=\infty$. In our tests we saw changes in saturation levels with immobile ions at $\beta =10$, so it is possible that their conclusion is affected by comparison to a run with immobile ions instead of a different mass ratio, resulting in a deviation higher than ours.

In contrast to K18, we find that the ion-cyclotron resonance affects the saturated state of the HWI significantly. In figure 8 we plot the box-averaged magnetic-field fluctuation amplitude and heat flux as functions of time for two $\beta _{e,0}=100$ runs, one with ${m_i/m_e=1600}$ and the other with $m_i/m_e=100$. While the heat-flux accumulation and exponential stages are nearly identical between runs, the cyclotron-resonant run saturates much sooner than does the run with mass ratio $1600$, with a $25\,\%$ smaller fluctuation amplitude and $57\,\%$ higher heat flux. Our result that the ion-cyclotron resonance significantly alters the saturated heat flux relative to the Landau resonance is consistent with our simple explanation of the instability in § 2. Evidently, cyclotron-resonant ions can strongly damp the oblique whistlers that are necessary for scattering heat-flux-carrying electrons. We also conjecture that Landau-resonant ions preferentially damp the parallel whistlers and have a diminishing effect on the heat flux. An in-depth analysis of the interaction between whistlers and ions, however, is outside the scope of this work.

Figure 8. Box-averaged magnetic perturbation (a) and box-averaged heat flux (b) versus time for two $\beta _{e,0}=100$, $L_T=250\rho _{e,0}$ runs with $m_i/m_e=1600$ and $m_i/m_e=100$. When $\beta _e\sim m_i/m_e$, ions are in cyclotron resonance with the whistler waves. This results in both a reduction of $\delta B^2/B_0^2$ and an increase in the saturated heat flux by a factor of 2.

5 Results: II. An effective collision operator for heat-flux-driven whistler turbulence

In this section we present detailed methods and results for three different means of obtaining a collision operator from our simulations. Section 5.1 is devoted to a method based upon the Chapman–Enskog expansion. After stipulating a model collision operator that captures pitch-angle scattering in the frame of the whistlers, we derive the implied ‘Chapman–Enskog’ scattering frequency $\nu _{\rm CE}$ in § 5.1.1, characterize the equilibrium electron distribution function in § 5.1.2 and present our numerical findings in § 5.1.3. Section 5.2 focuses on the quasi-linear operator; we define the operator in § 5.2.1 and give our numerical results in § 5.2.3. Finally, we detail the Fokker–Planck method in § 5.3. We define the resulting operator in § 5.3.1 and discuss the relevant time scales, particle statistics and limitations in § 5.3.2. In §§ 5.3.3 and 5.3.4 we present our numerical results for the velocity-space Fokker–Planck operator in $(v,\xi$) coordinates.

For all of these methods, we leverage a key finding from § 4: that the energy of magnetic-field fluctuations in the saturated state of the HWI and the effective collisionality implied by the steady-state heat flux both scale proportionally with $\beta _e (\rho _e/L_T)$. We organize the various features of our model collision operators according to this dimensionless free parameter.

5.1 Chapman–Enskog pitch-angle-scattering operator

Our first method to obtain a collision operator for whistler turbulence relies on a Chapman–Enskog expansion of the electron kinetic equation. In the usual expansion, one solves for the collisional transport resulting from the free-energy gradients (usually velocity or temperature) in a fluid with a known scattering operator. We instead solve for the collision operator that self-consistently describes the transport and fluid gradients measured in our simulations. Using this method, one assumes the mathematical form of the collision operator and solves for the implied diffusion coefficient as a function of phase-space variables. This is simultaneously a strength and a limitation, as the diffusion coefficient is designed to explain all the observed transport given the model.

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 1966) and is the same as chosen by Drake et al. (2021). 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 2018), 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 Quasi-linear operator

Our next model collision operator for the saturated HWI is a quasi-linear operator that depends explicitly on the energy spectrum of the electromagnetic fluctuations. This dependence arises from averaging the electron kinetic equation over a few wave periods in time and wavelengths in space and separating the electron distribution function into its mean and fluctuating parts; one then finds that the time rate of change of the mean distribution function is due to the mean of the interaction of the perturbed fields with the perturbed distribution function. If one further assumes that the perturbed distribution function is the linear response to the fields (hence, the qualifier ‘quasi-linear’), one obtains an operator that depends only on the fields resonantly interacting with particles. A benefit of this method is a rich insight into the physics of the saturation mechanism. Its main disadvantage is that the quasi-linear operator contains a number of strict assumptions that are often suspect. For our results in particular, wave–particle resonances are likely to be broadened.

5.2.1 Definition

We adopt the electromagnetic quasi-linear operator given in Stix (1992) (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 1992; Barnes 1966). 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 Fokker–Planck method

The most direct and detailed method for calculating the HWI collision operator is the Fokker–Planck method. This method is particularly powerful because it calculates drag and diffusion directly from the statistics of the particle motion. Moreover, this method requires no assumptions about the physics underlying the drag and diffusive processes. That being said, it does require statistical assumptions; in particular, the Fokker–Planck method works well when there is a large separation between the autocorrelation time of the forcing and the collisional time scale. For a number of our runs, this scale separation is not especially large and results in some complications that are otherwise avoided by the other methods we have employed.

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 (1957) 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 2002). 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 1930) 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 1991), 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 Results: III. A model effective collision operator for saturated HWI

In § 5 we used our numerical simulation data to calculate three model collision operators describing the drag and diffusion of electrons as they interact with saturated HWI fluctuations. Here we attempt to unify these operators, highlighting their essential ingredients, proposing a model that captures them and testing whether that model can result in a heat flux that matches, within reason, the heat flux measured directly in the simulations. While doing so, we feature similarities and differences between existing model operators and emphasize what aspects of our model operator could benefit from further refinement.

6.1 Salient features of a model effective collision operator

First and foremost, all of our model operators show scattering rates that scale as $\beta _ev_{\textrm {th}e}/L_T$, with the exact rate depending on the model. We also demonstrated that pitch-angle scattering dominates over drag and diffusion in velocity (§ 5.3), and showed in § 5.2.1 that this is expected analytically from an electromagnetic quasi-linear operator for waves with highly subthermal phase velocity. Such an operator includes a drag term proportional to the phase velocity of the scattering waves, which we showed is equivalent to boosting the pitch-angle-scattering operator to the wave frame (§ A.2). Finally, we found that the consequent advection of scattered particles at the phase speed $v_{{w}}\sim v_{\textrm {th}e}/\beta _e$ is of the same order as the diffusive scattering; unsurprisingly then, we show in § 6.3 that each effect is responsible for approximately half of the observed heat flux.

The pitch-angle average of the pitch-angle-scattering frequency derived from each of our methods all show two power laws in $v$: one decreasing as a function of $v$ for subthermal electrons (physically, due to transit-time damping) and the other increasing like ${\sim }(v/v_{\textrm {th}e})^3$ for superthermal electrons (physically, due to cyclotron-resonant scattering). The pitch-angle dependence of the operators, however, is more complicated. The quasi-linear operator showed (a) a structure that was a function of $v_{\parallel }$ rather than of $v$ and $\xi$ separately, and (b) that the scattering frequency for electrons with $v_{\parallel }<0$ was an order of magnitude larger than for electrons with $v_{\parallel } >0$. For runs with $L_T/\rho _{e,0}=250$, the Fokker–Planck operator showed a largely symmetric structure, but with some evidence that the scattering was a function of $v_{\parallel }$. At larger $L_T/\rho _{e,0}$, however, the operator began to show more qualitative similarities to the quasi-linear operator, particularly with a resonant structure and a strong asymmetry between electrons with positive versus negative parallel velocities.

6.2 Basic model collision operator

The preponderance of evidence presented in this work points towards whistlers acting as pitch-angle scatterers in a frame moving at the whistler phase velocity

(6.1a)\begin{equation} C[f]=\frac{\partial }{\partial \xi} \left[\frac{1-\xi^2}{2}\nu(v,\xi) \left(\frac{\partial f}{\partial \xi} + v_{{w}}\frac{\partial f}{\partial v}\right)\right], \end{equation}

where we have measured the wave phase velocity in simulation to be

(6.1b)\begin{equation} v_{{w}} \simeq 0.23 v_{\text{th}e}/\beta_e . \end{equation}

Following § 6.1, we take the model scattering frequency to be given by the quasi-linear value, repeated here, with a few modifications:

(6.1c)\begin{equation} \frac{\nu_{\rm model}}{\varOmega_e}=2{\rm \pi} \sum_n \int_0^\infty\,{\rm d} (k_{{\parallel}}\rho_e)\, \delta\left(\frac{\omega(k_{{\parallel}}\rho_e)}{\varOmega_e}-k_{{\parallel}}\rho_e\frac{v_{{\parallel}}}{v_{\text{th}e}}+n\right)\frac{1}{L_x}\left|\frac{\varPsi_{n,k_{{\parallel}}}}{B_0}\right|^2. \end{equation}

In order to make the full quasi-linear expression more tractable for our model, we remove the explicit dependence of the quasi-linear scattering frequency on $k_\perp$. We accomplish this simplification by assuming that the oblique resonance condition can be approximated by the parallel resonance condition (5.18) and that the effect of integrating over oblique modes – specifically the reduction in resonant power for electrons interacting with the left-handed wave component – can be encapsulated by a bespoke numerical weighting of the parallel spectrum for each of the principle resonances $n=[-1,0,1]$, viz.

(6.1d)\begin{align} \frac{1}{L_x}\left|\frac{\varPsi_{n,k_{{\parallel}}}}{B_0}\right|^2& =\frac{\langle \delta B^2\rangle}{B_0^2}\left\{ \left[\left(\frac{k_{{\parallel}}\rho_e - 0.5}{0.4}\right)^3 + 1\right]^{{-}1}\delta_{n,-1}+ 0.1\left[\left(\frac{k_{{\parallel}}\rho_e - 0.5}{0.25}\right)^{2.5} + 1\right]^{{-}1}\delta_{n,0}\right.\nonumber\\ & \quad \left. + \,0.1\left[\left(\frac{k_{{\parallel}}\rho_e - 0.5}{0.5}\right)^3 + 1\right]^{{-}1}\delta_{n,1} \right\}. \end{align}

The integral in (6.1c) has bounds $[0,\infty ]$, reflecting the fact that whistler waves travel only down the temperature gradient ($k_\parallel >0$).

Before we discuss our model spectrum in more detail, the use of the parallel whistler dispersion relation in our model resonance condition requires some justification. The oblique cold plasma dispersion relation (4.5) differs from the parallel dispersion relation by a factor $1/\cos \theta$, where $\theta$ is again the angle between the wave vector and the magnetic field. For $|v_{\parallel }/v_{\textrm {th}e}|\sim 1$ and $\theta =0$, the wave frequency is a factor $1/\beta _e$ smaller than the (order unity) other terms and $k_{\parallel }\rho _e\sim 1$. Only when $\cos \theta \sim 1/\beta _e$ does the obliquity of the wave begin to affect the solution – and even then only by an order-unity factor. At the same time, $\beta _e\gg 1$ implies that waves are propagating nearly perpendicular to the mean field. Waves in this parameter regime are almost certainly kinetic Alfvén waves, not whistler waves, and so the resonance condition in (6.1c) is no longer applicable. More realistic (smaller) propagation angles will change the resonant parallel wavenumber by a factor smaller than unity, thus justifying our use of the parallel dispersion relation.

For each of our runs, we fit (6.1d) to the numerical magnetic-field spectra in saturation; a plot of the fit for run b40 can be found in figure 24 along with $L_x^{-1}|\varPsi _{n ,k_{\parallel }}/B_0|^2$ calculated from the numerical spectrum using (5.12b). Overall the fit is good, especially so for the $n=-1$ and $n=0$ spectra. However, an issue with our fit becomes immediately clear: $|\varPsi _{n,k_{\parallel }}/B_0|^2$ for $n=-1$ and $n=1$ should approach each other at small $k_{\parallel }\rho _e$. As it pertains to the heat flux implied by the model, however, this deviation makes little difference as it is the fluctuations at the high-$k_{\parallel }\rho _e$ end of the cascade that matter most. One might argue that, in terms of energy, including this high-$k_{\parallel }$ part of the cascade makes very little difference to the operator, since the spectra die away very quickly rightward of their peaks. However, the implied heat flux of a collision operator of the form (6.1a) depends on the inverse of $\nu (v,\xi )$ (see § 6.3, specifically (6.2) and (6.3)). The heat flux is therefore largely insensitive to regions in velocity space in the vicinity of $v_{\parallel }=v_{\textrm {th}e}$ where particles are scattering and is instead dominated by regions where they are not. The latter include two populations of particles: those with small enough $v_{\parallel }$ that they bounce, and therefore are trapped, in the whistler wave frame; and those that are still passing, but scatter inefficiently because they are resonant with low-energy fluctuations. Proper modelling of both the spectrum at low energies and the effect of resonance broadening and particle trapping is crucial to obtain a scattering frequency that accurately reproduces the observed heat flux.

Figure 24. Plot of $L_x^{-1}|\varPsi _{n ,k_{\parallel }}/B_0|^2$ from run b40 using (5.12b) in solid lines and model (6.1d) in dashed lines, for $n=[-1,0,1]$. The spectra calculated from simulation vary from run to run; the model presented here is a good qualitative fit across all runs.

In figure 25 we plot the model collision frequency (6.1c) normalized by ${\beta _{e,0}\rho _{e,0}\varOmega _{e,0}/L_T}$ for all ${\beta _{e,0}=40}$, $\boldsymbol {\nabla } p_0=0$ runs using the model spectrum (6.1d). Qualitatively, results match those of the quasi-linear operator in figure 11: scattering for electrons travelling up the temperature gradient is higher than that for electrons travelling down the gradient by a factor of ${\sim }10$ and there is a small scattering rate at $v_{\parallel }/v_{\textrm {th}e}\ll 1$ due to transit-time damping. The $n=0$ features in the model, however, are significantly sharper than those in the quasi-linear case, likely a result of using a continuous model for the spectrum rather than a discrete calculation from simulation.

Figure 25. Two-dimensional plots of the model pitch-angle collision frequency $\nu _\textrm {model}$ (6.1) for all $\beta _{e,0}=40$ runs. The scattering frequency is normalized to $\beta _{e,0}\rho _{e,0}\varOmega _{e,0}/L_T$ and dashed lines correspond to contours of constant $v_{\parallel }$. These results are qualitatively similar to the quasi-linear scattering frequency, shown in figure 11.

6.3 Implied heat flux for model collision operators

We calculate numerically the heat flux implied by the operators described in §§ 5.2.3–5.3.4 by first computing the implied perturbation to the distribution function via

(6.2)\begin{equation} \langle\, f_{e1}\rangle_{\phi} ={-}\int \,{\rm d}\xi\, \left[\frac{1}{\nu(w,\xi)}\left(w\nabla_\parallel f_{e0}+\frac{\nabla_{{\parallel}} p_e}{m_en_e}\frac{{\rm d} f_{e0}}{{\rm d} w}\right) +v_w\frac{{\rm d} f_{e0}}{{\rm d} w}\right], \end{equation}

using $f_{e0}=f_{{\rm M}e}$ and $\nu (v,\xi )$ from the appropriate model operator. The first term on the right-hand side of (6.2) is the perturbation of the distribution function implied by the diffusive flux, which has a scattering frequency $\nu (w,\xi )$, and the last term represents the perturbation implied by advection of particles at the whistler phase velocity $v_{{w}}$. The perturbed distribution is then used to compute the parallel electron heat flux according to

(6.3)\begin{align} q_{{\parallel} e}& =\int \,{\rm d}^3\boldsymbol{w}\,\frac{1}{2}m_ew^2w_\parallel \langle\, f_{e1}\rangle_\phi\nonumber\\ & ={-}\frac{m_e}{2}\int_0^\infty \,{\rm d} w\,w^5\left[\int\,{\rm d}\xi\,\xi\left(w\nabla_\parallel f_{e0}+\frac{\nabla_{{\parallel}} p_e}{m_en_e}\frac{{\rm d} f_{e0}}{{\rm d} w}\right)\int_0^\xi\frac{{\rm d}\xi'}{\nu(w,\xi')} +\frac{2v_{{w}}}{3}\frac{{\rm d} f_{e0}}{{\rm d} w} \right]. \end{align}

Before discussing the results of this calculation, we should note that (6.2) was derived assuming a steady state, i.e. $\partial f_e/\partial t=0$. This assumption may be violated in regions of phase space where our model scattering frequencies are small. Electrons in these regions must experience many collisions to approach a steady state, so $\partial f_e/\partial t\rightarrow 0$ on time scales much longer than $1/\nu (v,\xi )$. Including this effect would serve to reduce the strong variation in $\langle\, f_{e1}\rangle _\phi$; however, it likely will not alter the qualitative results that follow.

In figure 26 we plot normalized heat fluxes as a function of $\beta _{e,0}\rho _{e,0}/L_T$: in blue we show the box-averaged heat flux measured in the $\boldsymbol {\nabla } p_0=0$ runs; and in orange, red, purple and cyan we show the normalized heat flux implied, respectively, by advection at the whistler phase speed, by the quasi-linear operator, by the Fokker–Planck operator and by the model from Drake et al. (2021) (discussed in § 6.6). In black, green and grey we show, respectively, heat fluxes implied by our model (6.1), the same model but where we neglect the heat-flux contribution of electrons expected to be trapped (see § 6.4.1 for evidence of trapped electrons in one of our runs, as well as § 6.4.2 for heat flux calculations) and our semi-empirical model (6.10). Points marked by an ‘x’ are runs with $\beta _{e,0}=40$ and points marked by circles are runs with $L_T=125\rho _{e,0}$; the former are connected with a dotted line and the latter are connected with a dashed line. Run b40, which belongs to both groups, is denoted with an ‘x.’ The blue dashed line is the average of the measured heat flux across all runs. We find that the heat flux implied by advection is on average half of the total measured heat flux, suggesting that the advection and pitch-angle diffusion contribute to the heat flux equally.

Figure 26. Heat flux measured in simulations ($q_{\parallel \textrm {measured}}$) compared with the heat flux implied by advection at the whistler phase speed ($q_{v_{{w}}}$), quasi-linear operator ($q_{\parallel \textrm {QL}}$), Fokker–Planck operator ($q_{\parallel \textrm {FP}}$), Drake et al. (2021) ($q_{\parallel \textrm {Drake2021}}$), our model (6.1) ($q_{\parallel \textrm {model}}$), our model only including passing electrons (6.5) assuming the trapped-passing boundary (6.4) ($q_{\parallel \textrm {model},\xi _\textrm {crit}}$) and our semi-empirical model (6.10) ($q_{\parallel {\textrm {model},\textrm {SE}}}$). Points denoted by an ‘x’ are runs with $\beta _{e,0}=40$ and are connected by a dotted line. Runs with $L_T=125\rho _{e,0}$ are denoted with a circle (except for run b40) and are connected with dashed lines; the blue dashed line is the average measured heat flux across the runs.

Overall, we find that the heat flux implied by our Fokker–Planck operator agrees well with the heat flux measured in our simulations. While there is some deviation of the former from the latter, particularly for $\beta _{e,0}\rho _{e,0}/L_T=0.04-0.1$, the deviations appear to be random, i.e. not a function of scale separation. We therefore associate these fluctuations with some error in our Fokker–Planck measurement process, perhaps due to our choice of jump moment or limited velocity-space resolution. In contrast, the heat fluxes implied by both our quasi-linear operator and basic model generally increase with $\beta _{e,0}$, as evidenced by the increasing trend of $(q_{\parallel \textrm {QL}})$ and $(q_{\parallel \textrm {model}})$ for runs with $L_T/\rho _{e,0}=250$ as a function of $\beta _{e,0}\rho _{e,0}/L_T$ (all runs with $L_T/\rho _{e,0}=250$ are denoted by dots in figure 26 except for run b40, which is denoted by an ‘x’ at $\beta _{e,0}\rho _{e,0}/L_T=0.16$). We attribute this trend to transit-time damping becoming weaker at larger $\beta _{e,0}$. For $n=0$, the resonant parallel wavenumber is proportional to $\beta _ev_\parallel /v_{\textrm {th}e}$ (5.18). At a given parallel velocity, therefore, the resonant $k_\parallel \rho _e$ increases proportionally to $\beta _e$. The result, given the steep spectrum we observe, is that transit-time damping is only active over an ever-narrower region of phase space close to $v_\parallel =0$ as $\beta _e$ increases, the result being an increasingly larger heat flux.

While $q_{\parallel \textrm {QL}}$ and $q_{\parallel \textrm {model}}$ exhibit the same qualitative trend due to their similar construction, the two imply very different levels of heat flux. In both models, electrons scatter resonantly off magnetic-field fluctuations whose amplitudes scale proportionally to some negative power of $k_{\parallel }\rho _e$, resulting in scattering frequencies that decrease rapidly as $|v_{\parallel }/v_{\textrm {th}e}|$ decreases from $|v_{\parallel }/v_{\textrm {th}e}|\sim 1$. Physically, low-$v_{\parallel }$ electrons interact via exact resonance only with high-$k_\parallel$ waves, which do not have appreciable power; this implies a large perturbation to the distribution function in that region of velocity space, which in turn contributes to a large heat flux. The critical difference between our quasi-linear operator and our model is that our basic model spectrum does not include the noise floor (due to a finite number of PIC particles) present in the actual simulated spectra used to calculate the quasi-linear operator. For a direct comparison of these spectra for run b40, see figure 24; the spectrum measured from our simulations levels out for $k_{\parallel }\rho _e\gtrsim 10$ while the model spectrum continues toward $0$. We can therefore expect noise to contribute to the cyclotron resonances of our quasi-linear model, at least in run b40, when $|v_{\parallel }/v_{\textrm {th}e}|\lesssim 0.1$. Referring to parallel velocity averages of the quasi-linear operator shown in figure 12, $|v_{\parallel }/v_{\textrm {th}e}|\sim 0.1$ appears to be a plausible estimate for noise to affect the scattering rate; however, transit-time damping complicates this picture somewhat. The only outlier in figure 12 is run b100, which has the weakest transit-time damping and PIC noise relative to fluctuation amplitude, resulting in the high implied heat flux in figure 26. Even though our quasi-linear operator does imply heat fluxes close to our simulation measurements for $\beta _{e,0}\rho _{e,0}/L_T\lesssim 0.1$, the observed scaling of the implied heat flux with $\beta _e$ attributed to transit-time damping and our difficulties in disambiguating it from the effects of cyclotron scattering and PIC noise call into question the applicability of the basic model to explain the saturated HWI.

There is evidently a physical reality to the Fokker–Planck results that is not captured by the quasi-linear numerical result and our model. The incompatible scaling of transit-time damping rate with $\beta _e$ found from our numerical quasi-linear and simple models also strongly suggest this to be the case. In other words, there must be something scattering electrons with $v_{\parallel }/v_{\textrm {th}e}\ll 1$ that is not captured by resonant physics. This is a manifestation of perhaps the longest-standing issue in quasi-linear theory called the $90^{\circ }$ scattering problem; we discuss this issue in detail in § 6.4. One physical effect that can alleviate the $90^{\circ }$ problem is trapping in large-amplitude waves; see § 6.4.1 for evidence of trapped electrons in our runs. If we treat these trapped electrons as if they have no contribution to the perturbed distribution function, and therefore, to the model diffusive heat flux (see § 6.4.2 for details), then the model produces implied heat fluxes, labelled $q_{\parallel \textrm {model,}{\xi _\textrm {crit}}}$ and shown in green in figure 26, which are less than the whistler advection value, i.e. the diffusive heat-flux contributions are negative. We argue in § 6.4.2 that this is a specific consequence of our model operator, as the same methodology applied to other operators produce reasonable results. It is likely that a more self-consistent treatment of the trapped population is required. In addition to this model where trapped electrons have no effect on $f_{e1}$, we also consider a model where electrons have a minimal scattering rate that is agnostic of a physical explanation and purely motivated by our Fokker–Plank scattering frequency. We call this a semi-empirical model, and cover it in detail in § 6.4.5. The heat flux implied by the model is plotted in figure 26 with grey markers and labelled by $q_{\parallel \textrm {model,SE}}$. This model matches simulation results well; however, the lack of a physics-based explanation for the scattering floor should give pause to a reader who wishes to extrapolate this model to astrophysical scale separation. Finally, another way around the $90^{\circ }$ scattering problem is to ignore the pitch-angle dependence of the scattering frequency entirely, as in the model from Drake et al. (2021), which we discuss in § 6.6. As shown in figure 26 by the cyan markers, this model does reproduce the proper heat-flux scaling and is reasonably consistent with our measurements from simulation. However, there are concerns with this operator in addition to the omission of any pitch-angle dependence of the scattering frequency; we refer the reader to § 6.6 for details on these concerns.

6.4 Scattering across the gap

As discussed at the end of § 6.3, the heat flux implied by our model scattering frequency (6.1) is unphysically large because the scattering rate approaches $0$ geometrically as ${|v_{\parallel }/v_{\textrm {th}e}|\rightarrow 0}$. This is a variation on the well-known ‘$90^{\circ }$ problem’, which also afflicts theories of cosmic-ray transport in the Galactic and intracluster contexts. Just as in our formulation, the scattering rate for cosmic rays approaches $0$ for small parallel velocities because the limited power in resonant (in this case, Alfvén) waves. It is argued that some nonlinear mechanism is in fact required to scatter cosmic rays across $v_{\parallel }=0$ and isotropize the distribution (e.g. Shalchi 2009; Holcomb & Spitkovsky 2019). Options proposed in the literature include adiabatic mirroring (i.e. trapping; Jones, Birmingham & Kaiser 1978; Felice & Kulsrud 2001) and resonance broadening (Dupree 1966), among others. We see strong evidence for such trapping or mirroring in our simulations, which we detail in § 6.4.1. We then continue with a discussion of the repercussions of limiting our heat flux to only passing electrons (§ 6.4.2). Finally, we end § 6.4 with a discussion of how to improve our resonant scattering frequency calculations in § 6.4.3 and a discussion of resonance broadening in § 6.4.4.

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 (2001), 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 (2001) 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 1992),

(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. 1995). 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 1992), where (Roberg-Clark et al. 2016; Cai, Wu & Tao 2020)

(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. 2018) 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 1966)

(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. 1992). Particles can therefore only diffuse across discrete wave modes when the wave amplitudes are large enough that their resonances overlap (Chirikov 1960). Earlier work on the HWI by Roberg-Clark et al. (2016) 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. (2016), 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.

6.5 Alternative transport regimes

In the ICM the saturated field amplitude $\delta B/B_0$ is expected to be ${\sim }3\times 10^{-7}$ at the largest temperature-gradient length scales, leaving $\omega _{b}/\varOmega _e\sim 5\times 10^{-4}$ and $\xi _\textrm {crit}\simeq 3\times 10^{-7}$. Thus, the effect of trapping is expected to be small. In reality, there are Coulomb collisions in the ICM that scatter electrons at a rate ${\sim }\nu _{ei}$, which may help regulate the heat flux. Let us explore this possibility. Coulomb collisions in the ICM correspond to a mean free path

(6.11)\begin{equation} \lambda_{e,{\rm ICM}}\simeq 680 \left(\frac{T_e}{6\,{\rm keV}}\right)^2\left(\frac{10^{{-}2}\,{\rm cm}^{{-}3}}{n_e}\right)\,{\rm pc}. \end{equation}

For temperature-gradient length scales shorter than $\beta _e\lambda _{e,\textrm {ICM}}$, the Coulomb-collisional heat flux will exceed the HWI threshold and the instability will be active. If the HWI is only able to scatter resonant electrons, then there will be a range of $|v_{\parallel }/v_{\textrm {th}e}|\ll 1$ where scattering is dominated by the Coulomb collisions. However, by definition, the Coulomb scattering frequency is not fast enough to stabilize the HWI, so it is unlikely that background Coulomb collisions will cause the HWI to saturate in the manner described and observed in this work. Except perhaps in the case of a specific spectral index for magnetic-field fluctuations or for temperature-gradient length scales near the Coulomb mean free path, the HWI fluctuations must continue to grow until the heat flux reaches the marginally stable value.

Another possibility is that large-scale turbulence traps electrons and isotropizes the distribution at small parallel velocities (for a discussion of this effect in the cosmic-ray literature, see, e.g. Yan & Lazarian 2008; Xu & Lazarian 2018). Assuming the wavelengths of the turbulent fluctuations are significantly smaller than the temperature-gradient length scale, the problem here is in principle the same as for Coulomb collisions. Either the trapping is strong enough to suppress the heat flux below the HWI threshold, resulting in no instability, or the HWI is active and, if only able to scatter resonantly, must grow until the heat flux is regulated to the threshold. Whether it is realistic to expect the system to adapt to some hybrid state with Coulomb collisions or large-scale turbulence – and even what such a state would look like – remains to be seen. It is worth noting that our PIC simulations do have noise that is analogous to Coulomb collisions resulting from a finite number of PIC particles, which clearly affected the results of our quasi-linear model that was computed from the magnetic spectrum measured in our simulations. When we introduced a quasi-linear model derived from a model spectrum that does not incorporate PIC noise, the heat flux implied was significantly larger than that implied by the operator using the measured spectrum. However, because the PIC noise does not scale proportionally to $\beta _e\rho _e/L_T$, the quasi-linear model incorporating the noise did not imply heat fluxes that were consistent with HWI saturation.

6.6 Comparison to model in Drake et al. (2021)

Drake et al. (2021) proposed a model collision operator for the saturated HWI similar to the one presented in § 6.2, with the same differential form of the operator and a scattering frequency

(6.12)\begin{equation} \nu_{\rm Drake21} = 0.1\varOmega_e \left(\frac{\delta B}{B_0}\right)^2\left(\frac{v}{v_{\text{th}e}}\right)^{4/3}. \end{equation}

The power-law index is similar to the pitch-angle-average results from §§ 5.2.3 and 5.3.4 for $v/v_{\textrm {th}e}\gtrsim 1$; however, this result is a coincidence. Drake et al. (2021) argues for a $4/3$ power law by invoking quasi-linear scattering off of an assumed $k^{-7/3}$ spectrum of electron-magnetohydrodynamic turbulence; our results are direct calculations from the simulation. The model (6.12) also contains no $\xi$ dependence, as Drake et al. (2021) argued the whistlers ‘should be able to scatter through the whole range of pitch angle$\ldots$ because of the multiple anomalous resonances’ under a broad spectrum of oblique waves. While we have found that the anomalous cyclotron resonances play a significant part in the saturation of the instability, their contribution to the scattering rate – according to quasi-linear theory and our Fokker–Planck measurements at the largest $L_T/\rho _{e,0}$ – is a factor ${\sim }10$ smaller than from the usual resonances. We also find that the ‘anomalous’ resonances on their own cannot be responsible for scattering electrons across all pitch angles because of the issues mentioned in § 6.4.

7 Conclusion

Making use of an ensemble of electromagnetic PIC simulations across different scale separations $L_T/\rho _e$ and plasma $\beta _e$ parameters, we have demonstrated the ability of the HWI to limit the parallel electron heat flux to quasi-linear threshold levels ${\propto }\beta ^{-1}_e$ in high-$\beta$, weakly collisional, magnetized plasma. We also show this result to be robust with respect to the initial hydrostatic equilibrium. To investigate in detail how the saturated state of the HWI regulates the heat flux, we have utilized three different methods to quantitatively characterize the effects of the HWI fluctuations on the velocity-space dynamics of the electrons.

The first method assumes that the whistlers pitch-angle scatter electrons in a frame travelling with the whistler phase velocity and leverages a Chapman–Enskog expansion to solve for the scattering frequency given the electron distribution function measured in the simulations. Although we found this method to be too noisy to resolve the full pitch-angle dependence of the scattering frequency, we were able to resolve the speed dependence of its pitch-angle average. The second method is derived from the electromagnetic quasi-linear operator, which takes as its input the spectrum of magnetic-field fluctuations measured in our simulations. The resulting operator was rigorously shown to converge – in the limit of slow whistler phase velocity – to a pitch-angle-scattering operator expressed in the frame of the phase velocity. This operator makes evident the previously known fact that oblique whistlers are necessary for heat-flux-limiting interactions to occur with electrons that travel down the temperature gradient. Using this operator, we also showed that the implied scattering frequency of these electrons is an order of magnitude less than for electrons that travel up the temperature gradient. This asymmetry in the scattering frequency directly follows from the ratio of energies in the left- and right-handed components of oblique whistler waves, which resonate with electrons moving down and up the temperature gradient, respectively. The third and final method relied on the calculation of Fokker–Planck jump moments in velocity space from tracked particle data. We showed not only how such jump moments are to be obtained in curvilinear coordinates, but also how to obtain drag and diffusion coefficients from said jump moments by using a jump interval that is self-consistently obtained from the particle statistics.

The pitch-angle-scattering frequency from all of our methods show two power laws in $v$: one decreasing as a function of $v$ for subthermal electrons, due to transit-time damping at $|v_{\parallel }/v_{\textrm {th}e}|\ll 1$; and the other increasing for superthermal electrons, due to cyclotron-resonant scattering at $|v_{\parallel }/v_{\textrm {th}e}|\sim 1$. Unfortunately, the pitch-angle dependence of the method leveraging the Chapman–Enskog expansion was unable to be resolved due to finite-particle-number noise, and our quasi-linear and Fokker–Planck methods produced scattering frequencies that differed significantly from one another. The quasi-linear scattering frequency exhibits steep power laws in $|v_{\parallel }/v_{\textrm {th}e}|$, which results from the steep spectral slope of the magnetic energy. The Fokker–Planck scattering frequency is flat in comparison, only showing slight asymmetry in pitch angle at the smallest $\beta _{e,0}\rho _{e,0}/L_T$. The heat flux implied by the Fokker–Planck operator matches well with our simulations, although there is some random variation in the result that suggests our method could be improved. The quasi-linear operator, however, implies a heat flux that varies as a function of plasma $\beta$ in a way that is inconsistent with the HWI marginal stability criterion, due to transit-time damping for electrons with $|v_{\parallel }|/v_{\textrm {th}e}\ll 1$. We also introduce a simple model, which is motivated by the quasi-linear model but includes a model whistler spectrum. Our simple model implies heat fluxes that are qualitatively similar to the quasi-linear model, except that the implied heat flux is even higher in our simple model. Because a diffusive heat flux is proportional to the inverse of the scattering frequency, it is therefore sensitive to where $\nu$ is the smallest. The increased heat flux implied by the simple model suggests that the numerical quasi-linear result is sensitive to the level of PIC noise in the runs, which is not incorporated into our model spectrum. Even though our numerical quasi-linear operator does imply heat fluxes that are consistent with simulations producing small fluctuation amplitudes, these issues call into question whether a resonant quasi-linear approach can explain physically the saturation of the HWI.

Ultimately, though, it is the quasi-linear cyclotron resonance that explains why the marginal whistler scattering rate $\nu /\varOmega _e\sim \beta _e\rho _e/L_T$ is also proportional to the square of the whistler fluctuation amplitude, i.e. $\nu /\varOmega _e\sim \langle \delta B^2\rangle /B_0^2$, which we clearly show to be the case in our simulations. But there is the rub: electrons with $|v_{\parallel }/v_{\textrm {th}e}|\ll 1$ are cyclotron resonant with $k_{\parallel }\rho _e\gg 1$ whistler waves, which have vanishingly little power at high parallel wavenumbers. The problem of how exactly these $|v_{\parallel }/v_{\textrm {th}e}|\ll 1$ electrons are scattered in a way that is consistent with the simple cyclotron-resonant quasi-linear picture of the HWI is a vexing one. This is a variation on the well-known $90^{\circ }$ scattering problem in the cosmic-ray community, albeit with tighter constraints than is usually found in the literature. We commented on various ways to improve the quasi-linear scattering rate at low $|v_{\parallel }/v_{\textrm {th}e}|$, including using a more accurate treatment of the dispersion relation and incorporating various schemes for resonance broadening. The former would require an in-depth method to solve the plasma dispersion function for our distribution function and is beyond the scope of this work; the latter runs into the problem that resonances are broadened in a way that depends on the turbulent fluctuation amplitude. While detailed calculations of resonance-broadened scattering frequencies also lie beyond the scope of this work, we argue that, in general, such scattering frequencies will not be proportional to $\langle \delta B^2\rangle /B_0^2$ and will violate the basic picture of the HWI to which we are constrained by simulation results.

The best model we can currently produce is a simple quasi-linear scattering frequency derived from a model whistler spectrum, but with an empirically motivated floor that is proportional to the marginal whistler scattering frequency ${\propto }\beta _e\rho _e/L_T$. It is no surprise that this model does imply a heat flux consistent with our simulations; however, we stress that we do not have an explanation for the physical process that underlies such a floor, other than to say that whatever that process is, it must result in the same scaling. This conclusion limits the possibility that either Coulomb collisions or large-scale turbulence could play a role in HWI saturation. For either process to help saturate the instability effectively, they would need to scatter electrons at or near the instability threshold value; but, this also means that the plasma would be stable to the HWI in the first place and so the scheme is not logically consistent. Clearly, for the HWI to regulate the electron heat flux self-consistently, therefore, the effective collision operator must look closer to our relatively flat Fokker–Planck operator or to the model from Drake et al. (2021). Unfortunately, a rigorous statement as to how this operator works in detail is beyond the scope of the current work.

While our three methods to obtain an effective HWI collision operator did not converge numerically to an unambiguous result, we hope the procedures outlined in this paper find utility in characterizing the phase-space dynamics of other kinetic instabilities in future studies. In fact, a concurrent work by some of the authors uses similar techniques to this work in the characterization of a pressure-anisotropy-driven firehose instability in hybrid-PIC simulations. That work saw more success in particular with the method that leverages the Chapman–Enskog expansion, owing to a much better resolved particle distribution function. That study also points out limitations with the Fokker–Planck method in regards to the locality of individual particle jumps; a curious reader is encouraged to refer to Bott et al. (2024, in preparation) for an in-depth discussion. We also plan to apply many of the lessons learned here in the study of a corresponding ion-heat-flux instability, discovered by Bott et al. (2024), using hybrid-PIC simulations. The utility of these methods to obtain effective collision operators lies in their distillation of complex kinetic physics into relatively simple expressions, which can then be used in a small-parameter expansion for an even simpler fluid-like closure of the kinetic equations. Such a closure would be much less computationally expensive to simulate and could be incorporated into Braginskii-magnetohydrodynamic simulations of the ICM (e.g. Berlok et al. 2021; Perrone, Berlok & Pfrommer 2024) and simple models of the solar wind or low-luminosity black-hole accretion flows.