below $\omega_{\rm A}\delta B_{\perp}^{2}/B_{0}^{2}$ and flattening its
spectrum, perhaps to what is observed. However, the most straightforward
estimate for this effect – that magneto-immutability would lead to $\Delta p$
fluctuations of approximately constant dynamical importance across scale (see
§ 2.2.3) – would suggest $\delta\Delta p\propto B_{0}^{2}$ or
$\mathcal{E}_{\Delta p}\sim\mathcal{E}_{p_{\|}}\sim k_{\perp}^{-1}$, which
also does not agree with what is observed (§ 4.3). So, there is no obvious
phenomenological explanation for the measured spectral slopes within this
framework. As described above, the observation that $\delta p_{\perp}\ll\delta
p_{\|}$ can be naturally explained by postulating that the finite-amplitude
Alfvénic fluctuations are in perpendicular pressure balance with the magnetic-
field-strength fluctuations, $\delta p_{\perp}\sim\delta(B^{2})$ (as observed
in § 4.3). Because $\delta p_{\|}$ is instead driven by
$\hat{\bm{b}}\hat{\bm{b}}\,\mbox{\boldmath{$:$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}$
to $\delta p_{\|}\sim\beta\max(1,\omega_{\rm A}/\nu_{\rm
c})\delta(B^{2})\gg\delta(B^{2})$, this suggests that $\delta\Delta
p\sim\delta p_{\|}\gg\delta p_{\perp}$.

The other possibility – a compressive magnetosonic cascade – more naturally
explains some spectral features, but appears to disagree with other key
properties. Its most problematic aspect is that oblique CGL magnetosonic slow
waves satisfy $\delta p_{\|}/\delta p_{\perp}\approx(5\beta+6)/2$ (see §
A.1.2), a far larger ratio than what is observed in all three simulations
shown in § 4.3.171717Heat fluxes and collisions modify $\delta p_{\|}/\delta
p_{\perp}$ modestly, but the predicted ratio remains far too large to match §
4.3 in all relevant regimes. Furthermore, although the observed ratio $\delta
p_{\|}/\delta p_{\perp}$ increases modestly with $\beta$, the increase is far
slower than the predicted linear dependence, which, for example, would yield
$\mathcal{E}_{p_{\|}}/\mathcal{E}_{p_{\perp}}\sim 1000$ in the _CL100_
simulation. This suggests that at least some other component or complication
is necessary to decrease $\delta p_{\|}/\delta p_{\perp}$ well below the
linear prediction to match the simulations. On the other hand, the observed
spectral scalings of the rates of strain and pressure seem to be most
naturally explained via a compressive cascade: if one postulates
$u_{\|}\propto u_{\perp}\sim k_{\perp}^{-1/3}$ (assuming the compressive and
Alfvénic cascades scale similarly),
$\hat{\bm{b}}\hat{\bm{b}}\,\mbox{\boldmath{$:$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}$
should scale as ${\sim}k_{\|}u_{\|}\propto k_{\perp}^{1/3}$, giving
$\mathcal{E}_{\nabla_{\|}u_{\|}}\propto k_{\perp}^{-1/3}$, which is close to
what is seen in figure 10. Predicted scalings of the other rates of strain
($\nabla_{\|}u_{\perp}\propto k_{\perp}^{1/3}$, $\nabla_{\perp}u_{\|}\propto
k_{\perp}^{2/3}$, and $\nabla_{\perp}u_{\perp}\propto k_{\perp}^{2/3}$) are
also consistent with what is observed, being the same as what is expected for
passive slow-mode cascade in isothermal MHD (and indeed, the passive-$\Delta
p$ rates of strain exhibit similar scalings in figure 10). In this scenario,
pressure spectra would result predominantly from the passive mixing of larger-
scale fluctuations, as opposed to local driving by
$\hat{\bm{b}}\hat{\bm{b}}\,\mbox{\boldmath{$:$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}$,
which also naturally yields the observed scaling
$\mathcal{E}_{p_{\|}}\propto\mathcal{E}_{p_{\perp}}\sim k_{\perp}^{-5/3}$
through the inertial range where $u_{\perp}\sim k_{\perp}^{-1/3}$, as well as
explaining the similarity of the collisionless and weakly collisional
simulations. Another interesting observation is the long parallel scale of
$p_{\|}$ and $u_{\|}$ fluctuations (see $\ell_{\|}(\lambda)$ in figure 8b),
which would be expected if for some reason the frequencies of Alfvénic and
compressive fluctuations were matched around the outer scale.

Overall, we see a plausible agreement of the Alfvénic hypothesis with most
diagnostics, although its predictions remain too qualitative to say much of
substance. Its biggest flaw is the similarity of the collisionless and weakly
collisional rate-of-strain and pressure spectra, which argues against pressure
fluctuations being driven locally in scale by
$\hat{\bm{b}}\hat{\bm{b}}\,\mbox{\boldmath{$:$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}$.
A compressive cascade alone cannot explain the observed spectra because the
$\delta p_{\|}/\delta p_{\perp}$ ratio characteristic of linear slow-mode
fluctuations is too large to fit the data (particularly the dependence on
$\beta$); however, we cannot rule out their importance in some form (e.g., if
$\delta p_{\perp}$ was dominated by residual Alfvénic fluctuations), and
various other properties, particularly spectral scalings, seem to be more
naturally explained via the compressive-cascade hypothesis. Clearly, more work
is needed to make further progress.

###  4.4 Heating and energy fluxes

Perhaps the most interesting and macroscopically relevant impact of the
pressure-anisotropy reduction through magneto-immutability is the suppression
of viscous (pressure-anisotropy) heating that it entails. As discussed above,
the effect implies that viscous damping is suppressed compared to naïve
estimates for individual Alfvénic perturbations of similar amplitudes to the
observed turbulent fluctuations, leading to a larger turbulent-cascade
efficiency, with more energy transfered to small scales. Consequently, a
larger fraction of the heating will occur through processes mediated by the
small-scale turbulent cascade (e.g., Landau damping and nonlinear phase mixing
around the ion gyroscale; Schekochihin et al., 2009; Kunz et al., 2018), as
opposed to viscous damping of outer-scale eddies. These different heating
mechanisms could in turn could have important thermodynamic consequences,
modifying, for instance, the perpendicular-to-parallel ion heating fraction,
or the ion-to-electron heating fraction.

In this section, we quantify these ideas by measuring directly the turbulent
transfers and fluxes between scales. This allows us to measure the local
damping due to pressure anisotropy across each scale in the cascade, and to
determine the cascade efficiency by measuring the fraction of input energy
that is lost into heat near the outer scale. This will again highlight a large
difference between active- and passive-$\Delta p$ simulations, demonstrating
that magneto-immutable turbulence is effectively conservative (no damping)
below the outer scale.

####  4.4.1 Energy transfers and fluxes

As discussed in § 3.3.3, the turbulence energetics can be usefully quantified
using ‘transfer functions’ $\mathcal{T}^{\rm AB}_{q\rightarrow k}$, which
measure the energy transfer between shells in Fourier space due to different
types of nonlinear interactions. Each of these transfers can involve several
terms in the compressible system, so, as detailed in Grete et al. (2017), we
combine relevant compressive terms into the non-compressive terms in order to
reduce the number of quantities under consideration. The compressive terms are
all individually very small and do not show interesting features, although
they are necessary to include in order to respect energy conservation. Their
full forms, as we compute them, are as defined in Grete et al. (2017), but we
simply use the placeholder $\mathcal{C}_{\rm AB}$ to refer to them below. We
follow Grete et al. (2017) in our nomenclature; ${\rm AB}={\rm UU}$ refers to
the transfer between shells of the kinetic energy through
$\bm{u}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}+\mathcal{C}_{\rm
UU}$ (see 27), ${\rm AB}={\rm BU}$ refers to the transfer from magnetic to
kinetic energy through
$\bm{B}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\bm{B}+\mathcal{C}_{\rm
BU}$, ${\rm AB}={\rm UB}$ refers to the transfer from kinetic to magnetic
energy through
$\bm{B}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}+\mathcal{C}_{\rm
UB}$, ${\rm AB}={\rm BB}$ refers to the transfer between shells of the
magnetic energy through
$\bm{B}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\bm{B}+\mathcal{C}_{\rm
BB}$, ${\rm AB}={\Delta p{\rm U}}$ refers to the transfer from the kinetic to
thermal energy via the pressure anisotropy through
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,(\hat{\bm{b}}\hat{\bm{b}}\Delta
p)$ (see below), and ${\rm AB}={\rm FU}$ refers to the transfer to kinetic
energy from the external forcing.181818We do not plot the pressure-stress
transfer from $-\mbox{\boldmath{$\nabla$}}p_{\perp}$, ${\rm AB}={\rm PU}$,
because it is small and uninteresting, although it is included when computing
sums of all terms.  We compute transfers as a function of the isotropic $k$ as
opposed to $k_{\perp}$, but results versus $k_{\perp}$ are similar.

As introduced in A+22, the only important new term compared to Grete et al.
(2017) is that due to the pressure-anisotropy stress
$\mathcal{T}_{q\rightarrow k}^{\Delta p{\rm U}}$. There is some freedom in the
definition of this; A+22 use

$\mathcal{T}_{q\rightarrow k}^{\Delta
p\mathrm{U}}=\int\mathrm{d}\bm{x}\left\langle\sqrt{\rho}\bm{u}\right\rangle_{k}\,\mbox{\boldmath{$\cdot$}}\,{\rm
sign}({\Delta p})\frac{\sqrt{|\Delta
p|}}{B}\frac{\bm{B}}{\sqrt{4\pi\rho}}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\langle\sqrt{|\Delta
p|}\hat{\bm{b}}\rangle_{q},$ (33)

which has the advantage that the square of $\langle\sqrt{|\Delta
p|}\hat{\bm{b}}\rangle_{k}$ can be interpreted as the pressure anisotropy’s
contribution to the thermal energy (Kunz et al., 2015), but the disadvantage
of sign discontinuities that arise from $\sqrt{|\Delta p|}$. If we instead
interpret the pressure-anisotropy stress as a damping of kinetic energy (as
opposed to the transfer between energy reservoirs) a more natural definition
is

$\mathcal{T}_{q\rightarrow k}^{\Delta
p\mathrm{U}}=\int\mathrm{d}\bm{x}\left\langle\sqrt{\rho}\bm{u}\right\rangle_{k}\,\mbox{\boldmath{$\cdot$}}\,\frac{\bm{B}}{\sqrt{4\pi\rho}}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\left\langle\frac{\Delta
p}{B^{2}}{\bm{B}}\right\rangle_{q}.$ (34)

We have computed both versions, finding qualitatively similar results, but
will focus on (34), because it fits more naturally with our focus on pressure-
anisotropy damping. Note also that the influence of the sign-discontinuity
issue of 33 was mitigated in A+22 because the $\Delta p$ distribution was
mostly confined to negative values, but the spread of $\Delta p$ is somewhat
broader in our simulations (see figures 3 and 4).

Figure 11: Important terms in the net energy-transfer spectra, as defined in
§ 3.3.3 and 28. The top panels show the simulations with $\beta_{0}=10$,
$\nu_{\rm c}=0$ and the bottom panels show the simulations with
$\beta_{0}=100$, $\nu_{\rm c}=33$, with the upper and lower subpanels for each
showing the active- and passive-$\Delta p$ cases, respectively. Insets in the
plots of the two active-$\Delta p$ cases show a zoom of the grey-box region.
All transfers are normalised to $\varepsilon/(k/k_{f})$, which approximately
represents the local cascade rate (see text). The different colours show
different energy-transfer terms as labeled in the top panel:
$\mathrm{T}^{\Delta p{\rm U}}$ is the transfer from $\bm{u}$ to thermal energy
due to
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,(\hat{\bm{b}}\hat{\bm{b}}\Delta
p)$; $\mathrm{T}^{\rm UU}$ is momentum transfer via
$\bm{u}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}$;
$\mathrm{T}^{\rm BU}$ is the contribution to kinetic energy from magnetic
tension and pressure; $\mathrm{T}^{\rm BB}$ is the advection of $\bm{B}$ in
the induction equation; $\mathrm{T}^{\rm UB}$ is the contribution to magnetic
energy from field stretching; and $\mathrm{T}^{\rm FU}$ is the forcing
contribution. In the active-$\Delta p$ simulations below the outer scale of
the turbulence, we see little contribution of $\Delta p$ to the energy
transfer, which shows directly that the effect of the pressure-anisotropy
stress is minimised below the outer scale, even though the local interruption
number is still smaller than unity in this range (see figure 7). In contrast,
in the passive-$\Delta p$ simulations (bottom panel), the pressure anisotropy
that develops is such that, if it did feed back on the flow, it would cause
$\mathcal{O}(1)$ dissipation at all scales in the cascade (blue line), which
would completely damp the turbulence.

Figure 12: The energy fluxes (see 29), which measure the transfer of energy
across a particular $k$, for the two active-$\Delta p$ simulations with
$\beta_{0}=10$, $\nu_{\rm c}=0$ (top panel) and $\beta_{0}=100$, $\nu_{\rm
c}=33v_{\rm A}/L_{\perp}$ (bottom panel). We normalise to the forcing input
$\varepsilon$, which implies that the total flux would be $\Pi/\varepsilon=1$
across the full inertial range for a conservative cascade with no damping. As
also seen in § 4.4.1, although there is some energy loss due to $\Delta p$ at
the outer (forcing) scale, at smaller scales the flux remains constant,
showing that there is little energy damping through the inertial range.

In § 4.4.1, we plot the net energy transfers ${\rm T}^{\rm AB}(k)$ (see 28),
which measure the net contribution from each term to the rate of change of the
kinetic or magnetic energy spectrum at each $k$. For a conservative cascade,
these are zero, because the transfer into a given shell from larger scales is
balanced by the transfer out of the shell to smaller scales. Thus, ${\rm
T}^{\rm AB}(k)$ provides a direct measure of the damping of kinetic or
magnetic energy at each scale due to pressure anisotropy or other terms. To
plot the results clearly, we normalise ${\rm T}^{\rm AB}(k)$ by the
dimensional (Kolmogorov) estimate $\partial_{t}\mathcal{E}\sim
ku_{k}\mathcal{E}\sim\varepsilon(k/k_{0})^{-1}$ (where
$k_{0}=2\pi/L_{\perp}$); with this normalisation, a line that is constant and
nonzero with $k$ symbolises a contribution that is of approximately constant
importance compared a conservative cascade towards smaller scales. As above,
we show the _CL10_ and _B100_ simulations, comparing to the results from the
corresponding passive-$\Delta p$ simulations (lower subpanels). The thick
black lines, which show the sum of all terms (excluding ${\rm T}^{\Delta p{\rm
U}}$ for the passive runs), are approximately zero throughout the inertial
range, as expected, with no clear difference between the active and passive
runs. While the contribution from ${\rm T}^{\Delta p{\rm U}}$ is slightly
negative in both _CL10_ and _B100_ , indicating a slight damping of inertial-
range motions, it is small compared to the cascade rate. Contrast this to the
equivalent ${\rm T}^{\Delta p{\rm U}}$ from the passive simulations, which
gives an indication of what the damping would be in the absence of feedback
from $\Delta p$. Although this measurement is counterfactual – the cascade
could not have proceeded in the presence of such strong damping – it concisely
illustrates how motions driven by normal turbulence would be strongly damped,
at a rate comparable to the cascade rate, across all scales. This is not
entirely obvious for the collisionless case particularly, because the
contribution to the kinetic energy from
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,(\hat{\bm{b}}\hat{\bm{b}}\Delta
p)$ can be both positive and negative (it is negative definite in the
Braginskii-MHD regime, but not otherwise; see 10). The effect of magneto-
immutability all but eliminates such damping below the outer scale, leaving
only a modest damping of outer-scale motions for the weakly collisional case.

§ 4.4.1 shows similar information in a different form by plotting the
contributions to the turbulent energy flux (29). We show only the
active-$\Delta p$ runs because the comparison with passive simulations is less
interesting in this case. Unlike the net transfers (§ 4.4.1), individual
terms in the fluxes are not easily interpreted because they do not include the
diagonal $\mathcal{T}^{\rm AB}_{k\rightarrow k}$ transfer. While this diagonal
contribution necessarily vanishes for the total energy flux, it dominates the
transfer between different energy reservoirs, including the damping of the
flow by the pressure anisotropy (blue line). Thus the fact that $\Pi^{\Delta
p{\rm U}}$ is small and positive is not particularly relevant, while the
individual lines are of interest only insofar as they indicate a contribution
to the total flux (black lines). The most important feature that we observe is
that the total energy fluxes are nearly constant for $10\lesssim k\lesssim
200$, consistent with the result of § 4.4.1 that there is little pressure-
anisotropy damping through the inertial range (a slight decrease in $\Pi(k)$
is observed for _B100_ , indicating a slight damping). The value of the
$\Pi/\varepsilon$ in the inertial range is a direct measure of the cascade
efficiency, which, as expected, is less than unity because there is some
pressure-anisotropy damping near the forcing scales (the reduction of
${\simeq}25\%$ for C100 and ${\simeq}50\%$ for _B100_ is consistent with the
total inferred pressure anisotropy heating, which is discussed below)

####  4.4.2 Pressure-anisotropy heating fraction and cascade efficiency

The net energy transfer ${\rm T}^{\Delta p{\rm U}}(k)$ also gives a convenient
way to evaluate the total heating rate due to pressure anisotropy. Because the
pressure-anisotropy (viscous) heating is confined to the outer scales (§
4.4.1), this is also a measure of the cascade efficiency – the energy that is
available to cascade in the usual way to heat via small-scale collisionless
mechanisms. Its precise value depends on the forcing scheme, since it mostly
occurs at the outer scales where the forcing drives the flow, and could thus
differ more significantly in systems with more realistic turbulence generation
mechanisms (e.g., magneto-rotational turbulence; Kunz et al., 2016; Kempski et
al., 2019). Nonetheless, its dependence on parameters ($\beta$ and the
collisionality) is interesting to consider, as is a comparison to A+22.

We define

$\mathcal{H}_{\Delta p}=\int_{k_{0}}^{k_{\rm max}}{\rm d}k\,{\rm T}^{\Delta
p{\rm U}}(k),$ (35)

so $\mathcal{H}_{\Delta p}/\varepsilon$ measures the fraction of the energy
input that is converted into heat via pressure anisotropy up to wavenumber
$k=k_{\rm max}$. Because the grid-scale dissipation in our simulations is non-
physical – it is supposed to represent energy absorbed by the smaller scales,
which would eventually cascade to mostly dissipate around the ion gyro-radius
scale (Schekochihin et al., 2009) – it is most appropriate to choose $k_{\rm
max}$ to lie near the bottom of the inertial range. We set $k_{\rm max}=100$
as appropriate for the resolution $400\times 200^{2}$, but the choice hardly
affects the results for the active-$\Delta p$ simulations anyway, because
almost all of the viscous heating is confined near the outer scales.191919For
passive-$\Delta p$ cases, $\mathcal{H}_{\Delta p}$ increases with $k_{\rm
max}$ through the inertial range because there is more and more damping (see
§ 4.4.1).

Figure 13: Net pressure-anisotropy heating rate $\mathcal{H}_{\Delta p}$
normalised to the cascade rate $\varepsilon$ from all simulations. The cascade
efficiency, defined as the fraction of energy that participates in the cascade
to small scales, is $1-\mathcal{H}_{\Delta p}/\varepsilon$. This is computed
using the net transfer functions (§ 4.4.1) by summing up the contributions
from all $k\leq 100$, _viz.,_ $\mathcal{H}_{\Delta p}=\sum_{k\leq
100}\mathrm{T}^{{\rm U}\Delta p}(k)$, so as to exclude the grid scales. The
marker style denotes the simulation set (see § 3.2 and table 1): circles
denote the _lrB_ set at constant $\mathcal{I}$, crosses denote the _lrCL_ set
with $\nu_{\rm c}=0$, and stars denote the $\beta 16$ set. We see that viscous
(pressure-anisotropy) heating damps up to ${\simeq}45\%$ of the cascade in
collisional regimes ($\nu_{\rm c}/\omega_{\rm A}\gtrsim 1$), and less in
collisionless cases. This value is approximately constant for
$\mathcal{I}\lesssim 10$ and independent of $\beta$ and $\nu_{\rm
c}/\omega_{\rm A}$ (so long as $\nu_{\rm c}/\omega_{\rm A}\gtrsim 1$). The
small symbols show the same computation from the passive-$\Delta p$
simulations; in this case the heating rate is larger than unity, because
$\Delta p$ does not actually feed back on the flow. Such large values imply
that if $\Delta p$ had not been reduced by magneto-immutability, it would have
nearly completely damped the cascade.

We evaluate $\mathcal{H}_{\Delta p}$ for all of the _lrCL_ , _lrB_ , and
$\beta 16$ simulations (see table 1), time averaging across the steady state
of each. Results are shown in figure 13, plotted against the measured
interruption number, with marker colour indicating the collisionality and
marker style indicating the simulation set (_lrCL_ , _lrB_ , or $\beta 16$).
The CGL-LF (active-$\Delta p$) simulations, shown with large symbols, have
$\mathcal{H}_{\Delta p}$ between $\mathcal{H}_{\Delta p}\simeq 0.2\varepsilon$
for the collisionless simulations and $\mathcal{H}_{\Delta p}\simeq
0.45\varepsilon$ for the weakly collisional and Braginskii-MHD simulations at
$\mathcal{I}\sim 1$, with no significant dependence on $\beta$. In other
words, for this choice of forcing, the cascade efficiency is always above
${\simeq}50\%$, meaning most of the input energy can participate in an MHD-
like turbulent cascade. Contrast this with the small symbols, which show the
passive-$\Delta p$ simulations for _lrCL_ and _lrB_ , indicating what the
pressure-anisotropy heating rate would be in the absence of magneto-
immutability. The numerical values are of course meaningless in this case –
the situation is counterfactual and having $\mathcal{H}_{\Delta
p}>\varepsilon$ is clearly not possible – but, as above, it demonstrates that
the motions involved in standard MHD Alfvénic turbulence would be strongly
damped by pressure anisotropy were it not for the pressure-anisotropy feedback
modifying their structure. From the $\beta 16$ simulations, we see also that
$\mathcal{H}_{\Delta p}(\mathcal{I})$ remains relatively large up to large
$\mathcal{I}$, despite the effective driving of pressure anisotropy
fluctuations decreasing towards large $\mathcal{I}$ due to their larger
$\nu_{\rm c}$. This shows that as $\mathcal{I}$ is decreased below
${\simeq}10$, the suppression of viscous heating is sufficiently strong that
it balances the stronger driving of $\Delta p$ fluctuations that results from
lower $\nu_{\rm c}$ (otherwise $\mathcal{H}_{\Delta p}$ would continue
increasing below $\mathcal{I}\simeq 10$).

###  4.5 Effect of electron pressure

In all of the simulations presented so far, we have set $p_{e}=\rho T_{e}=0$,
thus ignoring the influence of electrons. These add an additional isotropic
pressure response that is worthy of exploration in case it interferes with
magneto-immutability once it starts to dominate over the anisotropic ion
pressure response. We explore this in figure 14 via the rate-of-strain
spectra, motivated by the finding above that these showed a clear difference
between CGL-LF and standard MHD (passive) simulations. We show the low-
resolution versions of the spectra for the _CL10_ case already shown in figure
10, comparing them with two additional simulations that include electrons with
$T_{e}/T_{i0}=1$ and $T_{e}/T_{i0}=5$ (here $T_{i0}$ is the initial ion
temperature). We see no important differences, even for $T_{e}\gg T_{i}$,
which further implies that the cascade efficiency is independent of
$T_{e}/T_{i}$. This is expected since the plasma is already almost
incompressible anyway, and an additional isotropic pressure should simply help
it be even more so. The simulations of S+19 and Kempski et al. (2019) used a
fully incompressible Braginskii-MHD model and found similar results.

Figure 14: Rate-of-strain spectra at $\beta_{0}=10$, $\nu_{\rm c}=0$,
comparing the effect of including an (isothermal) electron pressure with
different $T_{e}$. We show $T_{e}=0$ (solid line), $T_{e}/T_{i0}=1$ (dot-
dashed lines), $T_{e}/T_{i0}=5$ (dotted lines), and the passive-$\Delta p$
simulation (dashed lines). The plasma’s behaviour is almost identical to the
$T_{e}=0$ case, because the isotropic electron pressure has little effect once
motions are already nearly incompressible.

The isotropic fluid electron model is formally valid either when the plasma is
semi-collisional (e.g., the ICM) or when electrons are cold. For the case of
hot collisionless electrons, a separate electron pressure-anisotropy equation
must be solved, which will add an electron contribution to the pressure
anisotropy. Presumably this will generically act to enhance the influence of
pressure anisotropy, even at large scales, while also bringing in a plethora
of electron-scale kinetic instabilities driven by pressure-anisotropies, heat
fluxes, and other non-Maxwellian features (e.g., Gary & Nishimura, 2003;
Verscharen et al., 2022). However, because the Coulomb interspecies
temperature equilibration is extremely slow and (as far as we are aware) there
exist no plasma instabilities that act to decrease it (e.g., Zhdankin et al.,
2021), there is not an obvious connection between the ion and electron
thermodynamics. Thus, even though electron instabilities could strongly
influence the electron temperatures (e.g., by limiting pressure anisotropies
and heat fluxes; Verscharen et al. 2022), insofar as they remain very small-
scale and unable to interact directly with ions, we do not expect them to have
a strong influence on the large-scale ion dynamics other than through the
pressure-anisotropy stress.202020Perhaps, by enhancing the electron scattering
rate, electron instabilities could improve the validity of a collisional
electron closure. That said, there is clearly more work needed to understand
this complex physics better.

##  5 Discussion: uncertainties and observable effects

The numerical results presented herein have demonstrated how pressure
anisotropy feeds back on Alfvénic turbulent flows at high $\beta$, minimising
its own influence to reduce the variance of $B^{2}$ and the associated
parallel viscous heating. The basic effect is neither surprising nor
escapable, arising from the large pressure-anisotropy forces that are rapidly
generated with any change of $B$ and act to oppose this change, in the same
way that the large isotropic pressure forces that result from changes to
density render fluids incompressible. However, our more detailed conclusions
and quantitative features (e.g., the differences between collisional and
collisionless results) do depend on the assumptions of the CGL-LF with hard-
wall pressure-anisotropy limiters. In this section, we consider these in more
detail by reference and comparison to the hybrid-kinetic simulations of A+22,
particularly focusing on uncertainties relating to the interplay of magneto-
immutability and microinstability scattering. We argue that most results are
robust and compare well to A+22, which is promising more generally for the use
of CGL-LF models in the study of weakly collisional plasmas. We also discuss
various diagnostics that can be used in observations or kinetic simulations,
where the comparison to passive-$\Delta p$ simulations – so useful in our
analysis above – is not available.

###  5.1 Interplay between magneto-immutability and microinstabilities

Perhaps the largest uncertainty concerning the CGL-LF model relates to the
influence of particle scattering through microinstabilities. As discussed in
§ 2.3, the ‘hard-wall-limiter’ method used here and in previous works is
based on the idea that large-scale turbulent motions are extremely slow
compared to microinstability-saturation time scales. This suggests that
microinstabilities should instantaneously halt the growth of $|\Delta p|$ when
it is driven beyond the stability thresholds, then instantaneously disappear
again when the driving reverses and the plasma returns towards stability. In
its practical implementation, the hard-wall-limiter model achieves this by
taking $\nu_{\rm c}^{\rm lim}$ to be very large in unstable regions; however,
kinetic theory and simulations show that the saturated scattering rate should
actually be $\nu_{\rm c}\sim S\beta$, where
$S\sim\hat{\bm{b}}\hat{\bm{b}}\,\mbox{\boldmath{$:$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}$
is the shear driven by the turbulence. The effect on $\Delta p$ is similar –
both hard-wall limiters and scattering with $\nu_{\rm c}\sim S\beta$ act to
pin $\Delta p$ at the stability thresholds – but the differences could lead to
inaccuracies.212121Most importantly, the physics of particles moving spatially
between stable and unstable regions of the plasma can, at least in principle,
be captured through the heat fluxes in the CGL-LF model; but, with large
$\nu_{\rm c}^{\rm lim}$, the heat fluxes are completely suppressed in unstable
regions meaning such physics is missed. Nevertheless, while kinetic
simulations find that heat fluxes are well approximated by the collisional
estimate in firehose regions, in mirror regions, they seem to be nearly
completely suppressed (Kunz et al., 2020). Thus, it is likely that a more
accurate model ought to treat $\Delta p>0$ and $\Delta p<0$ regions
differently, meaning that the complete suppression of heat fluxes via a large
$\nu_{\rm c}^{\rm lim}$ in limited regions may be no worse than using
$\nu_{\rm c}\sim S\beta$ anyway. Furthermore, kinetic simulations have shown
that the assumption that microinstabilities saturate and decay instantaneously
compared to the turbulent motions is questionable. This is most acute for the
mirror instability, which can limit $\Delta p$ growth via particle trapping
through much of a macroscopic shearing time (Schekochihin et al., 2008; Kunz
et al., 2014; Rincon et al., 2015). Decaying firehose and mirror fluctuations
can also be extremely long lived, continuing to scatter particles even when
not actively driven by the large-scale shear (Squire et al., 2017b; Kunz et
al., 2020). Indeed, their decay rate, as it involves changing $B^{2}$, may
itself be limited by the need to avoid creating unstable pressure anisotropies
(Melville et al., 2016). Thus, it may not be appropriate to consider
scattering only in unstable regions – rather, the continual driving will
create a ‘soup’ of magnetic fluctuations that drive scattering across the
entire plasma, putting it in the weakly collisional or Braginskii regime (Ley
et al., 2023). This situation applies to almost any foreseeable kinetic
numerical simulation, including those of A+22, but is at odds with the
assumptions of our collisionless simulations, where we set $\nu_{\rm c}=0$ in
stable regions.

We are thus left with some uncertainty between two limiting scenarios: in the
first, scattering sites pervade the plasma; in the second, they exist only in
spatiotemporal locations where the plasma is unstable. Reality, for
astrophysically relevant scale separations, likely lies somewhere between
these extremes.222222The scale separation in the solar wind, which is
${\simeq}10^{4}$ between the outer scales and $\rho_{i}$, is definitely not
large enough to lie in the second regime, being foreseeably approachable with
kinetic simulations (e.g., Bott et al., 2021). Comparing the microinstability-
scattering estimate above, $\nu_{\rm c}\sim\beta S\sim\beta\omega_{\rm
A}\mathcal{M}_{\rm A}^{2}$, with the Alfvénic interruption estimates discussed
in § 2.2, we see that the $\nu_{\rm c}$ expected in the first scenario is
exactly the collisionality required to maintain $\mathcal{I}\simeq 1$. Indeed,
the correspondence is inherent in the estimate, since this $\nu_{\rm c}$ is
simply that which is necessary to suppress the production of $\Delta p$ to a
level where its influence becomes comparable to the magnetic forces. This
shows that if the first scenario applies, all collisionless plasmas
effectively become weakly collisional with $\nu_{\rm c}/\omega_{\rm
A}\sim\mathcal{M}_{\rm A}^{2}$ and $\mathcal{I}\gtrsim 1$, meaning our
$\mathcal{I}<1$ collisional and collisionless simulations are never formally
applicable (see further discussion in A+22). More generally, in both
scenarios, the scattering from microinstabilities and the pressure anisotropy
stress (magneto-immutability) must be of similar importance to the plasma’s
dynamics. Both effects limit $\Delta p$ to $|\Delta p|\lesssim B^{2}/4\pi$,
and scale in the same way with $\beta$ and turbulence amplitude. As a
corollary, there is no limit in which one or the other can be ignored, except
in the presence of mean pressure anisotropy (resulting in microinstabilities
but with no corresponding pressure-anisotropy stress; Bott et al., 2021).

In order to provide a more detailed theory, we must estimate the relative
importance of two effects that scale in the same way with key parameters. This
is difficult to do without kinetic simulations at asymptotically large scale
separation, currently not feasible. The kinetic simulations of A+22 exhibited
a box-averaged scattering rate consistent with $\nu_{\rm c}\simeq\beta S$ at
$\beta=16$, or $\mathcal{I}\simeq 1$ (see their figure 12).232323Their
$\beta=4$ simulation exhibited a larger scattering rate, which does not fit
with this scaling, but that simulation also had several other peculiarities
that may explain this discrepancy. While consistent with the first scenario
(‘scattering everywhere’), given the modest separation (a factor ${\simeq}20$
between the outer scales and $\rho_{i}$), this does not constitute strong
evidence against the second scenario in general, as there is no reasonable way
to separate stable and unstable regions in the later stages of their
simulations. This measured scattering rate, $\nu_{\rm c}\simeq\beta S$, also
shows that either our $\beta 16\nu 6$ or $\beta 16\nu 12$ simulation should be
directly comparable to the results of A+22 (see table 1). Indeed, for scales
above $\rho_{i}$, most diagnostics appear relatively similar (e.g., pressure
anisotropy spectra, the standard deviation of $B$ and $\Delta p$, and the
viscous heating rate), with perhaps the most obviously significant large-scale
difference being that the kinetic simulation drives itself to negative mean
pressure anisotropy $\langle\Delta p\rangle$, while our $\beta 16$ simulations
have $\langle\Delta p\rangle\approx 0$. A likely explanation for this lies in
the crudeness of the scattering process from $\nu_{\rm c}$ in the CGL-LF
model. In the kinetic simulation, negative pressure anisotropy builds over
time as a result of the preferentially parallel heating due to Landau damping,
seemingly reaching steady state once it triggers firehose scattering across
the box (see their figures 3 and 5). However, once particle scattering from
firehose fluctuations is triggered, rather than driving the plasma towards
isotropy, it should be expected to maintain the plasma near the firehose
threshold, allowing the system to maintain $\langle\Delta p\rangle<0$ even
with a relatively high scattering rate. In contrast, the scattering included
in our weakly collisional simulations with the CGL-LF model obviously drives
the system towards $\langle\Delta p\rangle=0$, making it much harder for the
system to maintain $\langle\Delta p\rangle<0$ in the weakly collisional case.

Given all of these uncertainties, our approach in this paper has been to
explore a range of options within the confines of the fluid model. We have
found that the pressure-anisotropy stress has a strong influence on the
plasma’s dynamics across all collisionality regimes, for $\mathcal{I}\lesssim
10$ (see, e.g., figures 13 and 16). Thus even if the first scenario applies at
astrophysically relevant scale separations, magneto-immutability should play
an important role in the plasma’s dynamics, allowing a nearly conservative
turbulent cascade to be set up below the forcing scales. Promisingly, A+22
measured a cascade efficiency of $55$–$60\%$ (meaning that $55$–$60\%$ of the
input energy is processed via the turbulent cascade to kinetic scales), which
compares very well to all of our weakly collisional and Braginskii-MHD
simulations with $\mathcal{I}\lesssim 10$ (all of which have similar
$\mathcal{H}_{\Delta p}$; see figure 13). Given the significant differences in
physics and forcing, the correspondence could hint at ${\simeq}50\%$ being a
quasi-universal (or at least minimum) cascade efficiency for weakly
collisional plasmas, although an exploration of more realistic forcing
mechanisms would be needed to confirm this. Overall, the use of the semi-
phenomenological fluid model has been both a blessing and a curse: while we
clearly miss very important kinetic effects, leading to the uncertainties
discussed above, the influence of different physics can be explored more
thoroughly by selectively studying different options, spanning wider parameter
ranges, and comparing with counterfactual scenarios.

###  5.2 Observable consequences and diagnostics

Much of the evidence for magneto-immutability that we presented in this paper
has relied on the direct comparison to the equivalent passive-$\Delta p$ or
MHD simulations, a diagnostic luxury that is clearly not available in
observations or even in kinetic simulations. In this section, we use the
$\beta 16$ simulation set, which scans between the collisionless and MHD-like
regimes by changing $\nu_{\rm c}$, to suggest various methods for diagnosing
magneto-immutability in simulations or spacecraft data. There are several
difficulties in this endeavour. First, the general effect of magneto-
immutability is to make the weakly collisional turbulent plasma look more
similar to MHD than it otherwise would; the corollary is that those
differences that do persist are rather subtle. Secondly, there are other
effects in MHD that tend to reduce the variance of $B$, which are unrelated to
pressure anisotropy but very important in the near-Sun (low-$\beta$) solar
wind especially. Thirdly, in either the solar wind or kinetic simulations, we
will often not know _a priori_ the effective particle scattering rate (see,
e.g., Hellinger et al., 2011; Coburn et al., 2022), making it hard to
disentangle the relative contributions of scattering and magneto-immutability
to the reduction in the variance of $|\Delta p|$ (see § 5.1 above).

####  5.2.1 Constant-$B$ states in the solar wind – imbalance and spherical
polarisation

_In-situ_ spacecraft measurements in the solar wind and magnetosphere provide
us with an unparalleled laboratory for studying fundamental plasma physics,
and especially collisionless plasma turbulence (Chen et al., 2016; Verscharen
et al., 2019). Regions with $\beta\gtrsim 1$ are regularly observed in the
bulk solar wind, with values as high as $\beta\sim 10^{2}$ in specific regions
(e.g., Cohen et al., 2017; Chen et al., 2021), making it a promising arena for
studying the physics of pressure anisotropy and magneto-immutability. However,
when working in this laboratory, we must recognise and consider other possible
physics at play. In this context, a key point to note is that constant-$B$
‘spherically polarised’ fluctuations, which are observed ubiquitously near the
Sun (e.g., Belcher & Davis, 1971; Kasper et al., 2019), are most likely not
related to pressure anisotropy and magneto-immutability. Here we discuss these
and other possible mechanisms for reducing the variance of $|\bm{B}|$, which
should be taken into account when interpreting observational data.

The robustness of spherically polarised states is related to the fact that any
perturbation that satisfies

$\frac{\delta\bm{B}}{\sqrt{4\pi\rho}}=\pm\bm{u},\quad B={\rm const},\quad
p_{\perp},p_{\|},\rho={\rm const}$ (36)

is an exact nonlinear solution of the Kinetic MHD system (1)–(5) (or even the
standard MHD equations). Any such solution propagates at velocity
$\mp\bm{B}_{0}\sqrt{\Uptheta/4\pi\rho}$, where $\bm{B}_{0}$ is the mean of
$\bm{B}$ and $\delta\bm{B}$ is the remainder (Barnes & Hollweg, 1974). Since
the solution (36) holds even if $|\delta\bm{B}|\gg|\bm{B}_{0}|$, these states
are the natural nonlinear generalisation of linear Alfvén waves to large
amplitudes. Their ubiquitous presence in the near-Sun solar wind is not
unexpected given that Alfvénic perturbations grow in amplitude as they
propagate outwards from the low corona due to the decreasing Alfvén speed
(e.g., Völk & Aplers, 1973; Hollweg, 1974b). Given that they are observed
consistently when $\beta\ll 1$, which necessarily implies $\mathcal{I}\gg 1$
for trans-Alfvénic motions, it seems unlikely that pressure anisotropy is
playing an important role in their formation and sustenance (but see Tenerani
& Velli 2020). Indeed, most aspects of their evolution can be well described
by isothermal MHD (Hollweg, 1974b; Vasquez & Hollweg, 1998; Squire et al.,
2020; Mallet et al., 2021).

Clearly, as exact nonlinear solutions, such states are not meaningfully
turbulent, although there do exist turbulent-like states with very small $B$
variance that are similar to (36) (Dunn et al., 2023). But in order for such
states to occur, the turbulence must be strongly _imbalanced_ , with
$|\bm{z}_{+}|=|\bm{u}+\bm{B}/\sqrt{4\pi\rho}|\gg|\bm{z}_{-}|=|\bm{u}-\bm{B}/\sqrt{4\pi\rho}|$
(or vice versa). This is the norm in the fast solar wind, but is less
prevalent in the slow wind – presumably turbulence is also generically less
imbalanced in other astrophysical plasmas with less ordered global structure.
In this paper, we have not considered imbalanced turbulence, taking
$|\bm{z}_{+}|\sim|\bm{z}_{-}|$ in all simulations, a limitation of our study
that should be relaxed in future work. At any rate, the basic message is that
there exist both MHD-related and pressure-anisotropy-related effects that
reduce $B$ fluctuations in turbulence. In order to diagnose the influence of
pressure anisotropy by considering the variance of $B$ between different
regions, the imbalance must be taken into account, lest one inadvertently
study the prevalence of nearly spherically polarised solutions instead.
Similarly, in highly imbalanced high-$\beta$ turbulence, both pressure
anisotropy and other reasons (e.g., wave growth in a plasma with a decreasing
Alfvén speed) could cause the plasma to develop constant-$B$ spherically
polarised fluctuations. Once formed, such fluctuations would not generate
significant turbulent pressure anisotropy, in which case magneto-immutability
might not play such an important role even for $\mathcal{I}\lesssim 1$.

Figure 15: The variation in magnetic-field strength, or spherical
polarisation, for all simulations, quantified in terms of $C_{B^{2}}$ given by
37. This measures the relative tendency for the components of $\bm{B}$ to be
correlated in order to reduce the variation in $|\bm{B}|$. We plot multiple
snapshots for each simulation, with the marker colour showing $\mathcal{I}$,
which changes modestly during simulations due to heating and random
fluctuations. Marker styles denote the simulation set and the small markers
show the equivalent passive-$\Delta p$ simulations, with the same styles as in
figure 13. We see a clear correlation between the magneto-immutability and the
interruption number, even within individual simulations, and without any clear
dependence on the collisionality regime. Further, even for $\beta\sim 1$,
$C_{B^{2}}$ is reduced significantly by pressure-anisotropy effects. However,
even the lowest values of $C_{B^{2}}$ observed here are larger than those
regularly observed in the near-Sun solar wind due to the predominance of
highly imbalanced, spherically polarised states 36.

To demonstrate more directly the reduction in the variance of $B$ across our
simulations, in figure 15 we plot

$C_{B^{2}}=\frac{(\delta B^{2})_{\rm rms}}{(\delta\bm{B}_{\rm
rms})^{2}}=\frac{\langle(|\bm{B}|^{2}-\langle|\bm{B}|^{2}\rangle)^{2}\rangle^{1/2}}{\langle|\bm{B}-\langle\bm{B}\rangle|^{2}\rangle},$
(37)

which was introduced by Squire et al. (2020) as a simple measure of the
relative tendency of the components of $\bm{B}$ to become correlated in order
to reduce the variance of $B$. This measure is mostly insensitive to the
amplitude of the turbulence (small-amplitude, linearly polarised Alfvén waves
have $C_{B^{2}}\approx 1/2$).242424A number of other reasonable choices exist
to quantify features of the spherical polarisation. The similar statistic
$C_{B}=(\delta B_{\rm rms})^{2}/(\delta\bm{B}_{\rm rms})^{2}$ decreases with
turbulence amplitude, meaning it is more appropriately interpreted as a total
magnetic compressibility (Chen et al., 2020). Another possibility, which may
be preferable because it is likely to depend less strongly on the
fluctuations’ intermittency, is $\widetilde{C}_{B^{2}}=(\delta B^{2})_{\rm
rms}/\langle|\bm{B}-\langle\bm{B}\rangle|^{4}\rangle^{1/2}$. This produces
very similar results to 37, aside from all points being reduced by a factor
${\approx}2$.  For reference, Parker Solar Probe (PSP) observations and highly
imbalanced expanding-box MHD simulations have $C_{B^{2}}$ in the range of
$0.1$ to $0.3$ (Squire et al., 2020; Johnston et al., 2022). We plot
$C_{B^{2}}$ from a number of snapshots of the $\beta 16$, _lrCL_ , and _lrB_
simulation sets, with each point coloured by the interruption number
$\mathcal{I}$ (this can vary modestly between snapshots in a given
simulation). As in figure 13, the MHD (passive-$\Delta p$) simulations are
plotted with small markers, showing the expected higher variance of $B$ for
otherwise similar conditions (see figures 3 and 4). The dependence on
$\mathcal{I}$ is obvious, with lower-$\mathcal{I}$ simulations generally
showing lower $C_{B^{2}}$ at otherwise identical parameters. We also see a
tendency for modestly lower $C_{B^{2}}$ at lower $\beta$, which is likely
related to the larger influence of magnetic pressure forces compared to
thermal pressure forces at low $\beta$ (Vasquez & Hollweg, 1998). However,
also of note is that all simulations, including those where pressure-
anisotropy forces are very strong (low $\mathcal{I}$), show values of
$C_{B^{2}}$ that are rather large compared to what is routinely observed in
the solar wind or imbalanced turbulence (Kasper et al., 2019; Squire et al.,
2020). It seems that magneto-immutability, while an important influence on the
evolution of $B$, can never drive balanced turbulence arbitrarily close to
keeping $B$ truly constant across the domain – imbalanced spherically
polarised solutions are much better at this task.

####  5.2.2 Spectral diagnostics

Figure 16:  The left panel shows the ratio of pressure spectra
$\mathcal{E}_{p_{\|}}/\mathcal{E}_{p_{\perp}}$, a simple diagnostic of
magneto-immutability that should be readily observable in the solar wind. We
show the simulations from the $\beta 16$ simulation set (see § 3.2 and table
1), with the line colour showing $\mathcal{I}$, and the black lines showing
the corresponding results for isothermal MHD. The dashed line shows the _lrB_
600, $\beta_{0}=600$ Braginskii-MHD simulation, to demonstrate that the
observed decrease in $\delta p_{\|}/\delta p_{\perp}$ is not due to the change
from the weakly collisional to the Braginskii-MHD regime at high $\mathcal{I}$
and $\beta_{0}=16$. The right panel shows the ratio of the spectra of $u_{\|}$
and $u_{\perp}$ in the parallel and perpendicular direction;
$\mathcal{E}_{\nabla_{\|}u_{\perp}}/\mathcal{E}_{\nabla_{\|}u_{\|}}$ is shown
with solid lines, and
$\mathcal{E}_{\nabla_{\perp}u_{\perp}}/\mathcal{E}_{\nabla_{\perp|}u_{\|}}$ is
shown with dashed lines (see figure 10 and § 3.3.1). The increase in
$\nabla_{\|}u_{\perp}/\nabla_{\|}u_{\|}$ with decreasing $\mathcal{I}$ shows
the change in the structure of the flow that results from the turbulence
becoming more magneto-immutable at smaller collisionality. The lack of the
same change in $\nabla_{\perp}u_{\perp}/\nabla_{\perp|}u_{\|}$ allows the
effect to be considered separately from $u_{\|}$ itself, which may depend on
the forcing and other parameters.

In § 4, we saw that, while kinetic- and magnetic-energy spectra remained
rather similar in magneto-immutable turbulence compared to MHD, the rate-of-
strain spectra and pressure anisotropy spectra did not. The $\beta 16$
simulation set provides a useful scan from low-$\mathcal{I}$ (small $\nu_{\rm
c}$) to large-$\mathcal{I}$ collisional MHD, while other parameters
($\mathcal{M}_{\rm A}$ and $\beta$) remain similar. In figure 16 we plot two
interesting ratios of spectra that clearly show the change in the structure of
the turbulence at low $\mathcal{I}$. In the left panel, we show the ratio of
parallel- to perpendicular-pressure spectra, which should be readily
observable with _in-situ_ solar-wind measurements. As discussed in § 4.3, it
becomes larger than unity in magneto-immutable turbulence because of
perpendicular pressure balance. We see the clear trend with $\mathcal{I}$: the
higher-collisionality simulations ($\nu_{\rm c}/\omega_{\rm A}\gtrsim 50$,
$\mathcal{I}\gtrsim 8$) all collapse down to $\delta p_{\|}\sim\delta
p_{\perp}$, because the pressure anisotropy ceases to have a strong dynamical
influence. In addition, the dashed line shows the _lrB_ 600 simulation, which
is in the Braginskii regime but with $\mathcal{I}\approx 0.4$. This verifies
that the observed change is not related to a transition between collisionality
regimes (collisionless, weakly collisional, or Braginskii MHD). A similar
diagnostic could be constructed from the slopes of the pressure spectra: in
low-$\mathcal{I}$ turbulence, $\mathcal{E}_{p_{\perp}}$ flattens with scale
while $\mathcal{E}_{p_{\|}}$ steepens, with the two spectra converging at
small scales (see § 4.3); in contrast, when unaffected by the pressure-
anisotropy feedback, the $p_{\perp}$ and $p_{\|}$ spectra have similar shapes.

Since we saw that the reduction in $\nabla_{\|}u_{\|}$ was one of the more
obvious consequences of the pressure anisotropy feedback, the right panel
compares the ratios of parallel-flow gradients
$\mathcal{E}_{\nabla_{\|}u_{\perp}}/\mathcal{E}_{\nabla_{\|}u_{\|}}$ (solid
lines) to those of perpendicular-flow gradients
$\mathcal{E}_{\nabla_{\perp}u_{\perp}}/\mathcal{E}_{\nabla_{\perp}u_{\|}}$
(dashed lines). While the latter ratio is almost constant as a function of
$\mathcal{I}$ and looks very similar to MHD (black lines), the former varies
by a factor of at least ${\sim}10$ as the plasma transitions to collisional
MHD. We also see that even at the highest collisionality ($\nu_{\rm
c}/\omega_{\rm A}=200$, $\mathcal{I}\approx 50$), there remain clear
differences in the flow structure compared to MHD.

In summary, these diagnostics could provide a useful way to understand and
measure the influence of pressure-anisotropy feedback in more realistic
simulations or observations.252525Rate-of-strain spectra are not possible to
measure with a single spacecraft, but could be computed from data taken by
multi-spacecraft constellations such as the upcoming NASA MidEx mission
HelioSwarm (Klein & Spence, 2021). While the exact values of these metrics
will depend on details (e.g., the forcing, $\mathcal{M}_{\rm A}$, and the
imbalance), we expect that the qualitative features should be robust,
especially the fact that
$\mathcal{E}_{\nabla_{\|}u_{\perp}}/\mathcal{E}_{\nabla_{\|}u_{\|}}\gg\mathcal{E}_{\nabla_{\perp}u_{\perp}}/\mathcal{E}_{\nabla_{\perp}u_{\|}}$
and $\delta p_{\|}\gg\delta p_{\perp}$. This may allow one, for example, to
diagnose and understand a transition to $\mathcal{I}\gg 1$ turbulence at small
scales (as seen, e.g., for _CL10_ in § 4.3).

##  6 Conclusions

In this paper, we have studied the influence of fluctuation-generated pressure
anisotropy ($\Delta p$) on magnetised plasma turbulence. The simulations and
theory are based on a CGL-Landau-fluid model, which is effectively drift
kinetics (Kulsrud, 1983) with heat fluxes approximated by a fluid closure that
captures linear parallel Landau damping (Snyder et al., 1997). This makes it a
reasonable model for collisionless plasma dynamics on scales far above the
proton gyroradius. Indeed, our results compare well to the recent hybrid-
kinetic simulations of A+22, which capture a much more comprehensive array of
kinetic processes at the price of having a limited inertial range. Our results
focus on the $\beta>1$ regime, where small changes of the magnetic-field
strength lead to large pressure-anisotropy stresses on the plasma. This regime
is relevant to many hot and diffuse astrophysical plasmas, including the
intracluster medium and regions of the solar wind.

We show numerically, and argue theoretically, that the primary effect of the
pressure-anisotropy stress on the fluid is to reduce its own influence via a
dynamical feedback. This makes high-$\beta$ weakly collisional turbulence
behave more like standard collisional MHD than would be naïvely expected,
aside from a few important differences. The effect is straightforwardly
understood by analogy with the origin of incompressibility in fluids, which
results from the large $-\mbox{\boldmath{$\nabla$}}p$ force that rapidly
opposes any flow with
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,\bm{u}\neq 0$
attempting to change the fluid’s pressure. With pressure and density tied
together, this feedback eliminates density fluctuations. Analogously, there is
a large force from
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,(\hat{\bm{b}}\hat{\bm{b}}\Delta
p)$ that rapidly opposes any flow with
$\hat{\bm{b}}\hat{\bm{b}}\,\mbox{\boldmath{$:$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}\neq
0$ attempting to change the pressure anisotropy; since $\Delta p$ and $B$ are
tied together, the feedback minimises magnetic-field-strength fluctuations in
the turbulence. In the same way that the Mach number of a hydrodynamic flow
quantifies the importance of pressure compared to inertial forces, the
‘interruption number’ $\mathcal{I}$ (introduced in S+19) quantifies the
importance of pressure-anisotropic forces compared to magnetic tension (see §
2.2.3). The most important consequence of this physics is to reduce viscous
heating from pressure anisotropy, confining it to a small range of scales near
the outer scale so that most of the cascade remains nearly conservative. This
increase in ‘cascade efficiency’ will in turn increase the fraction of heating
that occurs via kinetic processes near and below the ion gyroscale
(Schekochihin et al., 2009), thus influencing the thermodynamics of the
plasma.

Compared to the previous work on Braginskii-MHD Alfvénic turbulence (S+19;
Kempski et al., 2019), the results here explore plasmas in more realistic
regimes at more modest $\beta$ and lower collisionality. Our results highlight
the resilience of the pressure-anisotropy feedback mechanism: despite the much
more complex evolution of the pressure anisotropy in the regimes studied here,
particularly the effect of heat fluxes that act quasi-diffusively to smooth
out $\Delta p$, we find that the basic features of magneto-immutable
turbulence are robust. Our detailed analysis is aided by direct comparisons to
‘passive-$\Delta p$’ simulations, which are identical in setup but
artificially remove the feedback of the pressure anisotropy on the flow. This
allows us to compare, for example, the pressure-anisotropy distribution and
viscous heating to the counterfactual situation in which $\Delta p$ evolved in
standard MHD turbulence. A familiar illustration of this comparison is given
in figure 1, where we plot the classic ‘Brazil plot’ of the joint PDF of
$\beta$ and pressure anisotropy (e.g., Kasper et al., 2002; Hellinger et al.,
2006). In the CGL-LF simulations, the pressure anisotropy is naturally
constrained close to zero by magneto-immutability; in the passive-$\Delta p$
simulations, most of the volume sits instead at the microinstability
boundaries, enforced here by artificial limiters.

A primary result of this study is that magneto-immutable turbulence behaves
surprisingly similarly to standard MHD, illustrating the fundamental
robustness of the MHD turbulent cascade. Specifically, despite the system
having to rearrange itself to avoid dissipating most of the injected energy,
the simulations show that saturated fluctuation amplitudes are similar to MHD
for identical forcing, and the magnetic- and kinetic-energy spectra have
slopes close to the expected $k_{\perp}^{-3/2}$ or $k_{\perp}^{-5/3}$. This
has the effect of making the influence of pressure anisotropy rather subtle:
rather than it causing some dramatic, obvious modifications to the system,
magneto-immutability interferes and causes such turbulence to resemble MHD
rather closely. That said, there are important differences that enable this
behaviour, particularly that the system significantly suppresses the spectrum
of the parallel strain
$\hat{\bm{b}}\hat{\bm{b}}\,\mbox{\boldmath{$:$}}\,\mbox{\boldmath{$\nabla$}}\bm{u}$
(see figure 10), which reduces the magnitude of the $\Delta p$ and $B$
fluctuations (e.g., figures 4 and 15) and steepens the pressure spectra,
leaving a characteristic signature with $\delta\Delta p\sim\delta
p_{\|}\gg\delta p_{\perp}$ (§ 4.3). These changes are summarised in figure
16, which shows how turbulence transitions from being magneto-immutable (low
$\mathcal{I}$, blue) to collisional MHD (high $\mathcal{I}$, red and black),
demonstrating how such effects could be diagnosed and studied with _in-situ_
solar-wind observations (see § 5.2).

Macroscopically, the most important effect of magneto-immutability is to
suppress strongly the viscous heating that would otherwise almost completely
damp the turbulent cascade. This leaves an effectively conservative cascade
below the forcing scales that processes ${\approx}50\%$ to ${\approx}80\%$ of
the input energy (depending on the regime; see §§ 4.4.1 and 13). While the
particular values for the cascade efficiency will depend on the forcing
scheme, so will certainly differ for more realistic systems (e.g., those
driven by large-scale instabilities), the general point is that the magneto-
immutability modifies the conversion between mechanical and thermal energy in
the plasma. This could have interesting consequences for its macroscopic
thermodynamics. For example, important quantities like the fraction of
turbulent energy that heats ions versus electrons, or in the perpendicular
versus parallel directions, can depend on the cascade efficiency and viscous
heating. A specific application where such physics could have direct
consequence is the intracluster medium, in which turbulent heating is thought
to offset a large fraction of radiative losses to mitigate cluster cooling
flows (e.g., Churazov et al., 2004; Zhuravleva et al., 2014). Kunz et al.
(2010) suggested that, by forcing the plasma to adjust its parallel rate of
strain, pressure anisotropy could regulate turbulent heating to maintain
naturally the plasma’s thermal stability against Bremsstrahlung cooling (with
viscous heating ${\propto}B^{4}T^{-5/2}$ at fixed gas pressure increasing
faster than cooling ${\propto}T^{-3/2}$ as the temperature drops). While the
details of their scheme differ from the turbulent heating studied here
(specifically, Kunz et al. 2010 assumed that kinetic microinstabilities
provide the primary regulation mechanism by contributing to, and thereby
regulating, the total parallel rate of strain), the basic idea – that
suppression of the parallel rate of strain provides a strong constraint on the
heating rate – is analogous to magneto-immutability. Thus, similar thermal-
stability considerations may apply to magneto-immutable turbulence, with the
modification that it is the large-scale (as opposed to the small-scale)
adjustment to the parallel rate of strain that regulates the viscous heating.
At the same time, the mechanism would allow some of the injected energy (the
cascade efficiency) to pass into a turbulent cascade, as needed to explain
observations of small-scale fluctuations on scales below the Coulomb mean-free
path (Zhuravleva et al., 2019; Li et al., 2020). The effectiveness of such a
mechanism likely depends on how fluctuations are driven (the forcing), as well
as complications such as anomalous particle scattering (see § 5.1), so
further study in more realistic settings is clearly needed.

Due to the ad-hoc nature of some approximations that were needed for the CGL-
LF model, our study is beset with some serious (but, we think, ultimately not
game-changing) uncertainties. The most important one is the influence of
kinetic microinstabilities, e.g., the firehose and mirror instabilities, which
grow extremely rapidly compared to large-scale motions in regions of large
$|\Delta p|$. We use simplified hard-wall limiters (Sharma et al., 2006), but
it is clear from kinetic simulations that these could miss important effects
such as long-lived fluctuations that continue to scatter particles even in
stable regions. In this context, our study complements the recent hybrid-
kinetic simulations of A+22, which capture all such complexities (except the
electron physics), but suffer from unavoidable limitations related to low
scale separation. This complicates the exploration of the inertial-range
$k_{\perp}\rho_{i}\ll 1$ turbulence in the system, both because the
microinstabilities’ growth and saturation are not particularly fast compared
to the integral-scale motions, and because $\rho_{i}$-scale physics (e.g., the
damping of perturbations; Foote & Kulsrud, 1979) can start to interfere
directly with scales comparable to the size of the box. Nonetheless, our
comparable simulations that have an imposed collisionality to approximate the
effects of microinstabilites (see § 5.1), provide a good match to most
features observed in A+22. This includes most spectra of Alfvénic and
compressive quantities, $\Delta p$ and $B$ distributions (aside from a mean
$\langle\Delta p\rangle<0$ in A+22), and the cascade efficiency, which is
${\simeq}55\%$ in A+22 (matching effectively all of the weakly collisional and
Braginskii-MHD simulations; see figure 13). Thus, a secondary result of our
work is that the CGL-LF model provides an accurate and relatively simple way
to model weakly collisional plasmas in high-$\beta$ regimes. That said, there
remain significant caveats, and working to build a better bridge between the
regimes accessible via the CGL-LF and hybrid-kinetic approaches should be a
priority for future work.

###### Acknowledgements.

We thank Chris Chen, Ryan Davis, Sergey Komarov, Stephen Majeski, and Romain
Meyrand for useful and illuminating discussions, as well as the anonymous
referees who provided helpful comments that improved the presentation of this
article. Support for J.S. was provided by Rutherford Discovery Fellowship
RDF-U001804, which is managed through the Royal Society Te Apārangi. M.W.K.
was supported in part by NSF CAREER Award No. 1944972. Support for L.A. was
provided by the Institute for Advanced Study. The work of A.A.S. was supported
in part by grants from STFC (ST/W000903/1) and EPSRC (SP/R034737/1), as well
as by the Simons Foundation via a Simons Investigator award. EQ was also
supported in part by a Simons Investigator award from the Simons Foundation.
We wish to acknowledge the use of New Zealand eScience Infrastructure (NeSI)
high performance computing facilities as part of this research. New Zealand’s
national facilities are provided by NeSI and funded jointly by NeSI’s
collaborator institutions and through the Ministry of Business, Innovation &
Employment’s Research Infrastructure programme. We also wish to acknowledge
the generous hospitality of the Wolfgang Pauli Institute, Vienna, where these
ideas were discussed during several ‘Plasma Kinetics’ workshops.

##  Appendix A Numerical method and validation

In this appendix we discuss various features of the numerical implementation
of our equations. We first give the relevant equations and other results that
are needed for the conservative part of the equations (i.e., for the CGL
equations; § A.1), followed by those needed for the heat fluxes and
collisions (§ A.1.1). We then present some simple wave tests and discuss a
collection of numerical problems that arose for this system, speculating on
their possible causes. As part of this discussion, we present a more accurate
Riemann solver that works well for some problems, but has turned out to be too
unreliable for large-amplitude (large-$\delta B_{\perp}$) turbulence
simulations.

###  A.1 Conservative form

As is standard for finite-volume implementations, we solve 1–5 in the
conservative form

$\frac{\partial\bm{U}}{\partial t}+\frac{\partial\bm{F}_{x}}{\partial
x}+\frac{\partial\bm{F}_{y}}{\partial y}+\frac{\partial\bm{F}_{z}}{\partial
z}+\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,\bm{Q}=0,$ (38)

where

$\bm{U}=\begin{pmatrix}\rho\\\ m_{x}\\\ m_{y}\\\ m_{z}\\\ E\\\ \mathcal{A}\\\
B_{x}\\\ B_{y}\\\ B_{z}\end{pmatrix}=\begin{pmatrix}\rho\\\ \rho u_{x}\\\ \rho
u_{y}\\\ \rho u_{z}\\\
p_{\perp}+\frac{1}{2}p_{\|}+\frac{1}{2}B^{2}+\frac{1}{2}\rho|\bm{u}|^{2}+p_{e}\ln
p_{e}\\\ \rho\ln(p_{\perp}\rho^{2}/p_{\|}B^{3})\\\ B_{x}\\\ B_{y}\\\
B_{z}\end{pmatrix};$ (39)

$\bm{F}_{x}$, $\bm{F}_{y}$, and $\bm{F}_{z}$ are the CGL fluxes; and $\bm{Q}$
are the heat fluxes, which are defined below (for simplicity of notation, we
normalise $\bm{B}$ by $\sqrt{4\pi}$ through this section). In the limit of
$T_{e}=0$, the $p_{e}\ln p_{e}$ contribution to $E$ disappears. There is a
fundamental degeneracy in the choice of the second conserved thermodynamic
variable (in addition to energy), since both $\mu=p_{\perp}/\rho B$ and
$\mathcal{J}=p_{\|}B^{2}/\rho^{3}$ are passively advected by CGL dynamics. Our
choice of the logarithm of their ratio $\mathcal{A}$ was inspired by Santos-
Lima et al. (2014) and also by extensive numerical testing of solvers that use
$\rho\mu$ instead, which can become highly unstable for certain types of
problems. We suspect the superiority of $\mathcal{A}$ compared to $\mu$ is
related to the fact that $\mathcal{A}$ treats $p_{\perp}$ and $p_{\|}$ on
approximately equal footing, rather than assigning artificial importance to
one over the other (see further discussion below). Note that, in principle,
$\rho$ multiplied by any function of $\mu$ and/or $\mathcal{J}$ could be used
as the conserved variable, and it is possible that there exist other options
with better numerical properties.

The fluxes are

$\bm{F}=\begin{pmatrix}\rho u_{x}\\\ \rho
u_{x}^{2}+P^{\star}-B_{x}^{2}\Uptheta\\\ \rho u_{x}u_{y}-B_{x}B_{y}\Uptheta\\\
\rho u_{x}u_{z}-B_{x}B_{z}\Uptheta\\\
(E+P^{\star})u_{x}-B_{x}(\bm{u}\,\mbox{\boldmath{$\cdot$}}\,\bm{B})\Uptheta\\\
u_{x}\mathcal{A}\\\ 0\\\ B_{y}u_{x}-B_{x}u_{y}\\\
B_{z}u_{x}-B_{x}u_{z}\end{pmatrix},$ (40)

with the obvious permutation of $x$, $y$, and $z$ used to get $\bm{G}$ and
$\bm{H}$. Here $P^{\star}\doteq p_{\perp}+B^{2}/2+p_{e}$ and $\Uptheta\doteq
1+{\Delta p}/{B^{2}}$ is the anisotropy parameter.

####  A.1.1 Heat fluxes and collisions

A downside of using $\mathcal{A}$, as opposed to $\mu$, is that the heat
fluxes take a rather complex form that does not seem possible to write as a
total divergence
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,\bm{Q}$. In order to
maintain this property, thus simplifying the numerical implementation of the
heat fluxes, we transform $\mathcal{A}$ to $\mu$, add the heat fluxes to $e$
and $\mu$, then transform back to $\mathcal{A}$. For this purpose, the
$e$-component of the heat flux is $Q_{e}=\hat{\bm{b}}(q_{\perp}+q_{\|}/2)$,
while the $\mu$-component is $Q_{\mu}=\hat{\bm{b}}q_{\perp}/B$, where
$q_{\perp}$ and $q_{\|}$ are given by 7 and 8, respectively (all other
components of $\bm{Q}$ are zero). Further information is provided in § 3.1.

Collisions, including microinstabilities, are evaluated implicitly on the
primitive variables $p_{\perp}$ and $p_{\|}$ at the end of each global time
step. This is done via the exact solution of $\partial_{t}p_{\perp}=-(\nu_{\rm
c}/3)\,\Delta p$, $\partial_{t}p_{\|}=(2\nu_{\rm c}/3)\,\Delta p$ across time
$\delta t$:

$\displaystyle p_{\perp}(t+\delta t)=\frac{1}{3}(2+{\rm e}^{-\nu_{\rm c}\delta
t})p_{\perp}(t)+\frac{1}{3}(1-{\rm e}^{-\nu_{\rm c}\delta t})p_{\|}(t),$
$\displaystyle p_{\|}(t+\delta t)=\frac{2}{3}(1-{\rm e}^{-\nu_{\rm c}\delta
t})p_{\perp}(t)+\frac{1}{3}(1+2{\rm e}^{-\nu_{\rm c}\delta t})p_{\|}(t).$ (41)

The microinstability limiters are also implemented implicitly, via the first-
order in $\delta t$ solution of $\partial_{t}p_{\perp}=-(\nu_{\rm c}^{\rm
lim}/3)\,(\Delta p-\Lambda_{\rm MI}B^{2}/8\pi)$,
$\partial_{t}p_{\|}=(2\nu_{\rm c}^{\rm lim}/3)\,(\Delta p-\Lambda_{\rm
MI}B^{2}/8\pi)$, applied only in regions where $\Delta p$ lies beyond the
relevant instability threshold (here $\Lambda_{\rm MI}=-2$ and $\Lambda_{\rm
MI}=1$ for the firehose and mirror thresholds, respectively). The methods
relax $\Delta p$ back to isotropy or to the relevant instability threshold
accurately and with no constraint on the time step, meaning that adiabatic MHD
can be recovered by setting $\nu_{\rm c}\rightarrow\infty$, and true ‘hard
wall’ limiters are recovered by setting $\nu_{\rm c}^{\rm
lim}\rightarrow\infty$.

####  A.1.2 Dispersion relation

The HLL Riemann solver and its relatives (Toro, 2009), which are used to
approximate the solutions of 38 across cell boundaries, require the wave
speeds for the hyperbolic double-adiabatic system ($\bm{Q}=0,\,\nu_{\rm
c}=0$). Assuming a mode of wavenumber $\bm{k}=k\bm{\hat{x}}$ and linearizing
about the equilibrium, $\rho=\rho_{0}$, $\bm{u}=u_{x}\bm{\hat{x}}$,
$p_{\perp}=p_{\perp 0}$, $p_{\|}=p_{\|0}$, and $\bm{B}=(B_{x},B_{y},0)$, one
finds the set of eight eigenvalues (Baranov, 1970; Meng et al., 2012),

$\displaystyle k^{-1}\omega_{C}=u_{x},\qquad k^{-1}\omega_{\Delta p}=u_{x},$
(42) $\displaystyle k^{-1}\omega_{\rm
A}=u_{x}\pm\left(\frac{\Uptheta_{0}}{\rho_{0}}\right)^{1/2}B_{x},$ (43)
$\displaystyle k^{-1}\omega_{{\rm MS}\pm}=u_{x}$
$\displaystyle\pm\left(\frac{1}{2\rho_{0}}\right)^{1/2}\Bigg{\\{}2p_{\perp
0}+B^{2}+\hat{b}_{x}^{2}(2p_{\|0}-p_{\perp 0})+T_{e}\rho_{0}$
$\displaystyle\pm\bigg{[}\left(2p_{\perp
0}+B^{2}+\hat{b}_{x}^{2}(2p_{\|0}-p_{\perp
0})+T_{e}\rho_{0}\right)^{2}+4p_{\perp
0}^{2}\hat{b}_{x}^{2}(1-\hat{b}_{x}^{2})$ $\displaystyle\qquad\>\>-12p_{\perp
0}p_{\|0}\hat{b}_{x}^{2}(2-\hat{b}_{x}^{2})+12p_{\|0}^{2}\hat{b}_{x}^{4}-12p_{\|0}B^{2}\hat{b}_{x}^{2}\bigg{]}^{1/2}\Bigg{\\}}^{1/2},$
(44)

where $\Uptheta_{0}=1+(p_{\perp 0}-p_{\|0})/B^{2}$. Here $\omega_{\rm A}$ is
like the MHD shear-Alfvén wave, modified by the pressure anisotropy, while
$\omega_{C}$ and $\omega_{\Delta p}$ are two entropy-like waves. One of these
entropy waves ($\omega_{C}$) involves only density perturbations and no
$\Delta p$, like the MHD entropy mode, but is strongly damped by heat fluxes.
The other mode ($\omega_{\Delta p}$) involves balanced perturbations to
$p_{\perp}$, $p_{\parallel}$, $\rho$, and $B$ in general and becomes the
gyrokinetic non-propagating mode with the inclusion of heat fluxes (Howes et
al., 2006).

The CGL magnetosonic waves $\omega_{{\rm MS}+}$ and $\omega_{{\rm MS}-}$,
which correspond to the $+$ and $-$ on the second line of (44), respectively,
are compressive waves that also involve perturbations to the magnetic field.
In general, for $u_{x}=0$, $|\omega_{{\rm MS}+}|\geq|\omega_{{\rm MS}-}|$ and
$|\omega_{{\rm MS}+}|\geq|\omega_{\rm A}|$, but for $\beta\gtrsim 1$,
$|\omega_{{\rm MS}-}|>|\omega_{\rm A}|$, meaning that CGL ‘slow’ waves
propagate faster than shear-Alfvén waves, unlike in standard MHD. As discussed
in § 4.3.3, the properties of obliquely propagating CGL slow waves are
potentially of interest to explain the compressive features of the cascade. A
key characteristic is their perpendicular pressure balance, $\delta
p_{\perp}+\delta(B^{2})/8\pi\approx 0$ or $\delta
p_{\perp}\approx-(B_{0}^{2}/4\pi)\delta B_{\|}/B_{0}$ (just like the MHD slow
mode), which is derived by applying the standard RMHD ordering (Schekochihin
et al., 2009) to the perpendicular part of the momentum equation 2 (with
$T_{e}=0$). Then, using the linearised form of the CGL invariants
($d\mu/dt=d\mathcal{J}/dt=0$) yields $\delta
p_{\perp}/p_{0}-\delta\rho/\rho_{0}-\delta B_{\|}/B_{0}=0$ and $\delta
p_{\|}/p_{0}-3\delta\rho/\rho_{0}+2\delta B_{\|}/B_{0}=0$, which can be
combined with the $\delta p_{\perp}$ constraint to yield $\delta
p_{\|}=-(5+6/\beta)\delta B_{\|}/B_{0}$. Thus we see that $\delta
p_{\|}/\delta p_{\perp}=(5\beta+6)/2\gg 1$, a strong dominance of parallel
pressure, as observed in § 4.3 (although this prediction is much stronger
than that observed). Then, because $\delta p_{\|}$ provides the restoring
force for the wave and it is not constrained by $\delta p_{\perp}$ (unlike
MHD, where $p_{\perp}=p_{\|}$), the frequency scales with the thermal speed,
becoming $\omega_{{\rm MS}-}^{2}=(5/2)k_{\|}^{2}p_{0}/\rho_{0}$ for $\beta\gg
1$. In the presence of heat fluxes, this general structure is maintained,
although the wave becomes strongly damped at a rate that is approximately half
that of its propagation frequency. Weak collisions with $\nu_{\rm c}\lesssim
k_{\|}v_{\rm th}$ cause even stronger damping by weakly coupling together
$\delta p_{\perp}$ and $\delta p_{\|}$; the mode becomes non-propagating
around $\nu_{\rm c}\sim k_{\|}v_{\rm th}$ before reverting towards the MHD
slow mode for $\nu_{\rm c}\gg k_{\|}v_{\rm th}$ (see Majeski et al., 2023).

Full expressions for the eigenvectors for all angles are inconveniently
complex and not important for our purposes, so we do not list them here.

###  A.2 Linear-wave convergence

Figure 17: Propagation of oblique linear waves within the Athena++ CGL Landau-
fluid implementation in one spatial dimension. Background parameters are
$p_{\perp}=p_{\|}=5$, $\rho=1$, $\bm{u}_{0}=0$, and
$\bm{B}_{0}=(1,\sqrt{2},0.5)$, with $\bm{k}=2\pi\hat{\bm{x}}$ and
$k_{L}=|\bm{k}|$. The left panel shows the time evolution of $\rho
u_{y}(x_{m},t)$, where $x_{m}$ is chosen so that $\rho u_{y}(x_{m},0)$ is the
maximum of the sinusoidal initial condition. Results are normalised to the
initial amplitude of ${\simeq}10^{-4}$. The curves show the expected linear
solution, computed analytically from the dispersion relation (42)–(44), while
symbols show results from Athena++ with $N_{x}=1024$. Blue lines/symbols show
$\nu_{\rm c}=10$, with the solid line (plus symbols), dashed line (circle
symbols), dot-dashed line (star symbols), and dotted line (triangle symbols)
corresponding to the Alfvén, slow-like, fast-like, and entropy-like modes,
respectively. The red dashed curve (circle symbols) shows the slow mode at
$\nu_{\rm c}=10^{10}$ to illustrate how the system reverts to undamped
adiabatic MHD at high $\nu_{\rm c}$. The right panel shows the normalised
root-mean-squared error in the solution, as measured from the analytic
solution, as a function of spatial resolution $N_{x}$. Line styles and colours
are as in the left panel, with the addition of yellow curves for collisionless
modes ($\nu_{\rm c}=0$) and black for pure CGL (no heat fluxes;
$k_{L}\rightarrow\infty$). We see almost second-order convergence for most
solutions, except for those that are dominated by damping from heat fluxes, in
particular the collisionless/weakly collisional slow mode (blue and yellow
dashed lines; the former is overdamped $|{\rm Im}(\omega)|>|{\rm
Re}(\omega)|$, see left-hand panel).

Code testing was carried out with standard problems and methods, such as
examining the total-energy conservation, nonlinear Alfvén-wave dynamics
(Squire et al., 2016), and linear waves. Figure 17 illustrates the ability of
the solver to capture the diverse linear behaviour of the CGL-LF system. The
tests involve initialising the solver with a chosen linear eigenmode, computed
analytically from the linearised system, then comparing the frequency and
decay rate from the code to the analytic expectations. In the right-hand
panel, we show the convergence to the analytic solution with spatial
resolution for a range of different waves and parameters, scanning from the
collisionless system out to the adiabatic-MHD with $\nu_{\rm c}=10^{10}$. For
modes that are only weakly damped, meaning their dynamics are dominated by the
conservative (CGL) solver, we see almost second-order convergence. With more
strongly damped modes, for which the error becomes dominated by the non-
conservative heat flux and/or collisions in the solver, we see first-order
convergence above some resolution that depends on the mode in question (this
behaviour is clearest for the slow-like modes shown with dashed lines, because
these are more strongly damped).

###  A.3 Numerical problems and challenges

The CGL system seems to be particularly prone to numerical instability, which
hampered efforts to design more accurate numerical solvers. As part of the
code development effort for this work, we tested a wide variety of different
numerical options and solvers, all of which showed some issues in certain
situations. Our final choice of the piecewise parabolic method with an HLL
Reimman solver and the conserved variable $\mathcal{A}$ (see 39) was made by
trading off the various issues that arose in different numerical tests.
Examples include a pernicious numerical instability that acts on the
compressive components of the system at the smallest scales, creating flute-
like ($k_{\perp}\gg k_{\|}$) grid-scale modes that can grow sufficiently to
break up wave solutions in certain circumstances262626In particular, they are
strongest when they can directly align with the numerical grid, which can
occur when the magnetic field is nearly constant and grid aligned.; or
virulent numerical instabilities that arise in regions of small $|\bm{B}|$
when $\mu$ rather than $\mathcal{A}$ is used as the conserved variable. While
we do not fully understand the root causes of these issues, we speculate that
they may be worse than standard in MHD because the system involves two
different volumetrically conserved variables ($\rho$ and $\mu$ or
$\mathcal{A}$) that feed back on the momentum in different ways.

As part of this development process, we designed a new Riemann solver for the
CGL system that treats the Alfvénic-branch solutions separately, just like the
popular HLLD solver for MHD (Miyoshi & Kusano, 2005; Mignone, 2007). This has
significantly improved accuracy for problems involving predominantly Alfvénic
fluctuations. Although this made it an obvious candidate for our Alfvénic
turbulence simulations, the solver was ultimately found to cause issues in our
higher-resolution turbulence simulations, even though it worked quite well for
many simpler test problems. This type of behaviour is not uncommon, for
example, also occurring in HLLC and related hydrodynamic finite-volume
solvers, which can suffer from serious numerical issues if shocks develop
(see, e.g., Simon & Mandal, 2019, and references therein). For this reason,
for our main simulations, we reverted back to the simplest HLL solver, using
parabolic reconstruction to reduce the large numerical dissipation that is
inherent to this method. Nonetheless, in case the HLLD solver proves useful
for other problems in future works, we provide a brief description of it
below.

####  A.3.1 An HLLD CGL Riemann solver

Rather than proceeding in the most obvious way, by extending the HLLD method
of Miyoshi & Kusano (2005) to the double-adiabatic mode structure, we instead
base our solver on the isothermal HLLD solver of Mignone (2007). The method
thus computes the fluxes of thermodynamic variables using the simple HLL
method (Harten et al., 1983), while those of the transverse momentum and
magnetic field are obtained using a more accurate computation of the wave
structure. It enables shear-Alfvénic dynamics (or rotational discontinuities)
to be captured accurately, while avoiding various very complex and time-
consuming nonlinear solves that arise if one attempts to replicate the method
of Miyoshi & Kusano (2005)272727Indeed, following the logic of Miyoshi &
Kusano (2005) for the CGL system, one finds it necessary to solve for nearly
all inner states of the Riemann fan simultaneously in a complex system of
nonlinear equations. This leads to a slow and impractical solver..

The basis for the solver is to take $\rho$, $m_{x}$, $E$, $\mathcal{A}$, and
their fluxes as constant across the Riemann fan. In contrast, the ‘Alfvénic
variables’, $m_{y}$, $m_{z}$, $B_{y}$, and $B_{z}$, are allowed to vary across
the fan, being separated into three states $\bm{U}_{L}^{*}$, $\bm{U}_{c}^{*}$,
and $\bm{U}_{R}^{*}$ by the left and right Alfvén discontinuities with
velocities $S^{*}_{L}$ and $S^{*}_{R}$, as in isothermal MHD (Mignone, 2007).
We then use the consistency conditions for a five-state Riemann solver
separated by waves $S_{L}$, $S_{L}^{*}$, $S_{R}^{*}$, and $S_{R}$ (Mignone,
2007; Toro, 2009):

$\frac{(S_{L}^{*}-S_{L})\bm{U}^{*}_{L}+(S_{R}^{*}-S_{L}^{*})\bm{U}^{*}_{c}+(S_{R}-S_{R}^{*})\bm{U}^{*}_{R}}{S_{R}-S_{L}}=\bm{U}^{\rm
hll},$ (45)

for the conserved variables $\bm{U}$, and

$\frac{(S_{L}^{*-1}-S_{L}^{-1})\bm{F}^{*}_{L}+(S_{R}^{*-1}-S_{L}^{*-1})\bm{F}^{*}_{c}+(S_{R}^{-1}-S_{R}^{*-1})\bm{F}^{*}_{R}}{S_{R}^{-1}-S_{L}^{-1}}=\bm{F}^{\rm
hll}$ (46)

for the fluxes $\bm{F}$. Here $\bm{U}^{\rm hll}$ and $\bm{F}^{\rm hll}$ are
the HLLE states and fluxes:

$\displaystyle\bm{U}^{\rm
hll}\doteq\frac{S_{R}\bm{U}_{R}-S_{L}\bm{U}_{L}-\bm{F}_{R}+\bm{F}_{L}}{S_{R}-S_{L}},$
(47) $\displaystyle\bm{F}^{\rm
hll}\doteq\frac{S_{R}^{-1}\bm{F}_{R}-S_{L}^{-1}\bm{F}_{L}-\bm{U}_{R}+\bm{U}_{L}}{S_{R}^{-1}-S_{L}^{-1}}.$
(48)

Applying the assumption
$\bm{U}_{L}^{*}=\bm{U}_{c}^{*}=\bm{U}_{R}^{*}=\bm{U}^{*}$ to the non-Alfvénic
variables gives $\bm{U}^{*}=\bm{U}^{\rm hll}$, $\bm{F}^{*}=\bm{F}^{\rm hll}$,
as expected. It remains to determine $\bm{U}_{L}^{*}$, $\bm{U}_{R}^{*}$, and
their fluxes for the Alfvénic variables, after which $\bm{F}_{c}^{*}$ is
determined from 46.

We compute the jump conditions for the Alfvénic variables across the $S_{L}$
and $S_{R}$ waves following Mignone (2007):

$\displaystyle\bm{F}^{*}_{L}=\bm{F}_{L}+S_{L}(\bm{U}_{L}-\bm{U}_{L}^{*}),$
(49)
$\displaystyle\bm{F}^{*}_{R}=\bm{F}_{R}+S_{R}(\bm{U}_{R}-\bm{U}_{R}^{*}).$
(50)

This leads to

$\displaystyle\rho^{*}u_{y\alpha}^{*}=\rho^{*}u_{y\alpha}-B_{x}B_{y\alpha}\frac{u_{x}^{*}\Uptheta_{\alpha}-u_{x\alpha}\Uptheta^{*}+S_{\alpha}(\Uptheta^{*}-\Uptheta_{\alpha})}{(S_{\alpha}-S_{L}^{*})(S_{\alpha}-S_{R}^{*})},$
(51)
$\displaystyle\rho^{*}u_{z\alpha}^{*}=\rho^{*}u_{z\alpha}-B_{x}B_{z\alpha}\frac{u_{x}^{*}\Uptheta_{\alpha}-u_{x\alpha}\Uptheta^{*}+S_{\alpha}(\Uptheta^{*}-\Uptheta_{\alpha})}{(S_{\alpha}-S_{L}^{*})(S_{\alpha}-S_{R}^{*})},$
(52) $\displaystyle
B_{y\alpha}^{*}=\frac{B_{y\alpha}}{\rho^{*}}\frac{\rho_{\alpha}(S_{\alpha}-u_{\alpha})^{2}-B_{x}^{2}\Uptheta_{\alpha}}{(S_{\alpha}-S_{L}^{*})(S_{\alpha}-S_{R}^{*})},$
(53) $\displaystyle
B_{z\alpha}^{*}=\frac{B_{z\alpha}}{\rho^{*}}\frac{\rho_{\alpha}(S_{\alpha}-u_{\alpha})^{2}-B_{x}^{2}\Uptheta_{\alpha}}{(S_{\alpha}-S_{L}^{*})(S_{\alpha}-S_{R}^{*})}.$
(54)

Here $\alpha$ denotes either $L$ or $R$, $\rho^{*}=\rho^{\rm hll}$,
$u_{x}^{*}=F_{\rho}^{\rm hll}/\rho^{\rm hll}$, and $\Uptheta^{*}$ is the
anisotropy parameter ($\Uptheta\doteq 1+\Delta p/B^{2}$) across the fan (see
below). The Alfvénic wave speeds, $S_{L,R}^{*}$, are

$S_{L}^{*}=u_{x}^{*}-\sqrt{\frac{\Uptheta^{*}}{\rho^{*}}}|B_{x}|,\quad
S_{R}^{*}=u_{x}^{*}+\sqrt{\frac{\Uptheta^{*}}{\rho^{*}}}|B_{x}|.$ (55)

These are fixed by the jump conditions from the left and right regions to the
central regions after asserting that $\Uptheta^{*}$ should be constant (this
is in keeping with the assumption that only Alfvénic quantities change inside
the fan; see equations (25)–(28) of Mignone 2007). The numerical flux is then
computed as

$\bm{F}=\begin{cases}\bm{F}_{L},&0<S_{L},\\\
\bm{F}_{L}^{*},&S_{L}<0<S_{L}^{*},\\\
\bm{F}_{c}^{*},&S_{L}^{*}<0<S_{R}^{*},\\\
\bm{F}_{R}^{*},&S_{R}^{*}<0<S_{R},\\\ \bm{F}_{R},&S_{R}<0,\end{cases}$ (56)

where $\bm{F}_{L}^{*}$ and $\bm{F}_{R}^{*}$ are computed from 49 and 50, and
$\bm{F}_{c}^{*}$ from the consistency condition (46).

The firehose parameter $\Uptheta^{*}$ remains unspecified. Unfortunately, it
is not possible to have $\Delta p/B^{2}$, $E$, and $\mathcal{A}$ all constant
across the central region for general $B_{y}$ and $B_{z}$ that vary as in 53
and 54. We are thus forced to consider $B^{2}$ in $\Uptheta$ to be some
average across the fan. There appears to be no way to circumvent this
inconsistency, so long as we assume that non-Alfvénic variables are constant
across then fan, while Alfvénic components vary.282828It is worth noting that
the same inconsistency arises in the isothermal MHD solver of Mignone (2007),
although this is not discussed in his paper. Specifically, this solver assumes
that, within the fan, the total pressure $p_{T}=c_{\rm s,iso}^{2}\rho+B^{2}/2$
is constant (this is required so that $u_{x}$ is constant), while also taking
$\rho$ to be constant. These assumptions are clearly inconsistent if $B^{2}$
varies. A similar inconsistency arises in the adiabatic MHD solver of Miyoshi
& Kusano (2005) across the $S_{\alpha}^{*}$ waves. We are thus left with
choosing a suitable average for defining $\Uptheta^{*}$. We have experimented
with various possibilities for doing this, finding that the most robust method
seems to be to use the HLL average to compute $B^{2}$ (i.e., using $B_{y}^{\rm
hll}$ and $B_{z}^{\rm hll}$). This method has the advantage of being
straightforward and fast to execute, while also being conceptually consistent
with the idea that non-Alfvénic variables ($B^{2}$ in this case) should be
computed using HLL averages.

Special cases. There are two special cases that must be dealt with in the
solver. The first, which occurs if $S_{\alpha}^{*}\rightarrow S_{\alpha}$, is
handled identically to the Mignone (2007) solver by imposing a zero jump
across $m_{y}$, $m_{z}$, $B_{y}$, and $B_{z}$. As in isothermal MHD, there is
no issue with the limit $B_{x}\rightarrow 0$, which simply leads to
$S_{\alpha}^{*}=u^{*}$ and various simplifications of (51)–(54).

The other special case occurs as $\Uptheta$ approaches 0. For $\Uptheta<0$,
Alfvén waves are unstable, and it no longer makes sense to compute the wave
structure in the same way. We thus switch to a standard HLL solver for all
variables if $\Uptheta_{L}$, $\Uptheta^{*}$, or $\Uptheta_{R}$ is below some
threshold $\Uptheta_{\mathrm{thresh}}$. In practice,
$\Uptheta_{\mathrm{thresh}}\approx 0.1$ seems to give reasonable results
(recall that Alfvénic perturbations propagate very slowly for small $\Uptheta$
anyway). The switch to the HLL solver is carried out using the ‘anti-diffusion
control’ method described presently.

Anti-diffusion control. As mentioned above, the CGL HLLD method displayed
numerical instabilities, in particular small-scale oscillations in $\Uptheta$,
which became worse for more complex problems. These can often be mitigated by
the ‘anti-diffusion control’ method of Simon & Mandal (2019), which was
proposed to solve the ‘Carbuncle’ problem of the hydrodynamic HLLC solver. The
method involves continuously switching to the more diffusive HLL solver at the
advent of grid-scale oscillations. It is derived by rewriting 49 and 50 as

$\displaystyle\bm{F}^{*}_{L}=\bm{F}^{\rm hll}+S_{L}(\bm{U}^{\rm
hll}-\bm{U}_{L}^{*}),$ (57) $\displaystyle\bm{F}^{*}_{R}=\bm{F}^{\rm
hll}+S_{R}(\bm{U}^{\rm hll}-\bm{U}_{R}^{*}),$ (58)

then recognising the second term on each right-hand side as an ‘anti-
diffusive’ contribution to the flux (it corrects the HLL flux, making it less
diffusive). Multiplying this contribution by a factor $\omega_{\mathrm{AD}}$
satisfying $0\leq\omega_{\mathrm{AD}}\leq 1$, gives a method to interpolate
between the HLLD and HLL solvers. As a simple option, $\omega_{\mathrm{AD}}$
can be chosen as
$\omega_{\mathrm{AD}}=\exp(-\alpha_{\mathrm{AD}}\max|\Delta\Uptheta|)$, where
$\max|\Delta\Uptheta|$ is the maximum value of $|\Uptheta_{L}-\Uptheta_{R}|$
across the interface and its two neighbors perpendicular to the interface. The
free parameter $\alpha_{\mathrm{AD}}$ controls the strength of the anti-
diffusive effect, with $\alpha_{\mathrm{AD}}\approx 10$ providing a reasonable
trade off between accuracy and stability for a range of parameters and
problems. We also multiply $\omega_{\mathrm{AD}}$ by
$H(\Uptheta_{\mathrm{min}}-\Uptheta_{\mathrm{thresh}})$, where $H(x)$ is the
Heaviside step function and
$\Uptheta_{\mathrm{min}}=\min\\{\Uptheta_{L},\Uptheta^{*},\Uptheta_{R}\\}$, to
switch to the HLL solver as the system approaches the firehose instability
condition (see above). More information and stability analysis of the anti-
diffusion control method in the HLLC context can be found in Simon & Mandal
(2019).

## References

  * Arzamasskiy et al. (2019) Arzamasskiy, L., Kunz, M. W., Chandran, B. D. G. & Quataert, E. 2019 Hybrid-kinetic simulations of ion heating in Alfvénic turbulence. Astrophys. J. 879, 53.
  * Arzamasskiy et al. (2023) Arzamasskiy, L., Kunz, M. W., Squire, J., Quataert, E. & Schekochihin, A. A. 2023 Kinetic turbulence in collisionless high-$\beta$ plasmas. Phys. Rev. X 13, 021014.
  * Bale et al. (2009) Bale, S. D., Kasper, J. C., Howes, G. G., Quataert, E., Salem, C. & Sundkvist, D. 2009 Magnetic fluctuation power near proton temperature anisotropy instability thresholds in the solar wind. Phys. Rev. Lett. 103, 211101.
  * Baranov (1970) Baranov, V. B. 1970 Simple waves in a plasma with anisotropic pressure. Fluid Dynamics 5, 181.
  * Barnes (1966) Barnes, A. 1966 Collisionless damping of hydromagnetic waves. Phys. Fluids 9, 1483.
  * Barnes & Hollweg (1974) Barnes, A. & Hollweg, J. V. 1974 Large-amplitude hydromagnetic waves. J. Geophys. Res. 79, 2302.
  * Belcher & Davis (1971) Belcher, J. W. & Davis, Leverett, J. 1971 Large-amplitude Alfvén waves in the interplanetary medium, 2. J. Geophys. Res. 76, 3534.
  * Boldyrev (2006) Boldyrev, S. 2006 Spectrum of magnetohydrodynamic turbulence. Phys. Rev. Lett. 96, 115002.
  * Borovsky & Valdivia (2018) Borovsky, J. E. & Valdivia, J. A. 2018 The Earth’s magnetosphere: A systems science overview and assessment. Surveys in Geophys. 39, 817.
  * Bott et al. (2021) Bott, A. F. A., Arzamasskiy, L., Kunz, M. W., Quataert, E. & Squire, J. 2021 Adaptive critical balance and firehose instability in an expanding, turbulent, collisionless plasma. Astrophys. J. Lett. 922, L35.
  * Braginskii (1965) Braginskii, S. I. 1965 Transport processes in a plasma. Rev. Plasma Phys. 1, 205.
  * Brandenburg & Subramanian (2005) Brandenburg, A. & Subramanian, K. 2005 Astrophysical magnetic fields and nonlinear dynamo theory. Phys. Rep. 417 1\.
  * Chandran et al. (2010) Chandran, B. D. G., Li, B., Rogers, B. N., Quataert, E. & Germaschewski, K. 2010 Perpendicular Ion Heating by Low-frequency Alfvén-wave Turbulence in the Solar Wind. Astrophys. J. 720, 503.
  * Chandrasekhar et al. (1958) Chandrasekhar, S., Kaufman, A. N. & Watson, K. M. 1958 The stability of the pinch. Proc. R. Soc. A 245, 435–455.
  * Chen (2016) Chen, C. H. K. 2016 Recent progress in astrophysical plasma turbulence from solar wind observations. J. Plasma Phys. 82, 535820602\.
  * Chen et al. (2020) Chen, C. H. K., Bale, S. D., Bonnell, J. W., Borovikov, D., Bowen, T. A., Burgess, D., Case, A. W., Chandran, B. D. G., de Wit, T. D., Goetz, K., Harvey, P. R., Kasper, J. C., Klein, K. G., Korreck, K. E., Larson, D., Livi, R., MacDowall, R. J., Malaspina, D. M., Mallet, A., McManus, M. D., Moncuquet, M., Pulupa, M., Stevens, M. L. & Whittlesey, P. 2020 The evolution and role of solar wind turbulence in the inner heliosphere. Astrophys. J. Supp. 246, 53.
  * Chen et al. (2013) Chen, C. H. K., Bale, S. D., Salem, C. S. & Maruca, B. A. 2013 Residual energy spectrum of solar wind turbulence. Astrophys. J. 770, 125.
  * Chen et al. (2021) Chen, C. H. K., Chandran, B. D. G., Woodham, L. D., Jones, S. I., Perez, J. C., Bourouaine, S., Bowen, T. A., Klein, K. G., Moncuquet, M., Kasper, J. C. & Bale, S. D. 2021 The near-Sun streamer belt solar wind: turbulence and solar wind acceleration. Astron. Astro. 650, L3.
  * Chen et al. (2012) Chen, C. H. K., Mallet, A., Schekochihin, A. A., Horbury, T. S., Wicks, R. T. & Bale, S. D. 2012 Three-dimensional structure of solar wind turbulence. Astrophys. J. 758, 120.
  * Chen et al. (2016) Chen, C. H. K., Matteini, L., Schekochihin, A. A., Stevens, M. L., Salem, C. S., Maruca, B. A., Kunz, M. W. & Bale, S. D. 2016 Multi-species measurements of the firehose and mirror instability thresholds in the solar wind. Astrophys. J. 825, 1.
  * Chew et al. (1956) Chew, C. F., Goldberger, M. L. & Low, F. E. 1956 The Boltzmann equation and the one-fluid hydromagnetic equations in the absence of particle collisions. Proc. R. Soc. London A 236, 112.
  * Cho & Lazarian (2009) Cho, J. & Lazarian, A. 2009 Simulations of electron magnetohydrodynamic turbulence. Astrophys. J. 701, 236\.
  * Churazov et al. (2004) Churazov, E., Forman, W., Jones, C., Sunyaev, R. & Böhringer, H. 2004 XMM-Newton observations of the Perseus cluster - II. Evidence for gas motions in the core. Mon. Not. R. Astron. Soc. 347, 29.
  * Coburn et al. (2022) Coburn, J. T., Chen, C. H. K. & Squire, J. 2022 A measurement of the effective mean free path of solar wind protons. J. Plasma Phys. 88, 175880502.
  * Cohen et al. (2017) Cohen, R., Gerrard, A. J., Lanzerotti, L. J., Soto-Chavez, A. R., Kim, H. & Manweiler, J. W. 2017 Climatology of high-$\beta$ plasma measurements in Earth’s inner magnetosphere. J. Geophys. Res. Space Phys. 122, 711.
  * Cowie & McKee (1977) Cowie, L. L. & McKee, C. F. 1977 The evaporation of spherical clouds in a hot gas. I - Classical and saturated mass loss rates. Astrophys. J. 211, 135.
  * Cox (2005) Cox, D. P. 2005 The three-phase interstellar medium revisited. Ann. Rev. Astron. Astrophys. 43, 337.
  * Dong & Stone (2009) Dong, R. & Stone, J. M. 2009 Buoyant bubbles in intracluster gas: Effects of magnetic fields and anisotropic viscosity. Astrophys. J. 704, 1309.
  * Dunn et al. (2023) Dunn, C., Bowen, T., Mallet, A., Badman, S. & Bale, S. 2023 Effect of spherical polarization on the magnetic spectrum of the solar wind. arXiv e-prints p. arXiv:2305.09763.
  * EHT Collaboration et al. (2019) EHT Collaboration et al. 2019 First M87 event horizon telescope results. i. The shadow of the supermassive black hole. Astrophys. J. Lett. 875, L1.
  * Fabian (1994) Fabian, A. C. 1994 Cooling flows in clusters of galaxies. Ann. Rev. Astron. Astrophys. 32, 277.
  * Finelli et al. (2021) Finelli, F., Cerri, S. S., Califano, F., Pucci, F., Laveder, D., Lapenta, G. & Passot, T. 2021 Bridging hybrid- and full-kinetic models with Landau-fluid electrons. I. 2D magnetic reconnection. Astron. Astro. 653, A156.
  * Foote & Kulsrud (1979) Foote, E. A. & Kulsrud, R. M. 1979 Hydromagnetic waves in high beta plasmas. Astrophys. J. 233, 302.
  * Gary & Nishimura (2003) Gary, S. P. & Nishimura, K. 2003 Resonant electron firehose instability: Particle-in-cell simulations. Phys. Plasmas 10, 3571.
  * Goldreich & Sridhar (1995) Goldreich, P. & Sridhar, S. 1995 Toward a theory of interstellar turbulence. Strong Alfvénic turbulence. Astrophys. J. 438, 763.
  * Grete et al. (2017) Grete, P., O’Shea, B. W., Beckwith, K., Schmidt, W. & Christlieb, A. 2017 Energy transfer in compressible magnetohydrodynamic turbulence. Phys. Plasmas 24, 092311.
  * Gross & Krook (1956) Gross, E. P. & Krook, M. 1956 odel for collision processes in gases: Small-amplitude oscillations of charged two-component systems. Phys. Rev. 102, 593.
  * Hammett & Perkins (1990) Hammett, G. W. & Perkins, F. W. 1990 Fluid moment models for Landau damping with application to the ion-temperature-gradient instability. Phys. Rev. Lett. 64, 3019.
  * Harten et al. (1983) Harten, A., Lax, P. & Leer, B. 1983 On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review 25, 35.
  * Hasegawa (1969) Hasegawa, A. 1969 Drift mirror instability of the magnetosphere. Phys. Fluids 12, 2642.
  * Hellinger et al. (2011) Hellinger, P., Matteini, L., Štverák, Š., Trávníček, P. M. & Marsch, E. 2011 Heating and cooling of protons in the fast solar wind between 0.3 and 1 AU: Helios revisited. J. Geophys. Res. Space Phys. 116, A09105.
  * Hellinger et al. (2006) Hellinger, P., Trávníček, P., Kasper, J. C. & Lazarus, A. J. 2006 Solar wind proton temperature anisotropy: Linear theory and wind/swe observations. Geophysical Research Letters 33, L09101.
  * Hellinger & Trávníček (2008) Hellinger, P. & Trávníček, P. M. 2008 Oblique proton fire hose instability in the expanding solar wind: Hybrid simulations. J. Geophys. Res. Space Phys. 113, A10109.
  * Hollweg (1971) Hollweg, J. V. 1971 Nonlinear Landau damping of Alfvén waves. Phys. Rev. Lett. 27, 1349.
  * Hollweg (1974a) Hollweg, J. V. 1974a On electron heat conduction in the solar wind. J. Geophys. Res. 79 (25), 3845.
  * Hollweg (1974b) Hollweg, J. V. 1974b Transverse Alfvén waves in the solar wind: Arbitrary $k$, $v_{0}$, $B_{0}$, and $|\delta B|$. J. Geophys. Res. 79, 1539.
  * Howes (2010) Howes, G. G. 2010 A prescription for the turbulent heating of astrophysical plasmas. Mon. Not. R. Astron. Soc. 409, L104.
  * Howes et al. (2008) Howes, G. G., Cowley, S. C., Dorland, W., Hammett, G. W., Quataert, E. & Schekochihin, A. A. 2008 A model of turbulence in magnetized plasmas: Implications for the dissipation range in the solar wind. J. Geophys. Res. Space Phys. 113, A05103.
  * Howes et al. (2006) Howes, G. G., Cowley, S. C., Dorland, W., Hammett, G. W., Quataert, E. & Schekochihin, A. e. A. 2006 Astrophysical gyrokinetics: Basic equations and linear theory. Astrophys. J. 651, 590\.
  * Hunana et al. (2019) Hunana, P., Tenerani, A., Zank, G. P., Khomenko, E., Goldstein, M. L., Webb, G. M., Cally, P. S., Collados, M., Velli, M. & Adhikari, L. 2019 An introductory guide to fluid models with anisotropic temperatures. Part 1. CGL description and collisionless fluid hierarchy. J. Plasma Phys. 85, 205850602.
  * Johnston et al. (2022) Johnston, Z., Squire, J., Mallet, A. & Meyrand, R. 2022 On the properties of Alfvénic switchbacks in the expanding solar wind: Three-dimensional numerical simulations. Phys. Plasmas 29, 072902\.
  * Kasper et al. (2019) Kasper, J. C., Bale, S. D., Belcher, J. W., Berthomier, M., Case, A. W., Chandran, B. D. G., Curtis, D. W., Gallagher, D., Gary, S. P., Golub, L., Halekas, J. S. et al. 2019 Alfvénic velocity spikes and rotational flows in the near-Sun solar wind. Nature 576, 228.
  * Kasper et al. (2002) Kasper, J. C., Lazarus, A. J. & Gary, S. P. 2002 Wind/SWE observations of firehose constraint on solar wind proton temperature anisotropy. Geophys. Res. Lett. 29, 1839.
  * Kawazura et al. (2019) Kawazura, Y., Barnes, M. & Schekochihin, A. A. 2019 Thermal disequilibration of ions and electrons by collisionless plasma turbulence. Proc. Nat. Acad. Sci. 116, 771.
  * Kawazura et al. (2020) Kawazura, Y., Schekochihin, A. A., Barnes, M., TenBarge, J. M., Tong, Y., Klein, K. G. & Dorland, W. 2020 Ion versus Electron Heating in Compressively Driven Astrophysical Gyrokinetic Turbulence. Phys. Rev. X 10, 041050.
  * Kempski et al. (2019) Kempski, P., Quataert, E., Squire, J. & Kunz, M. W. 2019 Shearing-box simulations of MRI-driven turbulence in weakly collisional accretion discs. Mon. Not. R. Astron. Soc. 486, 4013\.
  * Klein & Spence (2021) Klein, K. & Spence, H. 2021 Helioswarm: Leveraging multi-point, multi-scale spacecraft observations to characterize turbulence. In EGU General Assembly Conference Abstracts, pp. EGU21–6812.
  * Kulsrud (1983) Kulsrud, R. M. 1983 MHD description of plasma. In Handbook of Plasma Physics (ed. R N Sagdeev & M N Rosenbluth). Princeton University Press.
  * Kunz et al. (2018) Kunz, M. W., Abel, I. G., Klein, K. G. & Schekochihin, A. A. 2018 Astrophysical gyrokinetics: turbulence in pressure-anisotropic plasmas at ion scales and beyond. J. Plasma Phys. 84, 715840201.
  * Kunz et al. (2022) Kunz, M. W., Jones, T. W. & Zhuravleva, I. 2022 Plasma Physics of the Intracluster Medium, pp. 1–42. Singapore: Springer.
  * Kunz et al. (2015) Kunz, M. W., Schekochihin, A. A., Chen, C. H. K., Abel, I. G. & Cowley, S. C. 2015 Inertial-range kinetic turbulence in pressure-anisotropic astrophysical plasmas. J. Plasma Phys. 81, 325810501.
  * Kunz et al. (2010) Kunz, M. W., Schekochihin, A. A., Cowley, S. C., Binney, J. J. & Sanders, J. S. 2010 A thermally stable heating mechanism for the intracluster medium: turbulence, magnetic fields and plasma instabilities. Mon. Not. R. Astron. Soc. 410, 2446.
  * Kunz et al. (2014) Kunz, M. W., Schekochihin, A. A. & Stone, J. M. 2014 Firehose and mirror instabilities in a collisionless shearing plasma. Phys. Rev. Lett. 112, 205003.
  * Kunz et al. (2020) Kunz, M. W., Squire, J., Schekochihin, A. A. & Quataert, E. 2020 Self-sustaining sound in collisionless, high-$\beta$ plasma. J. Plasma Phys. 86, 905860603.
  * Kunz et al. (2016) Kunz, M. W., Stone, J. M. & Quataert, E. 2016 Magnetorotational turbulence and dynamo in a collisionless plasma. Phys. Rev. Lett. 117, 235101.
  * Ley et al. (2023) Ley, F., Zweibel, E. G., Riquelme, M., Sironi, L., Miller, D. & Tran, A. 2023 A heating mechanism via magnetic pumping in the intracluster medium. Astrophys. J. 947, 89.
  * Li et al. (2020) Li, Y., Gendron-Marsolais, M.-L., Zhuravleva, I., Xu, S., Simionescu, A., Tremblay, G. R., Lochhaas, C., Bryan, G. L., Quataert, E., Murray, N. W., Boselli, A., Hlavacek-Larrondo, J., Zheng, Y., Fossati, M., Li, M., Emsellem, E., Sarzi, M., Arzamasskiy, L. & Vishniac, E. T. 2020 Direct detection of black hole-driven turbulence in the centers of galaxy clusters. Astrophys. J. Lett. 889, L1.
  * Lithwick & Goldreich (2001) Lithwick, Y. & Goldreich, P. 2001 Compressible magnetohydrodynamic turbulence in interstellar plasmas. Astrophys. J. 562, 279\.
  * Lovejoy & Schertzer (2012) Lovejoy, S. & Schertzer, D. 2012 Haar wavelets, fluctuations and structure functions: convenient choices for geophysics. Nonlinear Proc. Geophys. 19, 513.
  * Lynn et al. (2012) Lynn, J. W., Parrish, I. J., Quataert, E. & Chandran, B. D. G. 2012 Resonance broadening and heating of charged particles in magnetohydrodynamic turbulence. Astrophys. J. 758, 78.
  * Majeski et al. (2023) Majeski, S., Kunz, M. & Squire, J. 2023 Microphysically modified magnetosonic modes in collisionless, high-$\beta$ plasmas. J. Plasma Phys. 89, 905890303.
  * Mallet et al. (2015) Mallet, A., Schekochihin, A. A. & Chandran, B. D. G. 2015 Refined critical balance in strong Alfvenic turbulence. Mon. Not. R. Astron. Soc. 449, L77.
  * Mallet et al. (2021) Mallet, A., Squire, J., Chandran, B. D. G., Bowen, T. & Bale, S. D. 2021 Evolution of large-amplitude Alfvén waves and generation of switchbacks in the expanding solar wind. Astrophys. J. 918, 62\.
  * Marsch (2006) Marsch, E. 2006 Kinetic physics of the solar corona and solar wind. Living Rev. Solar Phys. 3, 1.
  * Melville et al. (2016) Melville, S., Schekochihin, A. A. & Kunz, M. W. 2016 Pressure-anisotropy-driven microturbulence and magnetic-field evolution in shearing, collisionless plasma. Mon. Not. R. Astron. Soc. 459, 2701.
  * Meng et al. (2012) Meng, X., Tóth, G., Sokolov, I. V. & Gombosi, T. I. 2012 Classical and semirelativistic magnetohydrodynamics with anisotropic ion pressure. J. Comp. Phys. 231, 3610.
  * Meyer et al. (2014) Meyer, C. D., Balsara, D. S. & Aslam, T. D. 2014 A stabilized Runge-Kutta-Legendre method for explicit super-time-stepping of parabolic and mixed equations. J. Comp. Phys. 257, 594.
  * Meyrand et al. (2019) Meyrand, R., Kanekar, A., Dorland, W. & Schekochihin, A. A. 2019 Fluidization of collisionless plasma turbulence. Proc. Nat. Acad. Sci. 116, 1185.
  * Meyrand et al. (2021) Meyrand, R., Squire, J., Schekochihin, A. A. & Dorland, W. 2021 On the violation of the zeroth law of turbulence in space plasmas. J. Plasma Phys. 87, 535870301.
  * Mignone (2007) Mignone, A. 2007 A simple and accurate Riemann solver for isothermal MHD. J. Comp. Phys. 225, 1427.
  * Mikhailovskii & Tsypin (1971) Mikhailovskii, A. B. & Tsypin, V. S. 1971 Transport equations and gradient instabilities in a high pressure collisional plasma. Plasma Phys. 13, 785.
  * Miyoshi & Kusano (2005) Miyoshi, T. & Kusano, K. 2005 A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. J. Comp. Phys. 208, 315\.
  * Müller & Grappin (2005) Müller, W. & Grappin, R. 2005 Spectral energy dynamics in magnetohydrodynamic turbulence. Phys. Rev. Lett. 95, 114502.
  * Nocera et al. (1986) Nocera, L., Priest, E. R. & Hollweg, J. V. 1986 Nonlinear development of phase-mixed Alfvén waves. Geophys. Astrophys. Fluid Dyn. 35, 111.
  * Pan & Johnsen (2017) Pan, S. & Johnsen, E. 2017 The role of bulk viscosity on the decay of compressible, homogeneous, isotropic turbulence. J. Fluid Mech. 833, 717.
  * Parashar et al. (2018) Parashar, T. N., Matthaeus, W. H. & Shay, M. A. 2018 Dependence of kinetic plasma turbulence on plasma $\beta$. Astrophys. J. Lett. 864, L21.
  * Parker (1958) Parker, E. N. 1958 Dynamical Instability in an Anisotropic Ionized Gas of Low Density. Phys. Rev. 109, 1874.
  * Passot et al. (2014) Passot, T., Henri, P., Laveder, D. & Sulem, P.-L. 2014 Fluid simulations of ion scale plasmas with weakly distorted magnetic fields. FLR-Landau fluid simulations. Europ. Phys. J. D 68, 207.
  * Passot et al. (2012) Passot, T., Sulem, P. L. & Hunana, P. 2012 Extending magnetohydrodynamics to the slow dynamics of collisionless plasmas. Phys. Plasmas 19, 082113.
  * Peterson & Fabian (2006) Peterson, J. R. & Fabian, A. C. 2006 X-ray spectroscopy of cooling clusters. Phys. Rep. 427, 1.
  * Quataert (2001) Quataert, E. 2001 Low radiative-efficiency accretion flows. In Probing the Physics of Active Galactic Nuclei (ed. Bradley M. Peterson, Richard W. Pogge & Ronald S. Polidan), Astronomical Society of the Pacific Conference Series, vol. 224, p. 71.
  * Quataert & Gruzinov (1999) Quataert, E. & Gruzinov, A. 1999 Turbulence and particle heating in advection-dominated accretion flows. Astrophys. J. 520, 248\.
  * Rincon et al. (2015) Rincon, F., Schekochihin, A. A. & Cowley, S. C. 2015 Non-linear mirror instability. Mon. Not. R. Astron. Soc. Lett. 447, L45.
  * Riquelme et al. (2015) Riquelme, M. A., Quataert, E. & Verscharen, D. 2015 Particle-in-cell simulations of continuously driven mirror and ion cyclotron instabilities in high beta astrophysical and heliospheric plasmas. Astrophys. J. 800, 27.
  * Rosenbluth (1956) Rosenbluth, M. N. 1956 The stability of the pinch. Los Alamos Sci. Lab. Rep. LA-2030.
  * Rosin et al. (2011) Rosin, M. S., Schekochihin, A. A., Rincon, F. & Cowley, S. C. 2011 A non-linear theory of the parallel firehose and gyrothermal instabilities in a weakly collisional plasma. Mon. Not. R. Astron. Soc. 413, 7\.
  * Roy et al. (2022) Roy, S., Bandyopadhyay, R., Yang, Y., Parashar, T. N., Matthaeus, W. H., Adhikari, S., Roytershteyn, V., Chasapis, A., Li, H., Gershman, D. J., Giles, B. L. & Burch, J. L. 2022 Turbulent energy transfer and proton-electron heating in collisionless plasmas. Astrophys. J. 941, 137.
  * Russell (2023) Russell, A. J. B. 2023 Nonlinear damping and field-aligned flows of propagating shear Alfvén waves with braginskii viscosity. arXiv e-prints p. arXiv:2303.11128.
  * Santos-Lima et al. (2014) Santos-Lima, R., de Gouveia Dal Pino, E. M., Kowal, G., Falceta-Gonçalves, D., Lazarian, A. & Nakwacki, M. S. 2014 Magnetic field amplification and evolution in turbulent collisionless magnetohydrodynamics: An application to the intracluster medium. Astrophys. J. 781, 84.
  * Schekochihin (2022) Schekochihin, A. A. 2022 MHD Turbulence: A Biased Review. J. Plasma Phys. 88, 155880501.
  * Schekochihin et al. (2009) Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G., Quataert, E. & Tatsuno, T. 2009 Astrophysical gyrokinetics: Kinetic and fluid turbulent cascades in magnetized weakly collisional plasmas. Astrophys. J. Supp. 182, 310.
  * Schekochihin et al. (2005) Schekochihin, A. A., Cowley, S. C., Kulsrud, R. M., Hammett, G. W. & Sharma, P. 2005 Plasma instabilities and magnetic field growth in clusters of galaxies. Astrophys. J. 629, 139.
  * Schekochihin et al. (2008) Schekochihin, A. A., Cowley, S. C., Kulsrud, R. M., Rosin, M. S. & Heinemann, T. 2008 Nonlinear growth of firehose and mirror fluctuations in astrophysical plasmas. Phys. Rev. Lett. 100, 081301.
  * Sharma & Hammett (2007) Sharma, P. & Hammett, G. W. 2007 Preserving monotonicity in anisotropic diffusion. J. Comp. Phys. 227, 123.
  * Sharma & Hammett (2011) Sharma, P. & Hammett, G. W. 2011 A fast semi-implicit method for anisotropic diffusion. J. Comp. Phys. 230, 4899.
  * Sharma et al. (2006) Sharma, P., Hammett, G. W., Quataert, E. & Stone, J. M. 2006 Shearing box simulations of the mri in a collisionless plasma. Astrophys. J. 637, 952.
  * Sharma et al. (2007) Sharma, P., Quataert, E., Hammett, G. W. & Stone, J. M. 2007 Electron heating in hot accretion flows. Astrophys. J. 667, 714\.
  * Simon & Mandal (2019) Simon, S. & Mandal, J. C. 2019 A simple cure for numerical shock instability in the HLLC Riemann solver. J. Comp. Phys. 378, 477\.
  * Sironi & Narayan (2015) Sironi, L. & Narayan, R. 2015 Electron heating by the ion cyclotron instability in collisionless accretion flows. I. Compression-driven instabilities and the electron heating mechanism. Astrophys. J. 800, 88.
  * Snyder et al. (1997) Snyder, P. B., Hammett, G. W. & Dorland, W. 1997 Landau fluid models of collisionless magnetohydrodynamics. Phys. Plasmas 4, 3974\.
  * Squire et al. (2020) Squire, J., Chandran, B. D. G. & Meyrand, R. 2020 In-situ switchback formation in the expanding solar wind. Astrophys. J. Lett. 891, L2.
  * Squire et al. (2017b) Squire, J., Kunz, M. W., Quataert, E. & Schekochihin, A. A. 2017 Kinetic simulations of the interruption of large-amplitude shear-Alfvén waves in a high-$\beta$ plasma. Phys. Rev. Lett. 119, 155101\.
  * Squire et al. (2016) Squire, J., Quataert, E. & Schekochihin, A. A. 2016 A stringent limit on the amplitude of Alfvénic perturbations in high-beta low-collisionality plasmas. Astrophys. J. Lett. 830, L25.
  * Squire et al. (2017a) Squire, J., Schekochihin, A. A. & Quataert, E. 2017 Amplitude limits and nonlinear damping of shear-Alfvén waves in high-beta low-collisionality plasmas. New J. Phys. 19, 055005.
  * Squire et al. (2019) Squire, J., Schekochihin, A. A., Quataert, E. & Kunz, M. W. 2019 Magneto-immutable turbulence in weakly collisional plasmas. J. Plasma Phys. 85, 905850114.
  * St-Onge et al. (2020) St-Onge, D. A., Kunz, M. W., Squire, J. & Schekochihin, A. A. 2020 Fluctuation dynamo in a weakly collisional plasma. J. Plasma Phys. 86, 905860503.
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J. & Felker, K. G. 2020 The Athena++ Adaptive Mesh Refinement Framework: Design and Magnetohydrodynamic Solvers. Astrophys. J. Supp. 249, 4.
  * Tenerani & Velli (2020) Tenerani, A. & Velli, M. 2020 Alfvénic fluctuations in the solar wind: nonlinearities and pressure anisotropy effects. Plasma Phys. Control. Fusion 62, 014001.
  * Toro (2009) Toro, E. 2009 Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Springer Berlin Heidelberg.
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S. & Werk, J. K. 2017 The circumgalactic medium. Ann. Rev. Astron. Astrophys. 55, 389\.
  * Vaidya et al. (2017) Vaidya, B., Prasad, D., Mignone, A., Sharma, P. & Rickler, L. 2017 Scalable explicit implementation of anisotropic diffusion with Runge-Kutta-Legendre super-time stepping. Mon. Not. R. Astron. Soc. 472, 3147.
  * Vasquez & Hollweg (1998) Vasquez, B. J. & Hollweg, J. V. 1998 Formation of spherically polarized Alfvén waves and imbedded rotational discontinuities from a small number of entirely oblique waves. J. Geophys. Res. Space Phys. 103, 335.
  * Vedenov & Sagdeev (1958) Vedenov, A. A. & Sagdeev, R. Z. 1958 Sov. Phys. Dokl. 3, 278\.
  * Verscharen et al. (2022) Verscharen, D., Chandran, B. D. G., Boella, E., Halekas, J., Innocenti, M. E., Jagarlamudi, V. K., Micera, A., Pierrard, V., Štverák, Š., Vasko, I. Y., Velli, M. & Whittlesey, P. L. 2022 Electron-driven instabilities in the solar wind. Front. Astron. Space Sci. 9, 951628.
  * Verscharen et al. (2016) Verscharen, D., Chandran, B. D. G., Klein, K. G. & Quataert, E. 2016 Collisionless isotropization of the solar-wind protons by compressive fluctuations and plasma instabilities. Astrophys. J. 831, 128\.
  * Verscharen et al. (2019) Verscharen, D., Klein, K. G. & Maruca, B. A. 2019 The multi-scale nature of the solar wind. Living Rev. Solar Phys. 16, 5.
  * Völk & Aplers (1973) Völk, H. J. & Aplers, W. 1973 The propagation of Alfvén waves and their directional anisotropy in the solar wind. Astrophys. Space Sci. 20, 267.
  * White et al. (2016) White, C. J., Stone, J. M. & Gammie, C. F. 2016 An extension of the Athena++ code framework for GRMHD based on advanced riemann solvers and staggered-mesh constrained transport. Astrophys. J. Supp. 225, 22.
  * Yang et al. (2017) Yang, Y., Matthaeus, W. H., Parashar, T. N., Haggerty, C. C., Roytershteyn, V., Daughton, W., Wan, M., Shi, Y. & Chen, S. 2017 Energy transfer, pressure tensor, and heating of kinetic plasma. Phys. Plasmas 24, 072306.
  * Yang et al. (2022) Yang, Y., Matthaeus, W. H., Roy, S., Roytershteyn, V., Parashar, T. N., Bandyopadhyay, R. & Wan, M. 2022 Pressure-strain interaction as the energy dissipation estimate in collisionless plasma. Astrophys. J. 929, 142.
  * Yoon et al. (1993) Yoon, P. H., Wu, C. S. & de Assis, A. S. 1993 Effect of finite ion gyroradius on the fire-hose instability in a high beta plasma. Phys. Fluids B 5, 1971.
  * Zhdankin et al. (2021) Zhdankin, V., Uzdensky, D. A. & Kunz, M. W. 2021 Production and persistence of extreme two-temperature plasmas in radiative relativistic turbulence. Astrophys. J. 908, 71.
  * Zhuravleva et al. (2019) Zhuravleva, I., Churazov, E., Schekochihin, A. A., Allen, S. W., Vikhlinin, A. & Werner, N. 2019 Suppressed effective viscosity in the bulk intergalactic plasma. Nature Astron. 3, 832.
  * Zhuravleva et al. (2014) Zhuravleva, I., Churazov, E., Schekochihin, A. A. et al. 2014 Turbulent heating in galaxy clusters brightest in X-rays. Nature 515, 85.