1. Introduction
Understanding the transport properties of plasmas in various physical regimes has long attracted the attention of plasma physicists. A key achievement of twentieth-century plasma physics was the derivation of transport models in classical, collisional plasmas using kinetic theory. Such plasmas are weakly coupled (
$n_e \lambda _{\textrm {D}}^3 \gg 1$
, where
$n_e$
is the number density of electrons and
$\lambda _{\textrm {D}}$
is the Debye length), non-degenerate (
$n_e \lambda _{\mathrm{th}e}^{3} \gg 1$
, where
$\lambda _{\mathrm{th}e}$
is the electron thermal de Broglie wavelength) and have thermal electron mean free paths
$\lambda _e$
that are much smaller than the macroscopic length scale
$L_{\mathrm{T}} = (-\boldsymbol{\nabla} _\parallel \ln {T_e})^{-1}$
over which the electron temperature
$T_e$
varies. For such plasmas, Spitzer & Härm (Reference Spitzer and Härm1953) developed a model in which the electron heat flux is given (in a hydrogen plasma) by
$q_e = -3.16(n_e v_{\mathrm{th}e} T_e) \lambda _e / L_{\mathrm{T}}$
, where
$v_{\mathrm{th}e}$
is the electron thermal velocity. For sufficiently large values of
$\lambda _e / L_{\mathrm{T}}$
, the heat flux no longer follows the Spitzer–Härm model, but is instead limited to some fraction
$\alpha _0$
of the so-called free-streaming value:
$q_e \leqslant \alpha _0 q_{\mathrm{fs}}$
(Malone, McCrory & Morse Reference Malone, McCrory and Morse1975; Bell, Evans & Nicholas Reference Bell, Evans and Nicholas1981; Li et al. Reference Li, Zhao, Li and Zheng2010), where
$q_{\mathrm{fs}} \sim n_e v_{\mathrm{th}e} T_e$
(an exact definition is given in (2.2)). A magnetic field of sufficient strength – that is, one that causes the electron Larmor radius
$\rho _e$
to become less than the electron mean free path
$\lambda _e$
– introduces anisotropy into the system by limiting transport across magnetic field lines, as has been described by Braginskii (Reference Braginskii1965). In such a system, the heat transport perpendicular to the magnetic field is suppressed by a factor
$(\rho _e/\lambda _e)^{2}$
compared with the Spitzer–Härm model, whilst the parallel heat flux remains unchanged.
For many decades, it has been suspected that these classical estimates of the parallel heat flux fail in weakly collisional and collisionless, high-
$\beta$
plasmas (
$\beta \gg 1$
, where
$\beta$
is the ratio of thermal to magnetic pressure). The observation of sharp, long-lasting temperature fluctuations in the intracluster medium of galaxy clusters is inconsistent with classical expectations (Zakamska & Narayan Reference Zakamska and Narayan2003; Markevitch et al. Reference Markevitch2003; Richard-Laferrière et al. Reference Richard-Laferrière2023). A primary mechanism to generate such long-lasting fluctuations, the Field instability (Field Reference Field1965), cannot produce them in the plasma outside cluster cores if the classical heat flux is present; but in the absence of thermal conduction the instability produces such fluctuations (see Choudhury, Sharma & Quataert Reference Choudhury, Sharma and Quataert2019). There are also significant differences between simulation predictions and experimental results in inertial confinement fusion (ICF) research, as has been outlined in Rinderknecht et al. (Reference Rinderknecht, Amendt, Wilks and Collins2018), and laser-plasma experiments (Meinecke et al. Reference Meinecke2022). Explaining these discrepancies is an outstanding research problem.
One of the most promising explanations is the presence of kinetic instabilities, which can play a significant role in modifying the transport properties of a plasma. Of the many kinetic instabilities characterised so far (see Bott, Cowley & Schekochihin Reference Bott, Cowley and Schekochihin2024), the whistler heat-flux instability (WHFI) has been identified as one of the most important because of its ability to reduce heat transport (Levinson & Eichler Reference Levinson and Eichler1992; Pistinner & Eichler Reference Pistinner and Eichler1998; Gary & Li Reference Gary and Li2000). Analytical, computational and experimental studies into the WHFI have been carried out in recent years in high-
$\beta$
plasmas. These prior studies have primarily focused on moderately high-
$\beta$
plasmas that satisfy
$\beta _{e} \ll L_{\mathrm{T}}/\rho _e$
, a regime particularly relevant to astrophysical plasmas with large separations between characteristic macroscopic and microscopic length scales. In moderately high-
$\beta$
plasmas, the WHFI saturates with
$\delta\kern-1.6pt B^2 / B_0^2 \sim \beta _{e} \rho _e / L_{\mathrm{T}} \ll 1$
, where
$\delta\kern-1.6pt B$
is the amplitude of the fluctuation in the magnetic field due to the whistler wave and
$B_0$
is the magnitude of the background magnetic field, and suppression of the parallel heat flux to approximately
${\sim } \beta _e^{-1}\ll 1$
times the free-streaming heat flux (Roberg-Clark et al. Reference Roberg-Clark, Drake, Reynolds and Swisdak2016, Reference Roberg-Clark, Drake, Reynolds and Swisdak2018; Komarov et al. Reference Komarov, Schekochihin, Churazov and Spitkovsky2018; Drake et al. Reference Drake, Pfrommer, Reynolds, Ruszkowski, Swisdak, Einarsson, Thomas, Hassam and Roberg-Clark2021; Meinecke et al. Reference Meinecke2022; Lopez, Bott & Schekochihin Reference Lopez, Bott and Schekochihin2025; Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025). This implies a weakening of the power-law dependence of the heat flux on the electron temperature, reducing it from
$T_e^{7/2}$
to
$T_e^{1/2}$
in collisional plasmas (
$\lambda _e \ll L_{\textrm {T}}$
) and from
$T_e^{3/2}$
(the free-streaming scaling) to
$T_e^{1/2}$
in collisionless plasmas. Recently, it has been demonstrated that such weakening also opens the possibility of generating temperature fluctuations across the intracluster medium (Choudhury & Bott Reference Choudhury and Bott2026a
).
The fact that the WHFI saturates with small-amplitude fluctuations in the moderately high-
$\beta$
regime implies the existence of an ultra-high-
$\beta$
regime in which
$\beta _{e} \gtrsim L_{\mathrm{T}}/\rho _e$
, and therefore
$\delta\kern-1.6pt B^2 / B_0^2 \gtrsim 1$
. Some research has been conducted into pressure-anisotropy-driven kinetic instabilities in ultra-high-
$\beta$
plasmas, such as the firehose and mirror instabilities (see Melville, Schekochihin & Kunz Reference Melville, Schekochihin and Kunz2016; Bott et al. Reference Bott, Kunz, Quataert, Squire and Arzamasskiy2025). These studies have shown that very different physics becomes relevant in ultra-high-
$\beta$
plasmas with large-amplitude fluctuations, and hence the scaling relationships discovered for the moderately high-
$\beta$
regime cannot be assumed to remain valid. Furthermore, an experiment with laser generated plasma has recently been carried out to investigate heat-flux suppression due to the WHFI in plasmas approaching the ultra-high-
$\beta$
regime (Vincent et al. Reference Vincent2026); however, a theoretical description of the WHFI in this parameter space does not yet exist. The fact that this is an active area of research necessitates an accurate, analytical model for heat transport in WHFI-unstable, ultra-high-
$\beta$
plasmas.
Such a model would be relevant to a range of astrophysical plasmas, including the reionised intergalactic medium of the early universe, due to the very high values of
$\beta _e$
expected to be present there (Barkana & Loeb Reference Barkana and Loeb2001). There has been much interest in seed magnetic fields in the early universe and how they became amplified to the levels seen in and around galaxies today (St-Onge & Kunz Reference St-Onge and Kunz2018; St-Onge et al. Reference St-Onge, Kunz, Squire and Schekochihin2020; Achikanath Chirakkara et al. Reference Achikanath Chirakkara, Federrath, Trivedi and Banerjee2021; Zhou et al. Reference Zhou, Zhdankin, Kunz, Loureiro and Uzdensky2022, Reference Zhou, Zhdankin, Kunz, Loureiro and Uzdensky2024; Hanebring et al. Reference Hanebring, Juno, Hakim, TenBarge and Pusztai2026). These seed fields could be primordial or Biermann-battery generated (Garaldi, Pakmor & Springel Reference Garaldi, Pakmor and Springel2021) and it has been shown that kinetic instabilities could play an intermediary role in amplifying the seed magnetic field before dynamo action amplifies it further (Schoeffler et al. Reference Schoeffler, Loureiro, Fonseca and Silva2016). The structure of the initial field before dynamo amplification has been shown to lead to measurable differences in present-day observations of field configurations (Vazza et al. Reference Vazza, Brüggen, Gheller, Hackstein, Wittor and Hinz2017; Mtchedlidze et al. Reference Mtchedlidze, Domínguez-Fernández, Du, Brandenburg, Kahniashvili, O’Sullivan, Schmidt and Brüggen2022), raising the possibility that such instabilities could have left identifiable markers. These magnetic fields have strong implications for galaxy formation (Brandenburg & Ntormousi Reference Brandenburg and Ntormousi2023), hence understanding how and to what extent the WHFI may have acted on them is important.
The combination of high
$\beta$
and smaller scale separation (
$L_{\mathrm{T}}/\rho _e \sim 10^1$
–
$10^{3}$
) in high-energy-density plasmas means they can also lie in the ultra-high-
$\beta$
regime. This is of particular note to ICF research, as the large temperature gradients and weakly collisional plasmas created in ICF implosions create conditions in which the WHFI can be expected to develop, providing a mechanism for reducing energy loss from the implosion hot-spot and leading to increased hot-spot temperatures and neutron yields (Walsh et al. Reference Walsh2022; Lopez et al. Reference Lopez, Bott and Schekochihin2025). Magnetic fields have been proposed as a means of suppressing cross-field heat transport through magnetic confinement of electrons (Chang et al. Reference Chang, Fiksel, Hohenberger, Knauer, Betti, Marshall, Meyerhofer, Séguin and Petrasso2011; Walsh et al. Reference Walsh2022, Reference Walsh2025b
); however, optimising an implosion for the presence of kinetic instabilities such as the WHFI also provides a mechanism for suppressing parallel transport.
It is clear, therefore, that the WHFI in ultra-high-
$\beta$
plasmas requires deeper study. In this paper, we carry out particle-in-cell (PIC) simulations with OSIRIS (Fonseca et al. Reference Fonseca2002), a massively parallel open-source PIC code, to answer several key questions about the WHFI and WHFI-regulated heat flux in such ultra-high-
$\beta$
plasmas. We present a brief theoretical background to the WHFI in § 2, followed by the derivation from kinetic theory of how the instability and its associated heat flux may behave in the ultra-high-
$\beta$
, large-amplitude regime if heat transport is dominated by diffusive (§ 2.3.1; see also Ryutov et al. (Reference Ryutov, Drake, Kane, Liang, Remington and Wood-Vasey1999)) or advective (§ 2.3.2) processes. Our simulations of the WHFI are described in § 3. We perform collisionless simulations with two spatial dimensions and three velocity dimensions (2D3V) to probe the transition from the
$\beta _e\lt L_{\mathrm{T}}/\rho _e$
regime to the
$\beta _e \gt L_{\mathrm{T}}/\rho _e$
regime in § 3.2. This is achieved by varying both
$\beta _{e}$
and
$\rho _e / L_{\mathrm{T}}$
independently. These simulations characterise heat flux non-locality and test the predictions made in § 2 for how magnetic-field fluctuations and the parallel heat flux saturate, as well as investigating cross-field transport. In § 3.3 we report simulations with one spatial dimension and three velocity dimensions (1D3V), carried out to confirm the 2D3V findings and to test whether parallel whistler waves are capable of regulating heat flux for sufficiently large values of
$\delta\kern-1.6pt B^2 / B_0^2$
. We discuss our results and main conclusions in § 4.
2. Theory of the WHFI
2.1. Linear theory
In high-
$\beta$
, magnetised (
$\lambda _e \gg \rho _e$
) plasmas, whistler waves propagating down a temperature gradient that is coaligned with the macroscopic magnetic field can grow through cyclotron-resonant interactions with electrons if the parallel heat flux exceeds a critical value:
$q_{e\|}\gtrsim q_{\textrm {fs}}/\beta _e$
(Komarov et al. Reference Komarov, Schekochihin, Churazov and Spitkovsky2018; Bott et al. Reference Bott, Cowley and Schekochihin2024). Physically, the instability arises in this scenario because the free energy associated with the temperature gradient becomes large enough to overcome cyclotron damping. In a weakly collisional plasma, this condition on the heat flux is equivalent to the lower bound on the plasma
$\beta$
: that is,
$\beta _e \gtrsim L_{\mathrm{T}}/\lambda _e$
(Drake et al. Reference Drake, Pfrommer, Reynolds, Ruszkowski, Swisdak, Einarsson, Thomas, Hassam and Roberg-Clark2021). If active, the fastest-growing WHFI modes are right-hand circularly polarised and oriented parallel to the temperature gradient (Bott et al. Reference Bott, Cowley and Schekochihin2024). If
$\lambda _e/L_{\mathrm{T}} \gtrsim \beta _e^{-1}$
, these modes have a characteristic parallel wavenumber
$k_{\|} \sim \rho _e^{-1}$
, and grow at a rate
$\gamma _{\textrm {WW}} \sim \varOmega _e/\beta _e$
, where
$\varOmega _e$
is the electron Larmor frequency. They also have a real frequency, with an associated phase velocity
$v_{\textrm {ph}} \sim v_{\mathrm{th}e}/\beta _e$
directed down the temperature gradient. Oblique whistler modes are also destabilised when
$\beta _e \gtrsim L_{\mathrm{T}}/\lambda _e$
, but their linear growth rate is smaller than the parallel modes by an order-unity factor. The linear theory of the WHFI has been discussed in detail elsewhere (see, e.g. Bott et al. Reference Bott, Cowley and Schekochihin2024), so we do not discuss it further here.
2.2. Moderately high-
$\beta$
, small-amplitude whistler regime
After several linear growth times, previous research has shown that, if
$\beta _e \ll$
$L_{\mathrm{T}}/\rho _e$
– a regime we refer to as moderately high-
$\beta$
– the WHFI saturates while whistler waves still have a small amplitude compared with the macroscopic magnetic field:
$\delta\kern-1.6pt B^2 / B_0^2 \sim \beta _e \rho _e/L_{\mathrm{T}} \ll 1$
(Komarov et al. Reference Komarov, Schekochihin, Churazov and Spitkovsky2018; Roberg-Clark et al. Reference Roberg-Clark, Drake, Reynolds and Swisdak2018; Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025). These studies of the WHFI in moderately high-
$\beta$
plasmas have proposed that the primary mechanism by which small-amplitude whistler waves can reduce thermal transport in a plasma is via resonant interactions with heat-carrying electrons, leading to pitch-angle scattering (Pistinner & Eichler Reference Pistinner and Eichler1998). Quasilinear theory predicts that the scattering frequency of resonant electrons due to these interactions with small-amplitude whistlers scales as
$\nu _{\mathrm{w}} \sim \varOmega _e \delta\kern-1.6pt B^2 / B_0^2 \sim \beta _e v_{\mathrm{th}e}/L_{\textrm {T}}$
(Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025). Scattering at this rate suppresses the anisotropy of the electron distribution function until the WHFI becomes marginally stable. Because whistler waves are right-hand polarised, the fastest-growing parallel whistler modes only resonate with, and therefore scatter, electrons moving up the temperature gradient, and as such are inefficient at reducing thermal conduction (Levinson & Eichler Reference Levinson and Eichler1992). Instead, oblique whistler modes are needed to scatter electrons because they have a left-hand polarisation component which can pitch-angle scatter copropagating electrons (Roberg-Clark et al. Reference Roberg-Clark, Drake, Reynolds and Swisdak2016). These oblique modes, in fact, become the fastest-growing due to the marginal anisotropy generated from the scattering of antiparallel electrons (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; Bott et al. Reference Bott, Cowley and Schekochihin2024; Choudhury & Bott Reference Choudhury and Bott2026b
).
In this regime, the parallel heat flux becomes independent of the temperature gradient in the saturated state (Pistinner & Eichler Reference Pistinner and Eichler1998),
\begin{align} q_{e\parallel } \sim \frac {v_{\mathrm{th}e}}{\nu _{\mathrm{w}}L_{\textrm {T}}} q_{\mathrm{fs}} \sim \frac {1}{\beta _e} q_{\mathrm{fs}} , \end{align}
as has been observed in simulations by Komarov et al. (Reference Komarov, Schekochihin, Churazov and Spitkovsky2018), Roberg-Clark et al. (Reference Roberg-Clark, Drake, Reynolds and Swisdak2018) and Yerger et al. (Reference Yerger, Kunz, Bott and Spitkovsky2025). Here,
\begin{align} q_{\mathrm{fs}} = \int _{v_\parallel \gt 0}^\infty \mathrm{d}^3\boldsymbol{v} \ \frac {1}{2} {m_e}v_\parallel \boldsymbol{v}^2 f_{\mathrm{M}e}(\boldsymbol{v}) \end{align}
is the electron free-streaming heat flux,
$m_e$
is the electron mass,
$\boldsymbol{v}$
is particle velocity and
$f_{\mathrm{M}e}$
is a Maxwellian electron distribution function. Two different explanations for why the parallel heat flux scales in this way have been proposed. Roberg-Clark et al. (Reference Roberg-Clark, Drake, Reynolds and Swisdak2018) argued that the whistler waves act as mobile scattering centres that propagate at the whistler phase velocity
$v_{\mathrm{ph}} \sim v_{\mathrm{th}e} / \beta _e$
, pitch-angle scattering electrons in that frame and limiting the advection speed of thermal energy to this value, resulting in an advective model of the heat flux. In contrast, Komarov et al. (Reference Komarov, Schekochihin, Churazov and Spitkovsky2018) argue that the instability saturates quasilinearly to a marginal state due to pitch-angle scattering of electrons in the frame of the background flow, limiting the distribution function anisotropy and resulting in a diffusive heat-flux model. Subsequent work by Yerger et al. (Reference Yerger, Kunz, Bott and Spitkovsky2025) showed that both advective and diffusive contributions to the heat flux were present and contributed roughly equally to the total heat flux.
Following Drake et al. (Reference Drake, Pfrommer, Reynolds, Ruszkowski, Swisdak, Einarsson, Thomas, Hassam and Roberg-Clark2021), an estimate for the perpendicular heat flux can be made using the characteristic transverse scale of whistler waves,
$\rho _e$
, to estimate the cross-field diffusion
$D_\perp \sim \rho _e^2 \nu _{\mathrm{w}}$
. These estimates, along with
$D_\parallel \sim v_{\mathrm{th}e}^2 / \nu _{\mathrm{w}}$
yields a perpendicular heat flux,
\begin{align} q_{e\perp } \sim \left (\frac {\nu _{\mathrm{w}}}{\varOmega _e}\right )^2 q_{e\parallel } \sim \left (\frac {\rho _e}{L_{\mathrm{T}}}\right )^2 \beta _e q_{\mathrm{fs}} . \end{align}
An analogous estimate can be made for the diamagnetic heat flux by assuming that this process has characteristic length and time scales determined by electron Larmor motion. This leads to a diffusion coefficient
$D_\wedge \sim \rho _e^2 \varOmega _e$
and heat flux,
\begin{align} q_{e\wedge } \sim \frac {\rho _e^2 \varOmega _e}{v_{\mathrm{th}e}^2} \nu _{\mathrm{w}} q_{e\parallel }\sim \frac {\nu _{\mathrm{w}}}{\varOmega _e}q_{e\parallel } \sim \frac {\rho _e}{L_{\mathrm{T}}}q_{\mathrm{fs}} . \end{align}
2.3. Ultra-high-
$\beta$
, large-amplitude whistler regime
In ultra-high-
$\beta$
plasmas with
$\beta _e \rho _e / L_{\mathrm{T}} \gtrsim 1$
, the scaling
$\delta\kern-1.6pt B^2 / B_0^2 \sim \beta _e \rho _e/L_{\mathrm{T}}$
discussed in § 2.2 implies that the characteristic amplitude of the magnetic field fluctuations satisfies
$\delta\kern-1.6pt B / B_0 \gtrsim 1$
. Thus, ultra-high-
$\beta$
plasmas are expected to contain large-amplitude whistler fluctuations. In this scenario, quasilinear theory is no longer valid, and so the saturation mechanism of the WHFI might be expected to be distinct in ultra-high-
$\beta$
plasmas. Large-amplitude whistler fluctuations will also interact with electrons differently than their small-amplitude counterparts, potentially causing the heat flux to depart from (2.1). Accordingly, a new theory of both WHFI saturation and the resulting heat flux is required in ultra-high-
$\beta$
plasmas.
To construct such a theory, we begin with the equations governing the collisionless WHFI, starting with the Vlasov equation for the electron distribution function
$f_e$
,
\begin{align} \frac {\partial {f}_e}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }{f}_e - \frac {e}{m_e} \left ( \boldsymbol{E}+ \frac {\boldsymbol{v} \times \boldsymbol{{B}}}{c} \right ) \boldsymbol{\cdot }\frac {\partial {f}_e}{\partial\kern-0.8pt \boldsymbol{v}} = 0 , \end{align}
together with Faraday’s law and the Maxwell–Ampère law from Maxwell’s equations:
\begin{align} \frac {\partial\kern-1.1pt \boldsymbol{B}}{\partial t} & = -c\boldsymbol{\nabla }\times \boldsymbol{E} ,\\[-10pt]\nonumber \end{align}
\begin{align} \boldsymbol{\nabla }\times \boldsymbol{B} & = \frac {4 \pi }{c} \boldsymbol{j} + \frac {1}{c} \frac {\partial\kern-1pt \boldsymbol{E}}{\partial t} . \end{align}
Here,
$t$
is time,
$\boldsymbol{E}$
is the electric field,
$e$
is the elementary charge,
$c$
is the speed of light and
$\boldsymbol{j}$
is the current.
To simplify the analysis, we assume that the ions are much colder than electrons (
$T_i \ll T_e$
), and so do not affect the dynamics of either the temperature gradient or the whistlers.Footnote
1
Motivated by the relationship between the heat flux, which is driven by macroscopic temperature gradients, and unstable whistler waves, we then decompose
$f_e$
into two parts: a ‘mean’ distribution function
$\tilde {f}_e$
that varies on the macroscopic length scale
$L_{\textrm {T}}$
over which the temperature varies, and a perturbation
$\delta\kern-1.6pt f_e$
to the distribution function that describes the whistler fluctuations. To compute
$\tilde {f}_e$
and
$\delta\kern-1.6pt f_e$
from
$f_e$
, we define a spatial averaging operator
$\langle \boldsymbol{\cdot }\rangle$
, taken over scales large compared with the whistler wavelength but small compared with
$L_{\textrm {T}}$
, such that any mean quantity
$\tilde {X}$
is unaffected by the average, while perturbed quantities
$\delta X$
vanish: namely,
$\langle \tilde X \rangle = \tilde X$
, while
$\langle \delta X \rangle = 0$
. Applying this averaging operator to the Vlasov equation yields the following pair of kinetic equations for
$\tilde {f}_e$
and
$\delta\kern-1.6pt f_e$
:
\begin{align} \frac {\partial\kern-1.7pt \tilde {f}_e}{\partial t} & + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {f}_e - \frac {e}{m_e} \frac {\boldsymbol{v} \times \boldsymbol{\tilde {B}}}{c} \boldsymbol{\cdot }\frac {\partial\kern-1.7pt \tilde {f}_e}{\partial\kern-0.8pt \boldsymbol{v}} = \frac {e}{m_e} \left \langle \left ( \delta\kern-1pt \boldsymbol{E} + \frac {\boldsymbol{v} \times \delta\kern-1pt \boldsymbol{B}}{c} \right ) \boldsymbol{\cdot }\frac {\partial }{\partial\kern-0.8pt \boldsymbol{v}} \delta\kern-1.6pt f_e \right \rangle ,\\[-7pt]\nonumber \end{align}
\begin{align} \frac {\partial }{\partial t} \delta\kern-1.6pt f_e & + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\delta\kern-1.6pt f_e - \frac {e}{m_e} \left ( \delta\kern-1pt \boldsymbol{E} + \frac {\boldsymbol{v} \times \delta\kern-1pt \boldsymbol{B}}{c} \right ) \boldsymbol{\cdot }\frac {\partial\kern-1.7pt \tilde {f}_e}{\partial\kern-0.8pt \boldsymbol{v}} - \frac {e}{m_e} \frac {\boldsymbol{v} \times \boldsymbol{\tilde {B}}}{c} \boldsymbol{\cdot }\frac {\partial }{\partial\kern-0.8pt \boldsymbol{v}} \delta\kern-1.6pt f_e\nonumber\\[5pt]& \quad = \frac {e}{m_e} \left (\delta\kern-1pt \boldsymbol{E} + \frac {\boldsymbol{v} \times \delta\kern-1pt \boldsymbol{B}}{c} \right ) \boldsymbol{\cdot }\frac {\partial }{\partial\kern-0.8pt \boldsymbol{v}} \delta\kern-1.6pt f_e - \frac {e}{m_e} \left \langle \left ( \delta\kern-1pt \boldsymbol{E} + \frac {\boldsymbol{v} \times \delta\kern-1pt \boldsymbol{B}}{c} \right ) \boldsymbol{\cdot }\frac {\partial }{\partial\kern-0.8pt \boldsymbol{v}} \delta\kern-1.6pt f_e \right \rangle . \end{align}
For large-amplitude WHFI-unstable modes in a high-
$\beta$
plasma, both the evolution equation (2.7b
) for
$\delta\kern-1.6pt f_e$
and the equations governing the perturbed electromagnetic fields can be simplified. The characteristic frequency
$\omega _{\textrm {ww}}$
of the electron-Larmor-scale whistler waves that are preferentially excited by the instability in high-
$\beta$
plasmas is much smaller than the electron cyclotron frequency:
$\omega _{\textrm {ww}} \sim \varOmega _e/\beta _e$
. A corollary is that these whistler waves are magnetically dominated, with
$|\delta\kern-1pt \boldsymbol{E}| \sim v_{\mathrm{th}e}|\delta {\boldsymbol{B}}|/(c \beta _e)$
, and the whistler phase velocity
$v_{\textrm {ph}} \sim v_{\mathrm{th}e}/\beta _e$
is much smaller than the electron thermal velocity. Due to the cumulative effect of interactions between electrons and the large-amplitude whistler waves, it follows that the mean distribution function is composed of an isotropic distribution
$\tilde {f}_{e}^{(0)}$
and an anisotropic, heat-flux-carrying component
$\tilde {f}_e^{(1)}$
. Extrapolating (2.1) to
$\beta _e \sim L_{\textrm {T}}/\rho _e$
, we find that
$\tilde {f}_e^{(1)} \sim (q_{e\|}/q_{\textrm {fs}})\tilde {f}_e^{(0)} \sim \tilde {f}_{e}^{(0)}/\beta _e$
is much smaller than
$\tilde {f}_{e}^{(0)}$
. Finally, we order
$\delta\kern-1.6pt f_e \sim \tilde {f}_e^{(1)}$
on the assumption that nonlinear wave–wave coupling could, a priori, be important if the whistler waves have a large amplitude. Bringing these considerations together implies the following ordering of parameters:
\begin{align} \frac {\omega _{\textrm {ww}}}{\varOmega _e} \sim \frac {c}{v_{\mathrm{th}e}} \frac {|\delta\kern-1pt \boldsymbol{E}|}{|\tilde {\boldsymbol{B}}|} \sim \frac {v_{\textrm {ph}}}{v_{\mathrm{th}e}} \sim \frac {\tilde {f}_e^{(1)}}{\tilde {f}_e^{(0)}} \sim \frac {\delta\kern-1.6pt f_e}{\tilde {f}_e^{(0)}} \sim \frac {\rho _e}{L_{\textrm {T}}} \sim \frac {1}{\beta _e} \ll \frac {|\delta\kern-1pt \boldsymbol{B}|}{|\tilde {\boldsymbol{B}}|} \sim k \rho _e \sim 1 . \end{align}
Under this ordering, (2.7b ) becomes, at leading order,
\begin{align} \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\delta\kern-1.6pt f_e - \tilde {\varOmega }_e \frac {\partial }{\partial \varphi } \delta\kern-1.6pt f_e - \frac {e}{m_e} \delta\kern-1pt \boldsymbol{E} & \boldsymbol{\cdot }\frac {\partial\kern-1.7pt \tilde {f}_e^{(0)}}{\partial\kern-0.8pt \boldsymbol{v}} - \frac {e}{m_e c}\boldsymbol{v} \times \delta\kern-1pt \boldsymbol{B} \boldsymbol{\cdot }\frac {\partial\kern-1.7pt \tilde {f}_e^{(1)}}{\partial\kern-0.8pt \boldsymbol{v}} \nonumber \\[4pt]& = \frac {e}{m_e c} \left ( \boldsymbol{v} \times \delta\kern-1pt \boldsymbol{B} \boldsymbol{\cdot }\frac {\partial }{\partial\kern-0.8pt \boldsymbol{v}} \delta\kern-1.6pt f_e - \left \langle \boldsymbol{v} \times \delta\kern-1pt \boldsymbol{B} \boldsymbol{\cdot }\frac {\partial }{\partial\kern-0.8pt \boldsymbol{v}} \delta\kern-1.6pt f_e \right \rangle \right ), \end{align}
where
$\tilde {\varOmega }_e \equiv e |\tilde {\boldsymbol{B}}|/m_e c$
is the frequency of Larmor motions around the macroscopic field
$\boldsymbol{\tilde {B}}$
, and
$\varphi$
is the gyroangle with respect to
$\boldsymbol{\tilde {B}}$
. Meanwhile, Faraday’s law for the perturbed electromagnetic field simplifies to
\begin{align} \frac {\partial }{\partial t} \delta\kern-1pt \boldsymbol{B}= -c\boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{E} , \end{align}
while the displacement current in the Maxwell–Ampére law can be neglected due to the whistlers’ low frequency:
\begin{align} \boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} = \frac {4 \pi }{c} \delta\kern-1pt \boldsymbol{j} \simeq -\frac {4 \pi e}{c} \int \mathrm{d}^3\boldsymbol{v} \ \boldsymbol{v} \, \delta\kern-1.6pt f_e . \end{align}
The perturbed electric field can be determined by taking the first moment of (2.9), self-consistently applying the ordering (2.8), and rearranging the resulting expression:
\begin{align} \delta\kern-1pt \boldsymbol{E} & = -\frac {1}{e \tilde {n}_e} \boldsymbol{\nabla }\boldsymbol{\cdot }\delta {\unicode{x1D64B}}_e + \frac {1}{4 \pi e \tilde {n}_e} \left (\boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \tilde {\boldsymbol{B}}\nonumber\\[8pt]& \quad + \frac {1}{4 \pi e \tilde {n}_e} \left [ \left ( \boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \delta {\boldsymbol{B}} -\left \langle \left ( \boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \delta {\boldsymbol{B}} \right \rangle \right ] . \end{align}
Here,
\begin{align} \delta {\unicode{x1D64B}}_e \equiv \int \mathrm{d}^3\boldsymbol{v} \ \boldsymbol{v} \boldsymbol{v} \, \delta\kern-1.6pt f_e \end{align}
is the perturbed electron pressure tensor; we note that we do not need to distinguish between the total particle velocity and the peculiar velocity
$\boldsymbol{w}$
, because under the ordering (2.8), the perturbed electron bulk-velocity (of order the whistler phase speed) is much smaller than the electron thermal velocity. Faraday’s law for the whistler perturbations becomes
\begin{align} \frac {\partial }{\partial t} \delta\kern-1pt \boldsymbol{B} & = \frac {c}{e \tilde {n}_e}\boldsymbol{\nabla }\times \left \{ \boldsymbol{\nabla }\boldsymbol{\cdot }\delta {\unicode{x1D64B}}_e - \frac {1}{4 \pi } \left (\boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \tilde {\boldsymbol{B}}\right.\nonumber\\[5pt]& \quad \left. - \frac {1}{4 \pi } \left [ \left (\boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \delta {\boldsymbol{B}} -\left \langle \left (\boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \delta {\boldsymbol{B}} \right \rangle \right ] \right \}. \end{align}
This is a closed set of equations that completely characterises the WHFI.
The WHFI is, in essence, driven by distribution-function anisotropy. This manifests in the evolution equation for the perturbed magnetic field via the electron pressure tensor
$\delta {\unicode{x1D64B}}_e$
. Equation (2.13), in turn, implies that
$\delta {\unicode{x1D64B}}_e$
is determined directly from the perturbed electron distribution function
$\delta\kern-1.6pt f_e$
. Thus, the way in which distribution-function anisotropy impacts
$\delta\kern-1.6pt f_e$
ultimately determines the behaviour of the instability. Examining the evolution equation (2.9) for
$\delta\kern-1.6pt f_e$
, each term can be interpreted physically as follows: the streaming of particles through the perturbation, the Larmor motion of electrons around the macroscopic magnetic field, the interaction with the perturbed electric field, driven by the effect of distribution-function anisotropy, and nonlinear interactions (the final two terms). The distribution-function anisotropy provides the free energy that drives the instability. For marginality to be satisfied, this driving term needs to be balanced by one of the other terms. In the moderately high-
$\beta$
regime, this balance is provided by cyclotron damping mediated by the perturbed electric field. Under the ordering (2.8), however, nonlinear interactions become comparable to or larger than cyclotron damping if
$\delta\kern-1.6pt f_e \gtrsim f_e^{(1)}$
. In the ultra-high-
$\beta$
regime, we therefore posit that the saturation mechanism for the large-amplitude WHFI instability is nonlinearity.
Under this assumption, we can deduce a scaling relationship between the perturbed magnetic field and the heat flux. To establish this, a relation between
$\delta\kern-1.6pt B$
and
$\delta\kern-1.6pt f_e$
can be found from Faraday’s law (2.14) by balancing the term responsible for the instability drive with the term responsible for nonlinear damping:
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\delta {\unicode{x1D64B}}_e \sim \frac {1}{4 \pi } \left [ \left (\boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \delta {\boldsymbol{B}} -\left \langle \left (\boldsymbol{\nabla }\times \delta\kern-1pt \boldsymbol{B} \right ) \times \delta {\boldsymbol{B}} \right \rangle \right ] . \end{align}
Estimating the magnitude of the electron pressure tensor as
$\delta {\unicode{x1D64B}}_e \sim (\delta\kern-1.6pt f_e/\tilde {f}_e^{(0)}) p_e$
, using (2.13), (2.15) can be rearranged to give
\begin{align} \frac {\delta\kern-1.6pt B^2}{B_0^2} \sim \beta _{e0} \frac {\delta\kern-1.6pt f_e}{\tilde {f}_e^{(0)}} . \end{align}
Here,
$p_e$
is the isotropic electron pressure and the subscript 0 denotes the initial value of
$\beta _e$
. This distinction matters in the large-amplitude-whistler regime because the magnetic field at saturation can differ substantially from its initial value. For example, at saturation the local plasma beta
$\beta _e$
and Larmor radius
$\rho _e$
[defined using the root-mean-square (r.m.s.) magnetic field
$B_{\mathrm{rms}}$
] are not equal to their initial values. Our assumption that nonlinear interactions cause saturation of the WHFI in ultra-high-
$\beta$
plasmas then implies that
$\delta\kern-1.6pt f_e \sim f_e^{(1)} \sim (q_{e\|}/q_{\textrm {fs}})f_e^{(0)}$
. We deduce that
\begin{align} \frac {\delta\kern-1.6pt B^2}{B_0^2} \sim \beta _{e0} \frac {q_{e\|}}{q_{\textrm {fs}}} . \end{align}
To deduce a scaling theory for both the amplitude of the whistler perturbations and the heat flux in ultra-high-
$\beta$
plasmas, we need a model of the heat flux to combine with (2.17). Motivated by the observation that, in the moderately high-
$\beta$
regime, the heat flux has both advective and diffusive components, we consider both possibilities here.
2.3.1. Diffusive heat-flux model
A diffusive heat flux can be constructed following Ryutov et al. (Reference Ryutov, Drake, Kane, Liang, Remington and Wood-Vasey1999). This model assumes that electrons follow a random walk with a step size equal to the length scale on which the magnetic field changes direction. For whistler fluctuations, this distance is
$\rho _e$
, giving a thermal diffusivity of
$\chi \simeq \rho _e v_{\mathrm{th}e}$
. The heat flux is then
\begin{align} q_e \simeq -n_e \rho _e v_{\mathrm{th}e} \boldsymbol{\nabla }T_e \sim \frac {\rho _e}{L_{\mathrm{T}}} q_{\mathrm{fs}}. \end{align}
Combining this with (2.17), a prediction for the perturbed magnetic energy is found:
\begin{align} \frac {\delta\kern-1.6pt B^2}{B_0^2} \sim \frac {\beta _{e0} \rho _e}{L_{\mathrm{T}}} \sim \left ( \frac {\beta _{e0} \rho _{e0}}{L_{\mathrm{T}}} \right )^{2/3} . \end{align}
Using this to evaluate
$\rho _e$
in (2.18) gives
\begin{align} \frac {q_{e}}{q_{\mathrm{fs}}} = \left (\frac {\rho _{e0}}{L_{\mathrm{T}}} \right )^{2/3} \beta _{e0}^{-1/3} . \end{align}
Assuming that the peak wavenumber is
$k \rho _e \sim 1$
as in the small-amplitude case, then it follows that the characteristic length scale of the whistler fluctuations is smaller than the Larmor radius of electrons in the initial magnetic field:
\begin{align} k \rho _{e0} \sim k \rho _{e} \frac {\delta\kern-1.6pt B}{B_0} \sim \left ( \frac {\beta _{e0} \rho _{e0}}{L_{\mathrm{T}}} \right )^{1/3} \gg 1 . \end{align}
To ensure this theory is self-consistent, the advective component of the heat flux must be much smaller than the diffusive component. The whistler phase velocity
$v_{\mathrm{ph}}$
limits the advective heat flux, and if we assume that the whistler dispersion relation is the same as in the small-amplitude regime, then
\begin{align} \frac {v_{\mathrm{ph}}}{v_{\mathrm{th}e}} \sim \frac {1}{\beta _{e0}} k \rho _{e0} \sim \left ( \frac {\rho _{e0}}{L_{\mathrm{T}}} \right )^{1/3} \beta _{e0}^{-2/3} , \end{align}
giving the ratio of the diffusive to the advective heat fluxes:
\begin{align} \frac {q_e / q_{\mathrm{fs}}}{v_{\mathrm{ph}} / v_{\mathrm{th}e}} \sim \left ( \frac {\beta _{e0} \rho _{e0}}{L_{\mathrm{T}}} \right )^{1/3} \gg 1 . \end{align}
Thus, we conclude that this theory is self-consistent; the rate of advection of heat by these large-amplitude whistler fluctuations is indeed small compared with the diffusion rate. It should be noted that the assumption that the whistler dispersion relation remains the same for large-amplitude whistlers may not be true, as nonlinear effects may alter their dispersion relation.
This model predicts several key differences from the small-amplitude-whistler regime. First, the suppression of heat flux should have a much weaker dependence on
$\beta _{e0}$
, since
$q_{e \parallel } / q_{\mathrm{fs}} \sim \beta _{e0}^{-1}$
in the small-amplitude regime and is only
$\sim \beta _{e0}^{-1/3}$
according to the Ryutov model. Instead, the temperature-gradient length scale
$L_{\mathrm{T}}$
now contributes significantly to the heat flux, whereas there was no length-scale dependence in the small-amplitude regime. Second, the characteristic scale of whistler fluctuations in the saturated state is smaller than that found in moderately high-
$\beta$
plasmas. Finally, the heat flux would also become approximately isotropic because the magnetic field is stochastic, rendering the notion of a direction parallel to the global magnetic field largely meaningless on plasma microscales.
2.3.2. Advective heat-flux model
Alternatively, if the advective heat flux is assumed to dominate, then the behaviour differs materially from that of the diffusive model. This advective dominance may arise if the diffusion of electrons across field lines is heavily suppressed due to the large amplitude of the magnetic field fluctuations, and instead the transport of thermal energy is limited by its advection velocity. In this case, the relevant advection velocity is the whistler phase velocity
$v_{\mathrm{ph}}$
, such that
\begin{align} \frac {q_{e \parallel }}{q_{\mathrm{fs}}} \sim \frac {v_{\mathrm{ph}}}{v_{\mathrm{th}e}} , \end{align}
and so
\begin{align} \frac {\delta\kern-1.6pt B^2}{B_0^2} \sim \beta _{e0} \frac {v_{\mathrm{ph}}}{v_{\mathrm{th}e}} . \end{align}
Suppose further that the phase velocity scales with some power
$\alpha$
of
$\beta _{e0}$
, i.e.
$v_{\mathrm{ph}} \sim \beta _{e0}^{\alpha } v_{\mathrm{th}e}$
. It then follows that
\begin{align} \frac {\delta\kern-1.6pt B^2}{B_0^2} \sim \beta _{e0}^{1+\alpha } . \end{align}
In the small-amplitude whistler regime, the phase velocity scales as
$v_{\mathrm{ph}} \sim \beta _{e0}^{-1}$
; if the phase velocity of large-amplitude whistlers were also to scale with
$\beta _{e0}^{-1}$
(i.e.
$\alpha = -1$
), then the WHFI would saturate with
\begin{align} \delta\kern-1.6pt B \sim B_0, \end{align}
and
\begin{align} q_{e \parallel } \sim \beta _{e0}^{-1} q_{\mathrm{fs}} . \end{align}
In this case, the characteristic wavelength of the whistler fluctuations in saturation is then simply of order the Larmor radius of electrons in the macroscopic magnetic field. However, we emphasise that this value of
$\alpha$
is not necessarily realised, because nonlinear effects associated with large-amplitude waves may alter the whistler dispersion relation.
As with the diffusive heat-flux model, for this theory to be self-consistent, the diffusive heat flux must be significantly smaller than the advective component:
\begin{align} \frac {v_{\mathrm{th}e}}{\nu _{\mathrm{eff}} L_{\mathrm{T}}} \ll \frac {v_{\mathrm{ph}}}{v_{\mathrm{th}e}} . \end{align}
For
$\alpha = -1$
, this requires
\begin{align} \frac {\nu _{\mathrm{eff}}}{ \varOmega _{e0}} \gg \frac {\beta _{e0}\rho _{e0}}{L_{\mathrm{T}}} \gtrsim 1 . \end{align}
Thus, the scattering of electrons by large-amplitude whistlers must strongly disrupt Larmor motion caused by the macroscopic magnetic field in some manner.
The advective heat-flux model predicts behaviour that differs substantially from the diffusive heat-flux model and is instead much closer to that found in moderately high-
$\beta$
plasmas. Unlike the diffusive model, the advective model predicts that neither the heat flux nor the saturated magnetic field strength depends on the temperature-gradient length scale. However, the extent to which heat flux is suppressed depends on how the phase velocity scales with
$\beta _{e0}$
for large-amplitude waves; if
$v_{\textrm {ph}}$
is still inversely proportional to
$\beta _{e0}$
, then the advective heat flux retains the scaling (2.1). Another important distinction between the two models is that, in the advective model, cross-field transport continues to follow the Braginskii scalings (Braginskii Reference Braginskii1965), whereas the diffusive model predicts an approximately isotropic heat flux.
The simulations presented in the following sections test these predictions and determine which of the two transport models is realised in the ultra-high-
$\beta$
regime.
3. Simulations
3.1. Numerical set-up
All simulations were initialised with isobaric initial conditions similar to those used in Yerger et al. (Reference Yerger, Kunz, Bott and Spitkovsky2025), with
\begin{align} n(x) = n_{0} \left (\frac {2}{1 + x /L_x} \right ), \qquad T(x) = T_{\mathrm{c}}\left (1+\frac {x}{L_x}\right ) . \end{align}
Here,
$n$
is the number density of electrons,
$T$
is the electron temperature,
$T_c$
is the electron temperature at the cold (left-hand) boundary and
$L_x$
is the length of the simulation box in the direction parallel to the temperature gradient (the
$x$
direction). Throughout this paper, the subscript
$0$
denotes an initial value at the hot (right-hand) boundary. The reference free-streaming heat flux is taken to be its value at the hot boundary of the simulations:
$q_{\mathrm{fs}} = 0.030 n_e m_e c^3$
. The electron temperature and density are equal to the ion temperature and density at all locations.
For all simulations, thermal-bath boundary conditions were used at the left- and right-hand walls, such that any particle leaving the simulation domain at that boundary is immediately reinjected at the same location with a thermal velocity randomly sampled from a half-Maxwellian distribution for the velocity component perpendicular to the wall, and Maxwellian distributions for the other velocity components. This allowed us to maintain boundary temperatures of
$T_{\mathrm{c}} = 10.2\,\mathrm{keV}$
and
$T_{\mathrm{h}} = 20.4\,\mathrm{keV}$
, respectively, such that
$\sqrt {T_{\mathrm{h}} / m_e} =0.2$
. The top and bottom boundary conditions were periodic. A realistic proton-to-electron mass of
$1836 m_e$
was used to simulate a proton–electron plasma and the grid spacing
$\Delta x$
was chosen such that
$\lambda _{\textrm {D}} = 1.5 \Delta x$
to minimise spurious particle heating. All other simulation parameters are listed in table 1. The number of particles per cell
$n_{\mathrm{ppc}}$
was chosen to ensure that numerical noise was negligible, giving
$n_{\mathrm{ppc}} =3600$
in 1D3V and
$n_{\mathrm{ppc}}=900$
in 2D3V. The only exception to this was for the 2D3V simulation with
$\beta _{e0}=400$
, where a smaller
$n_{\mathrm{ppc}}$
was chosen to make the simulation run time feasible.
$L_x = L_{\mathrm{T}0}/2$
and
$L_y = L_{\mathrm{T}0}/10$
(except for simulations with varying
$L_{\mathrm{T}0} / \rho _{e0}$
, for which
$L_y = 10 \rho _{e0}$
). For simulations with the magnetic field inclined by an angle
$\theta$
, both box dimensions were additionally divided by
$\cos {\theta }$
.
Initial parameters for the simulations performed in this study. Here
$\theta$
is the angle between the temperature gradient and the applied magnetic field,
$n_{\mathrm{ppc}}$
is the number of particles per cell,
$t_{\mathrm{end}}$
is the duration of the simulation and
$\omega _{\mathrm{p}0}$
is the plasma frequency at the hot boundary of the simulation.
Table 1. Long description
The table presents initial parameters for simulations, organized into rows and columns. It includes dimensions, angle theta in degrees, various lengths normalized by rho_e0, number of particles per cell (nppc), beta_e0, and simulation durations in terms of Omega_e0 and omega_p0. The table has 28 rows and 9 columns. Column headers are Dimensions, theta (degrees), Lx/rho_e0, Ly/rho_e0, L_T0/rho_e0, nppc, beta_e0, tend*Omega_e0 (x10^3), and tend*omega_p0 (x10^3). Each row lists specific values for these parameters, detailing the setup for different simulation scenarios.
For the electromagnetic-field boundary conditions, we implemented a masking function in OSIRIS that decreases the field strengths near the thermal bath boundaries following Umeda, Omura & Matsumoto (Reference Umeda, Omura and Matsumoto2001), as subsequently adopted by Komarov et al. (Reference Komarov, Schekochihin, Churazov and Spitkovsky2018) and Yerger et al. (Reference Yerger, Kunz, Bott and Spitkovsky2025). The masking function is of the form
\begin{align} f_{\mathrm{Mask}}(x, L_D, r) = \begin{cases} 1 - \left (r\dfrac {x-L_D}{L_D}\right )^2, & x \leqslant L_D,\\[9pt] 1, & L_D \lt x \lt L_x - L_D,\\[3pt] 1 - \left (r\dfrac {x-L_x+L_D}{L_D}\right )^2, & x \geqslant L_x-L_D, \end{cases} \end{align}
where
$r$
controls the strength of the damping and
$L_D$
is the length of the masking region. The values
$L_D = 2\rho _{e0}$
and
$r=0.25$
were chosen to minimise the fraction of whistler waves reflected at the boundary.
The ratio
$L_{\mathrm{T}0} / \rho _{e0}$
was varied by changing
$L_x / \rho _{e0}$
while keeping the temperature at each boundary the same. The initial plasma beta
$\beta _{e0} = 8\pi n_0 T_{\mathrm{h}} / B_0^2$
was varied by applying a uniform background magnetic field across the simulation domain. Decreasing
$B_0$
increases
$\beta _{e0}$
; however, it also increases
$\rho _{e0}$
, hence larger simulation domains are required for higher
$\beta _{e0}$
simulations, making them more computationally expensive. Because the simulations were initialised with an isobaric equilibrium,
$\beta _{e0}$
is constant throughout the simulation domain. Each simulation was run until the WHFI reached saturation, which also takes longer for high
$\beta _{e0}$
simulations. These increased computational demands limit the maximum value of
$\beta _{e0}$
that can be investigated.
The simulations were successfully benchmarked against the TRISTAN-MP simulations of Yerger et al. (Reference Yerger, Kunz, Bott and Spitkovsky2025) by comparing the time evolution of the magnetic-field strength and parallel heat flux for simulations with similar initial parameters.
3.2. The 2D3V results
3.2.1. Varying plasma beta
We first performed 2D3V simulations of the WHFI designed to access the ultra-high-
$\beta$
regime by increasing
$\beta _{e0}$
at fixed
$L_{\mathrm{T}0} / \rho _{e0}$
. Values of
$\beta _{e0} \in \{20, 40, 60, 80, 120, 180, 280\}$
were used, with
$L_{\mathrm{T}0} / \rho _{e0}=100$
and an external magnetic field applied parallel to the temperature gradient in the
$x$
direction. The magnetic field for a sample simulation (
$\beta _{e0}=180$
) is shown in figure 1, demonstrating the growth and saturation of the WHFI in the ultra-high-
$\beta$
regime. The waves initially form in the hotter-half of the simulation before propagating to the colder half. Superimposing trajectories of a sample of electrons on these field snapshots shows that some electrons become trapped in low-field regions of the plasma. This demonstrates that, in addition to resonant scattering, large-amplitude waves possess an alternative mechanism for suppressing the heat flux: variations in magnetic-field strength are large enough to form effective magnetic mirrors. The large perturbations are also capable of forming locally closed or strongly distorted magnetic field lines that trap electrons.
The magnetic fields and
$T_e$
from a 2D3V WHFI simulation with
$\beta _{e0} = 180$
during (a) exponential growth at
$t\varOmega _{e0} = 60$
and (b) nonlinear growth at
$t\varOmega _{e0} = 1069$
. The temperature increases along the
$x$
direction, and the background magnetic field is parallel to the
$x$
axis. Panel (a) shows the formation of waves with wavefronts oblique to the background magnetic field in the hotter half (right-hand side) of the simulation box and how they propagate to the colder half in (b). Black lines show trajectories of test electrons traced around the times corresponding to the background-field images, demonstrating magnetic trapping and reflections when they interact with whistler waves. Purple shaded regions show the domain of the small-
$L_{\mathrm{T}}$
region at both times.
Temporal evolution of the spatially averaged parallel heat flux (upper) and perturbed magnetic energy (lower) for 2D3V simulations with background magnetic fields parallel to the temperature gradient,
$L_{\mathrm{T}0} / \rho _{e0} = 100$
, and a range of
$\beta _{e0}$
values. Both quantities are spatially averaged over the small-
$L_{\mathrm{T}}$
region and the heat flux is normalised to the free-streaming value. A black line is drawn at
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle = 1$
to show when the perturbed magnetic energy exceeds that of the background field. Markers indicate the time at which saturation was determined to have been reached (the start of the saturated interval for each simulation).
Figure 2. Long description
All heat fluxes initially increase, before peaking and then rapidly decreasing. Plots for higher values of beta saturate at lower heat flux values. The perturbed magnetic energy initially increases very quickly, followed by a long period of slow growth and then saturation. The saturated value increases with beta, until beta=120. For larger values of beta, the saturated perturbed magnetic energy is similar.
The region in the simulation in which the magnetic energy is maximised is found to correspond strongly with the location of steepest temperature gradient, consistent with the formation of a transport barrier by the whistler waves and confirming that the fluctuations draw energy from the temperature gradient (see figure 1). A ‘small-
$L_{\mathrm{T}}$
region’ was defined at each time step as follows. We averaged the temperature over the
$y$
direction and smoothed the result, which gives the temperature profile shown in the lower panels of figure 1. We then calculated the local
$\rho _{e0} / L_{\mathrm{T}}$
at each
$x$
coordinate using the gradient of this temperature profile, which yielded a peaked profile once large enough whistlers are present. Taking the full width at half maximum of a Gaussian fit to this
$\rho _{e0} / L_{\mathrm{T}}$
profile then yielded the domain of the small-
$L_{\mathrm{T}}$
region. The small-
$L_{\mathrm{T}}$
region was constrained not to include the masking regions near the boundaries, and was found to propagate down the temperature gradient with the whistler waves. This defines the region of interest for exploring the WHFI, as it is here that the temperature gradient is steep enough to strongly destabilise whistler waves. The temporal evolution of the heat flux (normalised to the free-streaming value, defined for the particle distribution at the hot wall of the simulation) and perturbed magnetic energy are therefore both averaged over this small-
$L_{\mathrm{T}}$
region. The resulting time histories are shown in figure 2.
Because the simulation is initialised without any heat flux, a time interval of order
${\sim }L_x/v_{\mathrm{th}e}$
is required before the heat flux can develop. Once this has happened, the WHFI is destabilised and starts to grow, causing a corresponding decrease in
$\langle q_{e\parallel } \rangle$
, where angle brackets here denote a spatial average over the small-
$L_{\mathrm{T}}$
region. The initial decrease in the heat flux once the whistler waves are destabilised appears very sharp partly because this is when the small-
$L_{\mathrm{T}}$
region becomes well-defined and is no longer the entire spatial domain. Thus, the rapid change reflects both a physical decrease in the heat flux and a rapid change in the spatial region over which the heat flux is averaged. A similar effect affects the preliminary growth of the magnetic fields. Initially the wave growth is exponential; however, after some time, this growth slows down and continues at a greatly decreased rate until saturation is reached. Figure 2 shows that
$\delta\kern-1.6pt B^2 / B_0^2 \gt 1$
is reached, indicating that the WHFI-perturbed magnetic field becomes comparable to, or larger than, the background field, but only by an order-unity factor.
To characterise the phase velocity of the whistler waves, we calculate the dispersion relation for the different values of
$\beta _{e0}$
investigated. This is done by carrying out a two-dimensional Fourier transform of the magnetic field fluctuation
$B_y + \mathrm{i}B_z$
at every point in time. The mode of interest is than extracted (either
$k_y \rho _{e0} = 0$
for parallel modes or
$k_y \rho _{e0} = 1$
for oblique modes) and then a Fourier transform in time is carried out on this mode, to produce the spectrograms in figure 3. We find that, in the saturated state, the dispersion relation in this large-amplitude regime differs from that measured from simulations in the small-amplitude regime (Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025). To demonstrate this, the dispersion relations for representative simulations (
$\beta _{e0} \in \{60, 280\}$
) are plotted in figure 3. These show that the dispersion relation is linear, following
$\omega / \varOmega _{e0} = \hat {v}_{\mathrm{ph}} k_\parallel \rho _{e0}$
, where
$\hat {v}_{\mathrm{ph}} = v_{\mathrm{ph}} / v_{\mathrm{th}e}$
is a fitting parameter which is the phase velocity normalised by the electron thermal velocity. The substantial decrease in slope between
$\beta _{e0}=60$
and
$\beta _{e0}=280$
shows qualitatively that
$\hat {v}_{\mathrm{ph}}$
decreases with increasing
$\beta _{e0}$
for both parallel and oblique modes.
Spectrograms in wavenumber-frequency space of
$B_y + \mathrm{i}B_z$
for (a,b)
$\beta _{e0} = 60$
and (c,d)
$\beta _{e0} = 280$
in the saturated state. The plots are taken from magnetic field data over the entire spatial domain, excluding the boundary masking regions, and over the last
$350 \varOmega _{e0}^{-1}$
of simulation time. Plots for both parallel (
$k_y \rho _{e0}=0$
) and oblique (
$k_y \rho _{e0}=1$
) wavenumbers are shown in (a), (c) and (b), (d), respectively. The linear best fit for the pixels with maximum intensity at each wavenumber is also shown in each panel, following
$\omega / \varOmega _{e0} = \hat {v}_{\mathrm{ph}} k_\parallel \rho _{e0}$
. For parallel modes, the fits are restricted to
$k_\parallel \rho _{e0} \gt 0$
since parallel modes only resonantly interact with electrons moving antiparallel to the background magnetic field, leading to much higher power in this region.
The Fourier spectra of the whistler waves are shown in figure 4 for both parallel and oblique modes for
$\beta _{e0}=60$
and
$\beta _{e0}=280$
, which appear very similar to those observed in the small-amplitude regime. Notably, the spectra continue to demonstrate
$k_\parallel ^{-4}$
scaling, as has been observed by Roberg-Clark et al. (Reference Roberg-Clark, Drake, Reynolds and Swisdak2018), Yerger et al. (Reference Yerger, Kunz, Bott and Spitkovsky2025) and Choudhury & Bott (Reference Choudhury and Bott2026b
). The Fourier spectra integrated over
$k_y$
are also shown. For
$\beta _{e0} = 60$
, the integration over all modes reveals some unexpected phenomena: a steeper gradient to the spectrum, a more pronounced peak at
$k_\parallel \rho _{e0} \sim 1$
, and more power in the high-wavenumber tail. This may indicate cyclotron resonance playing a role near the transition between the small- and large-amplitude regimes.
The peak parallel wavenumber was defined as the wavenumber at which the wave energy is maximal (i.e. the wavenumber at which
$k_\parallel E_B(k_\|)$
peaks). This value is also plotted in figure 4, for all values of
$\beta _{e0}$
simulated. The peak parallel wavenumber remains around
$k_\parallel \rho _{e0} \sim 1$
for all values of
$\beta _{e0}$
. This may seem to contradict the expectation that
$k_\parallel \rho _{e} \sim 1$
; yet, since
$\delta\kern-1.6pt B / B_0$
does not exceed order unity,
$k_\parallel \rho _{e} \sim k_\parallel \rho _{e0}$
and consequently the expectation remains valid. The peak perpendicular wavenumber
$k_{y, \mathrm{peak}} \rho _{e0}$
was calculated in the same way and is also plotted in figure 4, also remaining around
$k_y \rho _{e0} \sim 1$
.
Fourier spectra of
$B_y + \mathrm{i}B_z$
for 2D3V simulations with
$L_{\mathrm{T}0} / \rho _{e0} = 100$
and a parallel macroscopic magnetic field with (a)
$\beta _{e0} = 60$
and (b)
$\beta _{e0} = 280$
. Shown are spectra for parallel (
$k_y \rho _{e0} = 0$
) and oblique (
$k_y \rho _{e0} = 1$
) modes, as well as spectra integrated over all
$k_y$
. The spectra are normalised such that
$\int \mathrm{d}k_\parallel \ E_B (k_\parallel ) = \langle \delta\kern-1.6pt B_y^2 +\delta\kern-1.6pt B_z^2 \rangle / B_0^2$
, the perturbed magnetic energy perpendicular to the applied field for those modes. As in figure 3, the data are taken over the entire spatial domain, excluding the boundary masking regions, and the last
$350 \varOmega _{e0}^{-1}$
of each simulation. The
$k_\parallel ^{-4}$
scaling typical of whistler waves is shown over its relevant region. (c) The peak parallel wavenumber for both oblique and parallel modes for all of the
$\beta _{e0}$
values tested, where the peak wavenumber (the wavenumber at which
$k_\parallel E_B(k_\|)$
peaks) in the last
$350 \varOmega _{e0}^{-1}$
of each simulation and the mean of these values is plotted. Error bars show the standard deviation of this value. The same is plotted in (d) but for the peak perpendicular wavenumber for modes with
$k_\parallel \rho _{e0} \in \{1,1.5\}$
.
Because the large-amplitude whistler waves have an approximately linear dispersion relation, the phase velocity no longer depends on the peak wavenumber, and instead
$v_{\mathrm{ph}} = \hat {v}_{\mathrm{ph}}v_{\mathrm{th}e}$
. For an advective heat flux, the perturbed magnetic energy is expected to follow
$\delta\kern-1.6pt B^2 / B_0^2 \sim \beta _{e0} v_{\mathrm{ph}} /v_{\mathrm{th}e}$
and the heat flux is expected to follow
$q_{e\parallel } / q_{\mathrm{fs}} \sim v_{\mathrm{ph}} /v_{\mathrm{th}e}$
. The scalings
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle$
and
$\langle q_{e\parallel }/q_{\mathrm{fs}} \rangle$
with
$\beta _{e0}$
are plotted in figure 5, where these quantities have been spatially averaged over the small-
$L_{\mathrm{T}}$
region. Theoretical predictions based on the measured values of
$v_{\mathrm{ph}} / v_{\mathrm{th}e}$
at each
$\beta _{e0}$
, together with (2.24) and (2.25), are also shown. Predictions inferred from both parallel (
$v_{\mathrm{ph}\parallel }$
) and oblique (
$v_{\mathrm{ph}\angle }$
) modes are plotted and shown to be very similar.
The theoretical predictions made in § 2.3.2 for the advective heat-flux model, rather than the diffusive heat-flux model outlined in § 2.3.1, are a close match to the simulation results. In particular,
$q_{e\parallel }/q_{\mathrm{fs}}$
agrees with
$v_{\mathrm{ph}} / v_{\mathrm{th}e}$
for both oblique and parallel modes to within an order-unity prefactor. The more specific prediction (2.28), obtained by assuming that the large-amplitude whistler dispersion relation remains close to its small-amplitude form, also appears to hold:
$q_{e\parallel }/ q_{\mathrm{fs}} \sim \beta _{e0}^{-1}$
describes the heat flux results well. Given this scaling of the heat flux and phase velocity, it would be expected that
$\delta\kern-1.6pt B \sim B_0$
, and indeed for
$\beta _{e0}\rho _{e0}/L_{\mathrm{T}0} \gtrsim 1$
, the perturbed magnetic energy plateaus at a value corresponding to
$\delta\kern-1.6pt B \sim B_0$
(
$\delta\kern-1.6pt B_{\textrm {rms}} \approx 3.3 B_0$
). Because of the large computational expense of running such simulations, it was not possible to investigate larger values of
$\beta _{e0}$
with the same initial conditions; however, a plateau at a similar saturation value is also observed in the angled simulations, which investigated larger values of
$\beta _{e0}$
(see § 3.2.3). These results support the idea that heat transport becomes primarily advective when large-amplitude whistler waves are destabilised. We note that the power-law scaling of
$\delta\kern-1.6pt B^2 / B_0^2$
with
$\beta _{e0}$
in the moderately high-
$\beta$
regime is steeper than the scaling found in previous studies (Komarov et al. Reference Komarov, Schekochihin, Churazov and Spitkovsky2018; Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025). We believe that this is because our simulations were not designed to characterise the moderately high-
$\beta$
regime; the range in
$\beta _{e0}$
between the onset of heat-flux suppression by the WHFI and the transition to the ultra-high-
$\beta$
regime is too narrow for the
$\delta\kern-1.6pt B^2 / B_0^2 \propto \beta _e$
scaling to be resolved.
Scalings of (a) the heat flux and (b) the perturbed magnetic energy with
$\beta _{e0}$
in 2D3V simulations with
$L_{\mathrm{T}0} / \rho _{e0}=100$
, where the values on the
$y$
axes have been spatially averaged over the small-
$L_{\mathrm{T}}$
region. In (a), the best-fitting line of the form
$\langle q_{e\parallel }/q_{\mathrm{fs}} \rangle = A\beta _{e0}^{-1}$
, where
$A \approx 4.8$
is a fitting parameter, is shown by the black dotted line in the ultra-high-
$\beta$
region. The expected values calculated from the whistler phase velocities, found using the best fit to the dispersion relation at each
$\beta _{e0}$
, are also shown for both oblique (
$v_{\mathrm{ph}\angle }$
, corresponding to modes with
$k_y \rho _{e0} = 1$
) and parallel (
$v_{\mathrm{ph}\parallel }$
, corresponding to modes with
$k_y \rho _{e0} = 0$
) modes. The predictions in (b) have been multiplied by a best-fitting order-unity prefactor. The lowest
$\beta _{e0}$
simulation (
$\beta _{e0}=20$
) did not strongly destabilise the WHFI, and hence no dispersion relation could be found.
One important aspect of the advective heat-flux model proposed in § 2.3.2 is that electrons travelling from the hot wall to the cold wall are unable to transport thermal energy towards the cold wall faster than
$v_{\mathrm{ph}}$
via diffusion [cf. (2.29)]. Electron trajectories are shown in figure 6, showing that although a small fraction of particles pass through the small-
$L_{\mathrm{T}}$
region, the majority of particles are either reflected back to the hot wall or become trapped in the small-
$L_{\mathrm{T}}$
region and are advected at approximately the phase velocity. As such, the small-
$L_{\mathrm{T}}$
region acts as an effective transport barrier, as was suggested by figure 1. Indeed, of the particles tracked in figure 6, only 12 % made it through the transport barrier, and only 4.9 % made it through unaffected (defined as reaching the cold wall without travelling towards the hot wall at any point in their trajectory). The electrons that were unaffected by the transport barrier were near-thermal electrons with an average speed of
$v / v_{\mathrm{th}e} = 1.2 \pm 0.5$
and, as such, did not carry significant heat flux. This provides strong evidence that thermal energy cannot be transported down the temperature gradient faster than the whistler phase speed in such a system.
Electron trajectories for particles that were on the hot side of the small-
$L_{\mathrm{T}}$
region at
$t\varOmega _{e0}=370$
in the 2D3V simulation with
$\beta _{e0}=120$
and
$L_{\mathrm{T}0} / \rho _{e0} = 100$
. Only a subset of electrons that met these criteria had their trajectories plotted to avoid overcrowding the figure (coloured lines). However, the plotted median and interquartile range (IQR) of electron positions are computed using all of the electrons that met the criteria. The position of the small-
$L_{\mathrm{T}}$
region at each time is shown in blue. The whistler phase velocity
$v_{\mathrm{ph}}$
, calculated from the linear best fit to the dispersion relation for parallel modes from the same simulation, is also shown. The thermal velocity is also plotted for reference.
Figure 6. Long description
Coloured lines show electron trajectories. A cyan shaded region shows the small-LT region, which moves from the hot to cold wall at approximately the phase velocity, shown by a black dashed line. A solid black line and red shaded region shows the median electron position and inter-quartile range. Electron trajectories rarely cross the cyan region. The median electron position remains on the hot side of the cyan region and moves towards the cold wall at approximately the whistler phase velocity.
Temporal evolution of the spatially averaged parallel heat flux (upper) and perturbed magnetic energy (lower) for 2D3V simulations with background magnetic fields parallel to the temperature gradient and
$\beta _{e0}=120$
. Here
$L_{\mathrm{T}0} / \rho _{e0}$
is varied between 44 and 150 in order to increase
$\beta _{e0} \rho _{e0} / L_{\mathrm{T}0}$
past unity without changing
$\beta _{e0}$
. The heat flux is normalised to the free-streaming value. A black line is drawn at
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle = 1$
to indicate when the perturbed magnetic energy exceeds that of the background field.
Figure 7. Long description
4 lines are shown for both the evolution of heat flux and perturbed magnetic energy. All heat fluxes initially increase, before peaking and then rapidly decreasing. All plots saturate at similar values of around qe/qfs=0.04, except the red line for LT/ρe=150 which saturates with a slightly lower heat flux. The perturbed magnetic energy initially increases very quickly, followed by a long period of slow growth and then saturation. The saturated value is the same regardless of temperature gradient, at around 4.
3.2.2. Varying temperature gradient
To investigate whether the saturation state of the WHFI depends solely on
$\beta _{e0}$
in the ultra-high-
$\beta$
regime, we ran simulations with fixed
$\beta _{e0}=120$
and varied
$L_{\mathrm{T}0} / \rho _{e0}$
. Since a diffusive heat-flux component would depend on the temperature-gradient length scale, this provides a test of the extent to which diffusive heat transport contributes to the total heat flux.
The temporal evolution of
$\langle \delta\kern-1.6pt B^2 /B_0^2 \rangle$
and
$\langle q_{e\parallel } / q_{\mathrm{fs}} \rangle$
is plotted in figure 7. This figure shows that, for
$\beta _{e0} \rho _{e0} / L_{\mathrm{T}0} \gtrsim 1$
, the behaviour of each simulation is largely independent of
$L_{\mathrm{T}0} / \rho _{e0}$
. All simulations reach saturation at similar times with similar values of
$\langle \delta\kern-1.6pt B^2 /B_0^2 \rangle$
and
$\langle q_{e\parallel } / q_{\mathrm{fs}} \rangle$
. The only simulation that differs slightly has
$L_{\mathrm{T}0} / \rho _{e0} = 150$
, and hence is not in the ultra-high-
$\beta$
regime, so is not expected to exhibit the same behaviour. Although
$q_{e\parallel }/q_{\mathrm{fs}}$
is independent of
$L_{\mathrm{T}}/{\rho }_{e0}$
in both ultra-high and moderately high-
$\beta$
regimes, there is no a priori reason to expect it to saturate at the same value in both. Therefore, the slightly lower saturated heat flux in the
$L_{\mathrm{T}0} / \rho _{e0} = 150$
simulation is not entirely unexpected. The approximate independence of
$\langle \delta\kern-1.6pt B^2 /B_0^2 \rangle$
and
$\langle q_{e\parallel } / q_{\mathrm{fs}} \rangle$
from the temperature-gradient length scale is indicative of a change in scaling from the moderately high-
$\beta$
regime, where the saturation value of
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle$
depends on
$L_{\mathrm{T}0} / \rho _{e0}$
. Thus, these results suggest that diffusive contributions to heat transport are minimal. We note that the initial peak in the heat flux at
$t \varOmega _{e0} \lt 100$
occurs at different times for the different simulations, because changing the simulation box width changes the time required for the initial heat flux to be established.
We also find that
$L_{\mathrm{T}0} / \rho _{e0}$
does not significantly affect the dispersion relation or peak wavenumber of the whistler waves. This is shown in figure 8, where
$v_{\mathrm{ph}} / v_{\mathrm{th}e}$
and the peak wavenumber are both approximately independent of
$L_{\mathrm{T}0} / \rho _{e0}$
. This again agrees with the expectations of § 2.3.2 and indicates that the phase velocity depends primarily on
$\beta _{e0}$
.
Values of (a)
$v_{\mathrm{ph}} / v_{\mathrm{th}e}$
and (b) the peak wavenumber for simulations with
$\beta _{e0}=120$
and varying
$L_{\mathrm{T}0} / \rho _{e0}$
for both parallel (
$k_y \rho _{e0} = 0$
) and oblique (
$k_y \rho _{e0} = 1$
) modes. The peak wavenumber was found for each time step in the last
$350 \varOmega _{e0}^{-1}$
of each simulation and the mean of these values is plotted. Error bars show the standard deviation of this value.
Magnetic field and
$T_e$
at (a)
$t\varOmega _{e0} = 54$
, during exponential wave growth, and (b)
$t\varOmega _{e0} = 570$
, during nonlinear growth. Images are taken from a 2D3V simulation with an oblique background magnetic field and
$\beta _{e0} = 400$
. Purple lines show the magnetic field lines in the plane of view, while the colour plot shows the out-of-plane component. Black lines show test-electron trajectories sampled around the times of the field snapshots.
3.2.3. Cross-field transport
If the diffusive heat-flux component is indeed small compared with the advective component, perpendicular transport should remain strongly suppressed. To investigate whether the heat flux remains anisotropic for large-amplitude whistler waves, 2D3V simulations were initialised with conditions similar to those in § 3.2.1 but with the magnetic field inclined at
$45^{\circ }$
to the temperature gradient as described in table 1. This produces a diamagnetic heat flux through the Righi–Leduc effect and introduces a component of the temperature gradient perpendicular to the magnetic field, as in Choudhury & Bott (Reference Choudhury and Bott2026b
). To keep
$L_{\mathrm{T}0} / \rho _{e0} = 100$
, the thermal bath boundary conditions were kept at the same temperatures and the dimensions of the box were changed to
$L_x = 35 \rho _{e0}$
and
$L_y = 7.1 \rho _{e0}$
. Example magnetic fields during exponential whistler growth and nonlinear whistler growth are shown in figure 9 for
$\beta _{e0} = 400$
, showing how the whistler waves propagate along the magnetic field lines, at a
$45^{\circ }$
angle to the temperature gradient. The figure shows locally closed or strongly distorted field-line structures and magnetic mirrors, which effectively trap and reflect electrons. The electron trajectories remain closely tied to the magnetic-field lines, suggesting that perpendicular diffusion is weak.
Scalings of (a) the heat flux and (b) the perturbed magnetic energy with
$\beta _{e0}$
in 2D3V simulations with
$L_{\mathrm{T}0} / \rho _{e0}=100$
and angled magnetic fields, where the values on the
$y$
axes have been spatially averaged over the small-
$L_{\mathrm{T}}$
region and temporally averaged over the saturated interval. In (a), a linear best-fit of the form
$\langle q_{e\parallel }/ q_{\mathrm{fs}} \rangle = A\beta _{e0}^{-1}$
, where
$A \approx 4.6$
is a fitting parameter, is shown by the black dotted line in the ultra-high-
$\beta$
region.
Here
$\langle \delta\kern-1.6pt B^2/B_0^2 \rangle$
and
$\langle q_{e\parallel }/q_{\mathrm{fs}} \rangle$
scale with
$\beta _{e0}$
similarly to the case with magnetic fields applied parallel to the temperature gradient, as shown in figure 10. The beginning of the plateau in
$\langle \delta\kern-1.6pt B^2/B_0^2 \rangle$
seen in figure 5 is found to continue, suggesting that the instability saturates with
$\delta\kern-1.6pt B \sim B_0$
, and the best fit to the heat-flux scaling of
$q_{e\parallel }/q_{\mathrm{fs}} \sim 4.6\beta _{e0}^{-1}$
is very similar to that seen in figure 5. The similarities show that the scaling relationships observed in figure 5 and predicted in § 2.3.2 continue far into the ultra-high-
$\beta$
regime. In an advective heat-flux model, the parallel heat-flux and magnetic-energy scalings should not depend strongly on whether the background magnetic field is parallel or oblique to the temperature gradient. The agreement between the two sets of simulations therefore supports this model.
Temporal evolution of the spatially averaged (solid lines) and r.m.s. (dashed lines) (a)
$q_{e\wedge }$
and (b)
$q_{e\perp }$
, normalised to
$q_{\mathrm{fs}}$
. Averages are taken over the small-
$L_{\mathrm{T}}$
region. Panels (c) and (d) show how the cross-field transport, temporally averaged over the saturated interval, varies with
$\beta _{e0}$
, where
$1\sigma$
(standard deviation of the spatially averaged value over time steps in the saturated interval) error bars are plotted.
The temporal evolution of
$q_{e\wedge } / q_{\mathrm{fs}}$
and
$q_{e\perp } / q_{\mathrm{fs}}$
, as well as their r.m.s. values, are plotted in figure 11. Saturation is reached after around
$1000 \varOmega _{e0}^{-1}$
for each of the simulations. The average and r.m.s. cross-field transport values during the saturated interval are calculated and also plotted in figure 11. As the WHFI grows, the mean and r.m.s. diamagnetic heat flux increase in the ultra-high-
$\beta$
simulations. However, as the instability approaches saturation, the magnitude of the mean diamagnetic heat flux decreases again and begins to fluctuate strongly. The magnitude of these fluctuations decreases with increasing
$\beta _{e0}$
for
$\beta _{e0}\gt 180$
, such that the highest-
$\beta _{e0}$
simulations had relatively small mean and r.m.s. diamagnetic fluxes. The mean perpendicular heat flux remains very small throughout the simulation. This suggests that the initial increase in the mean diamagnetic heat flux before saturation is due to the Righi–Leduc effect, which is then suppressed once the whistler waves become sufficiently large. As there remains a non-zero diamagnetic heat flux due to the Righi–Leduc effect, it is possible that at sufficiently high
$\beta _{e0}$
the diamagnetic heat flux could become larger than the parallel heat flux.
The reason the r.m.s. heat flux is consistently so much higher than the mean cross-field heat flux can be seen clearly in figure 9. The electrons remain largely confined to moving along the magnetic-field lines. However, these field lines undergo large deviations from their initial direction, causing the electrons to follow the crests and troughs of the whistler perturbations. This means that at any given time, an electron is likely to be travelling at an oblique angle to the initial external field, but all of these fluctuations cancel out, leading to very small mean perpendicular transport. The fact that the mean perpendicular transport is very small regardless of
$\beta _{e0}$
rules out a diffusive, isotropic heat-flux model, especially given the perpendicular direction has a component down the temperature gradient. This behaviour is consistent with the theory outlined in § 2.3.2: the WHFI strongly affects parallel heat transport while leaving perpendicular transport suppressed.
3.3. The 1D3V results
Since the whistler waves in the ultra-high-
$\beta$
regime are large-amplitude fluctuations, scattering of electrons is no longer expected to depend solely on resonant processes. Therefore, the requirement that oblique waves be present to resonantly scatter heat-flux-carrying electrons no longer applies, and it is plausible that heat flux suppression in 1D3V simulations may be observed. We have therefore run 1D3V simulations similar to those described in § 3.2 to compare the results with the 2D3V case and test whether they also follow the advective heat-flux model from § 2.3.2.
3.3.1. Varying plasma beta
The 1D3V simulations were run with varying plasma beta, similar to those in § 3.2.1. Because 1D3V simulations are computationally simpler, we were able to investigate higher values of
$\beta _{e0}$
, namely
$\beta _{e0} \in \{40, 60, 80, 120, 180, 280, 400, 600, 1000\}$
(as described in table 1). An example of the magnetic field and electron temperature data is shown in figure 12, showing the formation of whistler waves and a steep temperature gradient in the same spatial region, similar to that observed in the 2D3V simulations.
The magnetic fields and
$T_e$
for a 1D3V simulation with
$\beta _{e0} = 280$
and
$L_{\mathrm{T}0} / \rho _{e0} = 100$
when the WHFI is undergoing nonlinear growth at
$t\varOmega _{e0} = 1617$
. The temperature increases along the
$x$
direction, and the background magnetic field is parallel to the
$x$
axis.
Temporal evolution of the spatially averaged parallel heat flux (upper) and perturbed magnetic energy (lower) for 1D3V simulations with background magnetic fields parallel to the temperature gradient and at a range of
$\beta _{e0}$
values. The heat flux is normalised to the free-streaming value. Both quantities are averaged across the small-
$L_{\mathrm{T}}$
region. A black line is drawn at
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle = 1$
to show when the perturbed magnetic energy exceeds that of the background field and markers are drawn on each of the lines to show where saturation was determined to have been reached (the start of the saturated interval). Particularly sharp gradients are typically caused by rapid changes in the width of the small-
$L_{\mathrm{T}}$
region.
Figure 13. Long description
All heat fluxes initially increase, before peaking and then rapidly decreasing. Plots for higher values of beta saturate at lower heat flux values. The perturbed magnetic energy initially increases very quickly, followed by a long period of slow growth and then saturation. The saturated value increases with beta, for all values of beta.
Figure 13 shows that very strong amplification of the perturbed magnetic fields (
$\delta\kern-1.6pt B ^2 / B_0^2 \gg 1$
) is seen in 1D3V, leading to significant suppression of parallel heat flux despite the absence of oblique whistler modes. The ability of parallel modes to affect both copropagating and counterpropagating electrons is confirmed by the trajectories of individual PIC macroparticles shown in figure 14. Electrons travelling both towards the hot wall and towards the cold wall reverse direction on time scales of tens of
$\varOmega _{e0}^{-1}$
, which is impossible in the small-amplitude-whistler regime where only particles travelling towards the hot wall can be resonantly scattered. Therefore, there must be other mechanisms that act to alter electron trajectories in this large-amplitude-whistler regime, suggesting the theory outlined in § 2.3.2 may apply to 1D3V simulations as well.
Electron trajectories for particles that were on the hot side of the small-
$L_{\mathrm{T}}$
region at
$t\varOmega _{e0}=1200$
in the 1D3V simulation with
$\beta _{e0}=120$
and
$L_{\mathrm{T}0} / \rho _{e0} = 100$
. Only a subset of electrons satisfying these criteria are plotted to avoid overcrowding the figure (coloured lines); however, the plotted median and IQR of electron positions are computed using all electrons that met the criteria. The median and IQR are highly variable due to the smaller number of tracked particles in 1D3V compared with 2D3V. The position of the small-
$L_{\mathrm{T}}$
region at each time is shown in blue. The thermal velocity and whistler phase velocity
$v_{\mathrm{ph}}$
are shown for reference.
Figure 14. Long description
Coloured lines show electron trajectories. A cyan shaded region shows the small-LT region, which moves from the hot to cold wall at approximately the phase velocity, shown by a black dashed line. A solid black line and red shaded region shows the median electron position and inter-quartile range. Electron trajectories rarely cross the cyan region.
Figure 14 confirms that the small-
$L_{\mathrm{T}}$
region acts as a transport barrier in 1D3V, as in 2D3V. Very few electrons pass through the small-
$L_{\mathrm{T}}$
region at all, with the majority being reflected back towards the hot wall. The electrons that make it through the transport barrier typically interact strongly with it, travelling much slower than the thermal velocity when inside the barrier. The fact that some electrons are able to pass through the barrier suggests that some diffusion occurs; however, the fact that this occurs for only a small fraction of the electrons is consistent with the assumption made in § 2.3.2 that diffusion should be negligible.
(a) Dispersion relation for a 1D3V simulation with
$\beta _{e0}=280$
and
$L_{\mathrm{T}0} / \rho _{e0}=100$
, computed using waves during the last
$2000 \varOmega _{e0}^{-1}$
of simulation time during the saturated interval and excluding data in the masking region. Two best-fitting curves have been plotted: a linear dashed line and a quadratic solid line. (b) A histogram showing the probability density of different phase velocities for the same simulation, weighted by the power at each phase velocity – i.e. the brightness in (a).
To determine whether the heat flux in 1D3V simulations is advective, we calculate the whistler phase velocity. The dispersion relation of whistler waves in the 1D3V simulations differs from the 2D3V simulations, leading to potential differences in heat flux between the two cases. This can be seen in figure 15, where the dispersion relation for the simulation with
$\beta _{e0}=280$
has been plotted. The dispersion of the waves has a quadratic dependence on
$k_{\parallel } \rho _{e0}$
when
$k_\parallel \rho _{e0} \lt k_{\parallel \mathrm{peak}} \rho _{e0}$
, and a linear dependence when
$k_\parallel \rho _{e0} \gt k_{\parallel \mathrm{peak}} \rho _{e0}$
, where the peak wavenumber
$k_{\parallel \mathrm{peak}} \rho _{e0}$
is easily identifiable by the region of maximum power. Thus, the dispersion relation is similar to the small-amplitude regime at small wavenumbers, and to the 2D3V large-amplitude regime at large wavenumbers.
The impact of the two different dispersion regimes on the phase velocity is shown in figure 15(b), where the power-weighted probability density of different phase velocities is plotted. It can be seen that the phase velocities are sharply peaked around one value, and that the presence of two different dispersion relations does not result in two different phase velocities being present. This is confirmed by noting that the region of highest power in figure 15(a) is near the intersection of the linear and quadratic fits. A Gaussian was fitted to the probability-density distribution of phase velocities, the mean of which is taken to be the phase velocity for whistler waves in this simulation. The same method was used for all 1D3V simulations to find their phase velocities.
Because in one dimension the perturbed magnetic field continues to grow beyond order-unity amplitudes, it is also possible to test the assumption that the peak wavenumber follows
$k \rho _{e} \sim k \rho _{e0} B_0 / \delta\kern-1.6pt B \sim 1$
. The Fourier spectra of whistler waves in 1D3V simulations, averaged over the final
$2000 \varOmega _{e0}^{-1}$
of the simulation, are plotted in figure 16(a). The spectra remain similar to those in both 2D3V and small-amplitude simulations, following
$k_\parallel ^{-4}$
for
$k_\parallel \rho _{e0} \lesssim 10$
, but the peak wavenumber increases with
$\beta _{e0}$
. To estimate
$\delta\kern-1.6pt B / B_0$
, we plot
$\sqrt {\langle \delta\kern-1.6pt B^2 / B_0^2\rangle }$
alongside the peak
$k \rho _{e0}$
in figure 16(b). The peak wavenumber was found using the same method as in § 3.2.1. It can be seen that once
$\sqrt {\langle \delta\kern-1.6pt B^2 / B_0^2\rangle } \gtrsim 2$
, the assumption that
$k \rho _{e} \sim k \rho _{e0} B_0 / \delta\kern-1.6pt B \sim 1$
applies well, again validating the proposed theory.
(a) Fourier spectra of
$B_y + \mathrm{i}B_z$
for 1D3V simulations with
$\beta _{e0} \in \{60,280, 1000\}$
and
$L_{\mathrm{T}0} / \rho _{e0}=100$
, averaged over the last
$2000 \varOmega _{e0}^{-1}$
of simulation time during the saturated interval and excluding data in the masking region. The spectra have been normalised such that
$\int \mathrm{d}k_\parallel \ E_B (k_\parallel ) = \langle \delta\kern-1.6pt B_y^2 +\delta\kern-1.6pt B_z^2 \rangle / B_0^2$
. The
$k_\parallel ^{-4}$
scaling typical of whistler waves is shown over its relevant region. (b) Plots of how the peak wavenumber (red) scales with
$\beta _{e0}$
, and how this compares with its expected values if the assumption
$k_\parallel \rho _e \sim 1$
is valid. Data are taken for 1D3V simulations with
$L_{\mathrm{T}0}/ \rho _{e0}=100$
. If
$k_\parallel \rho _e \sim 1$
, then
$k_\parallel \rho _{e0} \sim \delta\kern-1.6pt B / B_0$
, hence the actual quantity plotted with the purple error bars is
$\sqrt {\langle \delta\kern-1.6pt B^2 / B_0^2\rangle }$
as an estimate of this. All values have been averaged over the saturated interval and
$\sqrt {\langle \delta\kern-1.6pt B^2 / B_0^2\rangle }$
was spatially averaged over the small-
$L_{\mathrm{T}}$
region. The error bars come from the standard deviation of these quantities during the saturated interval (only the last
$2000 \varOmega _{e0}^{-1}$
of simulation time for the peak wavenumber).
Scalings of (a) the heat flux and (b) the perturbed magnetic energy with
$\beta _{e0}$
in 1D3V simulations for
$L_{\mathrm{T}0} / \rho _{e0}=100$
, where the values on the
$y$
axes have been spatially averaged over the small-
$L_{\mathrm{T}}$
region and temporally averaged over the saturated interval. A black dashed best-fitting line to the heat flux values is plotted. The expected values for both the perturbed magnetic energy and heat flux from the measured phase velocity are plotted, (multiplied by an order-unity constant for the perturbed magnetic energy). A best-fitting scaling relation for the perturbed magnetic energy has been plotted, based on the best-fitting scaling relation plotted in (a) and (2.26). The phase velocities and
$1\sigma$
error bars are calculated from the Gaussian fit to the phase velocity probability density function, as shown in figure 15(b).
Figure 17 shows the scalings for
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle$
and
$\langle q_{e\parallel } / q_{\mathrm{fs}} \rangle$
with
$\beta _{e0}$
.
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle$
and
$\langle q_{e\parallel } / q_{\mathrm{fs}} \rangle$
behave differently from their 2D3V counterparts, with
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle$
continuing to increase past
$\delta\kern-1.6pt B \sim B_0$
. Furthermore,
$\langle q_{e\parallel } / q_{\mathrm{fs}} \rangle \sim \beta _{e0}^{-0.5}$
, indicating weaker suppression than in two dimensions. Although the scaling is different from the 2D3V simulations, the results follow the proposed advective theory well. The heat flux and
$v_{\mathrm{ph}} / v_{\mathrm{th}e}$
agree reasonably well; their scalings with
$\beta _{e0}$
are very similar, and they quantitatively differ by only
$\sim$
30 %. The scaling
$\delta\kern-1.6pt B^2 / B_0^2 \sim \beta _{e0} v_{\mathrm{ph}} / v_{\mathrm{th}e}$
appears to hold within an order-unity factor as well. This supports the notion that the advective model does indeed apply in 1D3V simulations. A further aspect of figure 17 to note is the very steep dropoff in
$\delta\kern-1.6pt B^2 / B_0^2$
and enhancement of the heat flux for
$\beta _{e0}\lt 100$
, corresponding to simulations that are not in the ultra-high-
$\beta$
regime. This helps explain why significant heat-flux suppression in 1D3V simulations has not been observed in prior studies, as it only operates in the ultra-high-
$\beta$
regime.
3.3.2. Varying temperature gradient
As in § 3.2.2, the transition into the ultra-high-
$\beta$
regime was also studied in 1D3V by keeping
$\beta _{e0}$
constant and varying
$L_{\mathrm{T}0} / \rho _{e0}$
. Because 1D3V simulations are computationally cheaper,
$\beta _{e0}=180$
was used rather than
$\beta _{e0}=120$
. The temporal evolution of both the heat flux and perturbed magnetic energy is plotted in figure 18. As in the 2D3V case, all simulations behave very similarly, suggesting that the saturated WHFI becomes independent of
$L_{\mathrm{T}0}$
in this regime.
Temporal evolution of the spatially averaged parallel heat flux (upper) and perturbed magnetic energy (lower) for 1D3V simulations with background magnetic fields parallel to the temperature gradient and
$\beta _{e0}=180$
. Here
$L_{\mathrm{T}0} / \rho _{e0}$
is varied between 44 and 340 in order to vary
$\beta _{e0} \rho _{e0} / L_{\mathrm{T}0}$
without changing
$\beta _{e0}$
. The heat flux is normalised to the free-streaming value. A black line is drawn at
$\langle \delta\kern-1.6pt B^2 / B_0^2 \rangle = 1$
to show when the perturbed magnetic energy exceeds that of the background field.
Figure 18. Long description
Plots are shown for both the evolution of heat flux and perturbed magnetic energy for varying values of LT/ρe. All heat fluxes initially increase, before peaking and then rapidly decreasing. All plots saturate at similar values around 0.02. The perturbed magnetic energy initially increases very quickly, followed by a long period of slow growth and then saturation. The saturated value is the same regardless of temperature gradient, around 10.
Values of
$v_{\mathrm{ph}} / v_{\mathrm{th}e}$
for simulations with
$\beta _{e0}=180$
and varying
$L_{\mathrm{T}0} / \rho _{e0}$
. The mean phase velocity and
$1\sigma$
error bars were calculated in the same way as in figure 17.
Likewise, the dispersion relation appears to be independent of
$L_{\mathrm{T}0} / \rho _{e0}$
in both the linear and quadratic regions, as shown in figure 19. The phase velocity remains constant within the uncertainty, regardless of
$L_{\mathrm{T}0} / \rho _{e0}$
, as expected for a quantity that depends only on
$\beta _{e0}$
.
4. Discussion and conclusions
4.1. Key findings
In this paper, we have developed and tested a theory for the WHFI in ultra-high-
$\beta$
plasmas, defined by
$\beta _{e0}\rho _{e0}/L_{\mathrm{T}0}\gtrsim 1$
. In this regime, extrapolation of the moderately high-
$\beta$
theory predicts magnetic fluctuations with
$\delta\kern-1.6pt B/B_0\gtrsim 1$
, so the usual quasilinear, small-amplitude description is no longer self-consistent. We therefore proposed a large-amplitude saturation picture and heuristically derived two possible heat-flux closures: a diffusive model (§ 2.3.1) and an advective model (§ 2.3.2). Comparison with 2D3V and 1D3V PIC simulations shows that the advective model describes the simulations, whereas the diffusive model does not. Our main results are as follows.
-
(i) In 2D3V simulations, the magnetic energy in the whistler fluctuations can become larger than that in the background field but only by an order-unity factor, such that
$\delta\kern-1.6pt B \sim B_0$
is a reasonable approximation of the saturated state. Conversely, in 1D3V the amplification continues to increase over the range of
$\beta _{e0}$
investigated. -
(ii) The parallel heat flux is strongly suppressed relative to the free-streaming value. In 2D3V, we find
(4.1)while in 1D3V the suppression is weaker, approximately following
\begin{align} q_{e\parallel } \approx 4.7 \beta _{e0}^{-1} q_{\mathrm{fs}} \approx 2.65 \beta _{e0}^{-1} n_{e0} v_{\mathrm{th}e0}T_{e0}, \end{align}
$q_{e\parallel }/q_{\mathrm{fs}}\approx 0.3 \beta _{e0}^{-1/2}$
.
-
(iii) The heat flux is controlled by the whistler phase velocity, with
$q_{e\parallel }/q_{\mathrm{fs}}\sim v_{\mathrm{ph}}/v_{\mathrm{th}e}$
, as predicted by the advective model. -
(iv) The large-amplitude whistler fluctuations collectively form an effective transport barrier that traps or reflects most heat-carrying electrons. Thermal energy is therefore transported down the temperature gradient at approximately the whistler phase velocity rather than by spatial diffusion.
-
(v) Perpendicular transport remains strongly suppressed, even when the imposed temperature gradient has a component perpendicular to the background magnetic field. This is inconsistent with an isotropic diffusive heat-flux model.
-
(vi) In the ultra-high-
$\beta$
regime, the saturated heat flux and whistler phase velocity are approximately independent of
$L_{\mathrm{T}0}/\rho _{e0}$
at fixed
$\beta _{e0}$
within our numerical uncertainties, which is again consistent with the advective interpretation.
The finding that the WHFI-regulated heat flux is primarily advective in ultra-high-
$\beta$
plasmas, rather than diffusive, merits further discussion. The diffusive heat-flux model assumes that the electrons carrying most of the heat flux could undergo a random walk through the large-amplitude whistler fluctuations. Our simulations show otherwise; instead, their trajectories are such that the majority do not cross the transport barrier. Some electrons are deflected as they stream into the region of strongest fields, and subsequently stream away from the barrier. Others become trapped in local magnetic mirrors or locally closed-field line configurations. This implies that once the fluctuations become order-unity in amplitude, the transport mechanism changes: heat-flux suppression is no longer controlled primarily by quasilinear pitch-angle scattering, but by the motion of large-amplitude whistler structures that act as barriers to electron transport. Although an advective picture of WHFI-moderated heat transport has also been proposed for plasmas in the moderately high-
$\beta$
regime by Roberg-Clark et al. (Reference Roberg-Clark, Drake, Reynolds and Swisdak2018) and Drake et al. (Reference Drake, Pfrommer, Reynolds, Ruszkowski, Swisdak, Einarsson, Thomas, Hassam and Roberg-Clark2021), that picture differs substantively from the advective transport here. In those models, heat transport relies on cyclotron-resonant scattering of electrons by mobile scattering centres propagating at the whistler phase speed, a distinctive mechanism from the motion of large-amplitude structures that arise in the ultra-high-
$\beta$
regime considered here. In short, our results show that the transition to the ultra-high-
$\beta$
regime is not merely a quantitative extension of the moderately high-
$\beta$
WHFI, but instead involves a qualitative change.
4.2. Differences between 1D3V and 2D3V
Although the advective model captures the heat flux in both 1D3V and 2D3V simulations when the measured whistler phase velocity is used, the two geometries exhibit different dispersion relations, peak wavenumbers, and scalings with
$\beta _{e0}$
. This indicates that the dimensionality affects the nonlinear evolution of the whistler spectrum, even though the behaviour of the heat flux remains consistent with the advective model. Because the whistler waves are large-amplitude in the ultra-high-
$\beta$
regime, these differences are likely associated with nonlinear interactions among whistler modes, which need not proceed in the same way in 1D3V and 2D3V.
Fourier spectra of
$B_y + \mathrm{i}B_z$
at different times for (a) 1D3V
$\beta _{e0}=60$
, (b) 2D3V
$\beta _{e0}=60$
, (c) 1D3V
$\beta _{e0}=280$
and (d) 2D3V
$\beta _{e0}=280$
simulations. For the 1D3V plots, the Fourier spectra have been averaged over
$100 \varOmega _{e0}^{-1}$
around the time shown. For the 2D3V plots, the average was taken over
$40 \varOmega _{e0}^{-1}$
and the Fourier spectra are for both the parallel modes and for the spectra integrated over all
$k_y$
. No Fourier spectrum is shown in (b) for
$t \varOmega _{e0} = 2000$
, since the simulation did not run for that long. Fourier spectra have been normalised such that the integral of each spectrum equals
$\langle \delta\kern-1.6pt B_y^2 + \delta\kern-1.6pt B_z^2 \rangle / B_0 ^2$
, the perturbed magnetic energy perpendicular to the background magnetic field in the relevant modes at that time. All plots have the same axes limits for ease of comparison.
To interpret the difference in behaviour between the two types of simulations, Fourier spectra at different times during the simulation for different values of
$\beta _{e0}$
in both 1D3V and 2D3V are shown in figure 20. Several key differences are apparent. First, figure 20 shows that in one dimension, the peak wavenumber becomes larger for higher
$\beta _{e0}$
, which is not observed in 2D3V. This is related to the fact that the perturbed magnetic energy continues to increase with
$\beta _{e0}$
past
$\beta _{e0} \rho _{e0} / L_{\mathrm{T}0} \sim 1$
in 1D3V, decreasing the local electron Larmor radius
$\rho _e$
and hence increasing
$k_\parallel \rho _{e0}$
. Secondly, the fraction of power at larger wavenumbers is significantly higher in 1D3V than in 2D3V. Indeed, in 2D3V the power in the higher wavenumbers actually decreases in time initially, while it increases in 1D3V, suggesting that energy transfer from smaller towards larger wavenumbers is easier in 1D3V. Conversely, the proportion of power in small-wavenumber modes is much larger in the 2D3V simulations. This suggests that in two dimensions, large-wavenumber modes are suppressed and more energy is retained in, or transferred to, small-wavenumber modes. A natural interpretation is that this is due to the presence of oblique modes, which cannot form in 1D3V. The mechanism responsible for this redistribution is currently unclear and lies outside the scope of this paper. Finally, figure 20 shows that the height of the spectral peak is similar for the same
$\beta _{e0}$
, regardless of whether the simulation was run in 1D3V or 2D3V. This suggests that the continued growth of whistler waves beyond
$\delta\kern-1.6pt B \sim B_0$
in 1D3V is due to these high wavenumber modes, and that the peak wavenumber saturates with a similar power regardless of the simulation dimensionality.
To further investigate the nature of the nonlinear wave–particle interactions and how they differ between 1D3V and 2D3V simulations, we characterise the magnetic moment of particles during reflection events. For this characterisation, the velocities parallel to the background magnetic field of individual particles present in 1D3V and 2D3V simulations with
$\beta _{e0}=120$
were smoothed using a Savitsky–Golay filter (Savitzky & Golay Reference Savitzky and Golay1964) with window length corresponding to
$10 \varOmega _{e0}^{-1}$
to average out any Larmor motion. We then define a reflection event as the relative change in parallel velocity
$\Delta v_\parallel / v_{\parallel 0} \lt -1$
. The velocity of a particle between two reflection events
$v_{\parallel 0}$
was calculated by finding the extremum of the smoothed velocity between the two changes in sign, and averaging
$v_{\parallel }$
over
$10 \varOmega _{e0}^{-1}$
around that point. The change in velocity
$\Delta v_{\parallel }$
was calculated by finding the difference between two such consecutive values. The magnetic moment between two reflection events was calculated by taking the average magnetic moment
$\mu = m_e v_\perp ^2 / 2B$
over the same time interval as the velocity average, and the difference
$\Delta \mu$
was calculated in the same way.
We show the results of this analysis in figure 21, where the relative changes in magnetic moment and parallel velocity during reflection events are plotted. The magnetic moment is approximately conserved in 1D3V simulations, while it is not in 2D3V simulations. There are also clearly two distinct populations in the 1D3V figure, due to the different Doppler shifts experienced by electrons travelling towards the hot and cold wall, respectively, interacting with whistler waves. These plots suggest that, in 1D3V, large-reflection events conserve
$\mu$
well, consistent with magnetic-mirror-like reflection, whereas 2D3V simulations include other interactions, such as resonant scattering, which do not conserve
$\mu$
. As such, we conclude that the wave–particle interactions in 1D3V are distinct from those in 2D3V. This is, perhaps, unsurprising, as previous calculations have shown that resonant scattering of the heat-flux-carrying electrons cannot occur in 1D3V, whereas it can in 2D3V; however, this observation could explain the differences observed between the whistler dispersion relations for 1D3V and 2D3V simulations.
Probability-density maps for the relative changes in magnetic moment and velocity for particles that undergo reflection events, as defined in the text, in both (a) 1D3V and (b) 2D3V simulations. Both simulations have
$\beta _{e0}=120$
. In 1D3V, the magnetic moment is approximately conserved during these reflection events, whereas in 2D3V it is not.
The fact that the advective heat flux theory describes both 2D3V and 1D3V simulations raises the possibility that 1D3V simulations could be used as a reduced model for investigating the WHFI in ultra-high-
$\beta$
plasmas. This would substantially reduce the computational cost and make it possible to explore plasmas at higher
$\beta _{e0}$
compared with what can plausibly be investigated in 2D3V simulations. However, the quantitative behaviour in 1D3V differs from that in 2D3V;
$\delta\kern-1.6pt B^2 / B_0^2$
continues to increase beyond
$\delta\kern-1.6pt B \sim B_0$
, and the heat flux is less strongly suppressed. Thus, 2D3V simulations remain necessary for quantitatively reliable predictions of the WHFI in physical plasmas. Nevertheless, ultra-high-
$\beta$
1D3V simulations may be useful for quickly identifying parameter regimes in which interesting changes in behaviour occur, before more physically complete 2D3V simulations are performed.
4.3. Physical implications and applications
Perhaps the most important implication of the present work is that, despite the large electron mean free path, heat transport in the ultra-high-
$\beta$
regime appears to become effectively local once the WHFI saturates. Consequently, the advective heat-transport model supported by our simulations has a simple form that could be implemented in magnetohydrodynamic codes to better model collisionless, ultra-high-
$\beta$
plasmas. Indeed, the scaling relation found for parallel heat transport
$q_{e\parallel } \approx 4.7\beta _{e0}^{-1} q_{\mathrm{fs}}$
(where the prefactor
$4.7$
is the mean prefactor from the 2D3V simulations with angled and parallel magnetic fields) is surprisingly close to that found in Komarov et al. (Reference Komarov, Schekochihin, Churazov and Spitkovsky2018):
$q_{e\parallel } \approx 1.5 \beta _{e0}^{-1} q_0 \approx 3 \beta _{e0}^{-1} n_e v_{\mathrm{th}e} T_e \approx 5.3\beta _{e0}^{-1} q_{\mathrm{fs}}$
, despite the fact that the amplitude of the whistler waves becomes independent of
$\beta _{e0}$
in the ultra-high-
$\beta$
regime. Therefore, the same heat transport scaling with
$\beta _{e0}$
seems to be observed in all WHFI-unstable collisionless high-
$\beta$
regimes studied to date, regardless of the amplitude of the waves. This is a non-trivial result, as the very different mechanisms at play do not imply that this should be the case. For a typical front-side blow-off laser plasma with
$T_e \sim 1\,\mathrm{keV}$
and
$\beta _{e0}=100$
, the electron Larmor period is of the order
$\varOmega _e^{-1} \sim 1\times 10^{-13}\, \mathrm{s}$
. This paper has shown that the saturation time for the WHFI is
$\sim 10^3 \varOmega _e^{-1}$
, which corresponds to
$\sim$
0.1 ns in the laser plasma; however, the heat flux itself approaches a constant value more quickly, after just several hundred electron Larmor periods. This is short compared with the time scale for macroscopic evolution of the plasma, which is of the order of
$\sim$
1 ns. It should therefore be possible to implement a heat flux model into a magnetohydrodynamic code that follows
$q_{e\parallel } \approx 4.7 \beta _{e0}^{-1} q_{\mathrm{fs}}$
for regions of the plasma where the WHFI operates (i.e.
$\beta _e \lambda _e / L_{\mathrm{T}} \gg 1$
). The locality of the heat flux in this model makes it much more computationally manageable compared with current non-local methods (e.g. Luciani, Mora & Virmont Reference Luciani, Mora and Virmont1983; Schurtz, Nicolaï & Busquet Reference Schurtz, Nicolaï and Busquet2000).
The ICF plasmas are known to self-generate magnetic fields of order
$10$
–
$100$
MG, enough to render the plasma weakly magnetised with a Hall parameter greater than unity (Walsh et al. Reference Walsh, Chittenden, McGlinchey, Niasse and Appelbe2017; Walsh & Clark Reference Walsh and Clark2021; Frank & Bose Reference Frank and Bose2024, Reference Frank and Bose2025). Prior research has shown that even weakly collisional plasmas are kinetically unstable for high enough
$\beta _{e0}$
(Bott et al. Reference Bott, Cowley and Schekochihin2024; Lopez et al. Reference Lopez, Bott and Schekochihin2025); hence the WHFI is likely to be dynamically important in ICF plasmas. Understanding how this instability impacts the transport properties of ICF plasmas is therefore vital for accurate modelling and optimising experimental designs. For example, Walsh et al. (Reference Walsh2025a
) showed that, for premagnetised ICF implosions simulated with the classical Braginskii heat conduction model, the best-performing capsule design was distinct from the unmagnetised case, and ignition could be achieved at reduced drive energies. How WHFI-regulated heat transport affects the performance of different capsule designs remains to be seen. If the present scaling continues to hold for ICF-hotspot relevant parameters (
$\beta _e \sim 10^{3}$
–
$10^4$
), we would predict suppression of heat fluxes by up to an order of magnitude, suggesting that exploiting this instability could, in principle, improve hot-spot energy confinement and thereby increase neutron yield.
Outside of ICF research, these results may also be applicable to studies of the early universe in the collisionless, ultra-high-
$\beta$
plasma of the post reionisation intergalactic medium. Much research has been carried out in primordial plasmas to investigate different sources that may have contributed to the growth of magnetic fields observed in present-day galaxies (Turner & Widrow Reference Turner and Widrow1988; Vachaspati Reference Vachaspati1991; Ratra Reference Ratra1992; Widrow Reference Widrow2002; De Souza & Opher Reference De Souza and Opher2008; Durrer & Neronov Reference Durrer and Neronov2013; Subramanian Reference Subramanian2016; Vachaspati Reference Vachaspati2021; Jedamzik & Pogosian Reference Jedamzik and Pogosian2024). These plasmas are thought to be ultra-high-
$\beta$
and weakly collisional, with values of
$\beta _e$
as high as
$\beta _e \sim 10^{22}$
being possible (Barkana & Loeb Reference Barkana and Loeb2001; Bott et al. Reference Bott, Cowley and Schekochihin2024). As such, kinetic instabilities are likely to have been present, with potentially significant impacts on both transport properties and growth of the magnetic energy. In particular, this paper has shown that kinetic instabilities like the WHFI can significantly amplify magnetic energy. The 2D3V simulations showed significant magnetic field energy at scales significantly larger than the local electron Larmor radius, raising the possibility that the WHFI may act as an intermediary for transferring magnetic energy from electron to ion scales, before ion kinetic instabilities amplify the magnetic fields further.
Finally, the results of this paper have implications for laboratory-astrophysics experiments investigating anomalous heat conduction in weakly collisional plasmas (e.g. Vincent et al. Reference Vincent2026). In order for the results of such experiments to be interpreted correctly, it is important to understand whether the scaling relationships in the laboratory plasma and its astrophysical counterpart are the same. As this paper has shown, the transport mechanisms differs between plasmas with
$\beta _{e0} \rho _{e0} / L_{\mathrm{T}0} \gtrsim 1$
and those with
$\beta _{e0} \rho _{e0} / L_{\mathrm{T}0} \ll 1$
; hence any laboratory experiments investigating heat transport in astrophysical plasmas must ensure this parameter is in the appropriate regime for the phenomena they are studying.
4.4. Limitations and future work
While the present work establishes a theoretical framework for the WHFI and its regulation of heat transport in ultra-high-
$\beta$
plasmas, and demonstrates good agreement with kinetic simulations, several important questions remain unanswered. First, the simulations presented here are collisionless, which limits their direct applicability of the results to weakly collisional plasmas. Extending the present work to weakly collisional plasmas is therefore a natural next step. Previous research has shown that collisions damp, but do not necessarily prevent the formation of, kinetic instabilities (see Lopez et al. Reference Lopez, Bott and Schekochihin2025). It is therefore plausible that, in a weakly collisional plasma, the whistler waves would not attain the same amplitudes that they do in the collisionless case, and that the diffusive heat flux may become more significant due to the presence of collisions, although this remains to be investigated.
A second limitation is that the simulations presented here are 1D3V and 2D3V. Although the 2D3V simulations allow both parallel and oblique whistler modes, and can investigate cross-field transport in the simulation plane, they do not capture fully three-dimensional mode coupling or magnetic-field-line wandering. The latter effect might be expected to be important in ultra-high-
$\beta$
plasmas because the magnetic-field perturbations associated with the whistler fluctuations are comparable in size to the background magnetic field. Extending the present work to fully three-dimensional simulations is therefore important, though we caution that such simulations using current kinetic codes would have a high computational cost.
Another limitation is that we do not propose a complete theory that explains all of the quantitative properties of the WHFI in ultra-high-
$\beta$
plasmas that we observe in our simulations. For example, the exact process by which the whistler waves acquire a linear dispersion relation in the ultra-high-
$\beta$
regime is not studied here. Likewise, the physical origin of the difference between the dispersion relation in 1D3V simulations compared with 2D3V simulations is also not fully understood. This presents an interesting avenue for future research into nonlinear wave–wave and wave–particle interactions in high-
$\beta$
plasmas.
Another key question that remains unanswered is the upper limit on
$\beta _e$
at which the WHFI no longer regulates thermal transport. As unmagnetised plasmas are not susceptible to the WHFI, it should be possible to derive a bound for
$\beta _e$
, above which the heat flux is no longer suppressed as
$\beta _e$
increases further. In a collisionless plasma, it seems reasonable to assume that the WHFI can no longer operate when the Larmor radius becomes a finite fraction of the temperature length scale due to the insufficient scale separation between the characteristic whistler wavelength
$2\pi / k_{\textrm {ww}} \approx 2\pi / \rho _e$
and
$L_{\mathrm{T}}$
. This implies that an upper bound on
$\beta _e$
to destabilise the WHFI would be
$\beta _{e} \sim L_{\mathrm{T}}^2/d_e^2$
, where
$d_e = \beta _e^{-1/2} \rho _e$
is the electron skin depth. Investigating the exact nature of this upper bound for
$\beta _e$
is left to future work.
The WHFI, and kinetic instabilities more generally, remain comparatively unexplored in ultra-high-
$\beta$
plasmas, despite its relevance to high-energy-density plasmas, laboratory astrophysics and the early universe. In this paper, we have demonstrated that large-amplitude whistler fluctuations fundamentally alter the mechanisms supporting electron heat transport, replacing the quasilinear pitch of resonant pitch-angle scattering with one in which thermal energy is regulated by the propagation of large-amplitude magnetic fluctuations that act as effective transport barriers. Understanding the nonlinear plasma physics responsible for this transition, and understanding how robust it remains in weakly collisional plasmas, is essential for understanding and developing predictive heat-flux models in high-
$\beta$
plasmas encountered in both nature and the laboratory.
Acknowledgements
The authors would like to thank A. Fraser and A. Schekochihin for helpful and engaging discussions on this work.
Editor Luís O. Silva thanks the referees for their advice in evaluating this article.
Funding
R.D. was supported by the EPSRC and First Light Fusion under the AMPLIFI Prosperity Partnership – EP/X025373/1; P.P.C. and A.F.A.B. were supported by UKRI (grant number MR/W006723/1). Computing resources provided by STFC Scientific Computing Department’s SCARF cluster.
Declaration of interests
The authors report no conflict of interest.
θ
nppc
tend
ωp0
Te
βe0=180
tΩe0=60
tΩe0=1069
x
x
LT
LT0/ρe0=100
βe0
LT
⟨δB2/B02⟩=1
By+iBz
βe0=60
βe0=280
350Ωe0−1
kyρe0=0
kyρe0=1
ω/Ωe0=v^phk∥ρe0
k∥ρe0>0
By+iBz
LT0/ρe0=100
βe0=60
βe0=280
kyρe0=0
kyρe0=1
ky
∫dk∥ EB(k∥)=⟨δBy2+δBz2⟩/B02
350Ωe0−1
k∥−4
βe0
k∥EB(k‖)
350Ωe0−1
k∥ρe0∈{1,1.5}
βe0
LT0/ρe0=100
y
LT
⟨qe∥/qfs⟩=Aβe0−1
A≈4.8
β
βe0
vph∠
kyρe0=1
vph∥
kyρe0=0
βe0
βe0=20
LT
tΩe0=370
βe0=120
LT0/ρe0=100
LT
vph
βe0=120
LT0/ρe0
βe0ρe0/LT0
βe0
⟨δB2/B02⟩=1
vph/vthe
βe0=120
LT0/ρe0
kyρe0=0
kyρe0=1
350Ωe0−1
Te
tΩe0=54
tΩe0=570
βe0=400
βe0
LT0/ρe0=100
y
LT
⟨qe∥/qfs⟩=Aβe0−1
A≈4.6
β
qe∧
qe⊥
qfs
LT
βe0
1σ
Te
βe0=280
LT0/ρe0=100
tΩe0=1617
x
x
βe0
LT
⟨δB2/B02⟩=1
LT
LT
tΩe0=1200
βe0=120
LT0/ρe0=100
LT
vph
βe0=280
LT0/ρe0=100
2000Ωe0−1
By+iBz
βe0∈{60,280,1000}
LT0/ρe0=100
2000Ωe0−1
∫dk∥ EB(k∥)=⟨δBy2+δBz2⟩/B02
k∥−4
βe0
k∥ρe∼1
LT0/ρe0=100
k∥ρe∼1
k∥ρe0∼δB/B0
⟨δB2/B02⟩
⟨δB2/B02⟩
LT
2000Ωe0−1
βe0
LT0/ρe0=100
y
LT
1σ
βe0=180
LT0/ρe0
βe0ρe0/LT0
βe0
⟨δB2/B02⟩=1
vph/vthe
βe0=180
LT0/ρe0
1σ
By+iBz
βe0=60
βe0=60
βe0=280
βe0=280
100Ωe0−1
40Ωe0−1
ky
tΩe0=2000
⟨δBy2+δBz2⟩/B02
βe0=120
Open access
