Collisional whistler instability and electron temperature staircase in inhomogeneous plasma

Abstract

High-beta magnetised plasmas often exhibit anomalously structured temperature profiles, as seen from galaxy cluster observations and recent experiments. It is well known that when such plasmas are collisionless, temperature gradients along the magnetic field can excite whistler waves that efficiently scatter electrons to limit their heat transport. Only recently has it been shown that parallel temperature gradients can excite whistler waves also in collisional plasmas. Here, we develop a Wigner–Moyal theory for the collisional whistler instability starting from Braginskii-like fluid equations in a slab geometry. This formalism is necessary because, for a large region in parameter space, the fastest-growing whistler waves have wavelengths comparable to the background temperature gradients. We find additional damping terms in the expression for the instability growth rate involving inhomogeneous Nernst advection and resistivity. They (i) enable whistler waves to re-arrange the electron temperature profile via growth, propagation and subsequent dissipation, and (ii) allow non-constant temperature profiles to exist stably. For high-beta plasmas, the marginally stable solutions take the form of a temperature staircase along the magnetic field lines. The electron heat flux can also be suppressed by the Ettingshausen effect when the whistler intensity profile is sufficiently peaked and oriented opposite the background temperature gradient. This mechanism allows cold fronts without magnetic draping, might reduce parallel heat losses in inertial fusion experiments and generally demonstrates that whistler waves can regulate transport even in the collisional limit.

1. Introduction

X-ray observations of the diffuse plasma residing within galaxy clusters have revealed intricately structured temperature fields whose sharp gradients are inferred to have persisted far longer than classical transport theory predicts (Peterson & Fabian 2006; Markevitch & Vikhlinin 2007). Recent experiments at the National Ignition Facility with laser-produced hot, magnetised, high-beta plasmas also feature such anomalously structured temperature fields (Meinecke et al. 2022). Such structures would require a significantly reduced level of electron heat transport, which might be due to plasma microinstabilities. Indeed, it is now well established (Kunz et al. 2022) that the macroscopic transport properties of high-beta magnetised plasmas are significantly modified by small-scale instabilities such as firehose (Rosin et al. 2011; Kunz et al. 2014), mirror (Kunz et al. 2014; Komarov et al. 2016) or heat-flux whistler instabilities (Levinson & Eichler 1992; Pistinner & Eichler 1998; Komarov et al. 2018; Roberg-Clark et al. 2018; Drake et al. 2021; Yerger et al. 2025).

The heat-flux whistler instability is a particularly promising candidate because it is the fastest of all possible instabilities in the relevant parameter regime (Bott, Cowley & Schekochihin 2024) and quickly limits the parallel heat flux to a marginal value via electron scattering. However, this mechanism requires a resonant interaction with heat-carrying electrons that can be disrupted by collisions. This is not a problem for astrophysical plasmas whose magnetisation is typically of order $\mathcal {M} \doteq \Omega \tau _{ei} \sim 10^{12}$ , where $\Omega$ is the electron cyclotron frequency and $\tau _{ei}$ is the electron–ion collision time. However, it is not understood whether such an instability could persist in the more collisional laboratory analogues, whose magnetisation range is $\mathcal {M} \sim 10^{-2} {-} 10^{2}$ (Meinecke et al. 2022).

Recently a collisional mechanism for exciting whistler waves via anisotropic friction forces was identified (Bell et al. 2020). However, this initial analysis was restricted to the short-wavelength geometrical-optics limit, which prevented it from correctly describing all aspects of the long-wavelength fluid limit. This is problematic because in a typical laser-plasma experiment (Meinecke et al. 2022), the long-wavelength modes will be (i) the most easily observable in diagnostics and (ii) the first modes excited as the plasma heats up, and therefore the modes most capable of subsequently manipulating the plasma.

Here, we remove this shortcoming by deriving the Wigner–Moyal equations that govern the collisional whistler instability in a slab geometry, similar to what has been done in modelling drift-wave turbulence beyond the geometrical-optics approximation (Ruiz et al. 2016; Zhu et al. 2018; Tsiolis, Zhou & Dodin 2020). We find additional terms in the instability dispersion relation and growth rate due to gradients in the background plasma. We proceed to show that these additional terms can actually stabilise a non-constant temperature profile along the magnetic field, which is not possible if only the geometrical-optics approximation is used. Equivalently, these additional terms cause the instability to be damped at low temperatures, providing a mechanism for the instability to re-arrange the temperature profile into a marginally stable state via excitation at high temperature, propagation down the temperature gradient via Nernst advection and subsequent damping at low temperature. We proceed to show that the marginally stable temperature profile generically takes the form of a staircase where isothermal regions are insulated from each other by abrupt jumps in the temperature, which occur at the zeros of a certain function comprising an intricate combination of magnetic transport coefficients. These staircases can be in pressure balance, resembling the ubiquitous cold fronts in galaxy clusters (Markevitch & Vikhlinin 2007) but with the magnetic field no longer required to drape the front.

We then derive the back-reaction of the instability on the background temperature profile and show that, in the initial stages of the instability, the frictional work done by the instability actually cools the background plasma instead of heating it (as demanded by energy conservation). Moreover, when the gradient of the unstable whistlers’ amplitude is anti-aligned with the background temperature gradient along the magnetic field, the parallel heat flux can be reduced via the Ettingshausen effect, although this is more difficult to achieve in high-beta plasmas. If this mechanism can be reliably engineered, however, it might allow higher hotspot temperatures to be achieved in magnetised inertial fusion, since the parallel heat flux is the present limiting factor (Walsh et al. 2022).

2. Summary

Here, we first provide an executive summary that highlights the main definitions, discussions and results for each section, serving as an overall roadmap of the paper that can be consulted later for easy reference.

In § 3, the governing electron magnetohydrodynamics (MHD) equations and slab geometry are introduced, leading to the set (3.6). Fundamentally, the electron-MHD limit in slab geometry results in the simplification that only the electron temperature and the perpendicular components of the magnetic field (with respect to the single direction of inhomogeneity) have non-trivial time evolution; the electron density and the parallel magnetic field both remain constant in time. The perpendicular magnetic field is then expressed in the diagonalising eigenbasis for its evolution equation, leading to the simplified description in terms of mode amplitudes given in (3.12) and (3.13). This is all for a general, unspecified friction and heat flux; the remaining parts of this section (and the remainder of the paper) specialise to when the friction and heat flux are determined by the standard Chapman–Enskog expressions (3.14). The nine transport coefficients ( $\alpha _\parallel$ , $\alpha _\perp$ , $\alpha _\wedge$ , $\beta _\parallel$ , $\beta _\perp$ , $\beta _\wedge$ , $\kappa _\parallel$ , $\kappa _\perp$ , $\kappa _\wedge$ ) are all generally functions of the dimensionless magnetisation parameter $\mathcal {M}$ defined in (3.17). Ultimately, after inserting the same eigenmode decomposition into the expressions for the friction and heat flux, one arrives at the final set of working equations: the magnetic-field mode amplitudes are governed by (3.22) and the temperature evolution is governed by (3.26). Up to that point, all manipulations of the initial equations are exact and all non-linearities are retained.

In § 4, the dynamical equations are linearised about the small transverse field amplitude to obtain their dispersion relation and growth rate; since the temperature evolution involves terms quadratic or more in the mode amplitude, the two dynamical equations decouple in this limit. The temperature profiles are treated as fixed in time from the perspective of the magnetic field fluctuations, even though heat conduction is still present in the lowest-order temperature evolution equation; the validity of this approximation is outlined in § 6. To maintain an exact treatment (within the linear regime), the dispersion relation and growth rate are obtained as the Hamiltonian and non-Hamiltonian components of the generalised wave-kinetic equation (4.8) that governs the magnetic fluctuations. Specifically, the Hermitian and anti-Hermitian parts of the dispersion relation are given in (4.9), with the constituent terms defined in (4.1). This method of deriving the dispersion relation is chosen because it does not rely on making a short-wavelength approximation. The remainder of § 4 is therefore dedicated to defining the parameter space in which such a general treatment is necessary; it is found that the wavelength for the maximally growing mode becomes comparable to the background medium inhomogeneity length scale in the regime $\mathcal {M}^3 \beta _0 \ll 1$ , where $\beta _0$ is the plasma beta defined in (4.14).

In § 5, the expression for the growth rate (4.9b ) is examined in more detail to understand the condition for instability (5.2). Allowing for wavelengths comparable to the inhomogeneity scale of the medium introduces two additional stabilising mechanisms: one related to the gradient of the Nernst velocity that advects the perturbations and one to the curvature of the plasma resistivity that diffuses the perturbations. Simple heuristic descriptions of these two stabilising mechanisms are then presented briefly. The section concludes with a specific discussion for the simple case (5.5) when the plasma is in pressure balance with a linear temperature profile; in this case, it is shown that the additional stabilising terms cause a transition from unstable to stable behaviour analogous to the transition of the short-wavelength approximation from being valid to being violated discussed in § 4.

In § 6, the validity for the approximation that the temperature profile remain fixed in time for the linear stability analysis is assessed by comparing the instability growth rate (6.1) with the diffusion time (6.4) of the background temperature profile. This analysis accommodates arbitrarily shaped plasma profiles, with local length scales and curvatures defined in (6.2). For convenience, the growth rate (6.1) is also split into terms driven by temperature gradients, temperature curvatures, mixed temperature and density gradients, density gradients and density curvatures; the coefficient functions, defined in (6.3), are grouped together respectively as $F_{T,1}$ , $F_{T,2}$ , $F_{T,n}$ , $F_{n,1}$ and $F_{n,2}$ . Generally, it is found that the instability’s growth rate (6.5) normalised by the dynamical evolution time is only greater than unity in a sub-region of parameter space where both $\mathcal {M}$ and $\beta _0$ are large; this is in contrast to the short-wavelength approximation that predicts the instability growth to be dynamically relevant everywhere (cf. (6.9)). However, this result is nuanced because only regions where the instability is present are considered and, in fact, in much of the excluded parameter space there is no instability and perturbations are instead strongly damped (at a rate much larger than the diffusion time). The strong damping is due to the additional stabilising mechanisms that occur when the wavelengths become comparable to the inhomogeneity scale discussed in § 5.

In § 7, the marginally stable temperature profile (satisfying the condition for zero growth rate everywhere) is derived. Due to the additional stabilising terms in the expression for the growth rate, the solution involves non-trivial global structure. As discussed first in the general case, the governing equation (7.8) for the marginally stable $\mathcal {M}$ profile (a proxy for temperature) has the general structure of a nonlinear boundary-layer differential equation (7.12), due to the singular nature as $\beta _0$ becomes large of the coefficient $g(\mathcal {M})$ defined in (7.9) and (7.5). A variable transformation is performed to recast the nonlinear differential equation into a linear one that is then readily solved by successive integrations, yielding the solution (7.11) for the inverse function of the spatial profile $\mathcal {M}(z)$ . This solution is then shown to exhibit a staircase structure as a generic feature due to the boundary layers. The remainder of the section specialises the general theory to two situations: plasmas with constant pressure and plasmas with constant density. Both cases are qualitatively similar when viewed in a parameter space consisting of $\mathcal {M}$ and some measure of the plasma beta; for the constant-pressure case, this is simply $\beta _0$ , but for the constant density case, an effective plasma beta is defined in (7.29). The staircase is shown to have a single step, in rough correspondence to the transition between (dynamically relevant) instability growth and damping outlined in the previous § 6. The predicted staircase feature of the solution is then verified with direct numerical simulation of the governing differential equation along with a closed-form analytical approximation described in (7.23).

In § 8, the linear analysis is extended to a quasilinear study of the electron temperature response by restoring the lowest-order (quadratic) fluctuation terms in the temperature-evolution equation (8.1). Two different groupings of the various terms are presented depending on the physics to be emphasised, with definitions provided in (8.2). First, it is discussed how (as required by energy conservation) the growth of unstable perturbations also corresponds to a net cooling effect of the friction forces when the condition (8.6) is satisfied. Next, the total heat flux (including Ettingshausen terms generated by the instability) is analysed and found to be reduced compared with the standard conductive heat flux when the condition (8.9) is satisfied, although the reduction is generally modest. Lastly, the heat-flux reduction is analysed when the temperature has the marginally stable staircase profile derived in § 7. The extreme situation in which the heat flux is completely suppressed for the marginal-stability profiles is then considered to derive required conditions on the instability intensity profiles; the resulting expression (8.15) suggests that such high degrees of heat-flux suppression are unlikely to occur for high-beta plasmas within the validity of the quasilinear fluid model adopted here.

Finally, in § 9, the main results are summarised and directions for future work are outlined. Auxiliary calculations and review material are presented in the Appendices.

3. Governing fluid equations

3.1. Extended electron MHD equations in slab geometry

We are interested in the dynamics of electromagnetic oscillations in the whistler-frequency range in a collisional plasma. To allow an analytical description, let us consider for simplicity the extended electron MHD equations for a Lorentz plasma with stationary ions and isotropic pressure tensor:

(3.1a) \begin{align} \partial _{t} n &= 0, \end{align}

(3.1b) \begin{align} \partial _{t} {\textbf{B}} &= - c \nabla \times \left [ \frac { (\nabla \times {\textbf{B}}) \times {\textbf{B}}}{4 \pi e n} - \frac {\nabla (nT)}{e n } + \frac {{\textbf{R}}}{e n} \right ], \end{align}

(3.1c) \begin{align} \frac {3}{2} n \, \partial _{t} T &= \frac {c}{4 \pi n e} (\nabla \times {\textbf{B}}) \cdot \left ( \frac {3}{2} n \nabla T - T \nabla n + {\textbf{R}} \right ) - \nabla \cdot {\textbf{q}} . \\[6pt] \nonumber \end{align}

Here, all symbols have their usual meaning: $n$ and $T$ are the electron density and temperature, respectively, ${\textbf{B}}$ is the magnetic field, $c$ is the speed of light in vacuum, $e \gt 0$ is the absolute value of the electron charge, ${\textbf{R}}$ is the frictional force experienced by the electron fluid due to pitch-angle-scattering collisions with ions and ${\textbf{q}}$ is the electron heat flux. Expressions for ${\textbf{R}}$ and ${\textbf{q}}$ will be provided later in this section. Physically, the term within parenthesis in (3.1c ) comprises three distinct contributions: in order of appearance, they are (i) advection of temperature with the electron flow; (ii) compressional heating and (iii) frictional heating that can be related to Joule heating by replacing ${\textbf{R}}$ with ${\textbf{E}}$ via Ohm’s law – indeed, the terms inside square brackets in (3.1b ) are precisely ${\textbf{E}}$ .

The simplest set-up exhibiting the collisional whistler instability has a temperature gradient, density gradient, and wavevector of unstable perturbations all aligned with a mean background magnetic field (Bell et al. 2020). Hence, it can be adequately described by a one-dimensional ( $1$ -D) slab model in which the total magnetic field is given by

(3.2) \begin{align} {\textbf{B}}(t, z) = \begin{pmatrix} \widetilde {B}_x(t, z), \widetilde {B}_y(t, z), B_z(t) \end{pmatrix}^\intercal, \end{align}

with $B_z \gt 0$ being the mean field and $\widetilde {B}_x$ and $\widetilde {B}_y$ being fluctuating quantities associated with the instability. We also take $n$ and $T$ to be functions of $z$ and $t$ only. As discussed in Appendix A, this constraint means that the collisional whistler instability has no associated density or temperature fluctuations at the fundamental frequency. Since $B_z$ is independent of $z$ , one has automatically

(3.3) \begin{align} \nabla \cdot {\textbf{B}}(t, z) = 0 . \end{align}

One also has

(3.4) \begin{align} \nabla \times {\textbf{B}}(t, z) = \begin{pmatrix} - \partial _{z} \widetilde {B}_y(t, z) \\[4pt] \partial _{z} \widetilde {B}_x(t, z) \\[4pt] 0 \end{pmatrix} = \left(\begin{array}{c@{\quad}c} - \mathsf {J} & 0 \\ 0 & 0 \end{array} \right) \partial _{z} {\textbf{B}}, \end{align}

where we have introduced the skew-symmetric matrix

(3.5) \begin{align} \mathsf {J} = \left(\begin{array}{c@{\quad}c} 0 & 1 \\ -1 & 0 \end{array} \right) . \end{align}

Using the assumed form for the dynamical variables, the extended electron MHD equations (3.1) become

(3.6a) \begin{align} \partial _{t} n &= 0, \end{align}

(3.6b) \begin{align} \partial _{t} B_z &= 0, \end{align}

(3.6c) \begin{align} \partial _{t} \widetilde {\textbf{B}}_\perp &= \partial _{z} \left ( \frac { \Omega c^2}{\omega _p^2 } \mathsf {J} \, \partial _{z} \widetilde {\textbf{B}}_\perp + \frac {c}{e n} \mathsf {J} \, {\textbf{R}}_\perp \right ), \end{align}

(3.6d) \begin{align} \frac {3}{2} n \, \partial _{t} T &= - \frac {\Omega c^2}{\omega _p^2 B_z} {\textbf{R}}_\perp \cdot \mathsf {J} \cdot \partial _{z} \widetilde {\textbf{B}}_\perp - \partial _{z} q_z, \\[6pt] \nonumber \end{align}

where we have defined the local plasma frequency and the mean electron cyclotron frequency respectively as

(3.7) \begin{align} \omega _p^2(z) = \frac {4 \pi e^2}{m} n(z), \quad \Omega = \frac {e B_z}{m c} . \end{align}

Since the evolution equations for $n$ and $B_z$ are trivial, we shall omit them in the following analysis.

3.2. Eigenbasis projection

Further simplifications can be obtained by expanding $\widetilde {\textbf{B}}_\perp$ and ${\textbf{R}}_\perp$ onto the eigenbasis of $\mathsf {J}$ , viz. the circular polarisation vectors. These eigenvectors and their eigenvalues are given by

(3.8) \begin{align} {\check {e}}_\pm = \frac {1}{\sqrt {2} } \begin{pmatrix} 1 \\ \pm i \end{pmatrix}, \quad \mathsf {J} \, {\check {e}}_\pm = \pm i \, {\check {e}}_\pm, \end{align}

and satisfy the orthogonality condition

(3.9) \begin{align} {\check {e}}_\pm ^* \cdot {\check {e}}_\pm = 1, \quad {\check {e}}_\mp ^* \cdot {\check {e}}_\pm = 0 . \end{align}

Since $\widetilde {\textbf{B}}_\perp$ and ${\textbf{R}}_\perp$ are both real-valued vectors and since ${\check {e}}_+^* = {\check {e}}_-$ , the eigenbasis expansion takes the form

(3.10) \begin{align} \widetilde {\textbf{B}}_\perp = \epsilon B_z \frac { \psi \, {\check {e}}_+ + \psi ^* \, {\check {e}}_+^*}{2}, \quad {\textbf{R}}_\perp = \frac { \xi \, {\check {e}}_+ + \xi ^* \, {\check {e}}_+^*}{2}, \end{align}

where $\psi$ and $\xi$ are the complex scalar wavefunctions of $\widetilde {\textbf{B}}_\perp$ and ${\textbf{R}}_\perp$ , respectively, and $\epsilon \gt 0$ is a constant (assumed small) that parametrises the relative size of the fluctuations. Later, it will be shown that $\xi$ is of order $\epsilon$ as well. Note that since

(3.11) \begin{align} \begin{pmatrix} \widetilde {B}_x \\ \widetilde {B}_y \end{pmatrix} = \frac {\epsilon B_z}{\sqrt {2}} \begin{pmatrix} {\textrm {Re}} \, \psi \\ - \textrm {Im} \, \psi \end{pmatrix}, \quad \begin{pmatrix} R_x \\ R_y \end{pmatrix} = \frac {1}{\sqrt {2}} \begin{pmatrix} {\textrm {Re}} \, \xi \\ - \textrm {Im} \, \xi \end{pmatrix}, \end{align}

by introducing $\psi$ and $\xi$ , we have essentially traded two real-valued degrees of freedom for a single complex-valued degree of freedom.

Since ${\check {e}}_+$ is independent of $t$ and $z$ , one can readily show using orthogonality that the real-vector-valued evolution equation (3.6c ) for $\widetilde {\textbf{B}}_\perp$ is equivalent to the following complex-scalar-valued evolution equation for $\psi$ :

(3.12) \begin{align} i \partial _{t} \psi = - \partial _{z} \left ( \frac { \Omega c^2}{\omega _p^2} \partial _{z} \psi + \frac {c}{e n} \frac {\xi }{\epsilon B_z} \right ) . \end{align}

Similarly, the temperature equation (3.6d ) takes the form

(3.13) \begin{align} \frac {3}{2} n \, \partial _{t} T = \epsilon \frac {\Omega c^2}{\omega _p^2} \frac {\textrm {Im}\left ( \xi ^* \partial _{z} \psi \right )}{2 } - \partial _{z} q_z . \end{align}

3.3. Chapman–Enskog friction coefficients

Equations (3.12) and (3.13) are valid for any friction force, allowing one to study driven systems. For undriven systems, the friction force is determined by the plasma fluid variables themselves according to some closure. A common closure that we shall adopt here is provided by the Chapman–Enskog method (Helander & Sigmar 2002; Bott et al. 2024), which for the Lorentz collision operator, yields the following expressions for the friction force and for the heat flux (Epperlein 1984):

(3.14) \begin{align} {\textbf{R}} = \frac {n e c }{ \omega _p^2 \tau _{ei}} \mathsf {M}_\alpha \cdot \nabla \times {\textbf{B}} - n \mathsf {M}_\beta \cdot \nabla T, \quad {\textbf{q}} = - n \tau _{ei} v_t^2 \mathsf {M}_\kappa \cdot \nabla T - \frac {\Omega c^2}{\omega _p^2} \frac {n T}{B_z} \, \mathsf {M}_\beta \cdot \nabla \times {\textbf{B}}, \end{align}

where $\tau _{ei}$ is the electron–ion collision time (Epperlein & Haines 1986; Helander & Sigmar 2002) and $v_t$ is the thermal speed, defined respectively as

(3.15) \begin{align} \tau _{ei} = \frac {12 \pi ^2}{ \sqrt {2 \pi } } \frac { n v_t^3 }{ Z \omega _p^4 \log \Lambda }, \quad v_t = \sqrt {\frac {T}{m}} . \end{align}

The dimensionless resistivity ( $\alpha$ ), thermoelectric ( $\beta$ ) and conductivity ( $\kappa$ ) matrices are anisotropic with respect to the magnetic field, taking the form

(3.16) \begin{align} \mathsf {M}_\sigma = \sigma _\perp (\mathcal {M})\mathsf {I}_{3} + \Delta _\sigma (\mathcal {M}) \frac {{\textbf{B}} {\textbf{B}}}{|B|^2} \pm \sigma _\wedge (\mathcal {M}) \frac {\mathsf {B}_\wedge }{|B|}, \quad \sigma = \alpha, \beta, \kappa, \end{align}

where $\mathsf {B}_\wedge$ denotes the skew-symmetric hat-map matrix that enacts the cross-product $\mathsf {B}_\wedge \cdot {\textbf{v}} = {\textbf{B}} \times {\textbf{v}}$ (Zhang, Fu & Qin 2020) and $\Delta _\sigma \doteq \sigma _\parallel - \sigma _\perp$ is the anisotropy measure. In the last term, the minus sign applies to $\alpha _\wedge$ only. Also note that $\sigma _\perp$ and $\sigma _\wedge$ are both positive, but $\Delta _\alpha \le 0$ , while $\Delta _\beta \ge 0$ and $\Delta _\kappa \ge 0$ . The various transport coefficients are functions purely of the magnetisation parameter

(3.17) \begin{align} \mathcal {M} \doteq \frac {|B|}{B_z} \Omega \tau _{ei} &= \Omega \tau _{ei} \sqrt { 1 + \frac {\epsilon ^2}{2} |\psi |^2 } \nonumber \\ &\equiv \frac {3}{ 4 \sqrt {2 \pi } } \frac { \Omega \sqrt {m} }{ Z e^4 \log \Lambda } \frac {T^{3/2} }{n } \sqrt { 1 + \frac {\epsilon ^2}{2} |\psi |^2 }. \end{align}

In the remainder of the analysis, we shall use the rational interpolants for the various transport coefficients as functions of $\mathcal {M}$ developed by Lopez (2024). Their limiting forms as $\mathcal {M} \to 0$ and $\mathcal {M} \to \infty$ are listed in Appendix B.

Lastly, note that the slab geometry allows the relevant matrices to be constructed explicitly with a relatively simple form. Since

(3.18) \begin{align} \mathsf {B}_\wedge = \left(\begin{array}{c@{\quad}c@{\quad}c} 0 & - B_z & \widetilde {B}_y \\[3pt] B_z & 0 & - \widetilde {B}_x \\[3pt] - \widetilde {B}_y & \widetilde {B}_x & 0 \end{array} \right), \end{align}

the anisotropic transport matrices are given as

(3.19) \begin{align} \mathsf {M}_\sigma = \frac {1}{|B|^2} \left(\begin{array}{c@{\quad}c@{\quad}c} \Delta _\sigma \widetilde {B}_x^2 + \sigma _\perp |B|^2 & \Delta _\sigma \widetilde {B}_x \widetilde {B}_y \mp \sigma _\wedge |B| B_z & \Delta _\sigma \widetilde {B}_x B_z \pm \sigma _\wedge |B| \widetilde {B}_y \\[3pt] \Delta _\sigma \widetilde {B}_x \widetilde {B}_y \pm \sigma _\wedge |B| B_z & \Delta _\sigma \widetilde {B}_y^2 + \sigma _\perp |B|^2 & \Delta _\sigma \widetilde {B}_y B_z \mp \sigma _\wedge |B| \widetilde {B}_x \\[3pt] \Delta _\sigma \widetilde {B}_x B_z \mp \sigma _\wedge |B| \widetilde {B}_y & \Delta _\sigma \widetilde {B}_y B_z \pm \sigma _\wedge |B| \widetilde {B}_x & \Delta _\sigma B_z^2 + \sigma _\perp |B|^2 \end{array} \right) . \end{align}

3.4. Resulting equations for Chapman–Enskog friction forces

From (3.14), the perpendicular component of the Chapman–Enskog friction force in slab geometry is

(3.20) \begin{align} {\textbf{R}}_\perp = \frac {n e c}{\omega _p^2 \tau _{ei}} \left ( \frac {\alpha _\wedge B_z}{|B|} \mathsf {I}_{2} - \alpha _\perp \mathsf {J} - \frac {\Delta _\alpha }{|B|^2} \widetilde {\textbf{B}}_\perp \widetilde {\textbf{B}}_\perp ^\intercal \mathsf {J} \right ) \, \partial _{z} \widetilde {\textbf{B}}_\perp -\frac {n \, \partial _{z}T}{|B|} \left ( \frac {\Delta _\beta B_z}{|B|} \mathsf {I}_{2} + \beta _\wedge \mathsf {J} \right ) \widetilde {\textbf{B}}_\perp . \end{align}

Hence, the complex amplitude $\xi$ is given in terms of $\psi$ as

(3.21) \begin{align} \xi &= \epsilon B_z \frac {n e c}{\omega _p^2 \tau _{ei}} \left ( \frac {\alpha _\wedge B_z}{|B|} - i \alpha _\perp \right ) \partial _{z} \psi \nonumber \\ &\hspace {4mm}- \epsilon B_z \left [ \frac {n \, \partial _{z}T}{|B|} \left ( \frac {\Delta _\beta B_z}{|B|} + i \beta _\wedge \right ) - \frac {\epsilon ^2 B_z^2}{2} \frac {n e c}{\omega _p^2 \tau _{ei}} \frac {\Delta _\alpha }{|B|^2} \textrm {Im} \left ( \psi ^* \partial _{z}\psi \right ) \right ] \psi . \end{align}

Thus, $\xi = O(\epsilon )$ , as promised following (3.10). Hence, (3.12) becomes

(3.22) \begin{align} i \partial _{t} \psi &= - \partial _{z} \left [ \left ( G_d - i \eta \right ) \partial _{z} \psi \vphantom {\frac {}{}} \right ] + \partial _{z} \left [ \left ( u_\gamma - i v_N \right ) \psi \vphantom {\frac {}{}} \right ], \end{align}

where we have defined the following quantities:

(3.23a) \begin{align} G_d &= \frac { \Omega c^2}{\omega _p^2} \left ( 1 + \frac {\alpha _\wedge }{\mathcal {M} } \right ), \end{align}

(3.23b) \begin{align} \eta &= \frac { \Omega c^2}{\omega _p^2} \frac {|B|}{B_z}\frac {\alpha _\perp }{ \mathcal {M} }, \end{align}

(3.23c) \begin{align} u_\gamma &= \frac {\Delta _\beta }{m \Omega } \left ( \frac {B_z}{|B|}\right )^2 \partial _{z}T - \frac {\epsilon ^2}{2} \frac {\Omega c^2}{\omega _p^2} \frac {B_z}{|B|} \frac {\Delta _\alpha }{\mathcal {M}} \textrm {Im} \left ( \psi ^* \partial _{z}\psi \right ), \end{align}

(3.23d) \begin{align} v_N &= - \frac {\beta _\wedge }{m \Omega } \frac {B_z}{|B|} \partial _{z}T . \\[6pt] \nonumber \end{align}

Recall that $|B|$ and $\mathcal {M}$ depend on $|\psi |^2$ , so (3.22) still contains nonlinear effects.

The $z$ -component of Chapman–Enskog heat flux (3.14), $q_z$ , that appears in (3.13) can be shown to take the form

(3.24) \begin{align} q_z &= - \frac {n T \tau _{ei}}{m} \left [ \kappa _{\textrm {eff}} \, \partial _{z} T + \epsilon ^2 \frac {\Omega ^2 m c^2}{2 \omega _p^2} \left ( \frac {\beta _\wedge }{2 \mathcal {M}} \, \partial _{z} |\psi |^2 + \frac {B_z}{|B|} \frac {\Delta _\beta }{\mathcal {M}} \textrm {Im} \left ( \psi ^* \partial _{z}\psi \right ) \right ) \right ], \end{align}

where the effective conductivity is given as

(3.25) \begin{align} \kappa _{\textrm {eff}} = \frac {2 \kappa _\parallel + \epsilon ^2 \kappa _\perp |\psi |^2 }{ 2 + \epsilon ^2 |\psi |^2 } . \end{align}

Hence, (3.13) becomes

(3.26) \begin{align} \frac {3}{2} n \, \partial _{t} T &= \epsilon ^2 \frac {B_z^2}{8\pi } \left [ \eta |\partial _{z} \psi |^2 - u_\gamma \textrm {Im}\left ( \psi ^*\partial _{z} \psi \right ) - \frac {v_N}{2} \partial _{z} |\psi |^2 \right ] - \partial _{z} q_z, \end{align}

with $q_z$ given in (3.24).

In summary, the two main equations we shall be working with in the remainder of this paper are the evolution of the transverse field perturbations (describing the linear instability dynamics) governed by (3.22), and the evolution of the temperature profile (describing the back-reaction of the instability on the equilibrium dynamics) governed by (3.26). More will be said about the physical content of (3.26) in § 8, but as a brief preview, it is worthwhile to highlight now that only the first term within square brackets is manifestly positive (corresponding to resistive heating). The remaining terms, although still related to frictional forces, can take either sign depending on whether the instability is growing or damping, allowing the instability to re-distribute the temperature profile. As we shall discuss in later sections (beginning in the following section), this effect can be particularly prominent in the high-beta, low-magnetisation regime.

4. Dispersion relation, growth rate and breakdown of the short-wavelength approximation

Let us now restrict attention to the linear limit when $\epsilon \to 0$ . To lowest order, (3.26) is decoupled from (3.22), so we shall just consider the dynamics of (3.22) with prescribed stationary temperature and density profiles (the back-reaction of the instability on the temperature profile will be considered in § 8). Moreover, (3.22) maintains the same form when $\epsilon \to 0$ , but the constituent functions (3.23) become simply

(4.1a) \begin{align} G_d &= \frac { \Omega c^2}{\omega _p^2} \left ( 1 + \frac {\alpha _\wedge }{\mathcal {M} } \right ), \end{align}

(4.1b) \begin{align} \eta &= \frac { \Omega c^2}{\omega _p^2} \frac {\alpha _\perp }{ \mathcal {M} }, \end{align}

(4.1c) \begin{align} u_\gamma &= \frac {\Delta _\beta }{m \Omega } \partial _{z}T, \end{align}

(4.1d) \begin{align} v_N &= - \frac {\beta _\wedge }{m \Omega } \partial _{z}T, \\[6pt] \nonumber \end{align}

where $\mathcal {M} = \Omega \tau _{ei}$ . In this limit, the physical interpretation of these quantities becomes clearer: $G_d$ represents the familiar group-velocity dispersion for whistler waves but modified by friction-induced Hall effect (Davies et al. 2021), $\eta$ governs the resistive diffusion of magnetic-field perturbations, $u_\gamma$ is the cross-gradient Nernst advection velocity and $v_N$ is the standard Nernst advection velocity. In particular, the presence of $G_d$ explains our choice to call this instability the ‘collisional whistler instability’: in the absence of friction terms (i.e. with all $\alpha$ and $\beta$ coefficients set to zero), the dispersion relation (4.9a ), which we shall derive shortly, would be identical to the standard whistler dispersion relation in the electron-MHD limit (see, for example, Komarov et al. (2018)).

Figure 1. Plots of the normalised group-velocity dispersion $\widetilde {G_d} = \beta _0G_d/v_t r_L$ (see (4.1a )) for (a) a linear temperature profile (4.2) and (b) a Gaussian temperature profile (4.3). All normalisation quantities are defined with respect to $T_0$ , and $\mathcal {M}_0 = \mathcal {M}(T_0)$ .

Figure 2. Same as figure 1 but for the normalised resistivity $\widetilde {\eta } = \beta _0\eta /v_t r_L$ (see (4.1b )).

Figure 3. Same as figure 1 but for the normalised cross-gradient Nernst velocity $\widetilde {u_\gamma } = Lu_\gamma /r_L v_t$ (see (4.1c )).

Figure 4. Same as figure 1 but for the normalised Nernst velocity $\widetilde {v_N} = Lv_N/r_L v_t$ (see (4.1d )).

For illustration purposes, figures 1–4 show how $G_d$ , $\eta$ , $u_\gamma$ and $v_N$ vary in space for a linear temperature profile,

(4.2) \begin{align} T(z) = T_0 (1 + z/L), \end{align}

and for a Gaussian temperature profile,

(4.3) \begin{align} T(z) = T_0 \exp \left [-(z/L)^2\right ], \end{align}

both with an isobaric density profile $n \propto 1/T$ .

We shall now derive the linear dispersion relation and growth rate for the collisional whistler instability using two approaches – a phase-space-based approach and a more traditional approach presented in the appendix.

4.1. Derivation via Wigner–Moyal phase-space formulation

Here, we shall derive the phase-space analogue of (3.22) that governs the Wigner function of the complex mode amplitude $\psi$ , defined as

(4.4) \begin{align} \mathcal {W}(z, k_z, t) = \int \mathrm {d} s \, \psi ^*\left (z + \frac {s}{2}, t \right ) \psi \left (z - \frac {s}{2}, t \right ) \exp (i k_z s) . \end{align}

This will bring the Hamiltonian structure of (3.22) to light, allowing us then to extract the dispersion relation and growth rate immediately. The following derivation closely follows the presentation of Ruiz et al. (2016), so the reader is invited to consult that reference along with the brief summary provided in Appendix C for more detailed explanations of the steps involved. A more conventional derivation of the same dispersion relation and growth rate, based on a polar decomposition of the wavefield into an amplitude and phase, is presented in Appendix D.

To begin, let us introduce the state vector $|\psi \rangle$ whose spatial projection is given as $\langle z | \psi \rangle = \psi (z)$ . Let us also introduce the operators $\widehat {z}$ and ${\widehat {k}}_z$ whose action on state vectors is given respectively as $\langle z | {\widehat {z}} |\psi \rangle = z \psi (z)$ and $\langle z | {\widehat {k}}_z |\psi \rangle = - i \partial _{z} \psi (z)$ . Then, (3.22) can be viewed as the spatial projection of the Schrödinger equation

(4.5) \begin{align} i \partial _{t} |\psi \rangle = {\widehat {D}} |\psi \rangle, \end{align}

where the non-Hermitian Hamiltonian operator $\widehat {D}$ has the form

(4.6) \begin{align} {\widehat {D}} = {\widehat {k}}_z \left [ G_d({\widehat {z}}) - i \eta ({\widehat {z}}) \right ]{\widehat {k}}_z + {\widehat {k}}_z \left [ v_N({\widehat {z}}) + i u_\gamma ({\widehat {z}}) \right ] . \end{align}

By right-multiplying (4.5) by $\langle \psi |$ and subtracting its adjoint equation, one arrives at the evolution equation for the density operator ${\widehat {W}} \doteq |\psi \rangle \langle \psi |$ :

(4.7) \begin{align} i \partial _{t} {\widehat {W}} = {\widehat {D}} {\widehat {W}} - {\widehat {W}} {\widehat {D}}^\dagger . \end{align}

Finally, applying the Wigner–Weyl transform (WWT, see Appendix C) yields the Wigner–Moyal kinetic equation that governs the Wigner function (4.4) of the fluctuations:

(4.8) \begin{align} \partial _{t} \mathcal {W} = i \mathcal {W} \star \mathcal {D}^* - i \mathcal {D} \star \mathcal {W} \equiv 2 \, \textrm {Im} \left ( \mathcal {D}_H \star \mathcal {W} \right ) + 2 \, {\textrm {Re}} \left (\mathcal {D}_A\star \mathcal {W} \right ), \end{align}

where the Moyal product $\star$ is defined in Appendix C, and $\mathcal {D}$ is the dispersion function whose Hermitian and anti-Hermitian parts are

(4.9a) \begin{align} \mathcal {D}_H &= k_z v_N(z) + k_z^2 G_d(z) + \frac {1}{2} \partial _{z} u_\gamma (z) + \frac {1}{4} \partial _{z}^2 G_d(z), \end{align}

(4.9b) \begin{align} \mathcal {D}_A &= k_z u_\gamma (z) - k_z^2 \eta (z) - \frac {1}{2} \partial _{z}v_N(z) - \frac {1}{4} \partial _{z}^2 \eta (z) . \\[6pt] \nonumber \end{align}

Equation (4.8) states that $\mathcal {D}_H$ governs the Hamiltonian dynamics of the collisional whistler instability in phase space while $\mathcal {D}_A$ acts as the growth rate. This latter point is seen more easily by integrating (4.8) over $k_z$ (i.e. taking the lowest-order ‘fluid’ moment). This gives

(4.10) \begin{align} \partial _{t} \mathcal {I} = 2 \left \langle \mathcal {D}_A \right \rangle \mathcal {I} - \partial _{z} \left [ \left \langle \partial _{k_z} \mathcal {D}_H \right \rangle \mathcal {I} - \frac {1}{2} \partial _{z} (\eta \mathcal {I}) \right ], \end{align}

where moments of $\mathcal {W}$ have been defined as follows:

(4.11) \begin{align} \mathcal {I}(z) \doteq \int \frac {\mathrm {d} k_z}{2\pi } \mathcal {W}(z, k_z) \equiv |\psi (z)|^2, \quad \left \langle f(z) \right \rangle \doteq \frac { \int \mathrm {d} k_z f(z, k_z) \mathcal {W}(z, k_z) }{2\pi \mathcal {I}(z)} . \end{align}

Hence, if the flux at the boundary is negligible, then the total amount of energy contained within the fluctuations remains constant if $\left \langle \mathcal {D}_A \right \rangle = 0$ .

Both the Hermitian and anti-Hermitian parts (4.9) of the instability dynamics contain additional terms (the final two gradient terms) that are absent from previous treatments performed by Bell et al. (2020) based on the short-wavelength approximation. These terms arise because of the spatial variation in the plasma profiles and the consequent non-commutation with the differential operator $\partial _{z}$ , similar to how additional non-Hermitian terms arise when studying zonal flows (Ruiz et al. 2016). More will be said about the gradient terms in § 5.

4.2. Insufficiency of short-wavelength approximation

One might be tempted to drop the gradient drives when the equilibrium length scales are sufficiently long (i.e. to apply the short-wavelength approximation), but this is not always valid. From (4.9b ), we see that the wavevector for the fastest-growing mode at a given point $z$ is given by

(4.12) \begin{align} k_{z, {max}} = \frac {u_\gamma (z)}{2\eta (z)} . \end{align}

Since $u_\gamma \propto \partial _{z} T$ , (4.12) shows that the fluctuation wavelength may be comparable to, or even larger than, the temperature length scale $L_T = (\partial _{z} \log T)^{-1}$ , depending on the prefactor. Indeed, one has

(4.13) \begin{align} k_{z,\textrm {max}} L_T = \frac { \mathcal {M} \Delta _\beta }{4 \alpha _\perp } \beta _0 \sim \mathcal {M}^3 \beta _0, \end{align}

where the plasma beta is defined as

(4.14) \begin{align} \beta _0 = \frac {8\pi n T}{B_z^2}, \end{align}

and the final expression in (4.13) holds in the weakly magnetised limit $\mathcal {M} \ll 1$ . If $k_{z, {max}} L_T \gg 1$ , the additional gradient terms in $\mathcal {D}_A$ can be neglected and one recovers the growth rates of Bell et al. (2020). However, as shown in figure 5, this condition is not satisfied for a weakly magnetised plasma. In this parameter regime, the fastest-growing modes will have wavelengths comparable to the equilibrium scale, so describing them requires the Wigner–Moyal formalism employed here. Also, as we shall show in § 7, it is only by retaining the additional gradient terms in the growth rate that one can obtain non-trivial temperature profiles that are stable to the collisional whistler instability.

Figure 5. Region of parameter space where a geometrical-optics description of the collisional whistler instability is valid (green, $k_z L_T \gg 1$ ), questionable (yellow lined region, $k_z L_T \gtrsim 1$ ) and not valid (red crossed region, $k_z L_T \lt 1$ ). The regions are determined by the expression for $k_z L_T$ given by (4.13) using the transport coefficients of Lopez (2024).

Let us conclude this section with a brief discussion regarding the relevance of our analysis to current laser-plasma experiments. Consider the initial stages of an experiment such as that performed by Meinecke et al. (2022). In such an experiment, pressure balance is quickly established before self-generated magnetic fields have time to grow to appreciable strength (there are no imposed zeroth-order fields). Thus, we can view the initial phase of the experiment as residing within the upper-left corner of figure 5 (low $\mathcal {M}$ and high $\beta _0$ ). As time progresses, dynamo action causes magnetic fields to grow while maintaining constant pressure, so the system evolves to a higher $\mathcal {M}$ state along the trajectory $\beta _0 \sim \mathcal {M}^{-2}$ (if temperature is not constant during this time, then the evolution of $\mathcal {M}$ is even faster, following the shallower trajectory $\beta _0 \sim T^5 \mathcal {M}^{-2}$ ). Along such a trajectory, $k_{z, {max}} L$ will be an increasing function since the contours go as $\beta _0 \sim \mathcal {M}^{-3}$ (4.13). There is a subtlety, however, in that the maximum growth rate for the collisional whistler instability is negative when $k_{z, {max}} L \lesssim 1$ , as will be discussed later (see § 5 and also figure 9). This means that whistler waves will not be excited in the experiment until the plasma is magnetised enough so that $k_{z, {max}} L \sim 1$ , at which point, modes whose wavelengths are comparable with the gradient lengthscale will appear. Due to their early excitation, these modes will have the most time subsequently to manipulate the plasma evolution ¹ (see § 8), but due to their long wavelengths, they can only be accurately described by the Wigner–Moyal analysis performed here.

5. Linear stability condition

Let us consider the case when the whistler-intensity profile has some infinitesimally small (noise-level) initial value that is constant over space. Then, by integrating (D4) over all space, one can readily see that the growth rate for whistlers with a given $k_z = \partial _{z} \theta$ is governed by $\mathcal {D}_A(k_z, z)$ . Hence, the whistlers will be linearly unstable if the maximum growth rate is positive. Since (4.12) implies that

(5.1) \begin{align} \mathcal {D}_A\left [ k_{z, {max}}, z \right ] = \frac {u_\gamma ^2}{4\eta } - \frac {1}{2} \partial _{z}v_N - \frac {1}{4} \partial _{z}^2 \eta, \end{align}

linear instability requires that

(5.2) \begin{align} \frac {u_\gamma ^2}{\eta } \ge 2 \partial _{z}v_N + \partial _{z}^2 \eta . \end{align}

Note that the left-hand side of (5.2) is always positive; therefore, if gradients were neglected in (4.9b ) under the assumption that $k_z L_T \gg 1$ , corresponding to setting the right-hand side of (5.2) to zero, one would erroneously conclude that any non-constant temperature profile would be unstable, i.e. that $\partial _{z} T = 0$ is the only stable profile.

Figure 6. Evolution of a Gaussian pulse advected by an inhomogeneous velocity field $v_N(z) \gt 0$ (i.e. directed towards the right). The pulse spreads when $\partial _{z} v_N(z) \gt 0$ and compresses when $\partial _{z} v_N(z) \lt 0$ .

Instead, we see that two gradient-driven stabilisation mechanisms are present. The first term is the well-known compressional amplification that can result from a perturbation being advected by an inhomogeneous flow. As illustrated in figure 6, if the Nernst advection velocity is a decreasing function of the propagation direction ( $\partial _{z} v_N \lt 0$ ), then the flow can pile up and amplify the initial perturbation, otherwise, when $\partial _{z} v_N \gt 0$ , an initial perturbation will be spread out and stabilised. In the specific context of Nernst advection, this is a well-known mechanism for amplifying magnetic fields near the ablation fronts of laser-compressed fuel pellets (Nishiguchi et al. 1984, 1985). As seen in figure 4, $\partial _{z} v_N \lt 0$ tends to occur when the plasma is weakly magnetised, i.e. at small values of $\mathcal {M}$ .

Figure 7. Diffusion of a sinusoidal perturbation $f(x) = \sin (2x)$ by either a sinusoidal diffusion coefficient $\eta (x) = [2 + \sin (2x)]/10$ (solid) or a constant diffusion coefficient given by the maximum (dotted) or minimum (dashed) value of $\eta (x)$ .

Less well understood is the second stabilisation term due to ‘resistivity curvature’. As shown in figure 7, the diffusion of a perturbation is faster when $\partial _{z}^2 \eta \gt 0$ than homogeneous theory would predict (thus increasing the stability of the system against the perturbation), and the diffusion is slower when $\partial _{z}^2 \eta \lt 0$ (decreasing the stability of the system). Per figure 2, one generally has $\partial _{z}^2 \eta \gt 0$ in the vicinity of a hotspot, with $\partial _{z}^2 \eta$ becoming increasingly larger as the magnetisation level decreases.

In the limit when the geometrical-optics approximation is only weakly violated, this effect can be understood as the result of using the wavelength-averaged resistivity in place of the resistivity when determining the damping of a wave. Indeed, if we define the effective resistivity as

(5.3) \begin{align} \eta _{ {eff}}(z) = \frac {1}{2} \left [ \eta \left (z - \frac {1}{2k} \right ) + \eta \left (z + \frac {1}{2k} \right ) \right ], \end{align}

then in the limit that $k$ is still sufficiently large, a simple Taylor expansion yields

(5.4) \begin{align} \eta _{ {eff}}(z) \approx \eta (z ) + \frac {1}{4 k^2} \partial _{z}^2 \eta (z) . \end{align}

Thus, including the resistivity-curvature correction in (4.9b ) is equivalent to using $k^2 \eta _{ {eff}}$ as the dissipation term in growth rate of Bell et al. (2020).

An alternative explanation for the stabilisation by resistivity curvature can be formulated based on spectral leakage (i.e. the uncertainty principle), as depicted in figure 8. This figure shows the evolution of an initially sinusoidal perturbation as a heuristic diffusion operator is repeatedly applied. This heuristic diffusion operator acts as a low-pass filter for wavevectors larger than the diffusion scale $k_\eta \sim 1/\sqrt {\eta }$ ; accordingly, if a bilevel diffusion coefficient is used where one value of $\eta$ is much larger than the other, then the heuristic diffusion operator acts as a combined spectral and spatial filter with respect to the diffusive scale of the smaller value $k_{\eta _0}$ and the spatial domain of the larger value. Spectral leakage then enables the entire perturbation to decay away eventually, but at different rates when $\eta$ is a local minimum compared with a local maximum. Indeed, as seen by comparing the central peak at $t = \Delta$ and $t = 2\Delta$ in the figure, diffusion is increased when $\eta$ is concave up and diffusion is decreased when $\eta$ is concave down compared with the nominal diffusion one would expect if only the local value of $\eta$ was considered.

Figure 8. Evolution of a sinusoidal pulse subject to a discrete-time diffusion model (with step size $\Delta$ ) in which the diffusion acts as a low-pass filter with respect to wavevectors larger than $k_\eta \sim 1/\sqrt {\eta }$ ; accordingly, the spectral filter becomes a spatial filter where $\eta = \eta _\infty \to \infty$ .

To illustrate the impact of these additional stabilisation mechanisms, let us consider a plasma with a linear temperature profile $T \sim z/L_T$ and in pressure balance, so that $n \sim 1/T$ . Then, it can be shown (see (6.1)) that (5.2) becomes

(5.5) \begin{align} \frac { \beta _0 \Delta _\beta ^2(\mathcal {M}) }{\alpha _\perp (\mathcal {M}) } \ge \frac {15 \alpha _\perp (\mathcal {M}) }{\beta _0 \mathcal {M}^2} - \frac {15 \alpha _\perp '(\mathcal {M}) }{\beta _0 \mathcal {M}} + \frac {25 \alpha ''_\perp (\mathcal {M}) }{\beta _0} - 10 \beta _\wedge '(\mathcal {M}), \end{align}

where as a reminder, $\mathcal {M}$ is defined in (3.17) with $\epsilon = 0$ . Note that (5.5) is actually independent of the temperature gradient. Hence, when the right-hand side is sufficiently positive, there will be no unstable temperature gradients. It is clear that this will happen for weakly magnetised plasmas (small $\mathcal {M}$ ) or for low- $\beta _0$ plasmas due to the divergent denominator; ² this weakly magnetised regime will be discussed further in §§ 7 and 8.

At fixed $\beta _0$ , the same divergent denominator that ensures stability at low $\mathcal {M}$ causes the system to become unstable at high $\mathcal {M}$ . This implies the existence of a critical value $\mathcal {M}_{ {crit}}$ across which the transition from stable to unstable behaviour occurs. Hence, if one were to set up a simulation similar to Komarov et al. (2018), in which a linear temperature gradient is initialised across a plasma of length $L$ with $T(0)$ held fixed and $T(L)$ allowed to vary between simulations, one would see whistler waves beginning to grow once $T(L)$ exceeds a critical value corresponding to $\mathcal {M}_{ {crit}}$ . It might be tempting to conclude that the destabilisation is due to the temperature gradient exceeding a critical value, but subsequent simulations with increased box size $L$ and the same temperature difference should feature whistler waves continuing to be excited despite the temperature gradient being reduced.

6. Dynamical relevance of collisional whistler instability

The collisional whistler instability will grow on time scales determined by the maximum growth rate, denoted $\gamma _{ {whist}}$ . This is given by (5.1), which can be re-written as

(6.1) \begin{align} \gamma _{ {whist}} = \frac {3}{4} \left [ F_{T,1} + F_{T,2} \frac {L_T^2}{C_T} + F_{T,n} \frac { L_T }{ L_n } + F_{n,1} \left ( \frac { L_T }{ L_n } \right )^2 + F_{n,2} \frac {L_T^2 }{ C_n} \right ] \frac {\Omega r_L^2}{L_T^2}, \end{align}

where we have introduced the temperature length scale $L_T$ , the ‘temperature curvature’ $C_T$ , the density length scale $L_n$ and the ‘density curvature’ $C_n$ , as follows:

(6.2) \begin{align} L_T &= \frac {T}{\partial _{z} T}, \quad C_T = \frac {T}{\partial _{z}^2 T}, \quad L_n = \frac {n}{\partial _{z} n}, \quad C_n = \frac {n}{\partial _{z}^2 n}, \end{align}

and we have also introduced the auxiliary functions

(6.3a) \begin{align} F_{T,1} &= \frac {\beta _0 \mathcal {M} \Delta _\beta ^2 }{6 \alpha _\perp } + \beta _0 \, \partial _{\beta _0} F_{T,2} + \frac {3}{2} \mathcal {M} \, \partial _{\mathcal {M}} F_{T,2}, \end{align}

(6.3b) \begin{align} F_{T,2} &= \frac {2 \beta _\wedge }{3 } + \frac {\alpha _\perp - \mathcal {M} \alpha '_\perp }{\beta _0 \mathcal {M}}, \end{align}

(6.3c) \begin{align} F_{T,n} &= \beta _0 \, \partial _{\beta _0} F_{T,2} - \mathcal {M} \, \partial _{\mathcal {M}} F_{T,2} + \frac {3}{2} \mathcal {M} \, \partial _{\mathcal {M}} F_{n,2}, \end{align}

(6.3d) \begin{align} F_{n,1} &= - \mathcal {M} \, \partial _{\mathcal {M}} F_{n,2} - 2 F_{n,2}, \end{align}

(6.3e) \begin{align} F_{n,2} &= \frac {2 \alpha '_\perp }{3 \beta _0 } . \\[6pt] \nonumber \end{align}

Note that here and in what follows, we use $'$ to denote $\partial _{\mathcal {M}}$ for univariate functions of $\mathcal {M}$ . Also note that all the length scales in (6.2) are signed quantities.

Figure 9. Plots of the various drive terms for the collisional whistler instability, as defined in (6.1). Here, ‘WKB’ refers to the only drive term that survives the short-wavelength approximation: see (6.8). Importantly, note that the colour-bar axis can differ by orders of magnitude between plots.

Figure 9 shows the various drive terms as functions of $\mathcal {M}$ and $\beta _0$ . First, one notes that terms corresponding to the density-gradient drives ( $F_{T,n}$ , $F_{n,1}$ and $F_{n,2}$ ) are generally smaller (by approximately two orders of magnitude) than the terms corresponding to temperature-gradient drives ( $F_{T,1}$ and $F_{T,2}$ ). The temperature-gradient-drive terms are larger because they either contain a factor $\beta _0 \mathcal {M}$ that grows unbounded towards the upper right corner of parameter space, or a divergent factor $1/\beta _0 \mathcal {M}$ that grows unbounded towards the lower left corner. The former corresponds to the ‘WKB’ term $\beta _0 \mathcal {M} \Delta _\beta ^2/6 \alpha _\perp$ in (6.1), which is proportional to the maximum growth rate when no additional gradient terms are included in (4.9b ) (see (6.8)), while the latter is associated with the resistivity curvature. Thus, we expect that the instability dynamics will be largely independent of the density profile when a non-uniform temperature profile is present.

For the collisional whistler instability to be dynamically relevant, the whistler waves must grow faster than the time that it takes for the driving temperature inhomogeneity to diffuse away. The diffusion time $\tau _\kappa = (\partial _{t} \log T)^{-1}$ can be calculated from (3.24)–(3.26) with $\epsilon = 0$ :

(6.4) \begin{align} \tau _\kappa = \frac {L_T^2}{ \Omega r_L^2} \frac { 3 }{ \mathcal {M} \kappa _\parallel \left | 5 + 2 S \right | }, \end{align}

where $S = L_T^2/C_T$ describes the shape of the temperature profile. Hence, one has

(6.5) \begin{align} \gamma _{ {whist}} \tau _\kappa = \frac {9}{4} \frac { F_{T,1} + F_{T,2} S + F_{T,n} L_T L_n^{-1} + F_{n,1} \left ( L_T L_n^{-1} \right )^2 + F_{n,2} L_T^2 C_n^{-1} }{ \mathcal {M} \kappa _\parallel \left | 5 + 2 S \right | }, \end{align}

with dynamical relevance requiring $\gamma _{ {whist}} \tau _\kappa \gt 1$ . Note that when $S = - 5/2$ , the diffusion time becomes infinite ( $T \propto z^{2/7}$ yields a spatially constant heat flux), so the whistlers will be dynamically relevant anywhere there is a positive growth rate.

For simplicity, let us restrict attention to when the density profile is either isobaric or constant:

(6.6) \begin{align} L_n^{-1} = \begin{cases} - L_T^{-1} & \textrm {isobaric} \\ 0 & \textrm {constant} \end{cases}, \quad C_n^{-1} = \begin{cases} 2 L_T^{-2} - C_T^{-1} & \textrm {isobaric} \\ 0 & \textrm {constant} \end{cases} . \end{align}

Hence, one has

(6.7) \begin{align} \gamma _{\textrm {whist}} \tau _\kappa = \frac {9}{4 \mathcal {M} \kappa _\parallel \left | 5 + 2 S \right | } \times \begin{cases} F_{T,1} - F_{T,n} + F_{n,1} + 2 F_{n,2} + (F_{T,2} - F_{n,2}) S & \textrm {isobaric} \\ F_{T,1} + F_{T,2} S & \textrm {constant} \end{cases} . \end{align}

Figure 10. (a,b) Parameter space where the collisional whistler instability is dynamically relevant (green) for the specified density profile, as determined by comparing the peak growth rate $\gamma _{ {whist}}$ with the diffusion time $\tau _\kappa$ associated with standard parallel conduction (see (6.5)). The boundary of this region depends on the shape factor $S \doteq T T''/(T')^2$ , and is roughly symmetric about $S = -2.5$ (e.g. the boundaries for $S = -3$ and $S = -2$ are approximately the same). (c) The same, but when the temperature profile is constant and the instability is instead driven by a density inhomogeneity with shape factor S_n defined analogously to $S$ . (d) The same, but for the growth rate provided in (6.8), which is valid in the short-wavelength limit and is independent of the density profile.

The regions of parameter space where whistlers are dynamically relevant for isobaric or constant density profiles are shown in figure 10(a,b).

First of all, there is no visible difference between the results for an isobaric versus a constant density profile. This is because the density-drive terms in (6.1), $F_{n,1}$ and $F_{n,2}$ , are generally smaller than the principal temperature-drive term $F_{T,1}$ (figure 9). Second, not all of the parameter space is susceptible to whistlers even when the initial profile is diffusion-free ( $S = -2.5$ ). This is because the instability actually disappears for sufficiently low $\mathcal {M}$ and all whistler waves are instead strongly damped. This is in stark contrast to the prediction made with the short-wavelength asymptotic growth rate obtained by Bell et al. (2020):

(6.8) \begin{align} \gamma _{ {wkb}} = \frac {\beta _0 \mathcal {M} \Delta _\beta ^2 }{8 \alpha _\perp } \frac {\Omega r_L^2}{L_T^2} . \end{align}

The region of the instability’s dynamical relevance in this limit, which is shown in figure 10(d), is determined by the quantity

(6.9) \begin{align} \gamma _{ {wkb}}\tau _\kappa = \frac {3 \beta _0\Delta _\beta ^2 }{8 \kappa _\parallel \alpha _\perp } \frac { 1 }{ \left | 5 + 2 S \right | } . \end{align}

Since $\gamma _{ {wkb}} \ge 0$ , this approximation does not capture the strong damping that occurs at low magnetisation, instead predicting that all of the parameter space is susceptible to the collisional whistler instability.

Finally, one should note that there actually exists an instability even when the temperature is constant, driven instead by a density gradient. Indeed, setting $\partial _{z} T$ and $\partial _{z}^2 T$ both equal to zero in (6.1) yields

(6.10) \begin{align} \gamma _{ {whist}} = \frac {3}{4} \left ( F_{n,1} + F_{n,2} \frac {L_n^2 }{ C_n} \right ) \frac {\Omega r_L^2}{L_n^2} . \end{align}

Furthermore, since $T$ is constant, there is no diffusion so $\tau _\kappa$ is infinite; the density-gradient-driven collisional whistler instability will be dynamically relevant whenever $\gamma _{ {whist}} \ge 0$ , with $\gamma _{ {whist}}$ given now by (6.10). This region is shown in figure 10(c). Since $F_{n,2} \gt 0$ while $F_{n,1}$ has no definite sign, the region of dynamical relevance increases as $L_n^2/ C_n$ becomes increasingly positive. Eventually, as $L_n^2/ C_n \to \infty$ , the entire parameter space is susceptible to the density-gradient-driven instability.

That said, we should again emphasise that the density-gradient-driven instability has much smaller growth rates than the temperature-gradient-driven instability. Also, by using figure 9 to estimate $F_{T,1} \sim 10^5$ and $F_{n,1} \sim 10^2$ , (6.1) shows that only when $L_n^{-1} \gtrsim 30 L_T^{-1}$ will the density-gradient drives be important for the collisional whistler instability. We shall defer the study of such isothermal plasmas to future work, and instead consider either isobaric or constant-density plasmas to facilitate comparisons with Meinecke et al. (2022) or Bell et al. (2020), respectively. In these plasmas, the density-gradient drives play a negligibly small role.

7. Global marginally stable temperature profiles

It is interesting to consider what plasma profiles are marginally stable over an arbitrarily large spatial domain, since these can potentially correspond to the final states obtained after the collisional whistler instability has saturated quasilinearly. We shall obtain these global marginally stable states by considering the condition

(7.1) \begin{align} \gamma _{ {whist}} = 0, \end{align}

with $\gamma _{ {whist}}$ given by (6.1), as a differential equation governing $T(z)$ for a prescribed $n(z)$ , since $n$ does not evolve in time. Note that non-trivial (i.e. inhomogeneous) marginally stable profiles are only possible when the gradient terms are included in $\mathcal {D}_A$ ; if instead one were to consider marginally stable states with respect to $\gamma _{ {wkb}}$ given by (6.8), the answer would be simply a uniform temperature profile, regardless of the density profile.

We shall first discuss the general case before considering two special cases in detail. The first special case has the plasma density constrained by pressure balance, as occurs in astrophysical and recent experimental contexts (Markevitch & Vikhlinin 2007; Meinecke et al. 2022). The second special case will be the simpler situation in which the density is constant, corresponding to the analysis performed by Bell et al. (2020).

7.1. General theory

Let us consider monotonic profiles such that $\partial _{z} \mathcal {M} \neq 0$ everywhere, $\mathcal {M}$ having been defined in (3.17) but with $\epsilon = 0$ . Then, one can formally parametrise the inverse function $z(\mathcal {M})$ so that all functions can be considered functions of $\mathcal {M}$ . Suppose further that

(7.2) \begin{align} \partial _{\mathcal {M}} n[z(\mathcal {M})] \neq -\frac {n}{\mathcal {M}} . \end{align}

Then, one has $\partial _{\mathcal {M}} T \neq 0$ everywhere, so the composite map $z[\mathcal {M}(T)]$ can be formally constructed and all functions can be parametrised by $T$ instead of $z$ . One then has

(7.3) \begin{align} \partial _{z} n = n' (\partial _{T} \mathcal {M}) \partial _{z} T, \quad \partial _{z}^2 n = n'' ( \partial _{T} \mathcal {M})^2 (\partial _{z} T)^2 + n' (\partial _{T}^2 \mathcal {M}) (\partial _{z} T)^2 + n' (\partial _{T} \mathcal {M}) \partial _{z}^2 T, \end{align}

where $'$ again denotes $\partial _{\mathcal {M}}$ . Then, (7.1) can be recast in the form

(7.4) \begin{align} G_1(\mathcal {M}) (\partial _{z} T)^2 + G_2(\mathcal {M}) \partial _{z}^2 T = 0, \end{align}

where the two new auxiliary functions are defined as follows:

(7.5a) \begin{align} G_1(\mathcal {M}) &= \frac {2\pi n(\mathcal {M}) }{B_z^2} \mathcal {M} \frac { [\Delta _\beta (\mathcal {M})]^2 }{ \alpha _\perp (\mathcal {M}) } + G_2'(\mathcal {M}) \partial _{T} \mathcal {M}, \end{align}

(7.5b) \begin{align} G_2(\mathcal {M}) &= \beta _\wedge (\mathcal {M}) + \frac {B_z^2}{8\pi } \left \{ \alpha _\perp (\mathcal {M}) \left [ 1 + \mathcal {M} \frac {n'(\mathcal {M})}{n(\mathcal {M})} \right ] - \mathcal {M} \alpha _\perp '(\mathcal {M}) \right \} \frac {\partial _{T} \mathcal {M}}{n(\mathcal {M}) \mathcal {M}^2} . \\[6pt] \nonumber \end{align}

Using the chain rule,

(7.6) \begin{align} \partial _{z} T = (\partial _{\mathcal {M}} T) \partial _{z} \mathcal {M}, \quad \partial _{z}^2 T = (\partial _{\mathcal {M}}^2 T) (\partial _{z} \mathcal {M})^2 + (\partial _{\mathcal {M}} T) \partial _{z}^2 \mathcal {M}, \end{align}

along with the standard relations between the derivatives of inverse functions, viz.

(7.7) \begin{align} \partial _{x} y = \frac {1}{\partial _{y} x}, \quad \partial _{x}^2 y = - \frac {\partial _{y}^2 x}{(\partial _{y}x)^3}, \end{align}

we deduce that (7.4) can be written as a differential equation for $\mathcal {M}$ :

(7.8) \begin{align} \partial _{z}^2 \mathcal {M} = g(\mathcal {M}) (\partial _{z} \mathcal {M})^2, \end{align}

where we have defined

(7.9) \begin{align} g(\mathcal {M}) = \frac { G_2(\mathcal {M}) \partial _{T}^2 \mathcal {M} - G_1(\mathcal {M}) \partial _{T} \mathcal {M} }{G_2(\mathcal {M}) (\partial _{T} \mathcal {M})^2} . \end{align}

Clearly, (7.8) can be trivially satisfied by $\mathcal {M}$ = constant. To obtain a non-trivial solution, note that (7.8) possesses affine symmetry with respect to $z$ , i.e. $z \mapsto c_1 + c_2 z$ ; hence, any solution will have the general form $\mathcal {M}(c_1 + c_2 z)$ with $c_1$ and $c_2$ being the two integration constants. The fact that the two integration constants enter in this manner suggests that it will be simpler to solve for the inverse function $z(\mathcal {M})$ , since one expects $\log \partial _{\mathcal {M}} z$ to satisfy an equation of the form $\partial _{\mathcal {M}} \log \partial _{\mathcal {M}} z = f(\mathcal {M})$ . Indeed, by making use again of (7.7), the nonlinear differential equation (7.8) is recast as a linear differential equation

(7.10) \begin{align} z'' = - g(\mathcal {M}) z' . \end{align}

As this is now a first-order differential equation with respect to $z'$ , the solution to (7.10) can be obtained directly via two successive integrations as

(7.11) \begin{align} z(\mathcal {M}) = z(\mathcal {M}_1) + z'(\mathcal {M}_2) \int _{\mathcal {M}_1}^{\mathcal {M}} \mathrm {d} \mu \, \exp \left [ - \int _{\mathcal {M}_2}^\mu \mathrm {d} m \, g(m) \right ], \end{align}

where $\mathcal {M}_1$ and $\mathcal {M}_2$ are arbitrary values at which boundary conditions can be applied. One notes that (7.11) manifestly respects the affine symmetry of the original equation. For (7.11) to be physically relevant, though, it must be the case that $\mathcal {M}(z) \ge 0$ ; a sufficient condition to ensure positivity is that $g \sim A/\mathcal {M}$ with $A \gt 1$ as $\mathcal {M} \to 0$ , as shown in Appendix F. Realistic friction coefficients do indeed have this property, as we shall see in §§ 7.3 and 7.4.

7.2. Magnetisation staircases as a general class of solutions

Let us now consider a high- $\beta _0$ plasma. Although not obvious from (7.11), in this limit, $\mathcal {M}(z)$ , and thus $T(z)$ , naturally forms a staircase structure. To see this more easily, note that $G_1 \sim O(\delta ^{-1})$ and $G_2 \sim O(1)$ with respect to the small parameter $\delta \sim 1/\beta _0$ (see (7.5), or more simply, see that $F_{T,1}$ is significantly larger than any other coefficient in figure 9 when $\beta _0$ is large). Hence, (7.8) has the general abstract form

(7.12) \begin{align} \delta \, y''(z) = G(y) \left [ y'(z) \right ]^2, \end{align}

where $G$ is a nominally $O(1)$ function. It is well established (Bender & Orszag 1978) that the solutions to such equations can exhibit boundary layers when $\delta \to 0$ ; for (7.12), such boundary layers will occur where $G(y) = 0$ .

Away from the boundary layers, the ‘outer’ solution of (7.12) is approximately constant, i.e. $y \approx y_j$ for some $y_j$ . However, the small gradient of $y$ will eventually bring $G(y)$ sufficiently close to zero to trigger a rapid change in $y$ across the boundary layer to reach the next plateau region where $y \approx y_{j+1}$ . A staircase pattern thereby emerges whose steps are dictated by the root structure of $G$ (equivalently, the inflection points of $y$ ), with the widths $W$ of the steps set by $\delta$ as

(7.13) \begin{align} W \propto a^{G'(y_*)/\delta }, \end{align}

where $y_*$ is a root of $G(y) = 0$ and $a \gt 1$ is a constant that depends on boundary conditions. This behaviour is summarised as follows.

Conjecture. A staircase step forms in the solution of ( 7.12 ) for $|\delta | \ll 1$ when $G(y)$ traverses a root $y_*$ where $G(y_*) = 0$ and $G'(y_*)/\delta \lt 1$ . Therefore, a multi-step staircase forms when $G(y)$ has multiple roots.

Figure 11. Solutions of (7.12) when $G(y)$ takes the form shown in each respective inset. All solutions have $|\delta | = 0.01$ . The ‘arbitrary units’ (a.u.) designation on the $x$ -axis emphasises that, due to affine symmetry, there is formally no scale to the $x$ -dependence of the solutions. All solutions satisfy $y(0) = 0$ , with the other boundary condition $y'(0)$ , which simply controls the horizontal scale of the solution, adjusted for each case to fit the pertinent behaviour on the same axis. The functional forms for the plots shown in panels (a) and (b) can be derived analytically, as shown in Appendix G.

The basis for this conjecture is demonstrated in figure 11, which shows solutions of (7.12) for a polynomial $G(y)$ . The derivation of (7.13) and the analytical solution for certain special cases of $G(y)$ are presented in Appendix G. Let us now demonstrate the role that these staircase solutions play in determining the globally stable temperature profiles for isobaric and constant-density plasmas.

7.3. Isobaric density profile with Chapman–Enskog friction

Let us consider a situation when the density profile is set by pressure balance. This means that $n(T)$ is given as

(7.14) \begin{align} n(T) = \frac {\beta _0 B_z^2}{8\pi T} . \end{align}

One therefore has

(7.15) \begin{align} \mathcal {M}(T) = \left (\frac {T}{\tau } \right )^{5/2}, \quad T(\mathcal {M}) = \tau \mathcal {M}^{2/5}, \quad n(\mathcal {M}) = \frac {\beta _0 B_z^2}{8\pi \tau \mathcal {M}^{2/5}}, \end{align}

where we have introduced the magnetisation temperature

(7.16) \begin{align} \tau \doteq \left ( \frac { \sqrt {2} }{6 \sqrt {\pi }}Z e^2 m^{3/2} c^2 \beta _0 \Omega \log \Lambda \right )^{2/5} . \end{align}

Consequently,

(7.17) \begin{align} \partial _{T} \mathcal {M} = \frac {5}{2} \frac { \mathcal {M}^{3/5} }{ \tau }, \quad \partial _{T}^2 \mathcal {M} = \frac {15}{4} \frac { \mathcal {M}^{1/5} }{ \tau ^2 }, \quad n'(\mathcal {M}) = - \frac {2}{5} \frac {n}{\mathcal {M}} . \end{align}

The auxiliary functions (7.5) and (7.9) therefore take the following forms:

(7.18a) \begin{align} G_1 &= \frac {5}{2 \tau \mathcal {M}^{2/5}} \left \{ \beta _0 \mathcal {M} \frac { [\Delta _\beta (\mathcal {M})]^2 }{ 10 \alpha _\perp (\mathcal {M}) } + G_2'(\mathcal {M}) \mathcal {M} \right \}, \end{align}

(7.18b) \begin{align} G_2 &= \beta _\wedge (\mathcal {M}) + \frac { 3 \alpha _\perp (\mathcal {M}) - 5 \mathcal {M} \alpha _\perp '(\mathcal {M}) }{2 \beta _0 \mathcal {M}}, \end{align}

(7.18c) \begin{align} g(\mathcal {M}) &= \frac { 3 G_2(\mathcal {M}) - 2 \mathcal {M}^{2/5} \tau G_1(\mathcal {M}) }{5 G_2(\mathcal {M}) \mathcal {M}} . \\[6pt] \nonumber \end{align}

Figure 12. (a) Contour plot and (b) lineouts at select $\beta _0$ values for $g(\mathcal {M})$ when the friction coefficients are obtained using the Lorentz collision operator. The dashed black contour in panel (a) indicates the root set $g = 0$ across which a staircase step is expected to form when $\beta _0$ is large.

A contour plot of $g$ versus $\mathcal {M}$ and $\beta _0$ is shown in figure 12, along with lineouts along $\mathcal {M}$ for some values of $\beta _0$ . It is clear that $g(\mathcal {M})$ becomes large for large $\beta _0$ , as anticipated. Moreover, $g(\mathcal {M})$ has a single root corresponding to the single root of $F_{T,1}$ (figure 9), satisfying the criterion for a staircase to form.

Using known asymptotics (B1) of the Lorentz friction coefficients, it is straightforward to show that $\mathcal {M} g \to 8/5$ as $\mathcal {M} \to 0$ . Hence, the temperature profile is guaranteed to be positive everywhere (Appendix F). To obtain a simple analytical approximation for the solution to (7.11), it is reasonable to take

(7.19) \begin{align} g(\mathcal {M}) \approx \frac {8}{5 \mathcal {M}} \left ( 1 - \frac {\mathcal {M}}{\mathcal {M}_*} \right ), \end{align}

with $\mathcal {M}_*$ being the single root. Specifically, when $\beta _0 \gg 1$ , $\mathcal {M}_*$ can be approximately calculated using the $\mathcal {M} \ll 1$ limit of the Lorentz friction coefficients to be

(7.20) \begin{align} \lim _{\beta _0 \to \infty } \mathcal {M}_* = 8 \sqrt { \frac {\alpha _\parallel }{105 \beta _0} } \approx \frac {0.4}{\sqrt {\beta _0}} . \end{align}

Note, importantly, that a geometrical-optics description of this parameter regime is not valid because (4.13) predicts that $k_{z, {max}} L \sim \mathcal {M}_* \ll 1$ when $\mathcal {M}_*^2 \beta _0 \sim 1$ and $\mathcal {M}_* \ll 1$ . The continuation of the root line to small $\beta _0$ can be computed using the $\mathcal {M} \to \infty$ limit of the friction coefficients (B2) to give

(7.21) \begin{align} \lim _{\beta _0 \to 0} \mathcal {M}_* = \frac {2 \sqrt {6}}{\beta _\parallel \beta _0} \approx \frac {3.3}{\beta _0} . \end{align}

Since $\beta _0 \ll 1$ , no staircase is expected to form in this parameter regime. For arbitrary $\beta _0$ , a simple interpolation of the two limits can be used to obtain

(7.22) \begin{align} \mathcal {M}_* \approx \frac {3.3}{\beta _0} + \frac {0.4}{\sqrt {\beta _0}} . \end{align}

As shown in Appendix G, the solution to (7.11) can be computed analytically for $g(\mathcal {M})$ given by (7.19):

(7.23) \begin{align} \frac { z(\mathcal {M}) - z(\mathcal {M}_1) }{ z(\mathcal {M}_2) - z(\mathcal {M}_1) } = \frac { \gamma \left ( - \frac {3}{5}, - \frac {8 \mathcal {M}}{5 \mathcal {M}_*} \right ) - \gamma \left ( - \frac {3}{5}, - \frac {8 \mathcal {M}_1}{5 \mathcal {M}_*} \right ) }{ \gamma \left ( - \frac {3}{5}, - \frac {8 \mathcal {M}_2}{5 \mathcal {M}_*} \right ) - \gamma \left ( - \frac {3}{5}, - \frac {8 \mathcal {M}_1}{5 \mathcal {M}_*} \right ) }, \end{align}

Figure 13. Solution (7.11) for the marginally stable magnetisation $\mathcal {M} \propto T^{5/2}$ at various values of $\beta _0$ for Lorentz friction coefficients (Lopez 2024) and isobaric plasma (7.14). The boundary conditions are $z(\mathcal {M}_1) =0$ , $z(\mathcal {M}_2) = 1$ , $\mathcal {M}_1 = 0.1$ and $\mathcal {M}_2 = 1$ (solid) or $\mathcal {M}_2 = 0.4$ (dashed). The top plot uses the analytical approximation presented in (7.23) with $\mathcal {M}_*$ defined in (7.22), while the bottom plot is the numerically computed solution.

where $\gamma (s,z)$ is the lower incomplete Gamma function (Olver et al. 2010). Importantly, $\gamma (s,0)$ is divergent when $s \lt 0$ so $\mathcal {M}(z)$ is positive-definite. Figure 13 shows the solution (7.23) at different values of $\beta _0$ for two different boundary conditions using the approximation for $\mathcal {M}_*$ provided in (7.22). A step-function profile clearly develops as $\beta _0$ increases for both boundary conditions, demonstrating the robustness of the temperature staircase. For comparison, figure 13 also presents numerically computed solutions of (7.11). Overall, the analytical approximation is seen to capture all the salient features of the temperature profile, but underestimates the sharpness of the staircase step because the approximation (7.19) does not reproduce the correct gradient across the root, i.e. $g'(\mathcal {M}_*)$ .

7.4. Constant density profile with Chapman–Enskog friction

Let us now consider the simpler case of constant density:

(7.24) \begin{align} n(\mathcal {M}) = n_0 . \end{align}

One therefore has

(7.25) \begin{align} \mathcal {M}(T) = \left (\frac {T}{\tau _0} \right )^{3/2}, \quad T(\mathcal {M}) = \tau _0 \mathcal {M}^{2/3}, \end{align}

where the magnetisation temperature now takes the form

(7.26) \begin{align} \tau _0 \doteq \left ( \sqrt {\frac {2 m}{\pi }} \frac { Z e^2 \omega _p^2 \log \Lambda }{ 3 \Omega } \right )^{2/3} . \end{align}

Consequently,

(7.27) \begin{align} \partial _{T} \mathcal {M} = \frac {3}{2} \frac { \mathcal {M}^{1/3} }{ \tau _0 }, \quad \partial _{T}^2 \mathcal {M} = \frac {3}{4} \frac { \mathcal {M}^{-1/3} }{ \tau _0^2 }, \end{align}

and the auxiliary functions (7.5) and (7.9) take the following forms:

(7.28a) \begin{align} G_1 &= \frac {1}{2 \mathcal {M}^{2/3} \tau _0} \left \{ \beta _{\textrm {eff}} \mathcal {M}^{5/3} \frac { [\Delta _\beta (\mathcal {M})]^2 }{ 2 \alpha _\perp (\mathcal {M}) } + 3 G_2'(\mathcal {M}) \mathcal {M} \right \}, \end{align}

(7.28b) \begin{align} G_2 &= \beta _\wedge (\mathcal {M}) + \frac {3}{2} \frac { \alpha _\perp (\mathcal {M}) - \mathcal {M} \alpha _\perp '(\mathcal {M}) }{\beta _{\textrm {eff}} \mathcal {M}^{5/3}}, \end{align}

(7.28c) \begin{align} g(\mathcal {M}) &= \frac { G_2(\mathcal {M}) - 2 \mathcal {M}^{2/3} \tau _0 G_1(\mathcal {M}) }{ 3 G_2(\mathcal {M}) \mathcal {M} }, \\[6pt] \nonumber \end{align}

where we have defined the effective plasma beta

(7.29) \begin{align} \beta _{ {eff}} = \frac {8 \pi n_0 \tau _0}{B_z^2} \approx 2.61 \times 10^{6} \left ( Z \log \Lambda \right )^{2/3} \left ( \frac { n_0 }{ 10^{20}\,\textrm {cm}^{-3} } \right )^{5/3} \left ( \frac {B_z}{ 10^3 \,\textrm {G}}\right )^{-8/3} . \end{align}

Clearly, $\beta _{ {eff}}$ can be made large for realistic plasma parameters, so, provided that a root to (7.28c ) exists, a sharp magnetisation staircase is expected to form for constant density profiles as well.

Figure 14. Same as figure 12, but for constant density. The definition of $\beta _{ {eff}}$ is (7.29).

Figure 14 shows a contour plot of $g$ as a function of $\mathcal {M}$ and $\beta _{ {eff}}$ , along with lineouts along $\mathcal {M}$ for some values of $\beta _{ {eff}}$ . Analogously to figure 12, $g(\mathcal {M})$ becomes large for large $\beta _{ {eff}}$ and has a single root line. Therefore, the behaviour of the solution (7.11) will have the same qualitative features as those seen in figure 13, namely, a positive-definite magnetisation profile possessing a single staircase step across the root whose approximate interpolated form is

(7.30) \begin{align} \mathcal {M}_* \approx \frac {0.5}{\beta _{ {eff}}^{3/8}} + \frac {1.9}{\beta _{ {eff}}^{3/5}} \end{align}

(the two terms individually constitute the $\beta _{ {eff}} \to \infty$ and the $\beta _{ {eff}} \to 0$ limits of $\mathcal {M}_*$ , respectively). By comparing figures 14 and 12, however, we see that $g'(\mathcal {M}_*)$ is larger when the plasma density is constant instead of isobaric; hence, the staircase associated with figure 14 will be sharper than either the analytical approximation given by (7.23) or the numerical solution presented in figure 13.

8. Back-reaction on the background temperature profile

Having discussed at length the linear growth rate of the collisional whistler instability, let us now briefly investigate how the instability modifies the background temperature profile. To second order in $\epsilon$ , (3.26) can be written in two equivalent forms:

(8.1a) \begin{align} \frac {3}{2} n \, \partial _{t} T &= \epsilon ^2 B_z^2 \mathcal {S}_\alpha - \epsilon ^2 n \mathcal {V} \partial _{z} T - \partial _{z} q_z, \end{align}

(8.1b) \begin{align} \frac {3}{2} n \, \partial _{t} T &= - \partial _{z} \mathcal {Q}_z - \frac {\epsilon ^2 B_z^2}{8 \pi } \left \langle \mathcal {D}_A \right \rangle \mathcal {I} . \\[6pt] \nonumber \end{align}

The first form (8.1a ) emphasises the advection–diffusion dynamics involved in the temperature evolution, while the second form (8.1b ) emphasises the flow of energy throughout space and the transfer of energy from plasma to waves. Here, $\mathcal {S}_\alpha$ represents the heating source due to the work done by the resistive ( $\alpha$ ) friction force on the perturbed flow, $\mathcal {V}$ the wave-driven advection velocity due to the work done by the thermoelectric ( $\beta$ ) friction force on the perturbed flow, $\mathcal {Q}_z$ the modified heat flux due to the additional wave-driven Poynting flux contribution and $\left \langle \mathcal {D}_A \right \rangle \mathcal {I}$ the energy sink to excite fluctuations; their respective definitions are

(8.2a) \begin{align} \mathcal {S}_\alpha &\doteq \frac {\eta }{8\pi } \left ( \mathcal {I} \left \langle k_z^2 \right \rangle + \frac {1}{4} \partial _{z}^2 \mathcal {I} \right ) + \epsilon ^2 \frac { \Delta _\alpha c^2}{16 \pi \omega _p^2 \tau _{ei}} \mathcal {I}^2 \left \langle k_z \right \rangle ^2 \left ( 1 + \epsilon ^2 \frac {\mathcal {I}}{2} \right )^{-1} \nonumber \\ &= \frac {\eta }{8\pi } \left ( \mathcal {I} \left \langle k_z^2 \right \rangle + \frac {1}{4} \partial _{z}^2 \mathcal {I} \right ) + O(\epsilon ^2), \end{align}

(8.2b) \begin{align} \mathcal {V} &\doteq \frac {\Omega c^2}{2 \omega _p^2} \left ( \frac { \Delta _\beta \left \langle k_z \right \rangle \mathcal {I} }{ 1 + \epsilon ^2 \mathcal {I}/2 } - \frac { \beta _\wedge \, \partial _{z} \mathcal {I} }{ 2 \sqrt { 1 + \epsilon ^2 \mathcal {I}/2 } } \right ) \nonumber \\ &= \frac {\Omega c^2}{4 \omega _p^2} \left ( 2 \Delta _\beta \left \langle k_z \right \rangle \mathcal {I} - \beta _\wedge \, \partial _{z} \mathcal {I} \vphantom {\frac {}{}} \right ) + O(\epsilon ^2), \end{align}

(8.2c) \begin{align} \mathcal {Q}_z &\doteq q_z + \frac {\epsilon ^2 B_z^2}{16 \pi } \left ( v_N \mathcal {I} + \frac {\mathcal {I} \partial _{z} \eta - \eta \, \partial _{z} \mathcal {I} }{2} \right ), \end{align}

(8.2d) \begin{align} q_z &= - n T \left [ \frac {\tau _{ei}}{m} \left ( \kappa _\parallel - \frac {\epsilon ^2}{2} \Delta _\kappa \mathcal {I} \right ) \partial _{z} T + \epsilon ^2 \frac {\Omega c^2}{2 \omega _p^2} \left ( \Delta _\beta \left \langle k_z \right \rangle \mathcal {I} + \frac {\beta _\wedge }{2} \, \partial _{z} \mathcal {I} \right ) \right ] + O(\epsilon ^3) . \\[6pt] \nonumber \end{align}

The appropriate lowest-order expressions for $\eta$ , $v_N$ and $\mathcal {D}_A$ are given in (4.1) and (4.9b ). It is important to note that, when combined with (4.10), (8.1b ) manifestly conserves the total energy of the electron MHD system of equations (the appropriate expressions for energy conservation are presented in Appendix H).

8.1. Frictional cooling

As required by energy conservation, the growth of whistler waves due to friction implies that the friction must be cooling the temperature profile at the same time (at least volumetrically, i.e. neglecting fluxes). This is counterintuitive since friction is often considered a source of heating instead of cooling. Indeed, the frictional work due to resistivity ( $\mathcal {S}_\alpha$ ) is positive definite and therefore always a heat source. However, the advection velocity $\mathcal {V}$ of the temperature profile due to friction is not sign-definite and, depending on the signs of $\partial _{z} \mathcal {I}$ and $\left \langle k_z \right \rangle$ , it can be aligned with $\partial _{z} T$ and therefore be a cooling flow. ³

More quantitatively, let us suppose that the wave profile is given by a quasi-monochromatic (and also quasi-eikonal) field of the form ⁴

(8.3) \begin{align} \psi (z) = \sqrt {\mathcal {I}(0)} \, \exp \left [ \frac {z}{2 L_I} +i \int _0^z k_{ {max}}(\zeta ) \mathrm {d} \zeta \right ], \end{align}

where $L_I$ is the intensity gradient length scale. One can then calculate the instantaneous resistive heating rate (8.2a ) and thermoelectric advection velocity (8.2b ) as

(8.4) \begin{align} \mathcal {S}_\alpha = \left [ \left ( \frac {u_\gamma }{2\eta } \right )^2 + \left ( \frac {1}{2 L_I} \right )^2 \right ] \frac {\eta \mathcal {I}}{8\pi }, \quad \mathcal {V} = \left ( \frac {u_\gamma ^2}{\eta } - \frac {\beta _\wedge }{m \Omega } \frac {\partial _{z}T}{L_I} \vphantom {\frac {}{}} \right ) \frac {m c^2 \Omega ^2}{4 \omega _p^2} \frac {\mathcal {I}}{\partial _{z}T} . \end{align}

Hence, we see that (i) the resistive heating is manifestly positive-definite, as required and (ii) the thermoelectric advection velocity becomes a cooling flow, i.e. $\textrm {sign}(\mathcal {V}) = \textrm {sign}(\partial _{z} T)$ , when wave intensity and temperature gradients oppose each other, viz. when $L_I^{-1} \lt u_\gamma ^2 m \Omega L_T/(\eta \beta _\wedge T)$ . Furthermore, since the total frictional heating can be expressed as

(8.5) \begin{align} \epsilon ^2 B_z^2 \mathcal {S}_\alpha - \epsilon ^2 n \mathcal {V} \partial _{z} T &= \frac {\epsilon ^2 B_z^2}{32 \pi } \left ( \eta L_I^{-2} - 2 v_N L_I^{-1} - \frac {u_\gamma ^2}{\eta } \right ) \mathcal {I}, \end{align}

the total frictional heating will be negative when

(8.6) \begin{align} \frac { v_N - \sqrt {v_N^2 + u_\gamma ^2} }{\eta } \lt L_I^{-1} \lt \frac { v_N + \sqrt {v_N^2 + u_\gamma ^2} }{\eta } . \end{align}

If we take $\partial _{z} T \gt 0$ , then (8.6) can be equivalently written as

(8.7) \begin{align} - \frac {\mathcal {M} \beta _0}{2 \alpha _\perp } \left ( \sqrt {\beta _\wedge ^2 + \Delta _\beta ^2} + \beta _\wedge \right ) \lt L_T L_I^{-1} \lt \frac {\mathcal {M} \beta _0}{2 \alpha _\perp } \left ( \sqrt {\beta _\wedge ^2 + \Delta _\beta ^2} - \beta _\wedge \right ) . \end{align}

Interestingly, the condition (8.6) is satisfied for a constant intensity profile $L_I^{-1} = 0$ because the two endpoints of (8.6) necessarily have opposite signs (but note that the interval is not symmetric about zero since its centre is $v_N \neq 0$ ). This means that friction will cool the temperature profile in the early stages of the instability when whistlers grow from an initially homogeneous noise-level of fluctuations.

8.2. Reduced heat flux

Next, let us consider how the heat flux gets modified by the collisional whistler instability. For the quasi-eikonal field given by (8.3), the heat flux takes the form

(8.8) \begin{align} - \frac {q_z}{q_0} = \kappa _\parallel + \frac {\epsilon ^2}{2} \left ( \frac {\Delta _\beta ^2}{2 \alpha _\perp } + \frac {\beta _\wedge }{\mathcal {M} \beta _0} \frac {L_T}{L_I} - \Delta _\kappa \right ) \mathcal {I}, \quad q_0 = \frac {n T v_t^2 \tau _{ei}}{L_T} . \end{align}

Hence, we see that the net effect of the instability on the heat flux results from the competition of three terms. The first two $O(\epsilon ^2)$ terms are associated with the Ettingshausen effect, ⁵ which is the additional heat flux (beyond the standard enthalpy flux (Epperlein & Haines 1986)) carried by faster moving, less collisional electrons whose directional symmetry is broken with a mean flow. The first of these terms is always positive and therefore always enhances the heat flux; in contrast, the second term can reduce the heat flux when $L_I$ and $L_T$ are oppositely oriented and can even overcome the first term if $L_T L_I^{-1}$ is sufficiently negative (meaning that $\mathcal {I}$ is sufficiently sharply peaked):

(8.9) \begin{align} \frac {L_T}{L_I} \lt - \frac {\mathcal {M} \beta _0 \Delta _\beta ^2}{2 \alpha _\perp \beta _\wedge } \approx - 981 \beta _0 \mathcal {M}^4, \end{align}

where the final approximation is for $\mathcal {M} \ll 1$ . This heat-flux-reduction mechanism can be readily achieved in the high- $\beta _0$ , low- $\mathcal {M}$ regime in which the temperature staircases discussed in § 7 also form, since in this limit, the right-hand side of (8.9) goes to zero as $-25/\beta _0$ (see (7.20)). The third $O(\epsilon ^2)$ term in (8.8), which is always negative, is the reduction of the effective conductivity due to the transverse magnetic field perturbations generated by the instability causing the temperature gradient and the total magnetic field to become misaligned.

8.3. Marginally stable heat flux

Finally, let us suppose that the temperature profile is in a globally marginally stable state such that $\left \langle \mathcal {D}_A \right \rangle = 0$ (§ 7). Dynamically, one expects an arbitrary temperature profile to be driven towards such a state on the instability time scale, which can be faster than the conduction time scale (see figure 10). This is because, as shown in figure 9, the instability growth rate $\mathcal {D}_A$ is an increasing function of temperature (and conversely, the damping rate is a decreasing function of temperature). Energy conservation (8.1b ) then implies that higher temperatures are increasingly cooled by the instability (thereby decreasing $\mathcal {D}_A$ ), while lower temperatures are increasingly heated (thereby decreasing $-\mathcal {D}_A$ ) to create temperature plateaus separated by transition regions where the temperature profile has remained unchanged because $\mathcal {D}_A \approx 0$ initially. The result is a profile that has $\mathcal {D}_A = 0$ everywhere.

Figure 15. Regions (green) where the total heat flux $\mathcal {Q}_z$ (8.10) is reduced due to the presence of whistler waves generated by the collisional whistler instability at global marginal stability. This reduction is controlled by the length scale ratio $L_T L_I^{-1}$ and occurs when $ \Gamma \gt 0$ (8.11); these regions are shown in green above the correspondingly labelled line. For $L_T L_I^{-1} \lt -2$ , the entire plot range has a reduced total heat flux.

When $\left \langle \mathcal {D}_A \right \rangle = 0$ globally, the total frictional heating can be written as an energy flux (see (8.1b )), which combines with the heat flux to yield

(8.10) \begin{align} - \frac {\mathcal {Q}_z}{\kappa _\parallel q_0} = 1 - \epsilon ^2 \Gamma \mathcal {I}, \end{align}

where the flux-reduction factor is

(8.11) \begin{align} \Gamma = \frac {1}{2 \kappa _\parallel } \left [ \Delta _\kappa + \frac {5}{2} \frac {\alpha _\perp ' }{ \mathcal {M} \beta _0^2} - \frac {\beta _\wedge }{\mathcal {M} \beta _0} \left ( \frac {L_T}{L_I} + 1 \right ) - \frac {\alpha _\perp }{ \mathcal {M}^2 \beta _0^2 } \left ( \frac {L_T}{L_I} + \frac {3}{2} \right ) - \frac {\Delta _\beta ^2}{2 \alpha _\perp } \right ] . \end{align}

Here, we have imposed pressure balance for simplicity; the general expression is obtained by replacing $5/2$ with $L_T \partial _{z} \mathcal {M} /\mathcal {M}$ in the second term. As in (8.8), we see that the intensity gradient is capable of reducing the total heat flux. In this case, the two mechanisms available are the Ettingshausen mechanism discussed in § 8.2 and the frictional cooling discussed in § 8.1 (now coupled to the heat flux by imposing global marginal stability). Both are controlled by the length scale ratio $L_T L_I^{-1}$ . As shown in figure 15, making $L_T L_I^{-1}$ negative causes the heat flux to be reduced ( $|\mathcal {Q}_z| \lt \kappa _\parallel q_0$ ) over a large region of parameter space. ⁶

Indeed, the flux-reduction factor $\Gamma$ can be made arbitrarily large by making $L_T L_I^{-1}$ increasingly negative. To make this more quantitative, let us consider the weakly magnetised (small- $\mathcal {M}$ ) limit and approximate all transport coefficients by their lowest-order asymptotic limits (B1). Then one has

(8.12) \begin{align} \Gamma &\approx 124 \mathcal {M}^2 + \frac {1.34 }{ \beta _0^2} - \frac {0.362}{ \beta _0} \left ( \frac {L_T}{L_I} + 1 \right ) - \frac {0.011 }{ \mathcal {M}^2 \beta _0^2 } \left ( \frac {L_T}{L_I} + \frac {3}{2} \right ) - 1200 \mathcal {M}^4 \nonumber \\ &\approx \frac {1}{\beta _0} \left ( 21.7 - 0.423 \frac {L_T}{L_I} \right ) \approx \frac {\Delta _\kappa }{2\kappa _\parallel } - \frac {0.423}{\beta _0} \frac {L_T}{L_I}, \end{align}

where in the second line, (7.20) has been used to evaluate $\Gamma$ near the steepest point in the temperature staircase, taking $\beta _0$ to be large as well. Using this, we can place a bound on the required value of $L_T L_I^{-1}$ to achieve strong heat-flux reduction as follows. By simple rearrangement, (8.10) can be written as

(8.13) \begin{align} \epsilon ^2 \mathcal {I} = \frac {1}{\Gamma } \left ( 1 + \frac {\mathcal {Q}_z}{\kappa _\parallel q_0} \right ) \le 1, \end{align}

where the inequality ensures that the small-amplitude expansion is not grossly violated. If the heat flux is strongly reduced, then $\mathcal {Q}_z \approx 0$ and one would have

(8.14) \begin{align} \Gamma \ge 1, \end{align}

or, equivalently, using (8.12),

(8.15) \begin{align} - \frac {L_T}{L_I} \ge \frac {\beta _0}{0.423} \left ( \frac {\kappa _\parallel + \kappa _\perp }{2\kappa _\parallel } \right ) \gtrsim 2 \beta _0, \end{align}

since $\kappa _\parallel \le \kappa _\parallel + \kappa _\perp \le 2 \kappa _\parallel$ .

Thus, the collisional whistler instability is capable of reducing the electron heat flux in principle, somewhat similar to the collisionless whistler instability (Levinson & Eichler 1992; Pistinner & Eichler 1998; Komarov et al. 2018; Roberg-Clark et al. 2018). The persistence time $\tau$ of the temperature profile is then nominally lengthened by the same factor $\tau \sim \tau _\kappa /|1 - \epsilon ^2 \Gamma \mathcal {I}|$ until it is ultimately set by the persistence time of the intensity profile itself, which in turn is set by advection and refraction (i.e. evolving the $z$ and the $k_z$ dependence of $\mathcal {W}$ , respectively). The advection-limited persistence time is expected to be comparable to $\tau$ (and thus not a limiting factor) since the total Poynting flux that accounts for the intensity advection (the term in square brackets in (H5)) and the modified heat flux $\mathcal {Q}_z$ used to estimate $\tau$ only differ by subdominant terms. The refraction time scale, however, is difficult to estimate given that it inherently involves breaking the quasi-eikonal ansatz (8.3), requiring one to reconsider the full phase-space dynamics of $\mathcal {W}$ . This is beyond the scope of the present work.

Let us conclude this section with a brief word of caution regarding the bound obtained in (8.15). The inequality (8.15) represents a stricter requirement on $L_I$ compared with (8.9) because (8.15) requires second-order effects to become comparable to the lowest-order heat flux, whereas (8.9) results from comparing two second-order terms. Clearly, such an occurrence would also indicate that the perturbative approach underlying (8.2) has potentially become invalid; if $L_T L_I^{-1}$ is too large, even an infinitesimal intensity will grow quickly in space and become relatively large elsewhere. Therefore, our analysis here is merely a simple first attempt to understand what is required of $\mathcal {I}(z)$ to greatly reduce the heat flux assuming quasilinear theory holds for all $L_T L_I^{-1}$ . Future investigations can be conducted to see how nonlinear physics modifies this constraint. Conversely, it is likely that the heat-flux suppression by the collisional whistler instability will be a small effect for high-beta plasmas within the strict validity of our fluid slab model.

9. Conclusion

In this work, we have shown that the electron MHD equations with Braginskii friction in a 1-D slab geometry are unstable with respect to transverse magnetic perturbations. We call this instability the collisional whistler instability, since the dispersion relation contains the usual group-velocity dispersion of whistler waves. We show that for a large region of parameter space, the fastest-growing/least-damped whistler waves do not satisfy the geometrical-optics approximation. This necessitates using the Wigner–Moyal formalism to describe their dynamics, which we derive (§ 4). Extra terms are found in the instability growth rate involving gradients of the background plasma that would not be present had the geometrical-optics approximation been applied. The physical origin of these terms and their impact on the instability threshold are discussed in § 5.

In particular, we show that the extra stabilisation provided by the new terms allows for non-constant temperature profiles to emerge and persist (§ 7). These stable temperature profiles are expected to be established quickly, on the instability time scale, since the quasilinear damping of the instability on the background temperature drives the system to marginal stability ( $\mathcal {D}_A = 0$ globally). In the high-beta limit, these stable temperature profiles generically take the form of a staircase with affine symmetry (shifts and rescalings of the spatial coordinate). For simple density profiles, viz. constant or isobaric profiles, the staircase has a single step that occurs at low temperature where the plasma is effectively unmagnetised. More exotic density profiles can yield multi-step staircases: e.g. choosing a power-law density profile ( $n \propto \mathcal {M}^\sigma$ ) gives a temperature profile with multiple steps but only when the plasma beta is small, and as a consequence, the multi-step staircase is not ‘sharp’.

Finally, we discuss the back-reaction of the collisional whistler instability on the plasma temperature profile (§ 8). The instability is able to modify the temperature profile via frictional heating and Ettingshausen heat flux so that total energy is conserved. Interestingly, there exists a regime in which the instability cools the plasma via friction rather than heats it; this regime necessarily occurs in the initial stages of the instability. The Ettingshausen heat flux is also capable of cancelling a portion of the conductive heat flux when the intensity gradient of the collisional whistler instability is anti-aligned with the temperature gradient. In principle, the collisional whistler instability might be capable of strongly reducing the heat flux through these mechanisms, but for high-beta plasmas, a strong reduction is unlikely to occur in the manner envisioned here as this would require the wave-intensity profile to be essentially delta-shaped. Non-geometrical-optics behaviour, nonlinear effects or even synergistic interplay with kinetic microinstabilities (as these quickly modify fluid transport coefficients away from the standard Braginskii expressions used here) might relax this conclusion, but dedicated simulations are required to investigate this further. In this sense, our work here should be considered an initial investigation in which a number of simplifications have been made to elucidate the basic physics at play and to obtain initial estimates of the various parameter dependencies. More detailed follow-up investigations can now be performed in which these simplifications are sequentially relaxed.