Kinetic simulations of the interruption of large-amplitude shear-Alfv\'en waves in a high-beta plasma

TL;DR

This study uses hybrid-kinetic simulations to show that large-amplitude shear-Alfvén waves in high-beta collisionless plasmas are rapidly damped due to pressure anisotropy and firehose instabilities, preventing their support above a certain amplitude.

Contribution

It demonstrates the nonlinear interruption mechanism of shear-Alfvén waves in high-beta plasmas and highlights the role of velocity-space instabilities in plasma dynamics.

Findings

Large-amplitude waves decay rapidly above the amplitude limit.

Firehose modes generate long-lived fluctuations causing particle scattering.

Collisionless plasmas cannot support linearly polarized Alfvén waves beyond a certain amplitude.

Abstract

Using two-dimensional hybrid-kinetic simulations, we explore the nonlinear "interruption" of standing and traveling shear-Alfv\'en waves in collisionless plasmas. Interruption involves a self-generated pressure anisotropy removing the restoring force of a linearly polarized Alfv\'enic perturbation, and occurs for wave amplitudes $\delta B_{\perp}/B_{0}\gtrsim \beta^{\,-1/2}$ (where $\beta$ is the ratio of thermal to magnetic pressure). We use highly elongated domains to obtain maximal scale separation between the wave and the ion gyroscale. For standing waves above the amplitude limit, we find that the large-scale magnetic field of the wave decays rapidly. The dynamics are strongly affected by the excitation of oblique firehose modes, which transition into long-lived parallel fluctuations at the ion gyroscale and cause significant particle scattering. Traveling waves are damped more slowly, but are also influenced by small-scale parallel fluctuations created by the decay of firehose modes. Our results demonstrate that collisionless plasmas cannot support linearly polarized Alfv\'en waves above $\delta B_{\perp}/B_{0}\sim \beta^{\,-1/2}$. They also provide a vivid illustration of two key aspects of low-collisionality plasma dynamics: (i) the importance of velocity-space instabilities in regulating plasma dynamics at high $\beta$, and (ii) how nonlinear collisionless processes can transfer mechanical energy directly from the largest scales into thermal energy and microscale fluctuations, without the need for a scale-by-scale turbulent cascade.

Figures

Equations

$$\frac{\delta B_{\perp }}{B_0}\gtrsim {\beta }^{-1/2},$$

$$\frac{\partial f_i}{\partial t}+\mathbf{v}\cdot \frac{\partial f_i}{\partial \mathbf{x}}+\frac{q_i}{m_i}\left(\mathbf{E}+\frac{1}{c}\mathbf{v}\times \mathbf{B}\right)\cdot \frac{\partial f_i}{\partial \mathbf{v}}=0;$$

$$\mathbf{E}+\frac{1}{c}{\mathbf{u}}_i\times \mathbf{B}=-\frac{T_e\nabla n_i}{e{n}_i}+\frac{(\nabla \times \mathbf{B})\times \mathbf{B}}{4\pi q_i{n}_i}.$$

Full text

*Introduction.—*Shear-Alfvén (SA) fluctuations are fundamental to magnetized plasma dynamics Alfvén (1942); Cramer (2011); Ogilvie (2016). They are routinely observed in both laboratory and space plasmas Gekelman et al. (2011); Bruno and Carbone (2013), and are the basis for modern theories of magnetohydrodynamic (MHD) turbulence Goldreich and Sridhar (1995); Boldyrev (2006); Mallet and Schekochihin (2017). They are also uniquely robust among plasma waves, with large-scale linear dynamics that are nearly unmodified across both kinetic and fluid plasma models Cramer (2011).

The purpose of this Letter is to explore, using hybrid-kinetic simulations, an important exception to this robustness. We focus on linearly polarized large-scale SA waves above the “interruption limit” Squire et al. (2016, 2017),

[TABLE]

in a collisionless plasma. Here $\beta\equiv 8\pi p_{0}/B_{0}^{2}$ is the ratio of thermal to magnetic pressure, $B_{0}$ is a background magnetic field, and $\delta B_{\perp}$ is an Alfvénically polarized field perturbation. SA perturbations above the limit (1) rapidly transfer their mechanical energy from the largest scales to plasma microscales and thermal energy, without the help of a turbulent cascade. This paradigm is at odds with standard theories of Alfvénic turbulence in collisionless systems Schekochihin et al. (2009), and may be crucial for understanding turbulent energy dissipation in astrophysical plasmas ranging from the intracluster medium (ICM) Rosin et al. (2011); Peterson and Fabian (2006); Enßlin and Vogt (2006); Zhuravleva et al. (2014) to hot accretion flows Quataert (1998) and high-$\beta$ regions of the solar wind Bale et al. (2009); Bruno and Carbone (2013); Chen (2016); Bale et al. (2017).

The interruption of SA perturbations occurs due to the self-generation of pressure anisotropy, $\Delta p\equiv p_{\perp}-p_{\parallel}$ (where $p_{\perp}$ and $p_{\parallel}$ are the thermal pressures perpendicular and parallel to $\bm{B}$). Pressure anisotropy is created whenever $B=|\bm{B}|$ changes in a weakly collisional plasma. Further, if $\beta>1$, the anisotropic momentum stress $\nabla\cdot(\Delta p{\bm{B}}{\bm{B}}/B^{2})$ can be as important as, or even dominate over, the magnetic tension $\nabla\cdot(\bm{{B}}{\bm{B}})/4\pi$. This suggests that collisionless dynamics can differ from MHD predictions, even for large-scale perturbations satisfying $\lambda\gg\rho_{i}$, $\tau\gg\Omega_{i}^{-1}$ (where $\rho_{i}$ and $\Omega_{i}$ are the ion gyroradius and gyrofrequency, respectively).

“Interruption” occurs when a linearly polarized SA oscillation creates an anisotropy $\Delta p=-B^{2}/4\pi$, which offsets the magnetic tension and triggers the firehose instability on ion gyroscales Davidson and Völk (1968); Rosenbluth (1956); Yoon et al. (1993); Hellinger and Matsumoto (2000). Even at small wave amplitudes ($\beta^{-1/2}<\delta B_{\perp}/B_{0}\ll 1$), interruption is a nonlinear effect. We study this behavior using hybrid kinetics (kinetic ions, fluid electrons), in three velocity and two spatial dimensions (the latter is required to capture the 2-D oblique firehose instability, which is expected to evolve in similarly three spatial dimensions). We consider parallel standing and traveling SA waves in the large-scale-separation limit, $\lambda_{\mathrm{mfp}}>\lambda_{A}\gg\rho_{i}$ ($\lambda_{\mathrm{mfp}}$ and $\lambda_{A}$ are the ion mean-free path and SA wavelength), as relevant to many astrophysical systems (e.g., in the ICM $\lambda_{\mathrm{mfp}}\sim 10^{11}\rho_{i}$ Rosin et al. (2011)). Although the thresholds for the oblique and parallel firehose instabilities differ slightly Hellinger and Matsumoto (2000); Klein and Howes (2015), we organize our discussion around the latter ($\Delta p=-B^{2}/4\pi$) because of its importance for large-scale SA waves.

*Hybrid-kinetic method.—*By treating electrons as an isothermal massless fluid, the hybrid method removes electron kinetic scales, plasma oscillations, and light waves from the Vlasov-Maxwell equations, reducing simulation cost while retaining kinetic ion dynamics Byers et al. (1978); Hewett and Nielson (1978). The equations consist of (i) the Vlasov equation for the ion distribution function $f_{i}(\bm{x},\bm{v},t)$,

[TABLE]

(ii) Faraday’s law, ${\partial\bm{B}}/{\partial t}=-c\nabla\times\bm{E}$; and (iii), a generalized Ohm’s law,

[TABLE]

Here, $q_{i}$ and $m_{i}$ are the ion’s charge and mass, $\bm{E}$ is the electric field, $c$ is the speed of light, and $T_{e}$ is the electron temperature. The ion density $n_{i}(\bm{x})\equiv\int d\bm{v}\,f_{i}$ and bulk velocity $\bm{u}_{i}(\bm{x})\equiv\int d\bm{v}\,\bm{v}f_{i}$ are calculated from $f_{i}$, closing the system.

We use the second-order-accurate particle-in-cell (PIC) code, Pegasus (Kunz et al., 2014a). We employ the $\delta f$ method (Chen and Parker, 2003), which evolves $\delta f=f-f_{0}$ rather than $f$ itself, and take $f_{0}$ to be an isotropic Maxwellian. This reduces particle noise by $\sim\!(\delta f/f_{0})^{2}$, making it optimal for simulation of high-$\beta$ plasmas, where very small ($\ll 1/\beta$) deviations from a Maxwellian distribution must be accurately resolved.

*Simulation set up.—*We consider two initial conditions, which vary initially only on large scales. These are (i) a parallel standing SA wave initiated by a magnetic perturbation, and (ii) a parallel traveling SA wave. We focus on the standing wave because of its relevance to situations where $\langle dB/dt\rangle\neq 0$, e.g., Alfvénic turbulence ($\langle\cdot\rangle$ represents a spatial average). Although also important, we leave study of initial Alfvénic velocity perturbations to future work, due to the larger domains required to capture mirror instability dynamics Melville et al. (2016); Rincon et al. (2015). The initial ion distribution function is an isotropic Maxwellian with $T_{e}=T_{i}$, with a background magnetic field $\bm{B}=B_{0}\hat{\bm{x}}$ and $\beta_{i}=8\pi n_{i}T_{i}/B_{0}^{2}=100$. Our domains have width $L_{y}=50\rho_{i}$ and lengths up to $L_{x}=1000\rho_{i}$, to maximize scale separation between the SA wave and microscale dynamics. We use a spatial resolution of $\Delta x=0.3125\rho_{i}$ and $N_{\mathrm{ppc}}=4096$ particles per cell (ppc) for the two main simulations in this Letter. We initialize with $\lambda_{A}=L_{x}$ in the out-of-plane field, $\delta B_{z}=-\delta b\,B_{0}\cos(2\pi x/\lambda_{A})$, and, for the traveling wave, a corresponding velocity perturbation, $\delta u_{z}=\delta b\,v_{A}\cos(2\pi x/\lambda_{A})$. In each case, we take the wave amplitude $\delta b=0.5$, which is well above the interruption limit $\delta b_{\mathrm{max}}\approx 2\beta^{\,-1/2}$ Squire et al. (2016). Within the MHD model, these initial conditions would create continuing sinusoidal SA oscillations of period $\tau_{A}=2\pi/\omega_{A}=\sqrt{\beta_{i}}\lambda_{A}/\rho_{i}\Omega_{i}^{-1}$ (modified slightly by compressibility (Squire et al., 2017)).

Due to the wide range of time and space scales involved in this problem, careful numerical tests are crucial. In addition to previous Pegasus tests Kunz et al. (2014a), we tested the numerical parameters required to accurately propagate long-wavelength linear SA waves (with $\lambda_{A}/\rho_{i}=50$ to $1000$, $\delta b=0.05$). These tests demonstrated that high $N_{\mathrm{ppc}}\propto\lambda_{A}/\rho_{i}$ is required for large wavelengths, due to the build up of PIC noise over long simulation times. For production runs, $N_{\mathrm{ppc}}=4096$ was chosen based on these requirements. We also tested the convergence (with $N_{\mathrm{ppc}}$) of nonlinear standing waves at $\lambda_{A}/\rho_{i}=250$, and their dependence on $\lambda_{A}/\rho_{i}$ over the range $\lambda_{A}/\rho_{i}=125$ to $1000$. We observed broadly similar dynamics over this range.

*Shear-Alfvén standing wave.—*Figure 1 shows the spatiotemporal evolution of $\delta B_{z}$ for a standing SA wave with $\lambda_{A}/\rho_{i}=1000,\>\tau_{A}\Omega_{i}=10000$ ($\omega_{A}/\Omega_{i}=2\pi\times 10^{-4}$). The pictured snapshots are chosen to illustrate the distinct phases of nonlinear wave evolution. These are: (i) initial field decrease, which creates a negative anisotropy $\Delta p<-B^{2}/4\pi$, nullifying magnetic tension and triggering the firehose instability; (ii) eruption of oblique firehose modes Yoon et al. (1993); Hellinger and Matsumoto (2000); Kunz et al. (2014b) which push the wave back above $4\pi\Delta p/B^{2}=-1$; (iii) decay of oblique firehose modes into smaller-scale ($k_{\perp}=0$, $k_{\parallel}\rho_{i}\sim 1$) fluctuations that scatter particles and cause the large-scale $\delta B_{z}$ to decay; and (iv), dissolution of the wave into freely oscillating SA waves below the limit (1), which can oscillate freely. Of these stages, (iii) is notably different from the predictions of 1-D Landau-fluid (LF) models (Squire et al., 2016, 2017).

Figure 2 shows 1-D ($y$-averaged) wave profiles. Because of heat fluxes Squire et al. (2017), as $\Delta p$ decreases initially it is nearly homogenous in space [Fig. 2(a)]. This causes oblique firehose modes to erupt suddenly across the entire wave [Fig. 1(b)] at $t/\tau_{A}\approx 0.075$. These growing modes cause $\Delta p$ to increase Schekochihin et al. (2008) back into the stable regime $\Delta p>-B^{2}/4\pi$ by $t/\tau_{A}\approx 0.085$, where it stays until the SA wave decays.

The subsequent evolution of the oblique firehose modes controls the large-scale wave dynamics. If these (now stable) fluctuations scatter particles sufficiently strongly, $\delta B_{z}$ can decay with $\Delta p\approx-B^{2}/4\pi$; if they do not (e.g., if they are resonantly damped Foote and Kulsrud (1979); Quataert (1998); Melville et al. (2016)), $\delta B_{z}$ cannot decrease Squire et al. (2016). The firehose modes’ evolution is governed by $\Delta p$ Rosin et al. (2011), which varies in space. Near the wave nodes, where $S=|\nabla\bm{u}|\approx 0$ and $\delta B_{z}\approx 0$, $\Delta p$ is not driven by a large-scale $dB/dt$ and can freely decay Quest and Shapiro (1996); Schekochihin et al. (2008); Seough et al. (2015); Melville et al. (2016). Near the wave antinodes, where $S\sim\beta^{\,-1/2}\omega_{A}\approx 6\times 10^{-5}\Omega_{i}$ Squire et al. (2017) and $\delta B_{z}\neq 0$, $\Delta p$ is continuously driven by the decreasing field Matteini et al. (2006); Hellinger and Trávníček (2008); Kunz et al. (2014b); Melville et al. (2016); Riquelme et al. (2016).

Surprisingly, it is small-scale modes at the SA-wave nodes—the least firehose-unstable regions (with $4\pi\Delta p/B^{2}\approx-0.7$)—that cause the strongest particle scattering. This is illustrated by the background color in Fig. 2, which shows the effective ion collisionality $\nu_{c}/\omega_{A}$ as a function of space, measured by calculating the time it takes for $\mu$ to change by a factor of $1.2$ for 2048 sample ions 111We then rescale $\nu_{c}$ so that $\nu_{c}^{-1}$ measures factor-$e$ changes in $\mu$.. The scattering changes from being homogenous and weak at early times, to being stronger and localized around the SA-wave nodes at later times. This change is caused by the decay of oblique firehose modes into $k_{\parallel}\rho_{i}\sim 1,\,k_{\perp}\sim 0$ fluctuations [Fig. 1(c–d)], which scatter particles efficiently due to their small scale. These parallel modes are long lived (they are nonlinearly stabilized against cyclotron damping Gary and Saito (2003); Gary (2004)), as indicated by their presence after the large-scale SA wave has decayed and $\Delta p\sim 0$ [Fig. 1(d)].

Because $\omega_{A}\ll\nu_{c}\ll\Omega_{i}$, the plasma dynamics now resemble the Braginskii collisional limit Braginskii (1965) and the SA wave behaves as discussed in Squire et al. (2017). We illustrate the similarity in Fig. 2(b), which also shows $\delta B_{z}$ and $4\pi\Delta p/B^{2}$ for an SA wave governed by the Braginskii model (including heat fluxes; see Squire et al. (2017), App. B). The “humped” shape occurs because the perturbation splits into regions where $4\pi\Delta p\approx-B^{2}$ and $d\delta B_{z}/dt<0$ (around the antinodes), and regions where $4\pi\Delta p>-B^{2}$ and $\delta B_{z}=0$ (these spread from the nodes). The wave decay rate is determined by $\nu_{c}$ which is sufficiently large [$\nu_{c}/\omega_{A}\sim\beta(\delta B_{\perp}/B_{0})^{2}$] that the wave decays within $t/\tau_{A}<1$. We note parenthetically that the wave decay generates a $\delta B_{y}$ perturbation (see Fig. 1), although its origin is currently unclear.

As the large-scale SA wave decays, it heats the plasma. This process does not involve a turbulent cascade, but rather the direct transfer of large-scale mechanical energy into thermal energy. This heating is essentially viscous dissipation, with particle scattering from microscale fluctuations controlling the effective viscosity and making the process irreversible. In Fig. 3, we compare the measured $\partial_{t}E_{\mathrm{th}}$ with heating due to the SA wave decay. Although the agreement is not perfect due to spurious grid heating Birdsall and Langdon (1991) (tests at $\lambda_{A}/\rho_{i}=250$ show that this improves with ppc or reduced $\lambda_{A}/\rho_{i}$), the stages of wave decay are evident; e.g., the drop in $\partial_{t}E_{\mathrm{th}}$ as firehose fluctuations grow at $t/\tau_{A}\approx 0.075$, followed by heating as the large-scale $\delta B_{\perp}$ decays. Fig. 3 also shows that the overall energetics are well captured by considering only the large-scale dynamics, or by using the same effective collisionality as in Fig. 2(b). This supports closure models that approximate the effects of microinstabilities on large-scale dynamics without having to resolve the microscales.

*Shear-Alfvén traveling wave.—*Figures 4 and 5 illustrate the dynamics of the nonlinear SA traveling wave with $\lambda_{A}=250\rho_{i}$. The initial evolution differs from standing waves because $\langle dB/dt\rangle=0$ for an unperturbed traveling wave, implying that global (spatially constant) pressure anisotropy is created only as the wave decays Squire et al. (2016). The evolution broadly follows the expectations of Squire et al. (2016, 2017), proceeding in 4 stages: (i) the spatially dependent $dB/dt$ creates an anisotropy $\Delta p(x)\sim\beta^{1/2}\delta B_{z}^{2}\sin(2k_{\parallel}x)$; (ii) this $\Delta p$ damps the wave Squire et al. (2017); Lee and Völk (1973); Kulsrud (1978) causing $\langle B\rangle$ to decrease and thus $\langle\Delta p\rangle<0$; (iii) the wave consequently slows down, with $\delta u_{\perp}$ decaying faster than $\delta B_{\perp}$; and (iv) the wave excites oblique firehose modes, which subsequently scatter particles and cause the large-scale $\delta B_{\perp}$ to decay in a similar manner to the standing wave.

These stages are seen clearly in Figs. 4 and 5. In particular, note the global $\langle\Delta p\rangle<0$ that quickly develops [Fig. 5(b)] and the fast decay of $\delta u_{\perp}$ at early times [Fig. 5(c)]. By $t/\tau_{A}\approx 0.3$, once $\langle\Delta p\rangle$ has decreased sufficiently, the wave starts exciting oblique firehose modes. Unlike for standing waves, this occurs only in localized regions around the wavefronts (i.e., near where $\delta B_{z}=0$), because $|dB/dt|$ (and thus $|\Delta p|$) is largest in these regions [see, e.g., Fig. 4(c) at $x/\rho_{i}\approx 110$ and the shading in Fig. 5(b)]. Subsequently, Fig. 5(a) shows that the particle scattering is strongest behind the wavefronts. We interpret this as being due to the transition of oblique firehose modes into long-lived $k_{\parallel}\rho_{i}\sim 1$ fluctuations [see Fig. 4(d)], like in the standing wave. Because firehose modes are excited only briefly around the wavefronts, the scattering rate $\nu_{c}/\omega_{A}$ is lower than for the standing wave. Thus the final decay of the traveling wave’s $\delta B_{z}$ is slower than the standing wave’s, although it is qualitatively similar. At earlier times, the large-scale SA wave evolution matches well the predictions from a 1-D LF model at $\beta_{i}=100$ Squire et al. (2017) (dashed lines in Fig. 5).

*Discussion.—*We have presented hybrid-kinetic simulations of large-amplitude SA waves in a collisionless plasma. Our results demonstrate clearly the exceptional influence of microinstabilities on the large-scale ($\lambda_{A}\gg\rho_{i}$) dynamics of high-$\beta$ collisionless plasmas, illustrating how the evolution of self-excited oblique firehose modes controls the plasma’s fluid properties. The simulations also verify, using a realistic model with kinetic ions, that linearly polarized shear-Alfvénic perturbations do not exist in their linear wave form above the amplitude limit $\delta B_{\perp}/B_{0}\sim\beta^{\,-1/2}$ Squire et al. (2016). The SA wave dynamics depend strongly on how oblique firehose modes evolve as the plasma becomes stable ($\Delta p\gtrsim-B^{2}/4\pi$). We find that firehose fluctuations become parallel ($k_{\perp}=0$) and move to smaller scales ($k_{\parallel}\rho_{i}\sim 1$), surviving nonlinearly throughout the large-scale $\delta B_{\perp}$ decay and scattering particles at a high rate. These long-lived $k_{\parallel}\rho_{i}\sim 1$ modes cause SA standing-wave dynamics in a collisionless plasma to resemble those in a collisional (Braginskii) one Squire et al. (2017). The initial evolution of the traveling wave is effectively collisionless and matches analytic predictions Squire et al. (2017); however, after generating a global negative anisotropy and exciting firehose modes, its final decay resembles the standing wave. For both standing and traveling waves, the simulations provide an interesting example of direct transfer of energy from the largest scales to thermal energy and microscale fluctuations, without a turbulent cascade.

Our simulations cannot fully address what occurs at yet higher $\lambda_{A}/\rho_{i}$. This will depend on how oblique firehose modes decay and scatter particles, physics that is currently poorly understood. That said, it is clear that SA wave interruption provides a robust mechanism for dissipating energy directly from large-scale perturbations into heat and microinstabilities. Our results suggest that numerical models of weakly collisional high-$\beta$ plasmas would be better off damping large-amplitude SA waves, rather than letting them freely propagate. One concrete way to achieve this aim might be a LF model with pressure-anisotropy limiters Sharma et al. (2006) that enhance the collisionality to a rate that is determined by the large-scale Alfvén frequency. More work on developing and validating subgrid models of this kind is underway.

Given the strong deviations from MHD predictions, SA wave interruption could significantly impact the turbulent dynamics of weakly collisional plasmas in a variety of astrophysical environments Kunz et al. (2010). Some effects have already been observed in the $\beta\sim 1$ solar wind Bale et al. (2017). Other astrophysical plasmas—for instance the ICM, with $\beta\sim 100$ Enßlin and Vogt (2006); Rosin et al. (2011)—are likely to be more strongly affected by interruption, and work is underway to assess its impact on turbulence under such conditions.