Amplitude limits and nonlinear damping of shear-Alfv\'en waves in high-beta low-collisionality plasmas

This work, which extends Squire et al. [ApJL, 830 L25 (2016)], explores the effect of self-generated pressure anisotropy on linearly polarized shear-Alfv\'en fluctuations in low-collisionality plasmas. Such anisotropies lead to stringent limits on the amplitude of magnetic perturbations in high-beta plasmas, above which a fluctuation can destabilize itself through the parallel firehose instability. This causes the wave frequency to approach zero,"interrupting"the wave and stopping its oscillation. These effects are explored in detail in the collisionless and weakly collisional"Braginskii"regime, for both standing and traveling waves. The focus is on simplified models in one dimension, on scales much larger than the ion gyroradius. The effect has interesting implications for the physics of magnetized turbulence in the high-beta conditions that are prevalent in many astrophysical plasmas.


Introduction
In this paper, we derive and discuss stringent nonlinear limits on the amplitude of shear-Alfvén (SA) fluctuations in weakly collisional plasmas. The result, which was first presented in , is that collisionless linearly polarized SA waves -that is, low-frequency incompressible oscillations of magnetic field (δB ⊥ ) and velocity (u ⊥ ) perpendicular to a background field B 0 -cannot oscillate when where β ≡ 8πp 0 /B 2 is the ratio of thermal to magnetic pressure. Above this limit (or a related limit (2) in the weakly collisional regime), standing-wave fluctuations are "interrupted" before even a quarter oscillation, while traveling waves are heavily nonlinearly damped. In both cases, the magnetic field rapidly forms a sequence of zigzags -piecewise straight field line segments with zero magnetic tension -and evolves at later times with the magnetic energy far in excess of the kinetic energy (i.e., effectively in a near-force-free state).
What is the cause of such dramatic nonlinear behavior, even in regimes (δB ⊥ /B 0 1 for β 1) where linear physics might appear to be applicable? As we now explain, the effect depends on the development of pressure anisotropy -i.e., a pressure tensor that differs in the directions perpendicular and parallel to the magnetic field. In a magnetized plasma in which the ion gyro-frequency Ω i is much larger than the collision frequency ν c , a decreasing (in time) magnetic field leads to a decreasing pressure perpendicular to the magnetic field (p ⊥ ), while the parallel pressure (p ) increases. Such behavior originates in part from conservation of the particle's first magnetic moment µ = mv 2 ⊥ /2B, which suggests that p ⊥ /B should be conserved as B changes in a collisionless plasma. This anisotropy, ∆p ≡ p ⊥ − p < 0, provides an additional stress in the momentum equation that can neutralize the restoring effects of magnetic tension, even destabilizing the SA wave and triggering the parallel firehose instability if ∆p < −B 2 /4π (Rosenbluth, 1956;Chandrasekhar et al., 1958;Parker, 1958;Schekochihin et al., 2010).
Consider the ensuing dynamics if we start with ∆p = 0, but with a field that, in the process of decreasing due to the Lorentz force, generates a pressure anisotropy that would be sufficient to destabilize the wave. This is a nonlinear effect not captured in linear models of SA waves. As ∆p approaches the firehose limit, the magnetic tension disappears and the Alfvén frequency approaches zero, "interrupting" the development of the wave. Because the wave perturbs the field magnitude by δB 2 ⊥ , an amplitude δB ⊥ /B 0 β −1/2 is sufficient to generate such a ∆p in a collisionless plasma. As the field decrease is interrupted at the firehose stability boundary, the plasma self-organizes to prevent further changes in field strength, leading to the nullification of the Lorentz force through the development of piecewise-straight (and therefore, tension-less) fieldline structures. In addition, as this process proceeds, there is a net transfer of the mechanical energy of the wave to the plasma thermal energy due to "pressure-anisotropy heating," which occurs because of spatial correlations between the wave's self-generated pressure anisotropy and dB/dt.
A similar effect also occurs in the weakly collisional "Braginskii" regime (Braginskii, 1965). Here, collisions act to balance the anisotropy generation and SA waves cannot where ω A is the wave frequency and ν c the ion collision frequency (with ω A ν c required for the Braginskii equations to be valid). In addition, because a changing magnetic field is required to balance the collisional relaxation of ∆p, an "interrupted" wave slowly decays in time until its amplitude is below the limit (2), at which point it can oscillate. Although the details of the nonlinear dynamics differ from the collisionless regime, the dynamics in both regimes share some generic features, in particular the strong dominance of magnetic energy over kinetic energy after a wave is interrupted.
This myriad of applications has in turn led to intense study of the SA wave's basic properties across many plasma regimes (Cramer, 2011). The most relevant to our study here are several papers noting that linearly polarized SA waves are Landau damped nonlinearly at the rate ∼ ω A β 1/2 (δB ⊥ /B 0 ) 2 at high β (Hollweg, 1971b;Lee & Völk, 1973;Stoneham, 1981;Flå et al., 1989), although this rate is reduced by particle trapping effects at high wave amplitudes (Kulsrud, 1978;Cesarsky & Kulsrud, 1981;Völk & Cesarsky, 1982). This Landau damping has a similar form to the collisionless "pressure-anisotropy damping" that plays a key role in some of the effects described in this work. There have also been a wide variety of studies considering nonlinear effects due to parametric instabilities and compressibility (e.g., Galeev & Oraevskii, 1963;Hollweg, 1971a;Goldstein, 1978;Derby, 1978;Medvedev & Diamond, 1996;Medvedev et al., 1997;Del Zanna et al., 2001;Matteini et al., 2010;Tenerani & Velli, 2013), which have generally found large-amplitude SA waves to be unstable to parametric decay at low β, but with stability increasing as β approaches ∼ 1 (Bruno & Carbone, 2013). Our study here complements these previous works by showing that in the limit β 1, linearly polarized finite-amplitude SA waves in weakly collisional plasmas can be nonlinearly modified so strongly that they are unable to oscillate at all. Note, however, that circularly polarized SA fluctuations are unmodified by these effects because their magnetic field strength is constant in time.
The role of shear-Alfvén waves in magnetized turbulence deserves special emphasis: turbulence is fundamental to many areas of astrophysics and geophysics and may be significantly modified by the nonlinear amplitude limit.
The well-accepted phenomenology of Goldreich & Sridhar (1995) suggests that strong magnetized turbulence should be understood in terms of nonlinear interactions between SA wave packets, which cascade in such a way that their linear physics is of comparable importance to their nonlinear interactions (this is known as "critical balance"). Because of the resilience of SA waves to kinetic physics, it is often assumed -and patently true in some cases, e.g., the solar wind at β 1 -that Alfvénic cascades survive in collisionless plasmas (Schekochihin et al., 2009) even though naive estimates suggest the plasma viscosity is very large. † The nonlinear interruption of Alfvénic fluctuations above the amplitude δB ⊥ /B 0 ∼ β −1/2 may thus significantly alter our understanding of turbulence in weakly collisional plasmas at high β -conditions that occur, for example, in regions of the solar wind (Bale et al., 2009;Bruno & Carbone, 2013), the intracluster medium (ICM) ‡ (Rosin et al., 2011;Zhuravleva et al., 2014), and hot astrophysical disks (Balbus & Hawley, 1998;Quataert, 2001). The picture described above and in what follows suggests a limit on the amplitude (in comparison to a background field) of such turbulence, above which motions are quickly damped, leaving longer-lived magnetic perturbations in their wake. This paper, which extends the results of , is organized as follows. In Sec. 2, we present the Landau-fluid model (Snyder et al., 1997) used throughout this work to analyze nonlinear SA wave dynamics. This model is chosen as the simplest extension of MHD to weakly collisional plasmas with motions on scales that are large compared to the ion gyroradius. Given the model's relative simplicity in comparison to full Vlasov-Maxwell equations, particular focus is given to gaining qualitative understanding of various physical effects: the pressure anisotropy, collisions, and heat fluxes. Section 3 then contains a very brief description and definition of the two main physical effects -termed interruption and nonlinear damping -that form the basis for our results. We then treat Braginskii (Ω i ν c |∇u|) and collisionless (Ω i |∇u| ν c ) SA wave dynamics in Secs. 4 and 5, respectively. Because standing waves are primarily affected by the interruption effect, whereas traveling-wave dynamics are more naturally thought of in terms of nonlinear damping, we split each of these sections and separately discuss standing and traveling waves in each case. For all cases, we derive various scalings, amplitude limits, and damping rates, and describe the physics qualitatively with the aid of numerical examples. In Sec. 6, we discuss the importance of kinetic physics that is not included in our model, both due to the limitations of a † As recently argued by Verscharen et al. (2016) for the solar wind, large-amplitude compressive fluctuations may also play an important role in high-β turbulence, aiding in the isotropization of the distribution function. ‡ This is the hot plasma that fills the space between galaxies in clusters.
1-D domain of the Landau-fluid prescription for the heat fluxes. These considerations underscore the importance of future two-and three-dimensional kinetic simulations for further study of the effect. For an impatient reader, the summary of key results in Sec. 7 should be (mostly) comprehensible without reference to the main text.
Finally, the appendices deserve some mention here, being somewhat separate in character and content than the main text of the paper. In these, we derive the nonlinear wave equations asymptotically, both in the collisionless limit (Appendix A; we also consider the zero-heat-flux double-adiabatic equations there), and in the Braginskii regime (Appendix B). These calculations serve two main purposes. The first is to justify more formally many of the approximations in the main text. In this capacity, they may help comfort a reader who is skeptical of our arguments relating, e.g., to heat fluxes in collisionless waves. The second purpose is to derive explicitly various effects that are only heuristically derived in the main text, e.g., the damping rate for traveling waves. These calculations also provide a useful reference point for future fully kinetic studies that could account more formally for various effects not included in the Landau-fluid model.

Macroscopic equations for a weakly collisional plasma
Throughout this work, our philosophy is to consider the simplest modifications to macroscopic † plasma dynamics due to kinetic physics. We thus consider a two-species, fully ionized plasma, and assume that the pressure tensor is gyrotropic -i.e., invariant under rotations about the field lines -but can develop an anisotropy, viz., a different pressure parallel and perpendicular to the magnetic field lines. This approximation is generally valid for motions on spatiotemporal scales much larger than those relating to ion gyromotion. It leads to the following equations for the magnetic field and the first three moments of the plasma distribution function (Chew et al., 1956;Kulsrud, 1983;Schekochihin et al., 2010): Here Gauss units are used, u and B are the ion flow velocity and magnetic field, B ≡ |B| andb = B/B denote the field strength and direction, ρ is the mass density, ν c is the ion collision frequency, p ⊥ and p are the components of the pressure tensor perpendicular and parallel to the magnetic field, and q ⊥ and q are fluxes of perpendicular and parallel heat in the direction parallel to the magnetic field. Note that p ⊥ and p in Eq. (4) are summed over both particle species, while ρ and u in Eqs.
(3)-(5) are the ion density and flow velocity (although for kρ i 1 and m e /m i 1, they may equivalently be viewed as the total density and flow velocity). The pressure equations (6) and (7) should in principle be solved separately for each species; however, in this work we consider only the ion pressures, an approximation that may be formally justified by an expansion in the electron-ion mass ratio when the electrons are moderately collisional (see, e.g., Appendix A of Rosin et al., 2011). The double-dot notation used in Eqs. (6)-(7) meanŝ bb : ∇u ≡b ibj ∇ i u j =b · (b ·∇u). Note that nonideal corrections to the magnetic-field evolution, which are important for motions at scales approaching ρ i , are not included in Eq. (5) and will be ignored throughout this work (see, e.g., Schekochihin et al., 2010). We also define ∆ ≡ ∆p/p 0 with p 0 = 2p ⊥ /3 + p /3 (note that ∆p p 0 for β = 8πp 0 /B 2 1), the Alfvén speed v A = B 0 / √ 4πρ (with B 0 a constant background field), the sound speed c s = p 0 /ρ, parallel sound speed c s = p /ρ, and denote the ion gyroradius and gyrofrequency ρ i and Ω i , respectively. Although Eqs. (3)-(7) are derived directly from the Vlasov equation assuming kρ i 1 and ω/Ω i 1 (where k and ω are characteristic wavenumbers and frequencies of the system), the heat fluxes q ⊥, remain unspecified and must be solved for using some closure scheme (or the full kinetic equation) as discussed below.
2.1. The importance of pressure anisotropy at high β In a changing magnetic field, the termŝ (where d/dt is the Lagrangian derivative) in Eqs. (6) and (7) locally force a pressure anisotropy ∆ = ∆p/p 0 . Importantly, because this anisotropy generation depends onb rather than B, its dynamical influence increases as β increases (aside from the limiting effects of microsinstabilities; see below). Namely, the final term of Eq. (4) shows that ∆p has a strong dynamical influence (i.e., is comparable to the Lorentz force) when ∆p ∼ B 2 ; i.e., when ∆ ∼ β −1 . For β > 1, the pressure anisotropy generated by changing B will generally cause a stress that is stronger than the Lorentz force. It is also worth noting the importance of the spatial form of ∆p, which, as we shall show, can strongly influence the nonlinear dynamics. As will become clear below (Secs. 2.1.1 and 2.1.2), this spatial variation in ∆p depends on the balance between the drivingbb : ∇u and the other terms in Eqs. (6) and (7) (e.g., the heat fluxes or collisionality), so we should expect nonlinear wave dynamics to depend significantly on a particular physical regime.
In this work, we focus on two such regimes for the evolution of ∆p, neglecting compressibility for simplicity in both cases [this neglect is valid at β 1, δB ⊥ /B 0 1; see Appendix A.3 around Eq. (A.34) and Appendix B around Eq. (B.7)]. The first approximation is Braginskii MHD, which is valid in weakly collisional plasmas when Ω i ν c |∇u|; the second is collisionless (ν c = 0, or equivalently Ω i |∇u| ν c ), which we model using a simple Landau fluid (LF) closure for the heat flux.
2.1.1. Braginskii MHD. When collisions dominate (|∇u| ν c ), we may neglect ∂ t p ⊥ and ∂ t p in comparison to ν c ∆p in Eqs. (6) and (7). For β 1, these approximations also imply ∆p p 0 , leading to We have neglected q ⊥, for simplicity in deriving Eq. (9), although this is only valid in the limit δp ⊥, /p ⊥, |u|/c s (where δp ⊥, denotes the spatial variation in p ⊥, ; see Mikhailovskii & Tsypin, 1971;Rosin et al., 2011). † An expression for ∆p with heat fluxes included is derived in Appendix B [Eqs. (B.11) and (B.12)], where we also briefly discuss how the nonlinear SA wave dynamics are modified by the resulting different spatial form of ∆p. However, given the extra complexity of including this effect, we ignore the heat fluxes in the discussion of Braginskii dynamics in Sec. 4.
2.1.2. Collisionless plasma. The evolution of ∆ is strongly influenced by heat fluxes when ν c |∇u| and β 1. As a simple prescription, we employ a Landau fluid (LF) closure (Snyder et al., 1997;Hammett & Perkins, 1990;Hammett et al., 1992;Passot et al., 2012), which has been extensively used in the fusion community, and to a lesser degree for astrophysical applications (Sharma et al., 2006;Sharma et al., 2007). The heat fluxes are chosen to reproduce linear Landau damping rates, namely, where ∇ is the parallel gradient operator, while the parallel wavenumber |k | must be considered as an operator. In the regime of interest, ∆p p 0 and ν c = 0, with small perturbations to the magnetic field, the dynamical effect of q ⊥, can be easily understood. Equations (10) and (11) are These, combined withb · ∇q ⊥, q ⊥, ∇ ·b (valid for small perturbations to the background field), imply that the heat-flux contributions to the pressure equations (6) † For ∆p p 0 , this condition is approximately equivalent to ν c ∼ c s /λ mfp k c s (where λ mfp is the ion mean-free path). and (7) simplify to These terms, which model the Landau damping of temperature perturbations, † suppress spatial variation in p ⊥, over the particle crossing time τ damp ∼ (|k | c s ) −1 . This damping implies that if τ damp |∇u| −1 ∼ ω A , the k = 0 part of ∆ is suppressed by a factor of ∼ v A /c s ∼ β −1/2 compared to its mean. ‡ This leads us to the simple interpretation that the heat fluxes spatially average ∆p, by damping k = 0 components of the pressure perturbations, giving where · denotes the spatial average. The spatial form of the O(β −1/2 )(x) term generally follows the spatial variation ofbb : ∇u, and is calculated by asymptotic expansion in various regimes in Appendix A.3 and Appendix A.4 [see Eqs. (A.43) and (A.59); these calculations also justify more formally the spatial-averaging action of the heat fluxes derived heuristically above].

Microinstabilities
An important limitation of Eqs.
(3)-(7), which exists for both the Braginskii and LF closures, is their inability to capture correctly certain plasma microinstabilities. For our purposes, at high β, the most important of these are the firehose and mirror instabilities. Both of these grow fastest on scales approaching the Larmor radius, which are explicitly outside the validity of Eqs.
(3)-(7). Assuming ∆p p 0 , the firehose is unstable if ∆ < −2/β and comes in two flavors: the parallel firehose, which is present in fluid models and is the cause of SA wave interruption, and the oblique firehose (Yoon et al., 1993;Hellinger & Trávníček, 2008), which grows fastest at k ⊥ = 0, and is not correctly captured by Eqs. (3)-(7). The mirror instability is unstable if ∆ > 1/β and grows with k ⊥ k . Although the linear mirror instability is contained in the LF model (Snyder et al., 1997), its nonlinear evolution, which involves trapped particle dynamics (Rincon et al., 2015), presumably requires a fully kinetic model. It is worth noting that 1-D fully kinetic simulations would also not correctly include either the oblique firehose or mirror instabilities.
It has been common in previous literature (e.g., Sharma et al., 2006;Kunz et al., 2012;Santos-Lima et al., 2014) to model the effect of these instabilities in † Since the spatial variation in ρ will generally be similar to that of p ⊥, , the effective damping is less than what it would be if the variation in ρ were ignored in Eq. (13). However, because the spatial variation in T ⊥, is of the same order as that of p ⊥, , a damping of the pressure alone −c s |k |p ⊥ may be used for heuristic estimates. A full asymptotic calculation of the relative contributions of ρ and p ⊥, fluid simulations by applying "hard-wall" boundaries on ∆, limiting its value by the appropriate microinstability threshold. This is motivated by the fact that both in kinetic simulations and, it appears, in the observed solar-wind, microinstabilities act to limit the pressure anisotropy at its marginal values (see, for example, Hellinger et al., 2006;Bale et al., 2009;Kunz et al., 2014;Servidio et al., 2015). In addition, the enormous scale separation between the micro-and macroscales in many astrophysical plasmas implies that the effect of microscale instabilities on large-scale dynamics should be effectively instantaneous (Melville et al., 2016). Motivated by the fact that the parallel firehose instability is contained in fluid models, most of the numerical results in this work will (where appropriate) apply a limit on positive anisotropies (to model the action of the mirror instability), but not on negative anisotropies.
A more thorough discussion of microinstabities is given in Sec. 6, focusing in particular on the implications of previous kinetic results for SA wave dynamics and the possible changes that might result from a multi-dimensional fully kinetic treatment.

Energy conservation
Energy conservation arguments are used heavily throughout the paper, forming the basis for our estimates of traveling-wave damping rates in Secs. 4.2 and 5.2. With the kinetic, magnetic, and thermal energies defined as Eqs.
(3)-(7) conserve the total energy: A key difference compared to standard MHD arises in the evolution equation for the mechanical energy E mech = E K + E M : The final term in this equation describes the transfer of mechanical to thermal energy due to the presence of a spatial correlation between ∆p and B −1 dB/dt. In the (highβ) Braginskii limit, where ∆p ∝ B −1 dB/dt [see Eq. (9)], this term is always positive and represents a parallel viscous heating. In the collisionless case, it can in principle have either sign, although we shall see that for SA waves, there is a positive correlation between ∆p and B −1 dB/dt that leads to net damping of the waves. For later reference, the mean pressure anisotropy evolves according to

Shear-Alfvén wave dynamics
It is helpful to derive a simple wave equation that isolates the key features of linearly polarized shear-Alfvén waves and the influence of the pressure anisotropy. Although here the derivation is heuristic, with the aim of highlighting the key features of high-β SA dynamics, similar equations are derived asymptotically from the full LF system (3)- Our geometry is that of a background field B 0ẑ , with perturbations perpendicular toẑ and the wavevector k = k zẑ + k ⊥ . Since SA waves are unmodified by k ⊥ = 0 (the envelope is simply modulated in the perpendicular direction) and we analyze only linear polarizations, we assume x-directed perturbations that depend only on z and t, viz., Note that circularly polarized fluctuations are unaffected by the pressure-anisotropic physics because the field strength remains constant in time. Combining Eqs. (4) and (5) and neglecting compressibility, the field perturbation δb = δB ⊥ /B 0 satisfies where ∆ is given by Eq. (14) (collisionless closure) or Eq. (9) (Braginskii closure). In the absence of a background pressure anisotropy, Eq. (20) illustrates that linear longwavelength SA fluctuations are unmodified by kinetic effects. Similarly, fixing ∆ and linearizing in δb, the parallel firehose instability emerges because the coefficient of ∂ 2 z δb is negative for β∆/2 < −1.
In the following sections, we shall treat standing and traveling waves separately. While these differ only in their initial conditions, they can display rather different nonlinear dynamics. In the context of Eq. (20), a standing wave has initial conditions in either δb or u x /v A , viz., or a traveling wave involves initial conditions in both δb and u x /v A , viz., (for a wave traveling from left to right).

Numerical method
For all numerical examples, both of the Landau fluid equations and of various reduced equations (in the Appendices), we use a simple Fourier pseudospectral numerical method on a periodic domain. Standard 3/2 dealiasing is used, along with a k 6 hyperviscous diffusion operator in all variables, which is tuned so as to damp fluctuations at scales just above the grid scale. This is necessary with Fourier methods because there is little energy dissipation otherwise, and the energy can be spuriously reflected back from highk into lower-k modes. The only further approximation used in solving Eqs.
(3)-(11) is the identification of |k | in Eqs. (10)-(11) with |k z | (the |k | operator is nondiagonal in both Fourier and real space and thus somewhat expensive to evaluate). While this approximation is truly valid only for δb 1, various tests have shown that the exact form of the heat fluxes makes little difference; for example, the method of Sharma et al. (2006), which sets |k | = k L with k L a parameter, does not qualitatively modify the solutions presented here. Results shown in the figures throughout the text were obtained at a resolution N z = 512, but we see little modification of results at higher or lower resolutions.

Wave interruption and damping through pressure anisotropy
In this section, we explain the two key mechanisms that can lead to strongly nonlinear behavior of SA waves in high-β regimes. These are: (1) the nullification of the wave's restoring force (the Lorentz force) through the self-generated pressure anisotropy, which we term interruption, and (2) the channeling of wave energy into thermal energy due to spatial correlation of ∆ and dB/dt, which we term nonlinear damping.

Interruption
It is immediately clear from Eq. (20) that any time ∆(z) approaches −2/β, the solutions to Eq. (20) are fundamentally altered because the restoring force of the SA wave disappears (i.e., the coefficient of ∂ 2 z approaches zero). We term this effect "wave interruption," because the oscillation halts when this occurs.
In the Braginskii limit, with ∆ given by Eq. (9), the wave is thus interrupted when Here B = B 0 (1 + δb 2 ) 1/2 , which depends only on the current value of the field and its rate of change. By contrast, in the collisionless limit, with ∆ given by Eq. (14), the interruption occurs if a wave evolves so that 3 ln which is interesting for its explicit dependence on the initial conditions. Our derivations of amplitude limits in the following sections are simply applications of Eqs. (24) and (25).

Nonlinear damping
In the presence of a positive correlation between ∆p and B −1 dB/dt, the mechanical energy of the wave is converted to thermal energy at the rate [see Eq. (17)] We term this effect "nonlinear damping" because the fact that the B perturbation is proportional to δb 2 in a SA wave implies that the damping rate also scales with δb 2 .
In the collisionless case, there is no fundamental requirement that ∆p and B −1 dB/dt have a positive spatial correlation. Nonetheless, given that ∆p is driven by B −1 dB/dt [see Eqs. (6) and (7)], one might intuitively expect such a correlation for SA waves, and the calculations in Sec. 5.2 and Appendix A.3 show that this is indeed the case. Note, however, that its numerical value, and thus the wave damping rate, depends on the effect of the heat fluxes in smoothing ∆p [this is the O(β −1/2 )(x) part in Eq. (14)]. In the collisionless limit, there is also Landau damping of a nonlinear SA wave due to the spatiotemporal variation of the magnetic pressure (Hollweg, 1971b;Lee & Völk, 1973), which turns out to cause wave damping at a rate similar to the pressure anisotropy damping (neglecting particle trapping effects; Kulsrud, 1978).

Braginskii MHD -the weakly collisional regime
In this section we work out the behavior of SA waves in the Braginskii limit. As discussed in Sec. 2, Braginskii dynamics differ significantly from fully collisionless dynamics because the pressure anisotropy is determined by the current value of ∂ t B, rather than the time history of the magnetic field.

Standing waves
Starting from a finite-amplitude magnetic perturbation, a Braginskii standing SA fluctuation will be significantly modified (interrupted) if Eq. (24) is satisfied at some point during its decay. If we consider an unmodified standing wave with δb 0 1, then  (20) with the Braginskii closure (9), starting from a sinusoidal magnetic perturbation [initial conditions (21)], with amplitude δb 0 and some chosen ω A β/ν c (see Fig. 2). A red square indicates that an initial perturbation was interrupted before a half cycle (as in Fig. 2), while a blue circle indicates that the perturbation flipped polarity without interruption. The dashed line is δb 0 = 2.5(ω A β/ν c ) −1/2 . Note that in the incompressible limit, SA wave dynamics are determined entirely by δb 0 and the ratio ω A β/ν c , because the ν −1 c factor in ∆ in Eq. (9) multiplies β/2 in Eq. (20).
The wave will thus be significantly modified -i.e., interrupted -if ∆ −2/β at some point in space, which occurs if Above this limit, ∆p can remove the restoring force of the wave in regions where δb(z) = 0 (i.e., around the antinodes of the wave). As we show in Fig. 1, the limit (30) is well matched by numerical solutions. Figure 2 illustrates the dynamics of a standing wave above the limit (30). We solve the SA wave equation (20), using the Braginskii closure (9), which assumes only 1-D dynamics and incompressibility of the wave. † Note that within the incompressible limit, the dynamics are entirely determined by δb 0 and the ratio Although the nonlinear wave dynamics shown in Fig. 2 may appear quite bizarre, with angular field structures and sharp discontinuities in ∆p, many features can be straightforwardly understood by noting that if the field is to decrease significantly more slowly than in a linear SA wave, it must nullify the magnetic tension. This can be † Note that for the condition (30) to be met at the same time as the condition ν c ω A , required for the validity of the Braginskii equations, the system must be at very high β. Further, since u ∼ v A ∼ β −1/2 c s , the motions are very subsonic. It thus makes sense to assume incompressibility when studying Braginskii waves, and if one wished to study the lower-β, larger-δb 0 , limit between the collisionless and Braginskii regimes, it would be most sensible to solve the full LF equations (3)-(7) with the collisional relaxation terms included. More discussion, including the effects of heat fluxes, is given in Appendix B.1. Figure 2. Evolution of an initial magnetic perturbation δb = −0.5 cos(2πz) within the Braginskii model at ν c /ω A β −1/2 = 0.2. Panel (a) shows δb at t = 0 (black dotted line), δb at t = 0.6τ A (blue solid line), u ⊥ /v A at t = 0.6τ A (red solid line), and δb at t = 2τ A (black dashed line), which is after the amplitude has decreased below the interruption limit. Panel (b) illustrates the shape of the magnetic field lines in space at t = 0.6τ A (blue lines), with the shading showing where 4π∆p/B 2 = 0 (white) or −1 (gray). Panel (c) shows the anisotropy parameter, 4π∆p/B 2 , which is −1 at the firehose limit, at the same times as in panel (a) (the black dotted line shows t = 0.01τ A to illustrate the initial evolution). Note the velocity is much smaller than the magnetic field during the decay, and the lack of magnetic tension everywhere in the decaying wave.
achieved: (1) by keeping a pressure anisotropy at the firehose limit, which occurs in the δb(z) "humps" where the field is curved; or (2) by having straight field lines, which occurs where δb(z) is zero. Then, because the field must keep decreasing in order to maintain ∆ = −2/β (since ∆ ∼ ν −1 c dB/dt), it slowly decays in time, with the regions where the field is small reaching δb(z) = 0 first. Once the amplitude of the wave decays below the level at which it can sustain ∆ = −2/β (i.e., when the field at the wave antinode is δb 2 ∼ ν c /ω A β −1 ), the wave can oscillate freely again with an amplitude below the interruption limit (although δb is not sinusoidal because the final stages of the interrupted decay are nonsinusoidal). Note that throughout this decay process, the perturbation's magnetic energy dominates over the kinetic energy. This is because the pressure anisotropy stress cancels out the Lorentz force in the momentum equation, leading to a magnetic field that changes more slowly than that in a similar-amplitude Figure 3. Time t decay for an initial magnetic perturbation of the form δb(z) = δb 0 cos(k z z) to decay to an amplitude that is small enough that it can oscillate. Solid lines and symbols show the numerically measured t decay (normalized by ω −1 A ), while dashed lines of matching color show the theoretical prediction (33). The bright red circles are those points for which δb 0 is below the limit (30), meaning the wave is able to oscillate. The match with the theoretical prediction for t decay is surprisingly accurate, illustrating the usefulness of the simple arguments outlined in Sec. 4.1.
linear SA wave.
We can use these ideas to calculate the decay time of the field, t decay , as a function of β and the initial amplitude δb 0 . The idea is simply to ignore the spatial dependence of the solution, focusing on the antinode of the wave, where δb is maximal. The condition ∆ = −2/β is then which has the solution By solving for δb 2 = 0, we arrive at a prediction for the time for the interrupted field to decay to an amplitude at which it can oscillate: As shown in Fig. 3, this estimate agrees very well with numerically computed decay times (taken from calculations like that in Fig. 2) even quantitatively, illustrating the effectiveness of the simple dynamical model proposed above.
4.1.1. Standing waves with an initial velocity perturbation. It is worth briefly describing also the dynamics that one observes after initializing with a velocity rather than a magnetic perturbation. † Because in such a situation the field initially grows rather than decays, the ensuing dynamics depend on what occurs at positive pressure anisotropies, when ∆ > 1/β. Specifically, one expects growing mirror fluctuations (which are not captured in 1-D models) to act to limit ∆ at 1/β, and that this limiting action will be fast compared to ω A , so long as there is significant scale separation with the gyroscale (see Sec. 2.2; Kunz et al., 2014;Melville et al., 2016). If this limit on ∆ does not exist -i.e., if ∆ can grow without bound as the magnetic field grows in the standing wave -the extra magnetic tension arising from ∆ > 0 acts to reverse the fluctuation of the wave at low magnetic-field amplitudes, while strong nonlinear damping causes the wave to damp to an amplitude below the interruption limit (30) in less than a wave period (see Sec. 4.2). However, if the anisotropy is limited at positive values, this allows the field to grow to much higher amplitudes, viz., δb ≈ u ⊥0 /v A (for u ⊥0 /v A β −1/2 and β 1). When the field then starts decreasing again (once u ⊥ ≈ 0) it does so from an amplitude that is above the interruption limit, and thus behaves in effectively the same way as an initial purely magnetic perturbation. Thus, the dynamics -as long as growing mirror fluctuations act to limit positive pressure anisotropies -are similar to those for an initial static magnetic perturbation, and the limit on the amplitude of an initial velocity perturbation u ⊥0 /v A is similar to Eq. (30). In Sec. 6, we give a more detailed discussion of this physics. †

Traveling waves
With the Braginskii model, because ∆ ≈ ν −1 c B −1 dB/dt, the anisotropy and the rate of change of the magnetic field B −1 dB/dt are always strongly correlated. From Eq. (17), this implies the wave energy E wave = E mech = E K + E M is nonlinearly damped at the rate A sinusoidal traveling SA wave, B = B 0 (δb sin(k z z − ω A t), 0, 1), creates a changing magnetic field 1 B which, from Eq. (34), causes the wave to damp at the rate This is effectively a parallel viscous damping, which occurs because there is a component of u in the field-parallel direction due to the finite amplitude of the wave. Noting that E wave = (8π) −1 δb 2 B 2 0 , we conclude that Eq. (36) implies a wave damping rate 1 2E wave (37) † The reader may be puzzled that we are not applying such a "∆-limit" argument also to the negative anisotropies. The salient point is that ∆ can only be limited by firehose microinstabilities once ∆ < −2/β, by which point the magnetic tension has already been removed. More extensive discussion of this and related issues is given in Sec. 6.
In the top panels, blue lines illustrate δb and red lines show u ⊥ /v A , while the bottom panels show the anisotropy parameter 4π∆p/B 2 (this is −1 at the firehose limit). The left panels (a)-(b) show a wave with δb 0 = 0.5 -i.e., above the interruption limit -at t = 0 (dashed lines; we show t = 0.05 instead for the anisotropy parameter to illustrate the initial evolution) and t = 2τ A (solid lines) [the initial condition in panel (a) is sinusoidal; a reduced plot range is used to better see the later times]. The right panels (c)-(d) show a wave with δb 0 = 0.05 -i.e., below the interruption limit -at t = 0 [dashed lines; t = 0.05 in panel (b)], t = 0.8τ A (solid lines) (the earlier time is shown so as to fit the full evolution into one panel). In panel (c) hollow circles mark the same position on the wave front as it propagates. Above the interruption limit, there is such a strong nonlinear modification of the wave that it is effectively interrupted and stopped before it can propagate. In contrast, below the limit, the wave undergoes minor shape changes and slow nonlinear decay.
where δb max is the interruption limit, given by Eq. (30). For amplitudes above the interruption limit (30), Eq. (37) implies a damping rate greater than the frequency of the wave itself. In this case, the damping will cause such a strong nonlinear modification of the wave that it might be considered more accurately as an interruption, effectively stopping the wave. Indeed, the local pressure anisotropy will reach the firehose limit in regions where B −1 dB/dt < 0 (and the mirror limit where B −1 dB/dt > 0), so we should expect some of the arguments of Sec. 4.1 to apply here. A traveling wave (with the same parameters as the standing wave in Fig. 2) is illustrated in Figs. 4(a)-(b), showing how the wave is virtually stopped with a larger magnetic than kinetic energy, and an anisotropy that is similar to the standing wave [cf. Fig. 2(c)]. Thus, there is effectively an "interruption" of the same kind as for a standing wave.
Consider now a traveling wave that is well below the interruption limit, illustrated in Fig. 4(c)-(d). Such a wave exhibits much slower damping and moderate nonlinear modification to the wave shape, which should be expected because the mechanism causing the wave damping has nonlinear spatial variation in space. The more angular  (37). The agreement is reasonably good, although at larger decay rates, where the nonlinear effects are stronger, the wave reduces its damping rate by becoming more square.
structures that the wave develops act to reduce B −1 dB/dt over much of the wave, and thus reduce the damping rate somewhat.
The damping rate is measured quantitatively in the right-hand panel of Fig. 4, where the theoretical prediction (37) is compared to the rates measured in simulations. While our prediction agrees very well in the low-decay-rate limit, where the nonlinear modifications to the wave shape are small, it deviates as the wave damping increases and the pressure anisotropy causes more significant changes to the shape of the wave.
Finally, note that although, for consistency with the upcoming analysis of collisionless plasmas (Sec. 5), we have discussed standing waves and traveling waves separately, this distinction is less important for Braginskii dynamics. Indeed, we have seen that a Braginskii traveling wave above the amplitude limit (30) is effectively interrupted, because the anisotropy is so strong that it stops the wave. Analogously, a standing wave below the interruption limit will oscillate but will be nonlinearly damped at the rate (37), because there is still parallel variation in ∇u that is damped by the Braginskii viscosity. In the next section, we shall see that there is a stronger distinction between standing and traveling waves (and between interruption and nonlinear damping) for collisionless wave dynamics, because of the smoothing effect of the heat fluxes.

Collisionless waves
We now consider the dynamics of SA waves in a collisionless plasma at high β. Collisionless dynamics differ significantly from the Braginskii limit discussed in the previous section because, in a collisionless plasma, the pressure anisotropy remembers the time history of B 2 , rather than being set by the instantaneous value of dB/dt. This implies that, once the magnetic tension is removed when the anisotropy reaches the firehose limit, the field is not able to decrease; any further decrease in B would drive the plasma unstable. In contrast, in the Braginskii limit, maintaining ∆p at the firehose limit requires dB/dt < 0. In addition, the heat fluxes always play a significant dynamical role in collisionless plasmas at high β, acting to smooth ∆p. This leads to near perfect zig-zag magnetic field lines that minimize the spatial variation in B 2 .
Compared to the Braginskii model, which may be rigorously derived from the kinetic equations via a perturbative expansion (Braginskii, 1965), Landau-fluid closures are only heuristically motivated (see Sec. 2). Nonetheless, for the clarity of presentation throughout this section, we shall primarily focus our discussion on physics contained within the LF model, viz., large heat fluxes that result from particles streaming along field lines, with no particle scattering. A variety of other physical effects that may be important (e.g., particle trapping, or particle scattering by magnetic fluctuations due to microinstabilities) are discussed in Sec. 6.

Standing waves
As discussed in Sec. 2 and more formally justified in Appendix A [see Eq. (A.34) and related discussion], the primary effect of the heat fluxes is to damp all k = 0 components of ∆p, giving Since B(t)/B(0) decreases in time as a standing wave evolves, a wave will reach Assuming a sinusoidal initial perturbation δb 0 (x) = δb 0 cos(k z z), a SA wave is interrupted if This limit, including the 8/3 numerical coefficient, matches numerical simulations using the full LF model nearly perfectly . The dynamics of a perturbation that starts above the limit (40) are illustrated in Fig. 6, which shows a solution of the LF equations (3)-(11). Despite the bizarre appearance of the highly angular, zig-zag structures that develop here [see Fig. 6(c)], the main features can be relatively easily understood. Let us consider qualitatively the wave evolution in three phases: Approach to interruption. During the initial evolution of the wave, before the anisotropy reaches ∆ = −2/β, the spatial shape of the wave is largely unaffected by the developing anisotropy. This is because spatial variation of β∆ is Evolution of a shear-Alfvén standing wave in the collisionless LF model at β = 100, at an initial amplitude above the interruption limit (40), with initial condition δb = −0.5 cos(k z z). Panel (a) shows δb at t = 0, t = 0.15τ A , and t = 3τ A (dotted black, dashed black, and blue lines, respectively), and u ⊥ /v A at t = 3τ A (red line). Panel (b) shows the magnetic-field lines of the perturbation shown in (a) at t = 3τ A . Panel (c) shows the anisotropy 4π∆p/B 2 = ∆β/2 associated with the evolution shown in (a), at t = 0, t = 0.05τ A , t = 0.15τ A and t = 3τ A (dotted, dashed, dot-dashed, and solid lines, respectively). In stark contrast to the highly nonlinear behavior of collisionless waves shown here, an MHD perturbation at these parameters is almost perfectly linear. A detailed description of each phase of evolution is given in Sec. 5.1.
due to the pressure anisotropy β∂ 2 z (δb∆) [see Eq. (20)] is thus ∼ β∆∂ 2 z δb, which simply slows down the wave without modifying its spatial structure. The pressure anisotropy during this phase is shown as a dashed line in Fig. 6(b), while δb(z) looks very similar to the initial condition [dotted line in Fig. 6(b)] at these parameters.
Early nonlinear evolution. As the pressure anisotropy approaches ∆ = −2/β, the linear term in the wave evolution equation (20), 1 + β∆/2, becomes very small and then turns negative when the anisotropy overshoots the firehose limit. This overshoot has two effects: the first is to reverse the decrease in the magnetic field of the largest-scale mode (since the linear term has changed sign); the second is to cause small-scale firehose fluctuations to grow rapidly in the regions of low δb(z) (i.e., in the neighborhood of the wave nodes; the overshoot of the firehose limit is greatest at low fields because ∆ does not vary significantly in space). These growing small-scale modes act very quickly to return the anisotropy back to its marginal level. This process can be seen in the t = 0.15τ A curves in Fig. 6(a) and (c), which show the field and the pressure anisotropy just after the small-scale firehose modes have grown at the antinodes and returned the anisotropy to the marginal level. During this phase, the presence of a spatially varying nonlinearityi.e., the O(β −1/2 ) spatial variation in ∆ [Eq. (38)] and the (1 + δb 2 ) −1 nonlinearity arising from field strength variation [see Eq. (20)] -is crucial to the dynamics, because the linear term 1 + β∆/2 is small. Without these nonlinearities, there is no preferential location for the growth of firehose instabilities and the entire wave erupts in a sea of small-scale fluctuations. It is critical (but nontrivial) to account for such a nonlinearity in a reduced equation Late-time evolution. As the firehose modes push the anisotropy back to its marginal level, the smallest-scale fluctuations decay rapidly (Melville et al., 2016). Following a transient period during which the pressure anisotropy slowly oscillates around (and decays towards) ∆ = −2/β (Schekochihin et al., 2008), the system relaxes into a final state with regions of straight fields separated by sudden corners. Despite this state's bizarre appearance, the basic cause of the plasma's preference for such structures may be inferred from a rather simple physical argument (within the LF model). This argument follows from three important properties of the collisionless dynamics: (i) without particle scattering any decrease in the magnetic field will lead to a decrease in the pressure anisotropy towards more negative values †; (ii) the only way in which the magnetic-field strength an be constant in time is either for the anisotropy to be at the firehose limit or for the field lines to be straight (or both); (iii) the heat fluxes continue to remove the spatial variation in ∆ during the slow transient phase following the initial interruption. Property (i) tells us that the magnetic field cannot continue decreasing, as it did in the Braginskii model, without creating small-scale firehose fluctuations everywhere in the plasma. Then, if we assume that the plasma has reached some quasi-steady state with nonzero B, the plasma cannot be everywhere at the firehose limit and also have ∂ z B = 0, because the heat fluxes continue to flatten ∆p [see Fig. 6(b)]. Thus, in the absence of oscillatory behavior due to the Lorentz force, properties (ii)-(iii) together imply that that B is effectively also flattened by the heat fluxes, which in turn suggests δb(z) must be piecewise constant if it is nonzero. The result is the zig-zag field lines, shown in Fig. 6(c). Note that the double-adiabatic model, which neglects the heat fluxes and so lacks property (iii), does not produce constant B fields (see Fig. A1); instead, fields with curvature may be tensionless by being everywhere at the firehose limit but with a spatially varying ∆p, as in the case of Braginskii interruption (see Sec. 4.1).
It is worth noting that, within the LF model that we use [Eqs.
(3)- (7) with Eqs. (10)-(11)], the firehose fluctuations grow fastest at the smallest scale accessible † See Eq. (18) for the evolution of the spatially averaged anisotropy. The asymptotic scalings discussed in Appendix A.3 show that the compressional term is small at high β, while the heat-flux term, 3 q ⊥ ∇ ·b , relies on magnetic curvature and so will decrease with δb. Thus, the only effect able to cancel the creation of anisotropy through 3p 0 B −1 dB/dt is the collisional damping −3ν c ∆p .  Figure 7. Evolution of a SA standing wave in the collisionless LF model at β = 100 (same parameters as Fig. 6), with initial conditions in the velocity above the interruption limit u ⊥ /v A = −0.5 sin(k z z) [i.e., initial conditions (22)]. Panel (a) shows u ⊥ /v A at t = 0 and t = 3τ A (dotted and solid red lines respectively), and δb at t = 0.25τ A (when δb is at its maximum) and t = 3τ A (dashed and solid blue lines respectively). We limit ∆ ≤ 1/β to capture heuristically the anisotropy-limiting behavior of the mirror instability; this enables δb to reach amplitudes approaching that of the initial u ⊥ /v in the simulation, which is set by an artificial hyperdiffusion operator. In reality, this scale is set by the gyroradius, where the parallel firehose growth rate decreases due to finite-Larmor radius (FLR) effects (Davidson & Völk, 1968;Schekochihin et al., 2010). A more detailed study of FLR and other kinetic effects will be the subject of future work (see Sec. 6), but it is worth noting that we see very similar macroscopic dynamics independently of the numerical resolution (for resolutions N z 128). This suggests that the exact scale separation between the SA wave and the firehose fluctuations that erupt in the "early nonlinear evolution" phase is not important for the late-time large-scale evolution (so long as the scale separation is sufficiently large).
5.1.1. Standing waves with an initial velocity perturbation. The discussion above concerned the evolution of a wave starting from a magnetic perturbation. For an initial perturbation in the velocity, the anisotropy initially grows in the positive direction, effectively increasing the restoring force of the wave. Starting from ∆p = 0 but without a mechanism to limit ∆p at positive values, this increase of ∆p as B grows is exactly the same as its decrease after B has reached its maximum, so the system never reaches ∆ < 0. This results in nonlinear standing-wave oscillations with a frequency ω > ω A and u ⊥ /v A > B ⊥ /B 0 , † which decay in time because of pressure-anisotropy damping arising from the (small) spatial variation in ∆p (see Sec. 5.2 for details). However, the mirror instability (which is excited when ∆ > β −1 ) breaks this symmetry, allowing the magnetic field to grow in time while ∆ is fixed at ∆ ≈ β −1 . As the magnetic field starts decreasing again, the mirror modes that sustained ∆ ≈ β −1 presumably decay quickly (Melville et al., 2016; see also our Sec. 6), implying that the anisotropy starts decreasing from ∆ = β −1 and can reach negative values. Thus, the limit on u ⊥ /v A will be similar to Eq. (40), with perhaps a larger numerical prefactor to account for the fact that the magnetic-field decrease starts from a positive pressure anisotropy. The process is illustrated in Fig. 7, in which we artificially limit the anisotropy to ∆ ≤ β −1 . Although ∆ = β −1 > 0 when the magnetic perturbation reaches its maximum (at t ≈ τ A /4), ‡ the resulting final state is similar to Fig. 6.

Traveling waves
A traveling wave, with δb ∝ cos(k z z − ω A t), does not change B(t) as it evolves. The arguments developed for standing waves in the previous section thus no longer apply, since the O(1) part of Eq. (38) is zero. However, even though B −1 dB/dt = 0, B −1 dB/dt itself is large. The resulting ∆p -which is reduced by a factor ∼ β 1/2 by the heat fluxes -is correlated in space with B −1 dB/dt and thus damps wave energy into thermal energy at the rate (17), viz., As the wave is damped, the resulting decrease in B(t) causes ∆p to decrease also [according to Eq. (38)], slowing down the wave and eventually causing it to stop (interrupt) if the initial B is sufficiently large for ∆ to reach −2/β. The key point here † Using Eq. (38) and assuming that ∂ z ∆ = 0, it is straightforward to derive an equation for the amplitude δb of a (sinusoidal) perturbation. After normalizing time by ω A , this is with the initial condition δb(0) = 0, ∂ t δb(0) = δu 0 = u ⊥ (0)/v A . Equation (41) is an undamped Duffing equation, and for β 1/2 δu 0 1 has the approximate solution where sn denotes the Jacobi elliptic function. These solutions oscillate in time, with a maximum amplitude δb max < δu 0 and a frequency larger then ω A (this is 1 in these time units). There is no damping of the wave, which arises from the neglected spatial variation in the pressure anisotropy. ‡ The magnetic perturbation δB ⊥ /B 0 that results from u ⊥ is also slightly smaller than the initial u ⊥ /v A because of the larger restoring force when ∆ > 0. is that during the process of wave decay, there is no mechanism to isotropize the k = 0 component of the pressure, implying any decrease in B must also be accompanied by a decrease in ∆p. The process described above is illustrated in Fig. 8, which shows the solution of the LF equations at the same parameters as the standing-wave examples (Figs. 6 and 7). At early times, the pressure anisotropy [the dotted line in Fig. 8(b)] is a strong function of space † with ∆p = 0. This spatially periodic ∆p then damps the wave, as well as causing significant nonlinear modifications to its shape (which becomes more angular, reducing B −1 dB/dt). As is also clear in Fig. 8(b), this damping drives the mean anisotropy to negative values. This effectively reduces the Alfvén speed tõ v A = v A (1 + β∆/2) 1/2 , which causes the velocity to decay faster than the magnetic field [compare red and blue dashed and solid lines in Fig. 8(a)] because δu ⊥ /ṽ A = δB ⊥ /B 0 † Note that the smoothing effect of the heat fluxes is very strong here. Without it, the early time anisotropy shown in Fig. 8(b) would be much larger. In the example shown in Fig. 8, the early time anisotropy is below the mirror and firehose thresholds (|∆| β −1 ); however, for even larger amplitudes, δb β −1/4 the wave can cause |∆| β −1 in the regions where the field is changing fastest [this estimate results from |∆| max ∼ β −1/2 δb 2 ; see Eq. (A.43)], causing even stronger nonlinear modifications to the wave.
in a traveling wave andṽ A < v A . Although not shown in Fig. 8 to avoid clutter, the velocity continues to be damped faster than the field as ∆ approaches −2/β, and we are left with an angular, magnetically dominated final state that is similar to the final state of the standing-wave evolution (Figs. 6 and 7). 5.2.1. Landau damping versus pressure-anisotropy damping. It is worth noting that the Landau damping rate of a linearly polarized SA wave due to the (nonlinear) spatiotemporal variation of the magnetic pressure is similar to the pressure-anisotropy damping, causing the wave to be damped at the rate γ/ω A ∼ β 1/2 δb 2 [Hollweg, 1971b;Lee & Völk, 1973; see discussion around Eq. (48) below]. † This effect is different from pressure-anisotropy damping (but is still captured by the LF model), and can in fact also be included in models of wave propagation that do not include a pressure anisotropy (see, e.g., Medvedev & Diamond, 1996;Medvedev et al., 1997). The difference between the pressure-anisotropy damping and Landau damping can be explained as follows. Heat fluxes are necessary for the Landau damping of a nonlinear wave, because they directly damp out the pressure perturbation that arises from the magnetic-field-strength variation. In contrast, heat fluxes act to reduce the pressure-anisotropy damping, by smoothing the spatial variation in the pressure anisotropy (in other words, the Landau damping damps the pressure-anisotropy damping!). In the discussion below, our estimate of the damping rate is heuristic, so there is no need to work out these two effects separately; however, the Landau-damping effects are included in the calculation of the wave-decay rate given in Appendix A.3.

Semi-quantitative description of traveling wave evolution
We now analyze the traveling-wave decay process in more detail, deriving a simple ordinary differential equation (ODE) to describe the process of decay and interruption. We assume that the wave remains sinusoidal throughout its evolution, which, although far from quantitatively justified (see Fig. 8), allows one to construct an ODE that describes qualitatively how the nonlinear damping and mean anisotropy affect the magneticenergy decay. This is then used to derive the decay rates of kinetic and magnetic energy, as a function of β and δb. These turn out to match reasonably well the numerical LF solutions. Our first step is to work out the decay rate of the wave due to pressure-anisotropy and Landau damping. Because this relies on the LF prescription for the heat fluxes, Here we give a heuristic derivation so as to present a relatively simple description of the important physics. Defining ∆ =∆ + ∆ k , where∆ = ∆ and ∆ k is the spatially varying part of ∆ that arises from the perturbation δb(z) with wavenumber k, the important terms in the equation for ∆ k are [see Eq. (13)] where a 1 is an O(1) dimensionless coefficient, which depends on the details of a closure for the heat fluxes. We assume a monochromatic traveling wave of amplitude δb, giving where ω A (∆) = v A k(1 + β∆/2) 1/2 accounts for the slowing down of the wave as∆ becomes negative. Since k c s ∼ β 1/2 ω A , we may neglect the time derivative † on the left-hand side of Eq. (45), leading to ∆ k ∼ β −1/2 δb 2 1 + β∆ 2 (47) † The heat fluxes suppress spatial variation in ∆ k on the time scale Evaluating the integral (43) one finds, where a 2 ≈ √ π/8 is calculated in Appendix A.3 [Eq. (A.45)]. This expression is similar in form to the damping rate (36) in the Braginskii regime, but is reduced by β 1/2 due to the smoothing effect of the heat fluxes.
Armed with the energy damping rate (48), we now formulate an equation for the slow (compared to the sound propagation, ω ∼ k c s ) dynamics of the wave amplitude δb. For a sinusoidal wave, the magnetic energy is E M = δb 2 B 2 0 /16π, while the assumption that the wave remains traveling rather than standing (i.e., that it does not generate a global oscillation in time) gives E K = (1 + β∆/2)E M . Noting also that ∆ = 3/4(δb 2 − δb 2 0 ), Eq. (48) becomes This can be reformulated in the variables ζ = βδb 2 andt = ω A β 1/2 t as which has the benefit of being controlled by just one parameter ζ 0 = βδb 2 0 . A full analytic solution to Eq. (50) is intractable, but numerical solutions (not shown) match our expectations based on the qualitative discussion in Sec. 5.2 above. Specifically, above the interruption limit (ζ 0 1), u ⊥ decays much faster than δB ⊥ in time and δB ⊥ asymptotes to a constant nonzero value at late times, whereas below the interruption limit (ζ 0 1) the nonlinear damping more equally affects u ⊥ and δB ⊥ and there is simply a slow decay of both.
We now use Eq. (50) to derive the initial decay rates of u ⊥ and δB ⊥ . This is most easily done by linearizing Eq. (50) about ζ = ζ 0 , viz., letting ζ = ζ 0 (1 + δζ) and expanding in δζ. This gives Rewriting ζ in terms of δb, substitutingt = ω A β 1/2 t, then calculating the decay rates (44) gives Although these expressions appear rather complicated, they agree nicely with our intuitive picture described earlier. In particular, the decay transitions from a regime where γ B ≈ γ U below the interruption limit βδb 2 0 1, to one where γ B γ U (with γ B independent of δb 0 ) when βδb 2 0 1.
A comparison of the damping rates (52) to the full LF traveling-wave solutions is presented in Fig. 9, where we show the decay rates measured numerically for solutions starting with a sinusoidal traveling wave. We find good agreement with the damping rates at low δb 0 , when they are small, and qualitative agreement with the trends predicted by Eq. (52) at larger δb 0 . Note in particular that the decay rate of δB ⊥ changes from increasing to decreasing with β at high δb 0 , whereas the decay rate of u ⊥ does not. The quantitative agreement at high βδb 2 0 is lacking, and there are clear reasons for this discrepancy. First, there is our assumption that the wave remains sinusoidal, which is patently not true when βδb 2 0 > 1 (see Fig. 8). The strong nonlinear shape modifications that do occur early in the evolution presumably involve some exchange of energy between u ⊥ and B ⊥ in ways that are not included in our model. Secondly, the measurement of a decay rate is ambiguous for the strongly nonlinear βδb 2 0 > 1 solutions. For simplicity, we have fit the amplitude evolution from t = 0 to t = 2τ A to a decaying exponential function, but the decay rate can vary significantly over this range at βδb 2 0 > 1. We have explored a variety of methods for determining this initial decay rate of the wave, and although the quantitative results vary with method, the general properties and qualitative agreement with the predictions (52) are robust.

Fully kinetic and multi-dimensional effects
Throughout the preceding sections, we have primarily focused on physical effects contained within the simplest 1-D Landau fluid equations (3)-(11). Importantly, the mirror and oblique firehose instabilities are not included in this model, † because these grow at k ⊥ ∼ k on the Larmor scale and thus require 2-D or 3-D kinetic simulations to be resolved correctly. In this section, we discuss -based on previous fully kinetic theory and simulations -some possible effects of these microinstabilities on the global wave evolution, focusing on which aspects of the simple 1-D picture described above are robust, and which may be modified by the inclusion of this physics. We also discuss other kinetic effects that could modify our results, including FLR effects (these were neglected by assuming k ρ i 1), particle trapping (this is not contained with the LF prescription for the heat fluxes), heat-flux limits from the gyrothermal instability, and other scattering effects. Of course, this discussion is in no way intended to be a replacement for future fully kinetic theory and simulations in two or three dimensions; rather, its purpose is to motivate the design of such studies and provide some guidance for interpreting their results. † The linear mirror instability can be captured relatively accurately by the LF model that we use (if the equations are solved in 2 or 3 dimensions; see Sec. 8 of Snyder et al., 1997). However, given the importance of trapped particles in the nonlinear mirror evolution (Schekochihin et al., 2008;Kunz et al., 2014;Rincon et al., 2015;Melville et al., 2016), it seems quite unlikely that a LF model could correctly reproduce the pressure-anisotropy-limiting behavior of the mirror instability, although we know of no relevant study that tests this (see Sec. 6.1 for discussion).

Mirror instability
In our discussion of standing waves (Sec. 5.1), the mirror instability was invoked to justify a limit on positive pressure anisotropies when starting from a velocity perturbation. This in turn allowed the magnetic field to grow, reach its maximum, and then be interrupted in a similar way to an initial purely magnetic perturbation. Without this limiting effect, an initial velocity perturbation will create an oscillating wave [albeit not a linear SA wave because the restoring force is enhanced by the positive anisotropy; see Eqs. (41) and (42)]. Thus, although the mirror instability is not crucial for the interruption effect itself, its presence does imply that the effect cannot be significantly modified based on the initial conditions. We now discuss in more detail why it is reasonable to assume that the mirror instability should have this effect.
A variety of recent kinetic results (Kunz et al., 2014;Rincon et al., 2015;Hellinger & Trávníček, 2015;Melville et al., 2016) show that mirror fluctuations, which are unstable when ∆ β −1 and cause perturbations in the field strength δB, limit ∆p by trapping particles. Namely, as the macroscopic field grows and attempts to raise the pressure anisotropy, a larger and larger fraction of particles becomes trapped in the magnetic wells and "sees" a lower field. Thus, even though the volume-averaged field continues to increase, ∆ is maintained at the marginal level β −1 because a larger proportion of particles is trapped in the ever-deepening mirror wells. During this phase, the magnetic mirrors grow in time as |δB/B| ∼ (|∇u|t) 2/3 and there is very little particle scattering because their parallel scale is significantly larger than the Larmor radius. Further, since the mirrors only saturate and start scattering particles when |δB/B| ∼ 1, they should never saturate and cause significant particle scattering for any SA wave initial conditions with u x (0)/v A ∼ |∇u|τ A < 1. The numerical experiments and arguments of Melville et al. (2016) are particularly relevant to what happens as the magnetic field reaches its maximum and starts to decrease. For β Ω i /|∇u| -the "moderate-β", or large-scaleseparation, regime most relevant to our results -the mirrors should freely decay on time scales much shorter than τ A (the decay time is ∼ β/Ω i ), releasing their trapped particles and allowing the anisotropy to decrease towards the firehose limit. Although less well understood, it seems that in the opposite limit, β Ω i /|∇u|, the firehose limit is also quickly reached (see Melville et al., 2016, Sec. 3.2), probably because the smaller-scale firehose fluctuations are able to grow on top of the larger-scale decaying mirrors. Overall, it is thus reasonable to surmise that the mirror instability will effectively act as a passive limiter, ensuring ∆ β −1 but not strongly affecting large-scale wave dynamics. † It is worth reiterating a fundamental difference between the mirror and firehose limits for SA waves. At the firehose limit, the anisotropic stress nullifies the wave restoring force (i.e., the magnetic tension). In contrast, a plasma at the mirror limit merely feels a modestly stronger (factor 3/2) restoring force. This difference explains † As the scale separation (k ρ i ) −1 is reduced, the mirrors will presumably become less effective (see Kunz et al., 2014), allowing ∆ to overshoot β −1 before acting to limit the anisotropy. Thus very large domains (compared to ρ i ) are likely essential to see these effects in fully kinetic simulations.
why the firehose limit is of much greater importance than the mirror limit for SA wave dynamics.

Oblique firehose instability
The oblique firehose instability is not included in our model, both because it operates at k ⊥ ∼ k and because kinetic theory is required for its correct description. Linearly, oblique firehose fluctuations grow faster than parallel firehose fluctuations because of their smaller scale (Yoon et al., 1993;Hellinger & Trávníček, 2008), and are also seen clearly at larger amplitudes in nonlinear regimes (Kunz et al., 2014;Melville et al., 2016). Further, unlike mirror fluctuations, the firehose fluctuations will saturate and start scattering particles after t sat ∼ β 1/2 (|∇u|Ω) −1/2 (for β Ω/|∇u|; see Kunz et al., 2014;Melville et al., 2016); i.e., after a very short time set by the microphysics.
The most obvious question that arises is then whether the differences between oblique and parallel firehose dynamics will cause significant differences in the nonlinear interruption of SA waves, compared to 1-D models where only the parallel firehose exists. This remains unclear, and understanding such issues will require fully kinetic simulations in 2-D or 3-D with k z ρ i 1. Either scenario -that the oblique firehose does or does not significantly modify the SA wave dynamics -can be plausible. On the one hand, the oblique firehose may behave similarly to the parallel firehose in 1-D simulations: be strongly excited during the early phases of wave interruption, but then die away at later times because the pressure anisotropy is pushed back above the firehose limit. On the other hand, the enhanced particle scattering in kinetic oblique firehose fluctuations could possibly continue until the magnetic field decays completely, potentially leading to collisionless SA dynamics that more closely resemble a Alfvén wave in the Braginskii regime † (see Sec. 4). However, in either case, the presence of oblique firehose fluctuations cannot circumvent the interruption limit itself -they have the same instability threshold as the parallel firehose, so are not active until the wave restoring force has already disappeared.

Other kinetic effects
Here we outline several other possible kinetic effects. Unlike the mirror and oblique firehose modes discussed above, most of these effects could be studied using 1-D, but fully kinetic simulations.
Finite-Larmor radius effects. Although various FLR effects can be included in Landau fluid models (Goswami et al., 2005;Ramos, 2005;Passot et al., 2012;Sulem & Passot, 2015), the simple closure that we used here does not include these corrections. Assuming † In support of this idea, including an artificial hard-wall firehose limit at ∆ = −2/β in a standing-wave LF simulation leads to the wave being strongly nonlinearly modified and then rather quickly decaying to oscillate with an amplitude below the interruption limit. The effect on a traveling wave is less severe, because the wave remains above the firehose limit for much of its decay. large scale separation, the most obvious effect from such corrections in 1-D is the regularization of the small scales for the parallel firehose, which has its peak growth rate at k ρ i ∼ |∆ + 2/β| 1/2 (Davidson & Völk, 1968;Schekochihin et al., 2010). Since we found that collisionless wave-interruption dynamics did not depend significantly on numerical resolution (which effectively sets the fastest-growing firehose mode in our fluid model), it seems unlikely that the direct effect of this regularization will be particularly important to the large-scale interruption (but note that the requirement k ρ i 1 could be quite severe, because |∆ + 2/β| 1/2 1 and we need significant separation between the firehose modes and the wave). There could, however, be other effects that are of some significance. For example, FLR effects enable a new instability -the "gyrothermal instability" ) -which may act to limit the heat fluxes before the SA wave hits the interruption limit, in a similar way to how the firehose instability limits the pressure anisotropy (Rosin et al., 2011). Through the gyro-viscous terms (the off-diagonal elements of the pressure tensor) and the Hall effect, † FLR effects can also act to circularly polarize the wave, creating a B y perturbation from a spatially varying B x . However, this is presumably only directly important for the macroscopic wave when the scale separation is modest, or in regions with large gradients that form during nonlinear evolution. Simple extensions of the LF model (solved numerically throughout Sec. 5) to include gyro-viscous effects and/or the Hall effect (not shown) have illustrated that these terms cause only minor changes to the SA wave evolution, so long as k ρ i is sufficiently small.
Particle trapping. Since LF closures prescribe the heat fluxes based on linear Landaudamping rates, effects of particle trapping are not included in these closures and may provide an order-unity correction to the heat fluxes. In particular, trapping can be important whenever the bounce frequency ω b of particles approaches the frequency of large-scale motions (this is ∼ ω A for a SA wave). Given that particles with velocity v and parallel velocity v are trapped if ξ = v /|v| < ξ tr ∼ |δB/B 0 | 1/2 , while |δB/B 0 | 1/2 ∼ |δB ⊥ /B 0 |, a simple estimate for the bounce frequency is ω b ∼ β 1/2 δb ω A . Thus, trapping can be important if δb > β −1/2 ; i.e., for a wave above the interruption limit. Trapping has the effect of reducing the Landau damping rate of nonlinear SA traveling waves (O'Neil, 1965;Kulsrud, 1978) and presumably also modifies pressureanisotropy damping. These effects will be considered in detail in future work.
Other scattering effects. If the "corners" that develop in the magnetic-field lines (e.g., Figs. 6-8) are on the Larmor scale, unresolved by our LF closure, these may scatter particles. This would provide an interesting case where a plasma could set its own mean free path λ mfp ∼ k −1 based on the large-scale driving. If real, this effect would most significantly modify traveling-wave dynamics, because the square structures that develop before interruption (see Fig. 8) could possibly cause sufficient scattering to damp † The Hall effect becomes important when the global pressure anisotropy fast enough so that the wave decayed before reaching the interruption limit.
Overall, we would like to stress that although the details of wave interruption may be modified by the addition of other kinetic physics, our basic result -that weakly collisional SA waves cannot exist in their linear form above the limits (30) and (40) -is robust. The dominance of magnetic energy over kinetic energy is also a generic consequence of interruption, because in the approach to the firehose limit, the equipartition of energy in an Alfvén wave is modified by the decrease of magnetic tension. We find generic agreement on these points between the Landau fluid, Braginskii, and double-adiabatic models (q ⊥ = q = 0; see Appendix A.2). Many of our key results are thus quite insensitive to the form of the heat fluxes or particle scattering, relying purely on the physics of pressure-anisotropy generation in a changing magnetic field.

Conclusion
In this paper, we have explored the nonlinear "interruption" and damping of linearly polarized shear-Alfvén (SA) waves in weakly collisional plasmas. These effects, which arise due to the pressure anisotropy that is generated in the changing magnetic field of the wave, lead to a limit on the amplitude of propagating/oscillating SA waves in the collisionless regime: In the weakly collisional Braginskii limit, which applies when ν c ω A , propagating/oscillating SA waves are also limited in amplitude, to We summarize our main findings as follows: • Above the limit (53), collisionless SA waves are "interrupted" when their selfgenerated pressure anisotropy reaches the firehose boundary ∆p ≈ −B 2 /4π. At this boundary, the wave's restoring force (the Lorentz force) is cancelled by the anisotropy, and the magnetic energy dominates the kinetic energy because the effective Alfvén speed goes to zero.
• Due to the correlation between ∆p and B −1 dB/dt, there is a net transfer of wave energy to thermal energy of the plasma through "pressure-anisotropy heating" at the rate´dx ∆pB −1 dB/dt. This results in a nonlinear damping of the wave, even below the limit (53).
• Heat fluxes are always important in the high-β collisionless limit (because the thermal velocity is larger than v A ) and act to smooth the spatial dependence of the pressure anisotropy.
• In the collisionless limit, standing and traveling SA waves behave in qualitatively different ways because the spatial average of B decreases during a standing wave's evolution, whereas it does not for a traveling wave. Thus, while a standing wave above the limit (53) is interrupted within half a wave period, a traveling wave is first nonlinearly damped [at the rate ∼ ω A (δB ⊥ /B 0 ) 2 β 1/2 ], leading to a decreasing B and eventual interruption of the wave.
• The kinetic energy in a collisionless traveling wave is damped significantly faster than the magnetic energy for amplitudes approaching (or exceeding) the limit (53), and the magnetic energy can be a large fraction of its initial value when the wave interrupts. This occurs because, as the wave decays, the global decrease in ∆p reduces v A , which changes the ratio of u ⊥ and δB ⊥ (this also slows down the wave; see Fig. 8).
• Barring additional kinetic and higher-dimensional effects not contained within our Landau-fluid model (see Sec. 6), the outcome of wave interruption is the creation of a magnetically dominated state of nearly perfect zig-zag magnetic field lines (see Figs. 6 and 8) -i.e., a quasi-periodic pattern with spatially constant magnetic field strength. The emergence of this state may be understood by noting that it is the only state that has both zero magnetic tension and a spatially smooth pressure anisotropy along the field lines (because spatial variation in ∆p is damped by the heat fluxes).
• Wave interruption in the Braginskii limit involves a slow decay of the wave over the timescale t decay ∼ β/ν c (δB ⊥ /B 0 ) 2 , which occurs because a slowly changing B is necessary to maintain the anisotropy at the firehose limit. The characteristic fieldline structures (Fig. 2) differ from collisionless waves because the magnetic tension is zero if the anisotropy is at the firehose limit, even if there is spatial variation in B 2 .
• The pressure-anisotropy damping of SA waves is large in a Braginskii plasma because the spatially varying part of ∆p is comparable to its mean. For waves below the limit (54), this leads to the wave energy being damped at the rate ∼ ω 2 A /ν c (δB ⊥ /B 0 ) 2 β. † • The amplitude limits do not apply to circularly polarized SA waves because for these, dB/dt = 0.

Implications and applications
Given the ubiquity of shear-Alfvén waves in plasma physics (see Sec. 1), the stringent limits on their amplitude derived here may have interesting implications in a variety of hot, low-density (and therefore, weakly collisional) astrophysical plasmas. Although a detailed study of all these implications is beyond the scope of this work, it is worth commenting on magnetized turbulence in particular, given its importance in many subdisciplines of astrophysics. The salient point is that in well-accepted phenomenologies of strong magnetized turbulence (in particular, Goldreich &Sridhar, 1995, andextensions, e.g., Boldyrev, 2006), the physics of shear-Alfvén waves is critical at all scales in the turbulence. † A strong modification to SA wave physics would thus be expected to significantly modify the turbulent cascade. One might expect such modifications to be even stronger in the weak turbulence regime (Ng & Bhattacharjee, 1996;Galtier et al., 2000;Schekochihin et al., 2012) given the relative weakness of nonlinear interactions in comparison to the SA wave physics in such turbulence. More explicitly, turbulence in a weakly collisional high-β plasma may depend on the amplitude of its forcing. Since velocities are strongly damped when a wave is interrupted, it may behave as a fluid with Reynolds number 1 when u ⊥ /v A β −1/2 (and u ⊥ v A ; otherwise the waves are of such large amplitude that the turbulence would likely be in a dynamo regime, which is not particularly well understood even in MHD, but is less obviously Alfvénic). In contrast, for perturbations of amplitude u ⊥ /v A β −1/2 , the linear SA wave physics is mostly unaffected, and a standard Alfvénic cascade should develop. Since pressure-anisotropy heating is able to dissipate large-scale wave energy when wave interruption is important, a turbulent cascade may not be necessary for the plasma to absorb the energy injected by a continuous mechanical forcing (Kunz et al., 2010). While further study is necessary to understand this physics better, it is at least clear that the immediate interruption of SA fluctuations with amplitudes exceeding δB ⊥ /B 0 ∼ β −1/2 should significantly limit the application of MHD-based turbulence phenomenologies to high-β weakly collisional plasmas.

Future work
A first priority for future studies of wave interruption is the inclusion of the kinetic effects discussed in Sec. 6. Unfortunately (from a computational standpoint), due to the 2-D kinetic nature of the mirror and oblique firehose instabilities, significant progress in this endeavor requires kinetic simulations in two spatial dimensions and three velocity-space dimensions. Since firehose instabilities grow on scales somewhat above the gyroscale (see Sec. 6), if one hopes to study the asymptotic regime k ρ i 1, the required scale separation between the gyroscale and the SA wave is likely quite large. We thus expect a detailed kinetic study to be rather computationally expensive, although certainly feasible. That said, there will also be a variety of interesting insights to be gained from purely 1-D kinetics: for example, the role of the gyrothermal instability, particle scattering off magnetic discontinuities, and particle-trapping effects. Further, the behavior of SA waves with limited scale separation between k −1 and ρ i is also of interest physically, in particular for the solar wind, where observations easily probe turbulent fluctuations down to the ion gyroscales and below. A separate line of † A more common way to say this is that the cascade is in critical balance (Goldreich & Sridhar, 1995), which states that the linear (Alfvén) time is equal to the nonlinear turnover time at all scales in a strong MHD cascade.
investigation for future work involves applications of the amplitude limit, in particular to turbulence, as discussed in Sec. 7.1. This may be productively pursued using a 3-D Landau fluid code (as in Sharma et al., 2007;Sharma et al., 2006) or using Braginskii MHD.
Overall, although many questions remain, both the limit δB ⊥ /B 0 β −1/2 itself and the strong dominance of magnetic over kinetic energy are robust, appearing across a variety of models. Given the stringent nature of the amplitude limit and the interesting implications for high-β magnetized turbulence in weakly collisional plasmas, we anticipate a range of future applications to heliospheric, astrophysical, and possibly laboratory (Forest et al., 2015;Gekelman et al., 2016) plasmas. the spatially variation in ∆ plays a key role in forming square structures as the wave approaches the interruption limit. This motivates two separate expansions: the first is valid when 1 + β∆/2 ∼ O(1), the second is valid when 1 + β∆/2 ∼ O( 2 ) (i.e., when the wave has already evolved to be close to the interruption limit).
These difficulties motivate our arrangement of this appendix as follows. We start in Appendix A.2 by considering the double-adiabatic version of Eqs. (3)-(7), with q ⊥ = q = 0. This leads to a simple nonlinear wave equation that is free from the issues mentioned in the previous paragraph because there is large spatial variation in ∆. We then consider the initial wave evolution using the Landau-fluid closure in Appendix A.3, ordering u ⊥ /v A ∼ δB ⊥ /B 0 ∼ O( ) and 1 + β∆/2 ∼ O(1), which can be used to find the initial damping rate of a traveling wave. Finally we derive an equation for waves as they get very close to the interruption limit itself in Appendix A.4. This involves a spatially varying nonlinearity arising from both the field curvature and the spatial variation of ∆. Unfortunately, a closure problem prevents true asymptotic determination of the evolution of the spatial mean of ∆, although the expansion is still helpful for determining the residual spatial variation of ∆ and formulating a simple nonlinear wave equation that describes the wave's approach to zig-zag field-line structures. . This is because a nonlinearity with spatial dependence appears in the lowest-order wave equation, due to the spatial variation of ∆ being comparable to its mean. For this reason we start by outlining the procedure for the double-adiabatic equations, even though the neglect of the heat fluxes is not a valid approximation in the high-β limit. Our asymptotic ordering, motivated by our interest in solutions near the interruption limit with δB ⊥ /B 0 ∼ β 1/2 0 , is In addition, B z = 1 due to ∇ · B = 0. This leads to This says that p 2 has no spatial variation, p 2 = 0, expressing the parallel pressure balance.
Order 1 . The perpendicular velocity at O( ) satisfies , it is clear that we need an expression for ∆ 2 , and, for the system to be closed, this must depend only on u x1 and B x1 .
Order 2 . To calculate ∆ 2 , we require only p ⊥2 and p 2 , so may ignore the momentum (A.3) and induction (A.4) equations at this order. Noting thatbb : ∇u =b zbz ∂ z u z + b xbz ∂ z u x , the perpendicular and parallel pressure equations (A.5)-(A.6) become From the O( 0 ) parallel momentum equation (A.21), we know that p 2 = 0, so, using ∂ z u z2 = 0, we obtain Combining this with the perpendicular pressure equation (A.24) and assuming ∆ 2 (t = 0) = 0, we find  As should be expected because of the different spatial form of the nonlinearity, the nonlinear evolution at the interruption limit is quite different from the evolution when heat fluxes are included (cf. Fig. 6). In particular, since the heat fluxes no longer act to spatially smooth the pressure anisotropy, B 2 can vary in space even with the anisotropy everywhere close to the firehose limit. The model equation does a reasonably good job at capturing the main qualitative features of interruption and is nearly perfect for the initial wave interruption (compare dashed lines in each panel). It is worth noting that the relative spatial variation of p ⊥ compared to its mean (and the lack of variation in p ) can also be obtained by considering the compressible part of the CGL equations as a forced oscillator system (see , and this agrees with numerical solutions of the full equations (A.2)-(A.6).
Appendix A.3. Landau-fluid closure: Initial evolution Here, we repeat the calculation of the previous section but include the Landau-fluid prescription (A.7)-(A.8) for the heat fluxes. These act to smooth pressure perturbations on the sound-crossing timescale, leading to a pressure perturbation that is constant in space to lowest order, i.e., ∆ 2 = 0. The resulting equation for B x is thus not accurate when 1 + β 0 ∆/2 1 (i.e., in the approach to the interruption limit), because there is no spatially local nonlinearity to steepen the wave into zig-zag structures during the slow dynamics when β 0 ∆/2 ≈ −1. As well as motivating the use of a second expansion with 1+β 0 ∆/2 ordered small (this is done in the Appendix A.4), the expansion presented here is used to calculate the damping rate of a traveling wave due to the spatial correlation of B −1 dB/dt and ∆, as used in the arguments in Sec. 5.2. This is nonzero only at order O( 3 ), because ∆ 2 = 0.
We use the same ordering as the previous section, Eqs. (A.11)-(A.16). To lowest order, the heat fluxes (A.7)-(A.8) simplify significantly to while ∇ · (q ⊥b ) = 2 ∂ z q ⊥,2 + 3 ∂ z q ⊥,3 + O( 4 ) (and similarly for q ). These simplifications are tantamount to stating that the heat flows along the mean field (i.e., alongẑ) at the lowest two orders. As we did for the double-adiabatic calculation in Appendix A.2, let us go through each order of the expanded equations. These are where we also know from order O( ) that ∂ z u z2 = 0, p 2 = p 2 , and p ⊥2 = p ⊥2 This leads to the following wave equation for B x1 : which is the same as Eq. (20) in the main text, with Eq. (14) for the anisotropy and neglecting the δb 3 terms that arise from the magnetic curvature. The problems with using Eq. (A.38) to describe the wave near 1 + β 0 ∆/2 = 0 are immediately apparent: no matter how small one takes , there will always be a time for which the higher-order contributions from the magnetic curvature and spatial variation of ∆ play an important dynamical role. Indeed, numerical solutions to Eq. (A.38) stay perfectly sinusoidal until 1 + β 0 ∆/2 < 0, at which point small-scale (firehose) fluctuations grow rapidly. There is no tendency for the wave to become square. This issue will be resolved in Appendix A.4 through the use of a different ordering scheme.
Appendix A.3.1. Wave damping. A traveling wave, which satisfies B 2 x = const., propagates linearly, with no nonlinear modification, under Eq. (A.38). Although a continuation of the expansion to higher order is not very useful under this ordering, one can obtain an estimate of the lowest-order contribution to the damping of the wave energy into thermal energy that occurs for a traveling wave due to the spatial dependence of p ⊥3 . In the dimensionless variables (A.1), the kinetic-energy evolution Eq. (17) is The right-hand side of this equation includes compressional and pressure-anisotropy heating, which can cause the transfer of mechanical wave energy into thermal energy.
where we have split ∆ 3 into its mean and spatially varying parts, ∆ 3 and ∆ 3 . The compressional term p ∇ · u contributes only at order O( 5 ) and higher, because ∂ z p 2 = ∂ z p 3 = 0 [Eq. (A.31)] and ∂ z u z2 = 0 [Eq. (A.34)]. The O( 2 ) term on the right-hand side of Eq. (A.40) is zero for a traveling wave, because ∂ t B 2 x1 = 0, and similarly for the second of the order O( 3 ) terms. Since we are interested in the damping of a pure sine wave, we shall also ignore the first order O( 3 ) term, which is related to the development of shape changes in B x . This leaves us with (β 0 /2) ∆ 3 ∂ t (B 2 x1 ) in the right-hand side of Eq. (A.40). This term describes how the average spatial correlation of ∆ with B −1 dB/dt causes a net damping, even when the averages of ∆ and B −1 dB/dt individually are each zero.
To calculate ∆ 3 , consider the spatially varying part of Eqs. (A.35) and (A.36): Using p 3 = 0 [Eq. (A.31)], solving for ρ 3 , and inserting this solution into the p ⊥ equation gives Therefore, the third-order wave-damping rate is For a traveling sinusoidal wave B x1 = δb cos(z − t), by carrying out the spatial integrations using |k z | −1 B 2 x1 = cos(2z − 2t)/4, we find that the wave energy is damped at the rate ∂ t E wave = √ π 3 β 1/2 0 δb 4 , or, restoring dimensions, In this section, we derive a nonlinear wave equation to describe the final approach to wave interruption, which is not captured correctly by Eq. (A.38) because it lacks a spatially dependent nonlinearity. Specifically, the distance from marginality, 1 + β 0 ∆/2, is ordered as O( 2 ) even though 1 and β 0 ∆/2 are each O( 0 ). Since the previous expansion yielded the result that ∆ = 0 to lowest order [Eq. (A.34)], this may be considered as a re-ordering of the equations, which becomes valid when Eq. (A.38) loses its validity because 1 + β 0 ∆/2 1. Under the assumption of small 1 + β 0 ∆/2, we are also forced to assume u x ∼ B x and ∂ t ∼ , as should be expected. † The resulting wave equation (A.60) contains spatially dependent nonlinearities from both the magnetic curvature and the spatial variation of ∆. It thus contains the terms necessary to reproduce zig-zag field-line structures seen in solutions of the full LF equations (e.g., Fig. 6).
Given these considerations, our asymptotic ordering is modified from Eqs. (A.11)-(A.16) as follows: where f represents any variable [the change in the ordering of u z stems from the time derivative in the continuity equation (A.2)]. Before embarking on an order-by-order expansion, we can simplify our task significantly by noting that only every second term in each field expansion need be considered, viz., B x = B x1 + 3 B x3 + O( 5 ), p ⊥ = 1 + 2 p ⊥2 + 4 p ⊥4 + O( 6 ), ∆ = −2/β 0 + 4 ∆ 4 + O( 6 ) etc., for all fields. This is justified by the fact that in all equations, the expressions for order-n quantities depend † The reader may notice that the spatially varying part of ∆ was O( 3 ) in the expansion of the previous section [Eq. A.43], whereas here, by assuming 1 + β 0 ∆/2 ∼ 2 we are effectively ordering it to be O( 4 ). This apparent discrepancy is resolved by noting that as the solutions of Eq. (A.38) evolve towards 1 + β 0 ∆/2 = 0, the time derivatives (or equivalently u x ) also become one order smaller, meaning that the spatial variation of ∆ is pushed into ∆ 4 (recall that ∆ 3 was determined by the current value of ∂ t B 2 x1 , not its time history). Thus, the solutions of Eq. (A.38) will evolve into a regime where the expansion discussed in this section in valid.
only on order-(n + 2) quantities, because the terms that relate to modifications in B all contain B 2 x . † We must still work through all orders of the equations in , since some fields contain even powers of (e.g., u x , ∆) while others contain odd powers (e.g., B x , u z ).
Expandingb z and using ∆ = −2/β 0 + O( 4 ), this leads to The perpendicular induction equation appears at this order, giving as expected.
Order 3 . The perpendicular momentum equation, which forms the basis for our wave equation, appears at O( 3 ) and reads order O( 5 ). Noting that ∂ z u z3 = 0, the pressure equations (A.5)-(A.6) at O( 3 ) are then, We can obtain useful information from both the spatial average and the spatially varying which implies that ∂ t B 2 x1 ∼ 2 ∂ t B 2 x1 , i.e., that the spatial average of B 2 x1 varies in time more slowly than B 2 x1 itself. This can occur, for instance, if B x1 is increasing in some region and decreasing in another, as could occur in the approach to a square wave.
The spatially varying part of Eqs. (A.54) and (A.55) can be used to solve for ∆ 4 . We first require the heat fluxes (A.7)-(A.8), which are, † This may be inserted into Eqs. (A.52) and (A.53) to yield the nonlinear wave equation (A.60) Unfortunately, Eq. (A.60) still contains the undetermined mean pressure anisotropy ∆ 4 . While in principle, one can solve for ∆ 4 by considering the mean part of the pressure equations (A.5)-(A.6), the result contains B x3 , so ∆ 4 remains unknown. Of course, any attempt to subsequently solve for B x3 generates dependence on ∆ 6 , leading to a standard closure problem. Despite this issue, Eq. (A.60) remains useful for a number of reasons. First, the spatially dependent nonlinearities are interesting: because of the time derivative in ∆ 4 [the first term on the right-hand side of Eq. (A.59)], this term has a diffusive effect in Eq. (A.60), and it can dissipate wave energy into thermal energy. This is not the case for the the ∂ 2 z (B 3 x1 ) nonlinearity in Eq. (A.60), which arises from spatial variation in the field strength. Secondly, the exact form of ∆ 4 plays only a minor role, because it is simply a spatially constant number that must decrease as B 2 decreases † Note that the p ⊥ (1 − p ⊥ /p )∇ B/B part of q ⊥ in Eq. (A.7) has made an appearance at this order. (since ∂ t ∆ ∼ ∂ t B 2 ). Indeed, it has only one property that is key to the dynamics described by Eq. (A.60) -it must be able to approach 0 and become negative, so as to slow the linear dynamics and allow the nonlinear terms to dominate. Corrections (at order 4 ) to the exact point at which this zero crossing occurs will presumably not affect the dynamics of the wave strongly. For the purposes of exploring solutions to Eq. (A.60) numerically (see Fig. A2), we thus make the simple ansatz ∆ 4 = ∆ 4 (0) + 3 4 β 0 (B 2 x1 − B x1 (0) 2 ) . (A.61) Here ∆ 4 (0) is the initial anisotropy, which would arise through dynamics that satisfy Eq. (A.38). The second term comes from the form of ∆ 2 in Eq. (A.38) [but it is now of order O( 4 ) for the same reason that we obtained Eq. (A.56)], which is a clear choice that satisfies the requirements discussed above. As shown in Fig. A2, the nonlinear wave equation (A.60) with ∆ 4 given by Eq. (A.61) has solutions that are pleasingly similar to those of the full LF equations (cf. Figs. 6-8): we see a clear approach to zig-zag field lines, a much faster decay of the velocity in comparison to magnetic field, the eruption of small-scale firehose modes at early times for a standing wave, and the slowing of a traveling wave as the anisotropy approaches the firehose limit.

Appendix B. Asymptotic wave equations -Braginskii limit
In this appendix, we derive asymptotic wave equations in the Braginskii limit, with ω A ν c . As expected, this calculation is significantly simpler than the collisionless cases discussed Appendix A. We find two regimes, with the dynamics controlled by the parameterν c /β 1/2 0 (recall thatν c ≡ ν c /ω A ), which determines whether the effect of heat fluxes is important for the spatial form of the pressure anisotropy. In the first regime, whenν c β 1/2 0 , one recovers the Braginskii wave equation discussed in the main text [Eq. (20) with the closure (9)]; in the second, whenν c ∼ β 1/2 0 (orν c < β 1/2 0 ), one finds an anisotropy that becomes smoother in space as β 1/2 0 /ν c increases. The magnetic field and velocity dynamics in each regime are generally similar to each other, heuristically vindicating the neglect of the heat fluxes in the main text. Throughout this appendix, we use the definitions and dimensionless equations described in Appendix A.1.
To constrain further the ordering ofν c and β 0 individually, let us consider the basic scaling of the heat fluxes in the Braginskii regime. Ignoring -for reasons that will become clear momentarily -the effect of the collisionality on the heat fluxes, let us consider the pressure equations (A.5) and (A.6). Sinceν c must dominate over ∂ t , these equations, to lowest order, will be dominated by the terms β 1/2 0 ∂ z q ⊥, , B x ∂ z u x , andν c ∆. The balance between the collisions and heat fluxes is thus controlled byν c /β 1/2 0 : ifν c ∼ β 1/2 0 , the heat fluxes will enter at the same order as collisional isotropization in the pressure equations, while ifν c β 1/2 0 , the heat fluxes will simply cause minor (higher-order) corrections to the spatial form of ∆. Our final equations will thus depend on the ordering ofν c /β 1/2 0 . This leads us to two natural choices: the "high-collisionality regime," ν c ∼ O( 4 ), β 0 ∼ O( 6 ), and and "the moderate-collisionality regime," ν c ∼ O( 2 ), β 0 ∼ O( 4 ). † In the discussion above, we neglected to mention collisional modifications to the heat fluxes [see Eqs. (A.7)-(A.8)]. For the sake of qualitative discussion, this neglect is admissible because collisions reduce the heat fluxes at the same point,ν c ∼ β 1/2 0 , as † Note that the regime in which the heat fluxes dominate, viz.,ν c β 1/2 0 , will turn out to be a subset of the moderate-collisionality regime in Appendix B.3. There is thus no need to treat it separately. the heat fluxes become subdominant to SA wave dynamics. This can be seen from the form of the LF heat fluxes, Eqs. (A.7)-(A.8). Expanding these assuming small p and ρ perturbations, and ignoring numerical coefficients, one finds that the heat fluxes scale as Combined with the discussion of the previous paragraph, this shows that in the limit ν c β 1/2 0 , where the heat fluxes were not important in the pressure equation, the heat fluxes are even further reduced, making their neglect more valid than it would otherwise be. In contrast, whenν c ∼ β 1/2 0 , the heat fluxes are only moderately affected by collisionality. Thus, the effect of collisionality on the heat fluxes is simply to improve the validity of the high-collisionality ordering, while changing the moderate-collisionality results by O(1) numerical factors.  7) and (B.15), starting from a sinusoidal magnetic perturbation B x1 = −0.5 cos(2πz). In the high-collisionality case shown in panels (a)-(b) we takeν c = 5 4 = 625, β = 5 6 = 15625. In the moderatecollisionality case shown in panel (c)-(d) we takeν c = 5 2 = 25, β = 5 4 = 625, such thatν c = β 1/2 (we take ζ(ν c ) = (2π) −1/2 for simplicity). These parameters give the same interruption limit ν c /β = 0.2 in both cases. In the top panels, we show B x1 at t = 0 (dotted black line), B x1 at t = 0.6τ A (solid blue line), u x1 at t = 0.6τ A (solid red line), and B x1 at t = 2τ A (dashed black line), which is after the wave has decayed to below the interruption limit. In the bottom panels, we show the pressure anisotropy, β 0 ∆/2, at the same times. Although the pressure anisotropy profiles are quite different in each case [compare panels (b) and (d)], the dynamics of the magnetic perturbation, including the time taken for the wave to decay, are quite similar. This is because the parts of the wave where ∆ > −2/β have B x1 = 0 anyway (see discussion in text). Note that the sole difference between the calculation shown in Fig. 2 and that in panel (a) here is the (1 + δb 2 ) −1 field nonlinearity term, which is not included here because it is at higher asymptotic order. Because the spatially varying nonlinearity due to ∆p is larger than that due to δb 2 , this term makes little difference to the dynamics.
Put together, Eqs. (B.11) and (B.12), along with the perpendicular momentum equation (B.5), lead to the wave equation  Fig. B1(b)], when the heat fluxes significantly modify the pressure anisotropy, the dynamics are largely similar to the basic high-collisionality Braginskii limit discussed in the main text [ Fig. B1(a); see also Fig. 2]. The reason is for this is related to the nature of the Braginskii wave decay, as discussed in Sec. 4. Effectively, the dynamics of the decaying wave separate into regions where ∆ = −2/β 0 and the field has curvature, and regions where ∆ > −2/β 0 and the field has no curvature (i.e., where the perturbed field δb is zero). The primary effect of the heat fluxes is thus to decrease the anisotropy where B x1 is already zero anyway, causing only small modifications to the dynamics of the standing wave. It is worth noting, however, that because ∆p is smoothed more by the heat fluxes as ν c /β 1/2 0 decreases [i.e., the ratio of ∆ 4 (B.12) to ∆ 4 (B.11) decreases as ν c /β 1/2 0 decreases], the decay rate of a Braginskii traveling wave will be reduced by a factor between 1 and β 1/2 0 [compared to the estimate in the main text; see Eq. (37)], when the heat fluxes cause significant smoothing of the pressure anisotropy.