On hydromagnetic wave interactions in collisionless, high-β plasmas

We describe the interaction of parallel-propagating Alfvén waves with ion-acoustic waves and other Alfvén waves, in magnetized, high-$\beta$ collisionless plasmas. This is accomplished through a combination of analytical theory and numerical fluid simulations of the Chew–Goldberger–Low (CGL) magnetohydrodynamic (MHD) equations closed by Landau-fluid heat fluxes. An asymptotic ordering is employed to simplify the CGL-MHD equations and derive solutions for the deformation of an Alfvén wave that results from its interaction with the pressure anisotropy generated either by an ion-acoustic wave or another, larger-amplitude Alfvén wave. The difference in time scales of acoustic and Alfvénic fluctuations at high-$\beta$ means that interactions that are local in wavenumber space yield little modification to either mode within the time it takes the acoustic wave to Landau damp away. Instead, order-unity changes in the amplitude of Alfvénic fluctuations can result after interacting with frequency-matched acoustic waves. Additionally, we show that the propagation speed of an Alfvén-wave packet in an otherwise homogeneous background is a function of its self-generated pressure anisotropy. This allows for the eventual interaction of separate co-propagating Alfvén-wave packets of differing amplitudes. The results of the CGL-MHD simulations agree well with these predictions, suggesting that theoretical models relying on the interaction of these modes should be reconsidered in certain astrophysical environments. Applications of these results to weak Alfvénic turbulence and to the interaction between the compressive and Alfvénic cascades in strong, collisionless turbulence are also discussed.


Introduction
Until recently, our understanding of Alfvén waves (AWs) at scales much larger than plasma-kinetic scales has been largely agnostic of the rate of Coulomb collisions.Indeed, linear shear AWs and nonlinear torsional AWs do not change the form of the velocity distribution function of the particles, but rather define the moving frame in which any changes to it are to be measured (e.g., Schekochihin et al. 2009).As a result, it has been predicted that many defining characteristics of Alfvénic turbulence, such as its spectral indices or spatial anisotropies, are at most subtly affected by a lack of collisions.In introducing the concept of AW interruption however, Squire et al. (2016Squire et al. ( , 2017b) ) established that shear AWs in collisionless, high-β plasmas are susceptible to non-equilibrium effects that place constraints on their ability to propagate freely, complicating the Alfvénic dynamics in new and novel ways. 1 This is just one example of how even small deviations from local thermodynamic equilibrium (LTE), made possible by low collisionality, can have large dynamical consequences for the evolution of highβ systems.As a result, there has been increased interest in revisiting many problems in collisionless, high-β plasmas that are fundamentally connected to Alfvénic dynamics, such as turbulence (Bott et al. 2021;Arzamasskiy et al. 2023;Squire et al. 2023), the fluctuation dynamo (Santos-Lima et al. 2014;Rincon et al. 2015;St-Onge & Kunz 2018;St-Onge et al. 2020;Sironi et al. 2023;Zhou et al. 2023), and magnetic reconnection (Cassak et al. 2015;Alt & Kunz 2019;Winarto & Kunz 2022;Egedal et al. 2023).
In well magnetized and collisionless plasmas, the importance of these non-equilibrium effects can be parameterized by the product β∆, where ∆ .= p ⊥ /p ∥ − 1 is the normalized difference between thermal pressures across (⊥) and along (∥) the local magnetic-field direction.Aside from determining whether a plasma is unstable to various pressureanisotropy-driven microinstabilities, the quantity β∆ also quantifies the competition between the plasma's pressure anisotropy and the magnetic tension, with the characteristic propagation speed of magnetic disturbances being the effective Alfvén speed v A,eff .= v A 1 + β∆/2.When β∆ ∼ ±1, the restoring tension force can be significantly enhanced or even entirely eliminated, rendering Alfvénic fluctuations impotent.If β ≳ 10, then such circumstances are rather easy to achieve by even small deviations from LTE.
The collisionless high-β plasmas in which these effects appear are in fact pervasive throughout the universe.Indeed, from radiatively inefficient accretion flows onto supermassive black holes, to the hot and dilute intracluster medium of galaxy clusters, to certain regions of the solar wind at and beyond the Earth's orbit, it is often the case that Coulomb collisions are rare and magnetic fields, while strong enough to magnetize particles, are energetically weak (Kunz, Squire et al. 2019).Further, turbulence is a ubiquitous phenomenon within these systems, and it plays a crucial role in processes like plasma heating (e.g., Sharma et al. 2007;Kunz et al. 2018;Kawazura et al. 2019;Arzamasskiy et al. 2023) and angular-momentum transport (e.g., Sharma et al. 2006;Kunz et al. 2016;Bacchini et al. 2022;Sandoval et al. 2023).To describe these processes within high-β astrophysical environments accurately, some refinement of our understanding of their turbulent fluctuations is necessary.Therefore, with a focus on the role of β∆, this work revisits the following three tenets of Alfvénic turbulence that are known to hold in magnetohydrodynamics (MHD) and even β ∼ 1 collisionless gyrokinetics: (i) that co-propagating AWs do not interact with one another, (ii) that weak interactions between AWs do not modify their field-aligned wavenumbers, and (iii) that compressive fluctuations in the inertial range are only passively mixed by the Alfvénic cascade 2 .
The assertion that co-propagating AWs do not interact with each other in Alfvénic turbulence (point (i)) is a simple result of v A being relatively insensitive to nonlinearities in β ∼ 1 plasmas.This effect arises naturally from the MHD equations and governs wave interactions across both weak and strong turbulence.Specific to weak AW interactions however (point (ii)), one can effectively apply energy and momentum conservation to 1 β .= 8πp/B 2 is the ratio of the thermal pressure of the plasma particles, p, and the energy density stored in the magnetic field, B 2 /8π, or equivalently twice the ratio of the squares of the sound speed cs and the Alfvén speed vA, viz.β .= 2c 2 s /v 2 A .In this paper, 'high-β' means β ≫ 1. 2 We are primarily concerned with effects on the inertial range of Alfvénic turbulence where the fluctuations are spatially anisotropic with respect to the background magnetic field with k ∥ /k ⊥ ≪ 1 (e.g., Goldreich & Sridhar 1995).Within intervals of the cascade not exhibiting this predominantly perpendicular polarization (e.g., near the forcing scale(s)), strong interactions between slow and Alfvén waves can occur (e.g., Howes & Nielson 2013) pairs of interacting waves since they remain correlated for many linear times.Those constraints, combined with the fact that all Alfvénic fluctuations propagate at v A , prevent the generation of a cascade in k ∥ .Lastly, for point (iii), inertial-range compressive fluctuations being passively mixed by AWs results from the lack of overlap between the linearized eigenvectors of acoustic and Alfvénic fluctuations in β ∼ 1 plasmas when k ⊥ ≫ k ∥ .None of the perturbed quantities of the acoustic waves enter to linear order in the propagation of AWs.However, among other things, AWs distort the direction of the magnetic field along which acoustic waves linearly perturb the thermal pressure.As such, they are sensitive to the presence of Alfvénic fluctuations.Yet, all of these phenomena may be expected to change in collisionless high-β plasmas.
To understand the need for revisiting these wave interactions, consider the example of AW interruption (Squire et al. 2016(Squire et al. , 2017b)).In the interruption process, a transverse perturbation (δB ⊥ ) to the magnetic field begins to oscillate at the Alfvén frequency and initially decreases in amplitude.As a result, the field strength |B| ∼ B 0 + δB 2 ⊥ /2B 0 decreases nonlinearly.Conservation of the double adiabats p ⊥ /ρB and p ∥ B 2 /ρ3 then dictates that negative pressure anisotropy ∆ < 0 is produced.If the amplitude is large enough (δB ⊥ /B 0 ≳ 2β −1/2 ), a pressure anisotropy of ∆ ∼ −2β −1 can be produced by the changing |B|, although limited to ∆ = −2β −1 by the Alfvénic firehose instability threshold.However, the rate of anisotropy production is Alfvénic, while p ⊥/∥ are diffused by heat fluxes operating on sonic timescales (τ sonic ∝ L/c s ≪ L/v A , where c s ≡ (2p ⊥ /3 + p ∥ /3)/ρ is the sound speed).The anisotropy is then smoothed out by the heat fluxes and fills the domain with a nearly constant ∆ = −2β −1 .At this point, the magnetic tension that supports the AW is nullified and the fluctuation ceases propagation, hence the term 'interruption'. 3If an AW is capable of experiencing a dramatic change in its propagation due to its own nonlinearly generated anisotropy, then one would expect that it would similarly respond to the anisotropy generated by other fluctuations within a turbulent system.It is for this reason that we revisit wave-wave interactions within high-β collisionless plasmas.
In this paper, we focus on the dynamical and thermodynamical interactions between two waves: either ion acoustic waves (IAWs) and AWs, or two co-propagating AWs.In §2 we outline expectations for these interactions based off of analytical theory performed within the CGL-MHD framework with the help of simplified Landau-fluid heat fluxes (Chew et al. 1956;Snyder et al. 1997).Following this, in §3 we test these predictions using a finite-volume Riemann-based code that numerically solves the CGL-MHD equations using the Landau-fluid heat fluxes.Finally, §4 lays out the implications our results have for both weak and strong turbulence, and details some of the barriers that remain to establish a complete theory of Alfvénic turbulence in high-β collisionless plasmas.

Theory
The degrees of freedom associated with wave-wave interactions in a collisionless plasma are numerous, and so to focus our investigation we make several simplifying assumptions.First, we take the electrons to be much colder than the ions, so that we may set their temperature T e = 0. Certain high-β astrophysical plasmas are thought to produce these circumstances, such as radiatively inefficient accretion flows onto supermassive black holes (e.g., Quataert 2003).However, the results we present are easily generalizeable to nonzero electron temperature (see §4).Second, our focus is on the dynamics that occurs on wavelengths much larger than the ion Larmor radius and at frequencies much smaller than the ion gyrofrequency.Accordingly, we employ a fluid model that accounts for the adiabatic production of pressure anisotropy, its spatial redistribution by field-aligned heat fluxes, and its dynamical feedback on the plasma through the action of an anisotropic pressure, all while neglecting finite-Larmor-radius effects.These 'CGL-MHD' equations (Chew et al. 1956) are: Our notation is standard: ρ is the mass density, u is the bulk fluid velocity, B is the magnetic field, and p ⊥ and p ∥ are the components of the pressure tensor oriented perpendicular and parallel to the local magnetic-field direction b = B/B.Note that (2.1a) and (2.1c) may be used in (2.1d) and (2.1e) to obtain where d/dt .= ∂/∂t + u • ∇ is the comoving time derivative.These equations express conservation of the first and second adiabatic invariants in the frame of the flow, but for the conservative redistribution of p ⊥ and p ∥ by the field-aligned conductive flows of perpendicular and parallel heat, q ⊥ and q ∥ , respectively.Our third assumption is that these heat fluxes are given by the '3+1 model' of Snyder et al. (1997), which is tailored to capture the linear collisionless damping rates of small-amplitude fluctuations: where c s,∥ = p ∥ /ρ is the parallel sound speed.The parameter |k ∥ | is a characteristic field-aligned wavenumber for the fluctuations in B, p ⊥/∥ , and ρ (e.g., Sharma et al. 2006); it is a simplification of the magnitude of the magnetic-field-aligned gradient operator, which would otherwise be very costly to calculate exactly because it is defined with respect to the local magnetic field and involves global operations.Because these heat fluxes describe linear collisionless damping, inherent in their use is the assumption that any perturbations produced by the waves studied here are sufficiently small with respect to the background plasma that nonlinear effects, such as saturation of the collisionless damping, can be ignored.Accordingly, while their analytical simplicity is invaluable in closing the CGL-MHD system, we discuss in §2.3 some circumstances in which these heat fluxes become inaccurate.

Formulation of the problem
Armed with (2.1), we now constrain the waves whose interactions are the focus of this paper.First, we consider only waves that propagate parallel to the background magnetic field B 0 , such that k ⊥ = 0.This choice removes the nonlinear Alfvénic interaction that normally drives a cascade to small scales, meaning that any interaction seen is solely a consequence of the effects of pressure anisotropy discussed in this work.Admittedly, many problems of interest involve fluctuations that possess k ⊥ /k ∥ of order unity or even much greater, such as critically balanced Alfvénic turbulence (Goldreich & Sridhar 1995).However, as we demonstrate in §2.2, one of the most important aspects of these wave interactions concerns the timescales associated with each interacting wave.Because AWs and IAWs possess frequencies that are independent of k ⊥ , this facet of their interaction is captured regardless of k ⊥ (further discussion of the k ⊥ = 0 assumption can be found within §4).Second, we assume that no constant background pressure anisotropy ∆ 0 is present in the plasma.A non-zero background pressure anisotropy may be included with little effect to the results (so long as the plasma remains stable) by replacing the Alfvén speed v A .= B 0 / √ 4πρ 0 , where ρ 0 is the (assumed uniform) background density, with an effective Alfvén speed v A,eff .= v A 1 + β∆ 0 /2.Our chief interest in this paper is how an AW responds to the fluctuating pressure anisotropy that results from other waves -a channel of communication that is customarily ignored in studies of plasma turbulence.Therefore, to distill this physics from (2.1), we pursue an asymptotic ordering of the wave perturbations that isolates the Alfvénic response to an externally generated pressure anisotropy ∆ = ∆(t, r).First, the AWs subject to ∆ are taken to be small enough in amplitude that their own nonlinearities can be ignored on the timescales of interest here, viz.δB ⊥ /B 0 ∼ u ⊥ /v A ≪ 1.With k ⊥ = 0, the vectors δB ⊥ and u ⊥ can be oriented along any perpendicular coordinate without loss of generality, and so henceforth we refer only to δB ⊥ and u ⊥ .Second, we take ∆ to be large enough that its effect on the Alfvénic dynamics cannot be neglected, i.e. |β∆| ∼ 1.This ordering allows the stress associated with the pressure anisotropy, (p ⊥ −p ∥ ) bb , to be comparable to the Maxwell stress, BB/4π.Lastly, we allow for compressive fluctuations in the parallel velocity (u ∥ ) and density (δρ).Because k ⊥ = 0 implies δB ∥ = 0, these fluctuations are purely acoustic, propagating at a speed ∼c s that is much faster than v A .These assertions combined, the resulting maximal ordering is given by To complete the ordering, we also make assumptions regarding the time evolution of these quantities.In particular, we assume that the time derivatives of δB ⊥ and u ⊥ are Alfvénic (∂ t ∼ k ∥ v A ), while all other (compressive) quantities evolve on sonic timescales Without yet discussing the production/evolution of ∆ (to be addressed in §2.2 and §2.3), we may then derive the evolution equations for our ordered Alfvénic fluctuations.Applying (2.4) to the components of (2.1b) and (2.1c) that are perpendicular to the guide field, and retaining only the leading-order terms in ϵ, we find that (2.5b) The transparency of these equations can be improved somewhat by rewriting them in terms of Elsässer variables z ± = u ⊥ ± δB ⊥ / √ 4πρ 0 (Elsässer 1950): (2.6) The left-hand side of these Elsässer equations describe forward (−) and backward (+) propagating AWs.The right-hand side provides the nonlinear interaction with another fluctuation that perturbs the pressure anisotropically.This anisotropy can arise from linear acoustic modes, or even other large-amplitude AWs that generate ∆ nonlinearly.
We analyze both of these situations in sections 2.2 and 2.3 respectively.

Interaction between ion-acoustic and Alfvén waves
In this section, we describe the interaction between monochromatic IAWs and AWs in an otherwise uniform plasma with density ρ 0 and pressure p 0 = βB 2 0 /8π.This demands an understanding of the evolution of the ∆ that appears on the right-hand side of (2.6).In our high-β ordering (2.4), the IAW fluctuations are energetically dominant over the Alfvénic fluctuations and, given their compressive nature, are chiefly responsible for the pressure perturbations δp ∥ and δp ⊥ .This can be seen by applying (2.4) to the compressive components of (2.1), recalling that the time derivatives of compressive quantities are ordered as ∂ t ∼ kc s , and retaining only leading-order terms: (2.7d ) These equations describing ion-acoustic fluctuations are entirely linear, with no feedback from the Alfvénic fluctuations.Therefore, we may simply prescribe the pressure anisotropy in (2.6) being due to ion-acoustic fluctuations as where ω C ≈ (±1.47 + 0.46i)k C c s is the (complex) linear frequency of a damped ionacoustic wave.4Substituting this anisotropy (with +ω C to start) into the Elsässer equations (2.6) and Fourier transforming into k-space yields the following evolution equation for AWs at each wavenumber: (2.9) The right-hand side of (2.9) shows that both forward-and backward-propagating AWs are generated at different wavenumbers when an AW is subject to the pressure anisotropy of an acoustic wave.Although this interaction term is not precisely proportional to the , no change to either z + or z − will result if their initial amplitudes are exactly equal.
Assuming only an initial z − , which makes the acoustic and Alfvén waves counterpropagating, the dynamics at the initial Alfvén mode's wavenumber k A are simply (2.10) The leading-order wave-wave interaction can be found at k A + k C by evaluating (2.9) given (2.10): where we have defined , and thus a weak interaction in high-β plasmas.To obtain a stronger interaction, we focus on the amplitude of the reflected wave packet (for Re(ω C ) > 0) and maximize the absolute value of the z + coefficient.Doing so indicates that the strongest interaction occurs when where the subscripts 'r' and 'i' denote the real and imaginary parts of ω C .With β ≫ 1, equation (2.12) makes the absolute value of the z − coefficient ≈0.3z 0 β∆ 0 .For β∆ 0 ∼ 1, this is always an O(1) interaction, decreasing linearly with the amplitude of the acoustic mode.Perhaps unsurprisingly, the strongest interaction occurs when the frequencies are nearly matched, although this also means that there is effectively no change in wavenumber of the Alfvénic fluctuations, since In this sense, the strongest IAW-AW interaction resembles the frequency matching condition of the parametric decay instability (Sagdeev & Galeev 1969).Additionally, while the maximum z + coefficient occurs for ω C,r > 0, the situation with k C ∼ k A / √ β still yields an O(1) reflection if ω C,r < 0 (i.e., if the waves are co-propagating).
In order to apply these conclusions to physical systems, we must also consider the case where ∆ is purely real.This demands that the Fourier expansion of ∆ contains both +ik C and −ik C exponentials, which in turn create source terms on the right-hand side of (2.9) from both k − k C and k + k C .While this complicates the problem somewhat, equation (2.11) can still be used to obtain an order-of-magnitude estimate for the solution.The reason for this stems from the fact that the acoustic mode collisionlessly damps on the same timescale as it propagates.By the time an Alfvénic fluctuation at k A + k C grows because of the interaction to a similar amplitude as that of the original wave at k A , ∆ has largely vanished.The subsequent interaction of the k A + k C AW with the now nearly depleted ∆ to generate a k A + 2k C mode will be considerably weaker.As a result, the amplitude of the k A + 2k C fluctuation is likely to be small, and may not provide significant feedback to the evolution of the mode at k A + k C .This prediction is tested using numerical simulations in §3.
The conservation of certain quantities by this interaction may also be of interest.In particular, it may appear jarring that there is a nonlinear mechanism in (2.6) by which new Alfvénic fluctuations can be generated, but no corresponding sink present in (2.7).This would seem to imply that energy is not conserved if |z ± | 2 increase but nothing happens to the acoustic fluctuations.In this case, the apparent mismatch stems from the fact that, if u ∥ /c s ∼ u ⊥ /v A , then there is considerably more (∼β×) energy present in the acoustic fluctuations than in the Alfvénic ones.As a result, any energy given up by the acoustic fluctuations is higher order to themselves, but leading order to an AW with comparable amplitude.This can be demonstrated by going to third order in the compressive equations and re-establishing conservation of energy, a calculation presented in detail in Appendix A. 5 Why energy goes to the Alfvénic from the acoustic fluctuations is not as obvious.The reason becomes apparent however, when considering the same interaction without damping of the acoustic wave.With ω C,i = 0 in (2.11), the leading-order Alfvénic response would then have a purely periodic amplitude at each k.Essentially, the effects of a fluctuation with ∆ > 0 would be nearly wiped out by the subsequent equivalent negative anisotropy.Only the root-mean-square amplitude of the newly generated fluctuations would grow, and on a timescale much slower than the AW propagation.Instead, if the anisotropy decays, then the positive anisotropy generated during the first half of the acoustic period would be larger than the negative anisotropy generated during the second half-period, meaning that the difference in the nonlinear effects accumulates.Once the acoustic wave is decayed away, each k of the Alfvénic spectrum is left with only a steady-state amplitude and energy.

Interaction between two AWs
The other source of pressure anisotropy we consider for (2.6) is another, largeramplitude AW.In particular, we are interested in the ability of nonlinear AWs to generate pressure anisotropy, thereby modifying their own v A,eff in a manner that is dependent upon their amplitudes (similarly to the interruption process described in §1; Squire et al. 2016).We therefore choose to study the interaction between two co-propagating AW packets, which are initially isolated such that their self-modified v A,eff 's do not affect the other wavepacket.In this setup, one of the wave packets will be large enough in amplitude such that it generates a pressure anisotropy ∆ ∼ β −1 , thereby fulfilling the role of the IAW in §2.2.This will be called the 'primary' wave packet (with relevant quantities denoted by the subscript 'p').We first calculate the pressure anisotropy generated by this primary wave packet and demonstrate that its v A,eff is increased as a result.Given sufficient time then, the two initially separated wave packets will be able to interact if the trailing wave packet is of larger amplitude.To study the subsequent interaction, we consider the leading ('secondary') wave packet to satisfy the Alfvénic components of (2.4), thereby allowing us to employ (2.6) in solving for its nonlinear deformation, while ignoring the back reaction on the primary AW packet.
The mechanism for producing pressure anisotropy differs considerably between IAWs and shear-Alfvén waves.While pressure anisotropy is produced at linear order for the former, the latter change neither the magnetic-field strength nor the plasma density to linear order in the fluctuation amplitude, and so by adiabatic invariance have no associated linear pressure anisotropy.Indeed, equations (2.2) without the heat fluxes give6 d dt (2.13b) The magnetic-field strength changes only at second order in the fluctuation amplitude, according to B = B 2 0 + δB 2 ⊥,p (while δρ = 0).From this, it is clear that δB 2 ⊥,p /B 2 0 must be ≳1/β in order for the stress associated with the induced pressure anisotropy to be comparable to the Maxwell stress (the usual restoring force in an AW; Squire et al. 2016).
In AW packets of this amplitude, the addition of a competitive pressure-anisotropic stress dictates that their initial perturbation will be unable to propagate rigidly at v A as is the case when ∆ = 0. Instead, a strictly positive (owing to δB 2 ⊥,p > 0) nonlinear pressure anisotropy will be generated by the oscillating δB ⊥,p with a dominant wavenumber of k = 2k A .This cannot, however, last in the presence of the heat fluxes (2.4).The pressure anisotropy ∆ p is generated by the AW packet as it propagates at a rate of ∼ kv A , while it is diffused by the heat fluxes at a much faster rate of ∼ kc s .As a result, ∆ p is smeared out across the wave packet as it propagates, inheriting the spatial structure of the wave packet envelope rather than the mode itself.The end result is that shear AW packets initialized with an amplitude δB ⊥ /B 0 ≳ β −1/2 generate smooth, envelope-scale, positive pressure anisotropies.Accordingly, v A,eff > v A and such packets propagate faster than those with small (linear) amplitudes. 7ith a ∆ p profile resembling the envelope of the parent AW packet (rather than the oscillating wave within the envelope), there will be a peak in the magnitude of ∆ p at the center of the packet, where δB 2 ⊥,p is the largest.Through v A,eff = v A 1 + β∆ p (x)/2, this region of the packet will travel faster than the regions with lower δB 2 ⊥,p that occur at the trailing edge of the envelope.In time, an isolated AW packet should therefore steepen to form a shock-like wavefront, which propagates at an enhanced v A,eff from the background. 8It is this shocked, more rapidly propagating wave packet that can interact with smaller amplitude AW packets in front of it.In a collisionless plasma, however, the heat fluxes driven by a shock-like gradient in the pressure cannot be accurately described by (2.3).This complicates our effort to solve (2.1d) and (2.1e) for the ∆ p of the steepened packet.In this case, two alternative descriptions of q ⊥/∥ are likely to be more accurate.For one, the heat fluxes may simply be that of free streaming particles.Unless reflected by the mirror force, plasma particles would transit the shock front and change their perpendicular temperature adiabatically according to δT ⊥,p /T 0 ≈ δB 2 ⊥,p /2B 2 0 .Another possibility, which is likely for near-interruption AW packets in particular, is that the distribution of particle velocities resulting from an adiabatic response to δB ⊥,p would be unstable to kinetic microinstabilities.Larmor-scale magnetic fluctuations generated by the ion-acoustic or whistler heat-flux instabilities (Komarov et al. 2018;Roberg-Clark et al. 2018) could scatter particles in pitch angle, giving the heat fluxes a similar functional form as (2.3), but with a different effective diffusion rate.We choose to continue our analysis under this assumption of diffusive heat fluxes, both because it may apply in particular to near-interruption AW packets, and for the fact that our '3+1' Landaufluid heat fluxes (2.3) are also diffusive, allowing us to verify the results numerically for arbitrary diffusion coefficients by modifying |k ∥ | (appendix B).In this case, given that particle motions are thermal, the upstream plasma with its higher temperature is able to pass through the shock front and diffuse downstream, at a rate controlled by the diffusion coefficient.This leads to the development of a pressure anisotropy precursor, which stretches in front of the shock as it propagates.For a finite diffusion coefficient, this precursor decays away steadily as one moves further downstream of the shock.We derive this precursor analytically in appendix B for a general diffusion coefficient κ, finding that the functional form of the pressure anisotropy is a decaying exponential in a frame moving with the shock.What remains to be discussed is how the ∆ p precursor of a large AW packet affects smaller-amplitude AW packets that cannot outrun its increased propagation speed, v A,eff > v A .As the tail of a secondary wave packet begins to feel the ∆ p of the largeramplitude primary packet approaching from behind it, there will be a larger v A,eff at the back of the secondary wave packet than at the front.This will cause the secondary packet to compress, shortening its parallel wavelength.Unlike the IAW-AW interaction however, no backward propagating AW packet will be created.To see why, consider the evolution equation (2.6) for z + assuming only an initial z − : (2.14) The final term on the right-hand side represents a possible source or sink, depending on its sign.However, for ∂ x ∆ p < 0, the first term on the right-hand side damps z + exponentially along the wave characteristics.Given that the ∆-precursor of a large-amplitude z − packet is monotonically decreasing in x, z + fluctuations are always damped in this interaction.This simplifies the equation for z − , which can be written in terms of its energy density (2.15) Solving (2.15) using the method of characteristics yields where x(t, x 0 ) is the trajectory of the x 0 characteristic determined by the solution of (2.17) and x 0 (t, x) is the foot of the characteristic obtained by inverting x(t, x 0 ).Equation (2.16) states that the tail of the smaller-amplitude AW packet will steepen according to d t x, while growing exponentially along characteristics (and more rapidly so for steeper ∆ p ). Motivated by the derivation of the ∆ p precursor for diffusive heat fluxes in appendix B, we consider an exponentially decaying anisotropy (2.18) which propagates rigidly at v A,eff = v A 1 + β∆ p,0 /2. Figure 1(a) shows the characteristics that such an anisotropy would generate if ∆ p,0 = β −1 and l = 0.1L, where L is a characteristic domain length.The compression of characteristics indicates that the trailing end of the secondary mode develops a shorter wavelength over time.Figure 1(b) exhibits the gain in energy experienced by a Gaussian wave packet that is subjected to this ∆ p (x) after t = L/v A , for different values of the decay length l ∆ /L.As l ∆ /L becomes smaller, the gain in energy of the mode begins to increase.In general this trend of increasing energy gain will continue, but here the gain curve drops when l ∆ /L becomes small enough that the packets becomes too well separated to feel the effect of ∆ within a time t = L/v A (for verification of this, we show the gain curve for two identical secondary packets initially separated from the primary by δx = 0.4L and δx = 0.8L).Without the interference of other waves in this duo, the secondary packet compression will likely continue until the anisotropy of the primary wave damps away, or the wavelength of the secondary packet reaches dissipation scales.In a scenario where the large-amplitude wave interacts with many smaller-amplitude fluctuations, the compression (and therefore energization) of the smaller modes should be considered in the damping of the largeramplitude mode, likely enhancing the rate of decay.9 We do not describe in any detail the interaction between primary and secondary antipropagating AWs, however there is a notable asymmetry in the ∆-mediated interactions between co-and anti-propagating AWs.Anti-propagating interactions are considerably weaker, both because of the much shorter interaction time caused by the larger difference in propagation velocity, and because of the difficulty of establishing any sort of frequency matching with an exponential pressure anisotropy.With co-propagating AW packets, the primary packet modifies v A,eff nonlinearly, thus the difference in velocities between packets is proportional to δB 2 ⊥ /B 2 0 ; in the anti-propagating case, the difference in velocities is ≈2v A .Furthermore, due to the monotonically decreasing v A,eff associated with the ∆-precursor of the primary packet, the primary packet carries the co-propagating secondary packet with it and the two can interact until microphysical effects interfere.Alternatively, for an anti-propagating interaction, if one were to calculate the k-space dependence of the primary packet's ∆, then the corresponding version of (2.9) could be solved, but instead with numerous k C because ∆ would not be monochromatic.That said, unlike in the IAW-AW interactions, the fact that ∆ is not monochromatic makes it difficult to achieve any semblance of frequency matching, and thus the interaction cannot similarly be tuned to maximize its strength.

Numerical simulations
To test the predictions made in section 2, we perform fluid simulations of each wavewave interaction within the CGL-MHD framework (2.1).Our numerical approach utilizes a new Riemann solver implemented in a version of the finite-volume simulation code Athena++ (Stone et al. 2008) that evolves p ⊥ and p ∥ using CGL equations that are closed with Landau-fluid heat fluxes (2.3) (Squire et al. 2023).Unless otherwise stated, the field-aligned wavenumber |k ∥ | in the Landau-fluid heat fluxes is set equal to the initial dominant wavenumber of the primary pressure-anisotropy-generating fluctuation, viz., the IAW in the acoustic-Alfvén interaction, and the larger-amplitude AW in the Alfvén-Alfvén interaction.All simulations are performed with β = 400, with only linear perturbations being initialized for the fluid moments of each wave: u ⊥ and δB ⊥ for AWs; and u ∥ , δρ, and δp ⊥/∥ for acoustic waves.The interaction of acoustic and Alfvénic fluctuations is studied primarily using monochromatic plane waves, while Alfvén-Alfvén interactions are studied using wave packets to demonstrate how a larger-amplitude AW packet can 'catch up' with an initially separated, smaller-amplitude AW packet.An acoustic-Alfvén interaction simulation was also performed using wave packets to visually depict the generation of a reflected AW as well as the non-locality of the interaction with respect to wavenumber.All monochromatic plane-wave simulations use periodic boundary conditions, while all simulations involving wave packets make use of outflow boundary conditions (where the derivatives of all fluid fields are set to zero at the boundaries). 10The only spatial coordinate is x ∈ [0, 1], aligned with the background magnetic field and measured in dimensionless Alfvénic units (where v A = B 0 / √ 4πρ 0 = 1) such that the Alfvén crossing time of a domain of length L = 1 is L/v A = 1.AW packets are initialized according to δB ⊥ = αB 0 ψ(x) cos(kx) and u ⊥ = ±αv A ψ(x) cos(kx), with amplitude α and Gaussian envelopes ψ(x) having standard deviations of approximately 2π/k.For the AW-AW interaction, the primary packet amplitude is α = 4/ √ β, while the secondary packet amplitude is α = 1/β.In the IAW-AW interaction the acoustic wave is initialized with amplitude δu ∥ /v th = 2/3β, and the Alfvén wave is initialized with δB ⊥ /B 0 = δu ⊥ /v A = 1/2β.Because strong IAW-AW interactions at high β Figure 2: Landau-fluid CGL-MHD simulation of the interaction between the pressure anisotropy driven by an IAW packet (black) and an AW packet (orange) with matched frequencies.The IAW packet collisionlessly damps rapidly from the initial conditions (top frame), before nonlinearly deforming the short-wavelength AW packet (middle).After the IAW has decayed and propagated rightwards out of the domain, a backward-propagating AW packet has developed, with no apparent change to the wavenumber of the parent AW.
are very non-local in k space, AW fluctuations with wavelengths of as little as onehundredth of the domain length need to be resolved.Similarly, for an AW packet to remain well separated from the outflow boundaries but still propagate for several linear times, the wave packet envelope and its dominant wavelength must be much smaller than the domain length.To capture these high-k fluctuations adequately, a resolution of 9,216 cells is used for the AW-AW interactions and up to 12,288 cells for the IAW-AW interactions.Importantly, the amplitudes of all fluctuations are chosen such that their associated pressure anisotropies are within the bounds of the firehose and mirror instabilities, −2 ≲ β∆ ≲ 1 (Rudakov & Sagdeev 1958;Southwood & Kivelson 1993).

Interaction between ion-acoustic and Alfvén waves
To illustrate the problem being studied, we begin with a simple example of copropagating, quasi-monochromatic acoustic-and Alfvén-wave packets.We set the dominant wavenumber of the acoustic wave to k C = k A / √ β, such that the modes are nearly matched in their linear frequencies.Figure 2 shows the evolution of these two packets from the initial setup until the point when the IAW has passed out of the domain, leaving behind only the modified AW packet(s).The top panel, showing the initial conditions, demonstrates that the initial pressure anisotropy associated with the IAW (black) dips as low as β∆ ≈ −2, and does not overlap with the AW to any significant degree.As the acoustic wave propagates, it decays rapidly at a rate of ≃0.31ω C .Therefore, by the time it overlaps with the AW, the amplitude of its driven pressure anisotropy has decayed significantly (middle panel).Nonetheless, the modification of the local v A,eff can be seen in the tail of the AW as giving rise to a new, smaller-amplitude fluctuation.By the time the IAW has decayed away and propagated out of the simulation domain, two changes Figure 3: Energy contained within backward-propagating AWs after one Alfvén time of interaction between a forward-propagating monochromatic AW and a standing monochromatic IAW, determined numerically as a function of the ratio of the Alfvén and acoustic wavenumbers using Landau-fluid simulations (red pluses).An approximate analytical solution for the energy (3.1) is shown as the solid line.Overall good agreement between theory and simulation is found, demonstrating strong interaction when the wave frequencies are approximately matched.
have occurred to the initial wave packet.Most obvious is that, as predicted by (2.11), a smaller-amplitude, backward-propagating AW packet of roughly equal wavelength to the parent AW packet has been generated.This process therefore appears to mimic some features of the parametric decay instability of AWs.It distinguishes itself from parametric decay, however, by the more subtle increase in the amplitude of the initial forward-propagating AW packet.Instead of reflecting a portion of the initial wave packet and decreasing its amplitude, the IAW spawns these modifications to the AW packet by giving up an infinitesimal (with respect to itself, not the AW) portion of its own energy.In the end, the degree of Alfvénic imbalance, normalized by the total energy of the AW fluctuations, decreases by roughly 5%.
The interaction between the wave packets in figure 2 is expected to peak when the frequencies of the IAW and AW are roughly matched.However, the exact amplitude of the reflected fluctuation will depend on the minutia of the initial conditions, such as the packets' initial separation and the widths of their envelopes.We therefore move to the periodic plane-wave setup used to derive (2.9), so that we can assess the accuracy of (2.11) with fewer variable factors.In doing so, we choose to initialize a standing IAW with zero initial pressure anisotropy (otherwise the intensity of the interaction would be hidden by an immediate Alfvénic response to the nonzero anisotropy at t = 0).This is accomplished by perturbing u ∥ (x, t = 0) = αv th sin(k C x) alone, with an amplitude α = 2/3β such that the corresponding 'forced' anisotropy in (2.6) would be given by ∆(x) = β −1 sin(ω C,r t) cos(k C x) exp(−ω C,i t).We simulate the interaction for a range of k A /k C , all for an elapsed time of t f = 2πω −1 A .The energy E − = |z − | 2 contained within all z − fluctuations at the end of each simulation is plotted versus k A /k C in figure 3, compared against an analytic estimate based off of (2.11).The analytic estimate is simply obtained by accounting only for the initial AW within the source term on the right-hand side of (2.9).Separating the expression for the pressure anisotropy generated by the standing This explains the accuracy of (3.1), as fluctuations at k A ± 2k C are too weak to affect the solution dramatically.In panel (b), the change in energy imbalance is several orders of magnitude smaller than the increase in total energy of the Alfvénic fluctuations, demonstrating that the IAW-AW interaction decreases imbalance with respect to overall AW energy.
wave into its complex exponential parts, the approximation is given by where the ±ω C in parentheses indicates the direction of IAW propagation used when evaluating z(k A +k C ).As expected, both the prediction and the analytic estimate peak in interaction strength when k A /k C ∼ √ β/2 ∼ 10, and agree rather well.The reason for this agreement can be found within figure 4(a), which shows the energy spectra of the forwardand backward-propagating Alfvénic fluctuations at t = 2π(5ω A ) −1 and t = 2π(ω A ) −1 .In figure 4(a), the slow, diffusive k-space nature of the IAW-AW interaction can be seen.With the initial AW represented by the peak in E − (k), each Fourier component of the energy is several orders of magnitude smaller than the initial k A , with the sole exception being k A ±k C (denoted by two vertical dotted lines).Furthermore, the difference between E(k A ± k C ) at t = 2π(5ω A ) −1 and t = 2π(ω A ) −1 is much smaller than that for the higher/smaller k.Indeed, it is possible to show that, when energy cascades in both directions and k C /k A is infinitesimal, equation (2.9) exhibits diffusive-like behaviour in k-space, with a diffusion rate that is proportional to k 2 C .Therefore, our estimate (3.1) is supported by this observation of k A ± k C being the dominant newly generated wavenumbers within a single Alfvén time.The E ± additionally evolve nearly identically, demonstrated by the relatively small change in cross-helicity with respect to the total energy shown in figure 4(b).This is expected, given that the form of (2.9) suggests this interaction relates to the degree of imbalance at each wavenumber, and appears to effectively damp it over time.Figure 5: The magnetic field and pressure anisotropy profiles of the steepened AW packet (panel a), alongside the relationship between the perpendicular pressure and the magnetic perturbation amplitude at the front of the shock (panel b).The steepened wave packet is led by a shock in the magnetic-field profile, with a smoothed pressure anisotropy profile and ∆ precursor modifying the local Alfvén speed.The dip in ∆ behind the wave packet results from the initial conditions of the AW not including ∆, hence the magnetic-field strength decreases in this region as the packet propagates away.As this packet propagates, the approximate conservation of the T ⊥ /B adiabat sets the amplitude of the decaying anisotropy, and thus v A,eff , at the shock front.

Interaction between two AWs
The interaction of co-propagating AWs relies upon the differing propagation speeds that result when different amplitudes generate different local pressure anisotropies.Accordingly, all simulations presented in this section will be of initially isolated wave packets with outflow boundary conditions.Both wave packets are initialized with a dominant wavenumber of k A = 80π, with amplitudes of δB z /B 0 = −δu z /v A = 4/ √ β and δB y /B 0 = −δu y /v A = β −1 for the large-amplitude and small-amplitude waves, respectively.The magnetic perturbations have been assigned to different directions only for the purpose of visualizing them separately; this choice is not necessary for the interaction itself.The standard deviations of the Gaussian wave envelopes are 2π/k A for both packets, the centres of which are initially separated by a distance of ∆x = 0.11.
Before examining the interaction between these two wave packets, we first verify our predictions about the evolution of the primary wave packet.Figure 5(a) displays the primary AW packet after it has steepened to form a shock.As in Squire et al. (2016), the waveform (orange) has become square-like in order to minimize the change in δB 2 ⊥,p .Due to the rapid action of the heat fluxes, the pressure anisotropy (black) has assumed a shape similar to that of the wave envelope, with an additional precursor that extends in front of the magnetic perturbation.A negative dip in β∆ is present behind the tail of the wave packet, although it does not propagate with the packet; it actually results from the magnetic perturbation departing from its initial location, causing a net decrease in |B| at the location where it was initialized.This localized dip is akin to the negative anisotropy generated during the more conventional interruption process of a monochromatic AW (Squire et al. 2017a), as we did not initialize the wave packet with the pressure anisotropy that it generates shortly thereafter.Interruption does not occur here because the β∆ dip is not sufficiently negative to nullify the magnetic tension.The pressure anisotropy of the primary packet is evolved alongside the magnetic perturbation of the secondary packet.The packets are initially separated (top), yet over time the enhanced effective speed allow the primary packet to approach the secondary from behind (middle).By the end of the simulation, the effective wavenumber of the secondary packet has more than doubled from the steepening induced by the primary packet's ∆ precursor.magnetic-field strength.After an initial adjustment into the near-steady-state shocked AW packet, δT ⊥ closely follows the peak value of δ|B| = δB 2 ⊥,p at the shock front.This supports our prediction made in appendix B that the maximum value of the pressure anisotropy can be calculated via conservation of the double adiabats.
In §2.3, we predicted that the pressure-anisotropy-enhanced propagation speed of largeamplitude (primary) AW packets would allow them to 'catch up' to and then distort smaller-amplitude (secondary) AW packets.Evidence of this physics can be found in figure 6, which shows how the pressure anisotropy generated by the primary (black) interacts with the magnetic-field perturbation of the secondary (orange).In the top panel, shortly after the magnetic disturbance of the large-amplitude mode begins to propagate, there is a short-wavelength ∆.This pressure anisotropy is smoothed out rapidly by the action of heat fluxes, however, and leads to a smoothly increasing positive anisotropy perturbation whose front propagates at a slightly enhanced v A,eff > v A .The second panel at k A v A t = 9.9 illustrates the moment that this front catches up to the tail end of the smaller-amplitude AW; the overlap between the two fluctuations causes this tail to propagate slightly faster, thereby compressing the waveform of the secondary.As this process continues into the final frame, the small-amplitude wave packet has had its width reduced by roughly half, and its amplitude increased due to the growth of E − (x) (see (2.16)).Because the difference in Alfvén speeds between these two packets is relatively small, the cumulative deformation occurs slowly with respect to the Alfvén time of one of the packets, making this a weak interaction even though it leads to such a dramatic change in the structure of the wave packet envelope.This compression is expected to continue until the anisotropy of the large-amplitude packet has nonlinearly damped to be sufficiently small, or it is deformed by interaction with other wave packets.

Discussion: applications in turbulence
The analysis presented in § §2 and 3 suggests that modern theories of turbulence in highβ collisionless plasmas must account for fundamentally different wave-wave interactions.In weak MHD turbulence, energy and momentum conservation dictate that AWs are unable to modify their k ∥ through interaction (Montgomery & Turner 1981;Galtier et al. 2000).Yet, we have just demonstrated a weak interaction of two AWs with differing amplitudes where the parallel wavelength does in fact change significantly.In a similar vein, the compressive cascade of MHD turbulence relies on the fact that slow modes possess a linear timescale that never exceeds the Alfvénic timescale of the turbulence that mixes them.Further, these slow modes have no effect on the Alfvénic cascade and effectively decouple from it below the forcing scale (Lithwick & Goldreich 2001;Cho & Lazarian 2003;Schekochihin et al. 2009;Nazarenko & Schekochihin 2011;Howes & Nielson 2013).However, when slow modes are replaced by collisionless IAWs, these roles nearly reverse, with acoustic modes propagating too fast to be mixed effectively by the AWs, and in some cases even demonstrating an ability to reflect Alfvénic fluctuations.
Still, certain aspects of these waves that may be important in a turbulence setting are not considered in this study.In particular, because we focused solely on wave amplitudes for which the induced pressure anisotropies satisfy |∆| ≲ 1/β, we do not consider the effects of the firehose or mirror instabilities.Previous work has shown that AWs, IAWs, and other compressive fluctuations evolve quite differently if their self-generated pressure anisotropies are large enough in magnitude to excite these Larmor-scale instabilities (Squire et al. 2017a;Kunz et al. 2020;Majeski et al. 2023).As demonstrated from first principles by Arzamasskiy et al. (2023), the anomalous scattering that results from these instabilities is likely to influence the entire turbulent cascade, yielding yet another source of non-local energy transfer.
As mentioned in §2, this study was performed assuming cold electrons.However, generalizing to electrons with finite temperature is quite straightforward, as long as the electrons remain barotropic.With an electron pressure tensor that satisfies p e = p e I = (ρT e /m i )I, no additional pressure anisotropy is produced by either wave, and so AW propagation is entirely unaffected.The only modification is to the ratio of the real and complex parts of the IAW frequency ω C,i /ω C,r , which, unless the electrons partake in Landau damping of the acoustic mode, decreases with increasing T e .This would bolster the effectiveness of the IAW-AW interaction by extending the lifetime of the pressure anisotropy and allowing the interaction to occur over more than one linear time of the frequency-matched modes.
In §2.1, we noted the lack of k ⊥ as a key difference between these pressure-anisotropyinduced interactions and those found within standard Alfvénic turbulence.The nature of the pressure-anisotropic stress we highlight is significantly different from the Reynolds and Maxwell stresses, which dominate energy transfer in Alfvénic turbulence (Grete et al. 2017).Most importantly, β∆, being a scalar, will always affect v A,eff regardless of the ratio k ⊥ /k ∥ , doing so through the same nonlinearity that we retain in (2.6).Therefore, the effects of IAWs or AWs with β∆ ∼ 1 are expected to remain leading order when k ⊥ ̸ = 0. Furthermore, the linear timescales of IAWs and AWs are independent of k ⊥ , thus the frequency matching criterion of their interaction, and its spectral non-locality, are unlikely to change.
Given the further complications of microinstabilities and other high-β effects such as magneto-immutability (Squire et al. 2023), formulating an analytical model for turbulence in this parameter regime remains a challenging task.The intention of this work is primarily to hint at the physics underlying changes in the observed behaviour of collisionless turbulence at β ≫ 1, such as non-local energy transfer due to the pressureanisotropic stress (Arzamasskiy et al. 2023).To that end, we may still speculate upon some of the more qualitative consequences of the wave-wave interactions described herein.One of the most effective tools for describing strong turbulence is an appropriately reduced set of equations.For example, reduced MHD incorporates the concept of critical balance directly by ordering the linear frequency comparable to the nonlinear frequency, (Goldreich & Sridhar 1995;Schekochihin et al. 2009).Reinstating nonzero k ⊥ then, we can consider the reduced Alfvénic Elsässer equations that result from our 1D high-β ordering (2.4), modified to include k ∥ /k ⊥ ∼ δB ∥ /B 0 ∼ ϵ: where P total represents the combined perpendicular and magnetic pressures.Contained within the right-hand side of (4.1) are the various new effects discussed throughout this work that are absent in a collisional MHD model, or even a β ∼ O(1) gyrokinetic model.The first term on the right-hand side incorporates the ability of AWs to travel at different v A,eff , while both terms describe the deformation in z ± that results from wavewave interactions.Equation (4.1) does not distinguish between the source of β∆; be it generated by AWs with δB ⊥ /B 0 ∼ β −1/2 or IAWs with δρ/ρ ∼ β −1 , the Alfvénic cascade will be modified.In each case, however, the description of the source of fluctuations can be more difficult.For example, IAWs act non-locally to affect Alfvénic fluctuations in k-space, but the reverse may not be true for Alfvénic mixing of acoustic waves.In that case, there is not necessarily a single acoustic Elsässer equation that can be written at each scale, because multiple k ⊥ and k ∥ are inherent to the problem.Furthermore, in a collisionless turbulent cascade, β∆ cannot remain O(1) at every wavenumber throughout the inertial range.As all of the other quantities cascade, β will remain constant, so it is likely that interactions mediated by ∆ are a leading-order effect only in a specific portion of the inertial range where β∆ ∼ 1.One subtle conclusion of the above reduced equations is that pressure balance in these plasmas should be struck primarily between B 2 and the perpendicular pressure p ⊥ , rather than the isotropic pressure (Squire et al. 2023).The leading order of (4.1) implies that δp ⊥ /p 0 ∼ β −1 δB ∥ /B 0 ∼ ϵ 2 .In essence, the anisotropy generated should be dominated by the parallel pressure perturbation.This was found empirically in simulations of incompressibly driven turbulence using the same CGL-MHD Athena++ code employed here (Squire et al. 2023), as well as in the hybrid-kinetic simulations of Arzamasskiy et al. (2023).Additional effects may be uncovered by studying the weak turbulence resulting from these wave-wave interactions.Isolated wave packets are not typically considered within strong turbulence, however the long deformation times associated with weak turbulence allows waves to remain correlated over many local Alfvén-crossing times (Iroshnikov 1964;Kraichnan 1968).This also expands upon the range of scales over which β∆ exerts an influence by not requiring that the interaction remain strong.One effect that may be impactful to weak turbulence is the ability of AWs to modify each others' k ∥ .Current understandings of weak AW interactions involve resonance conditions describing threewave interactions (Galtier et al. 2000): terms to demonstrate that energy conservation is reestablished by including all terms of order ϵ 3 or higher.The total energy normalized to ρ 0 v 2 th,i is given by Already, it is clear from (A 1) that, according to the ordering (2.4), the Alfvénic contributions to the energy (terms involving u ⊥ and δB ⊥ ) are at most of order ϵ 3 , while the compressive contributions (the last three terms) are up to linear in ϵ.Taking the time derivative of (A 1), the rate of change of the total energy is given by We examine these terms in turn with the aim of establishing the lowest order at which they appear.First, we focus on those terms relating directly to the energy carried by the Alfvénic fluctuations, i.e., the top line of (A 2).Denote the characteristic Alfvénic and compressive wavenumbers by k A and k C , respectively.The first term is of order k C c s ϵ 4 (due to the continuity equation's compressive nature in ∂ t ), while the second and third terms are of order k A v A ϵ 3 .For all circumstances studied in this paper then (including frequencymatched fluctuations), the leading order of the second and third terms will dominate and be the main mechanism by which the Alfvénic fluctuations extract energy from the compressive ones.Therefore, neglecting the first term and inserting (2.5) for ∂ t u ⊥ and ∂ t δB ⊥ , the energy rate of change simplifies to For ω A ≈ ω C , the remaining Alfvénic term on the first line of (A 3) is of the exact same order as the kinetic energy terms on the bottom line, while, for k A ≈ k C , it is one half-order in ϵ smaller.In order to satisfy energy conservation, we must then include all compressive terms that are up to, but not including, those of order ϵ 4 .Before embarking on the compressive terms however, we note that the heat fluxes need not be included in (2.1e) and (2.1d) for the proof of energy conservation between Alfvén and acoustic waves.Their effect is to diffuse wave energy into the background, and not to facilitate communication between the Alfvénic and compressive fluctuations.We therefore proceed by evolving the pressure perturbations according to the double adiabats (2.13), which can be combined in the following manner: Noting that δB ∝ δB 2 ⊥ , which makes B −1 ∂ x B second order in ϵ, the final two terms cancel up to and including O(ϵ 3 ).We also use the fact that (A 4) will be integrated over space to swap the remaining x-derivative on the bottom line from p ⊥ /ρ to u ∥ ρ, picking up a sign change.Using the continuity equation to cancel the first two terms on the second line, no more of the second line remains.Normalization and integration of the first term on the right-hand side of (A 4) then gives dx 1 2 where (2.5a) has been used to substitute for ∂ t δB ⊥ .This term will cancel the first term on the right-hand side of (A 3), leaving us with the following expression for the energy rate of change: The fully nonlinear parallel momentum in this case is given by Multiplying (A 9) by u ∥ /ρ 0 v 2 th,i = u ∥ /2p 0 , the last two terms on its right-hand side become of order O(ϵ 4 ), while the first two cancel the remaining terms in (A 8).As a result, we have shown that, up to and including O(ϵ 3 ) terms, total energy is conserved by the interaction of Alfvén and ion-acoustic waves according to the ordering (2.4).This finding tells us that, even when the frequencies of the two modes are matched and k A ∼ √ βk C , the leading order for the dynamical equations describing each wave remains unchanged.Naturally, no terms neglected in deriving (2.5) are enhanced with respect to the leading order by letting k A increase by ϵ −1/2 .At the same time, Alfvénic feedback effects on the compressive fluctuations were of order ϵ 3 , meaning the leading order compressive equations (2.7) remain the same until ω A ∼ β 2 ω C .Nonlinear steepening will remain far more important in the frequency-matched regime of strong interaction, thus for (small) acoustic fluctuations of amplitude ∼β −1 , the assumption of linearity is quite robust.
Figure 7: Comparison of the analytically predicted and simulated decay lengths of the pressure anisotropy precursor of a steepened AW packet, for q ⊥/∥ given by the '3+1' heat fluxes (2.3).The calculated decay length λ T,⊥ is normalized by the dominant wavelength of the AW packet λ A .On the left, the decay length is varied with respect to the Landau wavenumber |k ∥ |, while on the right it is varied with respect to β. v A,eff evaluated at the shock front).Accordingly, the relation (B 3) becomes This gives the functional form of δT ⊥ leading up to the shock front (where δT ⊥,max is determined), and motivates our choice of exponentially decaying pressure anisotropy as a source for (2.6) used in §2.3.To test the accuracy of this estimate of the decay length of δT ⊥ , a series of Landau-fluid simulations were performed with varying choices of β and for |k ∥ | in (2.3), shown in figure 7. The setup of the wave packet in these runs is the same as that of the primary packet described in §3.2.The decay length of δT ⊥ is found in these simulations by waiting until the wave packet steepens to form a shock, then fitting an exponentially decaying function to the region directly in front of the peak of the magnetic-field perturbation.Compared with the theoretical estimate using the '3+1' heat fluxes, strong agreement is found.The actual magnitude of δT ⊥,max is determined by acknowledging that ∂ x δB 2 ⊥ is large in the region of the shock, and including its associated term at leading order in the equation left of the arrow in (B 4) provides Integrating this equation from the shock front to the distant downstream region where ∂ x δT ⊥ ≈ 0, we obtain the simple result that δT ⊥,max /T 0 ≈ (1/2)δB 2 ⊥,max /B 2 0 .Knowing that δp ∥ may be neglected, the same steps can be applied to (2.1e) to obtain a diffusion equation for δρ, which is found to decay on a length scale of (3/2)κ/v A,max instead of κ/v A,max .Similarly, δρ max is determined to be δρ max ≈ (1/3)δB 2 ⊥,max /B 2 0 .Combined, these yield the final result that ∆ max ≈ δp ⊥,max /p 0 ≈ (5/6)δB 2 ⊥,max /B 2 0 , and v A,max = v A 1 + β∆ max /2.

Figure 1 :
Figure 1: (a) Characteristic curves for an AW subjected to an exponentially decaying pressure anisotropy (2.18) with characteristic decay length l ∆ = 0.1L and phase speed v A,eff = √ 2v A .Characteristics originating near the left of the domain begin to converge by the time they reach the right side of the domain, indicating ∆-induced AW steepening.(b) Energy gain by the secondary AW packet after one Alfvén-crossing time for different l ∆ , demonstrating that a steeper ∆ profile leads to more rapid gain in energy by the secondary packet.

Figure 4 :
Figure 4: (a) Energy spectrum of forward-and backward-propagating Alfvénic fluctuations during interaction with a standing IAW.(b) Total change in wave energy (orange) and imbalance (black), normalized to the initial Alfvénic fluctuation energy E 0 , versus time.In panel (a), the IAW-AW interaction primarily generates new AW fluctuations at k A ± k C , exhibited by the steepness of the spectra outside of the interval k∈ [k A − k C , k A + k C ].This explains the accuracy of (3.1), as fluctuations at k A ± 2k C are too weak to affect the solution dramatically.In panel (b), the change in energy imbalance is several orders of magnitude smaller than the increase in total energy of the Alfvénic fluctuations, demonstrating that the IAW-AW interaction decreases imbalance with respect to overall AW energy.
Figure 5(b) demonstrates the relationship between the perpendicular temperature and the perturbation to the

Figure 6 :
Figure 6: Time slices of a Landau-fluid simulation of the the AW-AW packet interaction.The pressure anisotropy of the primary packet is evolved alongside the magnetic perturbation of the secondary packet.The packets are initially separated (top), yet over time the enhanced effective speed allow the primary packet to approach the secondary from behind (middle).By the end of the simulation, the effective wavenumber of the secondary packet has more than doubled from the steepening induced by the primary packet's ∆ precursor.
At this point, it is may become obvious that, with all of the remaining terms being compressive in nature, they must cancel because the nonlinear interaction of a wave with itself generally conserves energy in the absence of dissipative effects like heat fluxes or collisions.Regardless, to solidify the argument we will continue and demonstrate that there are no other Alfvénic feedback terms that appear at orders greater than ϵ 4 .Using the continuity equation to simplify the top line of (A 6) yields