2.1.1. Model equations and assumptions

The linear evolution of the NP mode at long wavelengths can be treated analytically in the drift-kinetic approximation, in which all relevant time and length scales are much larger than those associated with the particles’ gyromotion and the velocity distribution function of the particles is gyrotropic. We adopt this framework, and further simplify the calculation by treating the electrons as a massless, neutralizing, isothermal fluid having constant temperature $T_{\rm e}$.Footnote ² In this model, the velocity of magnetic-field lines, and equivalently the perpendicular fluid flow, is captured by the $\boldsymbol {E}{\boldsymbol {\times }}\boldsymbol {B}$ drift velocity $\boldsymbol {u}_\perp$. The perpendicular velocity peculiar to this drift, denoted by $\boldsymbol {w}_{\perp }$, then describes the perpendicular particle motion relative to the field lines and the fluid flow, under the constraint that the magnetic moment $\mu \doteq m_{\rm i}w_{\perp }^2/2B$ is conserved. The component of the particle velocity directed along the local magnetic-field direction is denoted by $v_\parallel$.

In what follows, we solve for the evolution of small perturbations $\delta f(t,\boldsymbol {r},v_\parallel,w_\perp )$ to a spatially uniform ‘background’ ion velocity distribution function $F_0(v_\parallel,w_\perp )$. The parallel ($\parallel$) and perpendicular ($\perp$) coordinate directions are fixed with respect to a uniform background magnetic field, $\boldsymbol {B}_0$. Assuming that spatial variations in the plasma are due only to a sinusoidal perturbation having wavenumbers $k_{\parallel }$ and $k_{\perp }$, the relevant equations in their linearized forms are the drift-kinetic Vlasov equation,

(2.1a) \begin{equation} \left( \frac{\partial }{\partial t} + {\rm i} k_\parallel v_\parallel \right) \left( \delta f + \frac{\delta B_\parallel}{B_0} \frac{w_\perp}{2} \frac{\partial F_0}{\partial w_\perp} \right) + \frac{e}{m_{\rm i}} \delta E_\parallel \frac{\partial F_0}{\partial v_\parallel} - {\rm i} k_{{\parallel}}\frac{\delta B_{{\parallel}}}{B_0}\frac{w^2_{{\perp}}}{2} \frac{\partial F_0}{\partial v_\parallel} = 0 ; \end{equation}

the force equation for the evolution of the drift velocity,

(2.1b)\begin{equation} \frac{{\rm d} u_{{\perp}}}{{\rm d} t} ={-}\frac{{\rm i} k_{{\perp}}}{m_{\rm i}n_0} \left(\delta p_{{\perp}{\rm i}} + T_{\rm e} \delta n \right) - {\rm i} k_{{\perp}}v_{\rm A}^2\frac{\delta B_{{\parallel}}}{B_0} + {\rm i} k_{{\parallel}}v^2_{\rm A} \frac{\delta B_{{\perp}}}{B_0}; \end{equation}

the ideal induction equation governing the parallel and perpendicular components of the perturbed magnetic field $\delta \boldsymbol {B}$,

(2.1c)\begin{equation} \frac{{\rm d} }{{\rm d} t} \frac{\delta B_{{\parallel}}}{B_0} ={-}{\rm i} k_{{\perp}} u_{{\perp}} \quad{\rm and}\quad \frac{{\rm d} }{{\rm d} t} \frac{\delta B_\perp}{B_0} = {\rm i} k_\parallel u_\perp ; \end{equation}

and a generalized Ohm's law for the parallel electric field,

(2.1d)\begin{equation} \delta E_{{\parallel}} ={-}{\rm i} k_\parallel \frac{T_{\rm e}}{e} \frac{\delta n}{n_0} . \end{equation}

The perturbed number density and perpendicular ion pressure are given by

(2.2a,b)\begin{equation} \delta n \doteq \int{\rm d}^3\boldsymbol{v} \, \delta f \quad{\rm and}\quad \delta p_{{\perp}{\rm i}} \doteq \int{\rm d}^3\boldsymbol{v} \, \frac{1}{2} m_{\rm i} w^2_\perp \delta f , \end{equation}

respectively, with ${\rm d}^3\boldsymbol {v} = 2{\rm \pi} w_\perp {\rm d} w_\perp {\rm d} v_\parallel$. The other symbols have their usual meanings: $e$ is the elementary charge, $m_{\rm i}$ is the ion mass, and $v_{\rm A} \doteq B_0/(4{\rm \pi} m_{\rm i} n_0)^{1/2}$ is the Alfvén speed given $B_0$ and a uniform background density $n_0$ (the zeroth moment of $F_0$). Note that $u_\perp$ is not an explicit moment of the perturbed distribution function, and must be evolved independently using (2.1b). This combination of the drift-kinetic equation with a fluid equation for the drift velocity and a frozen-in magnetic field is commonly referred to as ‘kinetic MHD’ (Kulsrud Reference Kulsrud1964, Reference Kulsrud1983).

At this point, we take $F_0$ to be a stationary, isotropic, Maxwell–Boltzmann distribution, $F_0=F_{\rm M}(v)$, with $\int {\rm d}^3\boldsymbol {v} F_{\rm M}(v)=n_0$ and $\int {\rm d}^3\boldsymbol {v}\, m_{\rm i}v^2 F_{\rm M}(v) = 3n_0 T_{{\rm i}0} \doteq 3 p_{{\rm i}0}$. This not only simplifies the analysis, but also ensures that the background distribution function itself is not kinetically unstable. Equation (2.1a) can then be readily integrated in time to obtain

(2.3)$$\begin{gather} \delta f(t,w_\perp,v_\parallel) = \delta f(0,w_\perp,v_\parallel)\exp\left({-{\rm i} k_\parallel v_\parallel t}\right) \nonumber\\ - \int^t_0{\rm d} t' F_{\rm M}(v) \exp\left({-{\rm i} k_\parallel v_\parallel (t-t')}\right) \left[ {\rm i} k_\parallel v_\parallel \frac{T_{\rm e}}{T_{{\rm i}0}} \frac{\delta n(t')}{n_0} - \frac{w^2_\perp}{v^2_{{\rm th},{\rm i}}} \frac{{\rm d} }{{\rm d} t'} \frac{\delta B_\parallel(t')}{B_0} \right] , \end{gather}$$

where $v_{{{\rm th}},{\rm i}} \doteq (2T_{{\rm i}0}/m_{\rm i})^{1/2}$ is the ion thermal speed. The first term on the right-hand side of (2.3) represents the parallel phase mixing of the initial perturbation by the free streaming of particles along the (unperturbed) magnetic field. If $\delta f(0,w_\perp,v_\parallel ) \propto F_{\rm M}(v)$, integrating this term and then completing the square shows that any moment of the initially perturbed distribution function will decay as $\exp [-(k_\parallel v_{{\rm th},{\rm i}} t/2)^2]$. The second term in (2.3) captures the self-consistent response of the plasma to the induced parallel electric field (${\propto }\delta n/n_0$) and the magnetic mirror force (${\propto }\delta B_\parallel /B_0$). It is this eigenmode response that we first calculate and discuss, before moving on to take the second moments of (2.3) and compute the time-dependent pressure anisotropy in § 2.1.3.

2.1.2. Eigenmode response for the NP mode

If we take the fluctuation amplitudes to be proportional to $\exp (-{\rm i}\omega t)$ with complex frequency $\omega$, the dispersion relation that results after combining (2.1) may be written as

(2.4) \begin{align} D(\zeta) \doteq \left( \omega^2 - k^2 v^2_{\rm A} \right) \left[ 1 + \frac{T_{{\rm i}0}}{T_{\rm e}} + \zeta Z(\zeta) \right] + k^2_\perp v^2_{{\rm th},{\rm i}} \, \zeta Z(\zeta) \left[ 1 + \frac{T_{{\rm i}0}}{T_{\rm e}} + \frac{1}{2} \zeta Z(\zeta) \right] = 0 , \end{align}

where $\zeta \doteq \omega /|k_\parallel |v_{{\rm th},{\rm i}}$ is the dimensionless phase speed and $Z(\zeta )$ is the plasma dispersion function. The first term in parentheses captures the combined restoring force of the magnetic pressure and tension, and indicates that we are examining magnetosonic modes. Indeed, setting the accompanying multiplicative term in square brackets to zero provides the dispersion relation for a Landau-damped IAW in the limit $(m_{\rm e}/m_{\rm i})^{1/2}\ll 1$. The final term in (2.4), proportional to $k^2_\perp v^2_{{\rm th},{\rm i}}$, couples these Alfvénic and acoustic responses; its presence can be traced back to the final term in (2.3) representing the mirror force, and thus introduces collisionless damping of the mode through transit-time damping.

To isolate the NP mode, we focus specifically on highly oblique wavenumbers ($k_{\perp } \gg k_{\parallel }$) and low frequencies ($\zeta \ll 1$). In this limit, the plasma dispersion function in (2.4) can be approximated as $Z(\zeta ) \approx {\rm i}\sqrt {{\rm \pi} }$, and we may simplify the dispersion relation further by neglecting terms of order $\zeta ^2$. The result is an approximate expression for the decay rate of the NP mode:

(2.5)\begin{equation} \zeta \simeq{-}\frac{{\rm i}}{\sqrt{\rm \pi}\beta_{{\rm i}0}} \frac{k^2}{k_{{\perp}}^2}, \quad{\rm where}\ \beta_{{\rm i}0} = \frac{8{\rm \pi} p_{{\rm i}0}}{B^2_0} = \frac{v^2_{{\rm th},{\rm i}}}{v^2_{\rm A}} . \end{equation}

For $\zeta \ll 1$ to be satisfied by (2.5), we require that $\beta _{{\rm i}0} \gg 1$, which aligns well with our interest in high-$\beta$ plasmas. Further properties of the NP mode can be found by taking moments of the kinetic equation (2.1a), such as the proportionalities between $\delta n$, $\delta p_{\perp,{\rm i}}$ and $\delta B_\parallel$:

(2.6a)\begin{gather} \frac{\delta n}{n_0} ={-} \zeta Z(\zeta) \left[ 1 + \frac{T_{\rm e}}{T_{{\rm i}0}}\left( 1 + \zeta Z(\zeta)\right)\right]^{{-}1} \frac{\delta B_\parallel}{B_0} \simeq{-}\frac{1}{\beta_0}\frac{k^2}{k^2_\perp} \frac{\delta B_\parallel}{B_0} , \end{gather}

(2.6b)\begin{gather} \frac{\delta p_{{\perp} {\rm i}}}{p_{{\rm i}0}} ={-}\frac{T_{\rm e}}{T_{{\rm i}0}} \frac{\delta n}{n_0} + 2\, \frac{\omega^2 - k^2 v^2_{\rm A}}{k^2_\perp v^2_{{\rm th},{\rm i}}}\frac{\delta B_\parallel}{B_0} \simeq \left(2+\frac{T_{\rm e}}{T_{{\rm i}0}}\right) \frac{\delta n}{n_0}, \end{gather}

where $\beta _0 \doteq \beta _{{\rm i}0} ( 1 + T_{\rm e}/T_{{\rm i}0})$. The latter equation implies approximate perpendicular pressure balance when $k_\parallel \ll k_\perp$, since then

(2.6c)\begin{equation} \delta p_{{\perp} {\rm i}} + \delta p_{\rm e} + \frac{\delta B^2}{8{\rm \pi}} \simeq{-}\frac{k^2_\parallel}{k^2_\perp}\frac{\delta B^2}{4{\rm \pi}} \ll \frac{\delta B^2}{4{\rm \pi}}. \end{equation}

Additionally, the parallel ion pressure perturbation is given by

(2.6d)\begin{equation} \frac{\delta p_{{\parallel}{\rm i}}}{p_{{\rm i}0}} ={-}\frac{T_{\rm e}}{T_{{\rm i}0}} \frac{\delta n}{n_0} - 2 \zeta^2 \left( 1+\zeta Z(\zeta) \right) \left( \frac{\delta B_\parallel}{B_0} + \frac{T_{\rm e}}{T_{{\rm i}0}}\frac{\delta n}{n_0}\right) \simeq{-}\frac{T_{\rm e}}{T_{{\rm i}0}} \frac{\delta n}{n_0} , \end{equation}

so that $\delta p_{\parallel {\rm i}} + \delta p_{\rm e} \simeq 0$. Equations (2.5) and (2.6) highlight some of the essential properties of the NP mode, namely, that it does not oscillate, that it decays slowly at high $\beta$, and that its perturbations to the magnetic-field strength and the density are anti-correlated.

The physical mechanism behind the damping rate is primarily transit-time magnetic pumping, in which Landau-resonant particles (technically, their guiding centres) that are trapped between large-scale magnetic mirrors formed by an oblique perturbation in the magnetic field extract energy from the mirror force. They experience net heating by betatron acceleration because the number of particles heated in regions where $|B|$ increases (lower $v_\parallel$ particles) is greater than the number of particles cooled where $|B|$ decreases (higher $v_\parallel$ particles). At higher plasma $\beta$, this difference is smaller, hence the $\beta ^{-1}$ dependence of the damping rate.Footnote ³

This type of collisionless damping is susceptible to nonlinear saturation, whereby the particles in the well explore the phase space available to them by $\mu$ conservation, phase-mixing out the original Maxwellian according to their differing bounce times and flattening the distribution function in the magnetic well to create a plateau around $v_\parallel \sim 0$. This effectively increases the plasma $\beta$ of the resonant particles and the damping rate weakens dramatically. Because of the slow nature of the NP mode's decay rate at high $\beta$, nonlinear saturation occurs comparatively rapidly, at a rate comparable to the bounce frequency of a thermal particle,

(2.7)\begin{equation} \varOmega_{{{\rm b}}} \doteq \frac{1}{2} k_{{\parallel}} v_{\text{th},\text{i}} \left|\frac{\delta B_\parallel}{B_0}\right|^{1/2} . \end{equation}

While most particles bounce at approximately this frequency, particles that are barely trapped bounce much slower due to their prolonged time spent traversing the edge of the magnetic well. As a result, the plateau forms inside-out, reaching a pitch-angle dependent maximum extent set by $|v_\parallel |/w_\perp < \sqrt {|B_{\max }|/|B_{\min }|-1}$. For $|\delta B_\parallel /B_0|\gtrsim \beta ^{-2}_{{\rm i}0}$, the bounce frequency in (2.7) will be larger than the decay rate in (2.5), and thus nonlinear saturation will be important. Because of our interest in plasmas having $\beta \gg 1$, even modes that may often be considered ‘linear’ in amplitude will thus decay by only a small amount before experiencing nonlinear saturation, the implication being that these structures should be long lived. That is, unless some process is able to erode the resonant plateau in the perturbed distribution function on a time scale ${\lesssim }\varOmega ^{-1}_{\rm b}$.

2.1.3. Generation of pressure anisotropy and triggering of the mirror instability

The eigenmode in (2.6) implies a dimensionless pressure anisotropy in the ions given by

(2.8)\begin{equation} \Delta_{\rm NP} \simeq 2 \left(1+\frac{T_{\rm e}}{T_{{\rm i}0}}\right) \frac{\delta n}{n_0} \simeq{-}\frac{2}{\beta_{{\rm i}0}}\frac{k^2}{k^2_\perp}\frac{\delta B_\parallel}{B_0} . \end{equation}

This suggests that, for $\delta B_\parallel /B_0 \sim 1$, the pressure anisotropy associated with the NP mode is sufficient to excite both the firehose and mirror instabilities, the former occurring in regions where $\delta B_\parallel > 0$, the latter occurring in regions where $\delta B_\parallel < 0$. There are two considerations that complicate this conclusion.

The first complication concerns the additional pressure anisotropy that is generated when the initial perturbation to the distribution function is anisotropically phase mixed by particles streaming freely along, but not across, the field lines. To see this effect, let us return to the time-dependent solution for the perturbed distribution function, (2.3), and suppose that, at $t=0$, the plasma hosts an isothermal, pressure-balanced perturbation with

(2.9)\begin{equation} \delta f(0,w_\perp,v_\parallel) = \frac{\delta n(0)}{n_0} F_{\rm M}(v) ={-} \frac{2}{\beta_0} \frac{\delta B_\parallel(0)}{B_0} F_{\rm M}(v) . \end{equation}

This initial condition guarantees that the pressure anisotropy that develops as the particles free stream and the plasma settles into the NP eigenmode is generated self-consistently and not put in by hand. Calculating the difference of the $(1/2) m_{\rm i} w^2_\perp$ and $m_{\rm i}v^2_\parallel$ moments of (2.3) with the initial condition (2.9) yields the following expression for the time-dependent pressure anisotropy:

(2.10)\begin{align} &\Delta_{\rm NP}(t) = 2 \left(\frac{k_\parallel v_{{\rm th},{\rm i}}t}{2}\right)^2 \exp\left({-(k_\parallel v_{{\rm th},{\rm i}}t/2)^2}\right) \left(1 + \frac{T_{\rm e}}{T_{{\rm i}0}}\right) \frac{\delta n(0)}{n_0} \nonumber\\ &\qquad + \int^t_0{\rm d} t' \exp\left({-[k_\parallel v_{{\rm th},{\rm i}}(t-t')/2]^2}\right) \frac{{\rm d} }{{\rm d} t'} \frac{\delta B_\parallel(t')}{B_0} \nonumber\\ &\qquad + 2 \int^t_0{\rm d} t' \left[ \frac{k_\parallel v_{{\rm th},{\rm i}}(t-t')}{2}\right]^2 \exp\left({-[k_\parallel v_{{\rm th},{\rm i}}(t-t')/2]^2}\right) \frac{{\rm d} }{{\rm d} t'} \left[ \frac{T_{\rm e}}{T_{{\rm i}0}} \frac{\delta n(t')}{n_0} + \frac{\delta B_\parallel(t')}{B_0} \right] . \end{align}

All terms involving the combination $k_\parallel v_{{\rm th},{\rm i}}t/2$ describe the damping effect of phase mixing on the moments of the perturbed distribution function due to the production of fine-scale structure along $v_\parallel$. As discussed by Kunz et al. (Reference Kunz, Squire, Schekochihin and Quataert2020, their (3.7)), the first term on the right-hand side of (2.10) captures a transiently produced pressure anisotropy resulting from the anisotropy of particle motion: as the magnetized particles free stream along, but not across, the field, the $w^2_\perp$ and $v^2_\parallel$ moments of $\delta f(0)$ phase mix differently. The integral terms in (2.10) capture the pressure anisotropy driven by adiabatic invariance as the mode is excited and then decays in time. It is this contribution to $\Delta _{\rm NP}(t)$ that includes the pressure anisotropy of the eigenmode, (2.8).

The integrals in (2.10) can be computed numerically (see Appendix A) and the pressure anisotropy $\Delta _{\rm NP}(t)$ determined for a given initial mode amplitude

(2.11)\begin{equation} \alpha \doteq \left|\frac{\delta B_\parallel(0)}{B_0}\right| . \end{equation}

The result is shown in figure 1(a) at a selection of values of $T_{\rm e}/T_{{\rm i}0}$. The initial rise in $\Delta _{\rm NP}$ is due to a combination of the anisotropic phase mixing of the initially perturbed density and the pressure anisotropy adiabatically produced as the system settles into the NP eigenmode. After approximately one thermal-crossing time of the mode's parallel wavelength, the eigenmode is established and the slow exponential decay of $\Delta _{\rm NP}$ seen in the figure reflects the Barnes damping of the mode. (The higher-frequency oscillations seen on top of this slow decay are caused by fast modes excited by the initial conditions and represent rapid oscillations about the approximately perpendicular pressure balance.) An approximate analytic solution for $\Delta _{\rm NP}(t)$ may be obtained in the limit of $\beta _{{\rm i}0}\gg 1$, $(k_\parallel /k_\perp )^2\ll 1$ and $T_{\rm e}/T_{{\rm i}0} \sim 1$ upon substituting the damped eigenmode (2.6a) into the time integrals in (2.10). The result is that

(2.12a)\begin{align} \Delta_{\rm NP}(t) &\simeq 2 \tau^2 {\rm e}^{-\tau^2} \left( 1 + \frac{T_{\rm e}}{T_{{\rm i}0}} \right) \frac{\delta n(0)}{n_0} - \left( {\rm erf}(\tau) - \frac{\tau}{\sqrt{\rm \pi}} {\rm e}^{-\tau^2} \right) \frac{2}{\beta_{{\rm i}0}}\frac{\delta B_\parallel(t)}{B_0} \end{align}

(2.12b)\begin{align} &={-} \left[ 2 \tau^2 {\rm e}^{-\tau^2} + {\rm e}^{{-}2{\rm i}\zeta\tau}\left( {\rm erf}(\tau) - \frac{\tau}{\sqrt{\rm \pi}} {\rm e}^{-\tau^2} \right) \right] \frac{2}{\beta_{{\rm i}0}} \frac{\delta B_\parallel(0)}{B_0} , \end{align}

where $\tau \doteq k_\parallel v_{{\rm th},{\rm i}} t/2$. The term in square brackets goes as ${\sim }2\tau ^2+\tau /\sqrt {{\rm \pi} }$ for early times, suggesting that the plasma would become mirror-unstable at a time $t_{\rm m} \sim (\sqrt {\alpha } k_\parallel v_{{\rm th},{\rm i}})^{-1}$, comparable to the inverse of the bounce frequency (2.7). With the mode then slowly decaying exponentially, the maximum value of the pressure anisotropy may be estimated by setting $\exp (-2{\rm i}\zeta \tau )\simeq 1$ and maximizing (2.12b) with respect to $\tau$. The result is a maximum pressure anisotropy ${\simeq }2.6 \alpha \beta ^{-1}_{{\rm i}0}$ (cf. (2.8)) occurring at $k_\parallel v_{{\rm th},{\rm i}} t\simeq 2.3$. The approximate solution (2.12) is traced by the dashed line in figure 1(a), and is a manifestly good description of the full solution.

Figure 1. (a) Solution of (2.10) using the method presented in Appendix A for the time-dependent root-mean-square pressure anisotropy of a linear NP mode with wavenumber $k_\parallel$ and dimensionless initial amplitude $\alpha \doteq \delta B_\parallel (0)/B_0$ for $\beta _{{\rm i}0}=16$ and various $T_{\rm e}/T_{{\rm i}0}$. The small oscillations present after the initial adjustment are due to fast waves generated as the isothermal, pressure-balanced initial condition settles into the NP eigenmode. The approximate analytic solution (2.12) is shown with the dashed line. (b) Maximum pressure anisotropy (divided by $\alpha$) versus $T_{\rm e}/T_{{\rm i}0}$; its values at $T_{\rm e}/T_{{\rm i}0}=1/2$, $1$ and $2$ are indicated.

The second complication when using (2.8) to determine the kinetic stability of the NP mode is related to how the mode perturbs the perpendicular and parallel plasma $\beta$ parameters that feature in the firehose and mirror instability thresholds. Using (2.6) and that $\delta B_\perp = - (k_\parallel /k_\perp )\delta B_\parallel$, one obtains

(2.13a)$$\begin{gather} \beta_{{\parallel}{\rm i}} \simeq \beta_{{\rm i}0} \left( 1 + 2 \frac{\delta B_\parallel}{B_0} + \frac{k^2}{k^2_\perp} \frac{\delta B^2_\parallel}{B^2_0} \right)^{{-}1} \left[ 1 - \frac{k^2}{k^2_\perp} \frac{1}{\beta_{{\rm i}0}}\frac{T_{\rm e}}{T_{{\rm i}0}}\left(1+\frac{T_{\rm e}}{T_{{\rm i}0}}\right)^{{-}1}\frac{\delta B_\parallel}{B_0} \right] , \end{gather}$$

(2.13b)$$\begin{gather}\beta_{{\perp}{\rm i}} \simeq \beta_{{\rm i}0} \left( 1 + 2 \frac{\delta B_\parallel}{B_0} + \frac{k^2}{k^2_\perp} \frac{\delta B^2_\parallel}{B^2_0} \right)^{{-}1} \left[ 1 - \frac{k^2}{k^2_\perp}\frac{1}{\beta_{{\rm i}0}}\left(2+\frac{T_{\rm e}}{T_{{\rm i}0}}\right)\left(1+\frac{T_{\rm e}}{T_{{\rm i}0}}\right)^{{-}1} \frac{\delta B_\parallel}{B_0} \right] . \end{gather}$$

The final terms in both of these expressions may be dropped in the limit of $\beta _{{\rm i}0}\gg 1$. Combining the result with (2.8) yields

(2.14)\begin{equation} \beta_{{\parallel}{\rm i}}\Delta_{\rm NP} \approx \beta_{{\perp} {\rm i}}\Delta_{\rm NP} \approx{-}2 \frac{\delta B_\parallel}{B_0} \left( 1 + 2\frac{\delta B_\parallel}{B_0} + \frac{k^2}{k^2_\perp} \frac{\delta B^2_\parallel}{B^2_0} \right)^{{-}1} . \end{equation}

Equation (2.14) indicates that it is impossible to produce a pressure anisotropy that is sufficiently negative to destabilize the plasma to the firehose. Regions in which $\Delta _{\rm NP}<0$ also have a reduced plasma $\beta$, and so the more negative the anisotropy becomes (for larger $\delta B_\parallel >0$), the further the firehose threshold (${\approx }-2/\beta _{\parallel {\rm i}}$) moves away. Indeed, minimizing the right-hand side of (2.14) for $\delta B_\parallel > 0$, the most negative value of $\beta _{\rm i}\Delta _{\rm NP}$ is found to be ${\approx } -(1+|k/k_\perp |)^{-1} > -1/2$. In contrast, the plasma in regions where $\delta B_\parallel <0$ that acquire a positive pressure anisotropy have an easier time of reaching the reduced mirror threshold (${\approx }1/\beta _{\perp {\rm i}}$). Setting the right-hand side of (2.14) to unity and solving for $\delta B_\parallel = -|\delta B_\parallel |$ then provides the following amplitude threshold for the NP mode to trigger the mirror instability:

(2.15)\begin{equation} \left|\frac{\delta B_\parallel}{B_0}\right| \gtrsim 0.3 \quad \text{(NP mode amplitude threshold)}. \end{equation}

When this criterion is satisfied, we anticipate regions of kinetically unstable plasma to be localized to where $\delta B_\parallel <0$ and to host ion-Larmor-scale mirrors.

With these predictions borne in mind, we now determine the spatial extent of these mirror-unstable regions and discuss how the mirrors growing within them evolve to regulate the pressure anisotropy.

2.1.4. Regulation of pressure anisotropy by the mirror instability

In § 2.1.3, we showed that the plasma where $\delta B_\parallel < 0$ becomes mirror-unstable at $t_{\rm m} \sim (\sqrt {\alpha }k_\parallel v_{{\rm th},{\rm i}})^{-1}$ if initialized from isothermal pressure balance. With $\alpha \gtrsim 0.3$ (i.e. when instability is possible), this time is comparable to the time scale over which the NP mode's pressure anisotropy is set up (see figure 1). We may then view the mirror instability as growing on top of an otherwise weakly decaying positive pressure anisotropy satisfying (2.14) with $\delta B_\parallel <0$. The maximum growth rate of the instability depends on how far the local pressure anisotropy ventures beyond the instability threshold, $\varLambda _{\rm m} \doteq \Delta - \beta ^{-1}_{\perp {\rm i}} > 0$. In the asymptotic limit $\beta _{\perp {\rm i}}\varLambda _{\rm m}\ll 1$, the maximum mirror growth rate and associated wavenumber are given by (Hellinger Reference Hellinger2007; A.F.A. Bott et al., in preparation)

(2.16) \begin{equation} \gamma_{\rm m}/\varOmega_{\rm i}\approx 0.07 \beta_{{\perp} {\rm i}}\varLambda_{\rm m}^2 , \quad k_{{\parallel},{\rm m}}\rho_{\rm i} \approx 0.2 \beta_{{\perp} {\rm i}}\varLambda_{\rm m} , \quad k_{{\perp},{\rm m}}\rho_{\rm i}\approx 0.6(\beta_{{\perp} {\rm i}}\varLambda_{\rm m})^{1/2}. \end{equation}

However, because of the sensitive dependence of the instability parameter $\beta _{\perp {\rm i}}\varLambda _{\rm m}$ on the NP mode amplitude (see (2.14)), with its value ranging from ${\sim }1$ to ${\sim }100$ for $\alpha \in [0.4,0.9]$, only very marginally unstable NP modes (viz., $\alpha \approx 0.3$) satisfy the ordering used to derive (2.16). The growth rate and wavenumber when $\beta _{\perp {\rm i}}\varLambda _{\rm m} \gtrsim 1$ can be obtained by a direct numerical solution of the linearized Vlasov–Maxwell equations for a bi-Maxwellian plasma, with (2.14) specifying the pressure anisotropy for a given NP mode amplitude $\alpha$ (A.F.A. Bott, private communication). The resulting growth rates and wavenumbers are shown versus $\alpha$ in figure 2(a). (For this figure, we used $\beta _{{\rm i}0}=16$ and $k_\perp /k_\parallel =4$, although the values shown are insensitive to either parameter as long as $\beta _{{\rm i}0}\gtrsim 10$ and $(k/k_\perp )^2\approx 1$.) As $\alpha$ increases, these quantities approach the empirical values

(2.17a–c)\begin{equation} \gamma_{\rm m}/\varOmega_{\rm i} \approx 0.2\varLambda_{\rm m}, \quad k_{{\parallel},{\rm m}}\rho_{\rm i} \approx 0.6 , \quad k_{{\perp},{\rm m}} \rho_{\rm i} \approx 1.2 . \end{equation}

Figure 2. (a) Perpendicular ($k_{\perp,{\rm m}}$) and parallel ($k_{\parallel, {\rm m}}$) wavenumbers of the fastest-growing mirror mode having growth rate $\gamma _{\rm m}$, all computed from linear Vlasov–Maxwell theory using the instability parameter $\varLambda _{\rm m}$ corresponding to an NP mode with $\alpha \doteq |\delta B_\parallel /B_0|$ and $k_\perp /k_\parallel =4$ in a $\beta _{{\rm i}0}=16$ plasma (see (2.18); these values are weakly dependent upon $\beta _{{\rm i}0}$ and $k_\perp /k_\parallel$ so long as $\beta _{{\rm i}0}\gtrsim 10$ and $(k/k_\perp )^2\approx 1$). Dotted lines correspond to the asymptotic expressions (2.16), valid for $\beta _{\perp {\rm i}}\varLambda _{\rm m}\ll 1$. (b) Predicted number of mirrors $N_{\rm m}$ within the $\delta B_\parallel < 0$ region of an NP mode having wavelength $\lambda _\parallel$ and amplitude $\alpha$ (see (2.20)).

For the mirror instability to be relevant to the linear evolution of the NP mode, two criteria must be satisfied. First, the mirror growth rate must be much larger than the rate at which the NP mode decays (2.5), i.e. $\gamma _{\rm m} \sqrt {{\rm \pi} }\beta _{\rm i} \gg k_\parallel v_{{\rm th},{\rm i}}$. This condition appears to be trivially satisfied in high-$\beta$ plasmas for unstable NP modes with parallel wavelengths $\lambda _\parallel \gtrsim 10^3\rho _{\rm i}$. The second criterion is that the mirror modes must actually fit inside the length of the region that is mirror unstable, viz. $2{\rm \pi} /k_{\parallel,{\rm m}}\lesssim \ell _{\rm mirror}$. We estimate $\ell _{\rm mirror}$ by asking where in the NP mode the quantity (2.14) is larger than unity:

(2.18)\begin{equation} \varLambda_{\rm m} \simeq \frac{1}{\beta_{{\rm i}0}} \left({-}1 - 4\frac{\delta B_\parallel}{B_0} - \frac{k^2}{k^2_\perp}\frac{\delta B^2_\parallel}{B^2_0} \right) \gtrsim 0. \end{equation}

Because the leading-order eigenvector components are all real, we can take $\delta B_\parallel = -\alpha B_0 \cos (k_{\parallel }x + k_{\perp }y)$ (as used in our simulations; see § 2.2). By courtesy of our assumption that $k_\perp \gg k_\parallel$, we have that $\delta B_\perp \ll \delta B_\parallel$, so the field lines are approximately straight everywhere and the paths taken by the trapped particles as they bounce are approximately parallel to $\boldsymbol {B}_0$. Then, taking the appropriate root of (2.18) to ensure that the inverse cosine is defined for mirror-unstable amplitudes, we find that the length of the mirror-unstable portion of the wave satisfies

(2.19)\begin{equation} \ell_{\text{mirror}} \approx \frac{\lambda_{{\parallel}}}{\rm \pi}\cos^{{-}1}\left( \frac{2-\sqrt{4-k^2/k^2_\perp}}{\alpha} \right) \doteq f_{\rm m}\lambda_\parallel . \end{equation}

For $\alpha \approx 0.3$–$0.9$ and $k_\parallel \ll k_\perp$, $f_{\rm m}\approx 0.1$–$0.4$. The number of maximally growing mirrors that can fit within $\ell _{\rm mirror}$ is then

(2.20)\begin{equation} N_{\rm m} \approx \frac{f_{\rm m}}{4} \left(\frac{k_{{\parallel},{\rm m}}\rho_{\rm i}}{2{\rm \pi}}\right)\left( \frac{\lambda_\parallel}{\rho_{\rm i}}\right) . \end{equation}

In writing (2.20), we have included an additional factor of ${\approx }1/4$ to account for the fact that the pressure anisotropy is not expected to be uniform within the mirror-unstable region and so the full extent of $\ell _{\rm mirror}$ is unlikely to be filled with mirrors of identical wavelengths; the bespoke factor of ${\approx }1/4$ was obtained empirically from examining the spatial extent of scattering mirrors formed in the hybrid-kinetic simulations of unstable NP modes presented in § 2.2. A further, and final, adjustment to $N_{\rm m}$ accounts for the fact that the ion-Larmor radius $\rho _{\rm i} \propto \sqrt {T_{\perp {\rm i}}}/B$ in the mirror-unstable region is larger than $\rho _{{\rm i}0}$, primarily because of the decrease in the local magnetic-field strength. Using (2.6) to express $\delta T_{\perp {\rm i}}$ in terms of $\delta B_\parallel$, and appending a multiplicative factor of ${\approx }1/2$ to $\alpha$ to account, if only approximately, for the effective reduction in $\alpha$ due to the non-uniformity of $\delta B_\parallel$ within the mirror-unstable region, we find that

(2.21)\begin{equation} \frac{\rho_{\rm i}}{\rho_{{\rm i}0}} \approx \left( 1 - \frac{\alpha}{2\beta_{{\rm i}0}} \frac{k^2}{k^2_\perp} \right)^{1/2} \left( 1-\alpha + \frac{k^2}{4k^2_\perp} \alpha^2 \right)^{{-}1/2} . \end{equation}

With $k_{\parallel,{\rm m}}\rho _{\rm i}$ taken from figure 2(a), we can assemble (2.18)–(2.21) to predict $N_{\rm m}$ for a given $\lambda _\parallel /\rho _{{\rm i}0}$, $\alpha$ and $k_\parallel /k_\perp$ of the NP mode at $\beta _{{\rm i}0}\gg 1$.

The result of this procedure is shown in figure 2(b) as the open circles. Note that the number of mirrors $N_{\rm m}$ is fairly independent of the NP mode amplitude for $\alpha \gtrsim 0.4$, with the consequence that several mirrors can fit within the mirror-unstable region of an NP mode with $\lambda _\parallel \sim 10^3\rho _{{\rm i}0}$. However, at the critical amplitude $\alpha \approx 0.3$, only one or perhaps two mirrors are predicted to fit if $\lambda _\parallel \sim 10^3\rho _{{\rm i}0}$. In this case, the mirror instability might be ineffective at regulating the pressure anisotropy.

In summary, we predict that an NP mode with $\alpha \gtrsim 0.4$ and $\lambda _\parallel \gtrsim 10^3\rho _{{\rm i}0}$ should be able to support a robust collection of mirror-unstable fluctuations.

2.1.5. Effective collisionality induced by the mirror instability

We now seek an estimate for the effective collision frequency instigated by these mirror-unstable distortions in the magnetic-field lines. For this, we follow the arguments of Newman (Reference Newman2020) for the pitch-angle diffusion of charged particles in regions of Larmor-scale magnetic irregularities. First, we conjecture that each encounter of an ion with the edges of a single mirror depletes the plasma's temperature anisotropy $\mathcal {A} \doteq \overline { w_{\perp }^2}/2 - \overline {v_{\parallel }^2}$ by a fraction $\chi$ (here, the overline indicates an average over the ion distribution function). Following Newman (Reference Newman2020), we identify $\chi$ with $(3/2)\sin ^2\vartheta$, where $\vartheta$ is the local deflection angle of the perturbed magnetic-field lines. We estimate $\sin ^2\vartheta \approx (\delta B_{\perp,{\rm m}}/B)^2 \approx (k_{\parallel,{\rm m}}/k_{\perp,{\rm m}})^2 (\delta B_{\parallel,{\rm m}}/B)^2$, and leverage prior results on the nonlinear evolution of the mirror instability showing that mirrors can grow to amplitudes $|\delta B_{\parallel,{\rm m}}/B| \approx 1/3$ before saturating through strong pitch-angle scattering (Kunz et al. Reference Kunz, Schekochihin and Stone2014a; Riquelme et al. Reference Riquelme, Quataert and Verscharen2015; Sironi & Narayan Reference Sironi and Narayan2015). The result is that

(2.22)\begin{equation} {\chi \approx 0.2 (k_{{\parallel},{\rm m}}/k_{{\perp},{\rm m}})^2} . \end{equation}

To obtain the effective collision frequency $\nu _{\rm eff}$, we then multiply $\chi$ by the number of Larmor-scale mirrors per unit time encountered by a typical particle. For an NP mode with amplitude $\alpha \gtrsim 0.4$, the criterion for a particle to be able to pass through the NP mode's enhancement in $|B|$ is $|v_\parallel |/w_\perp \gtrsim \sqrt {4/3}$. In other words, for a near-Maxwellian distribution of particle velocities, a majority of the particles will be confined to the trough of the NP mode where ion-Larmor-scale mirrors should be present, passing through this mirror-unstable region twice per bounce time. In this case, $N_{\rm m}$ scattering mirrors are encountered by each trapped particle every transit time $\Delta t \approx {\rm \pi}\varOmega ^{-1}_{\rm b}$. The average rate of change of the ion anisotropy is then

(2.23)\begin{equation} \frac{\Delta \mathcal{A}}{\Delta t} \approx{-}\frac{\chi}{\rm \pi} N_{\rm m} \varOmega_{{{\rm b}}} \mathcal{A} \doteq{-} \nu_{\text{eff}} \mathcal{A} , \end{equation}

where in the last equality, we have introduced the effective collision frequency $\nu _{\rm eff}$. Assembling (2.7) and (2.18)–(2.23), we find that

(2.24a)\begin{equation} \nu_{\rm eff} \approx 0.003 \,\mathcal{G} \varOmega_{{\rm i}0} , \end{equation}

where

(2.24b)\begin{equation} \mathcal{G} \doteq k_{{\parallel},{\rm m}}\rho_{\rm i} \,\left(\frac{k_{{\parallel},{\rm m}}}{k_{{\perp},{\rm m}}}\right)^2 \left( \alpha-\alpha^2 + \frac{k^2}{4k^2_\perp} \alpha^3 \right)^{1/2} \cos^{{-}1}\left( \frac{2-\sqrt{4-k^2/k^2_\perp}}{\alpha} \right) \end{equation}

is a function of only the amplitude and wavenumber obliquity of the NP mode.

Equation (2.24) states that the predicted $\nu _{\rm eff}$ is independent of the wavelength of the NP mode and increases with increasing $\alpha$, key features that are tested (and confirmed) in § 2.2.6. The predicted dependence of $\nu _{\rm eff}$ upon $\alpha$ at $\beta _{{\rm i}0}=16$ and $k_\perp /k_\parallel =4$ is shown in figure 3(a); the values shown are insensitive to either parameter as long as $\beta _{{\rm i}0}\gtrsim 10$ and $(k/k_\perp )^2\approx 1$. The predicted collision frequency drops gradually from $\alpha =0.9$ to $0.5$, and then falls sharply by more than an order of magnitude to $\nu _{\rm eff}\lesssim 10^{-5}\varOmega _{{\rm i}0}$ at $\alpha =0.3$. In figure 3(b), we plot the minimum parallel wavelength $\lambda _\parallel$ of an NP mode for which $\nu _{\rm eff}\Delta t\ge 1$, where $\Delta t = {\rm \pi}\varOmega ^{-1}_{\rm b}$. Such modes should host mirrors whose scattering frequency is comparable to the transit time. Note that, for $\alpha =0.3$, $\lambda _\parallel /\rho _{{\rm i}0}$ must be ${\gtrsim }10^5$ for the scattering frequency to be larger than the inverse transit time. It is worth bearing these numbers in mind when interpreting the simulation results presented in §§ 2.2.3 and 2.2.6.

Figure 3. (a) Predicted scattering frequency $\nu _{\rm eff}$ (see (2.24)) caused by the mirror instability for an NP mode with amplitude $\alpha$, using the values of $k_{\parallel,{\rm m}}\rho _{\rm i}$ in figure 2(a). (b) Minimum parallel wavelength $\lambda _\parallel$ of an NP mode for which $\nu _{\rm eff}\Delta t \ge 1$, where $\Delta t = {\rm \pi}\varOmega ^{-1}_{\rm b}$. Such modes should host mirrors whose scattering frequency is comparable to the transit time. The data in both panels correspond to $\beta _{{\rm i}0}=16$ and $k_\perp /k_\parallel =4$, although the values shown are insensitive to either parameter as long as $\beta _{{\rm i}0}\gtrsim 10$ and $(k/k_\perp )^2\simeq 1$.

2.1.6. Suppression of nonlinear saturation of the NP mode

Once $\nu _{\mathrm {eff}}$ becomes competitive with the bounce frequency, the induced scattering will isotropize the ion distribution function faster than the nonlinear saturation can maintain the plateau in $\delta f(v_\parallel )$ around $v_\parallel \sim 0$. In this case, the nonlinear saturation is suppressed and the NP mode should resume its decay at a rate comparable to (2.5) (Johnston Reference Johnston1971). At some point during this decay, the mode amplitude will pass below its critical threshold for triggering the mirror instability (2.15), and the mirror modes themselves will become short-wavelength decaying NP modes. Near the mirror-instability threshold, these short-wavelength NP modes decay very slowly, and so the associated magnetic-field-strength fluctuations will remain nonlinear for some time after the large-scale NP mode is no longer formally mirror unstable. We therefore conjecture that the NP mode will continue to decay until the mirror fluctuations (and their induced scattering) have had sufficient time to dissipate. Excepting perhaps the case of asymptotically long NP mode wavelengths, then, there should be some delay between when the NP mode passes below threshold and when its nonlinear saturation is re-established.

The preceding arguments imply that three distinct regimes exist for collisionless NP modes in high-$\beta$ plasmas. (i) When the mode amplitude satisfies $|\delta B_\parallel /B_0| < 0.3$, the associated pressure anisotropy is too small to trigger the mirror instability, and the mode experiences slow Barnes damping until the damping nonlinearly saturates as the distribution function flattens around $v_\parallel \sim 0$. These pressure-balanced structures are thus long-lived. (ii) When $|\delta B_\parallel /B_0|\gtrsim 0.3$, the pressure anisotropy triggers the mirror instability in regions where $\delta B_\parallel < 0$ and eventually introduces an effective collisionality that, for sufficiently large NP mode wavelengths, suppresses the maintenance of a nonlinear plateau. As a result, linear decay resumes until the NP mode decays back well below its amplitude threshold. (iii) Because the induced scattering rate (2.24) does not scale with the wavelength of the NP mode, one might expect a third fluid-like regime results at very long wavelengths when $\nu _{\mathrm {eff}} \gg k_\parallel v_{\mathrm {th},{\rm i}}$ and the collisionless damping is arrested altogether. We discuss the realizability of this third regime and speculate on its behaviour in § 2.2.6.

2.2.1. Method of solution and initial conditions

To test the theory presented in § 2.1 and explore the nonlinear evolution of a mirror-infested NP mode, we employ the hybrid-kinetic particle-in-cell code Pegasus++ (Kunz, Stone & Bai Reference Kunz, Stone and Bai2014b; Arzamasskiy et al., in preparation). Pegasus++ evolves the ion distribution function $f(t,\boldsymbol {r},\boldsymbol {v})$ using a collection of positively charged macro-particles that interact with the self-consistent electromagnetic fields $\boldsymbol {E}(t,\boldsymbol {r})$ and $\boldsymbol {B}(t,\boldsymbol {r})$, which are in turn evolved on a discrete mesh using Faraday's law and a generalized Ohm's law that includes the inductive electric field, the Hall effect and a thermoelectric field caused by pressure gradients in the (assumed massless) electron fluid. The latter ensures quasi-neutrality. For simplicity, we adopt an isothermal equation of state for the electrons with temperature $T_{\rm e} = T_{{\rm i}0}$. Both the interpolation of fields to the macro-particle locations and the deposition of the macro-particles’ phase-space information on the mesh are performed using second-order-accurate triangle-shaped stencils.

All simulations of the NP mode are performed on a two-dimensional mesh that is elongated in the direction of a mean magnetic field $\boldsymbol {B}_0 = B_0\hat {\boldsymbol {x}}$ and spans one full NP mode wavelength, $L_x \times L_y = \lambda _\parallel \times \lambda _\perp$. The latter ranges from $\lambda _\parallel = 1000\rho _{{\rm i}0}$ to $4000\rho _{{\rm i}0}$, with aspect ratios of either $\lambda _{\parallel }/\lambda _{\perp }=4$ or $8$. When varying these two dimensions, the transverse dimension is never smaller than $250\rho _{{\rm i}0}$, thereby guaranteeing sufficient scale separation between the NP mode and any ion-Larmor-scale instabilities. In all runs, the spatial resolution is $\Delta x = \Delta y \simeq 0.3\rho _{{\rm i}0}$ and the number of macro-particles per cell is either $10^4$ or $5\times 10^3$ (the latter used only in our largest simulations); these values are similar to those used in previously published Pegasus simulations of collisionless Alfvén waves (Squire et al. Reference Squire, Kunz, Quataert and Schekochihin2017a) and IAWs (Kunz et al. Reference Kunz, Squire, Schekochihin and Quataert2020) in firehose/mirror-susceptible plasmas. A digital low-pass filter is applied to the computed moments of the ion distribution function to reduce the impact of grid-scale, finite-particle-number noise on the evolution of the NP mode and the trajectories of the particles.

At $t=0$, we perturb the magnetic field using the vector potential

(2.25)\begin{equation} \boldsymbol{A}(x,y) ={-}\frac{\alpha B_0}{|k|} \sin(k_\parallel x + k_\perp y) \hat{\boldsymbol{z}}, \end{equation}

where $k_\parallel = 2{\rm \pi} /\lambda _\parallel$, $k_\perp = 2{\rm \pi} /\lambda _\perp$ and $\alpha$ is a dimensionless number quantifying the mode amplitude. To excite the NP mode, the associated change in the magnetic pressure,

(2.26)\begin{equation} \frac{\delta B^2}{8{\rm \pi}} ={-} \frac{\alpha B^2_0}{8{\rm \pi}} \cos(k_\parallel x + k_\perp y ) \left[ \frac{2k_\perp}{|k|} - \alpha \cos(k_\parallel x + k_\perp y ) \right] , \end{equation}

must be exactly balanced by a perturbation to the perpendicular pressure of the plasma (cf. (2.6c)). To keep the initialization of the latter relatively simple, we choose to begin not from an exact NP eigenmode but rather from an isothermal perturbation to the plasma density $\delta n$, in which case the perturbed perpendicular pressure is simply $\delta p_\perp = \delta n (T_{{\rm i}0} + T_{\rm e})$. Balancing this expression by (2.26) and solving for $\delta n$ leads to the initial ion distribution function

(2.27)\begin{equation} f(0,x,y,v) = F_{\rm M}(v) \left\{ 1 + \frac{\alpha}{\beta_0} \cos(k_\parallel x + k_\perp y) \left[ \frac{2k_\perp}{|k|} - \alpha \cos(k_\parallel x + k_\perp y) \right] \right\} . \end{equation}

In this case, the initial total (magnetic plus thermal) pressure in the simulation domain is constant and equal to $(B^2_0/8{\rm \pi} )(1+\beta _0)$; recall that $\beta _0 \doteq \beta _{{\rm i}0}(1+T_{\rm e}/T_{{\rm i}0})$. Starting from a pressure-isotropic plasma has the advantage that any pressure anisotropy that develops is generated self-consistently and not put in by hand. It is also consistent with the assumptions made to obtain the analytic solution for $\Delta _{\rm NP}(t)$, (2.10).

In all but two of our simulations, we set $\beta _{{\rm i}0} = 16$, a value large enough to allow comparison with the asymptotic expressions derived in § 2.1, but not so large that we cannot capture a full decay time of the linear NP decay rate. We vary $\alpha \in [0.1,0.8]$, spanning the predicted NP amplitude threshold for triggering the mirror instability (2.15). Special attention is paid to the case with $\lambda _\parallel = 2000\rho _{{\rm i}0}$, $\lambda _\perp = 500\rho _{{\rm i}0}$ and $\alpha =0.8$; we refer to this as our fiducial case. Two additional runs, one with $\beta _{{\rm i}0}=4$ and the other with $\beta _{{\rm i}0}=36$, both having $\alpha =0.8$, $\lambda _\parallel =1000\rho _{{\rm i}0}$ and $\lambda _\perp =250\rho _{{\rm i}0}$, were also performed.

Hereafter, $\langle \cdot \rangle$ denotes a spatial average taken over the entire domain; $\langle \cdot \rangle _{\rm m}$ denotes a spatial average taken over the mirror-unstable region of the NP mode; and $\langle \cdot \rangle _k$ denotes a spatial average taken over the $y$-direction while accounting for the changing position of the wavefront (so as to align all of the perturbed and unperturbed regions within the domain). The latter is referred to as a ‘wavefront average’; note that it leaves the $x$-coordinate (in the direction of $\boldsymbol {B}_0$) unchanged.

2.2.2. Overall evolution of the fiducial run

We begin our discussion of the Pegasus++ results by using the fiducial run to make contact with some of the predictions laid out in § 2.1. These predictions include the excitation and subsequent linear collisionless damping of the NP mode, its nonlinear saturation, the simultaneous generation of mirror-unstable pressure anisotropy in the regions of the mode where $\delta B_\parallel < 0$, and the resumption of linear damping following the pitch-angle scattering of trapped ions by the saturated Larmor-scale mirrors at a rate larger than the bounce frequency. Figure 4 illustrates these evolutionary phases by depicting the amplitude of the NP mode versus time. After a rapid adjustment from the isothermal pressure-balanced initial condition, the NP mode emerges and decays at the linear rate (black line) for approximately one bounce time, $\varOmega ^{-1}_{\rm b}$. Immediately thereafter, the decay stalls (blue line) as nonlinear saturation sets in. Figure 5 demonstrates that, meanwhile, the NP mode has produced a large, positive pressure anisotropy in the regions where $\delta B_\parallel < 0$ and almost zero pressure anisotropy elsewhere, consistent with the prediction (2.14) (dashed line). The mirror-unstable region (with $\langle \beta _{\perp {\rm i}}\Delta \rangle _k$ above the dotted line in figure 5) is seen to occupy ${\sim }40\,\%$ of the NP mode wavelength, consistent with (2.19) for $\alpha = 0.8$. It is in this mirror-unstable region that the magnetic field acquires moderate-amplitude, oblique fluctuations in its strength on ion-Larmor scales, which are clearly apparent in figure 6. The strongest fluctuations occupy roughly a quarter of the box length and acquire amplitudes comparable to that of the mean field. The associated distortions in the field lines ultimately scatter particles at a rate comparable to the bounce frequency (see figure 7 and the accompanying discussion in § 2.2.3). As a result, the NP mode amplitude enters a ‘suppressed saturation’ phase (figure 4, red line), during which the nonlinear plateau is eroded by the mirror-induced collisionality and the Barnes damping resumes.Footnote ⁴

Figure 4. Amplitude of the magnetic-field-strength perturbation of the NP mode versus time from the fiducial run, with the different phases of the predicted evolution labelled and colour-coded. The dashed line indicates the linear decay rate (2.5) of the NP mode in a pressure-isotropic plasma with $\beta _{{\rm i}0}\gg 1$. See § 2.2.2 for discussion.

Figure 5. Wavefront-averaged profiles of $\beta _{\perp {\rm i}}\Delta$ and $\beta _{\parallel {\rm i}}\Delta$ at $k_\parallel v_{{\rm th},{\rm i}}t=3.1$, when the pressure anisotropy is near its maximum value, compared against the theoretical predictions from the linear eigenmode (2.14), for $\alpha = 0.8$ and different NP mode wavelengths $\lambda _\parallel$ and $\lambda _\perp$. The fiducial run corresponds to the solid black line. Positive values of $\beta _{\rm i}\Delta$ far exceeding the mirror threshold occur in the regions where $\delta B_\parallel < 0$. Elsewhere, negative pressure anisotropy is compensated by a decrease in $\beta _{\rm i}$ to avoid exciting the firehose instability.

Figure 6. The $x$-component of the magnetic-field perturbation, filtered to remove wavenumbers associated with the $\alpha =0.8$ NP mode, at $k_\parallel v_{{\rm th},{\rm i}}t = 25$ in the fiducial run. By this time, the mirror instability is fully nonlinear, causing large-amplitude, small-wavelength deflections in the magnetic-field direction that pitch-angle scatter particles.

Figure 7. Effective collisionality $\nu _{\rm eff}$ caused by the mirror instability in the fiducial run with $\alpha =0.8$ and $\lambda _\parallel =2000\rho _{{\rm i}0}=4\lambda _\perp$. (a) Space–time diagram of $\langle \nu _{\rm eff}\rangle _k$ (colour). An illustrative particle trajectory is shown with the grey line, exhibiting resonant bouncing, followed by trapping within a mirror fluctuation and eventual scattering out of resonance with the NP mode. (b) Box-averaged (black) and maximum wavefront-averaged (red) collision frequencies versus time.

In the remainder of § 2.2, we examine these phases in more detail and their dependence on mode amplitude and scale separation, starting with the mirror-induced scattering and its impact on the NP mode's pressure anisotropy.

2.2.3. Effective collisionality: particle scattering and trapping

Figure 7 displays the evolution of the mirror-induced effective collisionality $\nu _{\rm eff}$ in the fiducial run, calculated following the method used by Kunz et al. (Reference Kunz, Schekochihin and Stone2014a, Reference Kunz, Squire, Schekochihin and Quataert2020), Melville, Schekochihin & Kunz (Reference Melville, Schekochihin and Kunz2016) and Squire et al. (Reference Squire, Kunz, Quataert and Schekochihin2017a). Namely, the individual magnetic moments of ${\sim }10^4$ particles are tracked and monitored for (both abrupt and accumulated) changes by at least a factor of $\kappa =1.2$ (as used by Kunz et al. Reference Kunz, Squire, Schekochihin and Quataert2020 to measure firehose/mirror-induced scattering in unstable IAWs). The time intervals $\tau$ between which these changes occur are stored, along with the locations at which the changes occurred, and a spatially dependent effective collision frequency $\nu _{\rm eff}$ is calculated from the mean scattering time $\langle \tau \rangle$ using $\nu _{\rm eff} \doteq (\ln \kappa )^2/\langle \tau \rangle$. This calculation was also performed using $\kappa \in [1.1,1.5]$, with no significant differences arising.

In figure 7(b), the box-averaged effective collisionality (black line) and maximum value of the wavefront-averaged effective collisionality (red line) are shown as functions of time. Both exhibit rapid growth during the initial phase of the mirror instability and then reach a quasi-steady state, with ${\rm max}(\langle \nu _{\rm eff}\rangle _k) \approx 0.0035\varOmega _{{\rm i}0} \approx 2.5\varOmega _{\rm b}$. We have found the time scale for the scattering rate to reach this steady state to be largely independent of the wavelength of the NP mode, although it increases somewhat with decreasing $\alpha$ because of the slower linear growth rate of the mirror instability. The space–time diagram of the wavefront-averaged collisionality shown in panel (a) indicates that the maximum value of $\nu _{\rm eff}$ is localized to the centre of the mirror-unstable region, with slightly smaller values occurring near this region's boundaries where the mirror amplitudes are smaller (cf. figure 6). A large fraction of the thermal plasma is subject to this collisionality, because the mode amplitude is large enough that most of the plasma particles are confined in the regions where $\delta B_\parallel < 0$. For example, when $\alpha =0.8$, particles whose pitch angles satisfy $|v_\parallel |/w_\perp \leq \sqrt {{\rm max}(B)/{\rm min}(B)-1} \approx 2.8$ would be mirror-confined in the absence of collisions. Outside of these regions, where the plasma is stable, the collisionality is very low; as a result, the box-averaged collisionality is more than a factor of 5 smaller than the maximum value. Figure 7(a) also shows the path of a single tracked particle as a grey line. The initial evolution demonstrates bouncing within the $\delta B_\parallel <0$ region. Once the mirror fluctuations reach nonlinear amplitudes, the particle is temporarily trapped within a growing mirror. Eventually, it scatters enough in pitch angle to become de-trapped and traverses the $\delta B_\parallel > 0$ region, breaking its resonance with the NP mode.

2.2.4. Evolution of pressure anisotropy

Figure 8(a) shows the evolution of the maximum of the wavefront-averaged $\Delta$ and $\beta _{\perp {\rm i}}\Delta$ in the fiducial run. Figure 8(b) depicts the growth of the root-mean-square amplitude of the mirror fluctuations, averaged over the mirror-unstable region where $\delta B_\parallel <0$ and normalized to the average ‘background’ (i.e. guide-field plus NP-mode) magnetic-field strength in this region. The fluctuations grow large enough to scatter particles and restore the linear decay of the NP mode, through which the pressure anisotropy decays. Indeed, $\langle \Delta \rangle _k$ is similar to the linear prediction (2.12), denoted here by the blue dashed line. Likewise, $\langle \beta _{\perp {\rm i}}\Delta \rangle _k$ is modelled well by (2.14) with the substitution $\delta B_\parallel /B_0 = -\alpha \exp (-{\rm i}\zeta k_\parallel v_{{\rm th},{\rm i}}t)$, where $\zeta$ is the linear eigenvalue (2.5). This expression is traced by the dashed red line in figure 8, where we have started the decay at $k_\parallel v_{{\rm th},{\rm i}} t =9$ and set $\alpha =0.68$ to account for the delay due to the (temporary) nonlinear saturation. At larger scale separations, we anticipate that faster pitch-angle scattering induced by the mirrors will be able to regulate the pressure anisotropy more efficiently than its linear decay, at which point the mode will no longer resemble the collisionless linear NP eigenmode (see § 2.2.6).

Figure 8. (a) Maximum of the wavefront-averaged $\Delta$ (solid blue line) and $\beta _{\perp {\rm i}}\Delta$ (solid red line) versus time in the fiducial run. The evolution of $\langle \Delta \rangle _k$ matches well the predicted linear evolution (blue dashed line), suggesting that the rapid reduction of $\beta _{\perp {\rm i}}\Delta$ is due mostly to the resumed decay of the NP mode and the decrease in $\beta _{\perp {\rm i}}$ caused by the growing mirror fluctuations. (b) Root-mean-square amplitude of the mirror fluctuations, averaged over the mirror unstable region and normalized to the average ‘background’ (i.e. guide-field plus NP-mode) magnetic-field strength in the mirror region. The growth of the mirror instability coincides with a drop in $\langle \beta _{\perp {\rm i}}\Delta \rangle _k$.

The growth of mirrors leads to modifications in the shape of the NP mode profile, as shown in figure 9. The evolution of the wavefront-averaged profile of $\beta _{\perp {\rm i}}\Delta$ in the fiducial run at $k_\parallel v_{{\rm th},{\rm i}}t = 3$, 6, 11 and 27 is shown. The profile in the region where the mirror instability is active has flattened, although the mode seems to remain close to the linear eigenmode, as evidenced by figure 8. The reduction in $\beta _{\perp {\rm i}}\Delta$ occurs considerably faster than the linear decay of $\Delta$ by itself, which highlights the importance of $\beta _{\perp {\rm i}}$ in achieving marginal stability. This reinforces the idea that the mirror fluctuations do not so much act directly on the anisotropy to achieve $\beta _{\perp {\rm i}}\Delta = 1$, but rather they enable the NP mode to decay and reduce both $\Delta$ and $\beta _{\perp {\rm i}}$ to achieve marginal stability more rapidly than would otherwise occur.

Figure 9. Temporal evolution of the wavefront-averaged profile of $\beta _{\perp {\rm i}}\Delta$. Four times are shown: just after the adjustment into the NP eigenmode during the initial decay phase (black line); an intermediate time during which the NP mode decay is saturated (blue line); after the mirrors become nonlinear and scatter particles fast enough to suppress the NP mode's saturation (red line); and later once $\beta _{\perp {\rm i}}\Delta$ has been reduced enough that the mirrors are marginally stable (grey line).

2.2.5. Suppression of nonlinear saturation and resumption of transit-time damping

The effects of nonlinear saturation and mirror-induced collisionality across a variety of NP mode amplitudes can be seen in figure 10. For reasons of computational cost, for these runs, we used $\lambda _\parallel =1000\rho _{\rm i}$ rather than the fiducial $2000\rho _{\rm i}$. A Fourier transform is used to select the magnitude of the box-wavelength perturbation to the background field (i.e. the amplitude of the NP mode); this quantity is plotted as a function of time. In panel (a), the initial phase of evolution is featured, at first demonstrating linear decay at a rate similar to the prediction (2.5) (shown by a black dashed line), approximately independent of $\alpha$. After roughly one bounce time (marked by dotted lines of matching colour), the decay begins to stall and the mode amplitude tends towards a constant value. This nonlinear saturation occurs at earlier times for larger mode amplitudes, trending with the $\alpha ^{-1/2}$ scaling of the bounce time (see (2.7)). At amplitudes $\alpha \gtrsim 0.4$, more than 90 % of the original mode amplitude is preserved by the nonlinear saturation, suggesting that large-amplitude collisionless NP modes at high $\beta$ can be rather long lived.

Figure 10. Amplitude of the magnetic-field-strength perturbation of the NP mode, normalized to its initial value, versus time for $\lambda _\parallel = 1000\rho _{{\rm i}0} = 4\lambda _\perp$ and different $\alpha$. (a) Early times, during which the NP mode nonlinearly saturates after approximately one bounce time ${\sim }\varOmega ^{-1}_{\rm b}$ (vertical dotted lines; see (2.7)). The dashed line indicates the linear decay rate (2.5). (b) Late times, showing suppression of nonlinear saturation for amplitudes $\alpha \ge 0.6$.

Figure 10(b) shows the behaviour of these modes over longer time scales. For amplitudes $\alpha \le 0.4$, nonlinear saturation remains and the linear decay rate is never again realized. By contrast, the larger values of pressure anisotropy in the $\alpha =0.6$ and $0.8$ NP modes produce mirror fluctuations with amplitudes large enough to interfere with the maintenance of the nonlinear plateau. These modes are then able to decay further and convert magnetic energy into particle energy through a balance between plateau generation and pitch-angle scattering. The linear decay rate is fully re-established at $\alpha = 0.8$. A slightly weaker decay rate is seen in the $\alpha =0.6$ case because of the slower mirror-induced scattering rate relative to $k_\parallel v_{{\rm th},{\rm i}}$; at the larger scale separation of $\lambda _\parallel =2000\rho _{{\rm i}0}$ (not shown), the full linear decay rate is re-established for $\alpha =0.6$. With the value of $\lambda _\parallel$ used in these runs being twice smaller than that in the fiducial run, it is notable that the time at which near-linear decay is restored by mirror-induced scattering is the same in units of $\varOmega _{{\rm i}0}$. At scale separations much larger than those we are able to simulate currently, we therefore anticipate the nonlinear plateau to be eroded almost instantly by rapid mirror growth and its associated particle scattering.

Our final piece of evidence that the nonlinear plateau is maintained at subcritical NP mode amplitudes and eroded at supercritical amplitudes is also the most direct. In figure 11, we show ion velocity distribution functions $f(v_\parallel,w_\perp )$ measured within the $\delta B_\parallel < 0$ region, with bi-Maxwellian fits subtracted, from two runs having $\lambda _\parallel =4\lambda _\perp$ and either $\alpha =0.4$ (figure 11a,c) or $0.8$ (figure 11b,d). These distribution functions were time-averaged over two intervals of duration $4(k_\parallel v_{{\rm th},{\rm i}})^{-1}$ centred approximately $k_\parallel v_{{\rm th},{\rm i}}t = 5.4$ (figure 11a,b) and $k_\parallel v_{{\rm th},{\rm i}}t = 21$ (figure 11c,d). In the $\alpha =0.4$ run, the distribution is reduced with respect to the bi-Maxwellian at high pitch angles where the ions are well trapped, which indicates flattening in the parallel distribution approximately $v_\parallel \sim 0$. This feature persists beyond $k_\parallel v_{{\rm th},{\rm i}}t = 21$, and is the cause of the stalled decay seen in figure 10. In the $\alpha =0.8$ run, the flattening observed in the early-time distribution function is removed later on, allowing transit-time damping to resume (see figure 4). In fact, a considerable enhancement in the phase-space density exists at higher $w_\perp$ near $v_\parallel =0$; we suspect that the resumed damping allows for further betatron heating of trapped particles, leading to a small population of non-thermal particles in the $w_\perp$ tail of the distribution. Note that the width of the flattened regions at early times increases dramatically with amplitude, as is expected from the trapping criterion $|v_\parallel |/w_\perp < \sqrt {|B_{\rm max}|/|B_{\rm min}|-1}$.

Figure 11. Ion velocity distribution functions $f(v_\parallel,w_\perp )$ measured within the $\delta B_\parallel < 0$ region, with bi-Maxwellian fits subtracted, from two simulations having $\lambda _\parallel =4\lambda _\perp$ and either $\alpha =0.4$ (a,c) or $0.8$ (b,d). Panels (a,b) and (c,d) correspond to a time $k_\parallel v_{{\rm th},{\rm i}}t=5.4$ and $=21$, respectively. The colour bar has been allowed to saturate for the purpose of showing detail. Dotted lines represent isocontours of total energy, $w_\perp ^2+v_\parallel ^2={\rm const}$.

2.2.6. Dependence on scale separation

The effective collision frequency predicted by (2.24) suggests that, if the initial NP mode amplitude and wavenumber obliquity were held constant, then increasing the wavelength of the mode should have no effect on the collision frequency. This can be recast as a more illustrative relationship between the thermal crossing time and the collision frequency, $\nu _{\mathrm {eff}}/(k_\parallel v_{{\rm th},{\rm i}}) \propto \lambda _{\parallel }$. Figure 12 shows the maximum value of the box-averaged effective collision frequency normalized to $k_\parallel v_{{\rm th},{\rm i}}$ for a few different NP mode wavelengths, wavenumber obliquities and amplitudes. The measured values exhibit good agreement with the proportionality expectation at both wavenumber obliquities. This evidence implies that, at yet longer wavelengths, the collision frequency will continue to increase with respect to the transit time. Note that the measured collisionality for $\alpha =0.6$ is approximately a factor of two smaller than for $\alpha =0.8$, in qualitative agreement with the prediction featured in figure 3(a) that the scattering should decrease with decreasing NP mode amplitude. The fact that the simulated NP mode with $\alpha =0.4$ and $\lambda _\parallel =1000\rho _{{\rm i}0}$ does not have its nonlinear saturation interrupted by mirrors is also consistent with the prediction in figure 3(b). Finally, the collisionality measured in the run having $\alpha =0.8$, $\lambda _\parallel /\rho _{{\rm i}0}=1000$ and $\beta _{{\rm i}0}=36$ (red diamond) is comparable to that in the otherwise-equivalent $\beta _{{\rm i}0}=16$ run (red circle), consistent with the theoretical expectation in (2.24) that $\nu _{\rm eff}$ should be independent of $\beta _{{\rm i}0}$ for $\beta _{{\rm i}0}\gtrsim 10$.Footnote ⁵

Figure 12. Maximum value of the measured mirror-induced effective collision frequency $\nu _{{\rm eff},{\rm max}}$ versus NP mode wavelength for $\beta _{{\rm i}0}=16$ at two different wavenumber obliquities and two different initial amplitudes (an additional run having $\alpha =0.8$, $\beta _{{\rm i}0}=36$ and $\lambda _\parallel /\rho _{{\rm i}0}=1000$ is also included). The predicted scaling $\nu /(k_\parallel v_{{\rm th},{\rm i}}) \propto \lambda _{\parallel }$ is shown (dashed black line), normalized to the fiducial case (red circle at $\lambda _\parallel =2000\rho _{{\rm i}0}$).

As conjectured in § 2.1.6, the linear scaling of $\nu _{\rm eff}/(k_\parallel v_{{\rm th},{\rm i}})$ with $\lambda _\parallel$ suggests a possible fluid-like regime at sufficiently long NP-mode wavelengths. To investigate this regime, if only approximately, we examine the linear decay rate of NP modes in the presence of a constant pitch-angle scattering rate, shown in figure 13. The details of how we determined this decay rate are given in Appendix B; note that the real part of the frequency is zero for all scattering rates, i.e. the mode remains non-oscillatory. On the left-hand side of the plot, the collision frequency is small and the collisionless NP mode is recovered; on the right-hand side, the collision frequency is large and the mode becomes the MHD entropy mode. The MHD entropy mode is similar to the kinetic NP mode in that it too has no real frequency, but in the fully collisional limit, it involves only a density perturbation. For the employed values of $k_\perp /k_\parallel =4$ and $\beta _{{\rm i}0}=16$, the transition between these two regimes occurs at $\nu \approx 3 k_\parallel v_{{\rm th},{\rm i}}$. Using an asymptotic expansion at high $\beta$ and $k\simeq k_\perp$, one can show that the transitional collisionality scales approximately as $\nu \sim (3/4)\sqrt {\beta _{\rm i}} k_\parallel v_{{\rm th},{\rm i}}$. With $\nu _{\rm eff}/\varOmega _{{\rm i}0} \sim (3$–$6)\times 10^{-4}$ for $\alpha \gtrsim 0.6$ (see figures 3 and 12), we estimate that the transition to the collisional regime requires a scale separation $\lambda _\parallel /\rho _{{\rm i}0}\gtrsim 10^4\sqrt {\beta _{{\rm i}0}}$. Under this condition, the mirror-induced scattering will both isotropize the pressure perturbation and prevent resonant particles from continuously sapping energy from the wave, thereby reducing the decay rate and morphing the collisionless NP mode into the MHD entropy mode (at least for as long as the mirrors continue to scatter particles faster than ${\sim }\sqrt {\beta _{\rm i}} k_\parallel v_{{\rm th},{\rm i}}$). Unfortunately, unless the scale separation is extremely large (e.g. $\lambda _\parallel /\rho _{{\rm i}0}\gtrsim 10^5$ for our parameters), the decay rate will not be much slower than in the $\nu =0$ case. In the absence of affordable numerical simulations to test this point,Footnote ⁶ we simply conjecture that at asymptotic wavelengths, the reduction in the decay rate would allow these NP structures to become long lived once again, much like their below-threshold, nonlinearly saturated counterparts.

Figure 13. Linear decay rate of the NP mode obtained from the Landau-fluid CGL-MHD equations (B1) (see Appendix B for details). The dimensionless (complex) frequency $\zeta \doteq \omega /(|k_\parallel |v_{{\rm th},{\rm i}})$ is computed numerically as a function of collisionality $\nu /(|k_\parallel |v_{{\rm th},{\rm i}})$ for $k_\perp = 4|k_\parallel |$, $\beta _{{\rm i}0} = 16$ and $T_{\rm e} = T_{{\rm i}0}$. Overlaid are red circles marking the maximum box-averaged scattering rates measured in our hybrid-kinetic simulations (see figure 12).

3.1.1. Model equations and assumptions

Collisionless fast magnetosonic waves are in many ways simpler than their non-propagating counterparts, particularly so if their wavevectors are nearly perpendicular to the background magnetic field, viz. $k_\perp \gg k_\parallel$. In this limit, collisionless damping is extremely weak, and magnetic tension plays virtually no role in the mode's propagation. In fact, for exactly perpendicular propagation ($k_\parallel =0$), Landau and Barnes damping are entirely absent at long wavelengths due to the limited cross-field transport of magnetized particles. In this case, no kinetic information to approximate these modes other than their pressure anisotropy is needed, and they can be described entirely within double-adiabatic MHD – a model that results from taking the first three fluid moments of the drift-kinetic system (see Appendix B) and dropping the heat fluxes. Setting $\boldsymbol {B} = B\hat {\boldsymbol {y}}$ and $\boldsymbol {\nabla }=\hat {\boldsymbol {x}} \,\partial /\partial x$, these equations are

(3.1a)$$\begin{gather} \frac{{\rm D} n}{{\rm D} t} ={-}n \frac{\partial u_\perp}{\partial x}, \end{gather}$$

(3.1b)$$\begin{gather}m_{\rm i} n \frac{{\rm D} u_\perp}{{\rm D} t} ={-} \frac{\partial }{\partial x} \left( p_{{\perp} {\rm i}} + p_{\rm e} + \frac{B^2}{8{\rm \pi}} \right), \end{gather}$$

(3.1c)$$\begin{gather}\frac{{\rm D} B}{{\rm D} t} ={-}B \frac{\partial u_\perp}{\partial x} , \end{gather}$$

(3.1d)$$\begin{gather}\frac{{\rm D} }{{\rm D} t}\left(\frac{p_{{\perp} {\rm i}}}{nB}\right) = 0 , \end{gather}$$

(3.1e)$$\begin{gather}\frac{{\rm D} }{{\rm D} t}\left(\frac{p_{{\parallel}{\rm i}} B^2}{n^3}\right) = 0, \end{gather}$$

where ${\rm D}/{\rm D}t \doteq \partial / \partial t + u_\perp \partial /\partial x$. Although the right-hand side of (3.1b) is independent of the parallel pressure, and so (3.1e) is not needed to close this set of equations, (3.1e) is nevertheless useful for calculating the fast-wave pressure anisotropy. As in § 2, we adopt a simple equation of state for the electrons, $p_{\rm e} = n T_{\rm e}$, with $T_{\rm e}$ being constant.Footnote ⁷

In what follows, we investigate analytically two features of fast-wave propagation in a collisionless, magnetized plasma, adopting the simple but illustrative case of $k_\parallel = 0$. First, we demonstrate that such waves nonlinearly steepen quicker in double-adiabatic MHD than they do in standard (pressure-isotropic) MHD, a direct consequence of the proportional relationship between $T_\perp$ and $B$ associated with $\mu$ conservation, (3.1d). Second, we show how the resulting pressure anisotropy can destabilize the plasma to both firehose and mirror instabilities. We then estimate the effective scattering frequency introduced into the plasma by these instabilities and discuss how the consequent regulation of the pressure anisotropy affects the characteristics of the fast wave.

Before proceeding, it is useful to linearize (3.1) to obtain the fast-wave dispersion relation and eigenmode. Perturbing the plasma with approximately a uniform background having density $n_0$, isotropic ion pressure $p_{{\rm i}0}$ and magnetic-field strength $B_0$, we find that

(3.2a,b)\begin{equation} \frac{\delta p_{{\perp},{{\rm i}}}}{p_{{\rm i}0}} = 2\frac{\delta B}{B_0} \quad \text{and} \quad \frac{\delta p_{{\parallel}{\rm i}}}{p_{{\rm i}0}} = \frac{\delta n}{n_0} = \frac{\delta B}{B_0} . \end{equation}

These equations state that the density and pressure anisotropy are positively correlated with the magnetic-field strength, with the parallel ion temperature remaining constant. The dispersion relation of this double-adiabatic (‘da’) fast wave is

(3.3)\begin{equation} \omega = k_\perp v_{\rm A} \sqrt{1 + \beta_{{\rm i}0} \left(1+\frac{T_{\rm e}}{2T_{{\rm i}0}} \right)} \doteq k_\perp v_{\mathrm{ms},\mathrm{da}}, \end{equation}

so that the bulk velocity $u_\perp = v_{{\rm ms},{\rm da}}(\delta B/B_0)$. For comparison, the dispersion relation of a fast wave in single-adiabatic (‘sa’) MHD is

(3.4)\begin{equation} \omega = k_\perp v_{\rm A} \sqrt{1 + \beta_{{\rm i}0}\left(\frac{\varGamma}{2} +\frac{T_{\rm e}}{2T_{{\rm i}0}} \right)} \doteq k_\perp v_{\mathrm{ms},\mathrm{sa}} , \end{equation}

where $\varGamma$ is the adiabatic index of the ions. The proportional relation between the magnetic-field strength and the density in the double-adiabatic model means that $v_{{\rm ms},{\rm da}}>v_{{\rm ms},{\rm sa}}$. This increase will play a role in allowing double-adiabatic fast waves to form shocks faster than single-adiabatic fast waves, especially so at high $\beta$.

3.1.2. Wave steepening in double- versus single-adiabatic MHD

For waves in which the perturbed quantities determine the wave propagation speed, steepening may result. Large-amplitude waves in particular generate significant differences in the propagation speed between the peaks and the troughs, a situation expected to occur in both double- and single-adiabatic MHD fast waves. In this section, we perform a series of manipulations to the system (3.1) to quantify this effect. Before proceeding, it is convenient to renormalize quantities using the Alfvén speed $v_{\rm A} = B_0/(4{\rm \pi} m_{\rm i} n_0)^{1/2}$ and the wavelength $\lambda$ as follows: $u_{\perp } = \tilde {u}_{\perp } v_{\rm A}$, $B = \tilde {B}{B_0}$, $n = \tilde {n} n_0$, $x = \tilde {x} \lambda$, $t = \tilde {t} \lambda /v_{\rm A}$ and $p_{\perp,{\rm i}} = \tilde {p}_{\perp {\rm i}} m_{\rm i}n_0 v_{\rm A}^2$. We also note that if the perturbations satisfy $\delta \tilde {n} = \delta \tilde {B}$ at $t=0$, then these two quantities will remain equal for all times (see (3.1a) and (3.1c)); we can then eliminate $\delta \tilde {n}$ in favour of $\delta \tilde {B}$.Footnote ⁸ Meanwhile, if $\delta \tilde {B}$ is small and its associated perturbations in $\tilde {p}_{\perp,{\rm i}}$ and $\tilde {n}$ are given by (3.2a,b), (3.1d) becomes

(3.5)\begin{equation} \frac{\partial }{\partial \tilde{t}} \left( \frac{\tilde{p}_{{\perp},{\rm i}}}{\tilde{n} \tilde{B}} \right) \approx{-} \tilde{u}_{{\perp}} \frac{\beta_{{\rm i}0}}{2} \frac{\partial (\delta\tilde{B})^2}{\partial \tilde{x}} \sim \mathcal{O}[(\delta \tilde{B})^3]. \end{equation}

Hence, to second order in $\delta \tilde {B}$, we may treat $\tilde {p}_{\perp,{\rm i}} = (\beta _{{\rm i}0}/2)\tilde {B}^2$ as the equation of state if the initial condition is an eigenmode.

Under these conditions, (3.1) may be combined to obtain the following system:

(3.6)\begin{equation} \frac{\partial }{\partial \tilde{t}} \left[ \begin{array}{c} \tilde{u}_{{\perp}} \\ \delta \tilde{B} \end{array} \right]+ \left[ \begin{array}{cc} \tilde{u}_{{\perp}} & 1+\beta_{{\rm i}0} \left(1 + \dfrac{T_{\rm e}/2T_{{\rm i}0}}{1+\delta \tilde{B}} \right)\\ 1+\delta\tilde{B} & \tilde{u}_{{\perp}} \end{array} \right] \frac{\partial }{\partial \tilde{x}} \ \left[ \begin{array}{c} \tilde{u}_{{\perp}} \\ \delta\tilde{B} \end{array} \right] = 0. \end{equation}

Defining $\boldsymbol {W} = [\tilde {u}_\perp, \delta \tilde {B}]^\mathrm {T}$, (3.6) can be rewritten as $\partial _{\tilde {t}} \boldsymbol {W} + \boldsymbol{\mathsf{A}}(\boldsymbol {W}) \partial _{\tilde {x}} \boldsymbol {W} = 0$, with $\boldsymbol{\mathsf{A}}(\boldsymbol {W})$ being the evolution matrix. By first finding the eigenvalues $l^{(i)}$ and left eigenvectors $\boldsymbol {L}^{(i)}$ of $\boldsymbol{\mathsf{A}}(\boldsymbol {W})$, this system can be solved via its characteristic equations, which are given by $\boldsymbol {L}^{(i)} \boldsymbol {\cdot } {\rm d}\boldsymbol {W} = 0$. These characteristic equations are obeyed along space–time trajectories following ${\rm d} \tilde {x}/{\rm d} \tilde {t} = l^{(i)}$. However, because our equation of state is only valid up to second order in the wave amplitude, we need only to retain those terms of first order in the evolution matrix, and hence in its eigenvalues. Therefore, we expand the characteristic equations to first order and integrate them to find that the combinations

(3.7)\begin{equation} \eta_{{\pm}} = \tilde{u}_\perp{\pm} \tilde{v}_{\mathrm{ms},\mathrm{da}} \delta \tilde{B} \end{equation}

are approximately constant along

(3.8)\begin{equation} \left( \frac{{\rm d} \tilde{x}}{{\rm d} \tilde{t}} \right)^\pm{=} \frac{\eta_+{+} \eta_-}{2} \pm \tilde{v}_{\mathrm{ms},\mathrm{da}} \left[ 1 + \frac{1+\beta_{{\rm i}0}}{4 \tilde{v}_{\mathrm{ms},\mathrm{da}}^3}(\eta_+{-}\eta_-) \right]. \end{equation}

These can be reformulated as two nonlinear advection equations,Footnote ⁹

(3.9)\begin{equation} \frac{\partial \eta_\pm}{\partial \tilde{t}} + \left( \frac{{\rm d} \tilde{x}}{{\rm d} \tilde{t}} \right)^\pm \frac{\partial \eta_\pm}{\partial \tilde{x}} = 0. \end{equation}

Note that if the initial conditions are those of the fast eigenmode (as previously assumed in the assertion that $\delta \tilde {B} = \delta \tilde {n}$ for all time), then $\eta _-=0$ for all time. We are then left with

(3.10)\begin{equation} \frac{\partial\eta_+}{\partial\tilde{t}} + \left[ \frac{\eta_+}{2} + \tilde{v}_{\mathrm{ms},\mathrm{da}} \left( 1 + \frac{1+\beta_{{\rm i}0}}{4\tilde{v}_{\mathrm{ms},\mathrm{da}}^3}\eta_+\right)\right]\frac{\partial\eta_+}{\partial\tilde{x}} = 0, \end{equation}

the solution of which for $\tilde {u}_{\perp }$ is given by the method of characteristics as

(3.11)\begin{equation} \tilde{u}_{{\perp}}(\tilde{t},\tilde{x}) = \delta\tilde{u}_{{\perp} 0} \left(\tilde{t},\tilde{x} - \tilde{v}_{\mathrm{ms},\mathrm{da}}\tilde{t} \left[ 1 + \delta\tilde{B}_0(\tilde{x}_{\rm i}) + \frac{1+\beta_{{\rm i}0}}{2\tilde{v}_{\mathrm{ms},\mathrm{da}}^2} \delta\tilde{B}_0(\tilde{x}_{\rm i}) \right] \right), \end{equation}

where the subscript ‘0’ denotes an initial value and $x_{\rm i}$ is the $x$-position of the source of a given characteristic. The time-dependent solution for an example large-amplitude, double-adiabatic fast wave is shown in figure 14. This solution is strictly valid only until a shock has formed, at a time that may be determined by evaluating the eigenvalue $l^+$ at the location $x_0$ where its derivative achieves its largest negative value:

(3.12)\begin{equation} t^{\rm da}_{\rm s} = \left[l^+(x_0)\right]^{{-}1} \approx \left[\alpha k_\perp v_{{\rm ms},{\rm da}} \left(1 + \frac{v_{\rm A}^2}{v_{\mathrm{ms},\mathrm{da}}^2} \frac{1+\beta_{{\rm i}0}}{2} \right)\right]^{{-}1} , \end{equation}

where $\alpha \doteq \delta \tilde {B}(0)$ is the initial fast-wave amplitude. This double-adiabatic (‘da’) shock-formation time is to be compared with the corresponding time in a single-adiabatic MHD plasma, in which $pn^{-\varGamma }=p_0n^{-\varGamma }_0$. The general problem of fast-wave steepening in MHD plasmas has been studied thoroughly under many conditions (Hada & Kennel Reference Hada and Kennel1985; Ödblom Reference Ödblom1998; Sujith Reference Sujith2005). Following an analogous process to that used for the double-adiabatic fast wave, we find the single-adiabatic (‘sa’) shock-formation time:

(3.13)\begin{equation} t_{\rm s}^{\rm sa} \approx \left[\alpha k_\perp v_{\mathrm{ms},\mathrm{sa}} \left(1 + \frac{v_{\rm A}^2}{v_{\mathrm{ms},\mathrm{sa}}^2} \frac{1 +\varGamma(\varGamma-1)\beta_{{\rm i}0}/2}{2} \right)\right]^{{-}1}. \end{equation}

Simplifying (3.12) and (3.13) at high $\beta$, and setting $T_{\rm e} = T_{{\rm i}0}$ and $\varGamma = 5/3$, yields

(3.14a,b)\begin{equation} k_\perp v_{\rm A}t_{\rm s}^\mathrm{da} \approx \frac{\sqrt{6}}{4\alpha \sqrt{\beta_{{\rm i}0}}} \quad \text{and} \quad k_\perp v_{\rm A}t_{\rm s}^\mathrm{sa} \approx \frac{12\sqrt{3}}{29\alpha \sqrt{\beta_{{\rm i}0}}} \simeq 1.17 k_\perp v_{\rm A}t^{\rm da}_{\rm s}. \end{equation}

The single-adiabatic shock-formation time is thus larger than the double-adiabatic shock-formation time. When $T_{\rm e}/T_{{\rm i}0} = 0$, their ratio reaches a maximum of ${\simeq }1.23$; for $T_{\rm e} \gg T_{{\rm i}0}$, it approaches unity. This increase is a consequence of the direct correlation between the magnetic-field strength and the perpendicular (ion) pressure in double-adiabatic MHD, which amplifies local changes in the mode propagation speed.

Figure 14. Approximate solution (3.11) to the fast-wave steepening problem with initial amplitude $\alpha = 0.3$ and $\beta _{{\rm i}0}=25$. The solution has just begun to form a shock, indicating a shock-formation time of $k_\perp v_{\rm A} t_{\rm s} \sim 0.4$.

3.1.3. Pressure anisotropy and its regulation by kinetic instabilities

By contrast with the NP mode, the fast wave generates a fluctuating pressure anisotropy as the wave propagates. At sufficiently large $\beta$, both firehose and mirror instabilities may therefore be triggered. With $\delta p_{\perp,{\rm i}}$ and $\delta p_{\parallel,\mathrm {i}}$ given by (3.2a,b), the amplitude threshold for triggering both firehose and mirror instabilities is

(3.15)\begin{equation} \left|\frac{\delta B}{B_0}\right| \gtrsim \frac{2}{\beta_{\rm i}} \quad\text{(fast-wave amplitude threshold)} . \end{equation}

At high $\beta$, this criterion can be satisfied for even small-amplitude fluctuations, justifying the use of the linear eigenvector and unperturbed $\beta _{\rm i}$ in determining the threshold.

To assess whether these micro-instabilities will be able to grow, we compare their linear growth rates with the linear frequency of the fast wave at high $\beta$, $\omega _{\rm fast}\sim k_\perp v_{{\rm th},{\rm i}}$. We adopt the maximal mirror growth rate from (2.16), and use the maximal oblique firehose growth rate $\gamma _{{\rm f}} \approx 0.3\varOmega _{\rm i}\varLambda ^{1/2}_{{\rm f}}$, where $\varLambda _{{\rm f}} \doteq |\Delta + 2/\beta _{\rm i}|$ (Yoon et al. Reference Yoon, Wu and de Assis1993; A.F.A. Bott et al., in preparation), both of which are appropriate for the near-threshold conditions we anticipate in our fast-wave simulations. Assuming $|\delta B/B_0|\gtrsim 2\beta ^{-1}_{\rm i}$, we find that

(3.16a,b)\begin{equation} \frac{\gamma_{\rm m}}{\omega_{\rm fast}} \sim 0.01 \beta^{{-}1}_{\rm i}\frac{\lambda_\perp}{\rho_{\rm i}} \quad{\rm and}\quad \frac{\gamma_{{\rm f}}}{\omega_{\rm fast}} \sim 0.1 \beta^{{-}1/2}_{\rm i} \frac{\lambda_\perp}{\rho_{\rm i}} , \end{equation}

where $\lambda _\perp = 2{\rm \pi} /k_\perp$ is the wavelength of the fast wave. It is immediately apparent from (3.16a,b) that, at high $\beta$, very large scale separation between the fast-wave wavelength and the ion-Larmor scale is necessary to allow enough time for mirror fluctuations to grow and become nonlinear. The scaling with $\beta _{\rm i}$ is much weaker for the firehose instability, and so there will exist wavelengths at which mirror regulation of the pressure anisotropy is effectively non-existent but the firehose regulation is rapid. For this reason our Pegasus++ simulations, which focus on $\beta _{{\rm i}0}=25$, require $\lambda _\perp \gg 10^3 \rho _{{\rm i}0}$ to realize both mirror and firehose regulation.

The unstable Larmor-scale fluctuations will ultimately grow to amplitudes at which the particles’ rate of pitch-angle scattering is sufficient to hold the pressure anisotropy at marginal stability. This rate may be estimated by calculating the pressure anisotropy driven by a small-amplitude fast wave in a weakly collisional plasma (following Braginskii Reference Braginskii1965) and asking what value of effective collisionality $\nu _{\rm eff}$ would be required to keep $|\Delta |\sim 2\beta ^{-1}_{\rm i}$. With the former given in the collisional regime by $\Delta \sim -(\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u})/\nu _{\rm eff} \sim (k_\perp v_{\rm ms}/\nu _{\rm eff}) |\delta B/B_0|$, the limiting collisionality is

(3.17)\begin{equation} \nu_{\rm eff} \sim k_\perp v_{\rm ms} \frac{\beta_{\rm i}}{2} \left|\frac{\delta B}{B_0}\right| . \end{equation}

Note its explicit dependence upon the scale of the fast wave, an indirect consequence of the pressure anisotropy of the fast wave being continuously driven by the fluctuating wave. This is very different from the case with the aperiodic NP mode, in which the pressure anisotropy – an essential feature of the mode's perpendicular pressure balance – actually decays in time through transit-time damping.

3.1.4. Viscous damping and collisional propagation

The estimate of the effective collisionality (3.17) suggests that, depending on the wave amplitude, one should see a variety of fast-wave behaviour. For example, if $|\delta B/B_0| \gg 2\beta _{{\rm i}0}^{-1}$, then the implied collisionality can be large enough to push the fast wave into the collisional Braginskii-MHD regime ($\nu \gg \omega$). Making the presently unjustified yet instructive assumption that this collisionality is distributed uniformly in space, the fast-wave dispersion relation at arbitrary $\nu$ can be obtained after including isotropizing collisional terms $-\nu \Delta p/nB$ and $\nu \Delta p B^2/n^3$ on the right-hand sides of (3.1d) and (3.1e), respectively, then linearizing the resulting system of equations. We find that

(3.18)\begin{equation} \omega^3 - i\nu\omega^2 - \omega k^2_\perp v_{{\rm ms},{\rm da}}^{2} +i\nu k^2_\perp v_{{\rm ms},{\rm sa}}^{2} = 0. \end{equation}

The numerical solution to (3.18) is shown in figure 15. In the collisionless limit $\nu \rightarrow 0$, one recovers propagation at the double-adiabatic fast speed; taking $\nu \rightarrow \infty$ returns propagation at the single-adiabatic fast speed. Viscous damping occurs at intermediate values of $\nu \sim \mathrm {Re}(\omega ) \sim k_\perp v_{{\rm th},{\rm i}}$ around the transition between the double- and single-adiabatic regimes, where the scattering rate is comparable to the wave's oscillation frequency. The damping rate is always small compared with the wave frequency. The dispersion relation (3.18) alongside the amplitude threshold (3.15) and the predicted effective collision frequency (3.17) imply three regimes for the behaviour of perpendicularly propagating fast modes in a high-$\beta$ plasma. For small amplitudes satisfying $|\delta B/B_0| < 2\beta _{{\rm i}0}^{-1}$, the mode propagates normally as a collisionless fast mode. It will steepen and eventually form a shock on the double-adiabatic shock time $t_{\rm s}^{\rm da}$. In the near-threshold regime where $|\delta B/B_0| \gtrsim 2\beta _{{\rm i}0}^{-1}$, the scattering rate from triggered mirror and firehose instabilities will not quite reach the value (3.17), though scattering is still expected to occur and result in some viscous damping. The wave will also steepen to form a shock, but only a fraction of the wavelength will be kinetically unstable and therefore the shock will occur on a hybrid of the double- and single-adiabatic shock times. Lastly, at amplitudes well above the threshold, the scattering rate should be given by (3.17). The viscous damping will be very weak, the wave will host firehose/mirror scattering sites throughout most of its wavelength and the shock time should be better represented by the single-adiabatic model.

Figure 15. Exact solution to the dispersion relation (3.18) for a $k_\parallel =0$ fast wave in a plasma having collision frequency $\nu$, $\beta _{{\rm i}0} = 25$ and $T_{\rm e}/T_{{\rm i}0}=1$.

Figure 16. Shock-formation time versus $\beta _{{\rm i}0}$ and $\alpha$ for a double-adiabatic fast wave computed from CGL-MHD simulations (lines) and predicted analytically using (3.12) (circles). The simulated waves are estimated to have formed a shock at the time when the rate of change of the maximum density gradient drops below half of its own peak value.

We now test these ideas using numerical simulations.

3.2.1. Method of solution and initial conditions

Due to the large scale separations needed to obtain asymptotic $\nu _{\mathrm {eff}}$ for both firehose and mirror fluctuations (§ 3.1.3), we use a combination of Pegasus++ and (much cheaper) Landau-fluid CGL-MHD simulations. All simulations initialize a $k_\parallel =0$ fast wave in an otherwise Maxwellian plasma using the collisionless eigenmode (3.2a,b), viz.,

(3.19)\begin{equation} \left. \begin{gathered} \boldsymbol{B}(0,x) = B_0\left[ 1 + \alpha\sin(k_\perp x) \right] \hat{\boldsymbol{y}} , \quad \boldsymbol{u}(0,x) = v_{{\rm ms},{\rm da}}\alpha \sin(k_\perp x) \hat{\boldsymbol{y}} , \\ \frac{n(0,x)}{n_0} = \frac{p_{{\parallel}{\rm i}}(0,x)}{p_{{\rm i}0}} = 1 + \alpha\sin(k_\perp x) , \quad \frac{p_{{\perp} {\rm i}}(0,x)}{p_{{\rm i}0}} = \left[ 1 + \alpha\sin(k_\perp x) \right]^2, \end{gathered} \right\} \end{equation}

where $k_\perp = 2{\rm \pi} /\lambda _\perp$ and $\alpha$ is a dimensionless number quantifying the mode amplitude. For the Pegasus++ runs, the mesh is two-dimensional and elongated in the propagation direction, with size $L_x\times L_y = \lambda _\perp \times 100\rho _{{\rm i}0}$. The size of the domain in the $y$ direction is large enough to capture all relevant firehose and mirror fluctuations. We set $\beta _{{\rm i}0}=25$ and $T_{\rm e}=T_{{\rm i}0}$; the slightly larger value of $\beta _{{\rm i}0}$, as compared with that used in the simulations of the NP mode ($\beta _{{\rm i}0}=16$), results in a shorter numerical integration time (and thus computational savings) without changing the physical character of the fast wave. The spatial resolution and the number of macro-particles per cell are the same as in the NP simulations (§ 2.2.1). In the manuscript, we only show results from a Pegasus++ run having $\lambda _\perp = 8000\rho _{{\rm i}0}$, corresponding to the largest domain size that we simulated. We found that this value of $\lambda _\perp /\rho _{{\rm i}0}$ was the minimum required for the mirrors to have time to grow and begin scattering particles before the wave oscillates and the sign of the driven pressure anisotropy reverses.

In the accompanying Landau-fluid simulations, the full system of CGL-MHD equations is solved using a new Riemann solver implemented in a version of the finite-volume Athena++ simulation code (Stone et al. Reference Stone, Gardiner, Teuben, Hawley and Simon2008) that includes Landau-fluid heat fluxes (J. Squire et al., in preparation). These equations are given in Appendix B; they reduce to (3.1) in our chosen geometry. For these runs, $\beta _{{\rm i}0}$ is varied between 1 and 100 to study the variance of the shock time. A ‘limiter’ collisionality $\nu _\textrm {lim}$ is set either to $0$ or to $\alpha \beta _{{\rm i}0} k_\perp v_{\textrm {ms},\textrm {da}}$, depending on whether the focus is on wave steepening and shock formation ($\nu =0$) or the effects of the instability-induced scattering. This anomalous scattering rate is active only within regions of the domain where the pressure anisotropy would be kinetically unstable, viz., where $\beta _{\rm i}\Delta \leq -2$ and $\beta _{\rm i}\Delta \geq 1$; elsewhere, it is zero. It serves to isotropize the plasma pressure where mirror or firehose fluctuations would otherwise do so in a kinetic system, by contributing a term proportional to $-\nu _\textrm {lim}\Delta p$ to the right-hand sides of the evolution equations for $p_\perp$ and $p_\parallel$.

As in § 2, $\langle \cdot \rangle$ denotes a spatial average taken over the entire domain, while $\langle \cdot \rangle _k$ denotes a spatial average performed along the wavefront (in this case, the $y$ direction).

3.2.2. Wave steepening and shock formation

Our first goal is to test the expression (3.12) for the shock-formation time $t_{\rm s}$. We perform a parameter survey by varying $\beta _{{\rm i}0}$ and the wave amplitude $\alpha$ using the CGL-MHD code with the micro-instability-limiting scattering turned off. At each time step in the simulation, the local density gradient (using a four-cell average) is calculated throughout the domain and its maximum value is recorded as a measure of the wave steepening. As a fast wave steepens, the growth rate of this maximum gradient increases until eventually the shock forms and the maximum gradient in the domain begins to plateau. We define the numerically calculated shock-formation time to be the time at which the rate of change of this maximum gradient drops below half of its own peak value. The resulting times are compared with (3.12) in figure 14. When testing the dependence on $\beta _{{\rm i}0}$ (blue, left), the perturbation amplitude is set to $\alpha =0.01$; when testing the dependence on amplitude (red, right), $\beta _{{\rm i}0} = 25$.

Overall, the agreement between (3.12) and the numerically calculated shock-formation times is quite good. Small variations occur due to differences in the rates at which the maximum gradients plateau and to minute fluctuations in the maximum value of the gradient after the shock is formed (this value does not necessarily reach a perfect steady state). Perhaps unsurprisingly, at high $\beta$ where $v_{\mathrm {ms},\mathrm {da}} \approx v_{\rm A}\sqrt {3\beta _{{\rm i}0}/2}$, the ratio of the wave-crossing time and the shock-formation time is $t_{\mathrm {cross}}/t_{{\rm s},\mathrm {da}} \approx 4\alpha /3$. This means the number of wavelengths propagated prior to forming a shock is dependent upon the mode amplitude only.

3.2.3. Generation of pressure anisotropy and triggering of kinetic instabilities

Prior to shock formation, the linearized fluctuations (3.2a,b) suggest that pressure anisotropy at a level capable of triggering both mirror and firehose instabilities will exist when the fast-wave amplitude satisfies $|\delta B/B_0| \gtrsim 2/\beta _{\rm i}$. For these supercritical amplitudes, the wavefront should carry with it rapidly growing firehose fluctuations and more slowly growing mirror fluctuations, as per (3.16a,b). To test this idea, we performed a large-scale Pegasus++ simulation, the parameters of which are described in § 3.2.1; the initial wave amplitude $\alpha = 0.1$ and $\beta _{{\rm i}0}=25$.

Figure 17(a) depicts the pressure anisotropy generated by the fast wave as it propagates through space at three different times ($k_\perp v_{\rm A}t = 0.0$, 0.08, 0.39; note that the aspect ratio of the plotted domain is far from unity, and that the mean magnetic field is in the $y$ direction). Initially, the positive and negative pressure anisotropies in the wave are equal in magnitude. Shortly thereafter, the (unstable) negative anisotropy is reduced significantly due to the rapid growth of the (primarily oblique) firehose instability. The positive pressure anisotropy does not show a comparable decrease, and in fact increases somewhat from its initial value. This is likely because the rapid change in the negative-anisotropy regions, which perturbs the wave and causes some deviation from the eigenmode, is not matched by a comparable regulation from the positive side because of the relatively slow mirror growth. Figure 17(b) zooms in on the corresponding magnetic-field fluctuations that emerge in two separate co-moving regions where the plasma is mirror unstable (left) or firehose unstable (right). To accentuate these fluctuations, the large-scale contribution from the fast wave has been removed. At $k_\perp v_{\rm A} t = 0.08$, oblique firehose fluctuations are strong and nonlinear; parallel firehose fluctuations are also present, though subdominant, in $\delta B_x$ (not shown). At this time, there is only a hint of mirror fluctuations emerging above the noise level. In the final frame ($k_\perp v_{\rm A} t = 0.39$) however, highly oblique mirror modes have grown to large amplitudes in the region encompassed approximately by $x/\rho _{{\rm i}0} \in [4000,5000]$. The scale separation achieved in this simulation ($L_x/\rho _{{\rm i}0}=8000$) was the minimum at which we could observe mirror fluctuations with strengths comparable to their firehose counterparts; increasing the scale separation further would come at considerable computational expense.

Figure 17. (a) Pressure anisotropy times the ion beta from a Pegasus++ simulation of a collisionless fast wave, showing that the compression and rarefaction of the magnetic-field lines generate oppositely signed anisotropies that move with the wavefront. Some sloshing due to firehose regulation of the negative pressure anisotropy causes an additional reversal of $\Delta$ in the final time frame. (b) Zoomed-in regions showing $\delta B_y$ and $\delta B_z$, with the contribution from the background fast wave removed. Recall that the mean field is oriented in the $y$ direction. In the left set of panels, the mirror instability, with its oblique orientation and dominance in $\delta B_\parallel =\delta B_y$, grows relatively slowly in the co-moving region of fast-wave compression from $k_\perp v_{\rm A}t = 0.08$ to $0.39$. The firehose instability in the right set of panels is predominantly oblique and exhibits rapid growth and saturation; smaller-amplitude parallel firehoses appear in $\delta B_x$ (not shown). These firehose fluctuations reside downstream of the mirrors, where the fast-wave $\delta B<0$.

3.2.4. Effective collisionality: particle scattering

Following § 2.2.3, the effective collisionality was determined for the fast wave shown in figure 17 by tracking thousands of ion macro-particles and measuring the frequency at which their $\mu$ changes statistically by a factor of 1.2 or more. Figure 18 depicts this scattering rate as a function of the position along the wave ($x/\rho _{{\rm i}0}$) and the time ($k_\perp v_{\rm A}t$). Sites of strong scattering are associated with the firehose modes, which appear more or less instantly and travel along with the trough of the wave. The trail of the scattering sites indicates that the trough of the wave moves at ${\approx }6v_{\rm A}$, as expected for a fast mode with $\beta _{{\rm i}0} = 25$. In this simulation, the rapid regulation of the pressure anisotropy by the firehose instability causes sloshing. The sloshing temporarily drives a higher positive pressure anisotropy, and therefore enhanced mirror growth, for a short period beginning at $kv_{\rm A}t \approx 0.4$. The measured scattering rate in the firehose-unstable regions is comparable to the predicted asymptotic scattering rate for a $\beta _{{\rm i}0}=25$ fast wave with $\alpha =0.1$ and $T_{\rm e}=T_{{\rm i}0}$, viz. $\nu _{\mathrm {eff}} \approx 16 k_\perp v_{\rm A}$ (see (3.17)). The mirror instability in this case also scatters particles at an average rate of a few times $k_\perp v_{\rm A}$, but these scattering sites are much less coherent and do not coincide with the peak in the positive pressure anisotropy. This delayed growth is a result of the limited achievable scale separation in our simulations, which only barely allows mirrors to grow to nonlinear levels within a fast-wave crossing time.

Figure 18. Space–time diagram of the effective collision frequency measured in a Pegasus++ fast wave. The simulation parameters are $\beta _{{\rm i}0} = 25$, $\alpha = 0.1$ and $T_{\rm e}/T_{{\rm i}0}=1$; using these numbers in (3.17) predicts $\nu _{\mathrm {eff}} \approx 16k_\perp v_{\rm A}$.

Figure 19. Wavefront-averaged $\beta _{\rm i}\Delta$ in the fast wave for the same time frames as figure 17. Pressure-anisotropy regulation from the firehose instability maintains $\beta _{\rm i}\Delta \gtrsim -1.4$, while the mirror fluctuations cause some distortion of the mode above $\beta _{\rm i}\Delta \approx 1$ but are unable to regulate fully the positive anisotropy to marginally unstable values. An increase in the rate at which positive pressure anisotropy is generated by the steepened wave and the asymmetry in the anisotropy's regulation by micro-instabilities causes an enhancement of the positive pressure anisotropy in the final time shown.

The effects of the induced scattering on the fast wave's pressure anisotropy are visible in figure 19, which shows $\langle \beta _{\rm i}\Delta \rangle _y$ at the same times as in figure 17. The negative anisotropy is regulated within a very short time by the firehose instability to a value close to the oblique threshold $\beta _{\rm i}\Delta \simeq -1.4$. This regulation persists, but is not matched on the mirror-unstable side. Some steepening has also occurred, as expected, but the positive anisotropy has not been driven down near marginal mirror stability. For mirror fluctuations to regulate the positive pressure anisotropy to marginal stability, they would need to grow faster with respect to the fast-wave crossing time; (3.16a,b) suggests that this could be achieved by increasing $\lambda _\perp /\rho _{{\rm i}0}$ even further (beyond $\lambda _\perp /\rho _{{\rm i}0} = 10^4$), or perhaps by decreasing $\beta _{{\rm i}0}$ (though in this case, the amplitude threshold (3.15) would increase, necessitating larger fast-wave amplitudes that would shock almost immediately). Unfortunately, such large scale separations become prohibitively expensive to simulate using Pegasus++, and so from this point onwards, we employ the CGL-MHD code with pressure-anisotropy limiters.Footnote ¹⁰

3.2.5. Viscous damping and collisional steepening

To study fast-wave behaviour at asymptotically large scale separations, we employ the Landau-fluid CGL-MHD code. These simulations are performed using a larger $\beta$ parameter than that used in the Pegasus++ run, $\beta _{{\rm i}0} = 100$ rather than 25, and with $\alpha = 0.2$. These parameters have the advantage that a large portion of the fast wave is initially above the threshold for instability while the wave remains somewhat linear in amplitude. As discussed in § 3.2.1, this code introduces a user-specified constant scattering rate in (and only in) the kinetically unstable regions of the plasma. We set this scattering rate according to (3.17) using the initial mode amplitude. In reality, this scattering rate should decay alongside the amplitude, and so our treatment will not precisely reproduce the results that would be obtained from a more rigorous kinetic calculation.

In figure 20, the propagation and nonlinear steepening of the CGL-MHD fast wave are presented. The top panel in figure 20(a) shows the bulk fluid velocity perpendicular to the background field at three different times, exhibiting steepening without a significant change in wave amplitude. This indicates that no significant viscous dissipation occurs on a time scale comparable to the shock-formation time scale (as predicted by (3.18)). The bottom panel shows the pressure anisotropy of the wave at the same times, multiplied by $\beta _{\rm i}$. The anisotropy is substantially reduced below what it would be in the absence of the limiting collisionality, particularly on the firehose-unstable side, although it is not perfectly regulated to the instability thresholds. In particular, a peak in the positive pressure anisotropy becomes prominent starting from $k_\perp v_{\rm A} t \approx 0.1$. This is a result of wave steepening, as the sharp gradient at the wavefront generates positive $\Delta$ much faster than the slow decline in the wake generates negative $\Delta$, as well as faster than our (constant) limiting collisionality is able to regulate. Figure 20(b) displays the evolution of the maximum absolute value of the density gradient from this run, alongside that from a comparable run with $\nu _\textrm {lim}=0$. On the abscissa is the simulation time normalized by the double-adiabatic shock time $t_{\rm s}^\textrm {da}$ (see (3.12)). We calculated the shock time for each run using the same detection method as in figure 16; these times, marked by filled circles in the figure, agree reasonably well with the predicted values of $t_{\rm s}^\textrm {da}$ and $t_{\rm s}^\textrm {sa}$ for the collisionless and collisional cases, respectively. The difference in steepening rate between the two runs can be interpreted as $\nu _\textrm {lim}$ forcing a more MHD-like, rather than collisionless, evolution in the fast wave. The collisional isotropization at the peaks of the wave (which are also the most rapidly moving regions) effectively changes the local adiabatic index of the ions, slowing down the steepening process and yielding better agreement with $t_{\rm s}^\textrm {sa}$ than with $t_{\rm s}^\textrm {da}$. In this sense then, all of the essential characteristics of large-amplitude, high-$\beta$, collisionless fast waves approach that of single-adiabatic MHD as a result of induced micro-instabilities.

Figure 20. (a) Propagation of an $\alpha =0.2$ fast wave with $\beta _{{\rm i}0}=100$ and $\nu _\textrm {lim}$ set by (3.17). The top panel shows wave steepening in the fluid velocity, with no noticeable viscous decay on the time scale of shock formation. The bottom panel shows regulation of the pressure anisotropy to near the mirror and firehose thresholds. A peak appears in $\beta _{\rm i}\Delta$ due to the rapid generation of positive pressure anisotropy in the steepening wavefront. (b) Maximum density gradient found within the domain of the same $\alpha = 0.2$, $\beta _{\rm i}=100$ fast wave, compared against an equivalent run with $\nu _\textrm {lim}=0$. The predicted shock times are labelled by $t_{\rm s}^\textrm {da}$ and $t_{\rm s}^\textrm {sa}$, and the shock times detected by the same method used for figure 14 are denoted by circular markers. The growth of the maximum gradient continues for a longer time in the single-adiabatic case than in the double-adiabatic case, indicating delayed shock formation.