2.1. Overview

The existence of qualitatively distinct states in firehose-susceptible, high- $\beta$ plasmas stems in part from properties of the instability in its linear stage. In this section, we therefore describe the linear theory of the firehose instability. Though the linear theory of firehose instabilities has been discussed extensively in prior research (which we review in § 2.2), previously reported results do not completely account for the instability’s properties in high- $\beta$ plasmas. We therefore report a new analytical and numerical linear study in this regime (§ 2.3). We find that oblique firehose modes are dominant for $\beta _i \gg 1$ , with parallel ion-Larmor-scale firehose modes always having a smaller growth rate, in contrast to plasmas with $\beta _i \sim 1$ . Furthermore, the value of the pressure anisotropy at which ion-Larmor-scale oblique firehose modes are destabilised ( $\Delta _i = p_{\perp i}/p_{\| i} - 1 \simeq -1.35/\beta _{\|i}$ ) is less negative than that for longer-wavelength firehose modes at fixed $\beta _i$ , and is similar to the threshold value in $\beta _i \sim 1$ plasma. Aided by analytic theory, we explain these results in §§ 2.4 and 2.5, respectively.

2.2. A review of the firehose instability’s linear theory

Although a comprehensive understanding of the linear theory of the firehose instability (including at kinetic scales) was only obtained in the last few decades, the instability itself was first identified well over sixty years ago. The first studies of the firehose instability (Chandrasekhar, Kaufman & Watson Reference Chandrasekhar, Kaufman and Watson1958; Parker Reference Parker1958; Vedenov & Sagdeev 1958) showed that the dispersion relation of long-wavelength Alfvén waves (i.e. those modes with frequency $\omega$ whose parallel and perpendicular wavenumbers satisfy $k_{\|} \rho _i \ll |\Delta _i + 2/\beta _{\|i}|^{-1/2} \lesssim 1$ and $k_{\perp } \rho _i \ll 1$ , respectively) is

(2.1) \begin{equation} \omega ^2 = k_{\|}^2 v_{\textrm {A}}^2 \left (1+\frac {\Delta _i \beta _{\|i}}{2}\right ) = k_{\|}^2 v_{\textrm {A,eff}}^2 . \end{equation}

These modes become linearly unstable if the ion-pressure anisotropy $\Delta _i$ satisfies $\Delta _i \lt -2/\beta _{\|i}$ (or, equivalently, if $v^2_{\textrm {A,eff}} \lt 0$ ). We identify this condition as the ‘fluid’ firehose instability threshold, and the resulting instability as the non-resonant (or fluid) firehose instability. The instability of these modes can be understood physically as follows: once the parallel ion pressure exceeds the perpendicular pressure by an amount equal to twice the magnetic energy, parallel pressure forces on an Alfvénic perturbation can overpower the restoring magnetic tension that, in a pressure-isotropic plasma, is responsible for the wave’s propagation.

It is immediately clear from (2.1) that shorter-wavelength perturbations grow more rapidly than longer-wavelength ones, implying that the scale of the fastest-growing firehose modes must be determined by finite-Larmor-radius (FLR) effects. For non-resonant parallel firehose modes, these FLR effects can be characterised analytically (Shapiro & Shevchenko Reference Shapiro and Shevchenko1963; Kennel & Sagdeev Reference Kennel and Sagdeev1967; Davidson & Völk Reference Davidson and Völk1968), with the parallel wavenumber $k_{\|,\mathrm{peak}} \sim |\Delta _i + 2/\beta _{\|i}|^{-1/2} \rho _i^{-1}$ at which peak growth occurs being determined by gyroviscosity, i.e. the off-diagonal components of the pressure tensor associated with agyrotropy in the distribution function (e.g. Schekochihin et al. Reference Schekochihin, Cowley, Rincon and Rosin2010). Whenever $|\Delta _i + 2/\beta _{\|i}| \ll 1$ , which is either achieved near threshold (that is, when $|\Delta _i + 2/\beta _{\|i}| \ll 1$ in plasmas with $\beta _i \sim 1$ , or whenever $|\Delta _i| \ll 1$ in high- $\beta _i$ plasmas), the wavelength of the fastest-growing non-resonant mode is much larger than the ion-Larmor scale.

More recent studies that solved the hot-plasma dispersion relation numerically for a bi-Maxwellian plasma discovered the existence of two kinetic variants of the firehose instability: the resonant parallel firehose instability (Gary et al. Reference Gary, Li, O’Rourke and Winske1998) and the oblique firehose instability (Yoon, Wu & de Assis Reference Yoon, Wu and de Assis1993; Hellinger & Matsumoto Reference Hellinger and Matsumoto2000). Modes of the resonant parallel firehose instability are destabilised by gyroresonant interactions with suprathermal ions having parallel velocities $v_{\|} = (\varpi + \varOmega _i)/k_{\|}$ (where $\varpi$ is the real frequency of the mode). These modes have a characteristic parallel wavenumber $k_{\|} \rho _i \sim 1$ when $\beta _{\|i} \sim \Delta _i \sim 1$ , are circularly polarised, right-handed and propagating. Although the instability is technically thresholdless (Sagdeev & Shafranov Reference Sagdeev and Shafranov1960; see also Appendix A.2), previous numerical studies found that such modes only attain growth rates $\gamma _{\| \textrm {f}}$ at ion-Larmor scales that are not infinitesimal fractions of the ion-Larmor frequency when $\Delta _i$ exceeds some $\beta _{\|i}$ -dependent threshold. For example, Matteini et al. (Reference Matteini, Landi, Hellinger and Velli2006) report that, in order for $\gamma _{\|\textrm {f}} \gtrsim 5 \times 10^{-3} \varOmega _i$ , one requires that $\Delta _i \lesssim -0.6 (\beta _{\|i}-0.63)^{-0.58}$ . By contrast, oblique firehose modes are non-propagating and linearly polarised, with $k_{\|} \rho _i \sim k_{\perp } \rho _i \sim 0.5$ . Studies with $\beta _{\|i} \gtrsim 1$ identified a threshold that scales with $\beta _{\|i}$ in the same way as the fluid firehose threshold, but with a less negative numerical prefactor: $\Delta _i \lesssim -1.4 \beta _{\|i}^{-1}$ (Hellinger & Matsumoto Reference Hellinger and Matsumoto2000, Reference Hellinger and Matsumoto2001). These conditions together imply that, when $\beta _{\| i} \sim 1$ , the resonant parallel firehose instability tends to dominate, but that the oblique firehose instability should become dominant when $\beta _{\| i} \gg 1$ .

The less negative values of the pressure anisotropy required for the resonant and oblique firehose instabilities to operate linearly at ion-Larmor scales have been considered and discussed extensively for collisionless, $\beta _{\|i} \gtrsim 1$ plasma similar to the solar wind (Hellinger et al. Reference Hellinger, Trávníček, Kasper and Lazarus2006; Matteini et al. Reference Matteini, Landi, Hellinger, Pantellini, Maksimovic, Velli, Goldstein and Marsch2007, Reference Matteini, Hellinger, Landi, Trávníček and Velli2012, Reference Matteini, Hellinger, Goldstein, Landi, Velli and Neugebauer2013). Several studies present results of direct relevance to high- $\beta$ plasmas. For example, in addition to identifying the existence of the resonant parallel firehose instability, Gary et al. (Reference Gary, Li, O’Rourke and Winske1998) characterise its linear threshold for $\beta _i \leqslant 10$ . Hellinger et al. (Reference Hellinger, Trávníček, Kasper and Lazarus2006) compute linear instability thresholds for both the resonant parallel and oblique firehose instabilities in a bi-Maxwellian plasma for $\beta _i \leqslant 30$ . However, complementary results for linear firehose instability thresholds in plasmas with larger $\beta _i$ have not been the focus of any previous published studies, nor have results for peak growth rates. We therefore report these results in the next section.

2.3. The kinetic firehose instability at $\beta _{i} \gg 1$

To determine the linear thresholds and growth rates of firehose-unstable modes as a function of wavenumber in a plasma with bi-Maxwellian ions and Maxwellian electrons,Footnote ² we solve their linear dispersion relation numerically. We take the electric field $\delta {\boldsymbol{E}}$ and magnetic field $\delta {\boldsymbol{B}}$ associated with such perturbations to be of the form

(2.2) \begin{equation} \delta {\boldsymbol{E}} \propto \exp [\textrm {i}({\boldsymbol{k}} \,{\boldsymbol{\cdot }}\, {\boldsymbol{r}} - \omega t)] , \quad \delta {\boldsymbol{B}} \propto \exp [\textrm {i}({\boldsymbol{k}} \,{\boldsymbol{\cdot }}\, {\boldsymbol{r}} - \omega t)] , \end{equation}

where $\boldsymbol{k}$ and $\omega$ are the wavevector and (complex) frequency of the perturbation. The dispersion relation of firehose perturbations having arbitrary $\boldsymbol{k}$ is the hot-plasma dispersion relation (Davidson Reference Davidson1983), which we provide for a plasma with arbitrary distribution functions in Appendix A.1. For a hydrogenic plasma with bi-Maxwellian ions (with parallel temperature $T_{\|i}$ and perpendicular temperature $T_{\perp i}$ ) and Maxwellian electrons (with temperature $T_e$ ), the dispersion relation simplifies to

(2.3) \begin{equation} \mbox{det}\Bigg \{\frac {k^2 \rho _i^2}{\beta _{\perp i}} \Bigg [\hat {{\boldsymbol{k}}}\hat {{\boldsymbol{k}}} - \unicode{x1D644}\left (1-\frac {\omega ^2}{k_\|^2 v_{\mathrm{th}\perp i}^2} \frac {v_{\mathrm{th}\perp i}^2}{c^2}\right )\Bigg ] + \widetilde {{\boldsymbol{\sigma }}}_{\textrm {bi-M}} \Bigg \}=0 , \end{equation}

where $v_{\mathrm{th}\perp i} \equiv \sqrt {2 T_{\perp i}/m_i}$ and $\widetilde {{\boldsymbol{\sigma }}}_{\textrm {bi-M}} = \widetilde {{\boldsymbol{\sigma }}}_{\textrm {bi-M}} (k_{\|} \rho _i,k_{\perp } \rho _i,{\omega /k_\| v_{\mathrm{th}\perp i}},{m_e/m_i},{T_e/T_{\|i}},\Delta _i )$ is a dimensionless rank-three tensor that can be written in terms of the plasma dispersion function and sums of modified Bessel functions (see Appendix A.1). To find the complex frequency $\omega$ of firehose-unstable modes at fixed values of $\beta _i$ , $m_e/m_i$ , $T_e/T_{\|i}$ , $v_{\mathrm{th}\perp i}/c$ and $\Delta _i$ , we choose values of $k_{\|} \rho _i$ and $k_{\perp } \rho _i$ at which such modes are expected to be realisable, and then solve for the roots ${\omega }/{k_{\|} v_{\mathrm{th}\perp i}}$ of (2.3). Numerically, this is carried out using the secant method, with the initial guesses inputted into the algorithm being determined by an analytical approximation to the hot-plasma dispersion relation that is valid when $\beta _i, \beta _e \gg 1$ (taken from Bott et al. Reference Bott, Cowley and Schekochihin2024).

Figure 2 shows the growth rate $\gamma = \mathrm{Im}(\omega )$ of firehose-unstable modes in a plasma with $\beta _{i} = 200$ as a function of $k_{\|} \rho _i$ and $k_{\perp } \rho _i$ for representative choices of the other parameters. Similarly to prior numerical studies in $\beta _i \gtrsim 1$ plasma, we observe that, as $\Delta _i$ is decreased from zero towards the (negative) instability thresholds, modes whose growth rates are not infinitesimally small first emerge only at ion-Larmor scales, at a critical value of the ion pressure anistropy, $\Delta _{\textrm {cr}} \simeq -1.35/\beta _{i}$ , that is less negative than the fluid firehose threshold $\Delta _{i} = -2/\beta _{i}$ at fixed $\beta _i$ (figure 2 $a$ ). Of these modes, the fastest-growing ones are oblique firehose modes ( $k_{\|} \rho _i \simeq 0.45$ , $k_{\perp } \rho _i \simeq 0.35$ ) with zero real frequency. As $\Delta _i$ is decreased further, the region of $(k_{\|} \rho _i,k_{\perp } \rho _i)$ space over which the firehose instability operates extends, with long-wavelength modes becoming unstable once $\Delta _i \lt -2/\beta _i$ (see figure 2 $b$ ).

Figure 2. ( $a$ ) Linear growth rate $\gamma$ of firehose-unstable modes as a function of parallel and perpendicular wavenumber for a range of different $\Delta _i$ at $\beta _{i} = 200$ , $m_i/m_e = 1836$ , $T_e = T_{\|i}$ and $v_{\mathrm{th}e}/c = 0.05$ . The growth rates are calculated on a $400^2$ grid in ( $k_{\|} \rho _i, k_{\perp } \rho _i$ ), with equal logarithmic spacing in both directions. ( $b$ ) Critical value of $\Delta _i$ below which firehose instability onsets, $\Delta _{\textrm {cr}}$ , as a function of parallel and perpendicular wavenumber at $\beta _{i} = 200$ . ( $c$ ) Peak growth rate $\gamma _{\textrm {peak}}$ of the firehose instability as a function of $(\Delta _i-\Delta _{\textrm {cr}}) \beta _i$ for a range of $\beta _i$ (solid lines). The dashed lines shows the semi-analytic result (2.4), the red dotted line shows the power-law scaling (2.8) that empirically is a good fit for moderately large $\beta _i$ . ( $d$ ) Plot of $\gamma _{\textrm {peak}}$ as a function of $\beta _i$ for a range of $|\Delta _i|\beta _i$ . The blue dashed line shows (2.4); the red dot-dashed line shows the analytic result (2.6); and the red dotted line shows the power law (2.7).

The growth rate $\gamma _{\textrm {peak}}$ of the fastest-growing firehose-unstable modes is an increasing function of $-(\Delta _i-\Delta _{\textrm {cr}}) \beta _{i}$ (see figure 2 $c$ ) and a decreasing function of $\beta _i$ at fixed $\Delta _i \beta _i$ (see figure 2 $d$ ). The particular scaling of $\gamma _{\textrm {peak}}$ with $\beta _i$ at fixed $\Delta _i \beta _i$ depends on the latter’s exact value. When $|\Delta _i - \Delta _{\textrm {cr}}| \beta _i \ll 1$ , we find that

(2.4) \begin{equation} \gamma _{\textrm {peak}} \approx 2.9 |\Delta _i-\Delta _{\textrm {cr}}| \varOmega _i \quad (\mathrm{for} \;|\Delta _i - \Delta _{\textrm {cr}}| \beta _i \ll 1). \end{equation}

When $\Delta _i = -2/\beta _i$ , a different scaling can be derived analytically (Bott et al. Reference Bott, Cowley and Schekochihin2024):

(2.5) \begin{equation} \gamma _{\textrm {peak}} \approx \frac {5^{1/2}}{2^{3/2}} (k_{\|} \rho _i)_{\mathrm{peak}}^2 \frac {\varOmega _i}{\beta _i^{1/2}} \quad (\mathrm{for} \;\Delta _i = -2/\beta _i) , \end{equation}

where the characteristic parallel wavenumber of the fastest-growing mode, $ (k_{\|} \rho _i)_{\mathrm{peak}}$ , is a weakly varying function of $\beta _i$ that is given by special mathematical functions. For values of $\beta _i$ that are very large ( $\beta _i \gg 10^{6}$ ), the following simple expression for $(k_{\|} \rho _i)_{\mathrm{peak}}$ (and therefore $\gamma _{\textrm {peak}}$ ) can be found through an asymptotic analysis (Bott et al. Reference Bott, Cowley and Schekochihin2024):

(2.6) \begin{align} (k_{\|} \rho _i)_{\mathrm{peak}} & \approx \left [\frac {2}{\log {(27 \pi \beta _i/10)}}\right ]^{1/2} , \nonumber \\[4pt] \gamma _{\textrm {peak}} & \approx \frac {5^{1/2}}{2^{1/2}} \frac {\varOmega _i}{\beta _i^{1/2} \log {(27 \pi \beta _i/10)}} & \quad (\mathrm{for} \;\Delta _i = -2/\beta _i, \; \beta _i \mathrm{\; very \; large}). \end{align}

By contrast, for values of $\beta _i$ that are only moderately large ( $\beta _i \in [10,10^4]$ ), it can be shown empirically via fitting to the direct numerical solution of the linearised Vlasov equation that $\gamma _{\textrm {peak}}$ at $\Delta _i = -2/\beta _{\|i}$ has, to a very good approximation, a simple power-law dependence on $\beta _{i}$ :

(2.7) \begin{equation} \gamma _{\textrm {peak}} \approx 0.3 \varOmega _i \beta _i^{-0.6} \quad (\mathrm{for} \;\Delta _i = -2/\beta _i, \; 1 \ll \beta _i \ll 10^{5}). \end{equation}

This result can be extended to $\Delta _i$ close to (but not exactly equal to) $-2/\beta _{\|i}$ , for which we find that the peak growth rate is related to the pressure anisotropy via a simple power-law scaling:

(2.8) \begin{equation} \gamma _{\textrm {peak}} \approx 0.4 |\Delta _i-\Delta _{\textrm {cr}}|^{0.6} \varOmega _i \quad (\mathrm{for} \;\Delta _i \approx -2/\beta _i, \; 1 \ll \beta _i \ll 10^{5}). \end{equation}

The validity of these asymptotic approximations is tested in figure 2( $c$ , $d$ ). Figure 2( $c$ ) confirms that, in the relevant regime, the expressions (2.4) and (2.8) are good approximations to the numerically determined growth rate as a function of $\Delta _i - \Delta _{\textrm {cr}}$ . Furthermore, figure 2( $d$ ) shows that the decrease of $\gamma _{\textrm {peak}}$ with increasing $\beta _i$ is primarily accounted for by the $\beta _i^{-0.6}$ dependence included in (2.7) for $\beta _i \ll 10^5$ . For quantitative agreement over a larger range of $\beta _i$ , an exact power-law fit is an oversimplification, as shown by the better agreement of the numerically determined growth rate with (2.6).

In summary, we find that, in plasma with $\beta _{\|i} \gg 1$ , the fastest-growing unstable modes are oblique firehose modes, and that these modes emerge at less negative pressure anisotropies ( $\Delta _i\lesssim -1.35/\beta _{\|i}$ ) than fluid firehose modes ( $\Delta _i \lt -2/\beta _{\|i}$ ). Furthermore, the resonant parallel firehose instability does not feature significantly in our numerical solution of the dispersion relation, seeming to imply that it is subdominant to the oblique instability in $\beta _{\|i} \gg 1$ plasma. We account for both of these findings in §§ 2.4 and 2.5, respectively.

2.4. Why the threshold of the oblique firehose instability is larger than $-2/\beta _{\parallel i}$

The numerical result that kinetic-scale oblique firehoses in a bi-Maxwellian, $\beta _i \gg 1$ plasma are destabilised at a less negative value of the pressure anisotropy can be elucidated by physical arguments and additional mathematical analysis.

The physical basis for a reduced threshold arises from modifications to the effective parallel pressure force acting on a magnetic-field perturbation when the thermal ion-Larmor radius is only a finite fraction of that perturbation’s wavelength. As explained in § 2.2, the fluid firehose instability is an instability of Alfvén waves in which, due to an excess of parallel pressure compared with perpendicular pressure, parallel pressure forces on the perturbed volume of plasma associated with the Alfvén wave become sufficiently large to overcome the restorative perpendicular pressure and magnetic-tension forces. When the scale of the perturbation is not much larger than $\rho _i$ , thermal ions are less well ‘tied’ to the field line’s trajectory because of their gyromotion. This results in an additional contribution to the net flux of perpendicular momentum into the perturbed volume of plasma, and therefore enhanced parallel pressure forces.

That the reduced threshold is a FLR effect can be proven analytically by taking advantage of the numerical observation that marginally unstable oblique firehose modes have no real frequency. Using this fact, we can derive a somewhat simplified (but still transcendental) equation for the threshold condition of the instability for a plasma with arbitrary ion and electron distribution functions. Via a subsidiary expansion in $k_{\|} \rho _i \sim k_{\perp } \rho _i \ll 1$ , we can then write down a simple expression for the threshold condition that includes the leading-order FLR corrections (see Appendix A.3). These corrections are proportional to high-order moments of the distribution function (specifically, fourth order or higher). For a plasma with bi-Maxwellian ions and Maxwellian electrons, we deduce that the threshold condition is

(2.9) \begin{equation} \Delta _i \left (1+\frac {3}{2} k_{\|}^2 {\rho }_i^2 - \frac {3}{8} k_{\perp }^2 {\rho }_i^2\right )+ \frac {2}{\beta _{\|i}} = \mathcal{O}\big(\Delta _i k^4 {\rho }_i^4\big) . \end{equation}

For any perturbation having $k_{\|} \lt 2 k_{\perp }$ , the condition (2.9) implies that the value of $\Delta _i$ required for instability is less negative than the fluid firehose threshold $\Delta _i = -2/\beta _{\|i}$ . For modes with the same wavevector as those oblique modes that we observe numerically to become unstable at $\Delta _i \approx \Delta _c \simeq -1.35/\beta _i$ (i.e. $k_{\|} \rho _i \simeq 0.45$ , $k_{\perp } \rho _i \simeq 0.35$ ), (2.9) implies $\Delta _i \approx -1.6/\beta _{\|i}$ , which is (to the order of accuracy of the subsidiary expansion) not too dissimilar to the numerically determined result. This agreement supports the conjecture that FLR effects are responsible for the observed weakening in the instability threshold.

2.5. Why the resonant parallel firehose instability is subdominant in $\beta _{\|i} \gg 1$ plasma

The apparent unimportance of the resonant parallel firehose instability when $\beta _{\|i} \gg 1$ , a finding consistent with previous numerical results (see § 2.2), can be proven analytically. We show in Appendix A.2 that, when $\Delta _i \simeq -1.35/\beta _{\|i}$ , the fastest-growing resonant parallel firehose modes (which, in contrast to plasma with $\Delta _i \sim \beta _{\|i} \sim 1$ , satisfy $k_{\|} \rho _i \ll 1$ ) have a growth rate that is exponentially small in $1/\beta _{\|i}$ , i.e. $\gamma _{\|\mathrm{f}} \sim \beta _{\|i}^{-1/2} \exp {(-0.74\beta _{\|i})}$ . By comparison, the peak growth rate $\gamma _\perp \textrm {f}$ of the resonant oblique firehose instability satisfies ${\gamma _{\perp \textrm {f}}} \sim \left |\Delta _i -\Delta _{\textrm {cr}}\right | {\varOmega _i}$ when $\Delta _i$ is close to the instability’s threshold anisotropy, $\Delta _{\textrm {cr}} \simeq -1.35/\beta _{\|i}$ (cf. (2.4)). Assuming that $|\Delta _i - \Delta _{\textrm {cr}}| \ll 1/\beta _{\|i}$ , it can be shown that $\gamma _{\perp \textrm {f}}$ greatly exceeds $\gamma _{\|\textrm {f}}$ when

(2.10) \begin{equation} (\Delta _{\textrm {cr}}-\Delta _i) \beta _{\|i} \gg \beta _{\|i}^{1/2} \exp {(-0.74 \beta _{\|i})}. \end{equation}

For $\beta _{\|i} \gtrsim 4$ , the right-hand side of (2.10) is at least an order of magnitude below unity. We conclude that in a high- $\beta _{\|i}$ plasma with an increasingly negative pressure anisotropy, the resonant oblique firehose instability becomes much faster growing than its parallel counterpart once $\Delta _i \lt \Delta _{\textrm {cr}}$ .

Perhaps more surprisingly, the resonant parallel firehose instability also becomes subdominant to the non-resonant (fluid) parallel firehose instability in $\beta _{\|i} \gg 1$ , bi-Maxwellian plasma at pressure anisotropies not much more negative than the fluid firehose threshold, $\Delta _i = -2/\beta _{\|i}$ . If $\Delta _i \lt -2/\beta _{\|i}$ , then the non-resonant parallel firehose operates at all parallel wavenumbers that satisfy (Schekochihin et al. Reference Schekochihin, Cowley, Rincon and Rosin2010)

(2.11) \begin{equation} k_{\|} \rho _{\|i} \lt 4 \left |\frac {1}{\beta _{\|i}} + \frac {\Delta _i}{2}\right |^{1/2} , \end{equation}

and the peak instability growth rate is

(2.12) \begin{equation} \gamma _{\|\textrm {f,nr}} \equiv \left |\frac {2}{\beta _{\|i}} + \Delta _i\right |\varOmega _{\textrm {i}} \quad \mathrm{at} \quad (k_{\|} \rho _i)_{\|\textrm {f,nr}} \equiv 2 \left |\frac {2}{\beta _{\|i}} + \Delta _i\right |^{1/2}. \end{equation}

If we assume that $|{2}/{\beta _{\|i}} + \Delta _i| \ll 2/\beta _{\|i}$ , it then follows that $\gamma _{\|\textrm {f,nr}} \gtrsim \gamma _{\|\textrm {f}}$ is equivalent to the conditionFootnote ³

(2.13) \begin{equation} |2+\Delta _i \beta _{\|i}| \gg \beta _{\|i}^{1/2} \exp {\left (-\frac {\beta _{\|i}}{2}\right )}. \end{equation}

This bound can be satisfied near marginality of the non-resonant firehose instability provided that $\beta _{\|i} \gtrsim 7$ . Thus, in stark contrast to plasmas with $\beta _{\|i} \sim 1$ (see e.g. Hellinger & Matsumoto Reference Hellinger and Matsumoto2001), the resonant parallel firehose instability is unimportant in $\beta _{\|i} \gg 1$ plasma.

The relative inefficacy of the resonant parallel firehose instability in high- $\beta$ plasma has a simple physical explanation. As was mentioned in § 2.2, the instability is driven by resonant wave–particle interactions: specifically, right-handed circularly polarised hydromagnetic waves drain energy from gyroresonant particles with parallel velocities $v_{\|} =(\omega + \varOmega _i)/k_{\|} \approx v_{\mathrm{th}i}/k_{\|} \rho _i$ . In a plasma with $\beta _i \sim \Delta _i^{-1} \gg 1$ , the gyroresonant particles have characteristic velocities $v_{\|} \sim \Delta _i^{-1/2} v_{\mathrm{th}i}$ that are much greater than the ion thermal velocity. This reveals why the growth rates of the unstable modes are very small: due to their long wavelengths, the hydromagnetic waves can only interact resonantly with suprathermal ions, of which there is only a small number compared with the thermal population. The stabilising action of cyclotron damping is weak on such modes, which in turn allows even a small anisotropy to be able to overcome this damping. However, for shorter-wavelength modes, cyclotron damping is simply too strong for the instability to operate. This conclusion is consistent with the findings of Matteini et al. (Reference Matteini, Landi, Hellinger and Velli2006), who presented evidence of distribution functions becoming less distorted by resonant interactions as $\beta _i$ was increased in one-dimensional expanding-box simulations of firehose-unstable plasma with $\beta _i \leqslant 10$ ; this finding was attributed to the particles that were resonant with parallel firehose modes being increasingly suprathermal.

While the resonant parallel firehose instability is generically unimportant in high- $\beta _i$ plasmas with bi-Maxwellian ion distributions, this conclusion does not necessarily hold for plasmas with non-bi-Maxwellian distributions. Indeed, we will show that right-handed circularly polarised modes can be destabilised by the distribution function that naturally arises during the nonlinear evolution of the oblique firehose instability. These ‘secondary parallel firehose modes’ are characterised and discussed in § 5.2.

3.1. Possible saturated states of the instability

Once firehose modes are linearly destabilised, they grow until they are able to backreact significantly on the pressure anisotropy that drives their growth. Previous analytical and numerical studies suggest that this backreaction causes a transition from exponential growth of the magnetic energy of the modes to secular, power-law growth (Schekochihin et al. Reference Schekochihin, Cowley, Kulsrud, Rosin and Heinemann2008; Rosin et al. Reference Rosin, Schekochihin, Rincon and Cowley2011). In turn, the secular-growth phase eventually transitions into saturation, with the magnetic energy no longer growing. Based on both previous studies (in particular, Melville et al. Reference Melville, Schekochihin and Kunz2016) and the results of this paper, we claim that there are three qualitatively distinct states – ultra-high-beta, Alfvén-inhibiting and Alfvén-enabling – that can be realised by the saturation of firehose instabilities in high- $\beta$ plasmas. Which of these states is realised depends on the relative magnitude of just two independent parameters: $\beta _i$ and $\tau \varOmega _i$ , where we formally define the macroscopic evolution time $\tau$ by

(3.1) \begin{equation} \tau \equiv \left |\left (\hat {{\boldsymbol{b}}}\hat {{\boldsymbol{b}}}-\frac {\unicode{x1D644}}{3}\right )\!\,{\boldsymbol{:}}\, \boldsymbol{\nabla } {\boldsymbol{u}}\right |^{-1} = \left | \frac {\mathrm{d}}{\mathrm{d}t} \log {\frac {B}{n^{2/3}}} \right |^{-1}. \end{equation}

We describe each of these states in §§ 3.1.1–3.1.3, respectively. To aid comparison between these states, table 1 summarises their key properties.

Table 1. Summary of typical values describing the ultra-high-beta, Alfvén-inhibiting and Alfvén-enabling states of a firehose-unstable plasma. Properties include the regulated pressure anisotropy in saturation $\Delta _{\textrm {sat}}$ , the particle-averaged effective collisionality $\nu _{\textrm {eff}}$ , the implied effective Braginskii viscosity $\mu _{\textrm {B,eff}}$ and the characteristic energy $\delta B^2/B_0^2$ of the firehose fluctuations. Whether or not the magnetic-energy spectrum of firehose fluctuations extends to wavelengths much greater than $\rho _i$ is also indicated. We note that $\delta B^2/B_0^2$ in the Alfvén-inhibiting state remains uncertain, because our study and that of Melville et al. (Reference Melville, Schekochihin and Kunz2016) obtain discrepant results (see discussion in § 3.1.2).

The one commonality between all states is the emergence of an effective collisionality $\nu _{\textrm {eff}}$ associated with the firehose fluctuations, which manifests as an additional isotropisation term in the Chew–Goldberger–Low (CGL) equations that describes the evolution of parallel and perpendicular pressures in magnetised plasmas (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956):

(3.2a) \begin{align} \frac {\mathrm{d} p_{\perp }}{\mathrm{d}t} = p_{\perp } \frac {\mathrm{d}}{\mathrm{d}t} \log {(n B)} -\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q}_{\perp } - q_{\perp } \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{b}} - \nu _{\textrm {eff}} (p_{\perp }-p_{\|}),\\[-10pt]\nonumber \end{align}

(3.2b) \begin{align} \frac {\mathrm{d} p_{\|}}{\mathrm{d}t} = p_{\|} \frac {\mathrm{d}}{\mathrm{d}t} \log {\frac {n^3}{B^2}}-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q}_{\|} + 2 q_{\perp } \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{b}} - 2 \nu _{\textrm {eff}} (p_{\|}-p_{\perp }) , \end{align}

where $\boldsymbol{q}_{\|}$ and $\boldsymbol{q}_{\perp }$ denote the parallel heat fluxes of parallel and perpendicular temperature, respectively. This effective collisionality in turn gives rise to an anomalous viscous stress tensor $\boldsymbol{\mathsf{\varPi }}$ . In a weakly collisional plasma ( $\nu _{\textrm {eff}} \ll \varOmega _i$ ), this tensor is approximated well by

(3.3) \begin{equation} \boldsymbol{\mathsf{\varPi }} \approx -\mu _{\textrm {B,eff}} \left (\hat {\boldsymbol{b}} \hat {\boldsymbol{b}}-\frac {\unicode{x1D644}}{3} \right )\left (\hat {\boldsymbol{b}} \hat {\boldsymbol{b}}-\frac {\unicode{x1D644}}{3} \right )\,{\boldsymbol{:}}\,{\boldsymbol{\nabla }} \boldsymbol{u} = \frac {\mu _{\textrm {B,eff}}}{\tau } \left (\hat {\boldsymbol{b}} \hat {\boldsymbol{b}}-\frac {\unicode{x1D644}}{3} \right )\! , \end{equation}

where

(3.4) \begin{equation} \mu _{\textrm {B,eff}} = \frac {p_i}{\nu _{\textrm {eff}}} \simeq -\frac {1}{2} (\Delta _i)_{\textrm {sat}} \beta _{\|i} \frac {B^2}{4 \pi } \,\tau \end{equation}

is the effective Braginskii viscosity.

3.2. Why an Alfvén-enabling state is realised when $\tau \gtrsim \tau _{\textrm {cr}}$

Although the reduced instability threshold for ion-Larmor-scale firehose modes had been identified previously, and while Alfvén-enabling states have been observed in simulations (see e.g. Hellinger & Trávníček Reference Hellinger and Trávníček2008; Bott et al. Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021), we are not aware of any existing theories explaining the physics underpinning the transition in high- $\beta$ plasma between the Alfvén-inhibiting and Alfvén-enabling states. We therefore outline such a theory here, based on our results from § 2.3.

In a plasma in which $\Delta _i$ is driven increasingly negative at a sufficiently slow rate, resonant oblique firehose modes can grow significantly and regulate the pressure anisotropy before $\Delta _i$ becomes negative enough for fluid firehose modes to be destabilised. More specifically, if the characteristic time scale ${\sim }\gamma _{\perp \textrm {f}}^{-1}$ over which the resonant oblique firehose modes grow linearly is much smaller than the time interval $\Delta t$ over which the pressure anisotropy would be driven by the macroscopic evolution from $\Delta _i \simeq -1.35/\beta _{\|i}$ to $\Delta _i = -2/\beta _{\|i}$ (i.e. $\gamma _{\perp \textrm {f}} \Delta t \gg 1$ ), then the growth of these modes will regulate $\Delta _i$ before $\Delta _i$ becomes $\leqslant -2/\beta _{\|i}$ . In this case, an Alfvén-enabling state will persist. If, by contrast, $\gamma _{\perp \textrm {f}} \Delta t \ll 1$ , then resonant oblique firehose modes will not have had the chance to grow before the plasma attains $\Delta _i\leqslant -2/\beta _{\|i}$ and (linear) Alfvén waves no longer propagate.

In the case when $\Delta _i$ is driven linearly in time over a time scale $\tau$ (i.e. $\Delta _i \approx -t/\tau$ , where $t = 0$ is defined as the time at which the pressure is isotropic), the condition for the Alfvén-enabling state to result is $\gamma _{\perp \textrm {f}} \Delta t ={0.65 \gamma _{\perp \textrm {f}} \tau }/{\beta _{i}} \gg 1$ . The transition time scale $\tau _{\textrm {cr}}$ is then the characteristic value of $\tau$ at which $\gamma _{\perp \textrm {f}} \Delta t \approx N_{\textrm {fold}}$ , where $N_{\textrm {fold}}$ is an order-unity factor equal to the number of e-folding times of the instability required for the resonant oblique firehoses to backreact on the plasma (in our simulations, we find $N_{\textrm {fold}} \approx 5)$ . This implies that

(3.5) \begin{equation} \tau _{\textrm {cr}}(\beta _{i}) \approx 1.5 N_{\textrm {fold}} \beta _{i} {\gamma }_{\perp \textrm {f}}^{-1} . \end{equation}

Because ${\gamma }_{\perp \textrm {f}}$ decreases with $\beta _{i}$ , we conclude that $\tau _{\textrm {cr}}$ monotonically increases as $\beta _{i}$ does. Then using the simplified expressions for the growth rate of resonant oblique firehose modes given in § 2.3, explicit expressions for $\tau _{\textrm {cr}}$ as a function of $\beta _{i}$ can be found in various different parameter regimes. When $\beta _ i \gg 10^{6}$ , substituting (2.6) into (3.5) gives

(3.6) \begin{equation} \tau _{\textrm {cr}}(\beta _{i}) \approx 0.9 N_{\textrm {fold}} \beta _{i}^{3/2} \log {\left (\frac {27 \pi \beta _i}{10}\right )} \varOmega _i^{-1}\quad \mathrm{(for \; very \; large \;} \beta _{i} \mathrm{)} . \end{equation}

By contrast, if $\beta _i$ satisfies $1 \ll \beta _{i} \ll 10^{5}$ , then we instead substitute in the empirical scaling (2.7) into (3.5) to obtain

(3.7) \begin{equation} \tau _{\textrm {cr}}(\beta _{\|i}) \approx 5 N_{\textrm {fold}} \beta _{i}^{1.6} \varOmega _i^{-1} \quad \mathrm{(for \;} 1 \ll \beta _i \ll 10^{5}\mathrm{)} . \end{equation}

In both regimes, the transition time scale is much greater than the ion gyroperiod. However, in plasmas where the ratio of the macroscopic evolution time scale $\tau$ to the ion gyroperiod is many orders of magnitude larger than $\beta _i$ – as is often the cases in astrophysical plasmas of interest (see § 7) – both (3.6) and (3.7) imply that $\tau \gt \tau _{\textrm {cr}}$ , with the consequence that such plasmas will always end up in an Alfvén-enabling state.

We can also make specific predictions for the relationship between the parameter $\tau$ and the minimum value $(\Delta _i)_{\textrm {min}}$ of the pressure anisotropy attained when $(\Delta _i)_{\textrm {min}}$ is close to $-2/\beta _{\|i}$ . For example, assuming that the first minimum of $\Delta _i$ is attained when oblique, kinetic-scale firehose fluctuations begin to modify the equilibrium – i.e. when $\gamma _{\perp \textrm {f}} \Delta t \approx N_{\textrm {fold}}$ – it follows that $(\Delta _i)_{\textrm {min}} \approx -1.35/\beta _i - N_{\textrm {fold}}/(\gamma _{\perp \textrm {f}}\tau )$ . Then, in the case of moderately large $\beta _i$ ( $1 \ll \beta _{i} \ll 10^{5}$ ), (2.8) for the peak growth rate of oblique firehose modes when $(\Delta _i)_{\textrm {min}}$ is close to $-2/\beta _{\|i}$ implies that

(3.8) \begin{equation} (\Delta _i)_{\textrm {min}} \approx \Delta _{\textrm {cr}} - 1.8 \left (\frac {N_{\textrm {fold}}}{\tau \varOmega _i}\right )^{0.625} . \end{equation}

This states that the difference $(\Delta _i)_{\textrm {min}} - \Delta _{\textrm {cr}}$ does not depend on $\beta _{\|i}$ in this parameter regime, instead being proportional to $(\tau \varOmega _i)^{-0.625}$ . Another corollary of (3.8) is that $(\Delta _i)_{\textrm {min}} \beta _i$ is not a function of $\tau$ , $\varOmega _i$ and $\beta _{i}$ independently, but rather only of the specific combination $\tau /\tau _{\textrm {cr}}$ .

3.3. Summary and the rest of this paper

In this section, we have summarised the possible states that can be realised in firehose-saturated, high- $\beta$ plasmas. These claims obviously require careful justification with recourse to nonlinear analytical studies and/or simulations. While such a study on the transition between ultra-high-beta and Alfvén-inhibiting states has already been completed by Melville et al. (Reference Melville, Schekochihin and Kunz2016), no prior study has been done of the analogous transition between the Alfvén-inhibiting and Alfvén-enabling states. We have carried out such a study and report its results in §§ 4–6. Readers who are happy to take such results on trust, and are instead keen to consider the relationship of our theory of firehose saturation with previous studies and the implications for astrophysical high- $\beta$ plasmas, are encouraged to skip forward to § 7.

4.1. Overview

While the existence of Alfvén-enabling and Alfvén-inhibiting states in firehose-susceptible plasmas can be predicted via the linear theory of the firehose instability, determining the equilibrium properties of these two states necessitates modelling the firehose’s nonlinear saturation across a range of different parameters (e.g. $\beta _{\|i}$ , $\tau$ ) as pressure anisotropy is driven by a plasma’s macroscopic evolution. This is most effectively done numerically. In this section, we first explain why so-called HEB simulations are particularly well suited to this purpose, and describe the method underpinning them. Then, we outline the results of a parameter study of numerous such simulations, characterising the time evolution of quantities such as the pressure anisotropy $\Delta _i$ , the effective Alfvén speed $v_{\textrm {A,eff}}$ and the magnetic-field strength $\delta B_{\textrm {f}}$ of the firehose fluctuations. This, in turn, allows us to determine the equilibrium thermodynamic and microphysical properties of the Alfvén-enabling and Alfvén-inhibiting states, respectively.

Our key finding is that, in expanding, high- $\beta _i$ plasmas with an expansion time $\tau _{\textrm {exp}}$ (see § 4.2.1) that satisfies $\tau _{\textrm {exp}} \gg \tau _{\textrm {cr}}(\beta _{\|i})$ (i.e the ‘asymptotic’ Alfvén-enabling state), the pressure anisotropy is regulated to a value $\Delta _i \simeq -1.6/\beta _{\|i}$ that is above the value $\Delta _i = -2/\beta _{\|i}$ at which Alfvén waves cease to propagate. By contrast, if $\beta _i \varOmega _i^{-1} \ll \tau _{\textrm {exp}} \ll \tau _{\textrm {cr}}(\beta _{\|i})$ (i.e. an Alfvén-inhibiting state), then $\Delta _i \simeq -2/\beta _{\|i}$ . We also show that the firehose fluctuations are qualitatively distinct in the two regimes. In the Alfvén-inhibiting state, a broad spectrum of magnetic fluctuations (including long-wavelength modes) is excited; in the Alfvén-enabling state, magnetic energy is primarily concentrated in fluctuations at ion-Larmor scales. In the latter case, there are two types of modes: oblique modes and parallel ion-Larmor-scale modes. The latter are not, in fact, resonant parallel firehose modes of the conventional type, but are instead a secondary instability associated with the (non-bi-Maxwellian) ion distribution that is created by resonant scattering of suprathermal ions by the oblique firehose modes. That the saturated value of $\Delta _i$ is somewhat more negative ( ${\simeq}-1.6/\beta _{\|i}$ ) than the linear threshold of the resonant oblique firehose instability in a bi-Maxwellian plasma ( ${\simeq}-1.35/\beta _{\|i}$ ) can also be attributed to the non-bi-Maxwellian form of the ion distribution function in saturation.

Finally, we characterise the velocity-space-averaged effective collisionality $\nu _{\textrm {eff}}$ for firehose-susceptible plasmas in both the Alfvén-enabling and the Alfvén-inhibiting states (see § 4.4). We confirm that, for all of our expanding-box simulations, $\nu _{\textrm {eff}} \sim \beta _{i}/\tau _{\textrm {exp}}$ , in agreement with previous shearing-box simulations of the firehose instability that are not in the ultra-high-beta regime (Kunz et al. Reference Kunz, Schekochihin and Stone2014a ; Melville et al. Reference Melville, Schekochihin and Kunz2016). We then provide quantitative estimates of the plasma’s effective parallel Braginskii viscosity $\mu _{\textrm {B}}$ in our simulations.

4.2. Simulation set-up

4.2.1. Why simulate an expanding plasma?

As was mentioned in the Introduction, a range of different macroscopic bulk-flow fluid motions – including shearing and expanding motions – can give rise to negative ion-pressure anisotropy ( $\Delta _i \lt 0$ ) in collisionless, magnetised plasma. To see this in more detail, let us assume that the parallel and perpendicular pressures evolve according to the double-adiabatic equations (i.e. the CGL equations (3.2b ) after dropping the heat fluxes and effective collisionality):

(4.1) \begin{equation} \frac {\mathrm{d}}{\mathrm{d}t} \log {\frac {p_{\perp }}{nB}} = 0 , \quad \frac {\mathrm{d}}{\mathrm{d}t} \log {\frac {p_{\|} B^2}{n^3}} = 0. \end{equation}

In any collisionless plasma governed by these equations whose initial temperature is isotropic (i.e. $T_{\|} = T_{\perp }$ at some time $t = 0$ ), the pressure anisotropy satisfies

(4.2) \begin{equation} \Delta _i = \frac {B^3/B_0^3}{n^2/n_0^2} - 1 , \end{equation}

where the subscript ‘ $0$ ’ is from here on used to denote the value of a quantity at $t = 0$ . Thus, the ion-pressure anisotropy will decrease in any double-adiabatic plasma in which $B^3/n^2$ decreases due to the plasma’s macroscopic evolution.

Of the motions that cause $B^3/n^2$ to decrease, a particularly advantageous one to choose for our purposes is that of spatially uniform expansion at a constant rate in one direction transverse to the mean magnetic field. There are a few different physical situations in which this type of expansion could arise: during the motion of a macroscopic, linearly polarised magnetosonic wave travelling perpendicularly to the background magnetic field; in certain regions of compressive turbulence; and, finally, the expansion of cylindrical, magnetised plasma – for example, generated by an exploding wire array in a laboratory astrophysics experiment. Assuming that expansion occurs at a rate $1/\tau _{\textrm {exp}}$ , where $\tau _{\textrm {exp}}$ is the expansion time, it follows that $B \propto n = n_0/(1+t/\tau _{\textrm {exp}})$ , and so

(4.3) \begin{equation} \Delta _i = \frac {n}{n_0} - 1 = -\frac {t}{t+\tau _{\textrm {exp}}}. \end{equation}

Beyond studying specific physical systems, there are three pragmatic reasons for this choice of motion in order to study firehose instabilities. First of these is the possibility of simulating such an expansion exactly via a coordinate transform method, allowing for a simulation domain that is both fixed and homogeneous to be used. Using a coordinate transform method (of which a shearing box is another example) maximises the effective separation between macro- and microscales for a fixed simulation domain size; it also allows for simulation-domain-averaged properties of the plasma (including the ion distribution function) to be used as a reasonable analogue for that plasma’s ‘equilibrium’ properties, minimising the uncertainty that could be introduced by macroscopic spatial variation of the plasma. The second reason is that, in contrast to a shearing-box simulation (e.g. Kunz et al. Reference Kunz, Schekochihin and Stone2014a ), the direction of the macroscopic magnetic field does not change as the motion proceeds, which simplifies comparing different times in the simulation. Finally, compared with other ‘simple’ motions, an expansion in a direction transverse to the background magnetic field gives rise to a comparably slow evolution of the pressure anisotropy over a fixed period of time. For example, a two-dimensional incompressible motion in which there is simultaneously expansion in one direction transverse to the background magnetic field and contraction in the parallel direction, causing the background magnetic-field strength to vary as $B = B_0/(1+t/\tau _{\textrm {exp}})$ , would give rise to a value of $|{\mathrm{d} \Delta _i}/{\mathrm{d}t}|$ that is initially three times larger than the analogous one-dimensional transverse expansion. As we show in § 4.3.1, accessing the Alfvén-enabling regime when $\beta _{i} \gg 1$ requires macroscopic evolution rates that are at least several orders of magnitude smaller than the ion-Larmor frequency $\varOmega _i$ ; because such simulations are expensive, choosing a motion that minimises the rate of change of $\Delta _i$ at a fixed time period is desirable. Motivated by these considerations, we choose in this paper to simulate plasmas expanding in a single transverse direction. In addition to its application to the specific physical systems mentioned earlier in this paragraph, we anticipate that the evolution and saturation of the firehose instability becomes insensitive to the specifics of the macroscopic motion driving it provided there is sufficient separation of relevant time scales; we revisit this assumption after describing our simulation results in § 7.

4.2.2. Hybrid expanding box simulations with Pegasus++

To capture all ion firehose instabilities correctly, the plasma’s ions (but not necessarily the electrons) must be modelled kinetically. We therefore choose to conduct HEB simulations. Although this approach and its implementation have been described elsewhere (e.g. Hellinger & Trávníček Reference Hellinger and Trávníček2005; Bott et al. Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021), we explain the method here for completeness. All HEB simulations reported in this paper were carried out using the PIC code Pegasus++ (Kunz et al. Reference Kunz, Stone and Bai2014b ; Bott et al. Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021).

In the HEB approach, one transforms from the locally co-moving frame of the expanding plasma to a co-moving frame in which the metric extends as the plasma expands, and then performs all subsequent calculations in this expanding frame. Denoting position in the co-moving, non-expanding frame by $\boldsymbol{r}$ , and in the co-moving, expanding frame by ${\boldsymbol{r}}'$ , the frame transformation is characterised by a matrix $\boldsymbol{\mathsf{\varLambda }}\equiv \partial {\boldsymbol{r}}/\partial {\boldsymbol{r}}'$ , with determinant $\lambda = \mathrm{det} \, \boldsymbol{\mathsf{\varLambda }}$ . For HEB PIC simulations using Pegasus++, we evolve two sets of equations: those describing the motion of ion macroparticles and those describing the evolution of electromagnetic fields. The former, which constitute evolution equations for the primed-frame positions ${\boldsymbol{r}}'_p=\boldsymbol{\mathsf{\varLambda }}^{-1}{\boldsymbol{r}}_p$ and velocities ${\boldsymbol{v}}'_p = \boldsymbol{\mathsf{\varLambda }}^{-1} {\boldsymbol{v}}_p$ of macroparticles, are given by

(4.4a) \begin{align} \frac {\mathrm{d} {\boldsymbol{r}}'_p}{\mathrm{d}t'} & = {\boldsymbol{v}}'_p ,\\[-10pt]\nonumber \end{align}

(4.4b) \begin{align} \frac {\mathrm{d}{\boldsymbol{v}}'_p}{\mathrm{d} t'} & = \frac {e}{m_i} \, \boldsymbol{\mathsf{\varLambda }}^{-2} \left [ {\boldsymbol{E}}'(t',{\boldsymbol{r}}'_p) + \frac {{\boldsymbol{v}}'_p}{c}\,{\boldsymbol{\times }}\,{\boldsymbol{B}}'(t',{\boldsymbol{r}}'_p)\right ] - 2 \boldsymbol{\mathsf{\varLambda }}^{-1}\frac {\mathrm{d}\boldsymbol{\mathsf{\varLambda }}}{\mathrm{d} t'}\,{\boldsymbol{v}}'_p , \end{align}

where the fields ${\boldsymbol{E}}'$ and ${\boldsymbol{B}}'$ are related to the physical electric field $\boldsymbol{E}$ and magnetic field $\boldsymbol{B}$ in the unprimed frame viaFootnote ⁴

(4.5) \begin{equation} {\boldsymbol{B}}' = \lambda \boldsymbol{\mathsf{\varLambda }}^{-1} {\boldsymbol{B}} \quad \textrm {and}\quad {\boldsymbol{E}}' = \boldsymbol{\mathsf{\varLambda }} {\boldsymbol{E}} . \end{equation}

To solve (4.4b ), Pegasus++employs a straightforward modification of the Boris push that groups the velocity-dependent non-inertial force with the ${\boldsymbol{v}}'_p \,{\boldsymbol{\times }}\,{\boldsymbol{B}}$ rotation. The fields ${\boldsymbol{B}}'$ and ${\boldsymbol{E}}'$ in turn satisfy modified versions of Faraday’s law and a generalised Ohm’s law, respectively:Footnote ⁵

(4.6a) \begin{align} \frac {\partial {\boldsymbol{B}}'}{\partial t'} & = -c\boldsymbol{\nabla }'\,{\boldsymbol{\times }}\,{\boldsymbol{E}}' , \end{align}

(4.6b) \begin{align} {\boldsymbol{E}}' & = -\frac {{\boldsymbol{u}}'}{c}\,{\boldsymbol{\times }}\,{\boldsymbol{B}}' - \frac {T_e}{en'}\boldsymbol{\nabla }' n' + [\boldsymbol{\nabla }'\,{\boldsymbol{\times }}\, (\boldsymbol{\mathsf{\varLambda }}^2{\boldsymbol{B}}')]\,{\boldsymbol{\times }}\,\frac {{\boldsymbol{B}}'}{4\pi e n'\lambda }. \end{align}

Fluid quantities in the primed frame are calculated in the usual way by taking moments of the primed-frame ion distribution function. This is done by summing up the (weighted) phase-space contributions from each ion macroparticle of shape $S$ centred on the phase-space position $({\boldsymbol{r}}'_p,{\boldsymbol{v}}'_p)$ to the phase-space position ( ${\boldsymbol{r}}',{\boldsymbol{v}}'$ ). For example, the primed-frame ion density $n'$ and primed-frame ion-flow velocity ${\boldsymbol{u}}'$ , which are related to their unprimed analogues via $n' = \lambda n$ and ${\boldsymbol{u}}' = \boldsymbol{\mathsf{\varLambda }}^{-1} {\boldsymbol{u}}$ , are computed via

(4.7) \begin{equation} n'({\boldsymbol{r}}') = \sum \nolimits _p \, S\big[{\boldsymbol{r}}'-{\boldsymbol{r}}'_p(t')\big] \quad \textrm {and}\quad {\boldsymbol{u}}'({\boldsymbol{r}}') = \frac {1}{n'} \sum \nolimits _p \, {\boldsymbol{v}}'_p \, S\big[{\boldsymbol{r}}'-{\boldsymbol{r}}'_p(t')\big] . \end{equation}

At any given time, physical variables can be computed directly from the primed-frame variables using the appropriate inverse coordinate transform.

For our simulations, we set

(4.8) \begin{equation} \boldsymbol{\mathsf{\varLambda }}(t) = \left (1+\frac {t}{\tau _{\textrm {exp}}}\right )\hat {{\boldsymbol{x}}} \hat {{\boldsymbol{x}}} + \hat {{\boldsymbol{y}}} \hat {{\boldsymbol{y}}} + \hat {{\boldsymbol{z}}} \hat {{\boldsymbol{z}}} , \end{equation}

where $\{\hat {{\boldsymbol{x}}},\hat {{\boldsymbol{y}}},\hat {{\boldsymbol{z}}}\}$ is a set of basis vectors of an orthogonal coordinate system in which $\hat {{\boldsymbol{z}}}$ is parallel to ${\boldsymbol{B}}(t = 0)$ (and remains parallel to the mean ‘guide’ field in the simulation domain throughout). In terms of the evolution of the side lengths $[L_x(t),L_z(t)]$ of a two-dimensional spatial domain, (4.8) gives

(4.9) \begin{equation} L_x(t) = L_{x0} \left (1+\frac {t}{\tau _{\textrm {exp}}}\right ) , \quad L_z(t) = L_{z0} . \end{equation}

We define the effective expansion time $\tau _{\textrm {exp,eff}}$ via

(4.10) \begin{equation} \tau _{\textrm {exp,eff}} \equiv \left (\frac {\textrm {d} }{\textrm {d} t}\log {\frac {B^3}{n^2}}\right )^{-1}= \left (\frac {\textrm {d} }{\textrm {d} t}\log {B}\right )^{-1} = t + \tau _{\textrm {exp}}. \end{equation}

We choose this definition for three reasons. Firstly, in the limit of small pressure anisotropy, it is a simple matter to show that $\mathrm{d} \Delta _i/{\mathrm{d} t} \approx \mathrm{d} (\Delta _i \beta _{\|i})/{\mathrm{d} t} \approx -1/\tau _{\textrm {exp,eff}}$ as the plasma expands. Secondly, in the saturated phase of the firehose instability, it can be shown that the box-averaged effective collisionality associated with firehose fluctuations is inversely proportional to $\tau _{\textrm {exp,eff}}$ (see § 4.4). Finally, although $\tau _{\textrm {exp,eff}}$ increases in time, the ion-cyclotron frequency $\varOmega _i \propto B = B_0 \tau _{\textrm {exp}}/\tau _{\textrm {exp,eff}}$ decreases in such a way that their product is constant: $\tau _{\textrm {exp,eff}} \varOmega _i = \tau _{\textrm {exp}} \varOmega _{i0}$ .

We ran numerous HEB simulations of this type with different values of $\tau _{\textrm {exp}}$ and $\beta _{i0}$ . We chose to perform these simulations in a 2.5-dimensional geometry – that is, particles move in three dimensions and the electromagnetic fields are three-dimensional, but spatial gradients are restricted to the two-dimensional $(x,z)$ plane – because such simulations can capture the relevant physics at significantly reduced computational costs. Periodic boundary conditions were applied in all spatial directions. Table 2 outlines the key parameters of all of the simulations reported in the paper. All simulations were initialised with equal parallel and perpendicular temperatures. The numerical resolution of the simulations was chosen ( $\Delta x = \Delta z = 0.26 \rho _i$ ) such that the characteristic wavenumbers of firehose modes would be sufficiently well resolved: the maximum wavenumber $k_{\textrm {max}}$ of modes that could be resolved in our simulations was $k_{\textrm {max}} = \pi /\Delta z \simeq 12.1 \rho _i^{-1}$ . We note that, as a result of this choice, for those of our simulations with $\beta _{i0} = 200$ , the wavenumber $k = d_i^{-1}$ at which fluctuations on the scale of ion skin depth $d_i$ exist was not resolved; however, we believe that this is an acceptable limitation, because in high- $\beta$ plasmas, the characteristic scale of firehose fluctuations is $\rho _i$ , not $d_i$ (see e.g. Bott et al. Reference Bott, Cowley and Schekochihin2024). Simulations were run until what appeared to be a saturated state was obtained. The large number $N_{\textrm {ppc}}$ of particles per cell used in these simulations is necessary in order to suppress the effect of numerical collisions (arising from the Poisson noise due to finite sampling of the ion distribution function) on both the evolution of the pressure anisotropy and the firehose instability itself. Even with such large values of $N_{\textrm {ppc}}$ , we find that numerical collisionality has a quantitative (but not a qualitative) effect on some of our results (see Appendix B).

Table 2. Parameters of all HEB simulations performed in this study. Here, $t_{\textrm {min}}$ is the time at which the (first) minimum of the pressure anisotropy is attained, $t_{\textrm {end}}$ is the time at which the simulation run was ended and $\tilde {\tau }_0$ and $\tilde {\tau }_{\textrm {eff}}$ are defined by $\tilde {\tau }_0 \equiv \tau _{\textrm {exp}} \varOmega _{i0}\beta _{i0}^{-1.6}/27$ and $\tilde {\tau }_{\textrm {eff}} \equiv \tau _{\textrm {exp,eff}} \varOmega _i({t_{\textrm {min}}})\beta _{i}({t_{\textrm {min}}})^{-1.6}/27$ , respectively. The empirical factor of $27$ is introduced so that runs with $\tilde {\tau }_{\textrm {eff}} \gtrsim 1$ are at all times in an Alfvén-enabling state (see § 4.3.1). In units of $\rho _i$ , all simulations were run with the same numerical resolution ( $\Delta x = \Delta z = 0.26 \rho _i$ ), and with the same initial side lengths of the simulation domain ( $L_{z0} = 1.5 L_{x0} = 300 \rho _i$ ).

4.3. Results

4.3.1. Box-averaged pressure anisotropy and effective Alfvén speed

The existence of both Alfvén-enabling and Alfvén-inhibiting states among the expanding plasmas we have simulated is illustrated in figure 3( $a$ ), which positions each simulation in a $[-\varOmega _i (\mathrm{d}\Delta _i/\mathrm{d}t)|_{t = 0})^{-1},\beta _{i0}] = [\tau _{\textrm {exp}}\varOmega _{i0},\beta _{i0}]$ phase space, and indicates whether the value of the ion-pressure anisotropy $(\Delta _i)_{\textrm {min}}$ attained at the first minimum (at some time $t = t_{\textrm {min}}$ ) is more (red) or less (blue) negative than $\Delta _i = -2/\beta _{\|i}$ . Qualitatively, at fixed $\beta _{i0}$ the simulated plasmas transition from being in an Alfvén-enabling state (for which $(\Delta _i)_{\textrm {min}} \gt -2/\beta _{\|i}$ ) to an Alfvén-inhibiting state (for which $(\Delta _i)_{\textrm {min}} \ll -2/\beta _{\|i}$ ) as the expansion time is decreased. Furthermore, the characteristic expansion time at which this transition occurs is (for the suite of simulations we have conducted) a monotonically increasing function of $\beta _{\|i}$ . More quantitatively, figure 3( $b$ ) shows that the scaling (3.7) of $\tau _{\textrm {cr}}$ with $\beta _{\|i}$ derived in § 2 (dashed line) is an excellent fit to the measured (effective) expansion time $\tau _{\textrm {exp,eff}}$ at which the instantaneous values of $\beta _i$ and $(\Delta _i)_{\textrm {min}}$ satisfy $(\Delta _i)_{\textrm {min}} \approx -2/\beta _{\|i}$ . For reference, in figure 3( $a$ ) we plot also the positions in the same $[-\varOmega _i (\mathrm{d}\Delta _i/\mathrm{d}t)|_{t = 0})^{-1},\beta _{i0}]$ phase space of previously published high- $\beta _i$ shearing-box simulations of the firehose instability (Kunz et al. Reference Kunz, Schekochihin and Stone2014a ; Melville et al. Reference Melville, Schekochihin and Kunz2016); these simulations all realised the Alfvén-inhibiting state in saturation, a finding consistent with their initialised parameters.

Figure 3. Phase-space maps of various simulations of high- $\beta _{i}$ firehose-susceptible plasmas, which indicate whether an Alfvén-enabling state ( $\Delta _i \gt -2/\beta _{\|i}$ ) is maintained at all times (blue points) or not (red points). In (a), we include the HEB simulations completed for this paper (denoted by ‘ $\times$ ’), as well as the shearing-box hybrid-kinetic simulations reported in Kunz et al. (Reference Kunz, Schekochihin and Stone2014a ) (‘ $+$ ’) and Melville et al. (Reference Melville, Schekochihin and Kunz2016) (‘ $\,{\boldsymbol{\cdot }}\,$ ’). The dashed line that provides an accurate delineation of consistently Alfvén-enabling states and other states is given by the equation $\tau _{\textrm {exp}} \varOmega _i = 27 \beta _i^{1.6}$ (cf. figure 1 and (3.7)).

A simple way to illustrate how the evolution of the pressure anisotropy and the effective Alfvén speed changes as the expanding plasma transitions from being Alfvén-inhibiting to Alfvén-enabling is to fix $\beta _{i0}$ , and compare the evolution of $\Delta _i \beta _{\|i}$ and $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ over time for a selection of increasing values of $\tau _{\textrm {exp}}$ . This comparison is made in figure 4. It is clear from figure 4 $(a)$ that, for $\Delta _i \gt \Delta _{\textrm {cr}} \simeq -1.35/\beta _{\|i}$ , the initial evolution of $\Delta _i \beta _{\|i}$ is independent of $\tau _{\textrm {exp}}$ , as predicted by (4.3).Footnote ⁶ However, once the oblique firehose instability is triggered, we see that for comparatively larger expansion times (e.g. black and red lines in figure 4 $a$ ), $-\Delta _i \beta _{\|i}$ stops increasing at smaller characteristic values of $t/\tau _{\textrm {exp}}$ , and attains a less positive value $-\Delta _i \beta _{\|i}$ at the time $t_{\textrm {min}}$ at which the first minimum of $\Delta _i$ , $(\Delta _i)_{\textrm {min}}$ , is attained. At times $t \gt t_{\textrm {min}}$ , $-\Delta _i \beta _{\|i}$ is regulated, eventually converging to an order-unity value in all of our simulations. For the largest values of $\tau _{\textrm {exp}}$ , we find that the pressure anisotropy is ultimately regulated to values $\Delta _i \approx -1.6/\beta _{\|i}$ . By contrast, for the comparatively smaller expansion times (e.g. cyan and blue lines), $\Delta _i \approx -2/\beta _{\|i}$ , characteristic of an Alfvén-inhibiting state. These saturated values of $\Delta _i$ imply that, for the simulations with comparatively smaller expansion times that we have run, the plasma attains an Alfvén-inhibiting state with $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2 \approx 0$ , while for the larger expansion times, $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2 \approx 0.2$ (see figure 4 $b$ ). For intermediate values of $\tau _{\textrm {exp}}$ , $(\Delta _i)_{\textrm {min}}$ drops below $\Delta _i = -2/\beta _{\|i}$ by an $\textit {O}(1/\beta _{\|i})$ value, but the ‘steady-state’ values of $\Delta _i$ that are subsequently attained imply that the state in these runs is, in saturation, Alfvén-enabling (albeit with a reduced value of $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ compared with runs in which $(\Delta _i)_{\textrm {min}} \gt -2/\beta _{\|i}$ ).

Figure 4. Time evolution of $(a)$ the firehose-instability parameter $-\Delta _i \beta _{\|i}$ and $(b)$ the squared and normalised effective Alfvén speed $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ for all of the D runs ( $\beta _{i0} = 50$ ). In $(a)$ , the dotted black line denotes the threshold $\Delta _i \beta _{\|i} = -1.35$ of the oblique firehose instability when $\beta _i \gg 1$ , while the dashed black line shows the fluid firehose threshold $\Delta _i \beta _{\|i} = -2$ . In $(b)$ , the dotted (dashed) black line denotes the corresponding value $v_{\textrm {A,eff}}^2 = 0.32 v_{\textrm {A}}^2$ ( $v_{\textrm {A,eff}}^2 = 0 v_{\textrm {A}}^2$ ) of the squared effective Alfvén speed at the threshold of the oblique firehose instability (fluid firehose threshold).

Figure 5. $(a)$ Values of the firehose-instability parameter $-\Delta _i \beta _{\|i}$ at the time $t_{\textrm {min}}$ at which the pressure anisotropy attains its first minimum, $(\Delta _i)_{\textrm {min}}$ , for all runs, as a function of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {min}})^{1.6}$ . The dotted (dashed) black line denotes the threshold $\Delta _i \beta _{\|i} = -1.35$ of the oblique firehose instability (the fluid firehose instability threshold, $\Delta _i \beta _{\|i} = -2$ ) when $\beta _i \gg 1$ ; the dotted-dashed purple line denotes (4.11). $(b)$ Values of the difference between $(\Delta _i)_{\textrm {min}}$ and the value $\Delta _{\textrm {cr}}$ at which the oblique firehose becomes unstable as a function of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {min}})^{1.6}$ , for all runs. The dotted-dashed purple line denotes (4.11).

Figure 6. $(a)$ Values of the firehose-instability parameter $-\Delta _i \beta _{\|i}$ at the time $t_{\textrm {sat}}$ at which the square of the perturbed magnetic-field strength $\delta B_{\textrm {f}}^2/B_0^2$ associated with the firehose fluctuations attains its maximum value, $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}}$ , for all runs as a function of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6}$ . The dashed black line denotes the fluid firehose instability threshold, $\Delta _i \beta _{\|i} = -2$ , when $\beta _i \gg 1$ . $(b)$ Values of the square of the effective Alfvén speed, $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ at $t = t_{\textrm {sat}}$ , for all runs as a function of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {min}})^{1.6}$ .

A key prediction of the theory outlined in § 2 is that the transition between Alfvén-enabling and Alfvén-inhibiting states in firehose-susceptible high- $\beta _{i}$ plasmas is a function of the parameter $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _i^{1.6}$ (in the limit where $1 \ll \beta _i \ll 10^{5}$ ). We test this prediction in figure 5 $(a)$ by plotting for each of our simulations the relationship between $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _i^{1.6}$ and $(\Delta _i)_{\textrm {min}} \beta _{\|i}(t_{\textrm {min}})$ . We see that the value of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _i^{1.6}$ is predictive of $(\Delta _i)_{\textrm {min}} \beta _{\|i}(t_{\textrm {min}})$ for all of our simulations, with the decreasing nonlinear relationship

(4.11) \begin{equation} (\Delta _i)_{\textrm {min}} \beta _{\|i}(t_{\textrm {min}}) \approx -1.35 - 5.1 \frac { \beta _{\|i}(t_{\textrm {min}})}{\left (\tau _{\textrm {exp,eff}} \varOmega _i\right )^{0.625}} \end{equation}

between the two parameters being a good fit to our data. This relationship is consistent with the prediction (3.8) that was based on the linear theory of the firehose instability (with $N_{\textrm {fold}} \simeq 5.4$ ). It follows that $(\Delta _i)_{\textrm {min}} \beta _{\|i}(t_{\textrm {min}}) \approx -2$ when $\tau _{\textrm {exp,eff}} \varOmega _i \approx 27 \beta _i^{1.6}$ . Furthermore, figure 5 $(b)$ shows that the power-law dependence of $(\Delta _i)_{\textrm {min}} - \Delta _{\textrm {cr}}$ on $(\tau _{\textrm {exp,eff}} \varOmega _i)^{-0.625}$ that was predicted by (3.8) is well satisfied.

The parameter $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _i^{1.6}$ also has a quasi-deterministic relationship with the values of $\Delta _i \beta _{\|i}$ and $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ in our simulations once the firehose instability has saturated. We illustrate this in figure 6 by plotting $-\Delta _i \beta _{\|i}$ and $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ at the time $t_{\textrm {sat}}$ at which the firehose fluctuations attain their peak magnetic-field strength; we denote the value of $\Delta _i$ attained at this time as $(\Delta _i)_{\textrm {sat}}$ . Figure 6 $(a)$ shows that, as $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6}$ increases from below unity to much greater values, $(\Delta _i)_{\textrm {sat}}\beta _{\|i}(t_{\textrm {sat}})$ increases monotonically from a value close to $-2$ to a less negative value of $\simeq$ $-1.6$ ; equivalently, $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ increases from being close to zero to $\simeq$ $0.2$ . For the simulations we have performed, we find that for $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6} \gtrsim 80$ , $(\Delta _i)_{\textrm {sat}}\beta _{\|i}(t_{\textrm {sat}})$ does not become less negative if $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6}$ is increased further still (and $v_{\textrm {A,eff}}^2/v_{\textrm {A}}^2$ does not increase). We infer from this that such a state is the ‘asymptotic’ Alfvén-enabling state for asymptotically large values of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6}$ . Given that the relevance of the parameter $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _i^{1.6}$ is derived entirely from the linear theory of the firehose instability, it is perhaps unsurprising that the correlation between $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6}$ and $(\Delta _i)_{\textrm {sat}}\beta _{\|i}(t_{\textrm {sat}})$ is indeed less strong than that between $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {min}})^{1.6}$ and $(\Delta _i)_{\textrm {min}}\beta _{\|i}(t_{\textrm {min}})$ ; however, the existence of any correlation at all suggests that the initial evolution of the firehose instability has a qualitative effect on the subsequent dynamics. Furthermore, the spread in values of $(\Delta _i)_{\textrm {sat}}\beta _{\|i}(t_{\textrm {sat}})$ at particular values of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6}$ is partially explained by the fact that, for plasmas in an Alfvén-enabling state, $\Delta _i\beta _{\|i}$ periodically fluctuates once the firehose instability has saturated; our chosen measure of $\Delta _i\beta _{\|i}$ in saturation is pointwise in time, and so does not account for this effect. Comparison with time-dependent phase-space plots of $[\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6},(\Delta _i)_{\textrm {sat}}\beta _{\|i}(t_{\textrm {sat}})]$ (not shown) supports this explanation, and also recovers the same general trend that is observed in figure 6 $(a)$ .

In summary, our simulation results confirm that the parameter $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _i^{1.6}$ is indeed a key metric for determining whether a firehose-susceptible high- $\beta _i$ plasma attains an Alfvén-inhibiting or Alfvén-enabling state once the firehose instability has saturated.

Figure 7. Two-dimensional visualisations of the out-of-plane component of the perturbed magnetic field in two simulations at $\beta _{i0} = 50$ representing $(a)$ an Alfvén-enabling state (run DVI) and $(b)$ an Alfvén-inhibiting state (run DIII). The perturbed field appears at three different times: in the linear, nonlinear and saturated phases of the firehose instability, respectively. We note that, as a fraction of the expansion time, the characteristic times at which the linear, nonlinear and saturated states are realised are longer in the Alfvén-inhibiting than Alfvén-enabling regime; this is because the firehose instability develops at a comparatively slower rate in this case when compared with the expansion rate.

4.3.2. Magnetic-field fluctuations

In addition to having distinct macroscopic properties – specifically, different equilibrium pressure anisotropies and effective Alfvén speeds – the Alfvén-enabling and Alfvén-inhibiting states are different microphysically. One manifestation of this is the nature of the firehose fluctuations that arise. Figure 7 visualises the out-of-plane (dominant) component of the perturbed magnetic field for two simulations having $\beta _{i0}=50$ but differing $\tau _{\textrm {exp}}\varOmega _{i0}$ , such that one realises an Alfvén-enabling state (figure 7 $a$ ; run DVI) while the other realises an Alfvén-inhibiting state (figure 7 $b$ ; run DIII). Initially, in both simulations oblique firehose fluctuations with characteristic wavenumber $k_{\|} \rho _i \sim k_{\perp } \rho _i \approx 0.45$ are destabilised. However, the magnitude of the magnetic-field perturbations in both the nonlinear regime and the saturated states is larger in the simulation that realises an Alfvén-inhibiting state (relative to an Alfvén-enabling state). In addition, oblique fluctuations occurring over a range of scales are much more prominent in the Alfvén-inhibiting state.

How the key parameters of the expanding plasma affect the characteristic amplitude of magnetic fluctuations can be most simply explored by considering the evolution of the box-averaged perturbed magnetic energy, $\delta B_{\textrm {f}}^2/B_0^2$ . Figure 8 $(a)$ shows the evolution of $\delta B_{\textrm {f}}^2/B_0^2$ in time at fixed $\beta _{i0}$ . The evolution of $\delta B_{\textrm {f}}^2/B_0^2$ in all of our simulations proceeds through four phases. First, there is a pre-firehose phase, in which the box-averaged magnetic-field strength of the fluctuations is simply that associated with random grid-scale fluctuations; next, a linear growth stage, during which the amplitude of firehose fluctuations grows exponentially; third, a nonlinear phase, in which the amplitude of fluctuations continue to grow, but no longer exponentially; finally, saturation. How $\delta B_{\textrm {f}}^2/B_0^2$ evolves qualitatively as a function of time during the nonlinear and saturated phase of the firehose depends on $\beta _{\|i}$ and $\tau _{\textrm {exp}} \varOmega _i$ . At sufficiently large values of $\tau _{\textrm {exp}}$ (at fixed $\beta _{\|i}$ ), $\delta B_{\textrm {f}}^2/B_0^2$ does not grow monotonically during the nonlinear phase, nor is it constant in the ‘saturated’ state (see especially the black line in figure 8 $a$ ). Instead, the magnetic energy oscillates around a mean value with a characteristic period that is much smaller than $\tau _{\textrm {exp,eff}}$ . These oscillations correlate with those seen in the pressure anisotropy in § 4.3.1, implying a direct link between the amplitude of the firehose fluctuations in saturation, and the regulation of the pressure anisotropy. For smaller values of $\tau _{\textrm {exp}}$ (again at fixed $\beta _{\|i}$ ), $\delta B_{\textrm {f}}^2/B_0^2$ does not oscillate in saturation. We also find that, as the expansion time $\tau _{\textrm {exp}}$ is decreased, the characteristic magnitude at which $\delta B_{\textrm {f}}^2/B_0^2$ attains its maximum, $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}}$ , increases.

Figure 8. $(a)$ Time evolution of the square of the perturbed magnetic-field strength $\delta B_{\textrm {f}}^2/B_0^2$ associated with the firehose fluctuations for all of the D runs ( $\beta _{i0} = 50$ ). $(b)$ Maximum value of $\delta B_{\textrm {f}}^2/B_0^2$ , $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}}$ , as a function of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}(t_{\textrm {sat}})^{1.6}$ , for all runs. The dash-dotted (dotted) line shows the relationship $\delta B_{\textrm {f}}^2/B_0^2 \simeq 1.6\beta _{\|i}(t_{\textrm {sat}})/\tau _{\textrm {exp,eff}} \varOmega _i$ ( $\delta B_{\textrm {f}}^2/B_0^2 \simeq 0.77 (200)^{0.3} \beta _{\|i}(t_{\textrm {sat}})^{0.2}/(\tau _{\textrm {exp,eff}} \varOmega _i)^{0.5}$ ).

Figure 9. Two-dimensional magnetic-energy spectra $E_B(k_{\|},k_{\perp })$ of the firehose fluctuations at a selection of different times during the firehose instability’s evolution: linear phase (far left), nonlinear phase (near left), and two times during the saturated state (near and far right). The top row corresponds to an Alfvén-enabling state (run DVI), while the bottom row corresponds to an Alfvén-inhibiting state (run DIII). The region circumscribed by the dashed line indicates the region of wavenumber space ( $k_{\|} \rho _i \geqslant 0.8, \; k_{\|} \gt k_{\perp }/\tan {(15^{\circ })} \approx 3.7 k_{\perp }$ ) that is used when calculating the magnetic energy of quasi-parallel firehose modes for figure 10.

Similarly to the pressure anisotropy and effective Alfvén speed, the specific value of $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}}$ , when renormalised by $\beta _{\|i}^{0.6}$ , can be predicted with a high degree of confidence by the parameter $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}^{1.6}$ for any given $\beta _{i}$ and effective expansion time $\tau _{\textrm {exp,eff}} \varOmega _i$ . This is demonstrated in figure 8 $(b)$ . However, the exact relationship between $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}} \beta _{\|i}^{0.6}$ and $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}^{1.6}$ is not simply a power law. For values of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}^{1.6}$ of order unity, $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}} \propto \beta _{\|i}/ \tau _{\textrm {exp,eff}} \varOmega _i$ , a prediction that arises from a naive quasi-linear scattering model (see § 5.4). However, a shallower power-law dependence arises for either sufficiently small or sufficiently large values of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}^{1.6}$ . That $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}}$ is not inversely proportional to $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}$ at sufficiently small values of the latter parameter is consistent with previous shearing-box simulations of firehose-susceptible high- $\beta$ plasma (Kunz et al. Reference Kunz, Schekochihin and Stone2014a ; Melville et al. Reference Melville, Schekochihin and Kunz2016); for example, Melville et al. (Reference Melville, Schekochihin and Kunz2016) found that $(\delta B_{\textrm {f}}^2/B_0^2)_{\textrm {max}} \approx 0.77 (\beta _i/\tau _{\textrm {exp,eff}} \varOmega _i)^{0.5}$ . Computing this formula for our $\beta _{i0} = 200$ runs, we find reasonable agreement for those of our runs with the smallest values of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}^{1.6}$ (figure 8 $b$ , dotted line). That the same also holds at sufficiently large values of $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}^{1.6}$ is a new finding, suggesting that the nature of the firehose modes present in this scenario is distinct.

To explore this possibility, figure 9 displays the evolution of the magnetic-energy spectrum, $E_B(k_{\|},k_{\perp })$ , corresponding to the fluctuations visualised in figure 7, with the top (bottom) row pertaining to the Alfvén-enabling (Alfvén-inhibiting) state. As expected, the magnetic-energy spectra are initially very similar, indicating oblique modes with $k_\parallel \rho _i\approx k_\perp \rho _i\approx 0.5$ . However, in the nonlinear phases of the instability, clear differences emerge. In the saturated Alfvén-inhibiting state (bottom row), a wide range of wavenumbers is excited (including fluctuations with characteristic wavelengths that are much larger than the ion-Larmor radius), and $E_B(k_{\|},k_{\perp })$ attains a quasi-steady state. By contrast, in the saturated Alfvén-enabling state (top row), the magnetic energy is primarily concentrated in two distinct populations of fluctuations whose scales are comparable to the ion-Larmor radius: oblique firehose modes and quasi-parallel modes (the latter circumscribed in the top-right panel by the dashed line). As is also clear from figure 8( $a$ ), the ‘saturated’ Alfvén-enabling state is not quasi-steady, but instead is quasi-periodic: while the spectrum of quasi-parallel modes does not change significantly, the spectrum of oblique firehose modes evolves periodically. In § 5.2, we argue that the quasi-parallel modes are associated with a secondary parallel firehose instability.

A simple way to illustrate the quasi-periodic behaviour of firehose-instability saturation in the Alfvén-enabling state is to examine the individual components of the perturbed magnetic energy, $\delta B_{\textrm {f}}^2/B_0^2$ , i.e. the component associated with the quasi-parallel modes, $\delta B_{\textrm {f,pl}}^2/B_0^2$ , and the component associated with the oblique modes, $\delta B_{\textrm {f,ob}}^2/B_0^2$ . These components are obtained by dividing the $(k_{\|},k_{\perp })$ plane into a quasi-parallel region and a non-quasi-parallel region (see figure 9), and then calculating the total magnetic energies residing within these two separate regions. Figure 10 shows the evolution of $\delta B_{\textrm {f}}^2/B_0^2$ , $\delta B_{\textrm {f,pl}}^2/B_0^2$ and $\delta B_{\textrm {f,ob}}^2/B_0^2$ for a selection of different simulations: specifically, a simulation of an asymptotic Alfvén-enabling state (figure 10 $a$ ), a marginal Alfvén-enabling state (figure 10 $b$ ) and an Alfvén-inhibiting state (figure 10 $c$ ). In the Alfvén-enabling states, we observe that $\delta B_{\textrm {f,ob}}^2/B_0^2$ oscillates quasi-periodically, with the magnitude of that oscillation being comparable to its mean value; $\delta B_{\textrm {f,pl}}^2/B_0^2$ also oscillates with a similar period, but with a comparatively smaller amplitude relative to its mean. As the parameter $\tau _{\textrm {exp,eff}} \varOmega _i/\beta _{\|i}^{1.6}$ decreases from large to small (left to right in figure 10), both the absolute and relative amplitudes of quasi-parallel and non-quasi-parallel modes change. This can be attributed to the distinct saturation mechanisms of the quasi-parallel and oblique firehose modes (see § 5.4). In Alfvén-inhibiting states, the saturated value of $\delta B_{\textrm {f,ob}}^2/B_0^2$ does not change on a period smaller than the expansion time. Furthermore, deviations from the maximum value of $\delta B_{\textrm {f,ob}}^2/B_0^2$ are much smaller than the maximum value itself, in contrast to the Alfvén-enabling states.

Figure 10. Time evolution of the square of the perturbed magnetic-field strength $\delta B_{\textrm {f}}^2/B_0^2$ (solid black line) associated with the firehose fluctuations, along with the analogous quantity $\delta B_{\textrm {f,pl}}^2/B_0^2$ for quasi-parallel fluctuations (solid red line) and $\delta B_{\textrm {f,ob}}^2/B_0^2$ for oblique fluctuations (solid blue line), for three different simulations: $(a)$ run CV ( $\tau _{\textrm {exp}} \varOmega _{i0} = 5 \times 10^4$ , $\beta _{i0} = 25$ ), $(b)$ run CII ( $\tau _{\textrm {exp}} \varOmega _{i0} = 5 \times 10^3$ , $\beta _{i0} = 25$ ) and $(c)$ run FIII ( $\tau _{\textrm {exp}} \varOmega _{i0} = 2 \times 10^3$ , $\beta _{i0} = 200$ ).

Computing $\delta B_{\textrm {f}}^2/B_0^2$ , $\delta B_{\textrm {f,pl}}^2/B_0^2$ and $\delta B_{\textrm {f,ob}}^2/B_0^2$ for all of our simulations in the Alfvén-enabling state gives a simple way to quantify – and thereby interpret – the oscillation period of the perturbed magnetic energy. In particular, for each of these simulations, we identify a period of the simulation in which firehose instabilities have saturated, and then directly calculate the period $\tau _{\textrm {osc}}$ between the maximum value of the perturbed magnetic energy and the next minimum value. The results of this analysis are shown in figure 11. We find that $\tau _{\textrm {osc}}$ is, indeed, much smaller than $\tau _{\textrm {exp}}$ for all of our simulations that attain Alfvén-enabling states. Furthermore, $\tau _{\textrm {osc}}/\tau _{\textrm {exp}}$ is, to a reasonable degree of approximation, inversely proportional to the square root of $\tau _{\textrm {exp}} \varOmega _{i0}$ (see figure 11 a), whilst being approximately independent of $\beta _i$ (see figure 11 b). This finding is consistent with the oscillation period being comparable in magnitude to the scattering rate of particles by the quasi-parallel modes which, in the Alfvén-enabling state, have the largest amplitude of all firehose-unstable modes (see § 6.4.2). This conclusion does not seem to depend on whether $\tau _{\textrm {osc}}$ is computed from $\delta B_{\textrm {f}}^2/B_0^2$ , $\delta B_{\textrm {f,pl}}^2/B_0^2$ or $\delta B_{\textrm {f,ob}}^2/B_0^2$ (see figure 11).

Figure 11. $(a)$ Numerically determined (half-)period $\tau _{\textrm {osc}}$ of oscillation of the perturbed magnetic energy $\delta B_{\textrm {f}}^2/B_0^2$ associated with all firehose modes (black), of the magnetic energy $\delta B_{\textrm {f,pl}}^2/B_0^2$ associated with parallel modes (red) and of the magnetic energy $\delta B_{\textrm {f,ob}}^2/B_0^2$ associated with oblique modes for all of our Alfvén-enabling simulations as a function of the expansion time. The dashed grey line shows the theoretical prediction $\tau _{\textrm {osc}} \propto \tau _{\textrm {exp}}^{1/2} \varOmega _i^{-1/2}$ . $(b)$ Same as $(a)$ , but as a function of $\beta _i$ .

4.3.3. Ion distribution functions

Another, more subtle manifestation of the distinct microphysics of Alfvén-enabling and Alfvén-inhibiting states can be seen by comparing the domain-averaged ion distribution functions $f(v_{\|},v_{\perp })$ arising in the two states. The time-dependent evolution of $f(v_{\|},v_{\perp })$ in representative Alfvén-enabling and Alfvén-inhibiting states during the linear, nonlinear and saturated stages of the firehose instability is shown in figures 12 and 13, respectively. As follows directly from double-adiabatic conservation laws (4.1), the ion distribution functions in all runs initially evolve to become bi-Maxwellian, with $T_{\|i} \approx T_{\|i0}$ and $T_{\perp i} \approx T_{\perp i0}/(1+t/\tau _{\textrm {exp}})$ ; indeed, figures 12( $a$ ) and 13( $a$ ) indicate little difference between $f_{\textrm {biM}}-f_{\textrm {M}}$ (left halves of these plots) and $f-f_{\textrm {M}}$ (right halves), where $f_{\textrm {biM}}$ is the bi-Maxwellian distribution with the parallel and perpendicular temperatures computed from $f(v_{\|},v_{\perp })$ , and $f_{\textrm {M}}$ the Maxwellian distribution function with its isotropic temperature computed from $f(v_{\|},v_{\perp })$ . However, once the firehose fluctuations acquire a sufficient magnitude to backreact on the ions, the distribution functions are no longer described well as bi-Maxwellians (see figures 12 c and 13 c). In saturation (figures 12 e and 13 e), the difference becomes even more pronounced.

Figure 12. Domain-averaged ion-distribution function $f(v_{\|},v_{\perp })$ in a simulation representative of an Alfvén-enabling state (run DVI) during the ( $a$ ) linear ( $t = 0.042 \tau _{\textrm {exp}}$ ), ( $c$ ) nonlinear ( $t = 0.07 \tau _{\textrm {exp}}$ ) and ( $e$ ) saturated ( $t = 0.375 \tau _{\textrm {exp}}$ ) stages of the firehose instability. The right half of each panel shows $f-f_{\textrm {M}}$ , where $f_{\textrm {M}}$ is a Maxwellian distribution with the same temperature as $f$ ; the left half of each panel shows $f_{\textrm {biM}}-f_{\textrm {M}}$ , where $f_{\textrm {biM}}$ is a bi-Maxwellian with the same parallel and perpendicular temperatures as $f$ . ( $b$ , $d$ , $f$ ) The non-Maxwellian component of the parallel ( $f(v_{\|})-f_{\textrm {M}}(v_{\|})$ , left panel) and perpendicular ( $f(v_{\perp })-f_{\textrm {M}}(v_{\perp })$ , right panel) distribution functions at the same times, respectively. Dashed lines denote the corresponding $f_{\textrm {biM}}$ .

Figure 13. Domain-averaged ion-distribution function $f(v_{\|},v_{\perp })$ in a simulation representative of an Alfvén-inhibiting state (run DIII) during the ( $a$ ) linear ( $t = 0.075 \tau _{\textrm {exp}}$ ), ( $c$ ) nonlinear ( $t = 0.2 \tau _{\textrm {exp}}$ ) and ( $e$ ) saturated ( $t = \tau _{\textrm {exp}}$ ) stages of the firehose instability. As in figure 12, the right half of each panel shows $f-f_{\textrm {M}}$ , and the left half of each panel shows $f_{\textrm {biM}}-f_{\textrm {M}}$ . ( $b$ , $d$ , $f$ ) The non-Maxwellian component of the parallel ( $f(v_{\|})-f_{\|\textrm {M}}$ , left panel) and perpendicular ( $f(v_{\perp })-f_{\perp \textrm {M}}$ , right panel) distribution functions at the same times, respectively. Dashed lines denote the corresponding $f_{\textrm {biM}}$ .

To characterise the departures from bi-Maxwellian distribution functions more carefully – and thereby identify the subtle differences between the Alfvén-enabling and Alfvén-inhibiting states – it is helpful to define one-dimensional distribution functions: the distribution function integrated over perpendicular and parallel velocities, $f(v_\parallel ) \equiv \int _{0}^{\infty } \textrm {d} v_{\perp } v_{\perp }\, f$ and $f(v_\perp ) \equiv \int _{-\infty }^{\infty } \textrm {d} v_{\|}\, f$ , respectively. The clearest non-bi-Maxwellian feature in the nonlinear phase of both states (figures 12 $d$ and 13 $d$ ) and in saturation (figures 12 $f$ and 13 $f$ ) is the comparatively more pronounced anisotropy of the distribution function at subthermal velocities. But the main difference between the two states is the distribution function of ions with suprathermal velocities: in the Alfvén-inhibiting state (figure 13 $f$ ), the distribution function is quasi-isotropic for all velocities $|v_{\|}| \gtrsim 1.25 v_{\mathrm{th}i}$ , whereas in the Alfvén-enabling state (figure 12 $f$ ), a significant anisotropy is retained at specific velocities that evolve periodically as a function of time.Footnote ⁷ The difference is challenging to discern from the distribution functions themselves, but can be more clearly seen by comparing the pitch-angle gradient of the distribution function (see figure 14). Figure 14(a) demonstrates that in the nonlinear and saturation phases of the Alfvén-enabling state, the pitch-angle gradient of the ion distribution is not close to zero for $|v_{\|}| \gtrsim 1.75 v_{\mathrm{th}i}$ , whereas the opposite is true for the Alfvén-inhibiting state. These features of the distribution function are directly related to properties of the effective collision operator associated with the firehose fluctuations (see § 6).

Figure 14. Pitch-angle gradient of the ion distribution function $f$ divided by $2 \tilde {v} f_{\textrm {M}}$ (solid red line), where $\tilde {v} = v/v_{\mathrm{th}i}$ , and $f_{\textrm {M}}$ is a Maxwellian distribution with the same temperature as $f$ , averaged over $v_{\perp }$ . The solid blue line is the analogous quantity, but calculated using $f_{\textrm {biM}}$ , the bi-Maxwellian distribution function with the same parallel and perpendicular temperatures as $f$ . The red-pink shading denotes the standard deviation of $(\partial f_i/\partial \xi )/(2 \tilde {v} f_{\textrm {M}})$ , determined from the range of $v_{\perp }$ over which the average is computed. $(a)$ Alfvén-enabling state (run DV) in the linear phase (left panel), nonlinear phase (middle panel) and saturation (right panel). $(b)$ Alfvén-inhibiting state (run DIII).

Figure 15. $(a)$ Values of the effective collisionality $\nu _{\textrm {eff}}$ measured directly in the simulations (solid lines) with Alfvén-enabling states (runs BIV, CIV and DVI). The expansion time in these simulations is $\tau _{\textrm {exp}} \varOmega _{i0} = 2 \times 10^{4}$ . The effective collisionalities predicted by the simple model (4.13) for each simulation are shown by the dashed lines, to which the curves asymptote at late times. $(b)$ Same as $(a)$ , but for three simulations (runs DIII, EI and FIII) with $\tau _{\textrm {exp}} \varOmega _{i0} = 2 \times 10^{3}$ and therefore Alfvén-inhibiting states.

4.4. Velocity-averaged collisionality and effective viscosity

Finally, we characterise the average collisionality $\nu _{\textrm {eff}}$ of all particles in our simulation. There are various approaches for measuring $\nu _{\textrm {eff}}$ in PIC simulations; we adopt that taken in Riquelme et al. (Reference Riquelme, Quataert and Verscharen2015) and Bott et al. (Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021), and calculate $\nu _{\textrm {eff}}$ via the rate of change of the simulation-domain-averaged first adiabatic invariant $\mu$ : $\nu _{\textrm {eff}} = \dot {\bar {\mu }}/\overline {(T_{\|i}-T_{\perp i})/B}$ . We adopt this measure because, in a plasma without collisionality, $\mu$ is well conserved, so its non-conservation is a clear signature of collisionality. More practically, this measure allows for a time-resolved estimate of the effective collisionality to be computed. Figure 15 shows $\nu _{\textrm {eff}}$ as a function of time for two representative sets of simulations, each at fixed $\tau _{\textrm {exp}}$ : figure 15 $(a)$ shows three simulations in the Alfvén-enabling regime with $\tau _{\textrm {exp}} \varOmega _{i0} = 2 \times 10^{4}$ , while figure 15 $(b)$ shows three Alfvén-inhibiting simulations with $\tau _{\textrm {exp}} \varOmega _{i0} = 2 \times 10^{3}$ . Qualitatively, it is clear that $\nu _{\textrm {eff}}$ increases with increasing $\beta _{i0}$ (blue to red) and decreasing $\tau _{\textrm {exp}}$ (left to right).

Similarly to Bott et al. (Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021), we can derive a theoretical estimate for $\nu _{\textrm {eff}}$ by using the firehose-collisionality-modified CGL equations (3.2b ). To derive an estimate for $\nu _{\textrm {eff}}$ , we make three simplifying assumptions: firstly, that heat fluxes are negligible (and so all terms proportional to $q_{\|}$ or $q_{\perp }$ in (3.2b ) can be ignored); secondly, that the dimensionless pressure anisotropy is small; and thirdly, that the expansion rate is much smaller than the effective collision rate. It follows from these three assumptions that (cf. Braginskii Reference Braginskii1965)

(4.12) \begin{equation} \nu _{\textrm {eff}} \Delta _i \simeq \frac {\textrm {d} }{\textrm {d} t} \log {\frac {B}{n^{2/3}}}. \end{equation}

Finally, noting that for a transversely expanding plasma, $B \propto n$ , we deduce from (4.12) that

(4.13) \begin{equation} \nu _{\textrm {eff}}^{\textrm {CGL}} \simeq \frac {1}{3 \Delta _i} \frac {\mathrm{d}}{\mathrm{d}t} \log {B} = -\frac {1}{3 \Delta _i \tau _{\textrm {exp,eff}}}. \end{equation}

This prediction is plotted in figure 15 (dashed lines). In the Alfvén-inhibiting regime, (4.13) compares very favourably to our numerical estimates of $\nu _{\textrm {eff}}$ in the saturated states of the simulations we show (figure 15 a). In the Alfvén-enabling regime, (4.13) agrees well with the numerical collisionality averaged over the saturated state, but does not capture significant time-dependent fluctuations (figure 15 b). Because $-\Delta _i \beta _{\|i} \sim 1$ in the saturated states of our simulations, it follows that $\nu _{\textrm {eff}} \sim \beta _{\|i}/\tau _{\textrm {exp,eff}}$ , as expected.

Figure 16. $(a)$ Values of the effective collisionality $\langle \nu _{\textrm {eff}} \rangle _{\textrm {sat}}$ measured directly in all simulations, averaged over the time interval between the time at which the firehose fluctuations attain their peak strength and the time at which the next local minimum is obtained. The dashed line indicates the effective (time-averaged) value $\nu _{\textrm {eff}} = \beta _{\|i}/6 \tau _{\textrm {exp,eff}}$ of the collisionality predicted in asymptotic Alfvén-inhibiting states, while the dotted line shows the value $\nu _{\textrm {eff}} \simeq 0.21 \beta _{\|i}/ \tau _{\textrm {exp,eff}}$ appropriate for asymptotic Alfvén-enabling states. $(b)$ Effective parallel Braginskii viscosity $\mu _{\textrm {B,eff}}$ associated with the collisionality measured directly in all simulations.

Turning to our complete set of runs, figure 16 $(a)$ shows the numerical estimates of the characteristic collisionality in the saturated state of all of our simulations. To account for the time variation of the collisionality in Alfvén-enabling states, we average it over a time interval in which the saturated state is realised. The effective collisionality is consistent across all of our simulations, but $\nu _{\textrm {eff}}$ does increase slightly as the critical parameter $\tau _{\textrm {eff,exp}} \varOmega _i/\beta _{\|i}^{1.6}$ increases. This trend follows directly from our prior result that, in saturation, $\Delta _i$ increases from $\Delta _i \simeq -2/\beta _{\|i}$ to $\Delta _i \simeq -1.6/\beta _{\|i}$ as $\tau _{\textrm {eff,exp}} \varOmega _i/\beta _{\|i}^{1.6}$ increases from small (i.e. plasma in an Alfvén-inhibiting state) to large (i.e. plasma in an Alfvén-enabling state). Based on these values and (4.13), it follows that we would expect $\nu _{\textrm {eff}} \simeq \beta _{\|i}/6 \tau _{\textrm {exp,eff}}$ in Alfvén-inhibiting states (dashed line in figure 16), while $\nu _{\textrm {eff}} \simeq 0.21 \beta _{\|i}/ \tau _{\textrm {exp,eff}}$ in asymptotic Alfvén-enabling states (dotted line). The prediction is realised in our simulations.

Having computed the domain-averaged collisionality, we can then determine the plasma’s effective parallel Braginskii viscosity $\mu _{\textrm {B,eff}}$ . By comparison with (3.4), it follows that in our HEB simulations,

(4.14) \begin{equation} \mu _{\textrm {B,eff}} \approx -\frac {3}{2} (\Delta _i \beta _{\|i})_{\textrm {sat}} \frac {B^2}{4 \pi } \tau _{\textrm {exp,eff}}. \end{equation}

This estimate agrees well with the value of $\mu _{\textrm {B,eff}}$ that is directly computed from our simulations (figure 16). That $\mu _{\textrm {B,eff}}{}$ is given by (3.4) is striking for two reasons: (i) in stark contrast to classical, strongly collisional plasmas, the plasma’s viscosity is dependent upon the magnetic-field strength and (ii) the viscosity coefficient decreases as the expansion rate increases, i.e. weakly collisional plasmas behave like non-Newtonian fluids.

5.1. Overview

While some of the results from our HEB simulations – for example, the particle-averaged effective collisionality, or the regulation of pressure anisotropy in Alfvén-inhibiting states – are consistent with the results of previous simulations of firehose instabilities (e.g. Hellinger & Trávníček Reference Hellinger and Trávníček2008; Kunz et al. Reference Kunz, Schekochihin and Stone2014a ; Riquelme et al. Reference Riquelme, Quataert and Verscharen2015; Melville et al. Reference Melville, Schekochihin and Kunz2016), other results are not, and so require further analysis and interpretation. Three findings about the Alfvén-enabling state in particular are unexpected, and so warrant additional investigation. First of these is the emergence of ion-Larmor-scale parallel firehose modes, which are specifically predicted not to be present by the linear theory of the firehose instability in a bi-Maxwellian plasma that was outlined in § 2.5. Secondly, the regulated pressure anisotropy $(\Delta _i)_{\textrm {sat}} \simeq -1.6/\beta _{\|i}$ in the Alfvén-enabling state does not correspond to the linear threshold $\Delta _{\textrm {cr}} \simeq -1.35/\beta _{\|i}$ for the oblique firehose instability in a bi-Maxwellian plasma (see § 2.3). The third unexpected finding is that the box-averaged perturbed magnetic energy, $\delta B_{\textrm {f}}^{2}/B_0^2$ , does not scale as $\delta B_{\textrm {f}}^{2}/B_0^2 \propto \beta _i/\tau _{\textrm {exp,eff}} \varOmega _i$ as might be naively anticipated, but instead has a weaker dependence. These findings are discussed in §§ 5.2–5.4, respectively.

5.2. Secondary parallel firehose instability in the Alfvén-enabling regime

A notable result from our simulations is the presence of ion-Larmor-scale parallel firehose modes in the Alfvén-enabling regime. The presence of such modes is, at first glance, inconsistent with the linear theory of the firehose instability in a bi-Maxwellian plasma with a negative pressure anisotropy (§ 2.5), which predicts that the resonant parallel firehose should be subdominant to oblique firehose modes in high- $\beta$ plasma. However, it can, in fact, be shown that these modes are not whistler/fast magnetosonic modes destabilised by the resonant parallel firehose instability (as would occur in the plasma with $\beta _i \sim 1$ ), but are instead a lower-frequency mode excited by a (newly identified) secondary instability associated with the non-bi-Maxwellian form of the distribution function. This form, presented in § 4.3.3, is caused by the backreaction of the oblique firehose modes on the otherwise bi-Maxwellian distribution function that is driven by the plasma’s expansion.

To understand this secondary parallel firehose instability better, it is helpful to describe qualitatively the types of parallel modes (and their growth) that the (high- $\beta$ ) plasma can support linearly as the expansion proceeds. Initially, at the start of the simulation, when the ion distribution function is Maxwellian, there are two types of forward-propagating parallel modes with $\beta _i^{-1/2} \ll k_{\|} \rho _i \lesssim 1$ : right-handed whistler/magnetosonic modes, which have characteristic real frequencies $\varpi \sim k_{\|}^2 \rho _i^2 \varOmega _i$ , and left-handed, ion-cyclotron modes which, in high- $\beta _i$ plasma, have $\varpi \sim \varOmega _i/\beta _i$ (Foote & Kulsrud Reference Foote and Kulsrud1979).Footnote ⁸ These two types of modes have very different characteristic frequencies because of their distinct physical mechanisms; the characteristic oscillation of the higher-frequency whistler/magnetosonic modes is supported by inertial and gyroviscous forces acting out of phase (with the action of the Alfvén restoring force being negligible), while for the lower-frequency ion-cyclotron modes, it is the out-of-phase action of the Alfvén restoring force and the gyroviscous force that gives rise to oscillatory dynamics. Despite their distinct mechanisms, both of these modes are damped ( $\gamma \lt 0$ ). As the plasma expands, the pressure anisotropy becomes increasingly negative, which changes the character of the ion-cyclotron mode; specifically, the real frequency of this mode becomes negative for $k_{\|} \rho _i \sim 1$ (though the mode remains damped). Because $k_{\|} \gt 0$ , this change of sign corresponds to initially forward-propagating ion-cyclotron modes becoming backward-propagating (and vice versa); in short, the initially left-handed forward-propagating ion-cyclotron mode becomes a type of right-handed (forward-propagating) mode that is qualitatively distinct from the whistler/magnetosonic mode. Physically, this change of handedness can be attributed to the Alfvénic restoring force being weakened by increasingly strong parallel pressure forces associated with the negative pressure anisotropy. The damping of these ion-cyclotron modes finally becomes growth once the oblique firehose fluctuations begin to backreact on the ion distribution function. These fluctuations, which have a characteristic parallel wavenumber $k_{\|} \approx 0.5 \rho _i^{-1}$ , efficiently scatter particles with a characteristic velocity $v_{\|} \approx v_{\mathrm{th}i}/(k_{\|} \rho _i)_{\textrm {ob}} \approx 2 v_{\mathrm{th}i}$ , and isotropise the distribution in a narrow $v_{\|}$ interval. This, in turn, enables the right-handed ion-cyclotron modes to extract energy from these same particles, and thereby grow.

With some effort, we can characterise the growth of the secondary parallel firehose modes (and their analogous damped modes in the initial stage of the simulation) analytically. For arbitrary background distribution functions ${f}_{s0}$ of species $s$ , the linear dispersion relation of parallel modes in a hot plasma is, neglecting the displacement current,

(5.1) \begin{align} D^{\pm } & = k_{\|}^2 c^2 + \sum _{s} \frac {\omega _{\mathrm{p}s}^2}{n_{s0}} \Bigg \{\pi \int _{C_{\textrm {L}}} \mathrm{d}v_{\|} \int _0^{\infty } \mathrm{d}v_{\perp } \frac {v_{\perp }^2}{k_{\|}v_{\|}-\omega \mp \varOmega _s} \nonumber \\ &\quad \times \left [k_{\|} \left (v_{\perp } \frac {\partial {f}_{s0}}{\partial v_{\|}}-v_{\|} \frac {\partial {f}_{s0}}{\partial v_{\perp }}\right )+\omega \frac {\partial {f}_{s0}}{\partial v_{\perp }}\right ]\Bigg \} = 0 , \end{align}

where $n_{s0}$ is the equilibrium number density of species $s$ , $\omega _{\mathrm{p}s}$ is the plasma frequency, $C_{\textrm {L}}$ is the usual Landau (‘L’) contour and we have assumed that $k_{\|} \gt 0$ . In a Maxwellian plasma, the ‘ $+$ ’ and ‘ $-$ ’ roots with $\omega \gt 0$ correspond to the whistler/magnetosonic modes and ion-cyclotron modes, respectively. Motivated by our simulation results, we further specialise to the ‘low-frequency’ ion-cyclotron modes with $k_{\|} \rho _i \sim 1$ , which satisfy $\omega \sim k_{\|} v_{\mathrm{th}i}/\beta _{i} \ll k_{\|} v_{\mathrm{th}i}$ . We also assume a Maxwellian electron population (as in the hybrid-kinetic simulations), and that the ion distribution function’s anisotropy is small compared with its characteristic magnitude:

(5.2) \begin{equation} \frac {v_{\mathrm{th}i}}{v}\frac {\partial {f}_{i0}}{\partial \xi } = v_{\mathrm{th}i} \left (\frac {\partial {f}_{i0}}{\partial v_{\|}}-\frac {v_{\|}}{v_{\perp }} \frac {\partial {f}_{i0}}{\partial v_{\perp }} \right ) \sim \frac {\omega }{k_{\|} v_{\mathrm{th}i} } f_{\mathrm{M}} \ll f_{\mathrm{M}} , \end{equation}

where $\xi \equiv v_{\|}/v$ is the pitch angle. Under these assumptions, simplified expressions can be derived for the real frequency $\varpi$ and growth rate $\gamma$ of these modes:

(5.3a) \begin{align} \frac {\varpi }{\varOmega _i} & \approx \pm \Bigg \{ \mathcal{G}(k_{\|} \rho _i)\left [k_{\|}^2 d_i^{2} + \left (\frac {\pi }{n_{i0}} \mathcal{P}\int _{-\infty }^{\infty } \mathrm{d} {v}_{\|} \int _0^{\infty } \mathrm{d}{v}_{\perp } \frac {k_{\|} {v}_{\perp }^3}{k_{\|}v_{\|} \mp \varOmega _i} \frac {1}{v} \frac {\partial {f}_{i0}}{\partial \xi }\right ) \right ] \nonumber\\& \quad - \left [\frac {\pi ^2}{n_{i0}} \int _0^{\infty } \mathrm{d}v_{\perp } v_{\perp }^{3} \left (\frac {1}{v} \frac {\partial {f}_{i0}}{\partial \xi }\right )\Bigg |_{v_{\|} = v_{\|\mathrm{res}}^{\pm }}\right ] \left [\frac {\sqrt {\pi }}{k_{\|} \rho _i} \exp {\left (-\frac {1}{k_{\|}^2 \rho _i^2}\right )} \right ] \Bigg \} \nonumber \\& \quad \times \left \{\left [\mathcal{G}(k_{\|} \rho _i)\right ]^{2} + \frac {\pi }{k_{\|}^2 \rho _i^2} \exp {\left (-\frac {2}{k_{\|}^2 \rho _i^2}\right )} \right \}^{-1}, \end{align}

(5.3b) \begin{align} \frac {\gamma }{\varOmega _i} & \approx \pm \left [\mathcal{G}(k_{\|} \rho _i)\right ]^{-1} \left [\frac {\pi ^2}{n_{i0}} \int _0^{\infty } \mathrm{d}v_{\perp } v_{\perp }^{3} \left (\frac {1}{v} \frac {\partial {f}_{i0}}{\partial \xi }-\frac {2 v_{\textrm {wv}}}{v_{\mathrm{th}i}^2} {f}_{\textrm {M}}\right )\Bigg |_{v_{\|} = v_{\|\mathrm{res}}^{\pm }}\right ] , \end{align}

where $d_i = \beta _i^{-1/2} \rho _i$ is the ion inertial length,

(5.4) \begin{equation} \mathcal{G}(k_{\|} \rho _i) = 1 + \frac {1}{k_{\|} \rho _i} \mathrm{Re} \, Z\left (\frac {1}{k_{\|} \rho _i}\right ) \end{equation}

is a special function related to the plasma dispersion function $Z(x)$ whose only root occurs at $k_{\|} \rho _i \simeq 1.08$ , $v_{\mathrm{wv}} \equiv \varpi /k_{\|}$ is the parallel phase velocity of the wave and $v_{\|\mathrm{res}}^{\pm } \equiv (\varpi \pm \varOmega _i)/k_{\|} \approx \pm \varOmega _i/k_{\|}$ is the parallel velocity of particles that are resonant with that mode. We note that, due to our assumed ordering, we have removed the whistler/magnetosonic branch, and so (5.3b ) describes just the real frequency and growth rate of (both forward- and backward-propagating) modes of the ion-cyclotron type. It follows that the damping or growth of such parallel modes depends upon the sign of the quantity

(5.5) \begin{equation} \mathcal{I}_{\pm } \equiv \frac {1}{v} \frac {\partial {f}_{i0}}{\partial \xi }-\frac {2 v_{\textrm {wv}}}{v_{\mathrm{th}i}^2} {f}_{\textrm {M}} \end{equation}

evaluated near the resonant velocity $v_{\|\mathrm{res}}^{\pm }$ . For $k_{\|} \rho _i \lt 1.08$ , growth occurs whenever $\pm \mathcal{I}_{\pm } \gt 0$ , and vice versa for $k_{\|} \rho _i \gt 1.08$ .Footnote ⁹

We can use (5.3a ) to evaluate $\varpi$ and $\gamma$ as the ion distribution function evolves from a Maxwellian via a bi-Maxwellian distribution to the non-bi-Maxwellian state associated with scattering by oblique firehoses. In a plasma with a bi-Maxwellian ion distribution, (5.3b ) simplifies considerably, because we have

(5.6) \begin{equation} \frac {v_{\mathrm{th}i}}{v}\frac {\partial {f}_{i0}}{\partial \xi } = -2 \Delta _i \frac {v_{\|}}{v_{\mathrm{th}i}} {f}_{\textrm {M}} , \quad \mathcal{I}_{\pm }|_{v_{\|} = v_{\|\mathrm{res}}^{\pm }} = \frac {2}{v_{\mathrm{th}i}}\left (\mp \frac {\Delta _i}{k_{\|} \rho _i} - \frac {v_{\textrm {wv}}}{v_{\mathrm{th}i}}\right ) {f}_{\textrm {M}} , \end{equation}

so that

(5.7a) \begin{align} \varpi & \approx \pm \Delta _i \pm \varOmega _i \frac {k_{\|}^2 \rho _i^{2}}{\beta _i} \mathcal{G}(k_{\|} \rho _i) \left \{\left [\mathcal{G}(k_{\|} \rho _i)\right ]^{2} + \frac {\pi }{k_{\|}^2 \rho _i^2} \exp {\left (-\frac {2}{k_{\|}^2 \rho _i^2}\right )}\right \}^{-1} ,\\[-10pt]\nonumber \end{align}

(5.7b) \begin{align} \gamma & \approx - \varOmega _i \frac {\sqrt {\pi } k_{\|} \rho _i}{\beta _i} \exp {\left (-\frac {1}{k_{\|}^2 \rho _i^2}\right )} \left \{\left [\mathcal{G}(k_{\|} \rho _i)\right ]^{2} + \frac {\pi }{k_{\|}^2 \rho _i^2} \exp {\left (-\frac {2}{k_{\|}^2 \rho _i^2}\right )} \right \}^{-1} . \quad \end{align}

In the plasma’s initial state, in which $\Delta _i = 0$ , the forward-propagating modes are indeed those associated with the ‘ $-$ ’ root, as expected, and so are left-handed, because the numerator of (5.7a ) is negative for $k_{\|} \rho _i \lt 1.08$ . These equations further imply that $\gamma \lt 0$ initially. In the bi-Maxwellian stage, (5.7a ) indicates that, for $\Delta _i \lt 0$ , the ‘ $-$ ’ mode transitions from being forward-propagating to backward-propagating at a smaller value of $k_{\|} \rho _i$ than for a Maxwellian distribution (and vice versa for the ‘ $+$ ’ mode). When $\Delta _i \lt -2/\beta _i$ , $\varpi \lt 0$ at all wavenumbers $k_{\|} \sim \rho _i^{-1}$ for the ‘ $-$ ’ mode, and $\varpi \gt 0$ for the ‘ $+$ ’ mode. However, both the ‘ $+$ ’ and ‘ $-$ ’ mode are still damped at this stage by ions with $v_{\|} \approx \pm \varOmega _i/k_{\|}$ . Finally, in the state with the non-bi-Maxwellian distribution, scattering by the oblique firehoses causes ${\partial {f}_{i0}}/{\partial \xi }|_{v_{\|} = v_{\|\mathrm{res}}}$ to decrease in magnitude near $v_{\|} \approx \pm 2 v_{\mathrm{th}i}$ , and $\varpi$ does not change its sign when these resonant particles start to be isotropised, because (5.3a ) implies that $\varpi$ is less sensitive than $\gamma$ to the value of $f_i$ at specific $v_{\|}$ . Once $\mathcal{I}_{\pm }$ reverse their sign for forward- and backward-propagating resonant parallel modes, respectively, it then follows that their growth rate becomes positive.

Figure 17. $(a)$ Two-dimensional magnetic-energy spectra of the firehose fluctuations in run CV at a selection of different times around the emergence of the parallel secondary firehose instability. $(b)$ Pitch-angle gradient of the ion distribution function $f$ divided by $2 \tilde {v} f_{\textrm {M}}$ (solid red line), where $\tilde {v} = v/v_{\mathrm{th}i}$ , and $f_{\textrm {M}}$ is a Maxwellian distribution with the same temperature as $f$ , averaged over $v_{\perp }$ , at the same times shown in $(a)$ . The solid blue line is the analogous quantity, but calculated using $f_{\textrm {biM}}$ , the bi-Maxwellian distribution function with the same parallel and perpendicular temperatures as $f$ . The dotted red and blue lines show $\tilde {v}_{\textrm {wv}} = {v}_{\textrm {wv}}/v_{\mathrm{th}i}$ calculated using a linear dispersion relation solver that finds the complex frequency of low-frequency modes with a given input numerical distribution function. $(c)$ Approximate real frequencies (red) and damping rates (blue) (cf. (5.3b )) of the ‘ $-$ ’ root, for $f$ (solid lines), $f_{\textrm {M}}$ (dotted lines) and $f_{\textrm {biM}}$ (dot-dashed lines) at the same times indicated in $(a)$ .

This evolution is illustrated using one of our simulations (run CV, an ‘asymptotic’ Alfvén-enabling simulation) in figure 17. Figure 17 $(a)$ shows the two-dimensional magnetic-energy spectrum at various times in the simulation around the time at which the parallel modes are observed; figure 17 $(b)$ shows the pitch-angle gradient of the ion distribution function $f_{i}$ at those same times; and figure 17 $(c)$ shows $\varpi$ and $\gamma$ of the ‘ $-$ ’ modes, which we calculate using the approximate expressions (5.3b ). The integrals in these expressions for $\varpi$ and $\gamma$ are evaluated numerically, taking as their input the numerical distribution function. We see that, in the initial stages of the growth of oblique firehose modes (figure 17 $a$ , left panel), when $f_{i}$ is still approximately bi-Maxwellian (and so the pitch-angle gradient of $f_i$ is well described by (5.6) – see figure 17 $b$ , left panel), parallel ‘ $-$ ’ modes with $k_{\|} \rho _i \gt 0.5$ have a negative sign, but are damped (figure 17 $c$ , left panel). However, concurrently with the emergence of parallel modes (figure 17 $a$ , middle panel), the ion distribution function becomes non-bi-Maxwellian (figure 17 $b$ , middle panel), and parallel ion-cyclotron modes become linearly unstable (figure 17 $c$ , middle panel), albeit over quite a narrow range of wavenumbers. For the fastest growing modes,

(5.8) \begin{equation} v_{\mathrm{th}i} \mathcal{I}_{-} = \frac {v_{\mathrm{th}i}}{v} \frac {\partial {f}_i}{\partial \xi }-\frac {2}{k_{\|} v_{\mathrm{th}i}} \varpi \left (-\frac {1}{k_{\|}\rho _i}\right ) {f}_{\textrm {M}} \ll \left |\frac {\varpi }{k_{\|} v_{\mathrm{th}i}}\right | {f}_{\textrm {M}} \sim \frac {1}{\beta _i} {f}_{\textrm {M}} , \end{equation}

and so, in contrast to the bi-Maxwellian, $f_i$ has the property that its pitch-angle gradient is approximately equal to twice the (normalised) phase velocity $v_{\textrm {wv}} \sim v_{\mathrm{th}i}/\beta _{i}$ of the linear modes that $f_i$ supports. In other words, the ion distribution function’s anisotropy is constrained by the parallel modes’ phase velocity. As the simulation progresses further, the unstable parallel modes tend to acquire slightly larger wavenumbers (figure 17 $a$ , right panel), with those that were initially unstable becoming forward-propagating, stable modes again (figure 17 $c$ , right panel).

In summary, scattering by the ion-Larmor-scale oblique firehose modes that initially arise due to the negative pressure anisotropy is responsible for the development of a non-bi-Maxwellian distribution function, which in turn is subject to an instability of right-handed parallel modes that would not be present if the distribution function were to have remained bi-Maxwellian. This secondary firehose instability could also explain the persistent parallel modes with $k_{\|} \rho _i \sim 1$ seen in regions of negative pressure anisotropy within the hybrid-kinetic simulation of a long-wavelength, large-amplitude Alfvén wave reported by Squire et al. (Reference Squire, Kunz, Quataert and Schekochihin2017).

5.3. Why $(\Delta _i)_{\textrm {sat}} \simeq -1.6/\beta _{\|i}$ in high- $\beta _i$ Alfvén-enabling states

One finding of our simulation results that was not anticipated from the linear theory of the firehose instability outlined in § 2 is that, in saturation, $\Delta _i \simeq -1.6/\beta _{\|i}$ . If, as we argue in § 5.4, the saturation of the oblique firehose instability can be described by quasi-linear theory, then it must be the case that the plasma attains a saturated state that is close to marginality with respect to the oblique firehose instability. However, we showed in § 2.3 that, in a bi-Maxwellian plasma, the oblique firehose instability’s threshold is given by $\Delta _{\textrm {cr}} \simeq -1.35/\beta _{\|i} \gt (\Delta _i)_{\textrm {sat}}$ . Naively, it might therefore be expected that oblique firehose modes should still grow, and, adopting the estimate (2.4) for these modes’ growth rate, will do so at a rate that far exceeds the rate of the plasma’s expansion, $\gamma _{\perp \mathrm{f}} \approx 0.7 \varOmega _i/\beta _{i} \gtrsim 19 \beta _{i}^{0.6}/\tau _{\textrm {exp,eff}}$ .

This seeming contradiction is resolved by the fact that the plasma’s ion distribution $f_i$ is not well modelled as a bi-Maxwellian distribution, but instead has a distinct form of anisotropy. More specifically, as was illustrated in § 4.3.3, the anisotropy of $f_i$ is concentrated at smaller characteristic values of $v$ compared with those of a bi-Maxwellian. This has the consequence of bringing the threshold of ion-Larmor-scale oblique firehose modes closer to the fluid firehose threshold. That the modified form of the anisotropy alters the oblique firehose instability’s threshold can be demonstrated mathematically by considering the leading-order FLR corrections to the fluid firehose threshold, which are computed for a general ion distribution function in Appendix A.3 (cf. (A.59)):

(5.9) \begin{equation} \frac {2}{\beta _{\|i}} + \Delta _i +k_{\|}^2 {\rho }_i^2 \mathcal{A}_{4i} - \frac {3}{16} k_{\perp }^2 {\rho }_i^2 \mathcal{B}_{4i} = \mathcal{O}\big(\Delta _i k^4 {\rho }_i^4\big), \end{equation}

where $\mathcal{A}_{4i}$ and $\mathcal{B}_{4i}$ are given by (cf. (A.55e ))

(5.10a) \begin{align} \mathcal{A}_{4i} & = -8 \pi \frac {T_{\|i}}{T_{i}} \frac {1}{n_{i} v_{\mathrm{th}i}^4} \int _{-\infty }^{\infty } \mathrm{d} {v}_{\|} \int _0^{\infty } \mathrm{d} {v}_{\bot } \, {v}_{\bot } {v}_{\|}^2 \left ({v}_{\|}^2-\frac {1}{4} {v}_{\perp }^2 \right ) {f}_{i} ,\\[-10pt]\nonumber \end{align}

(5.10b) \begin{align} \mathcal{B}_{4i} & = -4 \pi \frac {T_{\|i}}{T_{i}} \frac {1}{n_{i} v_{\mathrm{th}i}^4} \int _{-\infty }^{\infty } \mathrm{d} {v}_{\|} \int _0^{\infty } \mathrm{d} {v}_{\bot } \, {v}_{\bot }^3 \left ({v}_{\|}^2-\frac {1}{2} {v}_{\perp }^2 \right ) {f}_{i} . \end{align}

Inspecting the velocity-space integrands in (5.10b ) and comparing them with the analogous integrand for the pressure anisotropy,

(5.11) \begin{equation} \Delta _i = - 4 \pi \frac {T_{\|i}}{T_{i}} \frac {1}{n_{i} v_{\mathrm{th}i}^2} \int _{-\infty }^{\infty } \mathrm{d} {v}_{\|} \int _0^{\infty } \mathrm{d} {v}_{\bot } \, {v}_{\bot } \left ({v}_{\|}^2-\frac {1}{2} {v}_{\perp }^2 \right ) {f}_{i} , \end{equation}

it is clear that concentrating the anisotropy of a distribution function at smaller characteristic velocities will in general reduce the values of the ratios $\mathcal{A}_{4i}/\Delta _i$ and $\mathcal{B}_{4i}/\Delta _i$ . Thus, the distribution functions attained in the saturated state of the firehose instability simultaneously maintain comparatively larger values of $(\Delta _i)_{\textrm {sat}}$ than a bi-Maxwellian distribution and smaller values of $\mathcal{A}_{4i}$ and $\mathcal{B}_{4i}$ . Computing $\mathcal{A}_{4i}$ and $\mathcal{B}_{4i}$ directly for our ‘asymptotic’ Alfvén-enabling regime simulation (run CV), we find $\mathcal{A}_{4i} \simeq -1.6/\beta _{\|i}$ and $\mathcal{B}_{4i} \simeq -0.8/\beta _{\|i}$ ; setting $k_{\|} \rho _i \approx k_{\perp } \rho _i \simeq 0.5$ to match those of the dominant oblique firehose mode, (5.9) predicts that $(\Delta _i)_{\textrm {sat}} \approx -1.6/\beta _{\|i}$ . This agrees very well with its actual value in the simulation.

An outstanding question that follows naturally from our result is why, in prior $\beta _{\|i} \gtrsim 1$ simulations of firehose-susceptible plasmas (see e.g. Hellinger & Trávníček 2008; Hellinger et al. Reference Hellinger, Matteini, Landi, Franci, Verdini and Papini2019; Bott et al. Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021), it was found that $(\Delta _i)_{\textrm {sat}} \simeq -1.4/\beta _{\|i}$ , in closer agreement with the bi-Maxwellian threshold of the oblique firehose instability. The most plausible explanation of this (small) discrepancy pertains to the different linear characteristics of firehose instabilities at $\beta _{\|i} \gg 1$ versus $\beta _{\|i} \gtrsim 1$ . Specifically, as we demonstrated in § 2.5, when $\beta _{\|i} \gtrsim 1$ , the growth rate of resonant parallel firehose modes tends to be comparable to those of oblique modes. The presence of a saturated population of such modes, which would not be present in high- $\beta _i$ plasma, would be expected to affect the specific value of $\Delta _i$ attained in the saturated state. We note that, though the specific values of $(\Delta _i)_{\textrm {sat}}$ are distinct, both are such that the plasma still attains an Alfvén-enabling state.

5.4. The perturbed magnetic energy of firehose fluctuations in saturation: part I

It was shown in § 4.3.2 that the relationship between the perturbed magnetic energy associated with the firehose fluctuations in saturation and macroscopic plasma parameters is not simply a power law across all values of the key parameter $\tau _{\textrm {exp}} \varOmega _i/\beta _{\|i}^{1.6}$ , with a change occurring near the transition between the Alfvén-enabling and Alfvén-inhibiting states. This implies that the saturation physics in these two states must be distinct.

Such a conclusion is, at first glance, counter-intuitive. For the ion-Larmor-scale firehose modes that we observe in both the Alfvén-enabling and Alfvén-inhibiting states, which to a good approximation consist of perturbations to the direction of the magnetic field, a saturated state is most plausibly maintained via pitch-angle scattering at a rate sufficient to maintain near-marginality with respect to the firehose instability’s threshold. Assuming that the rate $\nu _{\textrm {eff}} \sim \beta _i/\tau _{\textrm {exp}}$ of pitch-angle scattering by (ion-Larmor-scale) fluctuations is related to their amplitude $\delta B_{\textrm {f}}/B_0$ by $\nu _{\textrm {eff}} \sim \varOmega _i \delta B_{\textrm {f}}^2/B_0^2$ – in effect, adopting a quasi-linear scattering model based on the assumption that $\delta B_{\textrm {f}} \ll B_0$ – we deduce that

(5.12) \begin{equation} \frac {\delta B_{\textrm {f}}^2}{B_0^2} \sim \frac {\beta _i}{\tau _{\textrm {exp}} \varOmega _i} . \end{equation}

We note that such a quasi-linear model should be self-consistent for any firehose-susceptible plasma in an Alfvén-enabling state, because ${\delta B_{\textrm {f}}^2}/{B_0^2} \sim {\beta _i}/{\tau _{\textrm {exp}} \varOmega _i} \lesssim \beta _i^{-0.6} \ll 1$ .

This argument, which provides testable predictions for the dependence of ${\delta B_{\textrm {f}}^2}/{B_0^2}$ on $\beta _i$ , $\tau _{\textrm {exp}}$ and $\varOmega _i$ , only partially accounts for the results of our numerical study. The scaling $\nu _{\textrm {eff}} \sim \beta _i/\tau _{\textrm {exp}}$ for the effective collisionality is indeed the same as reported in § 4.4. However, the scaling (5.12) only agrees for our simulations in Alfvén-inhibiting states, not Alfvén-enabling ones. We conclude that the argument must overlook aspects of firehose-instability saturation that affect the scaling of the perturbed magnetic energy.

In order to resolve this discrepancy, a more nuanced understanding of scattering of particles by both oblique firehose and secondary parallel firehose modes in Alfvén-enabling states – and how this leads to the saturation of both types of firehose instability – is required. We therefore characterise an effective ‘firehose collision operator’ in the next section.

6.1. Overview

One key property of the firehose instability in its saturated state is that it provides the plasma with an effective collisionality, $\nu _{\textrm {eff}}$ . Particles in the plasma experience this collisionality predominantly as pitch-angle scattering. In this section, we move beyond previous velocity-averaged estimates of this collisionality, and propose a model in the Alfvén-enabling state for the velocity-dependent pitch-angle scattering rate of particles with speeds of the order of the thermal speed. This allows us to construct a simple ‘effective firehose collision operator’, given by

(6.1) \begin{equation} \mathfrak{C}_{\textrm {f}}[f] = \frac {1}{2} \frac {\partial }{\partial \xi }\Bigg \{(1-\xi ^2)\nu _{\textrm {eff,pl}}(v \xi )\left [\frac {\partial f}{\partial \xi } - 2 \tilde {w} \tilde {v}_{\textrm {wv,pl}}(v \xi ) f_{\textrm {M}}\right ] + (1-\xi ^2)\nu _{\textrm {eff,ob}}(v \xi ) \frac {\partial f}{\partial \xi } \Bigg \} , \end{equation}

where

(6.2a) \begin{align} \nu _{\textrm {eff,pl}}(v_{\|}) & = \frac {0.15 v_{\mathrm{th}i}}{|{v}_{\|}|} \frac {\beta _i^{1/4} \varOmega _i^{3/4}}{\tau _{\textrm {exp}}^{1/4}} \exp {\left [-0.31(\tau _{\textrm {exp}}\varOmega _i)^{1/2} \left (\frac {v_{\mathrm{th}i}}{|{v}_{\|}|}-1.2\right )^2\right ]} \nonumber \\ & \quad + 0.09 H(k_{\|} \rho _i - 1.2) ({\tau \varOmega _i})^{-1/2} (k_{\|} \rho _i)^{-2.7} ,\\[-10pt]\nonumber \end{align}

(6.2b) \begin{align} v_{\textrm {wv,pl}}(v_{\|}) & = \mathrm{sgn}(v_{\|})\frac {v_{\mathrm{th}i}}{\beta _i}\left (4.9-2.9 \frac {|v_{\|}|}{v_{\mathrm{th}i}} \right ) ,\\[-10pt]\nonumber \end{align}

(6.2c) \begin{align} \nu _{\textrm {eff,ob}}(v_{\|}) & = \frac {1.4 v_{\mathrm{th}i}}{|{v}_{\|}|} \frac {\beta _i}{\tau _{\textrm {exp}}} \exp {\left [-13 \left (\frac {v_{\mathrm{th}i}}{|{v}_{\|}|}-0.75\right )^2\right ]} ,\\[-10pt]\nonumber \end{align}

where $H(x)$ denotes the Heaviside step function. We then compare the predicted properties of this collision operator with two different numerical diagnostics applied to our simulations, and confirm that the model collision operator accounts for both the characteristic anisotropy of the ion distribution function and the root mean square of the firehose fluctuations’ magnetic-field strength. In turn, this collision operator allows us to advance our qualitative understanding of the anomalous scaling of the perturbed magnetic energy in Alfvén-enabling states discussed in § 5.4.

6.2. An effective firehose collision operator

Beyond accounting for the saturated amplitude of firehose-unstable modes, there are two other motivations for investigating the velocity-dependent collisionality associated with firehose fluctuations. Firstly, it is the velocity dependence of effective collisions that determines the ion distribution function’s anisotropy, and thereby the specific saturation value of the pressure anisotropy at which further growth of firehose-unstable modes is inhibited. As discussed in § 2, long-wavelength firehose modes are insensitive to the form of the ion distribution function’s anisotropy, but kinetic-scale firehose modes are sensitive to it. Because it is these kinetic-scale modes that have the least stringent threshold for instability, the specific form of anisotropy is pertinent. Secondly, for certain other problems in astrophysical plasmas such as modelling cosmic-ray transport, understanding the effective collisionality of particles with specific velocities due to firehose fluctuations is a crucial component of the problem’s solution.

In general, characterising the effective collision operator associated with arbitrary electromagnetic fluctuations, which could cause slowing, parallel diffusion, and/or perpendicular diffusion of particles, is quite challenging. However, in the specific case of the effective collision operator associated with firehose fluctuations, various simplifying assumptions can be reasonably adopted. Based on the small amplitude of firehose fluctuations realised in the Alfvén-enabling state (figure 8 $b$ implies that the total magnetic energy of fluctuations satisfies $\delta B_{\textrm {f}}^2/B_0^2 \ll \beta _{\|i}^{-0.6} \ll 1$ ) and their broad spectra (see figure 9), we assume that the collision operator can be described by quasi-linear theory. Furthermore, we neglect the electric fields associated with the firehose fluctuations on the grounds that the electric contribution to the total Lorentz force is subdominant to the magnetic force; it follows from Faraday’s law that, for firehose fluctuations, $c \delta {\boldsymbol{E}}/|{\boldsymbol{v}} \,{\boldsymbol{\times }}\, \delta {\boldsymbol{B}}| \sim \omega /k v_{\mathrm{th}i} \sim 1/\beta _i \ll 1$ . Finally, we assume (based on our simulation results) that the magnetic-field perturbations caused by the firehose instability satisfy $\delta \boldsymbol{B} \approx \delta \boldsymbol{B}_{\perp }$ . Taking these assumptions together, the quasi-linear collision operator arising from magnetic fluctuations is simply a resonant pitch-angle-scattering operator that isotropises the distribution function in the frame moving at the (parallel) phase velocity $v_{\textrm {wv}} = v_{\mathrm{th}i} \tilde {v}_{\textrm {wv}}$ of the firehose modes (the wave frame) at a velocity-dependent scattering rate $\nu _{\textrm {eff}}(v_{\|}',v_{\perp })$ given by (e.g. Kulsrud & Pearce Reference Kulsrud and Pearce1969)

(6.3) \begin{equation} \nu _{\textrm {eff}}(v_{\|}',v_{\perp }) = \pi \frac {\varOmega _i^2}{v_{\|}'}\frac { \tilde {\mathcal{E}}_{B}\bigl ({\varOmega _i}/{v_{\|}'}\bigr )}{B_0^2/8\pi } , \quad \tilde {\mathcal{E}}_{B}(k_{\|}) \equiv \sum _{n \neq 0} n^2 \int \mathrm{d}^2 \boldsymbol{k}_{\perp } \, E_{B}(n k_{\|},k_{\perp }) \frac {[\textrm {J}_n(k_{\perp } v_{\perp }/\varOmega _i)]^2}{k_{\perp }^2 v_{\perp }^2/\varOmega _i^2} , \end{equation}

where the primes denote parallel velocities evaluated in the wave frame and $\textrm {J}_{n}(x)$ is the $n$ th-order Bessel function of the first kind.

If we also assume that both the anisotropy of the distribution function and $v_{\textrm {wv}}$ are small – more precisely, that $({\partial f_i}/{\partial \xi })/f_{\textrm {M}} \sim \tilde {v}_{\textrm {wv}} \sim 1/\beta _i \ll 1$ – it can be shown (see e.g. Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025) that the quasi-linear pitch-angle operator in the plasma’s rest frame has the following form:

(6.4) \begin{equation} \mathfrak{C}[f] = \frac {1}{2} \frac {\partial }{\partial \xi }\left \{(1-\xi ^2)\nu _{\textrm {eff}}(v,\xi )\left [\frac {\partial f}{\partial \xi } - 2 \tilde {w} \tilde {v}_{\textrm {wv}}(v,\xi ) f_{\textrm {M}}\right ] \right \} , \end{equation}

where we remind the reader that $\xi = v_{\|}/v$ is the pitch angle, $v \equiv \sqrt {v_{\|}^2 + v_{\perp }^2}$ is the particle speed and $\tilde {v}_{\textrm {wv}}(v,\xi )$ is the parallel phase velocity of the firehose modes with which specific particles having peculiar velocity ( $v,\xi$ ) are resonant. Note that if there are separate populations of modes with different characteristics that are responsible for scattering – as is the case in firehose-infested plasma in an Alfvén-enabling state, in which there are both oblique firehose and secondary parallel firehose modes – a collision operator associated with both populations is required.

Finally, to be able to write down simple expressions for $\nu _{\textrm {eff}}(v,\xi )$ and $\tilde {v}_{\textrm {wv}}(v,\xi )$ , we make one final assumption: that the fluctuations can be treated as being quasi-parallel in the sense that, for most particles, $v_{\perp }^2 \ll \varOmega _i^2/k_{\perp }^2$ . The assumption simplifies the sum in (6.3): all terms with $|n| \gt 1$ are then negligible, and the Bessel functions in the $n = \pm 1$ terms can be simplified using the identity $\textrm {J}_{\pm 1}(x) \approx \pm (x/2)(1-x^2/8+\ldots )$ for $x \ll 1$ . Under this final assumption, the effective pitch-angle scattering operator $\mathfrak{C}[f]$ associated with firehose modes in the Alfvén-enabling regime simplifies to

(6.5) \begin{equation} \mathfrak{C}[f] = \frac {1}{2} \frac {\partial }{\partial \xi }\Bigg \{(1-\xi ^2)\nu _{\textrm {eff,pl}}(v \xi )\left [\frac {\partial f}{\partial \xi } - 2 \tilde {v} \tilde {v}_{\textrm {wv,pl}}(v \xi ) f_{\textrm {M}}\right ] + (1-\xi ^2)\nu _{\textrm {eff,ob}}(v \xi ) \frac {\partial f}{\partial \xi } \Bigg \} , \end{equation}

where the velocity-dependent pitch-angle scattering rates $\nu _{\textrm {eff,pl}}$ and $\nu _{\textrm {eff,ob}}$ associated with the secondary parallel firehose modes and the oblique firehose modes, respectively, are now only functions of the parallel particle velocity $v_{\|} = v \xi$ ; they are directly related to the magnetic-energy spectra of the two firehose populations by

(6.6a) \begin{align} \nu _{\textrm {eff,pl}}(v_{\|}) & \simeq \frac {\pi }{2} \frac {\varOmega _i^2}{v_{\|}}\frac { E_{B,\mathrm{pl}}\!\left ({\varOmega _i}/{v_{\|}}\right )}{B_0^2/8\pi } ,\\[-10pt]\nonumber \end{align}

(6.6b) \begin{align} \nu _{\textrm {eff,ob}}(v_{\|}) & \simeq \frac {\pi }{2} \frac {\varOmega _i^2}{v_{\|}}\frac { E_{B,\mathrm{ob}}\!\left ({\varOmega _i}/{v_{\|}}\right )}{B_0^2/8\pi } . \qquad \end{align}

Here, $E_{B,\mathrm{pl}}(k_{\|})$ and $E_{B,\mathrm{ob}}(k_{\|})$ are the one-dimensional parallel magnetic-energy spectra of the secondary parallel and oblique firehose fluctuations, respectively, while

(6.7) \begin{equation} \tilde {v}_{\textrm {wv,pl}}(v_{\|}) = \frac {{v}_{\textrm {wv,pl}}(v_{\|})}{v_{\mathrm{th}i}} \simeq \varpi \bigl ({\varOmega _i}/{v_{\|}}\bigr ) \frac {v_{\|}}{v_{\mathrm{th}i}} \end{equation}

is an approximation (to leading order in the small parameter $1/\beta _{i}$ ) of the parallel phase velocity of the modes with which particles having parallel velocity $v_{\|}$ are resonant. Because the oblique firehose modes do not have a parallel phase velocity, the pitch-angle scattering operator associated with them is already in the plasma rest frame. The quasi-parallel assumption is reasonable for the secondary parallel firehose modes, but is less clearly appropriate for the oblique firehose modes. For the latter case, we estimate the error introduced in this approximation by using the numerical result that, in the saturated state of the firehose instability, $k_{\perp } \lesssim 0.5 \rho _i^{-1}$ . It follows that the magnitude of the first-order term in the Bessel function expansion is $k^2 v_{\perp }^2/8 \varOmega _i^2 \approx v_{\perp }^2/16 v_{\mathrm{th}i}^2$ . For particles with $v_{\perp } \lesssim 2 v_{\mathrm{th}i}$ (the majority of thermal particles), the error introduced by the approximation is therefore 25 % or less.

Thus we have constructed a simple model for the effective firehose collision operator that takes as its inputs two velocity-dependent scattering rates ( $\nu _{\textrm {eff,pl}}(v_{\|})$ and $\nu _{\textrm {eff,ob}}(v_{\|})$ ) and the parallel phase velocity $v_{\textrm {wv,pl}}(v_{\|})$ of the secondary firehose modes. The scattering rates are given directly by the one-dimensional parallel magnetic-energy spectra of oblique and secondary parallel firehose modes $E_{B,\mathrm{ob}}(k_{\|})$ and $E_{B,\mathrm{pl}}(k_{\|})$ , respectively, while $v_{\textrm {wv,pl}}(v_{\|})$ depends on the real frequency $\varpi (k_{\|})$ of the secondary firehose modes. Therefore, to compute the effective firehose collision operator, all that remains is to determine $E_{B,\mathrm{ob}}(k_{\|})$ , $E_{B,\mathrm{pl}}(k_{\|})$ and $\varpi (k_{\|})$ . We compute these functions numerically for all of the expanding-box simulations that we have conducted that attain Alfvén-enabling states. In order to obtain a time-averaged collision operator, for each simulation we choose a time interval during which the firehose instability has saturated, and then calculate averaged values of the oblique and parallel magnetic-energy spectra and the secondary firehose mode frequencies.Footnote ¹⁰

To calculate $E_{B,\mathrm{pl}}(k_{\|})$ , we first apply a mask to the total (time-averaged) magnetic-energy spectrum $E_{B}(k_{\|},k_{\perp })$ to isolate the secondary parallel firehose modes; this mask covers the same region of $(k_{\|},k_{\perp })$ space as that circumscribed by the white dashed line in figure 9. We then integrate the masked spectra over all perpendicular wavenumbers to obtain $E_{B,\mathrm{pl}}(k_{\|})$ . Then $E_{B,\mathrm{ob}}(k_{\|})$ is calculated by subtracting $E_{B,\mathrm{pl}}(k_{\|})$ from the total parallel one-dimensional magnetic-energy spectrum $E_{B}(k_{\|}) \equiv \int \textrm {d} k_{\perp } \,E_{B}(k_{\|},k_\perp )$ . We show $E_{B,\mathrm{pl}}(k_{\|})$ , $E_{B,\mathrm{ob}}(k_{\|})$ and $E_{B}(k_{\|})$ from three representative simulations in Alfvén-enabling states in figure 18 $(a$ – $c)$ , respectively.

Figure 18. $(a)$ One-dimensional (parallel) magnetic-energy spectrum $E_B(k_{\|})$ of all firehose fluctuations in the saturated, Alfvén-enabling state of run CV (solid red line). Also plotted are the magnetic-energy spectra of non-quasi-parallel fluctuations (blue solid line) and quasi-parallel ones (black solid line), as well as fits for these spectra discussed in the main text (dashed lines; (6.8b )). $(b)$ Same as $(a)$ , but for run DVI. $(c)$ Same as $(a)$ , but for run BIV.

Having calculated $E_{B,\mathrm{pl}}(k_{\|})$ and $E_{B,\mathrm{ob}}(k_{\|})$ numerically, we then fit both spectra with simple analytical functions of the form

(6.8a) \begin{align} E_{B,\mathrm{pl}}\!\left (k_{\|}\right ) & \simeq \frac {B_0^2}{8\pi }\left \{\frac { \bar {E}_{B,\mathrm{pl}}}{\sqrt {\pi }\Delta {k}_{\|,\textrm {pl}} \rho _i} \exp {\left [-\frac {(k_{\|}-{k}_\|, \textrm {pl})^2 }{\Delta {k}_{\|,\textrm {pl}}^{2}}\right ]} + H(k_{\|} - {k}_\|, \textrm {pl}) \frac {\bar {E}_{B,\mathrm{tail}}}{(k_{\|} \rho _i)^{p_{\textrm {tail}}}} \right \}\!, \qquad \quad \end{align}

(6.8b) \begin{align} E_{B,\mathrm{ob}}\!\left (k_{\|}\right ) & \simeq \frac {B_0^2}{8\pi }\frac {\bar {E}_{B,\mathrm{ob}}}{\sqrt {\pi }\Delta {k}_{\|,\textrm {ob}} \rho _i} \exp {\left [-\frac {(k_{\|} -{k}_\|, \textrm {ob})^2}{\Delta {k}_{\|,\textrm {ob}}^2}\right ]} , \qquad \end{align}

where ${k}_{\|, \textrm {pl}}$ ( ${k}_{\|, \textrm {ob}}$ ) is the wavenumber at which $E_{B,\mathrm{pl}}(k_{\|})$ ( $E_{B,\mathrm{ob}}(k_{\|})$ ) attains its maximum, $\Delta {k}_{\|,\textrm {pl}}$ ( $\Delta {k}_{\|,\textrm {ob}}$ ) is the characteristic width of the $k_{\|}$ interval over which $E_{B,\mathrm{pl}}(k_{\|})$ ( $E_{B,\mathrm{ob}}(k_{\|})$ ) extends and $\bar {E}_{B,\mathrm{pl}}$ ( $\bar {E}_{B,\mathrm{ob}}$ ) is the total energy in the secondary parallel (oblique) firehose fluctuations. We also find it necessary to model the high- $k_{\|}$ wavenumber of the distribution of secondary parallel firehose modes with a power-law tail (of amplitude $\bar {E}_{B,\mathrm{tail}}$ , and power-law index $p_{\textrm {tail}}$ ); although the magnetic energy associated with modes of such high wavenumbers is much smaller than modes with $k_{\|} \rho _i \sim 1$ , such modes nonetheless play a key role in determining the anisotropy of the ion distribution function in Alfvén-enabling states (see § 6.3.1), and so cannot be disregarded.

These particular functional fits are not derived analytically, but we find empirically that they describe the numerical spectra well. Practically, we first determine $\bar {E}_{B,\mathrm{pl}}$ and $\bar {E}_{B,\mathrm{ob}}$ by integrating each spectra, and then determine best-fit values to the other parameters, weighting the estimates by $E_{B,\mathrm{pl}}(k_{\|})$ and $E_{B,\mathrm{ob}}(k_{\|})$ , respectively. For determining the fits for the oblique firehose spectrum, we exclude all parallel wavenumbers $k_{\|} \lt 0.4 \rho _i^{-1}$ , because the spectrum of these longer-wavelength fluctuations is not well described by an analytic fit of the form (6.8b ), and such modes do not affect the anisotropy of thermal ions. We fit the power-law tail of the spectrum of parallel modes by first fitting the latter’s peak with the Gaussian analytic form, then subtracting this fit and the spectrum of the noise from the total spectrum, and fitting the power law to what remains. In figure 18 $(a$ , $b)$ , the good agreement between our fits of the form (6.8b ) to the one-dimensional magnetic-energy spectra of two representative simulations is illustrated.

Figure 19. $(a)$ Best-fit estimates for wavenumber parameters introduced in (6.6b ) for all of our Alfvén-enabling simulations as a function of the expansion time. $(b)$ Same as $(a)$ , but as a function of $\beta _i$ . $(c)$ Best-fit estimates for spectral amplitude parameters introduced in (6.6b ) for all of our Alfvén-enabling simulations as a function of the expansion time. $(d)$ Same as $(c)$ , but as a function of $\beta _i$ . $(e)$ Best-fit estimates for high-wavenumber tail of parallel modes introduced in (6.6b ) for all of our Alfvén-enabling simulations as a function of the expansion time. $(f)$ Same as $(e)$ , but as a function of $\beta _i$ .

The wavenumber parameters of our best-quality fits for all of our simulations of Alfvén-enabling states as functions of $\tau _{\textrm {exp}} \varOmega _{i}$ and $\beta _i$ are presented in figure 19. We find that all of the wavenumber parameters are approximately independent of both $\tau _{\textrm {exp}} \varOmega _i$ and $\beta _i$ , save for $\Delta {k}_{\|,\textrm {pl}} \rho _i$ , which has a weak dependence on $\tau _{\textrm {exp}} \varOmega _i$ :Footnote ¹¹

(6.9) \begin{equation} {k}_{\|, \textrm {pl}} \rho _i \approx 1.2 , \quad {k}_{\|, \textrm {ob}} \rho _i \approx 0.75 , \quad \Delta {k}_{\|,\textrm {pl}} \rho _i \approx \frac {1.3}{(\tau _{\textrm {exp}}\varOmega _i)^{0.25}} , \quad \Delta {k}_{\|,\textrm {ob}} \rho _i\approx 0.19 . \end{equation}

By contrast, both $\bar {E}_{B,\mathrm{ob}}$ and $\bar {E}_{B,\mathrm{pl}}$ do depend on $\tau _{\textrm {exp}} \varOmega _i$ and $\beta _i$ , with those relationships being well approximated by the following scalings:

(6.10) \begin{equation} \bar {E}_{B,\mathrm{pl}} \approx 0.3 \frac {\beta _{i}^{0.25}}{(\tau _{\textrm {exp}} \varOmega _i)^{0.5}} , \quad \bar {E}_{B,\mathrm{ob}} \approx 0.7 \frac {\beta _{i}}{\tau _{\textrm {exp}} \varOmega _i} . \end{equation}

Finally, for the high-wavenumber component of the secondary parallel firehose modes, we find that the power-law index is approximately independent of both $\tau _{\textrm {exp}} \varOmega _i$ and $\beta _i$ , but its amplitude has a comparable scaling to the peak amplitude of the secondary parallel firehose modes:

(6.11) \begin{equation} {p}_{\mathrm{tail}} \approx 2.7 , \quad \bar {E}_{B,\mathrm{tail}} \approx \frac {0.09}{(\tau _{\textrm {exp}} \varOmega _i)^{0.5}} . \end{equation}

We discuss possible theoretical justifications for these scalings in § 6.4.

To calculate the dispersion relation $\varpi (k_{\|})$ of the secondary firehose modes in our simulation, as well as to confirm that the oblique firehose modes are non-propagating, we compute the frequency-dependent magnetic-energy spectra $E_B(k_{\|},k_{\perp },\varpi )$ of the firehose fluctuations in saturation of our simulations of Alfvén-enabling states. Figure 20 $(a)$ shows $E_B(k_{\|},k_{\perp },\varpi )$ computed for a representative simulation at two fixed values of $k_{\perp }$ : for purely parallel modes ( $k_{\perp } = 0$ ) and for oblique modes with $k_{\perp }$ comparable to that of the oblique firehose modes.

Figure 20. $(a)$ Slice plots at fixed $k_{\perp }$ of the frequency-dependent magnetic-energy spectrum $E_B(k_{\|},k_{\perp },\varpi )$ of the firehose fluctuations in run CV, averaged over the saturated state. $(b)$ Fluctuation-energy-weighted average value of real frequency $\varpi _{\textrm {sim}}$ of the firehose fluctuations as a function of $k_{\|} \rho _i$ for several different runs that attain Alfvén-enabling states. The black line denotes the approximate fit (6.13).

For the parallel modes ( $k_\perp \rho _i=0$ ; figure 20 a, left panel), three distinct wave populations can be identified: at small wavenumbers ( $k_{\|} \rho _i \lt 0.4$ ), both left- and right-handed Alfvén modes, while at ion-Larmor scales, secondary parallel firehose modes. As expected, the latter do indeed have a non-zero real frequency. For the oblique modes ( $k_\perp \rho _i=0.4$ ; figure 20 a, right panel), we also observe three distinct populations: at long-wavelengths ( $k_{\|} \rho _i \lt 0.4$ ), shear Alfvén modes; just above ion-Larmor scales ( $k_{\|} \rho _i \in [0.4,0.9]$ ), zero-frequency oblique firehose modes; and a weak population of (propagating) oblique secondary parallel firehose modes for $k_{\|} \rho _i \gt 0.9$ . The presence of the long-wavelength modes in addition to the secondary parallel firehose and oblique firehose modes is perhaps surprising, because such long-wavelength modes are linearly damped at these levels of pressure anisotropy; we postulate that it is nonlinear coupling between secondary parallel firehose and oblique firehose modes that gives rise to them.

From $E_B(k_{\|},k_{\perp },\varpi )$ , we can then obtain a numerical estimate of the dispersion relation of the parallel secondary firehose modes as a function of $k_{\|}$ by taking a weighted mean:

(6.12) \begin{equation} \langle \varpi \rangle (k_{\|}) \equiv \frac {\int _0^{k_{\perp ,\textrm {max}}} \textrm {d} k_{\perp } \int _0^{\varpi _{\textrm {max}}} \textrm {d}\varpi \, \varpi E_{B}(k_{\|},k_\perp , \varpi )}{\int _0^{k_{\perp ,\textrm {max}}} \textrm {d} k_{\perp } \int _0^{\varpi _{\textrm {max}}} \textrm {d}\varpi \, E_{B}(k_{\|},k_\perp , \varpi )}. \end{equation}

Because the parallel secondary firehose modes are the dominant ones at $k_{\|} \rho _i \gt 0.9$ , the dependence of $\langle \varpi \rangle$ on $k_{\|}$ will correspond to their dispersion relation. This numerical estimate for several representative simulations is shown in figure 20 $(b)$ . We find that, for $k_{\|} \rho _i \in [0.7,1.7]$ , $\langle \varpi \rangle (k_{\|})$ is well approximated by the fit

(6.13) \begin{equation} \langle \varpi \rangle (k_{\|}) \approx \frac {\varOmega _i}{\beta _i} (4.9 k_{\|} \rho _i - 2.9). \end{equation}

The relatively narrow range of wavenumbers over which the spectrum of secondary parallel firehose modes exists means that the clearly unphysical parts of this fit to the dispersion relation (i.e. $k_{\|} \rho _i \ll 1$ , where $\varpi$ goes negative) are never used.

6.3. Testing the model collision operator

Having proposed an effective collision operator associated with firehose fluctuations in an Alfvén-enabling state, we now test whether this operator is consistent with two observables from our simulations: firstly, the velocity-dependent anisotropy of the distribution function for particles with speeds comparable to the thermal speed; secondly, Fokker–Planck coefficients calculated directly from the evolution of a subpopulation of tracked (macro)particles.

6.3.1. Velocity-dependent anisotropy of the distribution function

Because collision operators describe how specific collisional processes affect the distribution function, confirming that our proposed firehose collision operator accounts for the observed distribution function’s anisotropy is a natural test of our model. In the case of an expanding, high- $\beta _i$ plasma in an Alfvén-enabling state whose constituent particles have an effective collision rate $\nu _{\textrm {eff}}$ that satisfies $\tau _{\textrm {exp}}^{-1} \ll \nu _{\textrm {eff}} \ll \varOmega _i$ , the relationship between the anisotropy of the distribution function and the collision operator takes a simple form. This condition is expected to hold for most particles in firehose-susceptible plasmas that attain Alfvén-enabling states, because the velocity-averaged collisionality $\langle \nu _{\textrm {eff}} \rangle$ of particles satisfies $\langle \nu _{\textrm {eff}} \rangle \sim \beta _i/\tau _{\textrm {exp}} \gg 1/\tau _{\textrm {exp}}$ , while the pitch-angle scattering rate of even the most frequently scattered particles obeys the bound $\nu _{\textrm {eff}} \ll \varOmega _i$ .

To establish a relationship between the distribution function’s anisotropy and the effective collision operator, we employ a modified version of a mathematical technique used in classical transport theory of plasmas: the Chapman–Enskog expansion (e.g. Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025). This technique assumes that, in plasmas where the collision rate greatly exceeds the macroscopic evolution rate, the distribution function in the expanding plasma can be expanded in the form

(6.14) \begin{equation} f_i = f_{0i} + f_{1i} + \cdots\, , \end{equation}

where the first-order correction $f_{1i} \sim f_{0i}(\nu _{\textrm {eff}} \tau _{\textrm {exp}})^{-1}$ is asymptotically small compared with the leading-order term $f_{0i}$ . Simultaneously, the condition that $\nu _{\textrm {eff}} \ll \varOmega _i$ means that, over the evolution time scales of interest, $f_i$ is approximately gyrotropic. If the gyroaveraged kinetic equation satisfied by the distribution function is also expanded in the small parameter $(\nu _{\textrm {eff}} \tau _{\textrm {exp}})^{-1}$ , we find that, to leading order, $f_{0i}$ must satisfy $\mathfrak{C}_{\textrm {f}}[f_{0i}] = 0$ . Adopting our model firehose collision operator, and taking into account the Maxwellian initial condition of the distribution function in our simulations, this equation has the unique solution $f_{0i} = f_{\mathrm{M}i}$ .Footnote ¹² Considering the equation that arises to next order from the gyroaveraged kinetic equation, it follows that

(6.15) \begin{equation} \mathfrak{C}_{\textrm {f}}[f_{1i}] = \left [\left (\hat {{\boldsymbol{b}}}\hat {{\boldsymbol{b}}}-\frac {1}{3}\unicode{x1D644}\right )\,{\boldsymbol{:}}\,\unicode{x1D652}_i \right ] \tilde {v}^2 P_2(\xi ) f_{\mathrm{M}i} , \end{equation}

where $\unicode{x1D652}_i$ is the (traceless, symmetric) rate-of-strain tensor of the ion bulk flow, $P_2(\xi )$ is the Legendre polynomial of second degree and we remind the reader that $\tilde {v} \equiv v/v_{\mathrm{th}i}$ . In the case of plasma that is linearly expanding on a time scale $\tau _{\textrm {exp}}$ in one direction that is perpendicular to the background magnetic field, $\unicode{x1D652}_i = 2 \boldsymbol{\nabla } {\boldsymbol{u}} = -(2/\tau _{\textrm {exp}})\hat {{\boldsymbol{x}}}\hat {{\boldsymbol{x}}}$ , and so (6.15) becomes

(6.16) \begin{equation} \mathfrak{C}_{\textrm {f}}[f_{1i}] = -\frac {2}{3 \tau _{\textrm {exp}}} \tilde {v}^2 P_2(\xi ) f_{\mathrm{M}i}. \end{equation}

Now assuming that $\mathfrak{C}_{\textrm {f}}[f_{1i}] = \mathfrak{C}_{\textrm {f}}[f_{i}]$ takes the form given by (6.5), and integrating (6.16) from $\xi = -1$ to $\xi _0 \rightarrow \xi$ , we deduce that

(6.17) \begin{equation} \frac {1}{2 \tilde {v}}\frac {\partial f_{i}}{\partial \xi } \approx f_{\mathrm{M}i}(v) \left \{\frac {\nu _{\textrm {eff,pl}}(\tilde {v}_{\|})\tilde {v}_{\mathrm{wv}}(\tilde {v}_{\|})}{\nu _{\textrm {eff,pl}}(\tilde {v}_{\|})+\nu _{\textrm {eff,ob}}(\tilde {v}_{\|})} + \frac {\tilde {v}_{\|}}{3 \tau _{\textrm {exp}}[\nu _{\textrm {eff,pl}}(\tilde {v}_{\|})+\nu _{\textrm {eff,ob}}(\tilde {v}_{\|})]}\right \}\!. \end{equation}

Thus, we have established a simple relationship between the pitch-angle gradient of the distribution function and the functions $\nu _{\textrm {eff,pl}}(\tilde {v}_{\|})$ , $\nu _{\textrm {eff,ob}}(\tilde {v}_{\|})$ and $\tilde {v}_{\mathrm{wv}}(\tilde {v}_{\|})$ that characterise our model firehose collision operator.

Figure 21. $(a)$ Slice plot of the pitch-angle gradient of the ‘saturated’ ion distribution function $f_{i,\mathrm{sat}}$ divided by $2 \tilde {v} f_{\textrm {M}}$ in an Alfvén-enabling state (run CV). Here, $f_{i,\mathrm{sat}}$ is the domain-averaged ion distribution function $f(v_{\|},v_{\perp })$ time-averaged over the saturated period of the firehose instability. $(b)$ The same quantity as $(a)$ , but averaged over $v_{\perp }$ (solid red line); the average excludes the shaded region shown in $(a)$ which is negatively influenced by poor particle statistics. This is plotted with the theoretical prediction (6.17) for this quantity arising from our proposed collision operator (dashed line), and in the absence of any collisions (dotted line). The dark blue (grey) line shows the solutions of (6.19) at $t = 0.34 \tau _{\textrm {exp}}$ , including (excluding) the high-wavenumber power-law tail of secondary parallel firehose modes. Inset: complex frequency $\omega$ plotted against parallel wavenumber for linear low-frequency modes arising in a plasma with ion distribution function $f_{i,\mathrm{sat}}$ . The real frequency $\varpi$ (growth rate $\gamma$ ) is shown in red (blue). Also plotted is the observed real frequency $\varpi _{\textrm {sim}}$ of firehose fluctuations in the same run (red dashed).

Figure 21 provides a test of this relationship in the case of our asymptotic Alfvén-enabling run in its saturated phase. Firstly, figure 21 $(a)$ illustrates a key feature of (6.17): that $({\partial f_{i}}/{\partial \xi })(1/\tilde {v}f_{\mathrm{M}i})$ is approximately independent of $v_{\perp }$ , and primarily a function of $v_{\|}$ . Secondly, figure 21 $(b)$ compares the time- and $v_{\perp }$ -averaged value of $({\partial f_{i}}/{\partial \xi })(1/2\tilde {v}f_{\mathrm{M}i})$ (solid line) with the right-hand side of (6.17) (dashed line), where we first compute $\nu _{\textrm {eff,pl}}(\tilde {v}_{\|})$ , $\nu _{\textrm {eff,ob}}(\tilde {v}_{\|})$ and $\tilde {v}_{\mathrm{wv}}(\tilde {v}_{\|})$ assuming our quasi-linear model applies instantaneously, and then time-average the entire expression. The agreement is very strong, save for $|\tilde {v}_{\|}| \ll 0.5$ , supporting the claim that our model collision operator is appropriate for $|\tilde {v}_{\|}| \gtrsim 0.5$ .

To explain why reasonable agreement is not attained for comparatively small values of $\tilde {v}_{\|}$ , we note that, in deriving (6.17), we have implicitly assumed that the rate of anomalous scattering is large enough that, at the specified time $t$ at which the comparison is made, either $\nu _{\textrm {eff,pl}} t \gg 1$ or $\nu _{\textrm {eff,ob}} t \gg 1$ . As $\tilde {v}_{\|}$ is decreased from order-unity values to smaller ones, the amplitude of the increasingly high wavenumber firehose modes with which such particles are resonant decreases, leading to an ever-smaller scattering rate. Eventually, the collision rate decreases enough that $\nu _{\textrm {eff,ob}} t \ll \nu _{\textrm {eff,pl}} t \lesssim 1$ , at which point there is no expectation for (6.17) to hold. Indeed, if $\nu _{\textrm {eff,pl}} t \ll 1$ , one should expect the distribution function’s anisotropy to be consistent with continued double-adiabatic evolution; that is, in the absence of any collisions, the non-Maxwellian component of the distribution function should be given by

(6.18) \begin{equation} f_{1i} \approx \frac {2 t}{3 \tau _{\textrm {exp}}} \tilde {v}^2 P_2(\xi ) f_{\mathrm{M}i} , \quad \frac {1}{2 \tilde {v}}\frac {\partial f_{i}}{\partial \xi } \approx \tilde {v}_{\|} \frac {t}{\tau _{\textrm {exp}}} f_{\mathrm{M}i}. \end{equation}

This expression is plotted in figure 21 $(b)$ (dotted line) using the mean time of the ‘saturated’ interval over which $f_i$ is averaged; good agreement is found for $\tilde {v}_\parallel \lesssim 0.2$ . This implies that scattering of particles with small values of $v_{\|}$ compared with the ion thermal speed is indeed too infrequent to impact the distribution function anisotropy at such velocities.

We can further test this hypothesis by considering the evolution equation of the first-order correction $f_{1i}$ under the ordering $t \sim \tau _{\textrm {exp}}/\beta _{i} \sim \langle \nu _{\textrm {eff}} \rangle ^{-1}$ :

(6.19) \begin{equation} \frac {\partial f_{1i}}{\partial t} - \mathfrak{C}_{\textrm {f}}[f_{1i}] = \frac {2}{3 \tau _{\textrm {exp}}} \tilde {v}^2 P_2(\xi ) f_{\mathrm{M}i} . \end{equation}

It is clear that, taking the subsidiary limit $t \langle \nu _{\textrm {eff}} \rangle \gg 1$ recovers the steady-state solution (6.16), while the opposite limit $t \langle \nu _{\textrm {eff}} \rangle \ll 1$ returns adiabatic evolution, with $f_{1i}$ given by (6.18). We then solve (6.19) numerically, with the effective collision operator given by (6.1) and (6.2c ) when $\Delta _i \lt -1.35/\beta _i$ (i.e. when the oblique firehose is first destabilised). We integrate forward in time for the same duration as in our Pegasus++runs and compute the pitch-angle gradients. An illustrative comparison of the two results (dark blue versus red line) for run CV is shown in figure 21 $(b)$ ; in this (and other simulations) we find quantitative agreement, supporting our hypothesis.

We can also use numerical solutions of (6.19) to investigate the importance (or possible lack thereof) of the high-wavenumber power-law tail of secondary parallel firehose modes. If we remove the contribution of these modes from (6.2a ), and re-run our numerical solution of (6.19), we obtain the grey line in figure 21 $(b)$ . The resulting pitch-angle derivative of $f_{i,\textrm { sat}}$ matches the Pegasus++results well for $|v_{\|}| \gtrsim 0.6 v_{\mathrm{th}i}$ and for $|v_{\|}| \lesssim 0.2 v_{\mathrm{th}i}$ . For intermediate values of $v_{\|}$ , the numerical solution implies (erroneously) that the pitch-angle gradient of the distribution function should, for such values of $v_{\|}$ , be given by the double-adiabatic result (6.18). The reason that the double-adiabatic prediction is incorrect is simply that, if the high-wavenumber power-law tail of secondary parallel firehose modes is not modelled, then the scattering rate due to modes with $k_{\|} \rho _i \gtrsim 2$ implied by (6.2c ) is insufficient for the distribution function’s anisotropy to have been regulated in any meaningful way. We conclude that the high-wavenumber secondary firehose modes – which, as we argue in § 6.4.3, should be present physically – play a non-trivial role in determining the velocity-dependent anisotropy of the ion distribution function in saturation.