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. Reference Bell, Kingham, Watkins and Matthews2020). 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 Reference Helander and Sigmar2002; Bott et al. Reference Bott, Cowley and Schekochihin2024), which for the Lorentz collision operator, yields the following expressions for the friction force and for the heat flux (Epperlein Reference Epperlein1984):

(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 Reference Epperlein and Haines1986; Helander & Sigmar Reference Helander and Sigmar2002) 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 Reference Zhang, Fu and Qin2020) 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 (Reference Lopez2024). 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.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. (Reference Ruiz, Parker, Shi and Dodin2016), 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. (Reference Bell, Kingham, Watkins and Matthews2020) 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. Reference Ruiz, Parker, Shi and Dodin2016). 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. (Reference Bell, Kingham, Watkins and Matthews2020). 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 (Reference Lopez2024).

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. (Reference Meinecke2022). 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 evolutionFootnote ¹ (see § 8), but due to their long wavelengths, they can only be accurately described by the Wigner–Moyal analysis performed here.

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 Reference Bender and Orszag1978) 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.

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 Reference Lopez2024) 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. Reference Olver, Lozier, Boisvert and Clark2010). 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.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.Footnote ³

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

(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,Footnote ⁵ which is the additional heat flux (beyond the standard enthalpy flux (Epperlein & Haines Reference Epperlein and Haines1986)) 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.Footnote ⁶

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 Reference Levinson and Eichler1992; Pistinner & Eichler Reference Pistinner and Eichler1998; Komarov et al. Reference Komarov, Schekochihin, Churazov and Spitkovsky2018; Roberg-Clark et al. Reference Roberg-Clark, Drake, Reynolds and Swisdak2018). 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.