Adaptive Critical Balance and Firehose Instability in an Expanding, Turbulent, Collisionless Plasma

Many space and astrophysical plasmas are magnetized and weakly collisional, with the Larmor radii of the constituent particles being many orders of magnitude below their Coulomb mean free paths (e.g., Schekochihin & Cowley 2006). This feature results in a complex interplay between a plasma's macrophysical evolution (e.g., due to expansion, compression, or large-scale shear) and its microphysical response (e.g., departures from local thermodynamic equilibrium, triggering of kinetic instabilities) (Schekochihin et al. 2005; Kunz et al. 2014a; Hellinger & Trávníček 2015; Riquelme et al. 2015; Sironi & Narayan 2015; Squire et al. 2017; Kunz et al. 2020). This interplay becomes increasingly complex when that macrophysical evolution induces or accompanies a cascade of turbulent fluctuations down to microphysical scales, a situation thought to be ubiquitous in the solar wind, low-luminosity black hole accretion flows, and the intracluster medium (e.g., Alexandrova et al. 2013; Yuan & Narayan 2014; Simionescu et al. 2019).

In this Letter, we investigate to what extent the basic building blocks of strong, incompressible, Alfvénic turbulence—namely, the existence of a conservative cascade from large (injection) to small (dissipative) scales, the locality of interactions between turbulent fluctuations, and a scale-by-scale balance between the characteristic linear oscillation time of the fluctuations and their nonlinear interaction time known as "critical balance" (Goldreich & Sridhar 1995; Mallet et al. 2015; Schekochihin 2020)—survive when subject to microphysical constraints dictated by the kinetic evolution of a collisionless plasma. Theoretical work describing magnetized turbulence in weakly collisional or collisionless plasma, but adopting a pressure-isotropic background, suggests that these organizing principles endure, with a local, conservative, Alfvénic cascade extending from macroscopic scales down to the ion-Larmor scale (Schekochihin et al. 2009). However, the assumption of an isotropic background pressure is not always justified; instead, the pressure tensor is more naturally anisotropic with respect to the magnetic field, with the evolution of field-parallel and perpendicular pressures influenced by approximate adiabatic invariance of the charged particles. How this pressure anisotropy alters the properties of Alfvénic turbulence has been a question of particular interest in recent years (e.g., Klein & Howes 2015; Kunz et al. 2015, 2018; Markovskii et al. 2019).

To address this question, we use results from a hybrid-kinetic simulation in which strong Alfvénic turbulence is driven in a collisionless, magnetized plasma undergoing steady expansion transverse to a mean magnetic field. This expansion drives pressure anisotropy in the plasma through approximate adiabatic invariance. We find that, despite the consequent decrease in the characteristic linear frequency and Alfvén ratio of the fluctuations, the Alfvénic cascade adapts to maintain critical balance. Eventually the plasma becomes unstable to kinetic firehose instabilities, which grow rapidly on ion-Larmor scales, scatter particles, and thereby impede the further production of pressure anisotropy. Even in this state, critical balance of the Alfvénic cascade persists, with the majority of the turbulent motions remaining stable.

2.1. Why Expansion?

Of the various types of macroscopic evolution that a turbulent, collisionless magnetized plasma can undergo, there are two compelling reasons to consider expansion.

First, plasma expansion on a timescale ${\tau }_{\exp }$ much larger than the inverse cyclotron frequency ${{\rm{\Omega }}}_{s}^{-1}$ of each particle species s (∈ {e, i} for an electron-ion plasma) provides a natural way to drive temperature anisotropy, Δ_s ≡ T_⊥s /T_∥s − 1 ≠ 0, where T_⊥s (T_∥s ) is the field-perpendicular (-parallel) component of the temperature of species s. For example, as plasma expands transversely to a mean magnetic ("guide") field, mass and magnetic-flux conservation dictate that the mean number density n_s of each species s and the guide-field strength B_g satisfy ${n}_{s},{B}_{{\rm{g}}}\propto {L}_{\,\perp }^{-2}$, where L_⊥ is the characteristic transverse size of the plasma (taken to be much larger than the thermal Larmor radius ρ_s of each species; the characteristic parallel size L_∥ is held fixed). Combined with conservation of the first and second adiabatic invariants, viz. T_⊥s ∝ B_g and ${T}_{\parallel s}\propto {\left({n}_{s}/{B}_{{\rm{g}}}\right)}^{2}$ (Chew et al. 1956; hereafter, CGL), these scalings imply a decreasing T_⊥s while T_∥s remains approximately constant. Thus, if Δ_s = 0 initially, then it becomes increasingly negative. Simultaneously, the parallel plasma beta parameters, ${\beta }_{\parallel s}\equiv 8\pi {n}_{s}{T}_{\parallel s}/{B}_{{\rm{g}}}^{2}$, increase. That the combination β_∥s Δ_s grows increasingly negative has two important consequences. First, the effective Alfvén speed

$$\begin{eqnarray}&&{v}_{{\rm{A}},\mathrm{eff}}\equiv {v}_{{\rm{A}}}{\left(1+\sum _{s}\displaystyle \frac{{\beta }_{\parallel s}{{\rm{\Delta }}}_{s}}{2}\right)}^{1/2}\end{eqnarray} \tag{ 1 }$$

drops below the conventional Alfvén speed v_A, tending toward zero as ∑_s β_∥s Δ_s → −2 (at which point there is no energetic cost to bending the field). Thus, the effective tension in the magnetic-field lines is reduced, with Alfvén waves becoming unstable for ∑_s β_∥s Δ_s < −2 (the "fluid firehose" threshold; Chandrasekhar et al. 1958; Parker 1958). Concurrently, when β_∥s Δ_s ≲ −1 the plasma becomes unstable to the kinetic parallel and oblique firehoses (Yoon et al. 1993). For plasma with β_∥i ≈ 2–4 and Maxwellian electrons, the oblique firehose operates when ${{\rm{\Delta }}}_{i}\lesssim -1.4{\beta }_{\parallel i}^{-1}$ (Hellinger & Matsumoto 2000), while the growth rate γ_f of the (threshold-less) parallel firehose satisfies γ_f ≳ 10⁻³Ω_i for ${{\rm{\Delta }}}_{i}\lesssim -1.1{\beta }_{\parallel i}^{-1}$ (Matteini et al. 2006). Both effects prompt several questions, including whether critical balance persists during the expansion, how the kinetic instabilities interact with the Alfvénic turbulence, and whether the turbulent motions themselves become unstable and disrupt the cascade.

The second reason to consider the problem of expanding Alfvénic turbulence is its relevance to the solar wind. A parcel of solar-wind plasma initially located at a large distance R ≫ L_⊥, L_∥ from the Sun and moving radially outwards at speed v_sw will undergo (approximately linear) expansion on a characteristic timescale ${\tau }_{\exp }=R/{v}_{\mathrm{sw}}$ (e.g., Matteini et al. 2012). Expansion is thought to play an important role in various key physical processes in the solar wind, including plasma heating, the generation of turbulence, and kinetic physics such as the production of temperature anisotropy (Velli et al. 1989; Verdini & Velli 2007; Chandran & Hollweg 2009; Matteini et al. 2013; Chandran & Perez 2019). There have therefore been many complementary investigations of expanding plasmas in the solar-wind context (e.g., Grappin et al. 1993; Liewer et al. 2001; Matteini et al. 2006; Camporeale & Burgess 2010; Hellinger et al. 2015; Hellinger 2017; Hellinger et al. 2019; Squire et al. 2020).

2.2. Hybrid-kinetic Description of Expanding Alfvénic Turbulence

We adopt a hybrid-kinetic approach to solve for the multiscale dynamics of Alfvénic turbulence in a collisionless, expanding plasma. A nonrelativistic, quasi-neutral (n ≡ n_i = n_e ) plasma with kinetic ions (mass m_i , charge e) and massless, fluid electrons is threaded by a uniform magnetic field ${{\boldsymbol{B}}}_{{\rm{g}}}={B}_{{\rm{g}}}\hat{{\boldsymbol{z}}}$ and subjected to a random, time-correlated, solenoidal driving force F (t, r ) ⊥ B _g. This driving is the same as described in Arzamasskiy et al. (2019); it is designed to mimic the action of random inertial forces arising from an anisotropic cascade of turbulent fluctuations at scales larger than the simulation domain. The electrons are assumed to be pressure-isotropic and isothermal with temperature T_e = T_i0, the initial ion temperature. A fourth-order hyper-resistivity is used to remove magnetic energy at the smallest scales.

The subsequent evolution of this plasma is solved using the second-order–accurate, particle-in-cell code Pegasus++ (Arzamasskiy et al. 2021, in preparation), which is an optimized implementation of the algorithms detailed in Kunz et al. (2014b). Well-resolved 3D hybrid-kinetic simulations of Alfvénic turbulence are essential for modeling this problem, in particular for simultaneously capturing both the turbulent cascade above and below ion-Larmor scales and the physics of ion-firehose instabilities. That being said, our treatment of the electrons as an isothermal, isotropic fluid precludes any kinetic instabilities driven by electron temperature anisotropy (e.g., the electron firehose; Li & Habbal 2000). While the properties of inertial-range Alfvénic fluctuations and ion-scale firehose instabilities are not expected to be affected appreciably by electron kinetics, it remains an open question as to how electron anisotropy affects the sub-ion-Larmor cascade of kinetic Alfvén waves (KAWs; see Sections 3.6.2, 4.4, and 4.5 of Kunz et al. 2018). For now, we simply note that, in the near-Earth solar wind, the electrons' collisional age seems to control the electron temperature anisotropy (Salem et al. 2003) and the total temperature anisotropy at β ≳ 1 is dominated by protons (Chen et al. 2016). By modeling only a single ion species (protons), our simulations also preclude some other effects thought to be relevant in the solar wind, e.g., instabilities driven by drifting helium ions (Verscharen et al. 2019).

To model the expansion, Pegasus++ enacts a coordinate transform from a comoving, nonexpanding frame (position vector r ) to the comoving expanding frame (position vector ${\boldsymbol{r}}^{\prime}$) using the time-dependent (diagonal) Jacobian transformation matrix ${\boldsymbol{\Lambda }}(t)\equiv \partial {\boldsymbol{r}}/\partial {\boldsymbol{r}}^{\prime}$, as in the Hybrid Expanding Box (HEB) model of Hellinger & Trávníček (2005, Appendix A). Pegasus++ solves the following modified versions of Faraday's and Ohm's laws in the expanding frame for the magnetic field ${\boldsymbol{B}}^{\prime} \equiv \lambda {{\boldsymbol{\Lambda }}}^{-1}{\boldsymbol{B}}$ and the electric field ${\boldsymbol{E}}^{\prime} \equiv {\boldsymbol{\Lambda }}{\boldsymbol{E}}$:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{{\boldsymbol{B}}}^{{\rm{{\prime} }}}}{{\rm{\partial }}{t}^{{\rm{{\prime} }}}}=-c{{\boldsymbol{\nabla }}}^{{\rm{{\prime} }}}\times {{\boldsymbol{E}}}^{{\rm{{\prime} }}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\boldsymbol{E}}^{\prime} =-\displaystyle \frac{{\boldsymbol{u}}^{\prime} }{c}\times {\boldsymbol{B}}^{\prime} -\displaystyle \frac{{T}_{e}}{{en}^{\prime} }{\boldsymbol{\nabla }}^{\prime} n^{\prime} +({\boldsymbol{\nabla }}^{\prime} \times {\boldsymbol{B}}^{\prime} )\times \displaystyle \frac{{{\boldsymbol{\Lambda }}}^{2}{\boldsymbol{B}}^{\prime} }{4\pi {en}^{\prime} \lambda },\end{eqnarray} \tag{ 3 }$$

where the primed-frame number density $n^{\prime} \equiv \lambda n$ and ion-flow velocity ${\boldsymbol{u}}^{\prime} \equiv {{\boldsymbol{\Lambda }}}^{-1}{\boldsymbol{u}}$, $\lambda \equiv \det \,{\boldsymbol{\Lambda }}$, and $t^{\prime} =t$. These fields are used to update the simulation ion-particle positions ${{\boldsymbol{r}}}_{p}^{{\prime} }={{\boldsymbol{\Lambda }}}^{-1}{{\boldsymbol{r}}}_{p}$ and velocities ${{\boldsymbol{v}}}_{p}^{{\prime} }={{\boldsymbol{\Lambda }}}^{-1}{{\boldsymbol{v}}}_{p}$ via

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{{\boldsymbol{r}}}_{p}^{{\prime} }}{{\rm{d}}t^{\prime} }={{\boldsymbol{v}}}_{p}^{{\prime} },\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}{{\boldsymbol{v}}}_{p}^{{\prime} }}{{\rm{d}}t^{\prime} } & = & \displaystyle \frac{e}{{m}_{i}}\,{{\boldsymbol{\Lambda }}}^{-2}\left[{\boldsymbol{E}}^{\prime} (t^{\prime} ,{{\boldsymbol{r}}}_{p}^{{\prime} })+\displaystyle \frac{{{\boldsymbol{v}}}_{p}^{{\prime} }}{c}\,\times \,{\boldsymbol{B}}^{\prime} (t^{\prime} ,{{\boldsymbol{r}}}_{p}^{{\prime} })\right]\\ & & +{{\boldsymbol{\Lambda }}}^{-1}\displaystyle \frac{{\boldsymbol{F}}(t^{\prime} ,{{\boldsymbol{r}}}_{p}^{{\prime} })}{{m}_{i}}-2{{\boldsymbol{\Lambda }}}^{-1}\displaystyle \frac{{\rm{d}}{\boldsymbol{\Lambda }}}{{\rm{d}}t^{\prime} }\,{{\boldsymbol{v}}}_{p}^{{\prime} }.\end{array}\end{eqnarray} \tag{ 5 }$$

The final (velocity-dependent) term in Equation (5) is straightforwardly incorporated into the semi-implicit Boris algorithm for solving particle trajectories alongside the ${{\boldsymbol{v}}}_{p}^{{\prime} }\,\times \,{\boldsymbol{B}}^{\prime}$ rotation. Quantities in the nonexpanding frame are easily obtained ex post facto.

The expansion is taken to be perpendicular to $\hat{{\boldsymbol{z}}}$ and linear in time: ${\boldsymbol{\Lambda }}(t)=\hat{{\boldsymbol{z}}}\hat{{\boldsymbol{z}}}+\left(1+t/{\tau }_{\exp }\right)\left({\mathsf{I}}-\hat{{\boldsymbol{z}}}\hat{{\boldsymbol{z}}}\right)$, where ${\tau }_{\exp }$ is the expansion time and I the unit dyadic. Thus the perpendicular size of the simulated plasma increases in time as ${L}_{\,\perp }(t)={L}_{\,\perp 0}(1+t/{\tau }_{\exp })$, while the parallel size of the simulated plasma remains constant, L_∥(t) = L_∥0. (We denote any given quantity X evaluated at the start of the simulation by X₀.) Magnetic-flux conservation then gives ${B}_{{\rm{g}}}(t)={B}_{{\rm{g}}0}{\left(1+t/{\tau }_{\exp }\right)}^{-2}$. This prescription is physically relevant to the expanding solar wind at ≳ 0.1 au, on account of the solar wind's constant speed and radial direction at those distances (Verscharen et al. 2019), although our treatment of the mean magnetic field as radial is a simplifying assumption.

2.3. Physical Setup

At the start of the simulation (time t₀), N_ppc = 10³ simulation ion-particles per cell are drawn randomly from a stationary Maxwellian distribution with temperature T_i0 and number density n₀ and placed uniformly in an elongated 3D computational domain of size ${L}_{x}\times {L}_{y}\times {L}_{z}={\left(65{\rho }_{i0}\right)}^{2}\times 390{\rho }_{i0}$ containing 256² × 1536 cells. At this box size and resolution, the captured wavenumbers are initially in the ranges k_(x,y) ρ_i0 ∈[0.097, 12.37] and k_z ρ_i0 ∈ [0.016, 12.37]. The initial ion beta parameter is β_∥i0 = 2, representative of near-Earth conditions in the solar wind (Matteini et al. 2007). Prior to initiating expansion, steady-state Alfvénic turbulence is generated in the plasma by forcing the particles with an F (t, r ) having the correlation time τ_A0/2π, where ${\tau }_{{\rm{A}}0}\equiv {L}_{z}/{v}_{{\rm{A}}0}\approx 552{{\rm{\Omega }}}_{i0}^{-1}$ is the initial Alfvén-crossing time, ${v}_{{\rm{A}}0}\equiv {B}_{{\rm{g}}0}/{\left(4\pi {m}_{i}{n}_{0}\right)}^{1/2}$ is the initial Alfvén speed, and Ω_i0 ≡ eB_g0/m_i c is the initial ion-cyclotron frequency. The magnitude of the force is such that critical balance is maintained for the box-scale fluctuations: u_rms/v_A0 ≈L_⊥/L_∥, where u_rms is the rms turbulent velocity. Assuming a −5/3 power-law scaling for turbulent fluctuations on scales larger than the box, the inferred perpendicular wavenumber at which the energy of the turbulent fluctuations becomes comparable to that of the guide magnetic field is ${k}_{\perp }^{\mathrm{outer}}\sim {10}^{-3}{\rho }_{i0}^{-1}$, a comparable degree of separation to that observed in the fast, β ≳ 1 solar wind (Wicks et al. 2010). This initial nonexpanding phase of the simulation lasts for five Alfvén-crossing times until t = 0, so that ${t}_{0}=-5{\tau }_{{\rm{A}}0}\approx -2758{{\rm{\Omega }}}_{i0}^{-1}$. The turbulent magnetic fields at t = 0 are visualized in Figure 1(a).

Zoom In Zoom Out Reset image size

The plasma's expansion is then initiated as described in Section 2.2, with ${\tau }_{\exp }=10{\tau }_{{\rm{A}}0}\approx 5515{{\rm{\Omega }}}_{i0}^{-1}$. This expansion time is comparable to the inferred Alfvén-crossing time at the outer scale $\sim 1/{k}_{\perp }^{{\rm{outer}}}$of the turbulence. It also means that the turbulent heating time ${\tau }_{\mathrm{heat}}\sim (3/2){T}_{i}{L}_{\,\perp }/({m}_{i}{u}_{\mathrm{rms}}^{3})\gtrsim 5{\tau }_{\exp }$ in our simulation; as a result, the thermodynamic evolution of the plasma is dominated not by turbulent heating but rather by the approximately double-adiabatic expansion and the feedback from firehose instabilities.

As the plasma expands, strong Alfvénic turbulence is driven continuously such that ${u}_{\mathrm{rms}}(t)\approx {L}_{\,\perp }(t){v}_{{\rm{A}}}(t)/{L}_{\parallel }=\mathrm{const}$. In retrospect, our results suggest that a more realistic forcing prescription would maintain critical balance adaptively at the outer scale using v_A,eff instead of v_A. However, this prescription requires a priori knowledge of the temperature anisotropy's evolution to evolve v_A,eff(t) and, moreover, becomes problematic if v_A,eff were to approach 0. In practice, we find the only consequence of using v_A(t) to determine the forcing amplitude to be a slight excess of energy in the turbulent fluctuations at the outer scale.

The overall plasma evolution is summarized in Figure 2. Panel (a) shows that the box-averaged magnetic-field strength B(t) and density n(t) decrease in tandem once the expansion begins, with $B(t)/{B}_{{\rm{g}}0}\approx n(t)/{n}_{0}={\left(1+t/{\tau }_{\exp }\right)}^{-2}$. The rms field strength actually decreases slightly slower due to the growth of the turbulent Alfvénic fluctuations, δ B_rms, relative to B_g(t) as the plasma expands (Figure 2(b)). This growth, also evident in Figure 1, is caused by wave-action conservation and by the build-up of residual magnetic energy in the fluctuations from the reduced energetic cost of bending field lines in a plasma with Δ_i < 0. Namely, the Alfvén ratio ${r}_{{\rm{A}}}\equiv 4\pi {m}_{i}{{nu}}_{\mathrm{rms}}^{2}/\delta {B}_{\mathrm{rms}}^{2}$ becomes smaller than unity as the expansion proceeds, an effect that may be compensated by instead using the "kinetic normalization" ${r}_{{\rm{A}},\mathrm{eff}}\equiv {r}_{{\rm{A}}}{\left(1+{\beta }_{\parallel i}{{\rm{\Delta }}}_{i}/2\right)}^{-1}$ (Chen et al. 2013). The associated relation $\delta {B}_{\mathrm{rms}}/{B}_{{\rm{g}}}\approx {\left(1+{\beta }_{\parallel i}{{\rm{\Delta }}}_{i}/2\right)}^{-1/2}\,{u}_{\mathrm{rms}}/{v}_{{\rm{A}}}$, when combined with critical balance of the box-scale fluctuations, viz. ${u}_{\mathrm{rms}}\sim [{L}_{\,\perp }(t)/{L}_{\,\parallel \,}]{v}_{{\rm{A}},\mathrm{eff}}(t)\propto (1+t/{\tau }_{\exp }){\left[1+{\beta }_{\parallel i}(t){{\rm{\Delta }}}_{i}(t)/2\right]}^{1/2}$ (Figure 2(b), dashed red line), implies $\delta {B}_{\mathrm{rms}}/{B}_{{\rm{g}}}\propto (1+t/{\tau }_{\exp })$ (Figure 2(b), dashed blue line), a manifestly good fit to the data.

Zoom In Zoom Out Reset image size

Another key property of Alfvénic turbulence in an expanding collisionless plasma is the decreasing characteristic frequency of the fluctuations. This feature is demonstrated by Figure 2(c), in which the red pluses track the time evolution of the energetically dominant ("peak") oscillation frequency of the fluctuations, ω_peak. ⁵ While some decrease in ω_peak is caused by the decreasing Alfvén speed, ${v}_{{\rm{A}}}={v}_{{\rm{A}}0}{\left(1+t/{\tau }_{\exp }\right)}^{-1}$ (dashed blue line), it is mostly due to the reduction in the effective Alfvén speed caused by β_∥i Δ_i becoming increasingly negative. Indeed, the effective Alfvén frequency of the box-scale fluctuations, 2π v_A,eff/L_∥ (solid blue line), matches the data well.

The production of negative temperature anisotropy during the expansion is shown in Figure 2(d). During the initial phase, the parallel (blue line) and perpendicular (red line) ion temperatures evolve approximately double-adiabatically: T_⊥i (t) ≈ T_⊥i (0)[B(t)/B(0)] (dashed red line) and ${T}_{\parallel i}(t)\approx {T}_{\parallel i}(0){\left[n(t)/n(0)\right]}^{2}{\left[B(t)/B(0)\right]}^{-2}$ (dashed blue line). However, at $t\approx {t}_{{\rm{f}}}\equiv 0.4{\tau }_{\exp }$, an abrupt change in the evolution of T_⊥i (t) and T_∥i (t) occurs, and the double-adiabatic predictions no longer hold. This change is coincident with Δ_i decreasing sufficiently (and β_∥i increasing sufficiently—see Figure 2(e)) that Δ_i ≲ −1.4/β_∥i (see Figure 2(f)), at which point the plasma is unstable to kinetic firehose instabilities. Such firehose fluctuations, visually evident near the ion-Larmor scale in Figure 1(b), are characterized later in this section.

Figures 3(a) and (b) display 1D power spectra of the velocity (E_u ) and magnetic (E_B ) fluctuations at select times as functions of the perpendicular wavenumber k_⊥ normalized to the time-dependent ion-Larmor scale, ${\rho }_{i}\equiv {\left[2{T}_{\perp i}(t)/{m}_{i}\right]}^{1/2}/{{\rm{\Omega }}}_{i}(t)$. Their overall shapes are similar to those found in prior hybrid-kinetic simulations of nonexpanding, β_∥i ∼ 1 turbulence (e.g., Arzamasskiy et al. 2019): ${E}_{u}({k}_{\perp }),{E}_{B}({k}_{\perp })\propto {k}_{\perp }^{-5/3}$ in the inertial ("MHD") range, before steepening at k_⊥ ρ_i ≳ 1 due to finite-Larmor-radius effects. The "break point" at which this steepening occurs, ${\left({k}_{\perp }{\rho }_{i}\right)}_{\mathrm{break}}$ (blue curve, Figure 3(b) inset), decreases at a rate quantitatively consistent with theoretical expectations (Kunz et al. 2018, Section 3.6.4) that

$$\begin{eqnarray}&&{\left({k}_{\perp }{\rho }_{i}\right)}_{\mathrm{break}}\propto {\left(1+{\beta }_{\parallel i}{{\rm{\Delta }}}_{i}/2\right)}^{1/4}\,{\beta }_{\parallel i}^{-1/4}{\left({T}_{\perp i}/{T}_{\parallel i}\right)}^{1/2}\end{eqnarray} \tag{ 6 }$$

(red curve, Figure 3(b) inset). ⁶ These spectral features are maintained throughout the expansion, even for t ≳ t_f.

Zoom In Zoom Out Reset image size

Having provided evidence that various properties of the large-scale fluctuations adapt to the changing background pressure anisotropy in a manner consistent with critical balance, we now utilize the spectra in Figure 3 to show that critical balance is in fact maintained adaptively, scale by scale, as the plasma expands. We do so by computing the spectral anisotropy of the fluctuations using an approach proposed by Cho & Lazarian (2009) in which the characteristic parallel wavenumber k_∥(k_⊥) of magnetic-field fluctuations with perpendicular wavenumber k_⊥ is determined from their rms parallel lengthscale (see their Equation (34)). For fluctuations with a given k_⊥, this measure is most sensitive to the energetically dominant fluctuations with the largest k_∥, and so the approach can be used to determine the linear frequency ω_A ≡ k_∥ v_A,eff of these fluctuations and compare it with their nonlinear frequency ${\omega }_{\mathrm{nl}}\equiv {k}_{\perp }{\left[{k}_{\perp }{E}_{u}({k}_{\perp })+{v}_{{\rm{A}},\mathrm{eff}}^{2}{k}_{\perp }{E}_{B}({k}_{\perp })/{B}_{{\rm{g}}}^{2}\right]}^{1/2}$. In critically balanced turbulence, the turbulent energy is concentrated in a cone satisfying ω_A ≲ ω_nl, with the edge of the cone having ${k}_{\parallel }\propto {k}_{\perp }^{2/3}$ (Goldreich & Sridhar 1995).

The result of this calculation is shown at different times in Figure 3(c). At t = 0, the measured spectral anisotropy in the inertial range is consistent with the critical-balance scaling ${k}_{\parallel }\propto {k}_{\perp }^{2/3}$. As the expansion proceeds, this scaling is maintained as the overall degree of anisotropy decreases in tandem with the decreasing aspect ratio of the plasma. Furthermore, the inset shows that ω_A ≈ ω_nl scale by scale; thus critical balance holds adaptively. At t ≈ t_f, firehose modes (which, unlike the Alfvénic fluctuations, are not highly elongated in the field-parallel direction) emerge and bias slightly the calculated scaling of k_∥(k_⊥) in the inertial range. To mitigate this bias, a weight function is applied to the magnetic field that preferentially removes firehose-unstable regions before evaluating k_∥. Using this weight function, adaptive critical balance of the Alfvénic cascade is seen to persist. ⁷

In summary, no dramatic alterations to the fundamental nature of the Alfvénic turbulence are observed during expansion, even when kinetic-scale firehose modes are present. Importantly, there is no noticeable destabilization of the inertial-range Alfvénic cascade. This result is due to the efficient regulation of the box-averaged temperature anisotropy, which (as shown in Figure 2(f)) barely drops below Δ_i ≈ −1.4/β_∥i . While this value of Δ_i is negative enough to destabilize the plasma to kinetic firehose instabilities, it is above the "fluid" firehose instability threshold Δ_i = − 2/β_∥i below which ${v}_{{\rm{A}},\mathrm{eff}}^{2}\leqslant 0$ and Alfvén waves cease to propagate.

The character of the kinetic-scale firehose fluctuations can be ascertained by examining the 2D Fourier spectrum of the magnetic field E_B in $[{k}_{r}\equiv {\left({k}_{x}^{2}+{k}_{y}^{2}\right)}^{1/2},{k}_{z}]$ space. At $t\lesssim 0.32{\tau }_{\exp }$ (Figure 4(a), top), spectral power is concentrated in the region of (k_r , k_z ) space that satisfies k_z ≪ k_r , affirming the quasi-perpendicular nature of the Alfvénic cascade. By $t=0.4{\tau }_{\exp }$ (Figure 4(a), bottom), an additional region with spectral power is clearly visible, with its centroid located at (k_r ρ_i , k_z ρ_i ) ≈ (0.4, 0.3). We associate this power with growing oblique firehose fluctuations. ⁸ These fluctuations can be visualized by isolating the "firehose" part δ B_x,f of the magnetic field using a Fourier-space mask that filters out quasi-perpendicular modes; the region of (k_r , k_z ) space identified as the firehose part is indicated by the shaded region in Figure 4(a). While the Alfvénic turbulence does not evolve qualitatively during the time interval $t/{\tau }_{\exp }\in [0.32,0.40]$, the firehose fluctuations increase their amplitudes significantly (see Figure 4(b)). Figure 4(c), which shows the evolution of the magnetic energy of ion-Larmor-scale modes at different angles to the guide field, confirms that oblique firehose modes are unstable, with maximum growth rate comparable to that predicted by linear theory at β_⊥i ≈ 3.6 and Δ_i ≈ −0.4, viz. ${\gamma }_{{\rm{f}}\perp }\approx 0.02{{\rm{\Omega }}}_{i}\approx 120{\tau }_{\exp }^{-1}$ at (k_f⊥ ρ_i , k_f∥ ρ_i ) ≈ (0.4, 0.3). Parallel firehose modes (measured in region III of Figure 4(a)) are also unstable, but they have a significantly smaller amplitude than the oblique modes.

Zoom In Zoom Out Reset image size

The firehose fluctuations efficiently regulate the temperature anisotropy, even though their saturated rms magnetic-field strength is much smaller than that of the Alfvénic fluctuations at equivalent wavenumbers. They do so by pitch-angle scattering the ions so that the particles' first adiabatic invariants ($\mu \equiv {m}_{i}{v}_{\perp }^{2}/2B$, where v_⊥ is the peculiar perpendicular velocity) are no longer conserved (see Figure 4(d), orange line). The effective collisionality of this anomalous scattering, ν_c, may be estimated using the relation $\dot{\overline{\mu }}=-{\nu }_{{\rm{c}}}(\overline{{\rm{\Delta }}{T}_{i}/B})$, where the overline denotes a box average, ΔT_i ≡ T_⊥i − T_∥i , and $\dot{\overline{\mu }}$ is the rate of change of $\overline{\mu }\equiv \overline{{T}_{\perp i}/B}$. Figure 4(d) indicates that ${\tau }_{\exp }{\nu }_{{\rm{c}}}\ll 1$ for t < t_f (i.e., μ is approximately conserved pre-firehose), while ${\tau }_{\exp }{\nu }_{{\rm{c}}}\sim 1$ for t ≳ t_f (i.e., μ is significantly broken by the firehose fluctuations).

A simple model for ν_c may be constructed by adopting three assumptions: (i) that n(t)/n₀ ≈ B(t)/B_go; (ii) that contributions from heat fluxes and turbulent heating to the temperature anisotropy are negligible (the latter being because ${\tau }_{\mathrm{heat}}\gg {\tau }_{\exp };$ see Section 2.3); and (iii) that ${\beta }_{\parallel i}{{\rm{\Delta }}}_{i}\approx \mathrm{const}$ after t = t_f. Under these conditions, the CGL equations (including collisions) become ${\rm{d}}\mathrm{ln}({T}_{\perp i}/B)/{\rm{d}}t=-{\nu }_{{\rm{c}}}({T}_{\parallel i}/{T}_{\perp i}){{\rm{\Delta }}}_{i}$ and ${\rm{d}}\mathrm{ln}({T}_{\parallel i}{B}^{2}/{n}^{2})/{\rm{d}}t\approx {\rm{d}}\mathrm{ln}{T}_{\parallel i}/{\rm{d}}t\,=\,2{\nu }_{{\rm{c}}}{{\rm{\Delta }}}_{i}$. The third assumption then implies ${\nu }_{{\rm{c}}}\approx {\left(3{{\rm{\Delta }}}_{i}\right)}^{-1}\,{\rm{d}}\mathrm{ln}B/{\rm{d}}t\equiv {\nu }_{{\rm{c}}}^{{CGL}}$. The agreement between this model (Figure 4(d), blue line) and ν_c evaluated directly from the simulation is good, although ν_c fluctuates significantly. A direct calculation of the mean μ-breaking time of ∼ 10⁴ tracked particles, following Kunz et al. (2014a, 2020) and Squire et al. (2017), yields an effective collisionality $\simeq {\nu }_{{\rm{c}}}^{{CGL}}$ for t ≳ t_f. Setting ${\nu }_{{\rm{c}}}={\nu }_{{\rm{c}}}^{{CGL}}$ in the above equations leads to a simple equation for the parallel temperature, ${\rm{d}}\mathrm{ln}({T}_{\parallel i}/{B}^{2/3})/{\rm{d}}t\,=\,0$, so that ${T}_{\parallel i}(t)\approx {T}_{\parallel i}({t}_{{\rm{f}}}){\left[B(t)/B({t}_{{\rm{f}}})\right]}^{2/3}$. Further setting Δ_i ≈ −1.4/β_∥i yields ${T}_{\perp i}(t)\approx {T}_{\parallel i}(t)-1.4{B}_{{\rm{g}}}^{2}(t)/8\pi n(t)$. This model is plotted in Figure 2(d); given its simplicity, its agreement with the actual result is remarkable.

The regulation of temperature anisotropy can be elucidated further by considering PDFs of the simulation data in the (β_∥i , Δ_i ) phase space (e.g., Bale et al. 2009). Figure 5(a) shows these PDFs at different stages: at the expansion's start (t = 0), at t = t_f, and more than a full expansion time after t = 0; the phase-space trajectory of the PDF's average is indicated in the final panel by the black solid line. In all three cases, the relatively small dispersion in β_∥i and Δ_i is consistent with the small rms amplitude of the turbulent fluctuations. The temperature anisotropy clearly approaches the oblique firehose instability threshold Δ_i = −1.4/β_∥i (dashed line) and subsequently evolves along marginal instability.

Zoom In Zoom Out Reset image size

Despite the success of our collisionality model, the compartmentalization of all of the kinetic physics into an effective collision frequency hides some interesting emergent features in the ion distribution function f(v_∥, v_⊥). Figure 5(b) shows the difference between f and a Maxwellian distribution with the same temperature as f at three different times during the expansion (with all velocities normalized by the initial thermal speed ${v}_{\mathrm{th}i0}\equiv {\left(2{T}_{i0}/{m}_{i}\right)}^{1/2}$). For comparison, the difference between f and a bi-Maxwellian distribution with the same values of T_∥i and T_⊥i as f is also shown. Prior to the start of the expansion, the slight deficit of particles with (peculiar) parallel velocities v_∥ just below the Alfvén velocity v_A (Figure 5(b), left panel) is indicative of collisionless damping of the (kinetic) Alfvénic fluctuations. Once the expansion begins, these deviations are dwarfed by the expansion-driven temperature anisotropy (Figure 5(b), middle panel), which, on account of approximate double-adiabaticity, causes f to look like a bi-Maxwellian. However, by late times in the simulation, significant deviations from a bi-Maxwellian are evident (Figure 5(b), right panel), a finding seen in previous studies of the firehose instability (e.g., Hellinger 2017). In particular, the distribution function integrated over perpendicular velocities, $f({v}_{\parallel })\equiv {\int }_{0}^{\infty }{\rm{d}}{v}_{\perp }{v}_{\perp }f$, exhibits a flattened core (Figure 5(c)); the distribution function integrated over parallel velocities, $f({v}_{\perp })\equiv {\int }_{-\infty }^{\infty }{\rm{d}}{v}_{\parallel }\,f$, shows that the anisotropy of the distribution function at subthermal velocities is much more pronounced than in a bi-Maxwellian. These features can be attributed to resonant interactions between ions and the oblique firehose modes (Bott et al. 2021, in preparation).

That the nonlinear interactions between Alfvénic fluctuations adapt to satisfy critical balance, even as the characteristic linear frequency of those fluctuations is reduced by pressure anisotropy, is a vivid illustration of the complex interplay between velocity space and configuration space that is central to collisionless plasma physics. This interplay is made richer at β ≳ 1 by the emergence of ion-Larmor-scale firehose fluctuations, which establish a direct link between the microscales and macroscales by regulating the pressure anisotropy and thereby controlling the effective tension of magnetic-field lines. Despite the small-scale injection of magnetic energy by the firehose, those fluctuations are not sufficient in amplitude to contribute significantly to the magnetic power spectrum (at least perpendicular to the guide field). This finding should ease the concern expressed in Bale et al. (2009) that "these local [kinetic] instabilities ... may confuse the interpretation of solar wind magnetic power spectra." From the standpoint of the Alfvénic cascade, the most important (and potentially observable) roles played by the firehose are as a direct regulator of pressure anisotropy and an indirect mediator of adaptive critical balance and the transition to the KAW range.

The evolution of purely decaying, magnetized turbulence in an expanding, collisionless plasma with β_∥i ≳ 1 was recently investigated by Hellinger et al. (2019) using HEB simulations. In their setup, an isotropic spectrum of Alfvénically polarized waves (amplitude δ B_rms/B_g = 0.24) was initiated inside a cubic simulation domain with 512² × 256 cells spanning ${L}_{\,\perp }^{2}\times {L}_{\parallel }={\left(82{\rho }_{i0}\right)}^{3}$, before transverse expansion was introduced (${\tau }_{\exp }={10}^{4}{{\rm{\Omega }}}_{i0}^{-1}$) and the system evolved. The initial ion distribution function had nonzero temperature anisotropy, Δ_i0 = −0.25, with β_∥i0 = 2.4. Where there is overlap with their results, we find agreement: efficient regulation of the temperature anisotropy by kinetic firehose instabilities, persistence of a quasi-perpendicular Alfvénic cascade independent of firehose fluctuations, and distortion of the particle distribution function away from a bi-Maxwellian. There are, however, two important distinctions worth highlighting. First, because of the shape of the simulation domain (L_∥ ≤ L_⊥) in Hellinger et al. (2019), the Alfvénic fluctuations are likely not in critical balance. Alfvénic fluctuations in an MHD turbulent cascade become critically balanced for isotropic outer-scale fluctuations at a scale ${\lambda }_{\mathrm{CB}}\sim {L}_{\parallel }{\left(\delta {B}_{\mathrm{rms}}/{B}_{{\rm{g}}}\right)}^{3/2}$ (Schekochihin 2020); given the parameters in Hellinger et al. (2019), we estimate λ_CB ≈ 0.1L_∥ ∼ ρ_i , placing the entire inertial range in the weak-turbulence regime. Our demonstration of adaptive critical balance of strong Alfvénic turbulence when the distribution function is anisotropic (even unstably so) is one of our key results. Second, we followed the evolution of the turbulence for well over an expansion time, and so could confirm that the temperature anisotropy remains pinned to the kinetic firehose instability threshold as the expansion proceeds. This is an important result for solar-wind applications, because the expansion time there is comparable to the turnover time (and thus the characteristic decay time) of the outer-scale turbulent eddies.

Our conclusions may not hold for plasmas with much higher β_∥i than have been considered here. First, it is possible to show using linear theory that, if ${\tau }_{\exp }\lesssim 10{\beta }_{\parallel i}^{3/2}{\left(\mathrm{ln}{\beta }_{\parallel i}\right)}^{1/2}{{\rm{\Omega }}}_{i}^{-1}$, then Δ_i would not be regulated fast enough by the oblique firehose to remain > −2/β_∥i . In this case, ${v}_{{\rm{A}},\mathrm{eff}}^{2}$ would pass through 0 and the entire inertial-range Alfvénic cascade would be destabilized. For the value of ${\tau }_{\exp }$ used in our simulation, we expect this to occur for β_∥i ≳ 50. Second, negative pressure anisotropy driven by the Alfvénic fluctuations themselves can "interrupt" the fluctuations if $\delta {B}_{\mathrm{rms}}/{B}_{{\rm{g}}}\gtrsim {\beta }_{\parallel i}^{-1/2}$, by nullifying the restoring tension force and exciting a sea of scattering firehose fluctuations (Squire et al. 2017). An investigation of strong Alfvénic turbulence at such high beta is already underway.

A.F.A.B., M.W.K., and E.Q. were supported by DOE awards DE-SC0019046 and DE-SC0019047 made through the NSF/DOE Partnership in Basic Plasma Science and Engineering. Support for L.A. was provided by the Institute for Advanced Study. Support for J.S. was provided by Rutherford Discovery Fellowship RDF-U001804 and Marsden Fund grant UOO1727, which are managed through the Royal Society Te Apārangi. High-performance computing resources were provided by the Texas Advanced Computer Center at The University of Texas at Austin under grant number TG-AST160068 and the PICSciE-OIT TIGRESS High Performance Computing Center and Visualization Laboratory at Princeton University. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant OCI-1053575. This work benefited from useful conversations with Silvio Sergio Cerri and Alexander Schekochihin, and especially from contributions by Ryan Golant to the implementation of the expanding box in Pegasus++ while a Princeton University undergraduate in 2019.