2.1 The dynamical equations

Because we wish to apply our theory to systems in which the ion-Larmor radii are as much as ten orders of magnitude smaller than the collisional mean free paths, we choose to employ a collisionless model that neglects finite-Larmor-radius effects, saving micro-scale considerations like the mirror and firehose instabilities for the numerical analysis of § 3. The two most suited models for describing a plasma under these assumptions are then the drift-kinetic-MHD model of Kulsrud (Reference Kulsrud1983) and the CGL-MHD fluid approach of Chew et al. (Reference Chew, Goldberger and Low1956). While the drift-kinetic approach is more accurate in that it self-consistently determines the heat fluxes via moments of the guiding-centre distribution function, we instead employ the simpler CGL-MHD model. This, of course, places some limits on which applications of magneto-immutability can be studied using our simulation technique. For example, if investigating the manner in which cosmic rays are accelerated or diffused by the turbulent field and flow at large scales, then the CGL-MHD model may suffice. On the other hand, if one wishes to study the shape of the plasma distribution function beyond the fluid moments, or probe the microphysics of dissipation at sub-ion-Larmor scales, then the approach we take here will be insufficient. Additionally, CGL-MHD model by itself provides no closure for the form of the heat fluxes, and so for lack of anything markedly better we adopt the ‘Landau-fluid’ closure introduced by Snyder, Hammett & Dorland (Reference Snyder, Hammett and Dorland1997), which is designed to capture the effects of linear Landau damping on the fluid quantities. Fortunately, we find that the turbulence behaves in such a way that the exact form of the heat fluxes is relatively unimportant, suggesting that any errors introduced by applying these approximate heat fluxes are small (indeed, in § 3.6 we demonstrate that magneto-immutability occurs in a hybrid-kinetic simulation of high-$\beta$ Alfvénic turbulence). With that in mind, the CGL-MHD equations are given by (Chew et al. Reference Chew, Goldberger and Low1956)

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

(2.1b)\begin{gather} \rho \frac{{\rm d} \boldsymbol{u}}{{\rm d} t}={-}\boldsymbol{\nabla}\left(p_{{\perp}}+\frac{B^2}{8\pi}\right)+\boldsymbol{\nabla} \boldsymbol{\cdot}\left[\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}\left(p_{{\perp}} - p_\parallel{+} \frac{B^2}{4\pi}\right)\right] , \end{gather}

(2.1c)\begin{gather} \frac{{\rm d} \boldsymbol{B}}{{\rm d} t}= \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{B} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}, \end{gather}

(2.1d)\begin{gather} p_{{\perp}} \frac{{\rm d} }{{\rm d} t} \ln \frac{p_{{\perp}}}{\rho B} ={-} \boldsymbol{\nabla}\boldsymbol{\cdot}(q_{{\perp}}\hat{\boldsymbol{b}}) - q_{{\perp}}\boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{b}} , \end{gather}

(2.1e)\begin{gather} p_\parallel \frac{{\rm d} }{{\rm d} t}\ln \frac{p_\parallel B^2}{\rho^3} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}(q_\parallel\hat{\boldsymbol{b}}) + 2q_{{\perp}}\boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{b}} , \end{gather}

where $\boldsymbol {B}$ is the magnetic field, $\boldsymbol {u}$ is the ion flow velocity and $\rho$ is the mass density. Note that $p_{\perp }$ and $p_\parallel$ are defined with respect to the local magnetic-field direction $\hat {\boldsymbol {b}}=\boldsymbol {B}/B$. Here, ${\rm d}/{\rm d} t\doteq \partial /\partial t + \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the convective time derivative. Equations (2.1d) and (2.1e) demonstrate how, in the absence of heat fluxes, the quantities $p_{\perp }/\rho B$ and $p_{\|} B^2/\rho ^3$ are conserved in time along the flow of the plasma. Conservation of $p_{\perp }/\rho B$ and $p_{\|} B^2/\rho ^3$ is the collective result of individual particles conserving both their magnetic moment $\mu \doteq mw_{\perp }^2/2B$ and their parallel action (i.e. bounce invariant) $\mathcal {J} = \oint \mathrm {d}w_{\|} \,mw_{\|}$, respectively (where $\boldsymbol {w}\doteq \boldsymbol {v}-\boldsymbol {u}$ is the particle velocity peculiar to the fluid frame). Because of the plasma's strong magnetization, the flows of perpendicular-/parallel-thermal energy, $q_{\perp /\|}$, occur exclusively along the local magnetic-field direction. For these quantities we adopt the ‘3+1 model’ of Snyder et al. (Reference Snyder, Hammett and Dorland1997)

(2.2a)\begin{gather} q_{{\perp}}={-}\frac{v_{{\rm th} \|}}{\sqrt{\pi} \left|k_{\|}\right|}\left[\rho \boldsymbol{\nabla}_{\|}\left(\frac{p_{{\perp}}}{\rho}\right)-p_{{\perp}}\left(1-\frac{p_{{\perp}}}{p_{\|}}\right) \frac{\boldsymbol{\nabla}_{\|} B}{B}\right], \end{gather}

(2.2b)\begin{gather} q_{\|}={-}\frac{2v_{{\rm th} \|}}{\sqrt{\pi} \left|k_{\|}\right|} \rho \boldsymbol{\nabla}_{\|}\left(\frac{p_{\|}}{\rho}\right), \end{gather}

where $v_{{\rm th},\|}=\sqrt {2 p_{\|}/\rho }$ is the parallel-thermal speed, and $\boldsymbol {\nabla }_{\|} = \hat {\boldsymbol {b}} \boldsymbol {\cdot } \boldsymbol {\nabla }$ is the field-aligned gradient. The field-aligned wavenumber $|k_{\|}|$ in the denominators of (2.2) is meant to be representative of a characteristic parallel scale of the perturbations to $\rho$, $B$ and $p_{\perp /\|}$ (e.g. Sharma et al. Reference Sharma, Hammett, Quataert and Stone2006). This quantity is a stand in for the magnitude of the magnetic-field-aligned gradient operator, which is generally difficult to calculate as it must be evaluated along the exact perturbed fields at each time step. For the majority of the numerical simulations presented in § 3, we take $|k_\parallel |=4\pi /L_\parallel$ where $L_\parallel$ is the field-parallel outer scale; for the remainder of this section, however, the exact value of $|k_\parallel |$ is unimportant. Because the heat fluxes (2.2) are designed to capture only linear collisionless damping, it is implicitly assumed that the perturbations being studied are small enough with respect to the background plasma that nonlinear damping effects can be ignored. We also assume that the background plasma pressure is isotropic, and that the electrons are cold. Effects of a background anisotropy on inertial-range kinetic turbulence have been considered by Kunz et al. (Reference Kunz, Schekochihin, Chen, Abel and Cowley2015) and could be straightforwardly implemented within our model, however, they are not necessary for describing magneto-immutability and assessing its robustness. Similarly, electron pressure is not essential to the conclusions we will draw; as long as the electrons are sufficiently collisional and isothermal, a finite electron pressure has no qualitative effect on our results.Footnote ¹ Finally, the CGL-MHD equations’ inability to resolve $r_{L,{\rm i}}$-scales inevitably implies that the growth of the mirror and firehose instabilities and their interaction with the plasma particles are not properly captured. In our numerical simulations (§ 3), we model their effect on the plasma through an anomalous collisionality that isotropizes the pressure to marginally unstable values, but no effort is made to incorporate the effects of this micro-instability scattering into our analytic model for magneto-immutability.

2.2 Relevant properties of collisionless hydromagnetic waves

Developing an intuition of the linear modes and their interactions is particularly beneficial towards understanding the scale-by-scale transfer of energy in a turbulent MHD cascade. It is therefore important to review the linear modes of the Landau-fluid CGL system, so that we can make assumptions and formulate our ordering in light of these fundamental behaviours. Many of these waves are discussed in greater detail within Appendix B of Majeski et al. (Reference Majeski, Kunz and Squire2023); here, we restrict ourselves to aspects of the linear modes that are relevant to collisionless high-$\beta$ turbulence.

Our discussion begins with those modes that are least affected by the transition from collisional to collisionless regimes. In particular, linear collisionless Alfvén waves are identical to their collisional counterparts when their wavelengths far exceed $r_{L,{\rm i}}$, because neither density nor pressure perturbations are generated as they propagate. Meanwhile, the mode that experiences the most minor non-zero change is the fast mode, which has a nearly unchanged phase speed, although it is susceptible to collisionless damping depending on the propagation angle with respect to the background field, ($k_{\|}/k_{\perp }$). Consequently, most of the differences between high-$\beta$ CGL and collisional MHD turbulence do not originate from the behaviour of Alfvén or fast waves. The remaining collisionless hydromagnetic modes are non-propagating modes, ion-acoustic waves and kinetic entropy modes, which are most easily compared with the slow-magnetosonic and pressure-balanced entropy modes of collisional MHD.Footnote ² Because there is one more dynamical equation in CGL-MHD than in collisional MHD, there is one additional linear mode solution. However, the kinetic entropy mode is both static and heavily damped at high $\beta$, with no clear collisional counterpart (see, e.g. Majeski et al. Reference Majeski, Kunz and Squire2023), so it can be excluded from this comparison.

We are then left with the one-to-one comparisons of non-propagating modes with pressure-balanced entropy modes, and ion-acoustic waves with slow-magnetosonic waves. Both the non-propagating and pressure-balanced entropy modes have zero phase speed, with the entropy mode being fully static, and the non-propagating mode decaying due to transit-time damping at a rate $\gamma \sim |k_{\|}|v_{A}/\sqrt {\beta }$ that is small at high $\beta$. And while entropy modes satisfy perfect isotropic-pressure balance (allowing them to remain static), non-propagating modes exist in a state of approximate pressure balance between $p_{\perp }$ and $B^2$ when $k_{\|} \ll k_{\perp }$ (Majeski et al. Reference Majeski, Kunz and Squire2023). For these reasons, at high $\beta$ these modes are rather similar, with the important feature that they are both approximately static. The comparison of ion-acoustic and slow waves, on the other hand, yields far fewer similarities, especially at high $\beta$. The most essential difference lies in their dispersion relations, with the ion-acoustic wave having a frequency $\omega \sim k_{\|}v_{{\rm th}}$, while the slow-wave frequency scales as $k_{\|} v_{A}$ when $\beta \gg 1$ and $\omega \sim k_{\|}v_{\rm th}$ when $\beta \ll 1$. The ion-acoustic wave is also Landau damped at a rate $\gamma \sim |k_{\|}|v_{\rm th}$, meaning that at high $\beta$ all characteristic time scales of ion-acoustic waves are much shorter than those of slow (and Alfvén) waves. These disparate time scales are fundamental to how collisionless waves interact with an Alfvénic cascade. For example, in collisional MHD turbulence, a time scale separation exists between fast and Alfvén waves because fast waves have a phase speed proportional to $k$ rather than just $k_{\|}$ (with $k_{\|} \ll k$ in the inertial range of Alfvénic turbulence, see Goldreic & Sridhar Reference Goldreic and Sridhar1995). As a result, it is often argued that the rapid propagation of fast waves decouples them from any Alfvénic dynamics. In the next subsection (§ 2.3), we will argue that at $\beta \gg 1$, the ion-acoustic wave – the only other compressive ($\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} \neq 0$) wave – decouples from the Alfvénic dynamics as well because $v_{\rm th} \gg v_{A}$.

2.3 An asymptotic ordering for high-$\beta$ collisionless turbulence

Equations (2.1) are too complicated with which to work analytically, therefore we seek an asymptotic ordering that distils the physics responsible for magneto-immutability. Before doing so, however, it is instructive first to consider the asymptotic ordering used in the collisional reduced MHD (RMHD) employed by Zank & Matthaeus (Reference Zank and Matthaeus1992) and Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009, § 2), from which our ordering borrows heavily

(2.3)\begin{equation} \frac{u_{{\perp}}}{v_{A}} \sim \frac{u_{\|}}{v_{A}} \sim \frac{\delta B_{{\perp}}}{B_0} \sim \frac{\delta B_{\|}}{B_0} \sim \frac{\delta \rho}{\rho_0} \sim \frac{\delta p}{p_0} \sim \frac{k_{\|}}{k_{{\perp}}}\doteq \epsilon \ll 1 \sim \beta, \quad \partial_t \sim k_{\|} v_{A}. \end{equation}

Here, all of the perturbations belong to a given $k_{\perp /\|}$ ‘cone’ and interactions are taken to be local in $k$-space, meaning that, as perturbations become smaller through the cascade, so does the ratio $k_{\|}/k_{\perp }$. Arguably the most essential aspect of this ordering is that it enforces ‘critical balance’ of the Alfvénic fluctuations, viz. $k_{\perp } u_{\perp } \sim k_{\|} v_{A}$. This is a statement that the nonlinear time scale associated with eddy deformation by the Reynolds stress $(k_{\perp } u_{\perp })^{-1}$ is comparable to the linear propagation time scale of Alfvén waves $(k_{\|} v_{A})^{-1}$, as is expected for strong turbulence (Goldreic & Sridhar Reference Goldreic and Sridhar1995). As a result, the evolutionary time scale at each wavenumber is expected to be roughly $\partial _t \sim k_{\|} v_{A}$, since linear and nonlinear time scales are equivalent. This is also true of MHD slow modes, given that they propagate at $v_{A}$ or slower, allowing for the formation of a similarly strong compressive cascade. As mentioned in § 2.2, however, the enhanced propagation speed of fast modes prevents them from being effectively coupled to the Alfvénic motions, and as such they are ordered out of RMHD. Our ordering shares most of these assumptions, yet differs in a few key ways.

First, we incorporate the largeness of the background $\beta$ directly into the ordering, so that $\beta \sim \epsilon ^{-1}$. This may seem odd because $\beta$, being a background quantity, is constant across all scales even though $\epsilon$ becomes smaller as the fluctuations cascade anisotropically to larger $k$. However, we are concerned with plasmas in which perturbations near the outer scale often exceed $\beta ^{-1}$ in relative amplitude. As a result, it would be incorrect to eliminate terms of order $(u_{\perp }/v_{A})^2$ but keep those of order $u_{\perp }/\beta v_{A}$. Ordering $\beta \sim \epsilon ^{-1}$ allows us to retain both terms. We also assume that the Alfvénic component of our turbulence is critically balanced, as in collisional RMHD. This is not necessarily a given, however, there is evidence from the work of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) that immutability does not interfere with critical balance, a feature that we reproduce through simulations of our own in § 3. Next, we assume that, as a result of the plasma being collisionless, $\delta p$ is replaced in the ordering by $\delta p_{\perp }$ and $\delta p_{\|}$, with similar amplitudes to (2.3). This is a straightforward assumption, yet its implications are more nuanced. Doing so replaces slow-magnetosonic waves with ion-acoustic waves, and because we include $\beta$ in the overall ordering, this implies that nonlinear mixing of ion-acoustic waves by Alfvén waves is weak when local in $k$-space, since $k_{\|} v_{\rm th} \gg k_{\perp } u_{\perp }$. We will therefore assume that ion-acoustic waves, like fast waves, decouple from the Alfvénic turbulence. This assumption is supported by analytical calculations and simulations of the interaction of Alfvén and ion-acoustic waves in high $\beta$ plasmas, using the same CGL-MHD code employed in this work (Majeski & Kunz Reference Majeski and Kunz2024). Ordering out ion-acoustic waves therefore enables us to continue with the assertion that, for all variables, time derivatives are of order $\partial _t \sim k_{\|} v_{A}$.

Our final major assumption, which follows from those before, is that density fluctuations can be neglected in determining the leading-order dynamics (i.e. $\delta \rho /\rho _0 \sim \epsilon ^2$ at most). In Alfvénic turbulence, high values of $\beta$ naturally inhibit the magnitude of density fluctuations by reducing the sonic Mach number ($u/v_{\rm th}$) of the forcing. Furthermore, any forcing must be associated with time scales that are sonic or faster if it is to produce substantial compressive fluctuations, a rather rare occurrence in the systems with which we are concerned. As a result, density fluctuations are expected to be small in amplitude at the outer scale. We are able to apply this assumption throughout the inertial range because the only linearly compressive waves – ion-acoustic and fast modes – have been ordered out of the dynamics. This assumption of small $\delta \rho$ allows us to consider only perturbations to the temperatures $\delta T_{\perp /\||}$ rather than the pressures $\delta p_{\perp /\|}$.

The collisionless, high-$\beta$ ordering that results from the above considerations and assumptions is

(2.4)\begin{equation} \frac{u_{{\perp}}}{v_{A}} \sim \frac{u_{\|}}{v_{A}} \sim \frac{\delta B_{{\perp}}}{B_0} \sim \frac{\delta B_{\|}}{B_0} \sim \frac{\delta T_{{\perp}}}{T_0} \sim \frac{\delta T_{\|}}{T_0} \sim \frac{1}{\beta} \sim \frac{k_{\|}}{k_{{\perp}}} \doteq \epsilon , \quad \partial_t \sim k_{\|} v_{A}. \end{equation}

Note that the choice $\delta T_{\perp /\|}/T_0 \sim \beta ^{-1}$ means that $\beta \varDelta \sim 1$, or in other words, the anisotropic-pressure stress is present at the same order as the flow inertia.Footnote ³ This is important, because, in order to demonstrate that magneto-immutability avoids significant $\Delta p$ stress, $\Delta p$ must first be made large enough to disrupt the turbulence in the absence of immutability.Footnote ⁴ Some of the assumptions that yield (2.4), such as $\delta \rho /\rho _0$ being negligible, are only justified qualitatively, therefore we test them against numerical simulations in § 3 to establish confidence in the relevance of the ordering. We are now prepared to apply the ordering to (2.1) and obtain the high-$\beta$ reduced CGL-MHD equations.

2.4 The high-$\beta$ reduced CGL-MHD equations

In this section we formulate equations that describe the leading-order evolution for each of the fields present in the ordering (2.4). Doing so involves expanding these quantities in (fractional) orders of $\epsilon$, for example with $T_{\|}$

(2.5)\begin{equation} T_{\|} = T_0 + \delta T_{\|}^{(1)} + \delta T_{\|}^{(3/2)} + \delta T_{\|}^{(2)} + \cdots, \end{equation}

where the parenthetical superscript represents the order of each term in powers of $\epsilon$. Note that, beyond the leading-order perturbation, which is of order $\epsilon$, half-integer orders must be used for all higher-order perturbations. This is required by the heat fluxes, which are proportional to $v_{\rm th} \sim v_{A}/\sqrt {\epsilon }$.

Keeping in mind the expectation that density fluctuations are $\mathcal {O}(\epsilon ^2)$ at most and thus do not affect the leading-order dynamics, the continuity equation (2.1a) simplifies to ${\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0}$. If ordered according to (2.4), this becomes

(2.6)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = \boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}} + \boldsymbol{\nabla}_{\|}u_{\|} \approx \boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}}^{(1)} = 0. \end{equation}

The divergence-free condition for $\boldsymbol {B}$ similarly yields $\boldsymbol {\nabla }_{\perp } \boldsymbol {\cdot } \delta \boldsymbol {B}_{\perp }^{(1)} = 0$ to leading order. These conditions on $\boldsymbol {u}_{\perp }^{(1)}$ and $\delta \boldsymbol {B}_{\perp }^{(1)}$ assert that the perpendicular flow and magnetic-field perturbations are dominated by Alfvén waves, which are naturally incompressible. Given that large-scale collisionless Alfvén waves are linearly identical to those of collisional plasmas, it is no surprise that this aspect of the reduced equations is unchanged from standard RMHD. Continuing with the induction equation (2.1c), the leading-order parallel and perpendicular components become

(2.7a)\begin{gather} \frac{{\rm d} \delta B_\parallel^{(1)}}{{\rm d} t^{(0)}} = B_0\hat{\boldsymbol{b}}^{(0)} \boldsymbol{\cdot} \boldsymbol{\nabla} u_{\|}^{(1)} , \end{gather}

(2.7b)\begin{gather} \frac{{\rm d} \delta \boldsymbol{B}_{{\perp}}^{(1)}}{{\rm d} t^{(0)}} = B_0\hat{\boldsymbol{b}}^{(0)} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{{\perp}}^{(1)} , \end{gather}

where

(2.8a)\begin{equation} \frac{{\rm d} }{{\rm d} t^{(0)}} \doteq \frac{\partial }{\partial t} + \boldsymbol{u}_{{\perp}}^{(1)} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\perp}} \quad \mathrm{and} \quad \hat{\boldsymbol{b}}^{(0)} \doteq \hat{\boldsymbol{z}} + \frac{\delta \boldsymbol{B}_{{\perp}}^{(1)}}{B_0} , \end{equation}

with $\hat {\boldsymbol {z}}$ being the direction of the background magnetic field $\boldsymbol {B}_0$. Once again, these equations are equivalent to those obtained for the evolution of the magnetic field in standard RMHD. Equation (2.7a) states that leading-order changes to the magnetic-field strength are generated through field-aligned shear in $u_{\|}$, which we will see is reduced in magneto-immutable turbulence. That this has not been ordered out directly by (2.4) hints that magneto-immutability must be achieved through self-organization.

Next, we address the momentum equation. Beginning with the perpendicular direction, the first two orders are trivially

(2.9a)\begin{equation} \boldsymbol{\nabla}_{{\perp}} \delta T_{{\perp}}^{(1)} = 0 \quad \mathrm{and} \quad \boldsymbol{\nabla}_{{\perp}} \delta T_{{\perp}}^{(3/2)} = 0 , \end{equation}

as a result of $\beta$ being large. Given that $k_{\|} \ll k_{\perp }$ at all scales, these perturbations can have no parallel gradients either, thus $\delta T_{\perp } ^{(1)} = \delta T_{\perp }^{(3/2)} = 0$. The next two orders give perpendicular pressure balance

(2.10a)\begin{equation} \frac{\delta T_{{\perp}}^{(2)}}{T_0} + \frac{\delta \rho^{(2)}}{\rho_0} ={-}\frac{2}{\beta} \frac{\delta B_{\|}^{(1)}}{B_0} \quad \mathrm{and} \quad \frac{\delta T_{{\perp}}^{(5/2)}}{T_0} + \frac{\delta \rho^{(5/2)}}{\rho_0}={-}\frac{2}{\beta}\frac{\delta B_{\|}^{(3/2)}}{B_0}, \end{equation}

analogous to that of RMHD and observed within the simulations of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). Note that, although these equations involve the higher-order density perturbations, $\delta \rho ^{(2,5/2)}$ need not be determined independently as we will show that the perpendicular temperature perturbations $\delta T_{\perp }^{(2,5/2)}$ do not enter into the leading-order equations. The subsequent order of the momentum equation dictates the evolution of $\boldsymbol {u}_{\perp }$

(2.11)\begin{align} \rho_0 \frac{{\rm d} \boldsymbol{u}_{{\perp}}^{(1)}}{{\rm d} t^{(0)}} & ={-}\boldsymbol{\nabla}_{{\perp}} P_{\rm total}^{(3)} + \frac{B_0}{4\pi}\frac{\partial \delta \boldsymbol{B}_{{\perp}}^{(1)}}{\partial z} + \frac{1}{4\pi}\delta \boldsymbol{B}_{{\perp}}^{(1)} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\perp}} \delta \boldsymbol{B}_{{\perp}}^{(1)} \nonumber\\ & \quad -\frac{\beta}{8\pi} B_0\frac{\partial}{\partial z} \left(\delta \boldsymbol{B}_{{\perp}}^{(1)} \frac{\delta T_{\|}^{(1)}}{T_0} \right) - \frac{\beta}{8\pi}\delta \boldsymbol{B}_{{\perp}}^{(1)} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\perp}} \left(\delta \boldsymbol{B}_{{\perp}}^{(1)} \frac{\delta T_{\|}^{(1)}}{T_0} \right). \end{align}

Here, $P_{\rm total}^{(3)}$ represents the combined thermal and magnetic pressures evaluated at third order, which can be determined as a whole by enforcing $\boldsymbol {\nabla }_{\perp } \boldsymbol {\cdot } \boldsymbol {u}_{\perp } = 0$. Because of perpendicular pressure balance (2.9a,b), the leading-order pressure anisotropy is dominated by the contribution from $\delta T_{\|}^{(1)}/T_0$. The top line of (2.11) contains those terms that are already present in RMHD, while the bottom line incorporates feedback from the pressure anisotropy onto the perpendicular momentum. The contribution from pressure anisotropy remains because our ordering assumes ${\beta \varDelta \sim 1}$, although this holds different meanings at different scales. Near the outer scale, the anisotropy is dominated by fluctuations of amplitude $\varDelta \sim \beta ^{-1}$. As the cascade continues to smaller scales, however, fluctuations in $\varDelta$ become smaller in amplitude and no longer compete with the magnetic tension, thus dropping out of the leading order of (2.11). As a result, the pressure anisotropy that shows up in (2.11) is really only that of the largest scales, which the high-$k$ Alfvén modes feel as a nearly constant modified background $v_{A,{\rm eff}}$.

It is tempting to identify the pressure-anisotropy-related terms in (2.11) with magneto-immutability, however, (2.7a) suggests that magneto-immutability will primarily be mediated by $u_{\|}$. Therefore, the parallel component of (2.1b) must be examined, for which the leading two orders are simply

(2.12)\begin{equation} \hat{\boldsymbol{b}}^{(0)} \boldsymbol{\cdot} \boldsymbol{\nabla} \delta T_{\|}^{(1)} = \left(\hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \delta T_{\|}\right)^{(1)} = 0, \end{equation}

and

(2.13)\begin{equation} \hat{\boldsymbol{b}}^{(0)} \boldsymbol{\cdot} \boldsymbol{\nabla} \delta T_{\|}^{(3/2)} + \frac{\delta \boldsymbol{B}_{{\perp}}^{(3/2)}}{B_0} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\perp}} \delta T_{\|}^{(1)} = \left(\hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \delta T_{\|}\right)^{(3/2)} = 0. \end{equation}

These equations result directly from ordering out ion-acoustic waves. If the ordering were such that $\partial _t u_{\|} \sim \epsilon k_{\|} v_{\rm th}^2$, then the inertial term would be of the same order as the field-aligned temperature gradient (2.12), which would then be non-zero. Linearization of the resultant reduced equations would yield an ion-acoustic-like eigenfrequency proportional to $k_{\|}v_{\rm th}$, however, in our reduced equations no such eigenfrequency can be obtained. This is analogous to pressure balance in the perpendicular equations, which orders out fast modes by virtue of the assumption that the time derivative of $u_{\perp }$ is not proportional to $kv_{\rm th}$. Unlike (2.9a,b), (2.12) and (2.13) involve both parallel and perpendicular gradients and alignment with $\delta \boldsymbol {B}_{\perp }$. Hence, it is not necessarily true that $\delta T_{\|}^{(1)} = \delta T_{\|}^{(3/2)} = 0$ as is the case for $\delta T_{\perp }^{(1)}$ and $\delta T_{\perp }^{(3/2)}$, but only that the field-aligned gradients of $\delta T_{\|}$ are negligible. Also of note is that we have yet to make any mention of the heat fluxes, implying that this minimal variation of $\delta T_{\|}$ along field lines is dynamic, rather than diffusive.Footnote ⁵

For our purposes, there is no need to obtain the higher-order contributions to the parallel momentum equation, and so we continue applying our ordering by examining the evolution equations for the double adiabats. The first step in doing so is to consider the heat fluxes (2.2) in light of (2.9a,b), (2.12) and (2.13), which can be rewritten in the following more transparent manner:

(2.14a)\begin{gather} q_{{\perp}}={-}\frac{v_{{\rm th} \|}\rho }{\sqrt{\pi} \left|k_{\|}\right|}\left( \hat{\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla} T_{{\perp}} + \frac{\beta_{{\perp}}\varDelta}{2} \frac{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B}{4\pi\rho}\right), \end{gather}

(2.14b)\begin{gather} q_{\|}={-}\frac{2v_{{\rm th} \|}\rho}{\sqrt{\pi} \left|k_{\|}\right|} \left(\hat{\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla} T_{\|}\right). \end{gather}

Here, $\beta _{\perp } \doteq 8\pi p_{\perp }/B^2$ and $\beta _{\perp } \varDelta \sim 1$, so the right-most term in (2.14a) is smaller than $p_{\perp } \mathrm {d}_t \ln (T_{\perp }/B)$ by roughly $\epsilon ^{3/2}$, and can therefore be neglected. What remains are heat fluxes proportional to the field-aligned temperature gradients $\hat {\boldsymbol {b}} \boldsymbol {\cdot } \boldsymbol {\nabla } T_{\perp /\|}$. However, these gradients are zero not just to leading order, but also to next order thanks to (2.9a,b), (2.12) and (2.13). As a result, $q_{\perp /\|}$ can be neglected entirely to the leading two orders at which they would contribute to (2.1d) and (2.1e). Thus, the surprisingly simple outcome is that

(2.15)\begin{equation} \frac{{\rm d} \delta B_{\|}^{(1)}}{{\rm d} t^{(0)}} = \frac{{\rm d} \delta T_{\|}^{(1)}}{{\rm d} t^{(0)}} = 0. \end{equation}

This result is essentially the namesake of magneto-immutability. In order to satisfy both $\mu$ conservation and perpendicular pressure balance in a high-$\beta$ plasma, the flow must self-organize to enforce $(\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla } \boldsymbol {u})^{(1)}=0$, so that (2.7a) and (2.15) are in agreement. As a result, the magnetic field only experiences convective changes in its strength to leading order.Footnote ⁶ In describing magneto-immutability qualitatively, Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) also drew several comparisons with incompressibility in collisional MHD turbulence, and here we may draw another more direct comparison. In a collisional high-$\beta$ plasma, pressure balance requires the isotropic pressure perturbation $\delta p$ be reduced in order to match the pressure of the magnetic perturbation $\delta B_{\|}$. With conservation of the single adiabat ${\rm d}_t(p/\rho ^\gamma )=0$, the smallness of $\delta p$ then requires that ${\rm d}_t \rho ^{-\gamma }={\rm d}_t \delta \rho =0$. Thus both magneto-immutability and incompressibility in high-$\beta$ magnetized turbulence result from a combination of pressure balance and adiabatic invariance, only here the involvement of $B$ in the collisionless invariants necessitates magneto-immutability.

Provided proper initial and boundary conditions, we have successfully closed the system of equations. Collecting the full high-$\beta$, reduced CGL-MHD equations (and dropping the orders), we have

(2.16a)\begin{gather} \boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}} = \boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \delta \boldsymbol{B}_{{\perp}} = 0, \end{gather}

(2.16b)\begin{gather} \frac{{\rm d} \delta \boldsymbol{B}_{{\perp}}}{{\rm d} t} = B_0\hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{{\perp}}, \end{gather}

(2.16c)\begin{gather} \rho_0 \frac{{\rm d} \boldsymbol{u}_{{\perp}}}{{\rm d} t} ={-}\boldsymbol{\nabla}_{{\perp}} P_{\rm total} + \left( 1+\frac{\beta}{2}\frac{\Delta p}{p_0} \right)\left(\frac{B_0 \hat{\boldsymbol{z}} + \delta \boldsymbol{B}_{{\perp}} }{4\pi}\right) \boldsymbol{\cdot} \boldsymbol{\nabla} \delta \boldsymbol{B}_{{\perp}}, \end{gather}

(2.16d)\begin{gather} \frac{{\rm d} \Delta p}{{\rm d} t} = \frac{{\rm d} \delta B_{\|}}{{\rm d} t} = \hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} u_{\|} = \hat{\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \Delta p = 0 . \end{gather}

In (2.16d) and (2.16c), we have replaced $\delta T_{\|}/T_0$ with $-\Delta p/p_0$ without consequence, because $\Delta p/p_0 = \delta p_{\|}/p_0$ up to order $\epsilon ^{3/2}$. A homogeneous background pressure anisotropy $\Delta p_0$ can be straightforwardly included in these equations by appending it to the fluctuating $\Delta p$ in (2.16c). As in collisional RMHD, the divergence of $\boldsymbol {u}_{\perp }$ can be employed to determine $P_{\rm total}$, and the convective derivative and field-parallel gradient become

(2.17a)\begin{equation} \frac{{\rm d} }{{\rm d} t} = \frac{\partial }{\partial t} + \boldsymbol{u}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\perp}} \quad \mathrm{and} \quad \hat{\boldsymbol{b}} = \hat{\boldsymbol{z}} + \frac{\delta \boldsymbol{B}_{{\perp}}}{B_0}, \end{equation}

respectively. Note that the reduced order of $\boldsymbol {\nabla }_{\|} \delta T_{\|}$ compared with $\delta T_{\|}$ itself allowed the anisotropic pressure to be pulled outside of the gradient operator in (2.16c) (since $\boldsymbol {\nabla }_{\|} \Delta p$ becomes next order as well). In this sense, the pressure anisotropy is felt only as a modification to the effective Alfvén speed, thereby hindering its ability to interfere with the Alfvénic cascade as an anisotropic-pressure stress that could cause turbulent motions to damp into thermal energy. This is a central characteristic of magneto-immutability, and one of significant consequence. Without it, the turbulent cascade would cease long before reaching kinetic scales, dramatically modifying the transport properties of the plasma.

2.5 Features of the reduced system

We now discuss some properties of the reduced system (2.16), beginning with the quantities that its turbulent cascade conserves. In the same manner as Kunz et al. (Reference Kunz, Schekochihin, Chen, Abel and Cowley2015), this analysis is made easier by defining Elsässer variables that incorporate the modification to the Alfvén speed by the pressure anisotropy (Elsässer Reference Elsässer1950)

(2.18)\begin{equation} \boldsymbol{z}^{{\pm}} = \boldsymbol{u}_{{\perp}} \pm v_{A,{\rm eff}} \frac{\delta \boldsymbol{B}_{{\perp}}}{B_0}. \end{equation}

However, unlike in Kunz et al. (Reference Kunz, Schekochihin, Chen, Abel and Cowley2015), here, $v_{A,{\rm eff}}$ is understood to be a function of both position and time. It is a straightforward process to show that the reduced system (2.16) requires the cascades of Elsässer energies to obey

(2.19)\begin{equation} \frac{\partial W_{\rm AW}^{{\pm}}}{\partial t} \doteq \frac{\partial }{\partial t} \int\, \mathrm{d}^3 x \, \frac{1}{2} |\boldsymbol{z}^\pm|^2 ={\mp} \int\, \mathrm{d}^3 x \, \frac{P_{\rm total}}{\rho_0} \frac{\partial v_{A,{\rm eff}}}{\partial z}. \end{equation}

This is notably distinct from the Elsässer cascades of Alfvénic energy in collisional MHD and reduced kinetic-MHD (RKMHD), which satisfy $\partial _t W_{\rm AW}^\pm = 0$ individually (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Kunz et al. Reference Kunz, Schekochihin, Chen, Abel and Cowley2015). Here, only the total Alfvénic fluctuation energy $W_{\rm AW}^++W_{\rm AW}^-$ is conserved. Meanwhile, the cross-helicity $W_{\rm AW}^+-W_{\rm AW}^-$, which measures the imbalance in forward-propagating and backward-propagating Alfvén waves, is not conserved. The fact that $W_{\rm AW}^++W_{\rm AW}^-$ is conserved proves that the pressure anisotropy, having reduced-order parallel gradients, cannot effectively thermalize the energy contained within Alfvénic fluctuations (Squire et al. Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023, also see figure 7 here). That being said, the non-conservation of cross-helicity means that $\Delta p$ is still able to redistribute energy between forward- and backward-propagating waves (Majeski & Kunz Reference Majeski and Kunz2024). The other conserved quantities of the magneto-immutable cascade are rather trivial to detect from (2.16)

(2.20a)\begin{equation} \frac{\partial }{\partial t} \int\, \mathrm{d}^3x \,\frac{B_0 \delta B_{\|}}{4\pi} = 0 \quad \mathrm{and} \quad \frac{\partial }{\partial t} \int\, \mathrm{d}^3x\, \delta p_{\|} = 0, \end{equation}

where we have chosen to configure the constancy of $\delta B_{\|}$ in units of energy for consistency with the other conserved quantities. Equations (2.20a,b) represent, respectively, cascades of non-propagating modes and kinetic entropy modes. The reason for the respective associations of each conserved quantity comes down to the eigenvectors of each mode. Non-propagating modes, being in a state of near-perpendicular pressure balance, have $\delta B_{\|}/B_0 \gg \delta p_{\perp }/p_0$ at high $\beta$. However, they also exhibit $\delta p_{\perp } < \delta p_{\|}$, which is inconsistent with the fact that $\delta p_{\|} \gg \delta p_{\perp }$ in magneto-immutable turbulence. Thus, the passive cascade of $\delta p_{\|}$ can most easily be attributed to kinetic entropy modes, which vanish due to collisional damping in the MHD limit (Majeski et al. Reference Majeski, Kunz and Squire2023). To leading order, there is no conserved quantity for compressive fluctuations (namely those associated with $u_{\|}$), because ion-acoustic and fast waves have been ordered out of the dynamics.

To further highlight the features that make magneto-immutable turbulence unique, we can compare the reduced (2.16) with the standard RMHD system, which is obtained by applying (2.3) to the collisional MHD equations. The dynamics of the Alfvénic fluctuations $\delta \boldsymbol {B}_{\perp }$ and $\boldsymbol {u}_{\perp }$ is almost entirely unchanged from (2.16a)–(2.16c), with the only difference being that $\Delta p=0$ in collisional MHD (Zank & Matthaeus Reference Zank and Matthaeus1992; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). This suggests that the differences between the two models originate in the compressive cascade, a component of the turbulence that plays a passive role in collisional RMHD. The equations of $\delta B_{\|}$ and $u_{\|}$ in RMHD turbulence are (Schekochihin & Cowley Reference Schekochihin and Cowley2007)

(2.21a)\begin{gather} \frac{{\rm d} }{{\rm d} t}\frac{\delta B_{\|}}{B_0} = \frac{1}{1+2/(\gamma\beta)} \,\hat{\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla} u_{\|} \approx \hat{\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla} u_{\|} \quad (\beta \gg 1), \end{gather}

(2.21b)\begin{gather} \frac{{\rm d} }{{\rm d} t}\frac{u_{\|}}{v_{A}} = v_{A} \hat{\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla} \frac{\delta B_{\|}}{B_0}, \end{gather}

where $\gamma =5/3$ is the single-adiabatic index for a monatomic gas. These equations describe slow-magnetosonic waves, which are well known to comprise the passively advected compressive cascade in MHD turbulence, something that is clearly not present in (2.16). As discussed in the context of conserved quantities, the non-Alfvénic cascades can be attributed to non-propagating and kinetic entropy modes, with no truly compressive ($\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} \neq 0$) cascade of any kind present to leading order.

When constructing the ordering (2.4) under the assumption $\delta B_{\|}/B_0 \sim \beta ^{-1}$, we essentially guaranteed that realistic turbulent cascades eventually pass out of the parameter space within which (2.16) are strictly valid. Yet, these equations are nonetheless accurate even when $\delta B_{\|}/B_0 \ll \beta ^{-1}$ (so long as $\beta \gg 1$). To prove this, in Appendix A we show that the signatures of magneto-immutability, viz. $\boldsymbol {\nabla }_{\|}\Delta p$ and $\boldsymbol {\nabla }_{\|} u_{\|}$ suppression, can be recovered via a subsidiary high-$\beta$ ordering of the RKMHD equations, which assume that $\delta B_{\|}/B_0 \ll \beta ^{-1}$. Indeed, with the assumption of Landau-fluid heat fluxes we once again obtain the reduced system (2.16).Footnote ⁷ Note, however, that at these small scales where $\delta B_{\|}/B_0 \ll \beta ^{-1}$, such a lack of parallel gradients in $u_{\|}$ and $\Delta p$ need only be enforced passively, as they have already been suppressed by larger scales obeying the ordering (2.4). In this sense, magneto-immutability exists in a somewhat weaker, watchdog-like state. If $k_{\|}$ were somehow introduced to $u_{\|}$ or $\Delta p$ at small scales – for instance, by field-line wandering or magnetic reconnection (Meyrand et al. Reference Meyrand, Kanekar, Dorland and Schekochihin2019) – it would be suppressed by magneto-immutable self-organization. However, if these parallel gradients were not created at small scales, they simply would remain negligible with $\Delta p$ cascading passively. To contrast this with collisional RMHD or $\beta \sim 1$ RKMHD, any re-introduced parallel gradients in $u_{\|}$ and $\Delta p$ would be allowed to persist in the absence of dissipative effects. Unlike $u_{\|}$, the evolution of small-scale $u_{\perp }$ in high-$\beta$ collisionless turbulence becomes quite similar to that of collisional RMHD. The $u_{\perp }$ fluctuations only feel a non-local modification to $v_{A,{\rm eff}}$ produced by large scale $\Delta p$ satisfying $\beta \varDelta \sim 1$, thus $u_{\perp }$ and $\delta B_{\perp }$ effectively decouple from the $u_{\|}$ and $\delta B_{\|}$ cascades. Furthermore, the gradients in $v_{A,{\rm eff}}$ are negligible at these small scales, allowing the energetic coupling between $|\boldsymbol {z}^+|^2$ and $|\boldsymbol {z}^-|^2$ that appears in (2.19) to weaken. Thus, the Alfvénic cascade approaches that of RKMHD with a constant background anisotropy, as described by Kunz et al. (Reference Kunz, Schekochihin, Chen, Abel and Cowley2015).

Equations (2.16) leave us with several new testable predictions. First and foremost is that magneto-immutability, as a reduction of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla } \boldsymbol {u}$ and $\hat {\boldsymbol {b}} \boldsymbol {\cdot }\boldsymbol {\nabla } \Delta p$, is independent of scale. In the initial studies of Squire et al. (Reference Squire, Schekochihin, Quataert and Kunz2019) and Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023), the authors concluded that the pressure anisotropy being driven by turbulent fluctuations must be competitive with the background magnetic tension (i.e. $\beta \varDelta \sim 1$) in order for the suppression of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla } \boldsymbol {u}$ and $\hat {\boldsymbol {b}} \boldsymbol {\cdot }\boldsymbol {\nabla } \Delta p$ to occur. However, both of these signatures of magneto-immutability persist throughout the cascade simply as a result of $\beta \gg 1$, even when $\beta \Delta _k \ll 1$ at some large wavenumber $k$. Note, however, that because $\delta p_{\perp /\|}$ and $\delta B_{\|}$ are passively advected, the outer scale is important in that the compressive perturbations must be seeded with sufficiently large amplitude there.Footnote ⁸ The second testable prediction is that magneto-immutability is, to leading order, independent of the heat fluxes. This is owed to the dynamical, rather than diffusive, manner in which minimal field-parallel variation of $\delta T_{\|}$ is achieved, forcing the heat fluxes to contribute only to the next-order pressure anisotropy. Both magneto-immutability's independence of the heat fluxes, and the field-parallel spreading of $\delta T_{\|}$ in the absence of heat fluxes, can therefore be validated directly by well-constructed CGL-MHD simulations. Equation (2.15) also predicts that $\delta T_{\|}$ and $\Delta p$ behave as passive scalars, thereby adopting the statistics of the flow that nonlinearly mixes them (Biskamp Reference Biskamp2003). Although not addressed analytically, it is also implied that a sufficiently large scattering frequency in the pressure equations could disable magneto-immutability. It is essential that the $\delta p_{\perp }$ and $\delta p_{\|}$ evolve independently enough that $\delta p_{\perp }$ can be suppressed by pressure balance to give $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u} \approx 0$, while $\hat {\boldsymbol {b}}\boldsymbol {\cdot }\boldsymbol {\nabla } \delta p_{\|}$ is suppressed by the $u_{\|}$ momentum equation. This arrangement can still be preserved in the presence of a small amount of scattering ($\nu \ll k_{\|}v_{A}$), as the scattering term would remain next order in the $p_{\perp /\|}$ equations and the resultant reduced system would be unchanged. However, a scattering rate of the order of or larger than $k_{\|}v_{A}$ would force the pressures to evolve together, and the two magneto-immutability criteria would not be able to be met independently (see Appendix B for more). This point must be borne in mind when considering the effects of microinstabilities, because if their volume-filling fraction is large enough, then the consequent scattering may cause the plasma to become MHD-like.

3.1 Problem set-up and method of solution

To assess the claims set forth in § 2, we perform a suite of driven CGL-MHD turbulence simulations with a variety of parameters tuned to address each individual prediction or assumption. We employ a modified version of the Athena++ MHD code (Stone et al. Reference Stone, Tomida, White and Felker2020) that solves the CGL-MHD system (2.1) closed with the Landau-fluid heat fluxes (2.2). This code allows the parallel wavenumber $|k_{\|}|$ in the Landau-fluid heat fluxes (2.2) to be specified freely, and incorporates the effects of mirror and firehose instabilities (when excited) through a collisional closure. Unless capped by some maximal value, the Landau-fluid $q_{\perp /\|}$ could grow very large at small scales because of the numerical simplification that $|k_{\|}|$ is chosen to be constant. To prevent this from occurring, the heat fluxes are not allowed to surpass a maximal ‘free-streaming’ value of ${\approx }v_{\rm th}p_{\perp /\|}$ (Hollweg Reference Hollweg1974; Cowie & McKee Reference Cowie and McKee1977). Exact details of how this limitation is implemented in the code can be found in § 3.1.1 of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). The collisional microinstability closure uses a limiting scattering frequency $\nu _{\rm lim}$ (also specified by the user), which activates only in regions of the domain that exceed the mirror ($\beta \varDelta > 1$) or firehose ($\beta \varDelta < -2$) thresholds (e.g. Sharma et al. Reference Sharma, Hammett, Quataert and Stone2006). Once activated, the pressures $p_{\perp /\|}$ are driven back towards the instability thresholds, rather than to complete pressure isotropy, at a rate set by $\nu _{\rm lim}$. The ability to choose $|k_{\|}|$ and $\nu _{\rm lim}$ provides considerable freedom to explore how various collisionless effects change the behaviour of turbulence in this system. Further details about the CGL-MHD solver and microinstability closure can be found in Appendix A of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). As first described in § 1.3, we supplement these ‘active-$\varDelta$’ CGL-MHD simulations with a set of ‘passive-$\varDelta$’ simulations, which are performed in isothermal MHD but evolve the pressure anisotropy passively using (2.1d) and (2.1e) given the simulated MHD fields. Such passive-$\varDelta$ simulations are useful for comparing MHD-like turbulence with CGL-MHD turbulence.

Given that the particles’ Larmor scales are infinitesimally small in our model equations, the physical dimensions of our simulations are arbitrary. All simulations are performed in a fully periodic domain with dimensions $[L_x,L_y,L_z] = [1,1,2]$, where $L_z = L_{\|}$ is aligned with the background magnetic field $\boldsymbol {B}_0 = \hat {\boldsymbol {z}}$ and $L_x = L_y = L_{\perp }$ are the perpendicular dimensions. The ‘standard’ resolutions employed for these simulations are $n_{\perp } = 192$ and $n_{\|} = 2n_{\perp }$, however, higher- and lower-resolution simulations are performed (and explicitly referred to) for the sake of convergence and scale-dependence tests (see figure 11). Although the initial magnetic field $\boldsymbol {B}_0$ remains the same in all simulations, the initially isotropic thermal pressure $p_0 = \beta _0 B_0^2/8\pi$ is varied by choosing $\beta _0$ to be either 1, 10 or 100. For $|k_{\|}|$ and $\nu _{\rm lim}$, we make use of a set of ‘standard’ values in which $|k_{\|}| = 4\pi /L_{\|}$ represents the wavenumber of compressive fluctuations near the outer scale and $\nu _{\rm lim} = 10^{10} v_{A}/L_{\perp }$ yields a hard-wall limiter that prevents the pressure anisotropy from straying far beyond its microinstability thresholds. Specific instances in which $|k_{\|}|$ and $\nu _{\rm lim}$ are modified from their standard values are noted on a case-by-case basis. All simulations are run until a final time of at least $t_{f} = 10 L_{\perp }/v_{A}$ to ensure that a steady-state fluctuation level is achieved within the turbulence.

The turbulence is forced exclusively through sinusoidal perturbations to the flow velocity $\boldsymbol {u}$ using an Ornstein–Uhlenbeck correlated process (Uhlenbeck & Ornstein Reference Uhlenbeck and Ornstein1930), the strength of which is input numerically as the rate of change of the total kinetic energy in the domain $\mathrm {d}_t \mathcal {E}_{K}$. Our fiducial runs employ $\mathrm {d}_t \mathcal {E}_{K} = 0.32\rho _0 v_{A}^2L_{\perp }^3$ and a correlation time of $t_{\rm corr} = L_{\|}/v_{A}$. This choice of $t_{\rm corr}$ assumes Alfvénically correlated forcing, and the energy injection rate corresponds to an outer-scale magnetic perturbation amplitude of $\delta B_{\perp /\|} \approx B_0/2$ in steady state. For all simulations, the sinusoidal wavenumbers at which we force are limited to $k \in 2\pi /L_{\|} \times [1,3]$, over which the power distribution scales as $k^{-2}$. In this study, we frequently vary the mode of $\boldsymbol {u}$ forcing between Alfvénic and random forcing. Random forcing is as it sounds: $\boldsymbol {u}$ is perturbed randomly in all directions without any special conditions relating to the directions of $\boldsymbol {B}$ and $\boldsymbol {k}$ (although it is still time correlated). In Alfvénic forcing, however, we perturb only $\boldsymbol {u} \perp \hat {\boldsymbol {z}}$, and do so in a manner that enforces incompressibility, $\boldsymbol {\nabla }_{\perp } \boldsymbol {\cdot } \boldsymbol {u}_{\perp } = 0$, at the outer scale (below the outer scale, however, $\boldsymbol {u}$ does become slightly compressible as a result of the nonlinear amplitudes).

3.2 Numerical diagnostics

To analyse the simulations introduced in § 3.1, we make use of several numerical diagnostics, three of which we describe here because they are either used very frequently in our analysis or are particularly tailored to the subject of this work. These three are the Fourier spectra, the $\Delta p$ energy transfer function of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) and a novel scale-by-scale $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}$ alignment diagnostic based on the work of St-Onge et al. (Reference St-Onge, Kunz, Squire and Schekochihin2020).

The most frequently used diagnostic is a bin-averaged Fourier spectrum, defined for a field $\chi$ as

(3.1)\begin{equation} E_\chi(k) = \frac{1}{\delta k} \sum_{|\boldsymbol{k}|=k} |\mathcal{F} [\chi](\boldsymbol{k}) |^2 , \end{equation}

where $\mathcal {F} [\chi ](\boldsymbol {k})$ is the three-dimensional Fourier transform of $\chi$, and $k=\sqrt {k_x^2+k_y^2+k_z^2}$ is the magnitude of $\boldsymbol {k}$ (or $k_{\perp } = \sqrt {k_x^2+k_y^2}$ for $\boldsymbol {k}_{\perp }$). The spectra are calculated for $k$ bins, thus each bin contains contributions from several wavenumbers in the Fourier spectrum, and $\delta k$ represents the bin width in $k$ or $k_{\perp }$. Note that, in figure 6(b), when calculating the spectra of $\boldsymbol {\nabla }_{\perp /\|}u_{\perp /\|}$, the gradients of $\boldsymbol {u}$ are calculated with respect to the local (in configuration space) field, rather than using the same efficiency-motivated $\boldsymbol {k}_{\perp } = k_x\hat {\boldsymbol {x}} + k_y\hat {\boldsymbol {y}}$ assumption that the spectra employ.

Energy transfer functions illustrate the amount of turbulent kinetic energy transferred into or out of a given $k$-shell as a result of specific interaction terms in the model equations (e.g. Grete et al. Reference Grete, O'Shea, Beckwith, Schmidt and Christlieb2017). In this work, we are chiefly interested in the effects of pressure anisotropy on the turbulent cascade, thus the transfer function we make use of is designed to capture the energy removed from the flow $\boldsymbol {u}$ at a given $k_{\perp }$ due to the anisotropic pressure-stress $\boldsymbol {\nabla } \boldsymbol {\cdot } (\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\Delta p)$. We use the same definition given in Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023), where

(3.2)\begin{equation} \mathcal{T}_{\Delta p}(k_{{\perp}}) = \sum_{q_{{\perp}}} \int \,\mathrm{d}^3 \boldsymbol{x} \langle \sqrt{\rho}\boldsymbol{u} \rangle_{k_{{\perp}}} \boldsymbol{\cdot} \frac{\boldsymbol{B}}{\sqrt{4\pi\rho}} \boldsymbol{\cdot} \boldsymbol{\nabla} \left\langle \frac{\Delta p}{B^2}\boldsymbol{B} \right\rangle_{q_{{\perp}}}, \end{equation}

and $\langle \cdot \rangle _k$ indicates that the quantity has been Fourier transformed, filtered by wavenumber and returned to real space. There are other definitions for this transfer function (see Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023), however, as remarked by Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023), this particular formulation represents the $\Delta p$-stress as a damping of kinetic energy. This is well suited to the current study as we wish to understand how magneto-immutability permits the continuation of a turbulent cascade that would have otherwise been damped away near the viscous scale.

Finally, the $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}$ alignment diagnostic is a method for visualizing the effects of magneto-immutability's organization on the turbulent flow. It designed to calculate the cosine of an angle $\theta$ that is representative of the alignment between the rate-of-strain tensor $\boldsymbol {\nabla }\boldsymbol {u}$ and the (spatially) local magnetic-field direction $\hat {\boldsymbol {b}}$, on a scale-by-scale basis in $k_{\perp }$. Initially used by St-Onge et al. (Reference St-Onge, Kunz, Squire and Schekochihin2020) without separation by scale for a study of the incompressible fluctuation dynamo under the action of Braginskii viscous stresses, it is based off of analytical theory from Kazantsev (Reference Kazantsev1968). In the fluctuation dynamo, changes to the magnetic-field strength are mediated by $\boldsymbol {\nabla }\boldsymbol {u}$ through its three eigenvalues, which are associated with field-line compressing motions, field-line stretching motions and the incompressibility constraint (ensuring that the overall flow stretches as much as it compresses). The compressing and stretching motions result in changes to $|B|$, thus the angles between the compressing and stretching eigenvectors and the local magnetic-field direction dictate how efficient the flow is at changing the magnetic-field strength.Footnote ⁹ In practice, we must ensure that the eigenvectors and eigenvalues are real before dotting into $\hat {\boldsymbol {b}}$, therefore we diagonalize $(\boldsymbol {\nabla }\boldsymbol {u} + \boldsymbol {\nabla }\boldsymbol {u}^\mathsf {T})/2$ with $\mathsf {T}$ denoting the transpose. Once the eigenvectors are obtained, the cosine of their angle with $\hat {\boldsymbol {b}}$ is calculated, yielding a compressing $\cos \theta$ and a stretching $\cos \theta$ at each grid point within the domain. These cosines are then compiled into a probability distribution $\mathcal {P}(|\cos \theta |)$, representing the likelihood that the alignment angle cosine takes on a specific value at a given point in time. We have found that the compressing and stretching distributions are qualitatively identical for all simulations performed in this work, which is expected for Alfvénic turbulence (as opposed to the fluctuation dynamo studied in St-Onge et al. Reference St-Onge, Kunz, Squire and Schekochihin2020). For that reason, it suffices to only show one of the two, which we choose to be the stretching angle. To exhibit the $\mathcal {P}(|\cos \theta |)$ as a function of scale, we first Fourier transform $\boldsymbol {u}$, filter it using masks in $k_{\perp } = \sqrt {k_x^2+k_y^2}$ bins and then transform it back to real space before performing all of the aforementioned operations to obtain $\mathcal {P}(|\cos \theta |)$. The $\mathcal {P}(|\cos \theta |)$ values are normalized to unity at each individual $k_{\perp }$, rather than $\mathcal {P}$ being a distribution in both $|\cos \theta |$ and $k$, so as to avoid weighting larger wavenumbers less than smaller ones. When magneto-immutability is active, our expectation that the flow organizes to avoid changes in $|B|$ implies that we should see a small $|\cos \theta |$ between the eigenvectors of $\boldsymbol {\nabla }\boldsymbol {u}$ and $\hat {\boldsymbol {b}}$, whereas a non-immutable cascade should have $|\cos \theta | \sim 1$. Therefore, we can compare this cosine between passive- and active-$\varDelta$ simulations to detect whether the flow behaves in a different fashion so as to produce a cascade with minimal variation of $|B|$.

3.3 The effective equation of state

The most fundamental assumption underlying all of the conclusions of § 2 is that the plasma behaves in a sufficiently collisionless manner to be described by the CGL-MHD equations. For this to be true, the portion of the plasma having pressure anisotropy that is microphysically unstable, and therefore subject to an effective collisionality associated with particle scattering off Larmor-scale fluctuations (i.e. that with $|\varDelta | \gtrsim \beta ^{-1}$), must constitute a small fraction of the total volume. In their initial study of magneto-immutability in the CGL-MHD system, Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) found that magneto-immutability suppresses the overall level of fluctuations in the pressure anisotropy, thereby reducing the fraction of the plasma that is unstable. To establish the impact that this suppression has on the plasma's effective collisionality, we plot in figure 2(a) the relationship between the fluctuations in the parallel and perpendicular pressures and the (relatively small) fluctuations in the density, averaged in time over the interval $tv_{A}/L_{\perp } = [8,10]$. In an MHD-like plasma, both $\delta p_{\perp }$ and $\delta p_{\|}$ would scale with $\delta \rho$ with a slope of $5/3$, the single-adiabatic index characteristic of a collisional plasma (black dotted lines). Neither is particularly well aligned with the single-adiabatic index slope, with $\delta p_{\|}$ having a steeper slope, and $\delta p_{\perp }$ a shallower one. The significant differences between the pressure slopes suggest that the $\delta p_{\perp /\|}$ evolve in an independent manner as is expected for a predominantly collisionless plasma.

Figure 2. (a) Probability distributions of fluctuations in $p_{\|}$, $p_{\perp }$ and $\rho$, for $\beta =10$ Alfvénically driven turbulence. The black dotted lines represent a slope of $5/3$, the expectation for single-adiabatic, collisional MHD. Neither $\delta p_{\|}$ nor $\delta p_{\perp }$ appear to align with the MHD prediction or with each other. (b) The fraction of the domain whose pressure anisotropy lies beyond the microinstability thresholds as a function of time for all active-$\varDelta$ (solid) and passive-$\varDelta$ (dashed) simulations. The passive-$\varDelta$ simulations typically have ${\sim }50\,\%$ of their volume unstable, whereas the unstable fraction of the active-$\varDelta$ simulations typically remains below $10\,\%$.

This non-MHD equation of state is explained by figure 2(b), which compares the fraction of the domain that is unstable to micro-instabilities between active- (solid curves) and passive-$\varDelta$ (dashed curves) simulations as a function of time. In the steady state of the passive-$\varDelta$ simulations, nearly half of the domain lies beyond the instability thresholds, meaning that much of the plasma volume is experiencing the large scattering rate $\nu _{\rm lim}$, and a collisionless model would not be a good description for a significant portion of the turbulent dynamics. Conversely, the active-$\varDelta$ simulations are rarely beyond 10 % unstable, with an average falling closer to just 5 % in steady state. Importantly, this reduction in the unstable fraction appears to be independent of $\beta$, suggesting that higher values of $\beta$ are not likely to result in a significantly more collisional plasma. Infrequent excursions of the unstable fraction beyond 10 % do occur in all simulations, which could either be the result of intermittency or the randomness of the forcing.Footnote ¹⁰ Although generally short lived, these brief events do have the ability to change the statistics of the turbulence. However, their cause is unlikely to originate in the inertial range, because the fluctuation amplitudes there are small and satisfy the ordering (2.4). As we demonstrate in § 3.6, magneto-immutability is least active at the outer scales where the turbulence is being forced, and so these intermittent bursts are likely driven by large-scale motions. That being said, the anisotropic-pressure stress is inherently non-local in $k$-space, so these bursts’ effects on the turbulence are not limited to the largest scales. Therefore, all statistical measurements we report are obtained by averaging over a time interval of no less than $2L_{\perp }/v_{A}$, taken beyond $tv_{A}/L_{\perp } = 6$ to average over these bursts and ensure that the turbulence statistics have reached an approximate steady state.

3.4 Compressive forcing and ion-acoustic fluctuations

With the knowledge that the turbulence is taking place within a predominantly collisionless plasma, we surmise that the compressive wave fluctuations will behave not as collisional fast- and slow-magnetosonic modes, but rather as collisionless fast and ion-acoustic waves. As such, we evaluate the feasibility of an ion-acoustic wave cascade at high $\beta$, given our expectation that ion-acoustic waves are not effectively mixed and cascaded by the Alfvénic fluctuations in such plasmas. Conveniently, the only waves that actively make $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}\ne 0$ in a collisionless plasma are the fast wave and the ion-acoustic wave (non-propagating modes are pressure balanced and approximately incompressible). Thus, the presence of ion-acoustic waves in the inertial range of these turbulence simulations can be diagnosed by the compressive flow spectrum of $\hat {\boldsymbol {k}} \boldsymbol {\cdot } \boldsymbol {u}_k$, which is provided in figure 3(a,b).

Figure 3. Spectra of compressive velocity fluctuations vs perpendicular wavenumber. (a) The compressive spectra at $\beta =100$ are studied as a function of the type of forcing. When the correlation time is Alfvénic, random and Alfvénic forcing both produce similar results, with little energy in the compressive flow. If the correlation time is instead sonic, then an increase of two orders of magnitude is seen in the compressive flow energy. (b) Compressive spectra of randomly driven turbulence at different values of the plasma $\beta$. At $\beta =100$ the forcing is sonic, leading to significant compressive fluctuations at the outer scale, but the spectrum is extremely steep. As $\beta$ is reduced, the spectra become decreasingly steep, qualitatively approaching the passive-$\varDelta$ result (black dashed).

In figure 3(a), we explore the dependence of the compressive flow spectrum on the correlation time of the forcing. When the correlation time of the forcing is Alfvénic (as is expected for the astrophysical turbulence we are concerned with here), very little energy is present in compressive fluctuations, and the power-law index of the $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}$ spectrum is rather steep. This is true in both the randomly driven and Alfvénically driven set-ups, with very little difference between the two spectra. Only when the correlation time of the forcing is decreased to be sonic ($t_{\rm corr} \sim L_{\|}/v_{\rm th}$), can a substantial amount of energy enter into compressive fluctuations. This is likely because, in Alfvénically correlated forcing, the time scale associated with the randomization of the forcing is too slow for the forced wavenumbers to significantly drive ion-acoustic waves.

Indeed, the $\beta$ dependence of this conclusion is captured by figure 3(b), where we plot these compressive spectra for only randomly driven simulations at $\beta =1,10$ and 100, with the last being sonically correlated. In effect, the $\beta =1$ simulation is also sonically correlated because $v_{\rm th} = v_{A}$, but because a large difference between the Alfvén and ion-acoustic wave speeds is not present, the mixing of ion-acoustic waves by Alfvén waves is stronger, bringing the spectrum much closer to that of the dashed passive-$\varDelta$ simulation, performed at $\beta =100$. The spectrum of the $\beta =100$ CGL-MHD simulation, although sonically correlated to generate substantial outer-scale compressive fluctuations, has an extremely steep spectrum. This steep spectrum implies that very little of the compressive mode energy driven at the outer scale penetrates into the inertial range, supporting our claim that ion-acoustic waves may be ordered out of the dynamics. At large wavenumbers the spectrum does become less steep. However, given the very small amount of energy contained in $E_{\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {u}}$ at such high $k$, this decrease in steepness is likely due to nonlinearities from Alfvénic fluctuations. The weaker cascade of compressive modes at high $\beta$ is also evident from the percentage of flow energy contained in compressive fluctuations. This percentage in the randomly forced $\beta =100$ simulations is $1\,\%$ for the Alfvénically correlated passive run, $0.3\,\%$ for the sonically correlated active run and $0.007\,\%$ for the Alfvénically correlated active run. Clearly, sonic correlation is a requirement for any substantial amount of compressive flow to be generated, but the spectra in figure 3(b) show this does not guarantee a strong turbulent cascade of that flow.

The consequences of the lack of compressive modes in the turbulent inertial range can be seen in the statistics of the density fluctuations. Figure 4(a) shows the probability distribution of various values of the density within the simulation domain for each of the Alfvénically correlated simulations, both active (solid) and passive (dashed). As $\beta$ is increased, the density takes on fewer values that deviate significantly from the background, with the Alfvénically driven simulations being only slightly narrower than the randomly driven runs. These narrow distributions are likely a result of the Mach number becoming smaller, given that our simulations drive an approximately fixed amplitude of $u/v_{A}$ at the outer scale. How the level of these fluctuations depends on scale is more clear in the spectrum of density fluctuations, figure 4(b). Interestingly, even though these fluctuations are clearly ${\lesssim } \epsilon ^2 \rho _0$, making them too small to affect the dynamics to leading order, and ion-acoustic wave mixing by the Alfvénic cascade is weak, their spectra appear to follow a near $-5/3$ power law, as given by the dotted line. We expect that this is accounted for by the presence of kinetic entropy and non-propagating modes, which, given that they have no real frequency, should be well mixed by Alfvénic fluctuations (see §§ 2.2 and 2.4). As they do not feed back on the overall dynamics due to their small amplitudes, their spectrum is probably captured by passive advection ($\mathrm {d}_t\delta \rho =0$), hence the $-5/3$ index.Footnote ¹¹ The $\beta =100$ simulations do appear to show a break from the single power-law behaviour of the other simulations, although this is likely an effect of our choice of $\nu _{\rm lim}$, which is discussed in more detail within § 3.8. Note that, in figure 4(a,b), all but one of the passive simulations possess density fluctuations larger than their active counterparts; only the $\beta =100$ Alfvénically driven passive simulation possesses $\delta \rho$ as small as its active counterpart. The fact that the choice between random and Alfvénic forcing makes a difference in the passive simulations at $\beta =100$, but not the active simulations, highlights the separation of compressive time scales between MHD and CGL. In CGL, the type of forcing makes little difference unless it is sonically correlated so that it may excite ion-acoustic modes. On the other hand, the correlation time need not be sonic to excite MHD slow modes, so random forcing can in fact drive larger density fluctuations in the passive simulation.

Figure 4. (a) Probability distribution of $\rho$ within the domain for all Alfvénically correlated simulations, active (solid) and passive (dashed). At higher $\beta$, the decrease in Mach number leads to weaker density fluctuations being driven at the outer scale, and less overall variation in $\rho$. Randomly driven simulations exhibit somewhat increased variation in $\rho$, but the dominant parameter is $\beta$. (b) Density fluctuation spectra for all Alfvénically correlated simulations. While the overall amplitudes are decreased with increasing $\beta$, the spectra remain strong with spectral indices near $-5/3$. This may indicate that the density is passively advected in the absence of fast and ion-acoustic wave cascades. The apparent break in power-law behaviour of the $\beta =100$ spectrum is a result of our particular choice for $\nu _{\rm lim}$, which is discussed in more detail within § 3.8. In both panels (a,b), most all of the passive simulations exhibit larger density fluctuations, with only the $\beta =100$ Alfvénically driven passive run having $\delta \rho$ as small as its active counterpart.

3.5 Comparison with previous work

Certain aspects of our ordering (2.4) and reduced equations (2.16) explain features of high-$\beta$ CGL turbulence that were already observed in the Alfvénically driven simulations of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023), but not yet fully understood. Here, we reproduce some of the key results of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023), confirming that our simulations explore the same effects and qualifying the extent to which the compressibility of forcing matters in a system where ion-acoustic and fast modes are not effectively cascaded. All simulations shown within this section are performed at $\beta = 10$ with Alfvénically correlated forcing and standard resolution (as defined in § 3.1).

In figure 5(a), the kinetic and magnetic energy spectra from the randomly (solid) and Alfvénically (dash-dotted) forced simulations reveal inertial ranges that are close to the $-5/3$ power law expected in MHD turbulence. Although the individual spectral slopes deviate very slightly above or below this exact value, overall, there is little qualitative difference between the turbulence resulting from the two modes of forcing. This indicates that the pressure-anisotropy stress does not effectively remove energy from the cascade, at least not to the extent that would naïvely be expected when $\beta \varDelta \sim 1$ (Squire et al. Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). Upon careful inspection, the power-law index of the kinetic energy spectrum in the randomly driven simulation is steeper than that of the Alfvénically driven simulation, although the difference is small, likely originating from next-order effects not captured by (2.16). In figure 5(b), the field-perpendicular and parallel scales of the Alfvénic flow and magnetic-field perturbations are given. For both $u_{\perp }$ and $B_{\perp }$, the perpendicular direction is defined with respect to the average magnetic field of all points being used in the difference equation of the related structure functions (e.g. Chen et al. Reference Chen, Mallet, Yousef, Schekochihin and Horbury2011). The dash-dotted and solid coloured lines once again represent the Alfvénic and randomly driven simulations, while the black dotted line represents the $(l_{\|}/l_0) \sim (l_{\perp }/l_0)^{2/3}$ relationship predicted for critically balanced MHD turbulence. This appears to indicate that cascades of $\delta B_{\perp }$ and $u_{\perp }$ are critically balanced, another foundational assumption of our ordering (2.4). At the largest scales, some disagreement occurs, but this is to be expected. The definition of $B_{\perp }$ becomes vague in the presence of larger-amplitude fluctuations, and $u_{\perp }$ defined with respect to the averaged field will differ substantially from $u_{\perp }$ defined with respect to the local field at each point in the structure function. Additionally, larger scales may still be somewhat influenced by the forcing, which does not exclusively drive critically balanced fluctuations. Further details of how these characteristic eddy sizes are measured can be found within § 3.2.2 of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) or § 6.3 of Cho & Lazarian (Reference Cho and Lazarian2009).

Figure 5. (a) Kinetic energy spectra of $\beta =10$, Alfvénically correlated simulations with random (solid) and Alfvénic (dash-dotted) forcing. Only slight qualitative differences are visible between each type of forcing, reflecting the lack of a strong ion-acoustic cascade regardless of forcing. (b) Characteristic turbulent eddy sizes along and across the local magnetic field for Alfvénic variables. Both the randomly and Alfvénically driven simulations closely follow the scaling relationship predicted by critical balance in standard MHD, $(l_{\|}/l_0) \sim (l_{\perp }/l_0)^{2/3}$, represented by the black dotted line.

Figure 6. (a) Spectra of the perpendicular-thermal, parallel-thermal and magnetic pressures for randomly (solid) and Alfvénically (dash-dotted) driven $\beta =10$ simulations. The difference between $p_{\|}$ and $p_{\perp }$, combined with the rough equivalence of $B^2$ and $p_{\perp }$, reflects perpendicular pressure balance (Squire et al. Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). (b) The rate-of-strain spectra, showing suppression of $\boldsymbol {\nabla }_{\|} u_{\|}$ as predicted by (2.16). The most noticeable difference between forcing modes occurs here in the spectra of $\boldsymbol {\nabla }_{\perp } u_{\|}$.

The most easily verified predictions of (2.16) are that the parallel pressure dominates the pressure anisotropy, and that this results from the perpendicular pressure being balanced by magnetic pressure. In figure 6(a), the thermal and magnetic pressure spectra are shown, depicting both of these features.Footnote ¹² The perpendicular pressure, although larger at the outer scale, approaches and ultimately matches the magnetic pressure in the inertial range, as predicted by (2.10a,b); as stipulated in the reduced model, this empirical result requires the leading-order $\delta p_{\perp }$ perturbation to be $\mathcal {O}(\epsilon ^2)$. By contrast, fluctuations in $p_{\|}$ are much larger in amplitude, and all appear to follow the spectral index of $-5/3$ typical of passive advection, with some minor variation. The inertial ranges of these pressure spectra bear a striking resemblance to those of the $\beta =16$, hybrid-kinetic simulation of Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023, figure 8a). This is particularly important given that the dynamics of $p_{\perp }$ and $p_{\|}$ is generally reliant on the heat fluxes, which the CGL-MHD and hybrid-kinetic approaches model in very different ways. For such agreement to exist between the two models, it is likely that the heat fluxes are reduced by the field lines becoming nearly isotherms. This supports a key prediction of our reduced system (2.16) (see (2.12) and (2.13)). This suppression of $q_{\perp /\|}$ is also verified specifically for our CGL-MHD simulations within § 3.7. The $p_{\|}$ spectrum appears to have a slightly shorter inertial range than the $p_{\perp }$ and $B^2$ spectra, a feature which may be due to the fact that the heat fluxes are stronger for $p_{\|}$ than for $p_{\perp }$. Although the $q_{\perp /\|}$ are nominally ordered out by the reduction in $\boldsymbol {\nabla }_{\|} \delta T_{\perp /\|}$, it is possible that as the turbulence approaches grid scales magneto-immutability weakens somewhat from the effects of finite resolution. We find that in higher resolution simulations this steepening in $E_{p_{\|}}$ trends with the grid scale, therefore it is likely a numerical artifact and has no impact on our physical understanding of the $p_{\|}$ cascade. Figure 6(b) depicts the rate-of-strain spectra of the turbulent flow, broken up into field-perpendicular and parallel gradients, as well as perpendicular and parallel flows. As expected from (2.16), parallel gradients of $u_{\|}$ are dramatically suppressed with respect to all other elements of the rate-of-strain tensor (Squire et al. Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). Importantly, this includes $\boldsymbol {\nabla }_{\perp } u_{\|}$, which emphasizes that it is not simply that $u_{\|}$ itself is being reduced, only its gradients along the magnetic field. Between figures 5 and 6, the most apparent difference between random and Alfvénic driving exists in the spectra of $\boldsymbol {\nabla }_{\perp } u_{\|}$. All other gradients of the flow velocity are similar between each mode of forcing. However, the randomly driven simulation features a spectrum of $\boldsymbol {\nabla }_{\perp } u_{\|}$ that is roughly constant with $k_{\perp }$, unlike in the Alfvénically driven run where it is increasing with $k_{\perp }$, in accordance with standard MHD scalings (Squire et al. Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). This highlights the fact that, even though the forcing has a clear effect on $u_{\|}$, it has no apparent effect on the suppression of $\boldsymbol {\nabla }_{\|}u_{\|}$, given that it is equally suppressed for both random and Alfvénic forcing.

One of the most significant consequences of magneto-immutability is the reduction of the anisotropic-pressure stress via suppression of $\boldsymbol {\nabla }_{\|} \Delta p$. To compare the degree to which this is achieved in random vs Alfvénically forced turbulence, we plot the $\boldsymbol {\nabla }_{\perp /\|}\Delta p$ spectra and the $\Delta p$ transfer function $\mathcal {T}_{\Delta p}(k_{\perp })$ of our $\beta =10$ simulations in figure 7. In figure 7(a), the active-$\varDelta$ simulations exhibit clear suppression of $\boldsymbol {\nabla }_{\|} \Delta p$, especially with respect to $\boldsymbol {\nabla }_{\|} \Delta p$ from the Alfvénically driven, $\beta =10$ passive run. That being said, the difference between the ratio $\boldsymbol {\nabla }_\parallel \Delta p/\boldsymbol {\nabla }_{\perp } \Delta p$ in the passive and active runs, which is only a factor of ${\approx }2$, is more subtle. Although the reduced equations only explicitly constrain the parallel gradient of $\Delta p$ through (2.12), the overall production of pressure anisotropy is also reduced by suppressing changes in the magnetic-field strength. As a result, perpendicular gradients of $\Delta p$ are also smaller simply because of the smaller overall magnitude of $\Delta p$ fluctuations.

Figure 7. (a) Spectra of the field-parallel and -perpendicular gradients of the pressure anisotropy for $\beta =10$ simulations that are either active and Alfvénically forced, active and randomly forced or passive and Alfvénically forced. Both of the active-$\varDelta$ simulations show a significant decrease in $\boldsymbol {\nabla }_{\|} \Delta p$ with respect to the passive simulation, as well as greater difference between $\boldsymbol {\nabla }_{\perp } \Delta p$ and $\boldsymbol {\nabla }_{\|} \Delta p$. (b) Transfer rate into and out of the turbulent flow due to the anisotropic-pressure stress, normalized to the total cascade rate, for the same simulations as (a). Both active simulations show significant suppression of the pressure stress resulting from the reduction of $\boldsymbol {\nabla }_{\|}\Delta p$ seen in (a). The passive simulation predicts a stress that is capable of damping turbulent motions entirely, since $\mathcal {T}_{\Delta p} \sim \mathcal {T}_{\rm total}$ across the full inertial range.

The results of $\boldsymbol {\nabla }_{\|} \Delta p$ suppression are seen in the energy transfer due to the pressure anisotropy $\mathcal {T}_{\Delta p}(k_{\perp })$ on figure 7(b) (3.2). The value of $\mathcal {T}_{\Delta p}(k_{\perp })$ is normalized to the total energy transfer rate, which we approximate by the Kolmogorov cascade rate $\mathcal {E}_{K}/\tau _{k_0} \sim \mathcal {E}_{K} (2\pi u_{\rm rms}/L_{\perp })$, where $\tau _{k_0}$ is the outer-scale turnover time. For the passive run, $\mathcal {T}_{\Delta p}/\mathcal {T}_{\rm total}$ being $\sim 1$ means that if the turbulence cascaded according to MHD, the anisotropic pressure stress would remove all of the energy from the turbulent fluctuations, ending the inertial range as soon as it begins. However, with immutability, this pressure stress is much smaller, permitting the inertial range to be relatively conservative. Overall, very little difference exists between the randomly and Alfvénically driven simulations.

3.6 Organizing the turbulence

A key aspect of why magneto-immutability constitutes self-organization is that (2.4) does not explicitly order $\boldsymbol {\nabla }_{\|} u_{\|}$ as small; rather it becomes necessary for $\boldsymbol {\nabla }_{\|} u_{\|}$ to be suppressed in order to satisfy both perpendicular pressure balance and $\mu$ conservation. Given that $\hat {\boldsymbol {b}}$ is always $\mathcal {O}(1)$, and figure 6(b) demonstrates that $u_{\|}$ is significant, it can only be the alignment angle between the magnetic field and the flow rate-of-strain tensor that suppresses variations in the magnetic-field strength. To probe this organization, we analyse our simulations using the alignment angle diagnostic described in § 3.2. Examples of the scale-dependent $\mathcal {P}(|\cos \theta |)$ for both passive and active Alfvénically driven simulations at $\beta =10$ are shown in figure 8. Recall that the dynamics in the passive simulations is simply isothermal MHD, and so the passively evolved pressure anisotropy does not affect $\mathcal {P}(|\cos \theta |)$. In figure 8(a), it is clear that the peak of the distribution from the active run closely follows $\cos \theta \approx 0$ until dissipation scales are reached at large $k_{\perp }$. This indicates a near-complete misalignment of the rate-of-strain eigenvectors and $\hat {\boldsymbol {b}}$. By contrast, the passive (MHD) run distribution in figure 8(b) is peaked near $\cos \theta \approx 0.6$, which yields an $\mathcal {O}(1)$ dot product between $\hat {\boldsymbol {b}}$ and $\boldsymbol {\nabla } \boldsymbol {u}$, thereby permitting more change to $|B|$. In figure 9, we show the dominant alignment angle as a function of $k_{\perp }$ for all active (solid) and passive (dashed) simulations by tracking the peak of $\mathcal {P}(|\cos \theta |)$ across $k_{\perp }$. Without exception, all active simulations exhibit distributions of $\cos \theta$ that are concentrated between 0 and ${\approx }0.2$, while all passive simulations fall between $\cos \theta \approx$ 0.5 and 0.7. In some cases, this misalignment begins very close to the outer scale of the simulation where $u_{\perp }/v_{A}\sim 1/2$. This might not be expected given that our ordering (2.4) assumes $u_{\perp } /v_{A} \ll 1$, however, figures 6(a) and 7 indicate that it is not required for the fluctuation amplitude to be very small (e.g. $\lesssim 10\,\%$) before qualitative features of the reduced (magneto-immutable) system are detectable.

Figure 8. Probability distributions of the cosine of the angle between the rate-of-strain stretching eigenvector and the local magnetic-field direction, for $\beta =10$, Alfvénically driven turbulence. Distributions are calculated and normalized individually within each $k_{\perp }$ bin. The distributions of the compressing eigenvector cosines are qualitatively indistinguishable from their stretching counterparts for both the active and passive runs. The active-$\varDelta$ simulation yields $\cos \theta \approx 0$ throughout the inertial range, indicating that motions in the flow that would normally lead to magnetic-field growth are misaligned with $\hat {\boldsymbol {b}}$, rendering them incapable of significantly perturbing $|B|$. By contrast, the passive-$\varDelta$ simulation has its peak probability around $\cos \theta \approx 0.6$, which produces an $\mathcal {O}(1)$ dot product between the two vectors and allows significant changes in the magnetic-field strength.

Figure 9. Peak cosines of the $|\cos \theta |$ probability distributions vs $k_{\perp }$ for various Alfvénically correlated simulations, both passive (dashed) and active-$\varDelta$ (solid). All active simulations exhibit large misalignment in the inertial range, while all passive simulations show greater alignment, permitting larger changes to $|B|$. There is no clear trend with $\beta$ or forcing mode in the scale at which the active-simulation misalignment ends and begins to resemble the passive simulations. For most this occurs near numerical dissipation scales, which is equivalent for all simulations shown. The lack of misalignment at larger scales is likely due to the forcing not respecting magneto-immutability.

Many of the diagnostics shown elsewhere within this work and in Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) demonstrate magneto-immutability by comparing with passive simulations that possess exactly identical parameters, such as forcing and resolution. Unfortunately, for particle-in-cell simulations (or reality), there exist no corresponding passive-$\varDelta$ simulations that can so accurately represent the MHD equivalent turbulence. With the $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}$ alignment diagnostic, however, there is no need to compare with an exactly analogous passive-$\varDelta$ simulation to determine whether magneto-immutability is at work, as small values of $|\cos \theta |$ alone are sufficient indication. This allows us to search for magneto-immutability within a kinetic framework that self-consistently determines the heat fluxes, which the Landau-fluid heat fluxes can only approximate. In figure 10, we show the calculation of $|\cos \theta |$ for a $\beta =16$, Alfvénically driven hybrid-kinetic turbulence simulation initially described within Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023). The code used to perform this simulation was the Pegasus++ hybrid kinetic-ion fluid-electron particle-in-cell code (Kunz, Stone & Bai Reference Kunz, Stone and Bai2014). The simulation shown resolves scales both above and below the ion-Larmor scale, additionally incorporating a non-zero electron temperature $T_{e} = T_{i}$. The peak of the $|\cos \theta |$ distribution is traced by the red line, while a comparison with the inertial-range average of all of our passive simulations is provided by the dashed black line. The misalignment is distinctly stronger for the hybrid-kinetic simulation than for the passive-$\varDelta$ simulations, lasting until just past the ion-Larmor scale. While at scales satisfying $k_{\perp } r_{L,{\rm i}} \lesssim 0.4$ the misalignment is weaker, the peak value of the cosine is still significantly less than the MHD average, and the vast majority of the probability density falls well below $|\cos (\theta )|=0.62$. Such imperfect misalignment is also seen in our simulations in figure 9, and is perhaps not surprising given that these scales are closest to the forcing scales. The peak in alignment at $k_{\perp } r_{L,{\rm i}} \sim 0.4$ specifically is likely explained by the presence of oblique firehose modes, which Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023) note grow most rapidly at $k_{\perp } r_{L,{\rm i}} \sim 0.4$ and are expected to enhance $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}$. Those authors also reported significant viscous dissipation near the outer scale, perhaps further evidence that immutability is unable to coexist with the forcing, a fact reflected in the imperfect misalignment demonstrated at the largest scales within figure 10. This outer-scale heating can also be seen to some degree in our simulations in figure 7, as well as in most of the simulations of Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) (who used a different forcing scheme to that used here). Figure 18 of Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023) demonstrates that energy transfer from the viscous stress appears to be suppressed in the inertial range relative to the outer scale, instead being dominated by the Maxwell stress. Additional evidence hinting at the presence of magneto-immutability in Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023) includes the tendency of the turbulence to avoid the instability thresholds (their figure 5), and significant suppression of the spectrum of $p_{\perp }$ fluctuations in order to maintain pressure balance with the electron and magnetic-field contributions (their figure 9).

Figure 10. Distribution of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}$ alignment angle $\theta$ as a function of $k_{\perp }$ for the $\beta =16$, Alfvénically driven and correlated simulation of Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023). This simulation was performed within a $(120.5 r_{L,{\rm i}})^2 \times 241 r_{L,{\rm i}}$ box, at a resolution of $384^2 \times 768$. The peak alignment angle cosine (red line) is seen to be very close to 0 in the inertial range, in contrast to the average cosine measured from our passive-$\varDelta$ (isothermal MHD) simulations (black dashed line). As with our CGL simulations, this misalignment is weaker near the outer scale, where immutability struggles to coexist with the forcing.

Figure 11. (a) Kinetic energy spectra for Alfvénically driven $\beta =10$ simulations at 3 different resolutions across $\boldsymbol {B}_0$, all with $n_{\|}=2n_{\perp }$. The dotted black line represents a $k^{-5/3}$ power law. (b) Alignment angle cosines as a function of $k_{\perp }$ for each of the three resolutions, with vertical dotted lines marking the transition away from misalignment of $\boldsymbol {\nabla } \boldsymbol {u}$ with $\hat {\boldsymbol {b}}$. It appears that the transition trends with the dissipation scale of the turbulence, likely as a result of departure from the ordering (2.4).

In the CGL simulations depicted by figures 8 and 9, there is no obvious trend of the scale (in $k_{\perp }$) at which immutability ceases to affect $\cos \theta$; most simulations trend toward MHD-like alignment at roughly the same scale. To understand the meaning of this scale, we compare the alignment angle cosines with the kinetic energy spectra of the $\beta =10$, Alfvénically driven simulation at several different resolutions in figure 11. The resolutions are given for the coordinates perpendicular to $\boldsymbol {B}_0$ of each run, with $n_{\|} = 2n_{\perp }$. A clear trend with resolution exists, meaning that the alignment angle cosine only changes when the flow reaches dissipation scales and the turbulence no longer follows the ordering (2.4). Notably, this misalignment persists to smaller scales than those satisfying $\delta B_{\|}/B_0 \gtrsim \beta ^{-1}$, where the change in the magnetic-field strength due to individual turbulent fluctuations is large enough to generate $\beta \varDelta \gtrsim 1$ (in the absence of immutability). This implies that, at least as far as the resolutions we can probe go, magneto-immutability is not an outer-scale effect but rather persists throughout the inertial range regardless of how small $\beta \varDelta$ is at a given $k_{\perp }$, albeit in the somewhat modified sense discussed in § 2.5.

3.7 The role of heat fluxes

Not only is the $\Delta p$ stress suppressed by the reduction of $\boldsymbol {\nabla }_{\|} \delta T_{\|}$ and $\boldsymbol {\nabla }_{\|} \Delta p$, the heat fluxes are as well. Importantly, as predicted by (2.16), the $\boldsymbol {\nabla }_{\|}\delta T_{\perp /\|}$ reductions are not diffusive (i.e. caused by the strong heat fluxes at high-$\beta$), but rather dynamical, originating from the momentum equation. For that reason, we should be able to artificially enhance or suppress the heat fluxes in our simulations – for example, by adjusting the parameter $|k_\parallel |$ – and observe little effect on the signatures of magneto-immutability. In figure 12 we gather plots of $\cos \theta$ and the $\boldsymbol {\nabla }_{\|} \Delta p$ spectra, expressed as a function of $k_{\perp }$, for simulations where $|k_{\|}|$ is increased or decreased by a factor of 100. In all runs depicted, $\beta =10$ and the forcing is Alfvénic, with a corresponding passive simulation given as an orange dotted line, given that it was performed with the nominal Landau wavenumber of $|k_{\|}| = 2\pi$. The simulation with $|k_{\|}|$ increased by 100 (purple curve in each plot) is effectively double-adiabatic MHD because of how weak the heat fluxes are.

Figure 12. (a) Alignment angle cosine vs $k_{\perp }$ for $\beta =10$, Alfvénically driven simulations with $|k_{\|}|$ in (2.2) increased (purple), unmodified (orange) and decreased (blue) by a factor of 100. Given a change in heat-flux strength of $10^4$, the difference between strong (blue) and weak (purple) heat fluxes appears to have little effect on $\theta$. An orange dashed line represents the passive-$\varDelta$ equivalent of the unmodified heat-flux run. (b) The spectra of $\boldsymbol {\nabla }_{\|} \Delta p$ for each run. The large suppression of parallel gradients in the pressure anisotropy reflect the suppression of $\boldsymbol {\nabla }_{\|} T_{\|/\perp }$, as predicted by (2.16). This limits the ability of heat fluxes to play a role in the turbulence, and importantly does not rely on their strength, with even the $k_{\|} = 200\pi$ simulation, which is effectively double adiabatic (purple), showing a comparable reduction. The reduced $\boldsymbol {\nabla }_{\|} \Delta p$ also implies that the heat fluxes do not interfere in the avoidance of significant $\Delta p$-stress.

In figure 12(a), little variation is seen in the alignment angle cosine as a function of $k_{\perp }$, with only the double-adiabatic run being misaligned slightly further in $k_{\perp }$ than the others. Considering that the heat fluxes are $10^4$ times stronger for the blue curve than the purple curve, the extension of misalignment by less than a factor of 1.5 in $k_{\perp }$ is an extremely small difference. In the bottom panel, the field-parallel gradients of the pressure anisotropy are shown, featuring strong suppression of $\boldsymbol {\nabla }_{\|} \Delta p$ in the active-$\varDelta$ simulations as compared with the passive run. Although the active simulations are all more similar to each other than the passive one, parallel gradients of the double-adiabatic run are moderately larger, by a factor of $\sim 2$. The comparable reduction of $\boldsymbol {\nabla }_{\|} \Delta p$ for all values of $|k_{\|}|$ shown implies that the heat fluxes have little effect on whether or not immutability is able to effectively avoid $\Delta p$-stress. Note that the results of figure 10 suggest that that this heat-flux suppression in figure 12 extends to the kinetic simulations of Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023), where the $q_{\perp /\|}$ are not approximated via the Landau-fluid form. In order for $\hat {\boldsymbol {b}} \boldsymbol {\cdot } \boldsymbol {\nabla } u_{\|}$ to be relegated to next order in (2.15), the heat fluxes must be negligible, a fact also made clear by (A5) of Appendix A. Therefore, to observe the kind of misalignment seen in figure 10, it must be that the fully kinetic heat fluxes, like the Landau-fluid approximations, have little effect on the leading-order dynamics. Unfortunately, unless a closure for $q_{\perp /\|}$ is assumed this cannot be proven in general. Nonetheless, these results still further the favourable comparison between Landau-fluid CGL and kinetic simulations of high-$\beta$ collisionless turbulence.

3.8 The role of micro-instabilities

The only portions of this study that are not described by our ordered equations (2.16) are the role of microinstability limiters. Within the CGL model, dependence on this physics arises via our choice of $\nu _{\rm lim}$. As discussed in § 3.3, the sizes and distributions of micro-unstable patches are highly intermittent and difficult to predict; however, by varying $\nu _{\rm lim}$ we can obtain useful information about how they interact with the turbulence. For all other simulations outside of this section, the instability scattering rate is fixed to $\nu _{\rm lim} = 10^{10} v_{A}/L_{\perp }$, essentially providing a hard wall on the maximum possible pressure anisotropy. Given that the microinstabilities being accounted for grow at $r_{L,{\rm i}}$ scales, which are formally zero in the CGL-MHD system, it seems reasonable that they might scatter particles at a rate much faster than any of the dynamics being studied. For this reason, and for simplicity, such ‘hard-wall’ limiters have been employed frequently in past studies of pressure-anisotropic turbulence far above $r_{L,{\rm i}}$ scales (see e.g. Sharma et al. Reference Sharma, Hammett, Quataert and Stone2006; Santos-Lima et al. Reference Santos-Lima, de Gouveia Dal Pino, Kowal, Falceta-Gonçalves, Lazarian and Nakwacki2014; Squire et al. Reference Squire, Schekochihin, Quataert and Kunz2019, Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023). However, studies focused on the mirror and firehose instabilities in the absence of background turbulence have found that the scattering rate induced typically follows the relationship $\nu \sim \beta \hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}$ (Kunz et al. Reference Kunz, Stone and Bai2014; Riquelme, Quataert & Verscharen Reference Riquelme, Quataert and Verscharen2018), which in a turbulent environment takes on different values at different scales. Unfortunately, this is expensive to implement directly as it requires measuring the flow shear locally at each time step. To investigate the consequences of different choices for $\nu _{\rm lim}$, we perform Alfvénically driven simulations at $\beta =100$, where the relative forcing of $\Delta p$ is larger than at $\beta =10$, and vary $\nu _{\rm lim}$ between 20, 200 and $10^{10} v_{A}/L_{\perp }$.

The effects of these variations in $\nu _{\rm lim}$ on the kinetic spectra are shown in figure 13(a). Interestingly simulations with lower scattering rates, where the pressure anisotropy is allowed to make larger excursions beyond the microinstability thresholds, lead to less steep spectral indices that are much more in line with the $-5/3$ expectation of our ordering and a conservative cascade. On the other hand, the simulation with a hard-wall limiter $\nu _{\rm lim}$, as used in all other simulations throughout this work, exhibits a spectrum that is slightly steeper and appears to be dissipated earlier in $k_{\perp }$. This steepening implies an increase in damping from the $\Delta p$ stress, as some of the energy in the turbulent motions is being removed from the cascade.Footnote ¹³ The reason for this steepening and apparent dissipation can be seen in figure 13(b,c), which present probability distribution functions (p.d.f.s) of the measured values of $\beta \varDelta$ and $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u}$ in each run. These diagnostics essentially represent how effective magneto-immutability is at regulating the overall magnitude and production rate of pressure anisotropy. In the hard-wall-limited passive simulation (purple dashed), the distribution of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u}$ is less peaked in the stable regions ($-2<\beta \varDelta <1$) than any of the active runs, with values of $\beta \varDelta$ that are effectively pinned to the mirror and firehose limiters, and no local peak in $\beta \varDelta$ between the limiters. The hard-wall-limited active simulation is the least peaked of the active simulations at $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u} = 0$, also having no local peak in $\beta \varDelta$ between the instability limiters. By comparison, lower values of $\nu _{\rm lim}$ appear to better focus the distribution of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u}$ around 0, and at least in the case of $\nu _{\rm lim} =200v_{A}/L_{\perp }$, maximize the proportion of $\beta \varDelta$ that lies between the instability limiters. Note that, in the process of pushing $\Delta p$ between the instability thresholds, a slight background ($k=0$) pressure anisotropy develops in the $\nu _{\rm lim} = 200v_{A}/L_{\perp }$ run, with the most probable value of $\beta \varDelta$ in figure 13(b) being $\sim -0.5$ rather than 0. This is simply an effect of the asymmetry of the instability thresholds, and is expected to occur whenever the outer-scale amplitude is sufficiently large that fluctuations can reach both microinstability thresholds. If the simulation were initialized with a homogeneous, background (undriven) pressure anisotropy, the eventual steady state would still resemble that shown here, as the microinstability limiters alone set the location of this peak in the p.d.f. The consequences of the varying strengths of magneto-immutable self-organization are shown in figure 13(d), via the $\Delta p$-stress transfer functions of each simulation. It is clear that the $\nu _{\rm lim}=200v_{A}/L_{\perp }$ simulation, where magneto-immutability appears to enforce most strictly the suppression of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}$, experiences the least amount of viscous dissipation from the pressure anisotropy. The $\nu _{\rm lim}=20v_{A}/L_{\perp }$ run dissipates slightly more, and the hard-walled simulation experiences significant stress far into the inertial range. This likely explains why the hard-walled simulation exhibits kinetic energy spectra that are both steeper and shorter than the other two runs (panel a).

Figure 13. (a) Kinetic energy spectra for $\beta =100$, Alfvénically driven simulations with $\nu _{\rm lim} \in [20, 200, 10^{10}]v_{A}/L_{\perp }$. Significant spectral steepening is observed for the simulation with hard-wall limiters ($\nu _{\rm lim} = 10^{10}v_{A}/L_{\perp }$, the default in all other simulations), as well as a slightly earlier apparent dissipation scale. (b) Probability distribution of the values of $\beta \varDelta$ for each simulation. The hard-wall-limited simulation exhibits significant peaks with cutoffs near the mirror and firehose thresholds, while lower $\nu _{\rm lim}$ simulations extend beyond the cutoffs. The $\nu _{\rm lim} = 200 v_{A}/L_{\perp }$ simulation is the only run having a distribution with a global maximum between the instability thresholds. (c) Probability distribution of the values of $\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u}$ for each simulation, demonstrating that the width of the distribution reflects the strength of magneto-immutable organization. The distribution from the hard-wall-limited simulation is the closest of the active simulations to the passive, non-immutable simulation, with $\nu _{\rm lim} = 200v_{A}/L_{\perp }$ being the narrowest, suggesting that there is an intermediate value of $\nu _{\rm lim}$ that allows magneto-immutability to act most effectively. (d) The $\Delta p$-stress transfer functions for each value of $\nu _{\rm lim}$. There is clearly less dissipation resulting from pressure anisotropy at $\nu _{\rm lim}=200v_{A}/L_{\perp }$ than in any other run, with the hard-walled simulation by far experiencing the most viscous dissipation.

We hypothesize that this weakening of magneto-immutability occurs because, whenever a strong limiter scattering rate $\nu _{\rm lim}$ is activated, the magneto-immutable orderings of (2.1e) and (2.1d) are broken. In an MHD plasma with $\nu ^{-1}$ being much less than any dynamical time scale, this does not lead to any dissipation because the scattering drives the anisotropy to zero, and hence there is no anisotropic-pressure stress. However, when the scattering is induced by microinstabilities, it only drives the anisotropy to $\beta \varDelta \sim 1$ levels. This is visible in the hard-wall-limited curve of figure 13(b), where the distribution of $\beta \varDelta$ is sharply peaked at the instability thresholds (although the proportion of points $\mathcal {P}(\varDelta ){\rm d}\varDelta$ at the thresholds is still not large). As a result, $\varDelta$ still remains dynamically important, but it is no longer capable of organizing to avoid $\Delta p$-stress, because the collisional term makes the dynamical equations different from those that support magneto-immutability. Therefore, the more severely the ordering is violated by the instability limiters, the more the anisotropic-pressure stress is allowed to affect the cascade of turbulent energy. Interestingly then, the results of figure 13 imply that throughout this work, many of the signatures of magneto-immutability we detect would become even stronger if a lower (and more realistic) value of $\nu _{\rm lim}$ were used. That being said, it appears that if $\nu _{\rm lim}$ is too small, then the organization can again become somewhat less efficient. This drop in efficiency may be a result of the amplitude of $\beta \varDelta$ exceeding the values expected by the ordering, possibly upsetting the $\beta \varDelta \sim 1$ assumption.

While there is clearly a choice for $\nu _{\rm lim}$ that maximizes magneto-immutability, the question remains of whether such a choice is physically motivated or not. An instructive comparison to make in this pursuit is between our optimal $\nu _{\rm lim}$ and the effective scattering rate measured in the high-$\beta$, hybrid-kinetic simulations of Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023). Those authors found that microinstability scattering in their simulations appeared to follow the relationship $\nu \sim \beta (\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u})_L$, where the $L$ subscript denotes that the quantity is estimated at the outer scale. We can estimate this rate for the simulations of figure 13 from the root-mean-square velocity $\langle u\rangle _{\rm rms} \approx 0.75 v_{A}$ to find that $\nu \sim 235v_{A}/L_{\perp }$. Interestingly, this nearly coincides with the value of $\nu _{\rm lim}$ in our survey that is minimally disruptive to magneto-immutability.Footnote ¹⁴

(3.3)\begin{equation} \beta \frac{k_{\|}v_{A}}{\nu} \gg 1. \end{equation}

If we naively substitute the result of Kunz et al. (Reference Kunz, Stone and Bai2014) that, for pressure anisotropy driven by a flow shear, $\nu \sim \beta \hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u} \sim \beta k_{\|}u_{\|}$, then the condition for magneto-immutability subject to microinstability scattering is just that $u_{\|} \ll v_{A}$. This condition is automatically satisfied by virtue of our assumed ordering. The only effect that can then interrupt magneto-immutability in sub-Alfvénic turbulence is the forcing, hence why Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023) obtained a scattering rate based on the parallel rate-of-strain at the outer scale. Indeed, our figure 9 also demonstrates how magneto-immutable self-organization tends only to fail at the forcing scales for different values of and types of forcing. Thus, the value of $\nu _{\rm lim}$ in our simulations should also be estimated via $\nu \sim \beta (\hat {\boldsymbol {b}}\hat {\boldsymbol {b}} \boldsymbol {:} \boldsymbol {\nabla } \boldsymbol {u})_L$. In short, it is therefore likely that the most physically motivated choice of $\nu _{\rm lim}$ is also approximately that which best supports a strong, magneto-immutable cascade.

In considering our choice of $\nu _{\rm lim}$, it is worthwhile to also discuss alternative sources of anomalous scattering, such as that resulting from the ion-cyclotron instability or distinct versions of the firehose.Footnote ¹⁵ In Bott et al. (Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021) and Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023), the limiting negative pressure anisotropy is $-1.4/\beta$ corresponding to the oblique firehose instability, rather than the $-2/\beta$ employed in our simulations which represents the parallel (or ‘fluid’) firehose. If the limiters in our simulations activated at the oblique threshold, some minor quantitative results of the turbulence would change, but we do not expect any significant qualitative changes. In particular, the peak value of the p.d.f. of $\beta \varDelta$, which takes on $\approx -0.5$ in figure 13(b), would likely shift to $\approx -0.2$ due to there being a new centre between the mirror/firehose thresholds. Additionally, we expect a slightly larger fraction of microinstabilities given the more stringent threshold. To estimate this fraction, we calculated the average microinstability filling fraction with the oblique firehose threshold (rather than the parallel) for the final $\delta t=2L_{\perp }/v_{A}$ of the $\beta =100$, $\nu _{\rm lim}=200v_{A}/L_{\perp }$ simulation featured in figure 13. We find that $18.5\,\%$ of the domain is unstable with the oblique threshold, compared with $10.4\,\%$ with the parallel threshold. While significant quantitatively, this increase is not nearly sufficient to isotropize the entire plasma. Moreover, this estimated increase is almost certainly an overestimate, given that in this simulation no anomalous scattering ($\nu _{\rm lim}$) was turned on until the anisotropy reached $-2/\beta$. By comparison, the volume-filling fraction in the beta=16 hybrid-kinetic simulation of Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023, figure 10), which appears to obey the $-1.4/\beta$ threshold, was $17.9\,\%$. Finally, in their broad parameter survey, Squire et al. (Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2023) performed a simulation adopting the oblique firehose threshold, finding in their figure 4 that there was little qualitative difference from the simulations employing the parallel threshold.

With regards to the positive $\Delta p$ instabilities, Ley et al. (Reference Ley, Zweibel, Miller and Riquelme2024) found that a secondary ion-cyclotron instability can be triggered if the mirror threshold is exceeded substantially, causing the positive pressure anisotropy to be limited instead at ${\approx }0.5/\sqrt {\beta }$. However, that threshold is not likely to be relevant to our study for the following reasons. In order for this softer threshold to come into play, the pressure anisotropy must be able to overshoot the more restrictive $1/\beta$ mirror threshold by a substantial amount. Dedicated kinetic studies of the mirror instability in a driven system indicate that the overshoot is proportional to $(\hat {\boldsymbol {b}}\hat {\boldsymbol {b}}\boldsymbol {:}\boldsymbol {\nabla }\boldsymbol {u}/\varOmega _{i})^{1/2}$ (Kunz et al. Reference Kunz, Stone and Bai2014, where $\varOmega _{i} = eB/m_{{\rm i}}c$ is the ion-cyclotron frequency). In well-magnetized plasmas like the ICM or black-hole accretion flows, where the separation between the scales driving the pressure anisotropy (which are comparable to the outer scale when magneto-immutability is active) and the ion-Larmor scale is extremely large, this overshoot is predicted to be extremely small. In this case, the mirror instability should reach its saturated state and scatter particles to keep the pressure anisotropy bound within ${\approx }1/\beta$ before the pressure anisotropy can reach the ion-cyclotron threshold at ${\approx }0.5/\sqrt {\beta }$.