Pressure anisotropy and viscous heating in weakly collisional plasma turbulence

Abstract

Pressure anisotropy can strongly influence the dynamics of weakly collisional, high-beta plasmas, but its effects are missed by standard magnetohydrodynamics (MHD). Small changes to the magnetic-field strength generate large pressure-anisotropy forces, heating the plasma, driving instabilities and rearranging flows, even on scales far above the particles’ gyroscales where kinetic effects are traditionally considered most important. Here, we study the influence of pressure anisotropy on turbulent plasmas threaded by a mean magnetic field (Alfvénic turbulence). Extending previous results that were concerned with Braginskii MHD, we consider a wide range of regimes and parameters using a simplified fluid model based on drift kinetics with heat fluxes calculated using a Landau-fluid closure. We show that viscous (pressure-anisotropy) heating dissipates between a quarter (in collisionless regimes) and half (in collisional regimes) of the turbulent-cascade power injected at large scales; this does not depend strongly on either plasma beta or the ion-to-electron temperature ratio. This will in turn influence the plasma's thermodynamics by regulating energy partition between different dissipation channels (e.g. electron and ion heat). Due to the pressure anisotropy's rapid dynamic feedback onto the flows that create it – an effect we term ‘magneto-immutability’ – the viscous heating is confined to a narrow range of scales near the forcing scale, supporting a nearly conservative, MHD-like inertial-range cascade, via which the rest of the energy is transferred to small scales. Despite the simplified model, our results – including the viscous heating rate, distributions and turbulent spectra – compare favourably with recent hybrid-kinetic simulations. This is promising for the more general use of extended-fluid (or even MHD) approaches to model weakly collisional plasmas such as the intracluster medium, hot accretion flows and the solar wind.

1. Introduction

Many hot, diffuse plasmas in astrophysical environments are weakly collisional, with Coulomb mean free paths that are comparable to relevant macroscopic scales in the system. Canonical examples include the intracluster medium (ICM) in galaxy clusters (e.g. Fabian 1994; Peterson & Fabian 2006; Kunz, Jones & Zhuravleva 2022), hot accretion flows such as those observed by the Event Horizon Telescope (Quataert 2001; EHT Collaboration et al. 2019), the hotter phases of the interstellar and circumgalactic media (Cox 2005; Tumlinson, Peeples & Werk 2017) and the solar wind and magnetosphere (Marsch 2006; Borovsky & Valdivia 2018). In all of these examples, the plasma is extremely well magnetised, in that the system scales are much larger than the ion gyroradius $\rho _{i}$. Coupled with the plasmas’ large Coulomb mean free paths, this implies that magnetic fields are crucial in providing the plasma with an internal cohesion that allows it to behave more or less like a collisional fluid (Kulsrud 1983). However, despite the particle motion's allegiance to the local magnetic field, in many environments the magnetic fields are energetically subdominant as measured by $\beta = 8{\rm \pi} p/B^{2}$, where $p$ is a thermal pressure and $B^{2}/8{\rm \pi}$ is the magnetic energy density. Under such conditions, small relative changes to the magnetic field can easily cause particle distributions to become unstable, suggesting that non-equilibrium kinetic physics should play a key role in the dynamics.

Nonetheless, such conditions are commonly modelled using collisional magneto- hydrodynamics (MHD), usually for reasons of expediency, although with rigorous justification in some circumstances (e.g. Kulsrud 1983; Schekochihin et al. 2009; Kunz et al. 2015). It is the first purpose of this article to explain why this is usually not appropriate, as a result of the aforementioned kinetic physics; it is the second purpose to explain why, in the end, it is not so bad after all. The physics we explore is that of ‘pressure anisotropy’ (equivalently, temperature anisotropy), which is the difference between the thermal pressures in the directions perpendicular and parallel to the local magnetic field. The relevance of pressure anisotropy stems from it being the only non-isotropic piece of the pressure tensor that can survive on scales arbitrarily far above the particles’ gyroradii (Braginskii 1965), making it a key physical ingredient in the momentum balance of weakly collisional plasmas (compared with collisional MHD; Chew, Goldberger & Low 1956; Kulsrud 1983). At $\beta \gg 1$, a very small relative pressure anisotropy can lead to large bulk forces on the plasma. Furthermore, whenever the magnetic-field strength $B$ changes in time, $\Delta p$ is generated as a result of the conservation of single-particle adiabatic invariants. Put together, these two properties suggest that pressure-anisotropy stresses can vastly overwhelm the forces from the magnetic field or Reynolds stress in many high-$\beta$ environments.

Our study focuses on the role of pressure anisotropy in turbulence, specifically, in turbulence of the ‘Alfvénic’ variety that occurs when the system is threaded by a large-scale mean magnetic field. A particularly important example of such a system is the solar wind, which is our best natural laboratory for the study of Alfvénic turbulence under weakly collisional conditions. Such turbulence may also be relevant quasi-universally at smaller scales in systems without a mean magnetic field (Brandenburg & Subramanian 2005; Schekochihin 2022). We explore the structure and dynamics of weakly collisional Alfvénic turbulence, with a particular focus on their comparison with MHD and the implications for turbulent plasma heating. Pressure anisotropy provides another channel – effectively a viscosity – for damping injected mechanical energy into heat near the outer scales (e.g. Kunz et al. 2010; Yang et al. 2017). (Although this energy transfer can be reversible in some regimes, we will often refer to it as ‘viscous heating’ because at higher collisionalities the pressure-anisotropy stress increasingly resembles a parallel viscosity, with an associated heating rate that is positive definitive.) Because the pressure-anisotropy stress can easily overwhelm the magnetic tension and Reynolds stress in the high-$\beta$ limit (Squire, Quataert & Schekochihin 2016), simple estimates suggest that this viscous heating should completely dominate over other heating mechanisms. However, we show that, because of the dynamic feedback of the pressure anisotropy on the plasma flow, such heating is confined to a narrow range of scales near the forcing scales, with a sizeable fraction (typically, the majority) of the energy making it through to participate in a nearly conservative turbulent cascade. This phenomenon, which was termed ‘magneto-immutability’ by Squire et al. (2019) (hereafter S$+$19) who studied the effect in the context of the simplified ‘Braginskii MHD’ model (see also Kempski et al. 2019; St-Onge et al. 2020), involves the force from the pressure anisotropy acting to rearrange the turbulence in order to reduce the influence of itself. This rearrangement constrains the global variation of $B$, hence the term ‘magneto-immutability’ to describe the effect. Its physics is analogous to that of incompressibility, in which the pressure force that results from a fluid compression rapidly drives a flow that counters the compression, thereby eliminating such motions and minimising variations in the density. This explains our flippant declaration above about the article's purpose: simple estimates suggest pressure anisotropy should dominate the force balance and heating, modifying the dynamics compared with MHD; however, the effect of this modification is to minimise the pressure anisotropy's own influence, thus allowing a turbulent dynamics that looks broadly similar to MHD below the outer scales (albeit with a few important differences).

These ideas complicate the already complex story of turbulent heating in weakly collisional plasmas (e.g. Quataert & Gruzinov 1999; Howes 2010; Kunz et al. 2010; Kawazura et al. 2020; Meyrand et al. 2021; Yang et al. 2022; Arzamasskiy et al. 2023). A useful concept is the ‘cascade efficiency’: the fraction of energy available to cascade down to small scales and heat the plasma via different collisionless mechanisms (Schekochihin et al. 2009; Chandran et al. 2010; Arzamasskiy et al. 2019). In our simulations, between ${\simeq }20\,\%$ and ${\simeq }45\,\%$ of the energy is lost through pressure-anisotropy heating at large scales (giving a cascade efficiency of ${\gtrsim }55\,\%$ for our choice of large-scale forcing), independently of $\beta$ and the electron-to-ion temperature ratio. Magneto-immutability ensures that the remainder of the energy cascades nearly conservatively, presumably eventually heating at small scales in the way predicted by gyrokinetics (Howes et al. 2008; Schekochihin et al. 2009; Kunz et al. 2018). With different heating processes implying a different partition of turbulent energy between electron and ion heat, or between perpendicular and parallel heat, these apparently esoteric details of the turbulent structure can strongly influence the basic thermodynamics of the plasma. For instance, while it is well known that small-scale collisionless processes in high-$\beta$ plasmas generically heat ions more than electrons (Howes 2010; Parashar, Matthaeus & Shay 2018; Kawazura et al. 2020; Roy et al. 2022), large-scale pressure anisotropy can also heat collisionless electrons (Sharma et al. 2007), meaning that the cascade efficiency – and thus magneto-immutability – could directly control the ion-to-electron heating ratio. While we do not explicitly study electron physics, our study provides a useful foundation for distilling and quantifying the important physics, particularly through the idea that pressure-anisotropy heating is confined to a small range of large scales.

Our approach to the study of pressure anisotropy here is computational and simplified. We use the so-called Chew-Goldberger-Low-Landau-fluid (hereafter CGL-LF) model (Chew et al. 1956; Snyder, Hammett & Dorland 1997), which attempts to approximate collisionless heat fluxes in order to obtain a closed, simplified fluid model for the plasma dynamics on scales well above $\rho _{i}$. This model, which is implemented using new numerical methods into the Athena++ framework (White, Stone & Gammie 2016; Stone et al. 2020), is then used to explore a range of conditions as the key parameters of the problem are varied. Detailed comparisons with otherwise equivalent MHD simulations are used to diagnose the influence of the pressure anisotropy on the turbulence and heating. A downside of the simplified-fluid-model approach is that various ad hoc, loosely justified approximations are necessary. These approximations include those made to derive the fluid closure and its numerical realisation, and perhaps more importantly, the methods used to cope with the fast-growing $\rho _{i}$-scale firehose and mirror instabilities that will inevitably arise in real systems (Schekochihin et al. 2005; Bale et al. 2009; Kunz, Schekochihin & Stone 2014), but which cannot be captured correctly by drift kinetics (e.g. Rosin et al. 2011; Rincon, Schekochihin & Cowley 2015). The upside of the approach, by contrast, is that the model effectively assumes infinite scale separation between the outer scale and $\rho _{i}$, which is clearly not possible in kinetic models that explicitly resolve $\rho _{i}$, but is the appropriate limit for most astrophysical systems. With this in mind, a useful subsidiary purpose of the article is to assess the successes of this simplified model in reproducing the physics seen in the hybrid-kinetic simulations of Arzamasskiy et al. (2023) (hereafter A$+$22), which explicitly resolve $\rho _{i}$ and the associated chaos of kinetic-scale instabilities at the expense of a much reduced inertial range. We find relatively good agreement in general, despite the ad hoc nature of the fluid model and the interpretative difficulties associated with limited scale separation in the hybrid-kinetic simulations. In particular, we find similar cascade efficiencies (with a similar forcing scheme) and similar pressure-anisotropy distributions compared with A$+$22, with some caveats. Accordingly, the results of this work are quite promising for the use of CGL-LF approaches in modelling weakly collisional plasmas.

1.1. Outline

Although most of the basic theoretical concepts presented here have appeared in previous literature, we feel it is helpful to keep the majority of the discussion and notation self-contained to clarify the approximations and key concepts. Thus, § 2 presents an overview of the physics of pressure anisotropy and magneto-immutability, starting from the drift-kinetic model of Kulsrud (1983). Lacking any quantitative theory of the important effects, we focus on qualitative explanations for the behaviour and effect of the pressure anisotropy and heat fluxes in different regimes. We also define the ‘interruption number’ $\mathcal {I}$, which quantifies the expected strength of pressure-anisotropic effects in Alfvénic turbulence, analogously to the Reynolds number or Mach number in fluids. In § 3 we then outline the numerical methods and diagnostics, before presenting detailed simulation results in § 4. Leveraging the versatility of the simplified fluid model, a focus of the results is the direct comparison with ‘passive-$\Delta p$’ simulations. These are identical to the standard simulations but with the anisotropic pressure force artificially removed, thereby affording a direct comparison with a counterfactual situation where magneto-immutability does not exist. The comparison clearly demonstrates both the similarities and differences between our pressure-anisotropic Landau-fluid model and MHD in spectra, distributions and scale-dependent heating functions. We follow with discussions of the key uncertainties related to kinetic instabilities (§ 5.1) and of the distinction between magneto-immutability and the so-called ‘spherically polarised’ states measured in the solar wind (a complicating factor for interpreting spacecraft observations). We also summarise most salient differences between magneto-immutable and MHD turbulence (§ 5.2). A full list of the simulations with their important parameters is given in table 1. We conclude in § 6. An appendix presents various technical results related to the numerical finite-volume implementation in the Athena++ code (White et al. 2016; Stone et al. 2020), including linear-wave convergence tests and a new Riemann solver for the CGL system (albeit one that we ultimately did not use in the simulations).

Table 1. A list of all simulations used in this article. Key input parameters are $\beta _{0} = 8{\rm \pi} p_{0}/B_{0}^{2}$ and $L_{\perp }\nu _{c}/v_{A}$, where the subscript ‘$0$’ refers to an initial value ($\beta$ decreases as $B$ grows and the turbulence causes ion heating). The interruption number $\mathcal {I}$ is computed from these initial parameters using (2.19). The collisionality regime is specified via $\ell _{\perp }^{\rm WC}$ and $\ell _{\perp }^{\rm CL}$, which are the approximate scales below which motions transition into the weakly collisional regime and collisionless regime, respectively. These are computed by equating $\nu _{c}/\omega _{A}(\ell _{\perp })=\beta ^{1/2}$ and $\nu _{c}/\omega _{A}(\ell _{\perp })=1$, where $\omega _{A}$ is estimated by $(2{\rm \pi} v_{A}/L_{\|})(\ell _{\perp }/L_{\perp })^{-2/3}$, as for a non-aligned critically balanced cascade (Goldreich & Sridhar 1995; Schekochihin 2022). A turbulent eddy of scale $\ell _{\perp }$ is in the Braginskii-MHD regime for $\ell _{\perp }\gtrsim \ell _{\perp }^{\rm WC}$, in the weakly collisional regime for $\ell _{\perp }^{\rm WC} \gtrsim \ell _{\perp }\gtrsim \ell _{\perp }^{\rm CL}$ and in the collisionless regime for $\ell _{\perp }\lesssim \ell _{\perp }^{\rm CL}$ (see § 2.1.3). The ‘Passive?’ column indicates whether an otherwise identical passive-$\Delta p$ simulation was run for comparison purposes. All simulations have $L_{\|}=2 L_{\perp }$ and the energy-injection rate $\varepsilon = 0.16 v_{A}^{3}/L_{\perp }$.

2. Theoretical and phenomenological background

2.1. The evolution and effect of pressure anisotropy

In this section, we introduce the basic concepts necessary to understand the numerical results presented in § 4. The goal is to provide some intuitive understanding of the effects of pressure anisotropy in different plasma regimes, starting from the basic equations for a collisionless plasma on large scales. We begin by assuming that the plasma pressure tensor is gyrotropic but anisotropic – it is invariant under rotations about the magnetic field, but can differ in the directions perpendicular and parallel to the field. This is justified for ‘MHD-range’ scales, viz. those larger than the ion gyroscale in space and slower than the ion gyrofrequency in time. We also assume a single ion species and isothermal electrons, which allows us to drop the anisotropic component of the electron pressure and obtain single-fluid equations for the dynamics of the ions. Although this can be formally justified in the moderately collisional limit using the electron's larger collision frequency and fast thermal speed (Rosin et al. 2011), here the choice is primarily one of simplicity, since we wish to focus on the dynamical effects of pressure anisotropy on fluid motions. Most of our simulations in this work will assume cold electrons anyway, in order to better understand and diagnose the basic processes at play.

With these assumptions, the first three moments of the ion distribution function satisfy (Chew et al. 1956; Kulsrud 1983)

(2.1)\begin{gather} \partial_{t}\rho +\boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u} ) = 0, \end{gather}

(2.2)\begin{gather}\rho \left(\partial_{t}\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} \right)={-} \boldsymbol{\nabla}\left(T_{e}\rho+ p_{{\perp}} + \frac{B^{2}}{8{\rm \pi}}\right)+ \boldsymbol{\nabla} \boldsymbol{\cdot} \left[\hat{\boldsymbol{b}} \hat{\boldsymbol{b}}\left( \Delta p + \frac{B^{2}}{4{\rm \pi}}\right)\right], \end{gather}

(2.3)\begin{gather}\partial_{t} \boldsymbol{B} = \boldsymbol{\nabla} \times (\boldsymbol{u} \times \boldsymbol{B} ) , \end{gather}

(2.4)\begin{gather}\partial_{t}p_{{\perp}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (p_{{\perp}} \boldsymbol{u}) + p_{{\perp}} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{\cdot} (q_{{\perp}}\hat{\boldsymbol{b}}) +q_{{\perp}} \boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol{b}} = p_{{\perp}} \hat{\boldsymbol{b}}\hat{\boldsymbol{b}}:\boldsymbol{\nabla} \boldsymbol{u} - \tfrac{1}{3}\nu_{c}\Delta p, \end{gather}

(2.5)\begin{gather}\partial_{t}p_{{\parallel}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (p_{{\parallel}} \boldsymbol{u}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (q_{{\parallel}}\hat{\boldsymbol{b}}) -2q_{{\perp}} \boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol{b}} ={-}2 p_{{\parallel}} \hat{\boldsymbol{b}}\hat{\boldsymbol{b}}:\boldsymbol{\nabla} \boldsymbol{u} +\tfrac{2}{3}\nu_{c}\Delta p. \end{gather}

We use Gaussian-CGS units, $\boldsymbol {B}$ and $\boldsymbol {u}$ are the magnetic field and ion flow velocity (also approximately equal to the electron flow velocity for scales well above the ion gyroscale), $B\doteq | \boldsymbol {B} |$ and $\hat {\boldsymbol {b}}=\boldsymbol {B}/B$ denote the magnetic-field strength and direction, $\rho$ is the mass density, $\nu _{c}$ is the ion–ion collision frequency,¹ $p_{\perp }$ and $p_{\parallel }$ are the components of the ion-pressure tensor perpendicular and parallel to the magnetic field and $q_{\perp }$ and $q_{\parallel }$ are fluxes of perpendicular and parallel ion heat in the direction parallel to the magnetic field. The electron temperature $T_{e}$ is constant in time and space, because electrons are assumed isothermal (their pressure is $p_{e}=T_{e}\rho /m_{i}$, related to the ion density by quasi-neutrality assuming a single-ion-species plasma, and we neglect ion–electron collisions). Equations (2.1)–(2.5) are solved numerically in the conservative form given in Appendix A. Importantly, this avoids explicitly computing the parallel rate of strain $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$, which can introduce serious numerical errors in some situations (Sharma & Hammett (2011); S$+$19). We also define the Alfvén speed $v_{A}=B/\sqrt {4{\rm \pi} \rho }$, the total pressure $p_{0}\doteq 2p_{\perp }/3+p_{\|}/3$, the parallel sound speed $c_{s\|}\doteq \sqrt {p_{\|}/\rho }$, the plasma beta $\beta \doteq 8{\rm \pi} p_{0}/B^{2}$ (similarly $\beta _{\|}\doteq 8{\rm \pi} p_{\|}/B^{2}$), the pressure anisotropy $\Delta p\doteq p_{\perp }-p_{\|}$, the normalised pressure anisotropy $\varDelta \doteq \Delta p/p_{0}$ and the ‘anisotropy parameter’

(2.6)\begin{equation} \varTheta \doteq 1+\frac{4{\rm \pi} \Delta p}{B^{2}}, \end{equation}

which measures the relative change in the propagation speed of linear Alfvén waves (viz. $v_{A}$) due to the contribution from the pressure-anisotropy stress (see § 2.1.2).

In order to avoid solving a kinetic equation for the heat fluxes, we close (2.1)–(2.5) with the expressions

(2.7)\begin{gather} q_{{\perp}} ={-}\frac{2 c_{s\|}^{2}}{\sqrt{2 {\rm \pi}}c_{s\|} |k_{{\parallel}}|+\nu_{c}} \left[\rho\boldsymbol{\nabla}_{{\parallel}} \left(\frac{p_{{\perp}}}{\rho}\right) - p_{{\perp}}\left(1-\frac{p_{{\perp}}}{p_{{\parallel}}} \right)\frac{\boldsymbol{\nabla}_{{\parallel}} B}{B} \right],\end{gather}

(2.8)\begin{gather}q_{{\parallel}} ={-} \frac{8 c_{s\|}^{2}}{\sqrt{8 {\rm \pi}}c_{s\|} |k_{{\parallel}}|+(3{\rm \pi}-8)\nu_{c}} \rho \boldsymbol{\nabla}_{{\parallel}} \left(\frac{p_{{\parallel}}}{\rho}\right). \end{gather}

These forms, which are known as a ‘Landau-fluid’ closure, were devised by Snyder et al. (1997), in order to match the linear behaviour of the true kinetic system as closely as possible (Hammett & Perkins (1990); see also Passot, Sulem & Hunana 2012). The $k_{\|}$ and $\boldsymbol {\nabla }_{\|}=\hat {\boldsymbol {b}} \boldsymbol {\cdot }\boldsymbol {\nabla }$ in (2.7)–(2.8) are the parallel wavenumber and gradient operator, respectively, with the non-local gradient-like operator $\boldsymbol {\nabla }_{\|}/|k_{\|}|$ arising as a result of approximating the effects of collisionless damping within the fluid model (see § 3.1.1 for further discussion). We will use the form (2.7)–(2.8) both computationally and as a useful intuitive guide for understanding the physical effect of the heat fluxes.

2.1.1. Energy conservation

With $T_{e}=0$, (2.1)–(2.5) conserve the total energy $E_{K}+E_{M}+E_{\rm th}$, where

(2.9a--c)\begin{equation} E_{K} \doteq \frac{1}{2}\int {\rm d}\boldsymbol{x}\,\rho |\boldsymbol{u}|^{2},\quad E_{M} \doteq \frac{1}{8{\rm \pi} }\int {\rm d}\boldsymbol{x}\, |\boldsymbol{B}|^{2},\quad E_{\rm th} \doteq \int {\rm d}\boldsymbol{x}\left(p_{{\perp}} + \frac{p_{\|}}{2}\right). \end{equation}

The rate of change of the thermal energy is

(2.10)\begin{equation} \frac{{\rm d} E_{\rm th}}{{\rm d} t} = \int {\rm d}\boldsymbol{x}({-}p_{{\perp}}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} + \Delta p\, \hat{\boldsymbol{b}}\hat{\boldsymbol{b}}:\boldsymbol{\nabla} \boldsymbol{u} ) = \int {\rm d}\boldsymbol{x}\left(p_{\|}\frac{1}{\rho}\frac{{\rm d} \rho}{{\rm d} t} + \Delta p\frac{1}{B}\frac{{\rm d} B}{{\rm d} t}\right),\end{equation}

where ${\rm d}/{\rm d} t=\partial _{t}+\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the Lagrangian derivative and we used the continuity and induction equations to derive the second expression (see (2.13)). We see that a positive correlation between $\Delta p$ and ${\rm d} B/{\rm d} t$ drives net heating of the plasma. This is not necessarily guaranteed: in essence, the term is similar to the compressive term $p_{\|} {\rm d} \ln \rho /{\rm d} t$ and can mediate oscillatory transfer between mechanical and thermal energy through pressure anisotropy. However, when collisions dominate the evolution of $\Delta p$, its effect becomes that of a parallel viscosity and $\Delta p\,{\rm d} \ln B/{\rm d} t$ is guaranteed to be positive (see (2.16) and discussion below). For this reason, we will often refer to this term as ‘viscous heating’, even in the collisionless regime where it is more general. Note that, because ion-ion scattering must conserve energy, collisions only indirectly influence the thermal-energy evolution by changing the correlation between $\Delta p$ and ${\rm d} B/{\rm d} t$.

The case with $T_{e}\neq 0$ is addressed in Appendix A.1.

2.1.2. Wave behaviour

Equations (2.1)–(2.5) admit five types of linear-wave solutions: shear-Alfvén waves, two modified magnetosonic-like waves and two types of non-propagating entropy modes. These are discussed in more detail in Appendix A.1.2 (see also Hunana et al. 2019; Majeski, Kunz & Squire 2023), where we give the relevant ideal dispersion relations (computed from (2.1)–(2.5) with $q_\perp =q_\|=0$ and $\nu _c=0$) and use their non-ideal properties as a test of the numerical solvers. Here, we simply note that heat fluxes and/or collisions strongly damp all modes except for Alfvén waves, which, so long as $4{\rm \pi} \Delta p>-B^2$, propagate undamped with the modified speed

(2.11)\begin{equation} {v}_{A,{\rm eff}} = \sqrt{1+\frac{4{\rm \pi}\Delta p}{B^{2}}} v_{A} = \sqrt{{\varTheta}} \,v_{A}.\end{equation}

Linear Alfvénic modes are undamped because they do not involve any perturbation of $B^{2}$ and, therefore, do not create any pressure anisotropy. We see that for $4{\rm \pi} \Delta p<-B^{2}$, ${v}_{A,{\rm eff}}$ becomes imaginary, which is the fluid manifestation of the firehose instability.

2.1.3. The collisionless, weakly collisional and Braginskii-MHD regimes

It is helpful to examine the equation for the evolution of the pressure anisotropy, which may be obtained from (2.4) and (2.5)

(2.12)\begin{align} \frac{{\rm d} }{{\rm d} t}\Delta p & = (p_{{\perp}}+2p_{\|})\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}:\boldsymbol{\nabla} \boldsymbol{u} - (2p_{{\perp}}-p_{\|})\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} \nonumber\\ & \quad -\{\boldsymbol{\nabla}\boldsymbol{\cdot}[(q_{{\perp}}-q_{\|})\hat{\boldsymbol{b}}] + 3q_{{\perp}}\boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{b}} \} -\nu_{c}\Delta p. \end{align}

Each term on the right-hand side is grouped according to its physical effect, which we discuss in turn below.

1. (i) Changing field strength – the first term on the right-hand side of (2.12) captures the creation of pressure anisotropy through the parallel rate of strain $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$, which is related to changes in the magnetic-field strength through

   (2.13)\begin{equation} \frac{1}{B}\frac{{\rm d} B}{{\rm d} t} = \hat{\boldsymbol{b}}\hat{\boldsymbol{b}}:\boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\end{equation}
   (this equation is obtained from (2.3) after dotting it with $\hat {\boldsymbol {b}}$ and rearranging). Due to conservation of the collisionless adiabatic invariants $p_{\perp }/\rho B$ and $p_{\|}B^{2}/\rho ^{3}$, positive $\Delta p$ is created in regions of increasing field strength, while negative $\Delta p$ is created in regions of decreasing field strength.

2. (ii) Compressions – the second term is similar to the first but relates to changes in the plasma density. It is less important here, because of our focus on high-$\beta$ plasmas with low compressibility.

3. (iii) Heat fluxes – the third term involves the heat fluxes, which, although neglected in the so-called double-adiabatic model, are always large in the $\beta \gtrsim 1$ regime where $\Delta p$ has a dynamically important effect (e.g. Hunana et al. 2019). The effect of this term can be understood by examining the ‘Landau fluid’ form of the heat fluxes given in (2.7)–(2.8). When $\Delta p\ll p_{0}$ and the spatial variation of density and $\hat {\boldsymbol {b}}$ is small compared with that of the temperature, these take the approximate forms

   (2.14)\begin{gather} -\boldsymbol{\nabla}\boldsymbol{\cdot}(q_{{\perp}}\hat{\boldsymbol{b}} ) \approx{-}\hat{\boldsymbol{b}} \boldsymbol{\cdot}\boldsymbol{\nabla} q_{{\perp}} \approx \sqrt{\frac{2}{\rm \pi}} \boldsymbol{\nabla}_{\|}\left[\frac{c_{s\|}^{2}}{c_{s\|}|k_{\|}|+a_{{\perp}} \nu_{c}}\boldsymbol{\nabla}_{\|}p_{{\perp}}\right], \end{gather}
   (2.15)\begin{gather}-\boldsymbol{\nabla}\boldsymbol{\cdot}(q_{\|}\hat{\boldsymbol{b}} ) \approx{-}\hat{\boldsymbol{b}} \boldsymbol{\cdot}\boldsymbol{\nabla} q_{\|} \approx \sqrt{\frac{8}{\rm \pi}}\boldsymbol{\nabla}_{\|}\left[\frac{c_{s\|}^{2}}{c_{s\|}|k_{\|}|+a_{\|}\nu_{c}} \boldsymbol{\nabla}_{\|}p_{\|}\right], \end{gather}
   where $a_{\perp }$ and $a_{\|}$ are order-unity numerical factors (see (2.7) and (2.8)). We see that, broadly speaking, the heat fluxes act like a parallel diffusion operator on $p_{\perp }$ and $p_{\|}$ (albeit a nonlocal one if $c_{s\|}|k_{\|}|\gtrsim \nu _{c}$), thus reducing the spatial variation of $\Delta p$ along the magnetic field.

4. (iv) Collisions – the final term in (2.12) reduces the pressure anisotropy due to collisions, and, if large enough, causes (2.1)–(2.5) to reduce to standard, adiabatic MHD.

From this discussion, we see that the first two terms on the right-hand side of (2.12) are responsible for creating a pressure anisotropy from the motions of the plasma. By comparing the relative sizes of the other terms – ${\rm d}\Delta p/{\rm d} t$, the heat fluxes and the collisional term – one finds three regimes each with qualitatively different $\Delta p$ evolution. To distinguish these, it is helpful to consider $\Delta p$ structures of parallel scale $k_{\|}^{-1}$ and assume that $\boldsymbol {\nabla }_{\|} p_{\perp }\sim \boldsymbol {\nabla }_{\|} p_{\|}\sim \boldsymbol {\nabla }_{\|}\Delta p$, viz. that the relative variation in $p_{\perp }$, $p_{\|}$ and $\Delta p$ is of similar magnitudes (but note that it is not true that $p_{\perp }\sim p_{\|}\sim \Delta p$ because the variation in $p_{\perp }$ and $p_{\|}$ is small compared with their mean). Assuming Alfvénic motions, the time derivative scales as ${\rm d}\Delta p/{\rm d} t\sim \omega _{A}\Delta p$ where $\omega _{A}=k_{\|}v_{A}$ is the Alfvén frequency, the heat-flux term scales as ${\sim }[(k_{\|}c_{\rm s})^{2}/(|k_{\|}|c_{\rm s}+\nu _{c})] \Delta p \sim [\beta \omega _{A}/(\beta ^{1/2} + \nu _{c}/\omega _{A})] \Delta p$, and the collisional term is just $\nu _{c}\Delta p$. Restricting our discussion to $\beta \gtrsim 1$, the three regimes are:

1. (i) Collisionless, ${\nu _{c}\ll \omega _{A}}$ – if $\nu _{c}\ll \omega _{A}$, the collisional term cannot compete with ${\rm d}\Delta p/{\rm d} t\sim \omega _{A}\Delta p$ to reduce $|\Delta p|$ significantly and so can be neglected. Because $\nu _{c}/\omega _{A}\lesssim 1$, the heat-flux term scales as ${\sim }\beta ^{1/2}\omega _{A}\Delta p$, implying that heat fluxes are always important in this regime, smoothing $\Delta p$ in space rapidly compared with the Alfvénic motions that drive it.

2. (ii) Weakly collisional, ${\omega _{A}\ll \nu _{c}\lesssim \beta ^{1/2}\omega _{A}}$ – for $\nu _{c}\gtrsim \omega _{A}$, we can neglect ${\rm d}\Delta p/{\rm d} t$ in the balance of terms because collisions isotropise the pressure faster than $\Delta p$ can be created. However, if $\nu _{c}$ also satisfies $\nu _{c}\lesssim \beta ^{1/2}\omega _{A}$, the heat fluxes remain collisionless: they retain a similar form to the collisionless regime and thus have a similar effect on the dynamics, scaling as ${\sim }\beta ^{1/2}\omega _{A}\Delta p$, which remains larger than the collisional draining of $\Delta p$ (unless the parallel scales self-adjust; see Squire, Schekochihin & Quataert 2017b). This weakly collisional regime is thus a hybrid collisional–collisionless one: although motions in the plasma are slower than the collision time scales, heat fluxes remain strong and are governed by collisionless physics (Mikhailovskii & Tsypin 1971).

3. (iii) Braginskii MHD, ${\beta ^{1/2}\omega _{A}\ll \nu _{c}}$ – once $\nu _{c}\gg \beta ^{1/2}\omega _{A}$, in addition to neglecting ${\rm d}\Delta p/{\rm d} t$ compared with $\nu _{c}\Delta p$, we see that the heat fluxes take the collisional form, scaling as ${\sim }(\beta \omega _{A}^{2}/\nu _{c})\Delta p$. Further, unlike in the weakly collisional regime, the heat fluxes become subdominant to $\nu _{c}\Delta p$, since $\nu _{c}\gg \beta \omega _{A}^{2}/\nu _{c}$ for $\nu _{c}\gg \beta ^{1/2}\omega _{A}$. Thus at the same point that the heat fluxes take their collisional form, they become subdominant overall. In this regime, assuming incompressibility and $\Delta p\ll p_{\perp }\sim p_{\|}$, (2.12) takes the simple form

   (2.16)\begin{equation} \Delta p \approx \frac{3p_{0}}{\nu_{c}}\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}:\boldsymbol{\nabla} \boldsymbol{u} ,\end{equation}
   which, when inserted into (2.2), yields a parallel viscous stress. Given this is effectively the parallel viscosity of Braginskii (1965), we refer to this regime as ‘Braginskii MHD’.

In previous work (S$+$19), we explored the effect of pressure anisotropy in the Braginskii-MHD regime, because of the simplicity of its physics and computational implementation. However, when combined with the condition that the pressure anisotropy has a dynamically important influence on the turbulence, $\nu _{c}/\omega _{A}\lesssim \beta \delta B_{\perp }^{2}/B_{0}^{2}$ (see § 2.2), the constraint $\nu _{c}/\omega _{A}\gg \beta ^{1/2}$ requires $\beta$ to be very large in order for there to be sufficient range between $\beta ^{1/2}$ and $\beta \delta B_{\perp }^{2}/B_{0}^{2}$. Thus, for application to the solar wind and other hot astrophysical plasmas with $\beta \lesssim 100$, the collisionless and weakly collisional regimes are more relevant.

A detailed derivation of the effects described above for a single shear-Alfvén wave, including simplified asymptotic equations valid in each regime, is given in Squire et al. (2017b).

2.2. Magneto-immutability

In previous work (Squire et al. 2016, 2017a,b), we explored the idea that Alfvénically polarised magnetic fields or flow perturbations are ‘interrupted’ if their amplitude satisfies

(2.17)\begin{equation} \frac{\delta u_{{\perp}}}{v_{A}}\sim \frac{\delta B_{{\perp}}}{B_{0}} \gtrsim2 \beta^{{-}1/2} \begin{cases} 1, & \omega_{A}\gg \nu_{c} \\ \sqrt{{\nu_{c}}/{\omega_{A}}}, & \omega_{A}\ll \nu_{c} \end{cases} \end{equation}

(the first limit applies in the collisionless regime; the second applies in the weakly collisional and Braginskii-MHD regimes). The origin of the effect is straightforward: above the limit (2.17), the (nonlinear) perturbation of the magnetic-field strength drives the pressure anisotropy to the fluid firehose limit, $\Delta p = -B^{2}/4{\rm \pi}$, or $\varDelta = -2/\beta$. As can be seen from the final term of (2.2), this nullifies the plasma's magnetic tension (indeed this is the cause of the firehose instability), which thus robs the Alfvén wave of its restoring force (see (2.11)). The consequence, for a single linearly polarised Alfvén wave, is that the perturbation dumps most of its energy into plasma heating and/or magnetic-field perturbations that cease to evolve in time, rather than oscillating in the usual way (Squire et al. 2016, 2017a). Below the limit (2.17), waves slowly damp due to nonlinear Landau damping (e.g. Hollweg 1971) and/or nonlinear viscous damping (Nocera, Priest & Hollweg 1986; Russell 2023).

2.2.1. Alfvénic turbulence

Given that Alfvénic motions underlie magnetised plasma turbulence (Chen 2016; Schekochihin 2022), a natural question that arises is: what happens when large-scale random perturbations to the plasma are driven past the interruption limit (2.17)? Naïvely, one might expect that only motions below the limit (2.17) would be allowed, which would imply that the large-scale fluctuation energy would be limited to less than $\delta u_{\perp }^{2}\sim v_{A}^{2}/\beta$, with the remainder of the energy injected at large scales directly heating the plasma without creating smaller-scale motions and a turbulent cascade. Instead, we showed in S$+$19 that the turbulence rearranges itself mostly to avoid the motions that would create large pressure anisotropies in the first place, allowing the turbulent cascade to proceed in a way rather similar to MHD. It does this by reducing the variations in $B$ that would have driven large $\Delta p$, creating a turbulence in which $B$ varies significantly less than in turbulence where the pressure anisotropy does not back react on the plasma motions. This reduces the spread of $\Delta p$ produced by the turbulence, thus reducing the plasma heating done by these motions and increasing the ‘cascade efficiency’ (the proportion of energy that cascades to small scales).

That the plasma does this is not altogether surprising. Indeed, it is well known that collisionless effects damp out motions that involve variations in the magnetic-field strength, and that this damping is fast compared with Alfvénic time scales at high $\beta$ (Barnes 1966; Foote & Kulsrud 1979)². However, the effect does not simply involve a selective damping of those motions that involve variations in $B$, thereby leaving behind those motions that do not. Rather, there is a direct force on the plasma due to the final term $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ in (2.2), and this force opposes motions that would be strongly damped (those involving variation in $B$). The heating is thus significantly reduced compared with what would occur in the absence of this force. The origin of this behaviour is best understood by analogy to compressive motions and density fluctuations. It is well known that isotropic pressure forces in a fluid resist compressional flows with $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}\neq 0$: such a flow will generate a local pressure perturbation, which then (through the $-\boldsymbol {\nabla } p$ term) generates a force on the fluid that opposes the compressional motion. In a fluid where the thermal energy is large compared with other energies (such as a plasma at high $\beta$), this process is rapid compared with the time scales on which the flow or magnetic-field change, thus rendering the system effectively incompressible. Magneto-immutability involves a similar process, but with ‘magneto-dilational’ flows that have $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u} \neq 0$. Such flows generate a local pressure-anisotropy perturbation (see (2.12)), the feedback of which – through the pressure-anisotropy force $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ – opposes the motion that created $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ in the first place. If the generation of pressure anisotropy is fast compared with the Alfvénic time scales of the turbulence, these forces will render the system ‘magneto-immutable’, since motions with small $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ are those that minimise changes to the magnetic-field strength (see (2.13)).³

2.2.2. A reduced model for magneto-immutable turbulence?

The analogies between incompressibility and magneto-immutability lead one to speculate whether there could exist a reduced model – similar to incompressible hydrodynamics – that describes magneto-immutable turbulence. Incompressible fluid models are formulated by stipulating that, because the compressible back reaction happens so rapidly, the pressure force is just what it needs to be in order to enforce $\partial _{t}(\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u})= 0$. This lets one solve for $\boldsymbol {\nabla }^{2}p$ in terms of $\boldsymbol {u}$, thus closing the system. By analogy, for a magneto-immutable fluid model, we should solve for the $\Delta p$ in $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ that enforces $\partial _{t}(\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u} ) = 0$, which will ensure that the flow cannot generate magneto-dilations as it evolves. This immediately reveals a technical problem: while $\boldsymbol {\nabla }\boldsymbol {\cdot }$ is a simple operator, thus enabling straightforward solution of $p$, the combination $\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} :\boldsymbol {\nabla }$ is complicated, and solving for $\Delta p$ (which must be achieved at every time step for a numerical algorithm) becomes complex and expensive.

More generally, there is another key difference with incompressibility that argues against the utility of formulating a magneto-immutable fluid model: regions with large pressure anisotropies become unstable. The strong back-reaction force needed to suppress a large $|\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u} |$ in some region will require a large $|\Delta p|$, which will then grow small-scale instabilities, presumably tempering its influence. This effect is unavoidable because such instabilities are always triggered when $|\Delta p|\gtrsim B^{2}/4{\rm \pi}$ (see § 2.3), which is also the pressure anisotropy needed to start feeding back significantly on the flow. In contrast, in a large compression or rarefaction, the distribution function can remain isotropic and thus stable, so there is no similar effect for incompressibility. This subtle, but important, difference between incompressibility and magneto-immutability is discussed in more detail in §§ 2.4 and 5.1.

2.2.3. Interruption number

It will prove helpful to have a simple dimensionless number that can be used to quantify the expected influence of the pressure anisotropy on the flow. To construct this, we return to the idea that individual shear-Alfvén waves are ‘interrupted’ – meaning that they dissipate a large fraction of their energy into thermal energy within ${\lesssim } \tau _{A}$ – if their amplitude satisfies (2.17), or

(2.18)\begin{equation} \frac{\delta B_{{\perp}}}{B_{0}}\sim \frac{\delta u_{{\perp}}}{v_{A}}\gtrsim \mathcal{M}_{A}^{\rm int} \doteq \left[\max\left(1,\frac{\nu_{c}}{\omega_{A}}\right) \frac{1}{\beta}\right]^{1/2}. \end{equation}

Applying this limit to a random turbulent collection of fluctuations with root-mean-square (rms) amplitudes $\mathcal {M}_{A}\doteq \langle \delta u_{\perp }^{2}\rangle ^{1/2}/v_{A0} \approx \langle \delta B_{\perp }^{2}\rangle ^{1/2}/B_{0}$, we define the ‘interruption number’ to be

(2.19)\begin{equation} \mathcal{I} \doteq \left(\frac{\mathcal{M}_{A}^{\rm int}}{\mathcal{M}_{A}}\right)^{2}= \frac{\max\left(1,\dfrac{\nu_{c}}{\omega_{A}}\right)}{\beta\mathcal{M}_{A}^{2}}.\end{equation}

If $\mathcal {I}\lesssim 1$, the turbulence (unless otherwise restricted) should be of sufficiently large amplitude to generate $|\varDelta |\gtrsim 2/\beta$, which would be naïvely expected to damp out the energy faster than it cascades to small scales. In this sense, $\mathcal {I}$ can be interpreted similarly to the Reynolds number, capturing viscous-like effects due to the pressure anisotropy, with $\mathcal {I}\lesssim 1$ suggesting they dominate over the Alfvénic forces in the flow.

However, as described above, this expression ignores the feedback of the pressure anisotropy (magneto-immutability), which reduces the production of $\Delta p$ below the estimate used to derive (2.19). Thus, as shown by S$+$19, the turbulent cascade can in general proceed when $\mathcal {I}\lesssim 1$, meaning $\mathcal {I}$ is better thought of as an estimate of the importance of pressure-anisotropy effects, rather than whether the viscous damping dominates the dynamics (in other words, when $\mathcal {I}\lesssim 1$, the flow can rearrange itself so as to avoid the motions that would be strongly viscously damped). The purpose of our study is to understand the properties of turbulence in this $\mathcal {I}\lesssim 1$ regime and characterise how it heats the plasma.

Another way to think of the interruption number is by analogy to the Mach number of a compressible neutral fluid. Start by considering pressure-anisotropy production with the balance $\omega _{A} \Delta p \sim p_{0}\omega _{A} \delta B_{\perp }^{2}/B^{2}+ \nu _{c}\Delta p$ obtained from (2.12) (ignoring the heat fluxes) and giving the estimate $\Delta p \sim p_{0}\mathcal {M}_{A}^{2}/\max (1,\nu _{c}/\omega _{A})$ for the size of $\Delta p$ fluctuations that result from Alfvénic magnetic-field fluctuations. Because pressure anisotropy will feedback strongly on the plasma once $\Delta p\gtrsim B^{2}$, and the estimate above gives $\Delta p/B^{2}\sim \mathcal {I}^{-1}$, we see that $\mathcal {I}^{-1}$ is the ratio between the $\Delta p$ that is driven and the $\Delta p$ that would substantially change the flow. Analogously, in compressible hydrodynamics, isotropic pressure fluctuations $\delta p$ are related to density fluctuations $\delta \rho$ by $\delta p\sim c_{s}^{2} \delta \rho$, where $c_{s}$ is the sound speed. Pressure fluctuations will feed back strongly on the flow once they are comparable to the ram pressure $\delta p\sim \rho u^{2}$. Therefore, the ratio of the naturally generated pressure fluctuations to those needed to feed back on the flow is $\delta p/\rho u^{2}\sim \mathcal {M}^{-2}\delta \rho /\rho$, where $\mathcal {M}=u/c_{\rm s}$ is the Mach number. Thus, the isotropic hydrodynamic equivalent of $\mathcal {I}$ is $\mathcal {M}^{2}/(\delta \rho /\rho )$. However, in a completely unrestricted flow, $\delta \rho /\rho \sim 1$ (because the turnover rate ${\sim }ku$ at scale $k$ scales in the same way as $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}$), showing that $\mathcal {M}^{2}$ itself provides the relevant analogy with $\mathcal {I}$. This makes intuitive sense: $\mathcal {M}\sim 1$ marks the boundary between supersonic turbulence, which approaches the limit of the pressure-free Burgers equation, and incompressible turbulence, where $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}$ and $\delta \rho /\rho$ become strongly restricted by the feedback of the pressure force on the flow.

We can also define a local interruption number with respect to the scale-dependent amplitude of the eddies, assuming that the turbulence follows a standard Goldreich–Sridhar cascade (Goldreich & Sridhar 1995). This assumption is clearly highly questionable if the effects of pressure anisotropy are strong, but our simulations will show it to be nevertheless reasonable and it serves a useful purpose for basic estimates. Turbulent amplitudes scale as $\delta B_{\perp }^{2}\propto \delta u_{\perp }^{2}\propto k_{\perp }^{-2/3}$, where $k_{\perp }$ is the inverse scale of an eddy perpendicular to the local magnetic field. This implies that, in the collisionless regime where $\omega _{A} \varDelta \propto \omega _{A} \delta B_{\perp }^{2}/B_{0}^{2}$ (see (2.12) and (2.13)), the interruption number should be smallest at the outer scale and grow as $\mathcal {I}\propto k_{\perp }^{2/3}$. By contrast, in the Braginskii regime where $\nu _{c}\varDelta \propto \omega _{A} \delta B_{\perp }^{2}$, the critical balance scaling $\omega _{A}\propto k_{\|}\propto k_{\perp }^{2/3}$ implies that the interruption number is effectively constant across all scales above the collisionless transition where $\omega _{A}\gtrsim \nu _{c}$.

2.3. Microinstabilities, limiters and the choice of collisionality

A key physical effect that has been omitted in the discussion above is the influence of kinetic microinstabilities. Most important at high $\beta$ are the firehose and mirror instabilities (Rosenbluth 1956; Chandrasekhar, Kaufman & Watson 1958; Parker 1958; Vedenov & Sagdeev 1958; Hasegawa 1969), which are thought to be the fastest growing with the strongest back reaction on the large-scale plasma dynamics. In their simplest forms, they are triggered for $\varDelta \lesssim -2/\beta$ (firehose) or $\varDelta \gtrsim 1/\beta$ (mirror), with minor modifications to these limits from fundamentally kinetic effects (resonances and finite-Larmor-radius physics), particularly at moderate $\beta$ (Yoon, Wu & de Assis 1993; Hellinger et al. 2006). Although most features of the linear instability thresholds and growth rates are captured by (3.1a,b) and (2.5) with the Landau-fluid heat fluxes (2.7) and (2.8) (see Snyder et al. 1997), the detailed nonlinear saturation mechanisms, which involve particle scattering and trapping (Schekochihin et al. 2008; Kunz et al. 2014; Rincon et al. 2015; Melville, Schekochihin & Kunz 2016), are certainly not. Even more important is the separation of time scales that is inherent in how the microinstabilities feed back on the plasma: they grow and evolve on time scales comparable to the ion gyro-frequency, which is far faster than any motions related to the outer scale for any astrophysical system of interest. Thus, as far as the large-scale ($\ell \gg \rho _{i}$) plasma dynamics is concerned, microinstabilities should saturate and feed back on the plasma effectively instantaneously. This poses an extreme difficulty for fully kinetic simulations, which must attempt to determine which observed features are dependent on the necessarily modest scale separation, and which are robust to asymptotically large scale separations (e.g. Kunz et al. (2020); A$+$22).

While the details of firehose and mirror saturation are rather complex (Hellinger & Trávníček 2008; Schekochihin et al. 2008; Kunz et al. 2014; Riquelme, Quataert & Verscharen 2015; Sironi & Narayan 2015; Melville et al. 2016), roughly, the process involves the instabilities’ magnetic fluctuations scattering particles at the rate needed to maintain $\Delta p$ at its marginal level.⁴ This implies an additional microinstability-induced collisionality, $\nu _{c}\sim S\beta$, where $S\sim B^{-1}\,{\rm d} B/{\rm d} t\approx \hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ is the shearing rate, which operates in regions where $\Delta p$ is being driven beyond the mirror or firehose thresholds. Broadly speaking, such a picture has been commonly invoked to understand the solar-wind ‘Brazil plots’ (Kasper, Lazarus & Gary 2002; Hellinger et al. 2006; Bale et al. 2009), which show how the measured $\Delta p$ appears to be limited between the instability threshold values ($|\varDelta |\lesssim 1/\beta$; see figure 1): plasma that strays beyond the boundaries will be rapidly pushed back via scattering, thus maintaining only a small deviation from $\Delta p\approx 0$ at high $\beta$.⁵ Since the nominal effect of this scattering is simply to maintain $\Delta p$ at marginality, a simple phenomenological method to capture such effects in a fluid simulation is via the inclusion of ‘limiters’, which halt the growth of $\Delta p$ locally in space whenever it is driven past the firehose or mirror thresholds (Sharma et al. 2006). Physically, the approach assumes that (i) microinstabilities act quasi-instantaneously to return the plasma to marginal stability⁶ , and (ii) that the microinstabilities do not directly influence the plasma's evolution outside of the regions that are being driven unstable, either in space or time. Despite its clear shortcomings, which will be discussed in detail in § 5.1 after we present our computational results, the method is at least simple and well controlled, and we will use it throughout this work.

Figure 1. ‘Brazil plots’ – histograms of (local) $\beta$ versus $p_{\perp }/p_{\|}$ – formed by combining the data from all runs that have both active-$\Delta p$ and passive-$\Delta p$ simulations with the same parameters, as indicated by the ‘Passive?’ column of table 1 (for passive-$\Delta p$ simulations, the $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u} )$ force in (2.2) was artificially removed). The left panel combines all such active-$\Delta p$ simulations; the right panel combines all such passive-$\Delta p$ simulations. The colour (probability) scale is normalised differently for each $\beta$ to better show the main features, but is done so identically for the left and right panels. The insets zoom in on the right-hand regions of each plot. The clear differences between the active- and passive-$\Delta p$ simulations highlight the effect of the pressure-anisotropy feedback on the turbulence: active-$\Delta p$ cases maintain themselves primarily within the microinstability boundaries with small $|\Delta p|$, while passive-$\Delta p$ simulations (which have identical parameters otherwise, and similar turbulence amplitudes) have most of their volume artificially constrained by the hard-wall limiters. The basic conclusion is that the plasma employs dynamical pressure-anisotropy feedback, in addition to particle scattering (as more commonly discussed; Kasper et al. 2002; Hellinger et al. 2006; Bale et al. 2009), in order to limit its deviations from local thermodynamic equilibrium and stay within the kinetically stable region of parameter space.

2.4. The expected impacts of magneto-immutability

In this subsection we summarise the basic impact of pressure-anisotropy feedback (magneto-immutability) by comparison with the counterfactual situation where it does not exist. Because these effects tend to cause the plasma to revert to behaviour that more closely resembles the collisional (MHD) limit, they can be somewhat subtle and not easily diagnosed. Nonetheless, their influence on the heating processes and turbulent statistics can be strong and appreciation of it is needed to understand the behaviour of turbulent collisional high-$\beta$ plasmas.

As explained in § 2.2, the basic picture involves the plasma rapidly reacting to suppress ‘magneto-dilational’ flows with large $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$, which are those that would create large pressure anisotropies. The effect is very similar to incompressibility, if we substitute $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}$ with $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$, isotropic pressure $p$ with $\Delta p$ and the $-\boldsymbol {\nabla } p$ force with that from $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$. As a consequence:

1. (i) The standard deviation of $\Delta p$ will be suppressed, viz. it will be lower than if $\Delta p$ were driven by a turbulent flow with similar $\mathcal {M}_{\rm A}$ but that did not feel the force $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ (e.g. MHD turbulence).

2. (ii) The standard deviation of $B=|\boldsymbol {B}|$ will be suppressed in the same way as (i) (e.g. compared with a similar-amplitude MHD case). This is because suppressing $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ also suppresses changes of $B$ (fundamentally, it is the changing $B$ that drives $\Delta p$).

3. (iii) The primary influence on the flow statistics will be the suppression of magneto-dilations ($\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$) compared with MHD.

4. (iv) As a consequence, the net viscous-like heating of the plasma through pressure anisotropy will be suppressed by magneto-immutability (compared with the counterfactual situation where $\Delta p$ did not influence the flow). Since such heating is dominated by the outer scales in Alfvénic turbulence (see § 2.2.3), there will thus be more vigorous turbulence with a larger cascade efficiency, and therefore a larger fraction of heating will occur via kinetic processes at the smallest scales (Schekochihin et al. 2009). This may influence bulk thermodynamical properties such as the ion-to-electron heating ratio or parallel-to-perpendicular heating ratios (cf. Sharma et al. 2007; Howes et al. 2008; Kawazura, Barnes & Schekochihin 2019; Kawazura et al. 2020) or even the plasma's thermal stability (Kunz et al. 2010).

While we provide solid numerical evidence in support of each of these points (i)–(iv) in § 4, it is worth clarifying that these effects can never dominate all other dynamics, because there is no physical limit in which any of the aforementioned ‘suppressions’ becomes complete. The reason for this is discussed in § 2.3, with more detail to arrive in § 5.1: when pressure anisotropies become large, they cause plasma microinstabilities, which act to limit $\Delta p$ through scattering, trapping and/or microscale fields. Such effects always occur together with the direct $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ feedback of $\Delta p$ on the flow.⁷ This implies that even in the limit of extremely high $\beta$ and large amplitudes (small $\mathcal {I}$), $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ can never be arbitrarily strongly suppressed, because the $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ force that causes this suppression is also attenuated by the microinstabilities. This is the most important difference between magneto-immutability and incompressibility; the latter does not suffer the same fate because isotropic pressure changes (driven by $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}$) are neither kinetically unstable nor attenuated by particle scattering (see § 2.2.2). Similarly, related to point (iv) above, we still expect (and measure) significant viscous heating at the outer scale (a cascade efficiency below unity); but, the cascade efficiency is independent of $\mathcal {I}$ and larger than what would be expected without the $\Delta p$ feedback (in which case it would decrease continuously with $\mathcal {I}$, with turbulence amplitudes limited to ${\ll }\mathcal {M}_{A}^{\rm int}$; see (2.17)).

A simple, familiar visual illustration of the effects of magneto-immutability is shown in figure 1, where we combine most of the simulations run as part of this work into two ‘Brazil’ plots. These are two-dimensional (2-D) probability distribution functions (PDFs) of $\beta$ and temperature (pressure) anisotropy, illustrating the approximate magnitude of the deviation from pressure isotropy as a function of $\beta$. As mentioned earlier, its eponymous shape has usually been attributed to the action of microinstabilities scattering particles once $|\varDelta |\gtrsim \beta ^{-1}$ (the mirror and firehose thresholds; Kasper et al. 2002; Hellinger et al. 2006; Bale et al. 2009). The purpose of figure 1 is to compare simulations of (2.1)–(2.5) (a) with equivalent simulations where the dynamic feedback of $\Delta p$ has been artificially eliminated by removing $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ in (2.2) and assuming an isothermal ion-pressure response (‘passive-$\Delta p$’ simulations; see § 3.1.3). Both sets of simulations include the effect of microinstabilities via hard-wall limiters at the firehose and mirror boundaries, which prevent $\Delta p$ from straying beyond the relevant instability thresholds. We see a clear difference between the two cases, with $\Delta p$ distributions strongly dominated by regions at the microinstability boundaries (the artificial limiters) in the passive cases, but not when $\Delta p$ feeds back on the flow. This demonstrates that, in its efforts to remain close to local thermodynamic equilibrium, the plasma has two separate methods at its disposal. The first – particle scattering through microinstabilities – has been discussed by Kasper et al. (2002), Hellinger et al. (2006), Bale et al. (2009) and other subsequent works. The second – the dynamic feedback from the pressure anisotropy – should be of similar importance for maintaining near isotropy in most plasmas (see § 5.1), but has been largely ignored thus far.

3. Methods

3.1. The CGL Landau-fluid model

Our computational study is based on the ‘CGL-LF’ model of Snyder et al. (1997). This allows us to probe the effects discussed above without the complications of a true kinetic model (see § 2.3). The model solves (2.1)–(2.5) supplemented by the Landau-fluid closure for heat fluxes (described in § 3.1.1) and simple ‘hard-wall’ limiters on $\Delta p$ to approximate the effect of kinetic microinstabilities (§ 3.1.2).

We use the finite-volume Athena++ code (White et al. 2016; Stone et al. 2020), modified to solve (2.1)–(2.5) in the conservative form detailed in Appendix A. Briefly, this uses the total energy $E = p_{\perp }+p_{\|}/2+ B^{2}/8{\rm \pi} + \rho |\boldsymbol {u}|^{2}/2$ and the ‘anisotropy’ $\mathcal {A} \doteq \rho \ln (p_{\perp } \rho ^{2}/p_{\|}B^{3})$ as conserved variables, the latter following Santos-Lima et al. (2014). We found through extensive numerical tests that using $\mathcal {A}$ leads to a solver that is more numerically robust than if $p_{\perp }/B$ or $p_{\|}B^{2}/\rho ^{2}$ is used as the second conserved variable, particularly when there exists significant variation in $B$. Given that $\mathcal {A}$, and the equations themselves, become ill defined for $B\rightarrow 0$, we implement a numerical floor on $B$, below which the equations revert to standard adiabatic MHD. We use the piecewise parabolic method with an HLL Riemann solver (Toro 2009); although we developed and tested various other HLLD-inspired solvers (Miyoshi & Kusano 2005; Mignone 2007), these were found to be insufficiently robust for this study (see Appendix A.3.1).

The Riemann solver solves only the conservative part of (2.1)–(2.5), viz. that with $q_{\perp }=q_{\|}=0$ and $\nu _{c}=0$. We evaluate the heat fluxes in the form described below (see (3.2a,b)) using operator splitting with slope limiters to ensure numerical stability and monotonicity of the solutions (Sharma & Hammett 2007; Dong & Stone 2009). We use the RKL1 super-time-stepping algorithm of Meyer, Balsara & Aslam (2014) to allow for much larger time steps by stepping over the Courant–Friedrichs–Lewy (CFL) condition that results from the small-scale diffusive form of the computational heat fluxes (although heat-flux-related time-step constraints are still relatively severe at high $\beta$ and high resolution). Collisional terms are evaluated at the end of each global time step $\delta t$ from the exact solution of $\partial _{t}p_{\perp } = -\nu _{c}\Delta p/3$, $\partial _{t}p_{\|} = 2\nu _{c}\Delta p/3$ from $t$ to $t+\delta t$ (see (A4)). This method is implicit and numerically stable for any time step, and so has the property that adiabatic MHD can be easily recovered from the CGL-LF system by setting $\nu _{c}\gg \delta t^{-1}$ (which also causes the heat fluxes to become negligible). Similarly, if we choose $\nu _{c}\gtrsim \beta ^{1/2}\omega _{A}$ (see § 2.1.3), the method becomes a convenient, stable and computationally efficient way to include a Braginskii-viscous stress in the standard MHD system.

Turbulence is driven via a large-scale incompressible forcing term $\boldsymbol {F}$ added to the right-hand side of (2.2). This applies only in the directions perpendicular to the mean magnetic field in order to drive primarily Alfvénically polarised fluctuations, and consists of the eight largest-scale modes in the box evolved in time as an Ornstein–Uhlenbeck process with correlation time $t_{\rm corr}=\tau _{A}$ (where $\tau _{A}$ is the box-scale Alfvén time). At each time step, its amplitude is adjusted to enforce a constant energy-injection rate $\varepsilon = \langle \rho \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {F}\rangle$ (see, e.g. Lynn et al. 2012).

3.1.1. Heat fluxes

In the ‘Landau-fluid’ heat fluxes (2.7)–(2.8), the $k_{\|}$ and $\boldsymbol {\nabla }_{\|}=\hat {\boldsymbol {b}} \boldsymbol {\cdot }\boldsymbol {\nabla }$ terms contain $\hat {\boldsymbol {b}}$, which varies in space along with $c_{s\|}^{2}=p_{\|}/\rho$. This implies that that $2 c_{s\|}^{2}({\sqrt {2 {\rm \pi}}c_{s\|} |k_{\parallel }|+\nu _{c}})^{-1}$ and $8 c_{s\|}^{2}[\sqrt {8 {\rm \pi}}c_{s\|} |k_{\parallel }|+(3{\rm \pi} -8)\nu _{c}]^{-1}$ should rightly be considered operators that are not diagonal in either Fourier space or in real space, making them complex and computationally expensive to evaluate (see Snyder et al. 1997). For this reason, we use a simplified form, motivated by Sharma et al. (2006), in which $|k_{\|}|$ in the denominators of (2.7) and (2.8) is replaced by a constant $k_{L}$, which we take to be $2\times 2{\rm \pi} /L_{\|}$ to approximate $k_{\|}$ of larger-scale motions.⁸ However, this substitution also leads to the undesirable property that, due to large parallel gradients at small scales, the heat flux, which is now diffusive, can become ${\gg }\rho c_{s\|}^{3}$, the maximum possible heat flux (Hollweg 1974a; Cowie & McKee 1977). To mitigate this, we additionally limit $q_{\perp }$ and $q_{\|}$ to their maximum possible value deduced from (2.7) and (2.8), viz.

(3.1a,b)\begin{equation} q_{{\perp},\max} \doteq \sqrt{\frac{2}{\rm \pi}} c_{s\|} p_{{\perp}},\quad q_{\|,\max} \doteq \sqrt{\frac{8}{\rm \pi}} c_{s\|} p_{\|} \end{equation}

(the second term in (2.7) is ${\sim }\Delta \delta B_{\|}/B\ll 1$ times smaller than the first, so we ignore its contribution). The heat fluxes are thus computed as

(3.2a,b)\begin{equation} \tilde{q}_{{\perp}} = q_{{\perp}{L}}\frac{q_{{\perp},\max}}{q_{{\perp},\max}+|q_{{\perp}{L}}|},\quad \tilde{q}_{\|} = q_{\|{L}}\frac{q_{\|,\max}}{q_{\|,\max}+|q_{\|{L}}|},\end{equation}

where $q_{\perp {L}}$ and $q_{\|{L}}$ are evaluated using (2.7) and (2.8) with $|k_{\|}|=k_{L}$. We note that this approach effectively generalises the practice of limiting an MHD heat flux to be $|\boldsymbol {q}|\leq 5\phi \rho c_{s}^{3}$ with $\phi \approx 0.3$ (Cowie & McKee 1977),⁹ which is commonly used in MHD simulations (e.g. Vaidya et al. 2017).

The method (3.2a,b) is somewhat ad hoc, which led us to explore various other approaches in some detail. One possibility for some regimes is to evaluate $|k_{\|}|$ along the constant mean field using Fourier transforms (Passot et al. 2014; Finelli et al. 2021), thus effectively approximating $|k_{\|}|=|\hat {\boldsymbol {b}} _{0}\boldsymbol {\cdot }\boldsymbol {\nabla }|$, rather than $|k_{\|}|=k_{L}$, in (2.7) and (2.8). However, with extensive testing, we found this approach to be more prone to numerical instability and more computationally expensive than the simpler method described above, when the fluctuations are large compared with the background magnetic field (as is explored here). In any case, because heat fluxes will be strongly modified by scattering from microinstabilities (see § 3.1.2 below), an exact evaluation of the Landau-fluid form (2.7)–(2.8) is likely irrelevant for detailed agreement to nonlinear collisionless physics, so long as the method captures their general influence on the flow (see § 2.1.3).

3.1.2. Microinstability limiters

As discussed in § 2.3, although aspects of the firehose and mirror instabilities are captured by (2.1)–(2.5) with the closure (2.7)–(3.2a,b), their nonlinear saturation, which involves particle scattering and trapping, is not. We therefore artificially limit the pressure anisotropy to

(3.3)\begin{equation} -\varLambda_{\rm FH} \frac{B^{2}}{8{\rm \pi}} < \Delta p < \varLambda_{\rm M} \frac{B^{2}}{8{\rm \pi}}, \end{equation}

where $\varLambda _{\rm FH}$ and $\varLambda _{\rm FH}$ define the firehose and mirror instability thresholds, respectively.¹⁰ By default, we set $\varLambda _{\rm FH}=2$ and $\varLambda _{M}=1$, which describe the canonical versions of these instabilities, but our results do not depend strongly on this choice.¹¹ Computationally, the limiters work by applying a large scattering rate ($\nu ^{{\rm lim}}_{c}=10^{10}\tau _{A}^{-1}$ in all simulations) to any region outside of (3.3), which quickly (within one time step) reduces $\Delta p$ to lie on the relevant instability boundary (see Appendix A.1.1). We also apply this enhanced collisionality to the heat fluxes (2.7)–(2.8), thus strongly suppressing them in limited regions. This may be appropriate for mirror-limited regions, which show strong heat-flux suppression in simulations (Kunz et al. 2020), but it is likely much too strong in firehose-limited regions, which seem to be well described by the Braginskii estimate (Kunz et al. 2020), meaning it might be more appropriate to take $\nu _{c}^{\rm lim}\sim \beta \hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ (though this would be numerically complicated).

3.1.3. Passive-pressure-anisotropy simulations

For most simulations, in addition to the standard CGL-LF model (termed ‘active-$\Delta p$’ below), we have run an otherwise equivalent ‘passive-$\Delta p$’ simulation. The latter is identical to the former, except that the feedback of the pressure anisotropy into the momentum equation is artificially removed. Instead, we use an isothermal equation of state with the sound speed chosen to give the same value of $\beta$ as in the active-$\Delta p$ run. Note that these passive-$\Delta p$ runs still evolve the pressure anisotropy and include heat fluxes in an identical way to the standard CGL-LF model. In this way, the $p_{\|}$ and $p_{\perp }$ statistics can be compared directly to understand the feedback of the pressure-anisotropy stress on the flow.

In summary, a ‘passive-$\Delta p$’ simulation solves (2.1)–(2.5) using the same forcing, initial conditions and parameters as a standard (‘active-$\Delta p$’) run, but, in (2.2) we remove the $\Delta p$ term and set $p_{\perp }=T_{i0}\rho$ with $8{\rm \pi} \rho T_{i0}/B_{0}^{2} =\beta$ (a chosen initial parameter). This implies that the velocity and magnetic-field evolution in such a simulation is described by the isothermal MHD equations.

3.2. Study design

All simulations are run in a $L_{x}=L_{\|}=2L_{\perp }$ aspect ratio box of volume $V$, which has mean density $\rho =\rho _0$ and is threaded by a mean magnetic field $\boldsymbol {B}_{0}=B_0\hat {\boldsymbol {x}}$. The energy-injection (forcing) level $\varepsilon$ is set to $\varepsilon = 0.16 v_{A}^{3}/L_{\perp }$, which is chosen empirically to give $\mathcal {M}_{A}\sim \delta u_{\perp }/v_{A}\sim \delta B_{\perp }/B_{0}\sim 1/2$ in steady state for MHD, as needed for critical balance at the outer scale (from hereon, $v_{A}$ will refer to the mean-field Alfvén speed $v_{A}=B_0/\sqrt {4{\rm \pi} \rho _0}$). We initialise with isotropic pressure $p_{0}$, chosen to yield the desired initial $\beta$, $\beta _{0} = 8{\rm \pi} p_{0}/B_{0}^{2}$. Recall that, here and throughout, $\beta$ and $p_{0}$ refer to the ion contribution, and most simulations use $T_{e}=0$ by default in order to diagnose more easily the influence of pressure anisotropy. We do not use an explicit isotropic viscosity or resistivity in any simulations, relying on the grid to dissipate energy that reaches the smallest scales. Most simulations have a matching passive-$\Delta p$ run, which is set up as explained in § 3.1.3. With these default parameters, each simulation is specified by its initial ion $\beta$ and collisionality $\nu _{c}$ (in units of $v_{A}/L_{\perp }$; these are sometimes omitted for conciseness). These parameters then fix the expected interruption number from (2.19), assuming $\mathcal {M}_{A}$ is not strongly modified by the effects of pressure anisotropy. The full set of simulations is listed in table 1.

We focus particularly on three simulations to probe both the collisionless and weakly collisional regimes with $\mathcal {I}\lesssim 1$ (meaning that pressure-anisotropy feedback is a significant effect). These use $\beta _{0}=10,\ \nu _{c}=0$ (labelled CL10); $\beta _{0}=100, \nu _{c}=33v_{A}/L_{\perp }$ (B100); and $\beta _{0}=100,\ \nu _{c}=0$ (CL100). CL10 and B100, which have a numerical resolution of $1120\times 560^{2}$, are chosen to have similar $\mathcal {I}\simeq 0.4$ but explore two different, collisionless and Braginskii, regimes laid out in § 2.1.3. CL100 (with resolution $560\times 280^{2}$) is chosen to examine a situation with stronger pressure-anisotropy feedback, $\mathcal {I}\approx 0.04$. The extremely small time steps required for the stable evaluation of the heat-flux terms make high-resolution simulations quite costly in wall-clock time, and so CL10 and B100 are initialised from the saturated state of lower-resolution simulations (see below). We then run them for only a relatively short time of $t\approx L_{\perp }/v_{A}$, which is shorter than the outer-scale turnover time ${\simeq }2 L_{\perp }/v_{A}$, but longer than the time, $t\approx 0.1L_{\perp }/v_{A}$, that it takes the new smaller scales to reach turbulent steady state (we observe the evolution of the energy spectrum to ensure that this is the case). This means that these simulations are useful for exploring detailed properties of the turbulence (e.g. turbulence structure and energy transfers), but the largest scales (those above around a quarter of the box scale) are not properly statistically averaged.

In addition to these simulations, we have a large array of other ones at low resolution ($400\times 200^{2}$). These are designed to explore how the physics of magneto-immutability varies with plasma parameters. ‘lrCL’ simulations are all collisionless, changing $\beta _{0}$ from ${\lesssim }1$ to ${\approx }100$ in order to probe the dependence of the turbulence on the interruption number from $\mathcal {I}\approx 0.04$ to $\mathcal {I}\gtrsim 1$. In simulations with $\beta _{0}\lesssim 1$, the plasma is heated, thus increasing $\beta$, rather rapidly, limiting the time they sit near their initial parameters. The ‘lrB’ simulations probe the collisionless, weakly collisional and Braginskii-MHD regimes discussed in § 2.1.3 by keeping $\mathcal {I}\simeq 0.4$ approximately constant while changing $\beta _{0}$ and $\nu _{c}$, with $\beta _{0}$ ranging from $30$ to $600$ and $\nu _{c}\approx 0.33\beta _{0}L_\perp /v_{A}$ (see (2.19)). The ‘$\beta 16$’ simulations, which do not have a matching set of passive runs, all have $\beta _{0}=16$ and scan from the collisionless to the MHD regime in order to probe the approach to collisional MHD and compare more directly with the solar wind and/or kinetic simulations of A$+$22. Finally, we have carried out a number of simulations to test additional physical effects, including the effect of finite electron temperatures and different microinstability limiters.

3.3. Diagnostics

3.3.1. Energy and rate-of-strain spectra

The spectrum of a field $\phi$ is defined as

(3.4)\begin{equation} \mathcal{E}_{\phi}(k) = \frac{1}{\delta k}\sum_{|\boldsymbol{k}|= k}|\hat{\phi}(\boldsymbol{k})|^{2}, \end{equation}

where $\hat {\phi }(\boldsymbol {k})$ is the Fourier transform of $\phi$, $|\boldsymbol {k}|= k$ indicates a sum over all modes that fall in the relevant bin of $k= \sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$ or $k_{\perp } = \sqrt {k_{y}^{2}+k_{z}^{2}}$ and $\delta k$ is the width of the bin. The spectra are computed using a fine, logarithmically spaced $k$ or $k_{\perp }$ grid, removing the large-scale bins that do not contain any modes. Wavenumbers and spectra and are plotted in units of $L_\perp ^{-1}$ and $\phi ^2 L_\perp$ (for various $\phi$), respectively, which will be omitted for conciseness in most figures.

Useful measures of turbulence structure are the spectra of the parallel and perpendicular rates of strain. These are formed via

(3.5)\begin{equation} \left.\begin{gathered} \boldsymbol{\nabla}_{\|}u_{\|} \doteq \hat{b}_{i}\hat{b}_{j}\boldsymbol{\nabla}_{i}u_{j}, \\ (\boldsymbol{\nabla}_{\|}\boldsymbol{u}_{{\perp}})_{l} \doteq \hat{b}_{i}(\delta_{lj}- \hat{b}_{j} \hat{b}_{l}) \boldsymbol{\nabla}_{i}u_{j}, \\ (\boldsymbol{\nabla}_{{\perp}}u_{\|})_{l} \doteq (\delta_{li}- \hat{b}_{i} \hat{b}_{l})\hat{b}_{j}\boldsymbol{\nabla}_{i}u_{j}, \\ (\boldsymbol{\nabla}_{{\perp}}\boldsymbol{u}_{{\perp}})_{kl} \doteq(\delta_{ki}- \hat{b}_{i} \hat{b}_{k})(\delta_{lj}- \hat{b}_{j} \hat{b}_{l}) \boldsymbol{\nabla}_{i}u_{j}, \end{gathered}\right\} \end{equation}

with the spectrum of a vector or tensor (e.g. $\boldsymbol {\nabla }_{\|}\boldsymbol {u}_{\perp }$) computed by summing the spectra of each component. This gives a $\boldsymbol {\nabla }_{\perp }\boldsymbol {u}_{\perp }$ spectrum that is nearly identical to the dissipation spectrum ($\boldsymbol {\nabla } \boldsymbol {u}$ spectrum). Note that $\boldsymbol {\nabla }_{\|}{u}_{\|}$ with this definition is exactly what must be minimised to minimise the generation of pressure anisotropy.

We also consider pressure-anisotropy gradients defined in a similar way:

(3.6)\begin{equation} \left.\begin{gathered} \boldsymbol{\nabla}_{\|}\Delta p \doteq \hat{\boldsymbol{b}} \boldsymbol{\cdot}\boldsymbol{\nabla} \Delta p, \\ \boldsymbol{\nabla}_{{\perp}}\Delta p \doteq ({{{\textsf I}}}-\hat{\boldsymbol{b}} \hat{\boldsymbol{b}} )\boldsymbol{\cdot}\boldsymbol{\nabla}\Delta p. \end{gathered}\right\} \end{equation}

These will prove useful to quantify the reduction in the pressure-anisotropy stress $\boldsymbol {\nabla }\boldsymbol {\cdot } (\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ compared with $\Delta p$ itself.

3.3.2. Structure functions

To diagnose the 3-D structure of turbulent eddies, we use three-point second-order structure functions, conditioned on the angle between the point-separation vector and the local field and/or the local perturbation. For any field $\phi$, these are defined as

(3.7)\begin{equation} S_{2}[ \phi](\boldsymbol{\ell}) = \langle | \phi(\boldsymbol{x}+\boldsymbol{\ell}) - 2 \phi(\boldsymbol{x}) + \phi(\boldsymbol{x}-\boldsymbol{\ell})|^{2}\rangle, \end{equation}

with the average taken over all $\boldsymbol {x}$.¹² The separation vector $\boldsymbol {\ell }$ is conditioned on its angle to $\boldsymbol {B}_{\ell }=[\boldsymbol {B}(\boldsymbol {x}+\boldsymbol {\ell }) + \boldsymbol {B}(\boldsymbol {x}) + \boldsymbol {B}(\boldsymbol {x}-\boldsymbol {\ell })]/3$ to obtain parallel and perpendicular structure functions ($\ell _{\|}=\hat {\boldsymbol {b}}_{\ell }\boldsymbol {\cdot }\boldsymbol {\ell }$ and $\boldsymbol {\ell }_{\perp } = \boldsymbol {\ell }-\ell _{\|}\hat {\boldsymbol {b}}_{\ell }$ with $\hat {\boldsymbol {b}}_{\ell }\doteq \boldsymbol {B}_{\ell }/|\boldsymbol {B}_{\ell }|)$ (Cho & Lazarian 2009; St-Onge et al. 2020). We also condition $\boldsymbol {\ell }$ on its direction with respect to the local field and flow perturbations to study the alignment of the turbulent fluctuations (Boldyrev 2006; Chen et al. 2012; Mallet, Schekochihin & Chandran 2015).

Following the computation of $S_{2}$, a useful diagnostic of the turbulence structure is the scale-dependent anisotropy $\ell _{\|}^{\phi }(\ell _{\perp })$ for a given field $\phi$. This is computed by solving the equation $S_{2}[ \phi ](\ell _{\perp })=S_{2}[ \phi ](\ell _{\|})$ using numerical interpolation.

3.3.3. Energy-transfer functions

Energy-transfer functions are defined as in Grete et al. (2017) and A$+$22. The transfer function $\mathcal {T}_{q\rightarrow k}^\textrm {AB}$ measures the average transfer of energy from $k$-shell $q$ to shell $k$ due to the interaction labelled AB. Here, $k$ shells are defined by the Fourier-space filtering operation

(3.8)\begin{equation} \phi_{k}(\boldsymbol{x}) = \sum_{|\boldsymbol{k}|= k}\,\hat{\phi}(\boldsymbol{k})\exp({{\rm i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{x}}), \end{equation}

where $\hat {\phi }(\boldsymbol {k})$ is the Fourier transform of some field $\phi (\boldsymbol {x})$, and $|\boldsymbol {k}|= k$ represents those wavenumbers inside the logarithmic shell centred around $k$ (i.e. $\ln k - \frac 12 \textrm {d}\ln k \leq |\boldsymbol {k}| <\ln k+\frac 12 \textrm {d}\ln k$). Such a definition clearly satisfies the property $\phi (\boldsymbol {x}) = \sum _{k}\phi _{k}(\boldsymbol {x})$ and represents the part of $\phi$ centred around wavenumber $k$. The label AB relates to the influence of different terms in (2.1)–(2.5): e.g. kinetic energy can be transferred between shells through the Reynolds stress in the momentum equation $\rho \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}$, with

(3.9)\begin{equation} \mathcal{T}_{q\rightarrow k}^{\rm UU} ={-}\int {\rm d} \boldsymbol{x} \langle \sqrt{\rho}\boldsymbol{u}\rangle_{k}\boldsymbol{\cdot}[\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \langle \sqrt{\rho}\boldsymbol{u}\rangle_{q}],\end{equation}

or magnetic energy can be transferred to kinetic energy through $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {B}$. Further details of the specific transfer terms are given in § 4.4.

The full 2-D transfer functions can give useful information on the locality of the cascade, but can be difficult to interpret quantitatively. Two useful reductions are the net energy transfer

(3.10)\begin{equation} {\rm T}^{\rm AB}(k) = \sum_{q}{\mathcal{T}}^{\rm AB}_{q\rightarrow k},\end{equation}

and the flux

(3.11)\begin{equation} \varPi^{\rm AB}(k) = \sum_{q\leq k}\sum_{p>k}{\mathcal{T}}^{\rm AB}_{q\rightarrow p}.\end{equation}

The net transfer $\textrm {T}^\textrm {AB}(k)$ is the contribution of a given term $\textrm {AB}$ in (2.1)–(2.5) to the rate of change of the energy spectrum at a particular $k$ (for example $\textrm {T}^\textrm {UU}(k)$ is the net contribution to the kinetic-energy spectrum at $k$ from the term $\rho \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}$). In the inertial range in steady state, all terms should be zero because $\partial _{t}\mathcal {E}(k)=0$, unless there is a continual transfer of energy into or out of the particular shell or between terms (for example, a damping). The flux $\varPi ^\textrm {AB}(k)$ quantifies the transfer of energy across a particular $k$, such that the sum over all contributions AB measures the cascade rate. The interpretation of individual terms is less obvious, but gives interesting information about the dominant energy-transfer processes in the cascade (e.g. whether energy proceeds to smaller scales via transfers between $\boldsymbol {u}$ and $\boldsymbol {B}$, or due to transfers from larger scale $\boldsymbol {u}$ to smaller scale $\boldsymbol {u}$ directly).

4. Results

While the most astrophysically relevant consequences of magneto-immutability concern turbulent heating, it is necessary first to understand the details of how the field and flow self-organise to be magneto-immutable and how this characteristic manifests in the turbulence statistics. We thus start by diagnosing the basic effect of pressure-anisotropy feedback by comparing PDFs and Fourier spectra of various fluctuations in the active-$\Delta p$ simulations with those obtained in the passive-$\Delta p$ simulations (§§ 4.1 and 4.2). We then consider the changes to the flow structure that are necessary in order to enable these effects in § 4.3, before diagnosing the viscous heating and cascade efficiency using energy-transfer functions in § 4.4. Overall, the results show that, aside from a small range at the outer scale, magneto-immutability allows the system to set up a vigorous, nearly conservative cascade that is in most respects similar to standard MHD.

4.1. Reduction of the pressure anisotropy

A central result – that changes to the pressure anisotropy are suppressed by its feedback on the flow (magneto-immutability) – is illustrated in figures 1–4. Figure 1 shows the joint PDF of $\beta$ and $p_{\perp }/p_{\|}$ (the ‘Brazil plot’) for all simulations from the ‘lrCL’ and ‘lrB’ sets (18 in all; see table 1). We compare the set of CGL-LF simulations (left-hand panel) with the identical passive-$\Delta p$ simulations where the pressure-anisotropy feedback has been suppressed (right-hand panel). This diagnostic, although only qualitative (the total probability is normalised separately at each $\beta$ value for illustrative purposes, but done so identically for both simulation sets), clearly demonstrating the difference between the two. In particular, we see that the force on the flow from the pressure anisotropy causes most of the plasma to remain well within the microinstability limits. In contrast, the passive case, without $\Delta p$ feedback, has most of the volume stuck at the microinstability limits (the edges of the PDFs at high and low $p_{\perp }/p_{\|}$) because the turbulent fluctuations continuously push the plasma to positive and negative $\Delta p$. Similar results are demonstrated with a snapshot of the CL10 simulation in figure 2, where we again compare the active- and passive-$\Delta p$ simulations with otherwise identical parameters. While the perpendicular flows are of similar magnitudes, indicating similar turbulent fluctuation levels (a,b), the pressure-anisotropy variation (c,d) is much smaller in the presence of pressure-anisotropy feedback. As a consequence, the variation in $B=|\boldsymbol {B}|$ is also suppressed: both $\Delta p$ and $B$ are driven by $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$, so suppressing one suppresses the other (more fundamentally, $\Delta p$ is driven by changes in $B$ and $\rho$).

Figure 2. Visualisation of $u_{z}/v_{A}$ (a,b), $4{\rm \pi} \Delta p/B^{2}$ (c,d) and $B^{2}=|\boldsymbol {B}|^{2}$ (e, f) through an $x$–$y$ slice of the CL10 simulation with $\beta _{0}=10$ and $\nu _{c}=0$ ($\mathcal {I}\approx 0.4$). (a,c,e) show the active-$\Delta p$ run (which solves (2.1)–(2.5)), and (b,d, f) show the passive-$\Delta p$ run, in which the $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ force was artificially removed. While the fluctuations in $u_{z}$ look rather similar (i.e. the turbulence has a similar amplitude), magneto-immutability significantly reduces the variation in $B$ and $\Delta p$ in the active-$\Delta p$ run. Note that $-1<4{\rm \pi} \Delta p/B^{2}<1/2$ is enforced by the limiters in both runs.

Figure 3. Left panel: PDFs of $4{\rm \pi} \Delta p/B^{2}$ in the lrB simulation set, which explores the collisionless, weakly collisional and Braginskii regimes at fixed $\mathcal {I}\approx 0.4$ (see § 3.2). Active-$\Delta p$ simulations are shown with solid lines and passive-$\Delta p$ simulations with dashed lines. The percentages shown on the left and right sides of the panel indicate the proportions (in volume) of each simulation that lie at the firehose and mirror thresholds, respectively (colours match the line styles; the larger numbers in the upper part of the plot are for the passive simulations). The effect of pressure anisotropy is relatively strong; only a few per cent of the box lies at the mirror/firehose thresholds in the active-$\Delta p$ cases, while ${\sim }50\,\%$ of the box lies beyond these thresholds when $\Delta p$ has no effect on the flow. Right panel: PDFs of $B=|\boldsymbol {B}|$ for the same simulations. In active-$\Delta p$ cases, we see a factor-of-a-few decrease in the variance of $B$ (hence, ‘magneto-immutability’), which does not depend strongly on the regime of collisionality.

Figure 4. Same as in figure 3 but for the lrCL series of collisionless simulations with varying interruption number (see § 3.2). For very small interruption numbers (e.g. in the simulation with $\beta _{0}=100$, $\nu _{c}=0$, $\mathcal {I}\simeq 0.04$), the proportion of the box at the mirror/firehose thresholds is larger (${\simeq }5\,\%$) because it is naturally driven to much larger values by the dynamics (cf. the passive simulation shown with the black dashed line). We also include a simulation with $\varLambda _\textrm {FH}=1.4$ for comparison, to explore the possible effect of kinetic oblique firehose modes limiting $\Delta p$ before the fluid firehose threshold (see § 3.1.2).

Figures 3 and 4 demonstrate the same ideas quantitatively, showing PDFs of $4{\rm \pi} \Delta p/B^{2}$ and $B$ from the lrB simulation set in figure 3 (each with similar $\mathcal {I}\approx 0.4$) and from the collisionless lrCL simulation set in figure 4 (covering $\mathcal {I}\approx 0.04$, at $\beta _{0}=100$, to $\mathcal {I}>1$, at $\beta _{0}=1$). For each of the pressure-anisotropy PDFs (left-hand panels), we indicate the volume fraction of the box that sits at the mirror or firehose thresholds as a percentage (label colours match those of the curves, and the passive-$\Delta p$ cases are the larger numbers listed in the upper part of each panel). The PDFs have qualitatively different shapes between the active- and passive-$\Delta p$ simulations, with pressure-anisotropy feedback causing $\Delta p/B^{2}$ to peak near zero in the active-$\Delta p$ runs, as opposed to being nearly flat through the stable regions in the passive runs. As a consequence, the fraction of the plasma that sits at the microinstability thresholds is strongly reduced due to the pressure-anisotropy feedback on the flow. This fraction does increase modestly with decreasing interruption number (increasing $\beta _{0}$ in figure 4), but remains very small given that $\Delta p/B^{2}$ is effectively forced across a wider range with decreasing $\mathcal {I}$. We do not see any significant change with $\nu _{c}$ across the lrB simulations (figure 3), all of which have $\mathcal {I}\approx 0.4$, showing that the collisionality regime (collisionless, weakly collisional or Braginskii MHD) is of subsidiary importance. One exception to this is the mean negative $\Delta p$ in the collisionless cases, which becomes modestly more negative with increasing $\beta$ in collisionless simulations (figure 4), but is not seen in the weakly collisional or Braginskii regimes because the bulk $\nu _{c}$ also depletes any mean $\Delta p$ that might otherwise develop. As discussed in A$+$22, this feature is likely related to turbulent heating occurring predominantly in the parallel direction via Landau damping (as approximated by the Landau-fluid closure). In figure 4, we also test the effect of a reduced firehose limit $\varLambda _\textrm {FH}=1.4$, which is likely a better approximation for scattering from kinetic oblique firehose modes at very large scale separations $L_{\perp }/\rho _{i}$ (see § 3.1.2; Bott et al. 2021). This increases the proportion of the domain at the limiter thresholds (as expected, since the range between them is smaller) and decreases the mean $|\Delta p|$. The latter effect suggests that the asymmetry of the microinstability limits (i.e. the fact that $\varLambda _\textrm {FH}>\varLambda _{M}$) also contributes to the mean pressure anisotropy, but modest changes to $\varLambda _\textrm {FH}$ do not seem to lead to any other qualitatively important changes, so we will not consider this further.

4.1.1. Turbulent energy

While the results shown in figures 1–4 clearly demonstrate the significant impact of the pressure-anisotropy feedback on the statistics of $B$ and $\Delta p$, it is important to confirm that this does not occur purely as a result of the viscous damping of those motions that would otherwise drive significant $\Delta p$. Because all simulations are driven with the same energy-injection rate $\varepsilon$, the simplest diagnostic of this is the turbulent fluctuation energy, which is shown in figure 5 for the lrCL runs (we plot $E_{K\perp } + E_{M\perp } \doteq \int \textrm {d} \boldsymbol {x}\ (\rho |\boldsymbol {u}_\perp |^2 + |\boldsymbol {B}_\perp |^2/4{\rm \pi} )/2$, which includes only the fluctuating $y$ and $z$ components of $\boldsymbol {u}$ and $\boldsymbol {B}$). A fluctuation energy that decreased with decreasing $\mathcal {I}$ would imply a cascade efficiency that decreased with $\mathcal {I}$, as would be naïvely expected if fluctuations were increasingly strongly damped above the interruption limit (see § 2.2). Instead, we see in figure 5(a) that such a decrease is quite small (compare, e.g. the active and passive cases at $\beta _{0}=10$), with a similar steady-state energy reached even at $\beta _{0}=100$ (the lowest $\mathcal {I}$ that we explore). Thus, we see that there is only modest viscous damping of pressure-anisotropy-generating fluctuations, a feature that we explore in greater depth below. Note that, in contrast, a single linearly polarised shear-Alfvén wave with similar amplitude would be strongly damped under these conditions because its fluctuations are confined to the plane set by the initial conditions and so it cannot rearrange itself to avoid generating large pressure anisotropies (Squire et al. (2016); S$+$19).

Figure 5. Time evolution of the total perpendicular energy (a) and the Alfvén ratio $r_{A}=\sqrt {\varTheta E_{M\perp }/E_{K\perp }}$ (b) for the lrCL set of collisionless simulations with varying interruption number. Note that $r_{A}$ defined in this way (with $\varTheta$ computed as a volume average) accounts for the effect of the mean pressure anisotropy on Alfvénic fluctuations (see (2.11)). The dotted lines in each panel show the ‘expected’ values assuming $\delta u_\perp /v_{A}=\delta B_\perp /B_0=1/2$: $E_{K\perp } + E_{M\perp }=V\rho _0 v_{A}^2/4$ and $r_\textrm {A}=1$. We see that the bulk properties of collisionless turbulence are rather similar to standard MHD, even for $\mathcal {I}\ll 1$.

Figure 5(b) shows the Alfvén ratio, $r_{A} = \sqrt {\varTheta E_{M\perp } /E_{K\perp }}$, for some of the simulations. This should be approximately unity for Alfvénic turbulence. The time-dependent box-averaged mean anisotropy parameter $\varTheta$ is included in $r_{A}$ because an Alfvén wave satisfies $\delta B_{\perp }=\delta u_{\perp }\sqrt {4{\rm \pi} \rho /\varTheta }$ in the presence of a mean pressure anisotropy $\Delta p$ (see § 2.1.2; we set $\varTheta =1$ for the passive simulation). Figure 5 shows that, after accounting for the change in the Alfvén ratio due to variation of the mean $\Delta p$ with $\mathcal {I}$ (see figure 4), the bulk properties of the turbulence are rather similar to MHD, with just a slight excess of magnetic energy (also observed in MHD; see, e.g. Müller & Grappin 2005; Chen et al. 2013).

4.2. Alfvénic turbulence structure

We now diagnose the structure of the Alfvénic fluctuations in more detail, with the goal of understanding the similarities and differences between turbulence in high-$\beta$, weakly collisional plasmas and turbulence in standard MHD (Schekochihin 2022). We will see that the basic statistics of the flow and magnetic field are surprisingly similar to MHD. More detailed diagnostics, which focus on compressive fields and components of the rate-of-strain tensor (see § 4.3), are needed to reveal the changes to the turbulent flow structure that enable the suppression of $\Delta p$ discussed in § 4.1.

4.2.1. Turbulent energy spectra at $\mathcal {I}\lesssim 1$

Figure 6 shows the kinetic (blue lines) and magnetic (red lines) spectra obtained in the CL10, B100 and CL100 simulations, again comparing active-$\Delta p$ simulations (solid lines) with the passive-$\Delta p$ simulations (dashed lines; note that this case is just isothermal MHD for velocity and magnetic fluctuations). The purpose of this comparison is to exhibit the extremely similar spectra in each case, demonstrating that vigorous turbulence survives even when the expected damping from pressure anisotropy is strong across a wide range of scales starting at the outer scale (as is the case for all three simulations shown here since $\mathcal {I}<1$). A careful examination reveals minor differences caused by the pressure-anisotropy feedback, in particular a slight steepening of the kinetic-energy spectrum compared with MHD. In all cases the magnetic spectra exhibit a scaling that is flatter than ${\sim }k_{\perp }^{-5/3}$ and closer to ${\sim }k_{\perp }^{-3/2}$ (shown with the dashed lines), although the exact slopes vary somewhat across the inertial range. We shall study the eddies’ 3-D structure and dynamic alignment in § 4.2.2.

Figure 6. Kinetic (blue) and magnetic (red) energy spectra for simulations CL10, B100 and CL100, as labelled. Solid lines show CGL-LF (active-$\Delta p$) simulations, and dashed lines show isothermal MHD (passive-$\Delta p$) simulations. The dotted black lines indicate slopes of $k_{\perp }^{-3/2}$, $k_{\perp }^{-5/3}$ and $k_{\perp }^{-2}$ for comparison. We see spectra broadly consistent with a $k_{\perp }^{-3/2}$ scaling, although velocity spectra are steepened modestly by the effect of pressure anisotropy in the active-$\Delta p$ runs (however, they are clearly flatter than the ${\propto }k_{\perp }^{-2}$ spectrum observed in the hybrid-kinetic simulations of A$+$22). This demonstrates that vigorous turbulence, broadly resembling MHD, is maintained even when the feedback from the pressure anisotropy is strong ($\mathcal {I}<1$ in all cases). Note that the cutoff of the spectra at relatively larger scales seen in CL100 is simply a result of that run's lower numerical resolution.

Similar information in a different form is shown in figure 7. We integrate the spectra to obtain the scale-dependent amplitudes

(4.1a,b)\begin{equation} \delta B_{{\perp}}(k_{{\perp}})^{2} = \frac{2}{V} \int_{k_{{\perp}}}^{\infty} {\rm d} k_{{\perp}}' \mathcal{E}_{M\perp}(k_{{\perp}}'),\quad \delta u_{{\perp}}(k_{{\perp}})^{2} = \frac{2}{V{\rho}_{0}} \int_{k_{{\perp}}}^{\infty} {\rm d} k_{{\perp}}'\, \mathcal{E}_{K\perp}(k_{{\perp}}'), \end{equation}

and then use these to compute the scale-dependent interruption number (see (2.19))

(4.2)\begin{equation} \mathcal{I}(k_{{\perp}}) = \left( \beta \frac{\delta B_{{\perp}}}{B_{0}}\frac{\delta u_{{\perp}}}{v_{A}}\right)^{{-}1} \max\left(1, \frac{\nu_{c}}{\omega_{A}}\right)\end{equation}

(we ignore the effect of the mean $\varTheta$ on $v_{A}$ and $\omega _{A}$, as it does not make a significant difference for these qualitative estimates). Here $\omega _{A}$ should be considered a function of scale as per critical balance, so we use $\omega _{A}(k_{\perp }) = k_{\perp } \sqrt {\delta u_{\perp }\delta B_{\perp }/{\rho }_{0}^{1/2}}$. We see that, while the estimated local interruption number is slightly larger for the CGL-LF simulations, because the turbulence amplitudes are modestly smaller, $\mathcal {I}(k_{\perp })$ nonetheless remains ${\lesssim }1$ throughout a wide portion of the turbulent cascade in each case. Since the damping rate of a linearly polarised Alfvénic perturbation is comparable to its frequency for $\mathcal {I}\lesssim 1$ (i.e. for amplitude $\delta B_{\perp }/B_{0}\gtrsim \beta ^{-1/2} \max [1,(\nu _{c}/\omega _{A})^{1/2}]$), but the turbulent fluctuations exceed this amplitude throughout most scales of the simulations, the implication is that the fluctuations must be avoiding those motions that would cause strong damping in order for the cascade to proceed as observed. Also of note is the difference between the collisionless and weakly collisional regimes, as shown by the effectively flat $\mathcal {I}(k_{\perp })$ for the B100 simulation up to $k_{\perp }\approx 100$. This occurs because when $\omega _{A}\ll \nu _{c}$ (the weakly collisional and Braginskii MHD regimes; see § 2.1.3), the increasing frequency of the motions towards smaller scales cancels their decreasing amplitudes, conspiring to keep $\omega _{A}\delta u_{\perp } \delta B_{\perp }$ approximately constant through the inertial range (see § 2.2.3). This breaks down once $\omega _{A}\gtrsim \nu _{c}$, which happens here for $k_{\perp }L_{\perp }\gtrsim 100$, around the expected scale based on the parameters (see table 1). Below this, $\mathcal {I}$ increases at smaller scales since the turbulent eddies are fast enough to be effectively collisionless.

Figure 7. ‘Local’ interruption number defined by (4.2) (see also § 2.2.3). We show the three active-$\Delta p$ simulations whose spectra are given in figure 6, comparing with the passive (MHD) runs with dashed lines. In the B100 simulation ($\beta _{0}=100$, $\nu _{c}=33v_{A}/L_{\perp }$), $\omega _{A}$ depends on scale as per critical balance, $\omega _{A}=k \sqrt {\delta u_{\perp }\delta B_{\perp }/\rho _{0}^{1/2}}$, which causes a slower increase in $\mathcal {I}$ with scale. Clearly, $\mathcal {I}(k_{\perp })\lesssim 1$ across a reasonable range in each simulation, indicating that much larger pressure anisotropies would be created without dynamic feedback from $\Delta p$ (as is the case for the passive simulations).

4.2.2. Three-dimensional anisotropy and alignment

Many previous works have studied the 3-D statistical structure of eddies in MHD turbulence, which has important implications for the cascade rate and intermittency (see Schekochihin (2022), and references therein). Starting with Boldyrev (2006), these have argued that eddies become increasingly ‘dynamically aligned’ towards smaller scales, evolving into elongated sheet-like structures satisfying $\ell _\|\gg \xi \gg \lambda$ (where $\ell _\|$, $\xi$ and $\lambda$ are the eddy's correlation lengths in the field-parallel direction, in the direction of the turbulent fluctuation, and in the mutually perpendicular direction, respectively). Motivated by its utility as a sensitive diagnostic of the Alfvénic fluctuations’ structure, in figure 8 we illustrate the 3-D anisotropy of the CL10 simulations, computed using three-point structure functions (§ 3.3.2). In the left panel, we show the directionally conditioned structure functions of $\boldsymbol {Z}^{\pm }_{\perp } = \boldsymbol {Z}^{\pm } - \hat {\boldsymbol {b}}(\hat {\boldsymbol {b}}\boldsymbol {\cdot }\boldsymbol {Z}^{\pm })$. Here, $\boldsymbol {Z}^{\pm } = \boldsymbol {u}\pm \boldsymbol {v_{A}}$ and $\boldsymbol {v}_{A}$ is the local Alfvén speed $\boldsymbol {v}_{A,\textrm {eff}}=\boldsymbol {B}\varTheta ^{1/2}(4{\rm \pi} \rho )^{-1/2}$, with $\varTheta$ and $\rho$ evaluated using their local three-point average ($\varTheta =1$ for the passive simulation). The colours show different directions of the separation vector $\boldsymbol {\ell }=(\lambda,\xi,\ell _{\|})$, with $\ell _{\|}$ applying to separations within $15^{\circ }$ of the local $\hat {\boldsymbol {b}}$, $\xi$ to separations that are perpendicular (${>}45^{\circ }$) to $\hat {\boldsymbol {b}}$ and are within $15^{\circ }$ of the direction of the local increment of $\boldsymbol {Z}^{\mp }_{\perp }$¹³ and $\lambda$ to separations that are perpendicular to both $\hat {\boldsymbol {b}}$ and the $\boldsymbol {Z}^{\mp }_{\perp }$ increment (note that we use $\boldsymbol {Z}^{\mp }$ to define the directions of $\boldsymbol {Z}^{\pm }$ because nonlinear interactions are controlled by the opposite Elsässer variable). We see results that are quite similar to standard MHD turbulence (dashed lines), with $S_{2}$ flatter in the $\lambda$ direction than in the $\xi$ direction, and broadly consistent with (though modestly steeper than) the dynamically aligned MHD-turbulence scalings $S_{2}\sim \lambda ^{1/2}$, $S_{2}\sim \xi ^{2/3}$, $S_{2}\sim \ell _{\|}^{1}$ indicated by the dotted lines (Boldyrev 2006; Mallet et al. 2015). The right panel compares the parallel–perpendicular anisotropy of eddies measured from different quantities in the turbulence. As seen in the left panel, the Alfvénic eddies are relatively similar to isothermal MHD, with $\ell _{\|}\sim \ell _{\perp }^{1/2}$ as expected for an aligned cascade. However, while in MHD the compressive quantities $u_{\|}$ and $p = T_{i}\rho$ seem to scale as $\ell _{\|}\sim \ell _{\perp }^{2/3}$, suggesting such fluctuations are unaligned and passively mixed by the critically balanced Alfvénic cascade (Lithwick & Goldreich 2001; Schekochihin et al. 2009; Chen et al. 2012), the compressive fields in CGL-LF turbulence scale quite differently, with $p_{\|}$ and $u_{\|}$ having nearly constant anisotropy $\ell _{\|}\simeq 5\ell _{\perp }$ throughout the entire box. This provides the first hint of the differences caused by magneto-immutability, which will be discussed in more detail in the next section (§ 4.3).

Figure 8. Structure functions in the $\beta _{0}=10$, $\nu _{c}=0$ simulations. The left panel shows (3-point) second-order structure functions of $\boldsymbol {Z}_{\perp }^{\pm }$ (see § 3.3.2) with increments taken along the field (blue; $\ell _{\|}$), parallel to the local fluctuation in the opposite Elsässer variable $\boldsymbol {Z}_{\perp }^{\mp }$ ($\xi$; red), and perpendicular to the local fluctuation in $\boldsymbol {Z}_{\perp }^{\mp }$ ($\lambda$; black). Solid and dashed lines show the active-$\Delta p$ and passive-$\Delta p$ cases, respectively, with all length scales normalised to $L_\perp$. The scalings are very similar in both active and passive runs, and are broadly consistent with the expected $\lambda ^{1/2}$, $\xi ^{2/3}$ and $\ell _{\|}^{1}$ scalings of dynamic alignment (shown with dotted black lines). The right-hand panel shows the measured $\ell _{\|}(\ell _{\perp })$ for the same simulation, which is computed from the second-order structure functions of a variety of different quantities as labelled (this provides a measure of the average shape of an eddy for different variables). While $\boldsymbol {Z}_{\perp }^{\pm }$ and $p_{\perp }$ show the expected scale-dependent anisotropy $\ell _{\|}\sim \ell _{\perp }^{1/2}$ of a dynamically aligned cascade, $p_{\|}$ and $u_{\|}$ have nearly constant anisotropy $\ell _{\|}\sim 5\ell _{\perp }$ throughout the inertial range. This differs from MHD (dashed lines), where $u_{\|}$ and $p$ seem to follow $\ell _{\|}\sim \ell _{\perp }^{2/3}$, as expected for non-aligned eddies in a critically balanced cascade.

4.3. The rearrangement of flows and fields in magneto-immutable turbulence

Figures 6–8 show that the CGL-LF system can sustain random, turbulent-like motions across a wide range of scales, even when those motions might be expected to be strongly damped by the effects of pressure anisotropy. The properties of the kinetic and magnetic fluctuations appear broadly similar to standard MHD, including relatively detailed measures such as the scale-dependent anisotropy and eddy-alignment intermittency. In this section, we focus on the differences compared with MHD, understanding how the turbulent cascade rearranges itself due to the pressure-anisotropy stress. We first present the numerical results, then speculate on theoretical explanations in § 4.3.3.

4.3.1. Pressure statistics

Let us first examine the spectra of pressure fluctuations, which are shown in figure 9, again comparing active- and passive-$\Delta p$ for the CL10, B100 and CL100 simulations. Here, unlike for the kinetic- and magnetic-energy spectra, we see substantial differences due to the feedback of $\Delta p$ on the flow. The most obvious feature is the significant reduction and steepening of the pressure spectra in all active-$\Delta p$ runs. While the passive-$\Delta p$ simulations have quite flat spectra – approximately ${\propto }k_{\perp }^{-1}$ in B100 and yet flatter in the collisionless simulations – all three active cases exhibit pressure spectra ${\propto }k_{\perp }^{-5/3}$ across a reasonable range near the outer scales. Notably, the latter is in agreement with the hybrid-kinetic simulations of A$+$22, providing evidence that similar physics operates in their more realistic simulations. The steepening of pressure spectra due to the pressure-anisotropy feedback is clear evidence of the modified flow structure in the active-$\Delta p$ cases, despite their similar velocity spectra.

Figure 9. Spectra of $\Delta p$ (purple lines), $p_{\perp }$ (green lines), $p_{\|}$ (yellow lines) and the magnetic pressure $B^{2}/8{\rm \pi}$ (grey lines) for the same simulations as in figure 6, and again with the passive-$\Delta p$ simulations shown with dashed lines. The dotted lines indicate the spectral scalings $\mathcal {E}\propto k_{\perp }^{-1}$ and $\mathcal {E}\propto k_{\perp }^{-5/3}$ (see § 4.3.3 for discussion). We see a significant difference between active- and passive-$\Delta p$ cases, with $\mathcal {E}_{\Delta p}\sim \mathcal {E}_{p_{\|}}\gg \mathcal {E}_{p_{\perp }}$ observed most clearly in the active-$\Delta p$ simulations (and, to a lesser degree, collisionless cases). The active-$\Delta p$ spectra bear encouraging resemblance to those seen in the hybrid-kinetic simulations of A$+$22.

Another interesting feature of the active-$\Delta p$ pressure spectra is that $\mathcal {E}_{\Delta p}\sim \mathcal {E}_{p_{\|}}\gg \mathcal {E}_{p_{\perp }}$, while in the passive runs the difference between $\mathcal {E}_{p_{\|}}$ and $\mathcal {E}_{p_{\perp }}$ is far less pronounced (particularly for B100, where $\mathcal {E}_{p_{\|}}\sim \mathcal {E}_{p_{\perp }}$, as might be naïvely expected).¹⁴ This feature was also clearly observed in the hybrid-kinetic simulations of A$+$22. As discussed in more detail in § 4.3.3 and in A$+$22, it likely results from perpendicular pressure balance between $\boldsymbol {\nabla } p_{\perp }$ and the magnetic pressure $\boldsymbol {\nabla } B^{2}/8{\rm \pi}$, while $p_{\|}$ is driven more strongly by $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$, thus dominating $\Delta p$. The basic feature seems generic to anisotropic collisionless (CGL) dynamics, occurring because $\delta p_{\perp }$ remains highly constrained by perpendicular pressure balance, while $\delta p_{\|}$ is not, affecting the plasma only through the parallel force (see Appendix A.1.2 for a simple linear calculation to illustrate the effect). We provide evidence for the scenario in figure 9 by plotting the spectra of the magnetic pressure $\mathcal {E}_{B^{2}/8{\rm \pi} }$ (grey lines), which matches $\mathcal {E}_{p_{\perp }}$ very well below the forcing scales. In the CL10 simulation, where the feedback of the pressure anisotropy is only important at larger scales because $\mathcal {I}(k_{\perp })$ increases to ${\gg }1$ at small scales (see figure 7), the feature reverses to $\mathcal {E}_{p_{\|}}\simeq \mathcal {E}_{p_{\perp }}> \mathcal {E}_{\Delta p}$ for $k_{\perp }\gtrsim 200$ (there is no similar reversal for the passive-$\Delta p$ simulation). This is presumably the signature pressure-anisotropy feedback becoming subdominant (i.e. magneto-immutability no longer being important) once the driving of $\Delta p$ via $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ becomes subdominant to standard MHD effects (the perpendicular pressure balance). We see hints that similar behaviour would occur in B100 and CL100, but at smaller scales than in CL10 (as expected), which unfortunately are unresolved in these simulations. Overall, the general similarity of the active-$\Delta p$ spectra, but not the passive-$\Delta p$ spectra, to the large-scale hybrid-kinetic results of A$+$22 is encouraging for the applicability of the CGL-LF model.

4.3.2. Rate-of-strain spectra

Given that the pressure anisotropy is driven by the plasma motions, the pressure-anisotropy feedback must be driving important changes in the flow and magnetic-field structure, despite having only a small influence on the kinetic-energy spectrum (figure 6). We diagnose this in figure 10 via the spectra of the local rates of strain, as defined in § 3.3.1. These reveal the expected substantial suppression of $\boldsymbol {\nabla }_{\|}u_{\|}=\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ by the pressure-anisotropy feedback, as needed to suppress the turbulently driven fluctuations in $\Delta p$ and $B$. Of some interest is the comparison of $\boldsymbol {\nabla }_{\|}u_{\|}$ with $\boldsymbol {\nabla }_{\perp }u_{\|}$ and $\boldsymbol {\nabla }_{\|}\boldsymbol {u}_{\perp }$, with neither of the latter two significantly suppressed. This demonstrates that the system does not reduce $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ purely via a reduction in field-parallel flows, or via the reduction of field-parallel gradients, but through the combination of the two. We see broadly similar features in the collisionless (CL10) and weakly collisional (B100) regimes, with the approximate spectral scalings $\mathcal {E}_{\boldsymbol {\nabla }_{\perp }u_{\perp }}\propto k_{\perp }^{1/3}$, $\mathcal {E}_{\boldsymbol {\nabla }_{\perp }u_{\|}}\propto k_{\perp }^{1/3}$, $\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\perp }}\propto k_{\perp }^{-1/3}$ and $\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\|}}\propto k_{\perp }^{-1/3}$ for both the active- and passive-$\Delta p$ simulations. As discussed below (§ 4.3.3), the relatively flat $\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\|}}$ scaling may be a signature of some form of cascade of compressive fluctuations, since a $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ spectrum that is dominated by linearly polarised Alfvénic fluctuations would be naïvely expected to scale as ${\sim }k_{\perp }^{-1}$. The most interesting difference between the collisionless and weakly collisional spectra is the closer-to-constant offset between the active and passive $\boldsymbol {\nabla }_{\|}u_{\|}$ spectra in the weakly collisional case (i.e. the modestly steeper $\boldsymbol {\nabla }_{\|}u_{\|}$ spectrum in B100), which may be a result of the pressure-anisotropy feedback being approximately constant across scales when $\omega _{A}<\nu _{c}$, rather than being dominated by the outer scale (see figure 7).

Figure 10. Rate-of-strain and pressure-gradient spectra in the CL10 and B100 simulations (left and right panels, respectively). Dashed lines show the equivalent passive-$\Delta p$ runs, while the colours indicate the specific strain direction considered. Details and definitions are described in § 3.3.1 (see (3.5) and (3.6)); the dotted black lines indicate the ${\sim }k_{\perp }^{1/3}$, ${\sim }k_{\perp }^{-1/3}$ and $k_{\perp }^{-1}$ scalings (see text). While $\boldsymbol {\nabla }_{\perp }u_{\perp }$, $\boldsymbol {\nabla }_{\perp }u_{\|}$ and $\boldsymbol {\nabla }_{\|}u_{\perp }$ are each relatively similar in active- and passive-$\Delta p$ runs, $\boldsymbol {\nabla }_{\|}u_{\|}=\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$, which is responsible for the creation of $\Delta p$, is markedly reduced when $\Delta p$ is active. Note that neither $\boldsymbol {\nabla }_{\|}u_{\perp }$ nor $\boldsymbol {\nabla }_{\perp }u_{\|}$ is similarly reduced, showing that the effect is not just a reduction in $u_{\|}$ or in $\boldsymbol {\nabla }_{\|}$. The bottom panels show that $\boldsymbol {\nabla }_{\|}\Delta p$ is suppressed more than $\boldsymbol {\nabla }_{\perp }\Delta p$ in active-$\Delta p$ turbulence, which is a signal that the pressure-anisotropy stress is reduced beyond just the reduction in the variance of $\Delta p$, viz., $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)/\Delta p_\textrm {rms}$ is smaller in active-$\Delta p$ than in passive-$\Delta p$ turbulence.

The lower panels in figure 10 show spectra of gradients of the pressure anisotropy. As is clear from the $\boldsymbol {u}$ evolution (2.2), it is not the pressure anisotropy itself that feeds back on the flow, but only its divergence

(4.3)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}(\hat{\boldsymbol{b}} \hat{\boldsymbol{b}} \Delta p) = \hat{\boldsymbol{b}} \boldsymbol{\nabla}_{\|}\Delta p - \hat{\boldsymbol{b}} \Delta p \boldsymbol{\nabla}_{\|}\ln B + \Delta p \hat{\boldsymbol{b}} \boldsymbol{\cdot}\boldsymbol{\nabla}\hat{\boldsymbol{b}} .\end{equation}

Computing each of these terms separately, one finds that $\boldsymbol {\nabla }_{\|}\Delta p$ is larger than the other terms, and there are no significant correlations between the different terms. Thus, if an active-$\Delta p$ run exhibits a lower ratio of $\boldsymbol {\nabla }_{\|}\Delta p$ to $\boldsymbol {\nabla }_{\perp }\Delta p$ than an equivalent passive-$\Delta p$ run, this signifies that the effect of pressure anisotropy on the flow, $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$, is reduced beyond what would be assumed by just looking at the variance of $\Delta p$ or $\boldsymbol {\nabla }_{\|}u_{\|}$. We see from figure 10 that this is indeed the case, viz. there is a secondary mechanism – the turbulence reduces $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ beyond just $\Delta p$ (or $\boldsymbol {\nabla }\Delta p$) – which will further suppress the damping of plasma motions by pressure anisotropy.

We have examined a wide variety of correlations between the different terms that drive the pressure-anisotropy evolution (2.12). For example, one might imagine that a positive correlation between $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ and $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}$ would aid in reducing $\Delta p$ generation, particularly at lower $\beta$ (see (2.13)); or similarly for correlations between the heat fluxes and other terms in the $\Delta p$ evolution (2.12). However, directly comparing such measures between the active- and passive-$\Delta p$ simulations has not yielded any further effects that are worthy of note, so we do not discuss these here.

4.3.3. Possible theoretical interpretation

In this section we speculate on possible phenomenologies that could explain the observed pressure and rate-of-strain spectra in the inertial range below the forcing scales. We examine two qualitatively different possibilities: the first is that the pressures and $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ are dominated by the residual magnetic-field variation that arises from a collection of linearly polarised Alfvénic fluctuations; the second is that they are dominated by a cascade of compressive fluctuations, which is set up around the forcing scales as the system rearranges itself. In the first scenario, the $\Delta p$ and $B^{2}$ spectra are the residual local ‘left-overs’ that cannot quite be eliminated by magneto-immutability; in the second, they are diagnosing another type of cascade or mixing of larger-scale fluctuations, similar to the passive slow-mode cascade in reduced MHD or gyrokinetics (Schekochihin et al. (2009), Kunz et al. (2015) and Meyrand et al. (2019); however, unlike these gyrokinetic cases, the compressive and Alfvénic cascades could lack individually conserved invariants, meaning the distinction between them may be less precise than implied above).¹⁵ Motivated by the observation in figure 9 that $\mathcal {E}_{p_\|}\gg \mathcal {E}_{p_\perp }$, we will imagine the compressive cascade to consist of oblique slow-magnetosonic-like modes,¹⁶ which, because they maintain perpendicular pressure balance, satisfy $\delta p_{\perp }\ll \delta p_{\|}$. This property causes the mode's frequency to scale as $k_{\|}v_\textrm {th}$, rather than $k_{\|}v_{A}$ (as for the MHD slow mode), as well as leading to significant damping by heat fluxes and/or collisions. Further details are given in Appendix A.1.2 and Majeski et al. (2023).

Let us examine predictions for the rate-of-strain and pressure spectra for each possibility in turn, starting with the idea that finite-amplitude Alfvénic fluctuations dominate. In this case, a simple estimate for the spectrum of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ comes from (2.13), using $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u} \approx (1/2) B_{0}^{-2} \,\textrm {d}(\delta B_{\perp }^{2})/\textrm {d} t \sim \omega _{A} \delta B_{\perp }^{2}/B_{0}^{2}$. For a critically balanced Goldreich–Sridhar cascade ($\delta B_{\perp }\propto k_{\perp }^{-1/3}$, $\omega _{A}\propto k_{\perp }^{2/3}$), this suggests that $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u} \sim k_{\perp }^{0}$, or a $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ spectrum $\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\|}}\propto k_{\perp }^{-1}$, which is clearly significantly steeper than observed (cf. the steepest dotted line in figure 10). However, this prediction is less rigorous than it seems because ignores the influence of magneto-immutability: larger $\delta B_{\perp }$, with larger local $\mathcal {I}(k_{\perp })$, will be more affected by $\Delta p$ forces, thus reducing $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ below $\omega _{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 (figure 9). 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 figure 9). Because $\delta p_{\|}$ is instead driven by $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ to $\delta p_{\|} \sim \beta \max (1,\omega _{A}/\nu _{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 Appendix A.1.2), a far larger ratio than what is observed in all three simulations shown in figure 9.¹⁷ 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 {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ should scale as ${\sim }k_{\|}u_{\|}\propto k_{\perp }^{1/3}$, giving $\mathcal {E}_{\boldsymbol {\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 ($\boldsymbol {\nabla }_{\|}u_{\perp }\propto k_\perp ^{1/3}$, $\boldsymbol {\nabla }_{\perp }u_{\|}\propto k_\perp ^{2/3}$ and $\boldsymbol {\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 {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {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 {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {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 with 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 transferred 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}^\textrm {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}_\textrm {AB}$ to refer to them below. We follow Grete et al. (2017) in our nomenclature; $\textrm {AB}=\textrm {UU}$ refers to the transfer between shells of the kinetic energy through $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u} + \mathcal {C}_\textrm {UU}$ (see (3.9)), $\textrm {AB}=\textrm {BU}$ refers to the transfer from magnetic to kinetic energy through $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {B} + \mathcal {C}_\textrm {BU}$, $\textrm {AB}=\textrm {UB}$ refers to the transfer from kinetic to magnetic energy through $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}+\mathcal {C}_\textrm {UB}$, $\textrm {AB}=\textrm {BB}$ refers to the transfer between shells of the magnetic energy through $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {B}+\mathcal {C}_\textrm {BB}$, $\textrm {AB}={ \Delta p U}$ refers to the transfer from the kinetic to thermal energy via the pressure anisotropy through $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ (see below) and $\textrm {AB}=\textrm {F U}$ refers to the transfer to kinetic energy from the external forcing.¹⁸ 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 with Grete et al. (2017) is that due to the pressure-anisotropy stress $\mathcal {T}_{q\rightarrow k}^{ \Delta p U}$. There is some freedom in the definition of this; A$+$22 use

(4.4)\begin{equation} \mathcal{T}_{q\rightarrow k}^{\Delta p\mathrm{U}} = \int\mathrm{d}\boldsymbol{x} \left\langle\sqrt{\rho}\boldsymbol{u}\right\rangle_{k} \boldsymbol{\cdot} {\rm sign}({\Delta p}) \frac{\sqrt{|\Delta p|} }{B}\frac{\boldsymbol{B}}{\sqrt{4{\rm \pi} \rho}} \boldsymbol{\cdot} \boldsymbol{\nabla}\langle\sqrt{|\Delta p|}\hat{\boldsymbol{b}}\rangle_{q},\end{equation}

which has the advantage that the square of $\langle \sqrt {|\Delta p|}\hat {\boldsymbol {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

(4.5)\begin{equation} \mathcal{T}_{q\rightarrow k}^{\Delta p\mathrm{U}} = \int\mathrm{d}\boldsymbol{x} \left\langle\sqrt{\rho}\boldsymbol{u}\right\rangle_{k} \boldsymbol{\cdot} \frac{\boldsymbol{B}}{\sqrt{4{\rm \pi} \rho}} \boldsymbol{\cdot} \boldsymbol{\nabla}\left\langle \frac{\Delta p}{B^{2}}{\boldsymbol{B}}\right\rangle_{q}.\end{equation}

We have computed both versions, finding qualitatively similar results, but will focus on (4.5), because it fits more naturally with our focus on pressure-anisotropy damping. Note also that the influence of the sign-discontinuity issue of (4.4) 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).

In figure 11, we plot the net energy transfers $\textrm {T}^\textrm {AB}(k)$ (see (3.10)), 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, $\textrm {T}^\textrm {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 $\textrm {T}^\textrm {AB}(k)$ by the dimensional (Kolmogorov) estimate $\partial _{t}\mathcal {E}\sim k u_{k}\mathcal {E}\sim \varepsilon (k/k_{0})^{-1}$ (where $k_{0}=2{\rm \pi} /L_{\perp }$); with this normalisation, a line that is constant and non-zero 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 with the results from the corresponding passive-$\Delta p$ simulations (lower subpanels). The thick black lines, which show the sum of all terms (excluding $\textrm {T}^{\Delta p 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 $\textrm {T}^{\Delta p U}$ is slightly negative in both CL10 and B100, indicating a slight damping of inertial-range motions, it is small compared with the cascade rate. Contrast this to the equivalent $\textrm {T}^{\Delta p 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 $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ can be both positive and negative (it is negative definite in the Braginskii-MHD regime, but not otherwise; see (2.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.

Figure 11. Important terms in the net energy-transfer spectra, as defined in § 3.3.3 and (3.10). The top panels show the simulations with $\beta _{0}=10$, $\nu _{c}=0$ and the bottom panels show the simulations with $\beta _{0}=100$, $\nu _{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 labelled in the top panel: $\mathrm {T}^{\Delta p U}$ is the transfer from $\boldsymbol {u}$ to thermal energy due to $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$; $\mathrm {T}^\textrm {UU}$ is momentum transfer via $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}$; $\mathrm {T}^\textrm {BU}$ is the contribution to kinetic energy from magnetic tension and pressure; $\mathrm {T}^\textrm {BB}$ is the advection of $\boldsymbol {B}$ in the induction equation; $\mathrm {T}^\textrm {UB}$ is the contribution to magnetic energy from field stretching; and $\mathrm {T}^\textrm {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 shows similar information in a different form by plotting the contributions to the turbulent energy flux (3.11). We show only the active-$\Delta p$ runs because the comparison with passive simulations is less interesting in this case. Unlike the net transfers (figure 11), individual terms in the fluxes are not easily interpreted because they do not include the diagonal $\mathcal {T}^\textrm {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 $\varPi ^{\Delta p 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 figure 11 that there is little pressure-anisotropy damping through the inertial range (a slight decrease in $\varPi (k)$ is observed for B100, indicating a slight damping). The value of the $\varPi /\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)

Figure 12. The energy fluxes (see (3.11)), which measure the transfer of energy across a particular $k$, for the two active-$\Delta p$ simulations with $\beta _{0}=10$, $\nu _{c}=0$ (a) and $\beta _{0}=100$, $\nu _{c}=33v_{A}/L_{\perp }$ (b). We normalise to the forcing input $\varepsilon$, which implies that the total flux would be $\varPi /\varepsilon =1$ across the full inertial range for a conservative cascade with no damping. As also seen in figure 11, 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.

4.4.2. Pressure-anisotropy heating fraction and cascade efficiency

The net energy transfer $\textrm {T}^{\Delta p 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 (figure 11), 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, Stone & Quataert 2016; Kempski et al. 2019). Nonetheless, its dependence on parameters ($\beta$ and the collisionality) is interesting to consider, as is a comparison with A$+$22.

We define

(4.6)\begin{equation} \mathcal{H}_{\Delta p} = \int_{k_{0}}^{k_{\rm max}} {\rm d} k \,{\rm T}^{\Delta p U}(k),\end{equation}

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_\textrm {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_\textrm {max}$ to lie near the bottom of the inertial range. We set $k_\textrm {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.¹⁹

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 _{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 _{c}$ (otherwise $\mathcal {H}_{\Delta p}$ would continue increasing below $\mathcal {I}\simeq 10$).

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 (figure 11) by summing up the contributions from all $k\leq 100$, viz. $\mathcal {H}_{\Delta p}=\sum _{k\leq 100} \mathrm {T}^{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 _{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 _{c}/\omega _{A} \gtrsim 1$), and less in collisionless cases. This value is approximately constant for $\mathcal {I}\lesssim 10$ and independent of $\beta$ and $\nu _{c}/\omega _{A}$ (so long as $\nu _{c}/\omega _{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.

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 _{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, Uzdensky & Kunz 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.²⁰ 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 with 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 with 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 with 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 with 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 _{c}^\textrm {lim}$ to be very large in unstable regions; however, kinetic theory and simulations show that the saturated scattering rate should actually be $\nu _{c}\sim S\beta$, where $S\sim \hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {u}$ is the shear driven by the turbulence. The effect on $\Delta p$ is similar – both hard-wall limiters and scattering with $\nu _{c}\sim S\beta$ act to pin $\Delta p$ at the stability thresholds – but the differences could lead to inaccuracies.²¹ Furthermore, kinetic simulations have shown that the assumption that microinstabilities saturate and decay instantaneously compared with 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. 2017a; 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 _{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.²² Comparing the microinstability-scattering estimate above, $\nu _{c}\sim \beta S \sim \beta \omega _{A} \mathcal {M}_{A}^{2}$, with the Alfvénic interruption estimates discussed in § 2.2, we see that the $\nu _{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 _{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 _{c}/\omega _{A}\sim \mathcal {M}_{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{\rm \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 _{c}\simeq \beta S$ at $\beta =16$, or $\mathcal {I}\simeq 1$ (see their figure 12).²³ 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 _{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 _{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 with 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 with 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 _{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, Chen & Squire 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, Klein & Maruca 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 $|\boldsymbol {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

(5.1a--c)\begin{equation} \frac{\delta\boldsymbol{B}}{\sqrt{4{\rm \pi} \rho}} ={\pm}\boldsymbol{u},\quad B={\rm const},\quad p_{{\perp}},p_{\|},\rho={\rm const}\end{equation}

is an exact nonlinear solution of the Kinetic MHD system (2.1)–(2.5) (or even the standard MHD equations). Any such solution propagates at velocity $\mp \boldsymbol {B}_{0}\sqrt {\varTheta /4{\rm \pi} \rho }$, where $\boldsymbol {B}_{0}$ is the mean of $\boldsymbol {B}$ and $\delta \boldsymbol {B}$ is the remainder (Barnes & Hollweg 1974). Since the solution (5.1a–c) holds even if $|\delta \boldsymbol {B}|\gg |\boldsymbol {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, Chandran & Meyrand 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 (5.1a–c) (Dunn et al. 2023). But in order for such states to occur, the turbulence must be strongly imbalanced, with $|\boldsymbol {z}_{+}|=|\boldsymbol {u}+\boldsymbol {B}/\sqrt {4{\rm \pi} \rho } |\gg |\boldsymbol {z}_{-}|=|\boldsymbol {u}-\boldsymbol {B}/\sqrt {4{\rm \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 $|\boldsymbol {z}_{+}| \sim |\boldsymbol {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$.

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

(5.2)\begin{equation} C_{B^{2}}=\frac{(\delta B^{2})_{\rm rms}}{ (\delta \boldsymbol{B}_{\rm rms})^{2}}= \frac{\langle(|\boldsymbol{B}|^{2}-\langle|\boldsymbol{B}|^{2}\rangle)^{2}\rangle^{1/2}}{\langle|\boldsymbol{B}- \langle\boldsymbol{B}\rangle|^{2}\rangle},\end{equation}

which was introduced by Squire et al. (2020) as a simple measure of the relative tendency of the components of $\boldsymbol {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$).²⁴ For reference, Parker Solar Probe 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 with 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 with 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.

Figure 15. The variation in magnetic-field strength, or spherical polarisation, for all simulations, quantified in terms of $C_{B^{2}}$ given by (5.2). This measures the relative tendency for the components of $\boldsymbol {B}$ to be correlated in order to reduce the variation in $|\boldsymbol {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 (5.1a–c).

5.2.2. Spectral diagnostics

In § 4, we saw that, while kinetic- and magnetic-energy spectra remained rather similar in magneto-immutable turbulence compared with 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 _{c}$) to large-$\mathcal {I}$ collisional MHD, while other parameters ($\mathcal {M}_{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 _{c}/\omega _{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 lrB600 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 figure 9); in contrast, when unaffected by the pressure-anisotropy feedback, the $p_{\perp }$ and $p_{\|}$ spectra have similar shapes.

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 lrB600, $\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}_{\boldsymbol {\nabla }_{\|}u_{\perp }}/\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\|}}$ is shown with solid lines, and $\mathcal {E}_{\boldsymbol {\nabla }_{\perp }u_{\perp }}/\mathcal {E}_{\boldsymbol {\nabla }_{\perp |}u_{\|}}$ is shown with dashed lines (see figure 10 and § 3.3.1). The increase in $\boldsymbol {\nabla }_{\|}u_{\perp }/\boldsymbol {\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 $\boldsymbol {\nabla }_{\perp }u_{\perp }/\boldsymbol {\nabla }_{\perp |}u_{\|}$ allows the effect to be considered separately from $u_{\|}$ itself, which may depend on the forcing and other parameters.

Since we saw that the reduction in $\boldsymbol {\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}_{\boldsymbol {\nabla }_{\|}u_{\perp }}/\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\|}}$ (solid lines) with those of perpendicular-flow gradients $\mathcal {E}_{\boldsymbol {\nabla }_{\perp }u_{\perp }}/\mathcal {E}_{\boldsymbol {\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 _{c}/\omega _{A}=200$, $\mathcal {I}\approx 50$), there remain clear differences in the flow structure compared with 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.²⁵ While the exact values of these metrics will depend on details (e.g. the forcing, $\mathcal {M}_{A}$, and the imbalance), we expect that the qualitative features should be robust, especially the fact that $\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\perp }}/\mathcal {E}_{\boldsymbol {\nabla }_{\|}u_{\|}}\gg \mathcal {E}_{\boldsymbol {\nabla }_{\perp }u_{\perp }}/\mathcal {E}_{\boldsymbol {\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 figure 9).

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-LF 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 the collisionless plasma dynamics on scales far above the proton gyroradius. Indeed, our results compare well with 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 $-\boldsymbol {\nabla } p$ force that rapidly opposes any flow with $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {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 $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {\boldsymbol {b}} \hat {\boldsymbol {b}} \Delta p)$ that rapidly opposes any flow with $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {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 with inertial forces, the ‘interruption number’ $\mathcal {I}$ (introduced in S$+$19) quantifies the importance of pressure-anisotropic forces compared with 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 with 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 with‘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 with 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 {\boldsymbol {b}}\hat {\boldsymbol {b}}:\boldsymbol {\nabla } \boldsymbol {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 }$ (figure 9). 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 figures 11 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 with 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 with 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 microinstabilities (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.