2.1. Plasma parameters

The plasma parameters play a critical role in determining the nature of the turbulent fluctuations and their interaction with the underlying particle velocity distributions. For simplicity of notation, we assume here a fully ionized, hydrogenic plasma of protons and electrons with bi-Maxwellian equilibrium velocity distributions characterized by separate parallel and perpendicular temperaturesFootnote ⁴ for each species, $T_{\parallel,s}$ and $T_{\perp,s}$. In this idealized case, the equilibrium ion and electron number densities are equal, $n_i=n_e$.

The first and most important plasma parameter is the parallel ion plasma beta, $\beta _{\parallel,i} = 8 {\rm \pi}n_i T_{\parallel i}/B^2$, which characterizes the ratio of parallel thermal pressure to the magnetic pressure in the plasma and dominantly controls the phase speeds of different linear wave modes in a magnetized plasma. Note that with definitions for the parallel ion thermal velocity $v_{t\parallel,i}^2 = 2 T_{\parallel, i}/m_i$ and the Alfvén velocity $v_A^2 = B^2/(4 {\rm \pi}n_i m_i)$, the parallel ion plasma beta can be alternatively expressed as $\beta _{\parallel,i} = v_{t\parallel,i}^2 /v_A^2$. The definitions of all characteristic plasma quantities and useful conversion formulas are collected in table 2. The second parameter is the ion-to-electron parallel temperature ratio $\tau _\parallel =T_{\parallel,i}/T_{\parallel, e}$. The third and fourth parameters describe the ion and electron temperature anisotropies, $A_i=T_{\perp,i}/T_{\parallel,i}$ and $A_e=T_{\perp,e}/T_{\parallel,e}$. The fifth parameter characterizes the collisionality of the plasma by comparing the mean free path scale for electron–electron collisions relative with the parallel driving scale $k_{\parallel 0} \lambda _{{\rm mfp},e}$.

Table 2. Definitions of characteristic plasma quantities in cgs units and useful conversion relations, where $\beta _{\parallel,i}$, $\tau _\parallel$, $A_i$ and $\mu$ are defined in table 3. $^*$Note that the linear wave damping rate is defined unconventionally here such that $\gamma >0$ corresponds to wave damping, while $\gamma <0$ corresponds to wave growth; this choice eliminates negative signs in the numerous formulas for damping rates presented in this paper.

Although for a hydrogenic plasma of protons and electrons the ion-to-electron mass ratio $\mu =m_i/m_e$ has a fixed physical value of 1836, this ratio is included as an additional independent parameter in the calculations below since a reduced mass ratio is often used in numerical simulations of plasma turbulence for reasons of computational efficiency. Note also that, once $\mu$ and $\tau _\parallel$ are specifiedFootnote ⁵ , the single dimensionless parameter $k_{\parallel 0} \lambda _{{\rm mfp},e}$ is sufficient to specify all possible charged particle collision rates, since $\nu _{ee}/(k_{\parallel 0} v_{t\parallel,e}) = 1/(k_{\parallel 0} \lambda _{{\rm mfp},e})$ and $\nu _{ei} \sim \nu _{ee}$, $\nu _{ii} \sim \mu ^{-1/2} \tau _\parallel ^{-3/2} \nu _{ee}$ and $\nu _{ie} \sim \mu ^{-1}\tau _\parallel ^{-3/2}\nu _{ee}$.

Therefore, $(\beta _{\parallel,i},\tau _\parallel, A_i, A_e, k_{\parallel 0} \lambda _{{\rm mfp},e})$ are five key dimensionless plasma parameters that govern the properties of plasma turbulence, along with the fixed mass ratio $\mu$. These plasma parameters are summarized in table 3.

Table 3. For the general case of homogeneous turbulence in a fully ionized, hydrogenic plasma with bi-Maxwellian equilibrium velocity distributions, these are the fundamental dimensionless parameters of kinetic plasma turbulence, including plasma parameters and turbulence parameters. Note that the Boltzmann constant is absorbed to give temperature in units of energy.

Note that, if the plasma contains multiple ion species, possibly each with multiple charge states, the number of parameters increases dramatically, requiring a minimum of five new dimensionless parameters for each new species $s$: (i) mass ratio $m_s/m_i$, (ii) number density ratio $n_s/n_e$, (iii) charge ratio $q_s/q_e$, (iv) parallel temperature ratio $T_{\parallel, s}/T_{\parallel, i}$ and (v) temperature anisotropy $T_{\perp,s}/T_{\parallel,s}$. Partially ionized plasmas will similarly require additional dimensionless parameters to be characterized. Since the intention here is not to treat kinetic plasma turbulence in full generality, but rather to describe how one can use the Pi theorem and turbulence scaling theories to minimize the number of relevant dimensionless parameters, the remainder of this paper will focus on the case of a fully ionized, hydrogenic plasma of protons and electrons with bi-Maxwellian equilibrium velocity distributions. Refinements of turbulent heating models to include minor ions, such as the $\alpha$ particles that comprise a small fraction of ions in space and astrophysical plasmas, can be constructed following the same principles presented here.

2.2. Turbulence parameters

The turbulence parameters characterize the scale, amplitude and properties of the turbulent driving though five additional parameters. For simplicity, the focus here is restricted to only the case of turbulence in a spatially homogeneous system in which the scale length of gradients in the equilibrium quantities are much larger than the turbulent outer scale of the system. Dimensionless ratios of length scales are notationally simplified by characterizing the length scales of the turbulent fluctuations in terms of their perpendicular and parallel wavenumbers. Below quantities at the outer scale of the turbulence, commonly denoted the driving scale or energy injection scale, are denoted by the subscript ‘0’.

The driving of the turbulence can be characterized by three dimensionless parameters. The first parameter is the perpendicular wavenumber of the driving-scale fluctuations normalized by the ion Larmor radius, $k_{\perp 0} \rho _i$. The second parameter is the wavevector anisotropy of the fluctuations at the driving scale, $k_{\parallel 0} / k_{\perp 0}$. The third parameter is the nonlinearity parameter that characterizes the amplitude of the driving, $\chi _0 = (k_{\perp 0} \delta B_{\perp 0} \theta _{0})/(k_{\parallel 0} B_0)$, where $\delta B_{\perp 0}$ is the amplitude of the perpendicular magnetic field fluctuations, $B_0$ is the equilibrium magnetic field magnitude and $\theta _0$ is the angle of dynamic alignment at the driving scale in radians (Boldyrev Reference Boldyrev2006). In the next section, we demonstrate how modern scaling theories for anisotropic magnetohydrodynamic (MHD) turbulence can be used to combine these three dimensionless parameters into a single parameter, the isotropic driving wavenumber, $k_0 \rho _i$. Note also that, in principle, the dynamic alignment angle $\theta _0$ at driving scale is an additional dimensionless parameter (measured in radians), but dynamic alignment is typically understood to be a phenomenon that develops self-consistently due to the turbulence as it cascades to smaller scales, so we neglect this parameter at the outer scaleFootnote ⁶ , generally taking $\theta _0 \sim 1$ rad.

The nature of the turbulent fluctuations at the driving scale also plays an important role in determining the ultimate fate of the turbulent energy. The fourth turbulence parameter is the imbalance parameter, measuring the ratio of the amplitudes in upward-to-downward Elsasser fields, $Z_0^+/Z_0^-$, where the Elsasser fields at the outer scale are defined by $\boldsymbol {Z}_0^\pm = \delta \boldsymbol {U}_{\perp, 0} \pm \delta \boldsymbol {B}_{\perp, 0}/\sqrt {4 {\rm \pi}n_i m_i}$ (Elsasser Reference Elsasser1950), and where $\delta \boldsymbol {U}_{\perp, 0}$ and $\delta \boldsymbol {B}_{\perp, 0}$ are the perpendicular perturbed plasma flow velocity and perpendicular perturbed magnetic field of the turbulent fluctuations at the driving scale. Formally, $\delta \boldsymbol {U}_{\perp, 0}$ and $\delta \boldsymbol {B}_{\perp, 0}$ are vectors in the plane perpendicular to the equilibrium magnetic field $\boldsymbol {B}_0$, but to calculate the ratio $Z_0^+/Z_0^-$, we simply take the amplitudes of the Elsasser fields, $Z_0^\pm =|\boldsymbol {Z}_0^\pm |$. Physically, Alfvén wave packets propagating up the mean magnetic field have $Z^+=0$ and $Z^-\ne 0$, and those propagating down the field have $Z^+\ne 0$ and $Z^-= 0$. The Elsasser fields effectively characterize the wave energy flux in each direction along the magnetic fieldFootnote ⁷ of the typically dominant Alfvénic fluctuations in heliospheric plasma turbulence (Tu & Marsch Reference Tu and Marsch1995; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Bruno & Carbone Reference Bruno and Carbone2013). Since the nonlinear interactions that mediate the turbulent cascade in Alfvénic turbulence at MHD scales arise only between counterpropagating and perpendicularly polarized Alfvén waves (Kraichnan Reference Kraichnan1965; Sridhar & Goldreich Reference Sridhar and Goldreich1994; Goldreich & Sridhar Reference Goldreich and Sridhar1995; Ng & Bhattacharjee Reference Ng and Bhattacharjee1996; Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000; Howes & Nielson Reference Howes and Nielson2013; Howes Reference Howes2015a), the imbalance plays an important role in characterizing the turbulence. The importance of $Z_0^+/Z_0^-$ under sufficiently imbalanced conditions, such as occurs in the inner heliosphere where wave energy fluxes are dominantly anti-sunward, has been demonstrated clearly with the recent discovery of a new phenomenon known as the helicity barrier in plasma turbulence (Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021; Squire et al. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2022; Squire, Meyrand & Kunz Reference Squire, Meyrand and Kunz2023).

A fifth turbulence parameter is the ratio of the energy in compressible fluctuations to that in incompressible fluctuations, $E_{{\rm comp}}/E_{{\rm inc}}$. Since a magnetized plasma supports two distinct compressible wave modes, the fast and slow magnetosonic waves (differing by whether the thermal and magnetic pressure perturbations associated with the wave act in conjunction or in opposition), one could further characterize the compressible component of the turbulent fluctuations by the additional ratio of the energy in the fast magnetosonic fluctuations to the total energy of compressible fluctuations $E_F/E_{{\rm comp}}$, where the total compressible energy is the sum of the contributions from fast-mode and slow-mode fluctuations, $E_{{\rm comp}}=E_F+E_s$. The ratio $E_F/E_{{\rm comp}}$ is not included here as one of the fundamental parameters based on a study of the inertial range in solar wind turbulence at 1 AU that showed negligible energy in fast-mode mode fluctuations, or $E_F/E_{{\rm comp}} \ll 1$ (Howes et al. Reference Howes, Bale, Klein, Chen, Salem and TenBarge2012a). If significant energy exists in fast-mode fluctuations in other turbulent astrophysical plasmas, $E_F/E_{{\rm comp}}$ may need to be included.

Therefore, $(k_{\perp 0} \rho _i,k_{\parallel 0}/k_{\perp 0},\chi _0,Z_0^+/Z_0^-,E_{{\rm comp}}/E_{{\rm inc}})$ are five key dimensionless turbulence parameters that govern the properties of plasma turbulence. These turbulence parameters are summarized in table 3.

If the restriction to spatially homogeneous turbulent plasmas is relaxed, the large-scale context of the turbulent plasma can also impact the nonlinear dynamics of the turbulence and the evolution of the turbulent fluctuations. In the solar wind, the spherical expansion of the solar wind governs the evolution of the plasma equilibrium temperatures with heliocentric radius, possibly triggering kinetic temperature anisotropy instabilities that can mediate non-local transfer of energy to small scales, as depicted in figure 1(b). In rotating plasmas, such as found in the solar convection zone, the turbulence may be driven with helical motions that can initiate the growth of magnetic field energy through a dynamo effect (Parker Reference Parker1955; Glatzmaier Reference Glatzmaier1984, Reference Glatzmaier1985; Parker Reference Parker1993; Ossendrijver Reference Ossendrijver2003). Differential rotation in astrophysical accretion disks, which can trigger plasma turbulence through the magnetorotational instability (Balbus & Hawley Reference Balbus and Hawley1991; Hawley & Balbus Reference Hawley and Balbus1991; Balbus & Hawley Reference Balbus and Hawley1998), has been recently shown to generate turbulence with a dominant fraction of its energy in compressible fluctuations (Kawazura et al. Reference Kawazura, Schekochihin, Barnes, Dorland and Balbus2022). The inclusion of self-gravitational effects in the plasma, which play an important role in governing the dynamics of the multi-phase interstellar medium in molecular clouds and star-forming regions, adds significant additional complexity to the nature of the turbulent cascade. Finally, kinetic instabilities arising in the foot and ramp regions of collisionless shocks (Brown et al. Reference Brown, Juno, Howes, Haggerty and Constantinou2023) can generate turbulent fluctuations that are swept into the downstream region. In all of these astrophysically relevant cases, additional dimensionless turbulence parameters will be required to characterize properly the turbulent fluctuations and their evolution.

2.3. What happened to the Reynolds number?

Readers familiar with studies of MHD turbulence may find what appear to be two glaring omissions in the proposed list of plasma and turbulence parameters in table 3 since it does not include the Reynolds number, $Re = L U_0/\nu$, or the magnetic Reynolds number, $Re_M = L U_0/\eta$. The Reynolds number is a dimensionless quantity that characterizes the ratio of the amplitude of the convective term $\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {U}$ to that of the viscous diffusion term $\nu \nabla ^2 \boldsymbol {U}$ in the MHD momentum equation; the magnetic Reynolds number characterizes the ratio of the amplitude of the inductive term $\boldsymbol {\nabla } \times (\boldsymbol {U} \times \boldsymbol {B})$ to that of the resistive diffusion term $\eta \nabla ^2 \boldsymbol {B}$ in the MHD induction equation.

The kinematic viscosity $\nu$ and resistivity $\eta$ are transport coefficients that are rigorously defined by an extension of the Chapman–Enskog procedure for magnetohydrodynamic systems (Spitzer Reference Spitzer1962; Grad Reference Grad1963; Braginskii Reference Braginskii1965) in the limit of small but finite mean free path relative to the gradient scale, $\lambda _{{\rm mfp}}/L \ll 1$. As illustrated by figure 11, however, for many space and astrophysical plasma systems of interest, the scales at which energy is removed from the turbulence (beginning at the small-scale end of the inertial range) are weakly collisional, with $\lambda _{{\rm mfp}}/L \gtrsim 1$ or $\lambda _{{\rm mfp}}/L \gg 1$. In these complementary limits of moderate to large mean free path, relevant to space and astrophysical plasmas, the procedure to calculate the viscosity and resistivity ceases to be valid. In fact, the standard Laplacian viscosity and resistivity terms are likely to be poor approximations to the physical mechanisms that are believed to remove energy from weakly collisional plasma turbulence (see § 5.1 for a list of proposed damping mechanisms). Since it is not possible to define clearly the viscosity or resistivity in a weakly collisional system, it is also not possible to define the Reynolds number or magnetic Reynolds number in weakly collisional, kinetic plasma turbulence. Thus, we are forced to choose alternative dimensionless numbers that characterize the plasma and turbulence in weakly collisional systems, so one of the primary aims of this paper is to propose the specific set of fundamental dimensionless plasma and turbulence parameters that are suitable for characterizing kinetic plasma turbulence, presented here in table 3.

3.1. Exploiting turbulent scaling theories to reduce the number of fundamental dimensionless parameters

Given the large number of fundamental dimensionless parameters that characterize weakly collisional plasma turbulence, shown in table 3, it is imperative to reduce the number of important parameters to make progress in developing predictive models of turbulent heating. Adopting idealized approximations, such as maintaining isotropic equilibrium velocity distributions for the plasma, enables one to eliminate some of the parameters in table 3. Another approach is to exploit the MHD turbulent scaling theories described in § 3 to combine several parameters into a new single parameter, an approach presented here.

In general, the driving of the plasma turbulence depends on at least three independent parameters: (i) the perpendicular wavenumber at which the turbulence is driven $k_{\perp 0} \rho _i$, (ii) the anisotropy of the driven fluctuations with respect to the local magnetic field direction $k_{\parallel 0}/k_{\perp 0}$ and (iii) the nonlinearity parameter which characterizes the strength of the driving, $\chi _0 \equiv k_{\perp 0} \delta B_{\perp 0} \theta _0/(k_{\parallel 0} B_0)$. Recall that quantities evaluated at the driving scale are given the subscript ‘0’. Fortunately, for turbulence driven at large scales (relative to the ion kinetic scales), as is typical in most space or astrophysical plasma environments, MHD scaling theory can be used to predict the properties of the turbulence upon reaching the kinetic scales, where kinetic plasma processes act to dissipate the turbulence. In fact, as detailed here, it is possible to combine these three dimensionless parameters into a new, single dimensionless parameter: the isotropic driving wavenumber, $k_0 \rho _i$.

MHD turbulent scaling theories can be used, for a given set of driving parameters ($k_{\perp 0} \rho _i$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$), to predict the properties of the turbulent fluctuations (e.g. anisotropy, amplitude) upon reaching the ion kinetic scales, here taken to be the perpendicular wavenumber $k_\perp \rho _i = 1$. The isotropic driving wavenumber $k_0 \rho _i$ is the unique wavenumber that produces the same conditions at $k_\perp \rho _i = 1$ as if the turbulence were driven isotropically at that wavenumber with $k_{\perp 0} = k_{\parallel 0} \equiv k_0$ and in a condition of critical balance with $\chi _0=1$.

The key procedure for determining the isotropic driving wavenumber $k_0 \rho _i$ is illustrated in figure 2. Panel (a) shows the perpendicular wavenumber magnetic-energy spectrum $E_{B_\perp }(k_\perp )$ driven weakly with $\chi _0=0.1$ at perpendicular wavenumber $k_{\perp 0} \rho _i=10^{-5}$, where a weak turbulent cascade leads to a steep weak MHD turbulence spectrum with $E_{B_\perp } \propto k_\perp ^{-2}$ (Sridhar & Goldreich Reference Sridhar and Goldreich1994; Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000; Boldyrev & Perez Reference Boldyrev and Perez2009), followed by a flatter strong MHD turbulence spectrum with $E_{B_\perp } \propto k_\perp ^{-(5+\alpha )/(3 +\alpha )}$ (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Boldyrev Reference Boldyrev2006). The cascade reaches the perpendicular ion kinetic scale $k_\perp \rho _i = 1$, and at higher wavenumbers the cascade continues as strong kinetic Alfvén wave (KAW) turbulence with a significantly steeper magnetic-energy spectrum with $E_{B_\perp } \propto k_\perp ^{-2.8}$ (Alexandrova et al. Reference Alexandrova, Saur, Lacombe, Mangeney, Mitchell, Schwartz and Robert2009; Kiyani et al. Reference Kiyani, Chapman, Khotyaintsev, Dunlop and Sahraoui2009; Sahraoui et al. Reference Sahraoui, Goldstein, Robert and Khotyaintsev2009; Chen et al. Reference Chen, Horbury, Schekochihin, Wicks, Alexandrova and Mitchell2010; Howes et al. Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011b; Sahraoui et al. Reference Sahraoui, Huang, Belmont, Goldstein, Rétino, Robert and De Patoul2013). When the collisionless damping of the KAW turbulent cascade begins to diminish the rate of nonlinear energy cascade, the spectrum begins to fall off exponentially in the weakly dissipating KAW turbulence (WDKT) regime (Howes, Tenbarge & Dorland Reference Howes, Tenbarge and Dorland2011a), ultimately terminating the turbulent cascade at scales $k_\perp \rho _e \gtrsim 1$.

Figure 2. Diagram illustrating how to determine the isotropic driving wavenumber $k_0\rho _i$ in terms of the perpendicular driving wavenumber $k_{\perp 0} \rho _i$, the driving anisotropy $k_{\parallel 0}/k_{\perp 0}$ and the nonlinearity parameter of the driving, $\chi _0=k_{\perp 0} \delta B_\perp \theta _0/(k_{\parallel 0} B_0)$.

In figure 2(b) is presented the characteristic path of the anisotropic cascade of energy through $(k_\perp,k_\parallel )$ wavevector space for the different regimes of the turbulent cascade. There is no parallel cascade in the weak MHD turbulence regime (Sridhar & Goldreich Reference Sridhar and Goldreich1994; Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000), followed by a scale-dependent anisotropy given by $k_\parallel \propto k_\perp ^{2/(3+\alpha )}$ in the strong MHD turbulence regime (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Boldyrev Reference Boldyrev2006). The turbulence transitions to the strong KAW regime at $k_\perp \rho _i \sim 1$, with the anisotropy becoming stronger at $k_\perp \rho _i > 1$, scaling as $k_\parallel \propto k_\perp ^{1/3}$ (Cho & Lazarian Reference Cho and Lazarian2004; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a, Reference Howes, Tenbarge and Dorland2011a; TenBarge & Howes Reference TenBarge and Howes2012; TenBarge, Howes & Dorland Reference TenBarge, Howes and Dorland2013). As the turbulence weakens due to collisionless damping of turbulent fluctuations in the WDKT regime, the parallel cascade is predicted to cease once again (Howes et al. Reference Howes, Tenbarge and Dorland2011a).

The diagrams in figure 2 provide guidance to understand how MHD turbulence scaling theories can be used to combine ($k_{\perp 0} \rho _i$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$) into a single new parameter, the isotropic driving wavenumber, $k_0 \rho _i$. If the turbulence is driven at critical balance ($\chi _0=1$) and isotropically ($k_{\parallel 0}/k_{\perp 0}=1$), then the isotropic driving wavenumber is simply the same as the perpendicular driving wavenumber, $k_0 \rho _i=k_{\perp 0} \rho _i$. Therefore, it is necessary to use MHD turbulence scaling theory to determine $k_0 \rho _i$ only in the case that the turbulence is driven weakly with $\chi _0<1$ or anisotropically with $k_{\parallel 0}/k_{\perp 0}\ne 1$.

For weakly driven turbulence with $\chi _0<1$, MHD turbulence theory predicts that the parallel cascade is inhibited and that the perpendicular wavenumber magnetic-energy spectrum of the turbulence scales as $E_{B_\perp } \propto (\delta B_\perp )^2 /k_\perp \propto k_\perp ^{-2}$. Therefore, in terms of the driving amplitude $\delta B_{\perp 0}$ and perpendicular driving wavenumber $k_{\perp 0}$, the amplitude of the turbulence as a function of $k_\perp$ is given by $\delta B_\perp (k_\perp )=\delta B_{\perp 0} (k_{\perp 0}/k_\perp )^{1/2}$. Given the lack of parallel cascade in weak MHD turbulence, we fix $k_\parallel = k_{\parallel 0}$, and obtain the following scalingFootnote ⁸ for the nonlinearity parameter as a function of $k_\perp$:

(3.8)\begin{equation} \chi (k_\perp) = \chi_0 \left(\frac{k_{{\perp}} \rho_i}{k_{{\perp} 0} \rho_i}\right)^{1/2}. \end{equation}

Magnetohydrodynamic turbulence theory predicts that the nonlinearity parameter will increase through the weak turbulent cascade until it reaches a value of unity, $\chi \sim 1$, at which point the turbulence has reached a state of critical balance in the strong turbulence regime. From that point, theory predicts that an anisotropic cascade will ensue (transferring energy to smaller scales more rapidly in the perpendicular than the parallel direction), with the parallel and perpendicular cascades balanced to maintain a nonlinearity parameter of unity, $\chi \sim 1$. Therefore, the key point in the cascade is the perpendicular wavenumber at which the turbulence first reaches $\chi =1$, given by

(3.9)\begin{equation} k_{{\perp}} \rho_i|_{\chi=1}= \frac{1}{\chi_0^2} k_{{\perp} 0} \rho_i .\end{equation}

An example of this weakly turbulent, strictly perpendicular cascade in $(k_\perp,k_\parallel )$ space is depicted in figure 2(b), where turbulence driven weakly with $\chi _0 =0.1$, cascades to $k_{\perp } \rho _i= 100 k_{\perp 0} \rho _i$, at which point the turbulence becomes strong, with $\chi =1$.

To obtain the isotropic driving wavenumber $k_0 \rho _i$ from this point in wavevector space $(k_{\perp } \rho _i|_{\chi =1}, k_{\parallel 0})$, we need to extend the cascade backward along the critical balance line (red dashed lines) until it intersects with isotropy $k_\parallel =k_\perp$ (dotted line). From the MHD scaling theories in § 3, for a critically balanced, strong turbulent cascade that is driven strongly and isotropically at $k_0\rho _i$, the anisotropy as a function of $k_\perp$ scales as

(3.10)\begin{equation} \frac{k_\parallel}{k_\perp} = \left(\frac{k_0 \rho_i}{k_\perp \rho_i}\right)^{{(1+\alpha)}/{(3 +\alpha)}}. \end{equation}

The anisotropy at the point $(k_{\perp } \rho _i|_{\chi =1}, k_{\parallel 0})$ is given by $k_\parallel /k_\perp = (k_{\parallel 0}/k_{\perp 0}) \chi _0^2$. Substituting this expression for $k_\parallel /k_\perp$ on the left-hand side of (3.10) and substituting (3.9) into the denominator on the right-hand side of (3.10), we can solve for the isotropic driving wavenumber $k_0 \rho _i$, obtaining the results

(3.11)\begin{equation} k_0\rho_i = k_{{\perp} 0} \rho_i \left(\frac{k_{{\parallel} 0}}{k_{{\perp} 0}}\right)^{{(3+\alpha)}/{(1+\alpha)}} \chi_0^{{4}/{(1+\alpha)}}.\end{equation}

In figure 2(b), the isotropic driving wavenumber $k_0 \rho _i$ (red circle) is the shown where the critical balance line (solid black and red dashed lines) intersects the isotropy line, $k_\parallel =k_\perp$ (dotted line). Note that this reduction of the three driving parameters to $k_0\rho _i$ is valid only in the limit that the turbulent cascade becomes strong before it reaches the ion kinetic scales at $k_\perp \rho _i \sim 1$.

The relation (3.11) provides the means to combine the three turbulence driving parameters ($k_{\perp 0} \rho _i$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$) into the single new effective parameter, the isotropic driving wavenumber $k_0 \rho _i$. The properties of the turbulent fluctuations at perpendicular scales smaller than the point where the turbulent cascade becomes strong, or $k_\perp \rho _i> k_{\perp } \rho _i|_{\chi =1}$, will be the same for the original case with driving parameterized by ($k_{\perp 0} \rho _i$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$) as for turbulence driven strongly and isotropically at $k_0 \rho _i$. In the development of a turbulent heating model, it is the properties of the turbulent fluctuations at the small, kinetic length scales that will determine the kinetic damping mechanism and the resulting particle energization, so reducing three turbulence driving parameters to a single parameter substantially reduces the dimensionality of the parameter space to be modelled.

It is worthwhile noting that the evolution of the dynamic alignment angle with perpendicular wavenumber $\theta (k_\perp )$ is not well modelled by $k_0 \rho _i$. The development of current sheets at very small scales due to dynamic alignment, and the possibility for the interruption of the turbulent cascade by magnetic reconnection (Boldyrev & Loureiro Reference Boldyrev and Loureiro2017; Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a, Reference Loureiro and Boldyrevb; Mallet et al. Reference Mallet, Schekochihin and Chandran2017a, Reference Mallet, Schekochihin and Chandranb; Walker et al. Reference Walker, Boldyrev and Loureiro2018), will depend on this scaling of $\theta (k_\perp )$. But, from a practical point of view, this issue in modelling $\theta (k_\perp )$ while reducing ($k_{\perp 0} \rho _i$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$) to $k_0 \rho _i$ is likely to be a problem only in rare cases, for the following reasons. First, the scaling of dynamic alignment angle $\theta \propto k_\perp ^{-1/4}$ in B06 is rather weak, scaling as the $-1/4$ power, so only when the MHD inertial range spans many orders of magnitude in the perpendicular length scale will the development of current sheets likely have a significant impact on the turbulent dynamics. Second, dynamic alignment appears to arise only in strong MHD turbulence, so it seems reasonable to take $\theta _0 \sim 1$ rad (roughly equivalent to isotropy in the perpendicular plane) at the perpendicular scale $k_{\perp } \rho _i|_{\chi =1}$ where the turbulence becomes strong. One may then estimate the evolution of that dynamic alignment angle $\theta (k_\perp )$ from that point for the cases in which it impacts the evolution. Finally, although the development of current sheets in plasma turbulence is well established (Matthaeus & Montgomery Reference Matthaeus and Montgomery1980; Meneguzzi et al. Reference Meneguzzi, Frisch and Pouquet1981; Borovsky Reference Borovsky2008; Uritsky et al. Reference Uritsky, Pouquet, Rosenberg, Mininni and Donovan2010; Osman et al. Reference Osman, Matthaeus, Greco and Servidio2011; Zhdankin et al. Reference Zhdankin, Boldyrev, Mason and Perez2012, Reference Zhdankin, Uzdensky, Perez and Boldyrev2013), it is not clear that the scaling theory in B06 indeed properly describes this development; an alternative explanation for the development of current sheets arising self-consistently from Alfvén wave collisions exists (Howes Reference Howes2015a, Reference Howes2016; Verniero et al. Reference Verniero, Howes and Klein2018), although a specific scaling for this alternative mechanism has not yet been proposed.

3.2. Efficiently modelling the small-scale end of the turbulent inertial range

The multiscale physics of astrophysical plasma turbulence presents a particular challenge for numerical modelling that covers the entire dynamic range of the turbulent cascade while simultaneously capturing the kinetic mechanisms that govern the damping of the turbulence under weakly collisional plasma conditions, in particular since the physics of plasma turbulence is inherently three-dimensional (Howes Reference Howes2015b). For the development of turbulent heating models, it is critical to investigate the physical mechanisms that remove energy from the turbulence and energize the plasma particles. Fortunately, since the inertial range of turbulence is defined as the range of scales over which the turbulent energy cascades to smaller scales but the dissipation is negligible, as depicted in figure 1(a), it is not strictly necessary to model the entire inertial range in order to capture the physics of the turbulent dissipationFootnote ⁹ .

One less computationally costly approach to capture numerically the mechanisms of turbulent damping is to choose to model directly only the small-scale end of the inertial range, preserving as much of the resolved dynamic range of the simulation as possible to capture the dissipation mechanisms. To do so, one takes advantage of the scaling relations for turbulence in the MHD regime $k_\perp \rho _i \ll 1$, summarized in § 3, to ‘drive’ the turbulence at the domain scale of the simulation with turbulent fluctuations that have the properties dictated by the scaling theories.

As an example of how to set up simulations in this way, consider turbulence physically driven strongly and isotropically at an isotropic driving wavenumber of $k_0\rho _i =10^{-4}$, as depicted in figure 3. For a hydrogenic plasma of protons and electrons with a realistic mass ratio $m_i/m_e=1836$, the turbulent cascade extends at least down to the perpendicular scales of the electron Larmor radius $k_\perp \rho _e \sim 1$ (Sahraoui et al. Reference Sahraoui, Goldstein, Belmont, Canu and Rezeau2010; Howes et al. Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011b; TenBarge et al. Reference TenBarge, Howes and Dorland2013; Kiyani, Osman & Chapman Reference Kiyani, Osman and Chapman2015). Thus, a dynamic range of $L_0/\rho _e =k_{\perp \max }/k_0 = [(k_{\perp \max } \rho _e)/(k_0 \rho _i)] (\rho _i/\rho _e)$ $\propto (k_0 \rho _i)^{-1}(T_i/T_e)^{1/2}(m_i/m_e)^{1/2} \sim 4.3 \times 10^5$ in each of the three spatial dimensions is necessary when $T_i/T_e \sim 1$, well beyond the capabilities of any kinetic simulation code. If the simulation code can resolve a dynamic range of $10^3$ in each spatial dimension, we can choose a perpendicular domain scale in the middle of the physical inertial range at $k_{\perp D} \rho _i = 10^{-2}$ and use a reduced mass ratio $m_i/m_e=100$ so that $k_\perp \rho _e=1$ corresponds to $k_\perp \rho _i=10$ for $T_i/T_e=1$. Thus, the resolved perpendicular length scales of the simulation span $10^{-2} \le k_\perp \rho _i \le 10$, corresponding to $10^{-3} \le k_\perp \rho _e \le 1$ in terms of the electron Larmor radius $\rho _e$, as shown in figure 3.

Figure 3. Modelling the small-scale end of the turbulent inertial range using MHD scaling relations, choosing B06 scaling with $\alpha =1$, $m_i/m_e=100$, $k_0\rho _i =10^{-4}$, $k_{\perp D} \rho _i = 10^{-2}$.

To drive the turbulence at the perpendicular simulation domain scale $L_\perp =2 {\rm \pi}/k_{\perp D}$ with $k_{\perp D} \rho _i = 10^{-2}$, corresponding to domain scale in the middle of the turbulent inertial range, we use the scaling relations (3.1)–(3.4) to determine the appropriate turbulent scale-dependent anisotropy and amplitude of the fluctuations at the numerical domain scale for a physical inertial range driven strongly ($\chi _0=1$) and isotropically ($k_{\perp 0}= k_{\parallel 0}=k_{i 0}\equiv k_0$ and $\theta _0=1$ rad) at $k_0\rho _i =10^{-4}$. For this example, we will adopt the B06 scaling, choosing $\alpha =1$. The wavevector anisotropy and amplitude of the fluctuations at the domain scale can be expressed using scaling relations (3.1)–(3.4) to yield

(3.12)\begin{gather} k_{{\parallel} D} \rho_i = k_0 \rho_i \left( \frac{k_{{\perp} D}\rho_i}{k_0 \rho_i}\right)^{1/2} \end{gather}

(3.13)\begin{gather}k_{iD}\rho_i= k_0 \rho_i \left( \frac{k_{{\perp} D}\rho_i}{k_0 \rho_i}\right)^{3/4} \end{gather}

(3.14)\begin{gather}\theta_D = \theta_0 \left( \frac{k_{{\perp} D}\rho_i}{k_0 \rho_i}\right)^{{-}1/4} \end{gather}

(3.15)\begin{gather}\frac{\delta B_{{\perp} D} }{B_0} = \left( \frac{k_{{\perp} D}\rho_i}{k_0 \rho_i}\right)^{{-}1/4}. \end{gather}

The properties of the domain-scale fluctuations at $k_{\perp D} \rho _i = 10^{-2}$ determined by these scaling relations are illustrated in figure 3 for a turbulent cascade parameterized by dimensionless parameters $k_0\rho _i =10^{-4}$ and $\theta _0=1$ rad. The normalized value of the energy spectrum at $k_{\perp D} \rho _i = 10^{-2}$ is $(k_0/B_0)^2 E_{B_\perp }(k_\perp ) =10^{-3}$, consistent with the solution $\delta B_{\perp D}/B_0=0.32$ of (3.15) at the domain scale, given by the red circle in figure 3(a). The three-dimensional anisotropy of the wavevector of the turbulent fluctuations at the perpendicular domain scale $k_{\perp D} \rho _i = 10^{-2}$ is characterized by: (i) $k_{\parallel D} \rho _i =10^{-3}$ from (3.12), given by the blue circle in figure 3(b); and (ii) $k_{i D} \rho _i =3.2 \times 10^{-3}$ from (3.13), given by the green circle in figure 3(b). The dynamic alignment of the velocity and magnetic field fluctuations at the domain scale is given by $\theta _D=0.32$ rad from (3.14). Note that for these parameters, the scalings lead to a constant nonlinearity parameter as a function of perpendicular wavenumber $\chi (k_\perp ) = \chi _0=1$ at all scales $k_\perp \rho _i<1$, as dictated by the conjecture of critical balance in strong turbulence (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Mallet et al. Reference Mallet, Schekochihin and Chandran2015). These scaling calculations determine the nature of the turbulent fluctuations at the numerical domain scale, in the middle of the physical inertial range being modelled.

Since the largest perpendicular scale of the resolved dynamic range of the simulation at $k_{\perp D} \rho _i$ is constantly being fed turbulent cascade energy mediated by nonlinear interactions among turbulent fluctuations at scales slightly larger than the domain scale with $k_\perp < k_{\perp D}$, it is imperative that these turbulence simulations are driven, rather than simply simulating the decay of the energy turbulent fluctuations initialized at the domain scale. The nonlinear interactions among perpendicularly polarized, counterpropagating Alfvén waves drive the cascade of energy in Alfvénic plasma turbulence (Kraichnan Reference Kraichnan1965; Sridhar & Goldreich Reference Sridhar and Goldreich1994; Goldreich & Sridhar Reference Goldreich and Sridhar1995; Ng & Bhattacharjee Reference Ng and Bhattacharjee1996; Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000; Howes et al. Reference Howes, Drake, Nielson, Carter, Kletzing and Skiff2012b; Drake et al. Reference Drake, Schroeder, Howes, Kletzing, Skiff, Carter and Auerbach2013; Howes & Nielson Reference Howes and Nielson2013; Howes et al. Reference Howes, Nielson, Drake, Schroeder, Skiff, Kletzing and Carter2013; Nielson et al. Reference Nielson, Howes and Dorland2013; Howes Reference Howes2016; Verniero & Howes Reference Verniero and Howes2018; Verniero et al. Reference Verniero, Howes and Klein2018). Exploiting this physical insight, driving counterpropagating Alfvén waves, polarized in both dimensions perpendicular to the equilibrium magnetic field, has been shown to generate successfully a steady-state plasma turbulent cascade (Howes et al. Reference Howes, Dorland, Cowley, Hammett, Quataert, Schekochihin and Tatsuno2008b, Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011b; TenBarge et al. Reference TenBarge, Howes and Dorland2013; Told et al. Reference Told, Jenko, TenBarge, Howes and Hammett2015; Verniero & Howes Reference Verniero and Howes2018; Verniero et al. Reference Verniero, Howes and Klein2018). A finite-time correlated driving, such as an antenna driven by a Langevin equation (TenBarge et al. Reference TenBarge, Howes, Dorland and Hammett2014), mimics the unsteady nature of the turbulent cascade from larger scales expected of turbulent driving at a chosen scale within the turbulent inertial range.

To set up the simulation to describe optimally the domain-scale fluctuations, it has become common practice to utilize a numerical domain that is elongated along the direction of the equilibrium magnetic field, consistent with the anisotropic scaling of turbulent fluctuations at small scales, $k_\parallel /k_\perp \ll 1$. For the example shown in figure 3, individual Alfvénic fluctuations are predicted to have a perpendicular plane anisotropy given by $k_{iD}/k_{\perp D}=0.32$. But since these Alfvénic fluctuations can be polarized in either of the directions in the perpendicular plane, it is necessary for both perpendicular directions to be large enough to resolve the intermediate scale (the larger of the two perpendicular scales since $k_{iD}/k_{\perp D} < 1$). Thus, the simulation domain should be defined by $L_\parallel \times L_i^2$, where $L_\parallel = 2 {\rm \pi}\rho _i (k_{\parallel D} \rho _i)^{-1}$ and $L_i = 2 {\rm \pi}\rho _i (k_{i D} \rho _i)^{-1}$. Within this domain, counterpropagating Alfvénic fluctuations should be driven, polarized in both perpendicular directions, with wavevectors characterized by $(k_{\perp D},k_{i D},k_{\parallel D})$, with dynamic alignment characterized by $\theta _D$, and with amplitudes $\delta B_{\perp D}$. Note that setting up properly dynamically aligned Alfvénic fluctuations at the domain scale is a challenging task that will depend on the particular numerical implementation of the simulation code. Nonetheless, by driving the simulation appropriately within the inertial range, this approach enables the efficient simulation of the small-scale end of the inertial range consistent with the scaling properties of MHD turbulence, enabling kinetic simulations to focus on the physical mechanisms damping the turbulent fluctuations.

3.3. A reduced parameter space for kinetic plasma turbulence

Although the physics of plasma turbulence, its dissipation, and the resulting heating of the plasma species formally depends on all ten of the general plasma and turbulence parameters $(\beta _{\parallel,i},\tau _\parallel,A_i,A_e,k_{\parallel 0} \lambda _{{\rm mfp},e}; k_{\perp 0} \rho _i,k_{\parallel 0}/k_{\perp 0},\chi _0,Z_0^+/Z_0^-,E_{{\rm comp}}/E_{{\rm inc}})$ listed in table 3, focusing initially on the development of completely general turbulent heating models on a ten-dimensional parameter space is unlikely to be successful. Instead, first developing predictive turbulent heating models over a reduced parameter space is almost certainly a better strategy. As an initial exploration of Alfvénic turbulence and its kinetic damping mechanisms, we choose to focus on the simplified limit of isotropic temperatures, modifying our approach as needed for cases that violate these conditions.

For many of the calculations in this study, a reduced isotropic temperature case is adopted, assuming that the equilibrium velocity distribution for each species remains well described as an isotropic Maxwellian distribution. Note that kinetic temperature anisotropy instabilities generally bound the species temperature anisotropies in the vicinity of isotropy $A_s \sim 1$ (Kasper, Lazarus & Gary Reference Kasper, Lazarus and Gary2002; Hellinger et al. Reference Hellinger, Trávníček, Kasper and Lazarus2006; Štverák et al. Reference Štverák, Trávníček, Maksimovic, Marsch, Fazakerley and Scime2008; Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009; Štverák et al. Reference Štverák, Maksimovic, Trávníček, Marsch, Fazakerley and Scime2009; Verscharen et al. Reference Verscharen, Chandran, Klein and Quataert2016), although if these instabilities are triggered, particularly at $\beta _i \gg 1$, they can impact both the linear wave physics at play in the turbulent dynamics of critically balanced turbulence (Squire, Quataert & Schekochihin Reference Squire, Quataert and Schekochihin2016; Squire et al. Reference Squire, Kunz, Quataert and Schekochihin2017a; Squire, Schekochihin & Quataert Reference Squire, Schekochihin and Quataert2017b) and mediate the non-local transfer of energy directly from the large driving scales to the kinetic length scales (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023).

To define a reduced parameter space, we first restrict the temperature anisotropies for ions and electrons to unity, $A_i=1$ and $A_e=1$, so that $T_{\perp,i}=T_{\parallel,i}\equiv T_i$ and $T_{\perp,e}=T_{\parallel,e}\equiv T_e$. With this idealization, the ion plasma beta and ion-to-electron temperature ratio are redefined in terms of these isotropic temperatures $T_i$ and $T_e$, yielding the isotropic parameters $\beta _i=8 {\rm \pi}n_i T_{i}/B^2$ and $\tau =T_{i}/T_{e}$. In addition, for this reduced parameter space, we will also assume that the turbulent inertial range is sufficiently large, or the driving is sufficiently strong, that the dynamics at the perpendicular ion kinetic length scales $k_\perp \rho _i \sim 1$ is in a state of strong plasma turbulence with a nonlinearity parameter $\chi \sim 1$. Under these conditions, we can replace the three turbulent driving parameters $(k_{\perp 0} \rho _i,k_{\parallel 0}/k_{\perp 0},\chi _0)$ by the single isotropic driving wavenumber $k_0\rho _i$, as detailed in § 3.1.

We will also make two other simplifying assumptions for this reduced, isotropic temperature case. We will assume that the turbulence is balanced, with approximately equal wave energy fluxes up and down the local magnetic field $Z_0^+/Z_0^- \sim 1$, which is a necessary condition to apply the MHD turbulence scaling relations given by (3.1)–(3.4). Additionally, based on direct spacecraft measurements of turbulence in the heliosphere which find that incompressible Alfvénic fluctuations appear to dominate turbulent dynamics in space plasmas (Belcher & Davis Reference Belcher and Davis1971; Tu & Marsch Reference Tu and Marsch1995; Alexandrova et al. Reference Alexandrova, Carbone, Veltri and Sorriso-Valvo2008; Bruno & Carbone Reference Bruno and Carbone2013), as well as theoretical considerations that suggest that the dynamics of Alfvén waves dominates the turbulent cascade (Lithwick & Goldreich Reference Lithwick and Goldreich2001; Maron & Goldreich Reference Maron and Goldreich2001; Cho & Lazarian Reference Cho and Lazarian2003; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Howes Reference Howes2015a), we will also assume a negligible amount of energy in compressible turbulent fluctuations at the driving scale, $E_{{\rm comp}}/E_{{\rm inc}} \ll 1$.

Therefore, we propose the isotropic temperature case, a reduced parameter space for weakly collisional plasma turbulence in the limit of strong, balanced, Alfvénic turbulence, with a governing set of three dimensionless parameters, $(\beta _i,\tau ; k_{0} \rho _i)$, as shown in table 4. Although the isotropic temperature case eliminates a number of important physical effects such as temperature anisotropy instabilities, imbalanced turbulent wave energy fluxes and compressible fluctuations, it nonetheless represents an excellent starting point for the development of refined turbulent heating models. Such models can produce theoretical predictions that can be compared with kinetic numerical simulations and spacecraft observations of weakly collisional plasma turbulence, exemplified by the phase diagram for turbulent damping mechanisms shown in figure 12.

Table 4. Isotropic temperature case: reduced set of fundamental dimensionless parameters for isotropic Maxwellian equilibrium velocity distributions with $A_i=A_e=1$, collisionless conditions $k_{\parallel 0} \lambda _{{\rm mfp},e} \gg 1$, strong turbulence assumed at the ion kinetic scales, balanced turbulence $Z_0^+/Z_0^-\sim 1$, and negligible energy in compressible fluctuations $E_{{\rm comp}}/E_{{\rm inc}} \ll 1$. This results in a reduced set of plasma and turbulence parameters $(\beta _i,\tau ; k_{0} \rho _i)$.

In terms of the fundamental parameters of the isotropic temperature case, the MHD turbulent scaling relations (3.1)–(3.4) simplify to the following easy-to-use forms:

(3.16)\begin{gather} k_\parallel \rho_i = (k_0 \rho_i)^{(1+\alpha)/(3+\alpha)} (k_\perp \rho_i)^{2/(3+\alpha)} \end{gather}

(3.17)\begin{gather}k_i \rho_i = (k_0 \rho_i)^{\alpha/(3+\alpha)} (k_\perp \rho_i)^{3/(3+\alpha)} \end{gather}

(3.18)\begin{gather}\theta = \theta_0\left(\frac{k_0 \rho_i}{k_\perp \rho_i}\right)^{\alpha/(3+\alpha)} \end{gather}

(3.19)\begin{gather}\frac{\delta B_{{\perp} } }{B_0} = \left(\frac{k_0 \rho_i}{k_\perp \rho_i} \right)^{1/(3+\alpha)}. \end{gather}

5.1. Proposed damping mechanisms for kinetic turbulence

In weakly collisional space and astrophysical plasmas, the proposed kinetic damping mechanisms that remove energy from the turbulent fluctuations and transfer it to the particle species fall into three broad categories: (i) resonant wave-particle interactions, (ii) non-resonant wave–particle interactions and (iii) damping occurring in coherent structures.

Resonant mechanisms for the damping of the turbulent fluctuations in a magnetized plasma transfer energy to particles in regions of the velocity distribution that satisfy the resonance condition

(5.1)\begin{equation} \omega- k_\parallel v_\parallel{-} n \varOmega_s=0, \end{equation}

where $\varOmega _s$ is the cyclotron frequency for species $s$ and $|n|>1$ indicates the harmonics of the cyclotron frequency. The $n=0$ resonance is often denoted the Landau resonance and includes two damping mechanisms: (i) Landau damping is mediated by the parallel component of the electric field $E_\parallel$ doing work on the charged particles (Landau Reference Landau1946; Leamon et al. Reference Leamon, Matthaeus, Smith and Wong1998a, Reference Leamon, Smith, Ness, Matthaeus and Wongb; Quataert Reference Quataert1998; Leamon et al. Reference Leamon, Smith, Ness and Wong1999; Quataert & Gruzinov Reference Quataert and Gruzinov1999; Leamon et al. Reference Leamon, Matthaeus, Smith, Zank, Mullan and Oughton2000; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; TenBarge & Howes Reference TenBarge and Howes2013; Howes Reference Howes2015a; Li et al. Reference Li, Howes, Klein and TenBarge2016; Chen, Klein & Howes Reference Chen, Klein and Howes2019; Afshari et al. Reference Afshari, Howes, Kletzing, Hartley and Boardsen2021; Zhou, Liu & Loureiro Reference Zhou, Liu and Loureiro2023); and (ii) transit-time damping, also known as Barnes damping, is mediated by the magnetic mirror force associated with magnetic field magnitude fluctuations doing work on the magnetic moment of the charged particle gyromotion (Barnes Reference Barnes1966; Quataert Reference Quataert1998; Quataert & Gruzinov Reference Quataert and Gruzinov1999; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a; Howes, Huang & Felix Reference Howes, Huang and Felix2024). Both of these Landau-resonant collisionless damping mechanisms energize particles in the degree of freedom parallel to the local magnetic field. The $n \ne 0$ resonances yield cyclotron damping of the turbulent fluctuations by $E_\perp$, generally energizing particles in the two perpendicular degrees of freedom (Coleman Reference Coleman1968; Denskat, Beinroth & Neubauer Reference Denskat, Beinroth and Neubauer1983; Isenberg & Hollweg Reference Isenberg and Hollweg1983; Goldstein, Roberts & Fitch Reference Goldstein, Roberts and Fitch1994; Leamon et al. Reference Leamon, Matthaeus, Smith and Wong1998a; Gary Reference Gary1999; Isenberg, Lee & Hollweg Reference Isenberg, Lee and Hollweg2001; Hollweg & Isenberg Reference Hollweg and Isenberg2002; Isenberg & Vasquez Reference Isenberg and Vasquez2019; Afshari et al. Reference Afshari, Howes, Shuster, Klein, McGinnis, Martinovic, Boardsen, Brown and Huang2024).

Non-resonant damping mechanisms can remove energy from turbulent fluctuations through wave–particle interactions that do not require satisfying a particular resonance condition. Stochastic heating can arise when sufficiently large-amplitude turbulent electromagnetic fluctuations can effectively scatter the otherwise organized gyromotion of particles, leading to a heating of the particles in the perpendicular degrees of freedom (McChesney, Stern & Bellan Reference McChesney, Stern and Bellan1987; Chen, Lin & White Reference Chen, Lin and White2001; Johnson & Cheng Reference Johnson and Cheng2001; White, Chen & Lin Reference White, Chen and Lin2002; Voitenko & Goossens Reference Voitenko and Goossens2004; Bourouaine, Marsch & Vocks Reference Bourouaine, Marsch and Vocks2008; Chandran Reference Chandran2010; Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010, Reference Chandran, Dennis, Quataert and Bale2011; Bourouaine & Chandran Reference Bourouaine and Chandran2013; Chandran et al. Reference Chandran, Verscharen, Quataert, Kasper, Isenberg and Bourouaine2013; Klein & Chandran Reference Klein and Chandran2016; Arzamasskiy et al. Reference Arzamasskiy, Kunz, Chand Ran and Quataert2019; Martinović et al. Reference Martinović, Klein, Kasper, Case, Korreck, Larson, Livi, Stevens, Whittlesey and Chandran2020; Cerri, Arzamasskiy & Kunz Reference Cerri, Arzamasskiy and Kunz2021). Magnetic pumping energizes particles through a combination of two effects: (i) particles undergo the reversible, double-adiabatic evolution of their perpendicular and parallel velocities in a magnetic field with a time varying magnitude; and (ii) those particles are scattered by collisions during this evolution, introducing irreversibility and leading to a small net transfer of energy to the particles over a full oscillation cycle (Spitzer & Witten Reference Spitzer and Witten1953; Berger et al. Reference Berger, Newcomb, Dawson, Frieman, Kulsrud and Lenard1958; Lichko et al. Reference Lichko, Egedal, Daughton and Kasper2017; Lichko & Egedal Reference Lichko and Egedal2020; Montag & Howes Reference Montag and Howes2022). One final proposed non-resonant damping mechanism is the recently identified kinetic ‘viscous’ heating mediated by the effective collisionality associated with temperature anisotropy instabilities driven by large-scale fluctuations of the magnetic field magnitude in high beta plasmas (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023).

The final category involves the removal of turbulent fluctuation energy in coherent structures, such as particle energization arising from collisionless magnetic reconnection in current sheets that are found to arise naturally in plasma turbulence (Ambrosiano et al. Reference Ambrosiano, Matthaeus, Goldstein and Plante1988; Dmitruk, Matthaeus & Seenu Reference Dmitruk, Matthaeus and Seenu2004; Markovskii & Vasquez Reference Markovskii and Vasquez2011; Matthaeus & Velli Reference Matthaeus and Velli2011; Osman et al. Reference Osman, Matthaeus, Greco and Servidio2011; Servidio et al. Reference Servidio, Greco, Matthaeus, Osman and Dmitruk2011; Osman et al. Reference Osman, Matthaeus, Hnat and Chapman2012a, Reference Osman, Matthaeus, Wan and Rappazzob; Wan et al. Reference Wan, Matthaeus, Karimabadi, Roytershteyn, Shay, Wu, Daughton, Loring and Chapman2012; Karimabadi et al. Reference Karimabadi, Roytershteyn, Wan, Matthaeus, Daughton, Wu, Shay, Loring, Borovsky and Leonardis2013; Zhdankin et al. Reference Zhdankin, Uzdensky, Perez and Boldyrev2013; Dalena et al. Reference Dalena, Rappazzo, Dmitruk, Greco and Matthaeus2014; Osman et al. Reference Osman, Kiyani, Chapman and Hnat2014a, Reference Osman, Matthaeus, Gosling, Greco, Servidio, Hnat, Chapman and Phanb; Zhdankin, Uzdensky & Boldyrev Reference Zhdankin, Uzdensky and Boldyrev2015a, Reference Zhdankin, Uzdensky and Boldyrevb; Loureiro & Boldyrev Reference Loureiro and Boldyrev2017b; Mallet et al. Reference Mallet, Schekochihin and Chandran2017a, Reference Mallet, Schekochihin and Chandranb). It has also been argued that stochastic heating can also play a role in energizing particles in intermittent turbulence at the location of coherent structures where the fluctuation amplitudes are unusually large (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010; Xia et al. Reference Xia, Perez, Chandran and Quataert2013; Mallet et al. Reference Mallet, Klein, Chand Ran, Grošelj, Hoppock, Bowen, Salem and Bale2019).

Each of the kinetic damping mechanisms in these three categories (with helpful abbreviations) are summarized in table 5. Also included in table 5 are the primary dimensionless parameters upon which each damping mechanism depends, as explored in more detail in the following subsections.

5.2. A competition between turbulent cascade rates and damping rates

One of the most challenging aspects of modelling the damping of plasma turbulence is that the net energy transfer to each species arises from the competition between the turbulent cascade rate and the kinetic damping rates as a function of scale, as depicted by the arrows in figure 1(a). Nonlinear interactions among turbulent fluctuations – denoted Alfvén wave collisions (Howes & Nielson Reference Howes and Nielson2013) in the case of Alfvénic turbulence – mediate the transfer of energy from large to small scales that is dominantly local in scale space (Howes et al. Reference Howes, Tenbarge and Dorland2011a; Told et al. Reference Told, Jenko, TenBarge, Howes and Hammett2015), where energy is transferred locally from a wavenumber $k$ to a wavenumber $2k$, then on to $4k$ and so on. This is indicated by the black curved arrows in the figure. Because the physics of the cascade is dominantly local in nature (Howes et al. Reference Howes, Tenbarge and Dorland2011a; Told et al. Reference Told, Jenko, TenBarge, Howes and Hammett2015), scaling theory implies that the energy transfer rate at a given wavenumber $k$ is determined solely by the conditions at that scale. Note that, when analysed in Fourier space, the energy transfer to small scales in the turbulent cascade is the sum of the nonlinear interactions among many triads of plane-wave modes, with a wide spread of energy transfer rates about zero; the net energy transfer given by this sum at a given scale is much smaller in amplitude than the spread of all of the individual three-wave interaction rates, and the sum typically indicates a forward cascade of energy to small scales for turbulence in three spatial dimensions (Coburn et al. Reference Coburn, Smith, Vasquez, Forman and Stawarz2014, Reference Coburn, Forman, Smith, Vasquez and Stawarz2015). The turbulent cascade rate $\varepsilon$ used here specifically refers to this net energy transfer rate obtained by integrating over all contributing nonlinear interactions involving that scaleFootnote ¹⁰ .

Competing with the turbulent cascade rate at each scale is the damping rate due to collisionless interactions between the electromagnetic fields and the plasma particles, as shown by the blue arrows for the transfer of energy to ions and the magenta arrows for the transfer of energy to electrons in figure 1(a). The net damping rate of the entire turbulent cascade, resulting from the sum of the energization rates for each particle species $s$, can be computed by integration over all scales of the turbulent cascade. In a steady state, at each wavenumber $k_*$, the combination of three terms representing energy transfer into and out of wavenumber $k_*$ must balance to zero, as illustrated in figure 4: (i) the energy transfer from lower wavenumbers $k^-< k_*$ (larger scales), denoted by $\varepsilon (k^-)^{(k^-< k_*)}$; (ii) the energy transfer to higher wavenumbers $k^+>k_*$ (smaller scales), denoted by $\varepsilon (k_*)^{(k_*< k^+)}$; and (iii) the kinetic damping of the fluctuations at wavenumber $k_*$, given by $Q(k_*)$, which is the sum of the local (in scale) ion energization rate $Q_i(k_*)$ (blue) and electron energization rate $Q_e(k_*)$ (magenta). Complicating the determination of the partitioning of damped turbulent energy between ions and electrons $Q_i/Q_e$ is the fact that the kinetic damping rates due to ions and electrons may overlap as a function of scale (Howes et al. Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011b, Reference Howes, Tenbarge and Dorlanda; Told et al. Reference Told, Jenko, TenBarge, Howes and Hammett2015).

Determining this balance of cascade and damping rates as a function of scale, and ultimately using the terms in that balance to determine the partitioning of turbulent energy among species and degrees of freedom – parameterized by $Q_i/Q_e$, $Q_{\perp,i}/Q_{i}$ and $Q_{\perp,e}/Q_{e}$ – has been attempted in empirical turbulent cascade models (Pao Reference Pao1965; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a; Podesta, Borovsky & Gary Reference Podesta, Borovsky and Gary2010; Howes et al. Reference Howes, Tenbarge and Dorland2011a; Zhao, Wu & Lu Reference Zhao, Wu and Lu2013; Schreiner & Saur Reference Schreiner and Saur2017). An example of the modelling of the turbulent energy cascade and its kinetic dissipation is presented in figure 5 from the results of the weakened cascade model by Howes et al. (Reference Howes, Tenbarge and Dorland2011a). In addition to balancing the physics of the local cascade rate $\varepsilon (k_\perp )$ with the kinetic damping rates due to ions $Q_i(k_\perp )$ and electrons $Q_e(k_\perp )$ via the Landau resonance, this cascade model also includes the contributions to the local cascade rate from non-local interactions by large-scale shearing and small-scale diffusion. The model results shown here in figure 5 are from the same calculation as presented in figure 8 of Howes et al. (Reference Howes, Tenbarge and Dorland2011a), where we reproduce (a) the one-dimensional perpendicular magnetic-energy spectrum $E_{B_\perp } (k_\perp )$ and (b) the anisotropic cascade through wavevector space characterized by the characteristic parallel wavenumber as a function of the perpendicular wavenumber, $k_\parallel (k_\perp )$. This model for a plasma with $\beta _i=1$ and isotropic equilibrium Maxwellian velocity distributions with $T_i/T_e=9$ is driven weakly with $\chi (k_{\perp 0})=0.1$ at $k_{\perp 0} \rho _i=k_\parallel \rho _i=10^{-5}$, taking the model Kolmogorov constants to be $C_1=1.4$ and $C_2=1.0$. Note that an important avenue of research is to use kinetic numerical simulations or spacecraft observations to constrain the values of these Kolmogorov constants, as recently done by Shankarappa, Klein & Martinović (Reference Shankarappa, Klein and Martinović2023) using observations in the inner heliosphere by the Parker Solar Probe (Fox et al. Reference Fox, Velli, Bale, Decker, Driesman, Howard, Kasper, Kinnison, Kusterer and Lario2016).

Figure 5. Weakened cascade model results for a turbulent cascade driven weakly at $k_{\perp 0} \rho _i=k_{\parallel 0} \rho _i=10^{-5}$ with nonlinearity parameter $\chi (k_{\perp 0})=0.1$ in a plasma with isotropic equilibrium Maxwellian velocity distributions with $T_i/T_e=9$ and $\beta _i=1$. (a) One-dimensional perpendicular magnetic-energy spectrum $E_{B_\perp }(k_\perp )$, (b) anisotropic cascade through wavevector space given by $k_\parallel (k_\perp )$, (c) energy transfer rates to the ions $Q_i(k_\perp )$ (blue) and electrons $Q_e(k_\perp )$ (magenta) and (d) integrated energy damping rates by ions $\int _0^{k_\perp } d k'_\perp Q_i(k'_\perp )$ (blue) and by electrons $\int _0^{k_\perp } {\rm d} k'_\perp Q_e(k'_\perp )$ (magenta), all plotted as a function of the perpendicular wavenumber $k_\perp \rho _i$. The ranges of $k_\perp \rho _i$ of the full turbulent cascade include weak MHD turbulence (blue), the transition from weakto strong MHD turbulence (green), strong MHD turbulence (red), strong KAW turbulence (magenta) and WDKT (cyan) (Howes et al. Reference Howes, Tenbarge and Dorland2011a).

In figure 5(c) is plotted the rate of energy density transfer to the ions $Q_i(k_\perp )$ (blue) and to the electrons $Q_e(k_\perp )$ (magenta) as a function of $k_\perp \rho _i$, both normalized by the turbulent energy cascade rate at the driving scale $\varepsilon _0$. This shows that the ion kinetic damping occurs at the ion kinetic scales $k_\perp \rho _i \sim 1$, whereas the electron damping rises monotonically from the ion scales, peaking around $k_\perp \rho _i \sim 30$ for this case. To obtain the net partitioning of energy between ions and electrons $Q_i/Q_e$, it is necessary to integrate these energy damping rates over the entire cascade, as shown in figure 5(d): the turbulent energy cascade rate $\varepsilon (k_\perp )$ (black) is ultimately terminated at $k_\perp \rho _i \gg 1$ by the integrated energy damping rates by ions $\int _0^{k_\perp } {\rm d} k'_\perp Q_i(k'_\perp )$ (blue) and by electrons $\int _0^{k_\perp } {\rm d} k'_\perp Q_e(k'_\perp )$ (magenta).

The energy density transfer rates $Q_i(k_\perp )$ and $Q_e(k_\perp )$ that remove energy from the turbulent fluctuations and transfer it to the plasma particles in turbulent cascade models require as input the kinetic damping rates due to ions $\gamma _i(k_\perp )$ and electrons $\gamma _e(k_\perp )$. In the case of the weakened cascade model shown in figure 5 (Howes et al. Reference Howes, Tenbarge and Dorland2011a), the damping rates used are collisionless damping rates via the Landau ($n=0$) resonance for the Alfvén wave mode (and its extension as a KAW at scales $k_\perp \rho _i \gtrsim 1$) from the Vlasov–Maxwell linear dispersion relation (Stix Reference Stix1992). Previous studies have called into question whether resonant collisionless wave–particle interactions, such as Landau damping, can effectively damp turbulent fluctuations in the presence of the strong nonlinear interactions that mediate the turbulent cascade (Plunk Reference Plunk2013; Schekochihin et al. Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016). Using the field-particle correlation technique (Klein & Howes Reference Klein and Howes2016; Howes, Klein & Li Reference Howes, Klein and Li2017; Klein, Howes & TenBarge Reference Klein, Howes and TenBarge2017), several recent analyses of kinetic simulations (Klein et al. Reference Klein, Howes and TenBarge2017; Howes et al. Reference Howes, McCubbin and Klein2018; Klein et al. Reference Klein, Howes, TenBarge and Valentini2020) and spacecraft observations (Chen et al. Reference Chen, Klein and Howes2019; Afshari et al. Reference Afshari, Howes, Kletzing, Hartley and Boardsen2021) of plasma turbulence find clear velocity-space signatures of Landau damping. These findings indicate indeed that resonant wave–particle interactions do play a role, and possibly a dominant role (Afshari et al. Reference Afshari, Howes, Kletzing, Hartley and Boardsen2021), in the damping of space and astrophysical plasma turbulence, settling the question of whether Landau damping can effectively damp strong plasma turbulenceFootnote ¹¹ .

In quantifying how a particular damping mechanism competes with the nonlinear turbulent cascade, the key dimensionless quantity to evaluate is the ratio of linear damping rate due to nonlinear cascade rate $\gamma /\omega _{nl}$ (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a, Reference Howes, Tenbarge and Dorland2011a). Adopting the conjecture of critical balance between linear and nonlinear time scales in strong plasma turbulence $\omega \sim \omega _{nl}$ (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a, Reference Howes, Tenbarge and Dorland2011a; Mallet et al. Reference Mallet, Schekochihin and Chandran2015), the key dimensionless measure of the importance of any particular mechanism is the ratio $\gamma /\omega$, which can be estimated using the linear theory for each of the different damping mechanisms. Nonlinear kinetic simulations and spacecraft observations of weakly collisional plasma turbulence present opportunities to test whether the linear collisionless damping rates provide a quantitatively accurate, lowest-order estimate of the turbulent damping rates as a function of wavenumber. It is clear from the numerical and observational evidence, however, that Landau damping (Klein et al. Reference Klein, Howes and TenBarge2017; Howes et al. Reference Howes, McCubbin and Klein2018; Chen et al. Reference Chen, Klein and Howes2019; Klein et al. Reference Klein, Howes, TenBarge and Valentini2020; Afshari et al. Reference Afshari, Howes, Kletzing, Hartley and Boardsen2021; Zhou et al. Reference Zhou, Liu and Loureiro2023) and cyclotron damping (Afshari et al. Reference Afshari, Howes, Shuster, Klein, McGinnis, Martinovic, Boardsen, Brown and Huang2024) can play a significant role in the damping of astrophysical plasma turbulence.

5.3. The relevance of linear physics to strong plasma turbulence

The concept of critical balance in plasma turbulence (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a, Reference Howes, Tenbarge and Dorland2011a; Mallet et al. Reference Mallet, Schekochihin and Chandran2015) – which suggests that the nonlinear time scales of the turbulent cascade remain in approximate balance with the time scales of the linear wave physics, $\omega _{nl} \sim \omega$ – implies that the physics of the linear plasma response remains valid even in strong plasma turbulence. But, perhaps in analogy with the case of incompressible hydrodynamic turbulence in which no linear response exists (Howes Reference Howes2015a), the validity of linear physics has been called into question in the presence of strong plasma turbulence (Matthaeus et al. Reference Matthaeus, Oughton, Osman, Servidio, Wan, Gary, Shay, Valentini, Roytershteyn and Karimabadi2014). Although there exist numerous counterexamples in which linear physics properties have been shown to be relevant to strong plasma turbulence (Maron & Goldreich Reference Maron and Goldreich2001; Cho & Lazarian Reference Cho and Lazarian2003; Alexandrova et al. Reference Alexandrova, Carbone, Veltri and Sorriso-Valvo2008; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a, Reference Howes, Dorland, Cowley, Hammett, Quataert, Schekochihin and Tatsunob; Svidzinski et al. Reference Svidzinski, Li, Rose, Albright and Bowers2009; Sahraoui et al. Reference Sahraoui, Goldstein, Belmont, Canu and Rezeau2010; Howes et al. Reference Howes, TenBarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011b; Hunana et al. Reference Hunana, Laveder, Passot, Sulem and Borgogno2011; TenBarge et al. Reference TenBarge, Podesta, Klein and Howes2012; Chen et al. Reference Chen, Howes, Bonnell, Mozer, Klein and Bale2013), the question has persisted. A more recent careful study of kinetic numerical simulations and spacecraft observations has shown clearly that both small-amplitude and large-amplitude fluctuations, including variations interpreted to be coherent structures, preserve the linear properties of Alfvénic fluctuations (Grošelj et al. Reference Grošelj, Chen, Mallet, Samtaney, Schneider and Jenko2019), hopefully providing sufficient evidence to establish definitively the relevance of linear physics to strong plasma turbulence.

The key to predicting qualitatively which kinetic damping mechanisms are likely to dominate for a given instance of turbulence is to determine how the various collisionless damping mechanisms in § 5.1 depend on the fundamental dimensionless plasma and turbulence parameters. For a fully ionized proton and electron plasma with anisotropic (bi-Maxwellian) equilibrium velocity distribution functions (Stix Reference Stix1992; Quataert Reference Quataert1998; Quataert & Gruzinov Reference Quataert and Gruzinov1999; Swanson Reference Swanson2003; Klein & Howes Reference Klein and Howes2015), the complex eigenfrequencies determined by the linear Vlasov–Maxwell dispersion relation depend on the dimensionless parameters that characterize the nature of the turbulent fluctuations through the typical wavevector components in a plane-wave decomposition of the turbulent fluctuations, $k_\perp \rho _i$ and $k_\parallel \rho _i$, and on the plasma parameters through the parallel ion plasma beta $\beta _{\parallel,i}$, the parallel ion-to-electron temperature ratio $\tau _\parallel =T_{\parallel, i}/T_{\parallel, e}$, the temperature anisotropies of each species $A_i=T_{\perp,i}/T_{\parallel,i}$ and $A_e=T_{\perp,e}/T_{\parallel,e}$, and the ratio $v_{t\parallel,i}/c$. Thus, these dependencies can be expressed through the general relation

(5.2)\begin{equation} \omega =\omega_{{\rm VM}}^{{\rm aniso}}(k_\perp \rho_i, k_\parallel \rho_i, \beta_{{\parallel},i},\tau_\parallel,A_i,A_e,v_{t\parallel,i}/c). \end{equation}

In the limit of isotropic equilibrium velocity distributions for the ions and electrons – leading to the simplifications $A_i=1$, $A_e=1$ and $\tau _\parallel = \tau \equiv T_{i}/T_{e}$ – this set of dependencies is reduced to

(5.3)\begin{equation} \omega =\omega_{{\rm VM}}(k_\perp \rho_i, k_\parallel \rho_i, \beta_i,\tau,v_{ti}/c). \end{equation}

For both the isotropic and anisotropic equilibrium temperatures, in the non-relativistic limit $v_{ti}/c \ll 1$ typical of many space and astrophysical plasmas of interest, the linear Vlasov–Maxwell dispersion relation is practically independent of $v_{ti}/c$, so this parameter may effectively be dropped for the case of non-relativistic plasma turbulence.

In the limit of a sufficiently large inertial range of turbulence, the MHD turbulence scaling parameters given by (3.1)–(3.4) predict that turbulent fluctuations at the small-scale end of the inertial range – where most of the kinetic damping mechanisms are believed to become significant – will be anisotropic with $k_\parallel \ll k_\perp$ and small amplitude relative to the equilibrium magnetic field $|\delta \boldsymbol {B}| \ll B_0$ (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006, Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). In this limit, the turbulent frequencies of the incompressible Alfvénic fluctuations directly observed to dominate turbulence in space plasmas (Belcher & Davis Reference Belcher and Davis1971; Tu & Marsch Reference Tu and Marsch1995; Alexandrova et al. Reference Alexandrova, Carbone, Veltri and Sorriso-Valvo2008; Howes et al. Reference Howes, Bale, Klein, Chen, Salem and TenBarge2012a; Bruno & Carbone Reference Bruno and Carbone2013) have frequencies much smaller than the ion cyclotron frequency, $\omega /\varOmega _i = k_\parallel d_i \ll 1$, where $d_i=v_A/\varOmega _i=\rho _i/\sqrt {\beta _i}$ is the ion inertial length. In this limit, the turbulent fluctuations at the small-scale end of the inertial range satisfy the conditions for the applicability of gyrokinetic theory (Rutherford & Frieman Reference Rutherford and Frieman1968; Taylor & Hastie Reference Taylor and Hastie1968; Catto Reference Catto1978; Antonsen & Lane Reference Antonsen and Lane1980; Catto, Tang & Baldwin Reference Catto, Tang and Baldwin1981; Frieman & Chen Reference Frieman and Chen1982; Dubin et al. Reference Dubin, Krommes, Oberman and Lee1983; Hahm, Lee & Brizard Reference Hahm, Lee and Brizard1988; Brizard Reference Brizard1992; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Brizard & Hahm Reference Brizard and Hahm2007; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). Gyrokinetic theory can be rigorously derived from plasma kinetic theory under the limits of low frequency ($\omega /\varOmega _i \ll 1$) and anisotropic ($k_\parallel \ll k_\perp$) fluctuations in the non-relativistic limit ($v_{ti}/c \ll 1$). In the gyrokinetic approximation, the physics of the the cyclotron resonances and of the fast magnetosonic fluctuations – and their kinetic extension to whistler fluctuations (Howes, Klein & TenBarge Reference Howes, Klein and TenBarge2014) – are ordered out of the system. Gyrokinetics retains the physics of the $\boldsymbol {E} \times \boldsymbol {B}$ nonlinearity that mediates the turbulent cascade to small scales, full finite-Larmor-radius effects, and the collisionless damping associated with the $n=0$ Landau resonance. In this reduced limit, the complex eigenfrequencies of the linear gyrokinetic dispersion relation for isotropic equilibrium velocity distributions are linearly dependent on the parallel wavenumber, and have the significantly reduced dimensionless parameter dependency given by (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006)

(5.4)\begin{equation} \bar{\omega} \equiv \frac{\omega}{k_\parallel v_A} =\bar{\omega}_{GK}(k_\perp \rho_i, \beta_i,\tau). \end{equation}

Although the ion-to-electron mass ratio $\mu =m_i/m_e$ is not a parameter that varies in the physical world, many numerical studies of plasma turbulence choose to employ a reduced mass ratio for computational efficiency. Therefore, we choose to include the $\mu$ dependence of different proposed turbulent damping mechanisms in our calculations below to determine how the use of a reduced mass ratio may impact the partitioning of turbulent energy removed by different damping mechanisms. Note also that it is important to employ a sufficiently large mass ratio to achieve a separation between the ion and electron scales, with a previous study suggesting $m_i/m_e\ge 32$ is necessary to model with qualitative fidelity the distinct ion and electron responses to the electromagnetic fluctuations due to this scale separation (Howes et al. Reference Howes, McCubbin and Klein2018).

5.4. Landau damping and transit-time damping

To explore the two Landau-resonant collisionless damping mechanisms – Landau damping and transit-time damping – we take $n=0$ in the resonance condition in (5.1), obtaining the condition that particles with parallel velocities equal to the parallel phase velocity of the waves, $v_\parallel = \omega /k_\parallel$, may resonantly exchange energy with the electromagnetic waves.

For simplicity in this analysis, we will assume the isotropic temperature case for a proton–electron plasma with linear wave properties given by the isotropic linear Vlasov–Maxwell dispersion relation (5.3); extending these calculations to allow for anisotropic temperatures can be done simply following the same logical procedure and using the linear Vlasov–Maxwell dispersion relation for anisotropic temperatures given by (5.2). Furthermore, justified by observations in solar wind turbulence that incompressible Alfvénic fluctuations energetically dominate the turbulence (Belcher & Davis Reference Belcher and Davis1971; Tu & Marsch Reference Tu and Marsch1995; Alexandrova et al. Reference Alexandrova, Carbone, Veltri and Sorriso-Valvo2008; Bruno & Carbone Reference Bruno and Carbone2013), we will focus on the linear collisionless damping rates of the Alfvénic fluctuations, including both Alfvén waves in the MHD limit at $k_\perp \rho _i \ll 1$, and KAWs at $k_\perp \rho _i \gtrsim 1$. Consideration of the damping of compressible turbulent fluctuations, specifically fast and slow magnetosonic fluctuations at $k_\perp \rho _i \ll 1$ (Howes et al. Reference Howes, Bale, Klein, Chen, Salem and TenBarge2012a), is easily performed using the same procedure with the linear dispersion relations for these alternative wave modesFootnote ¹² .

Normalizing the resonant parallel velocity by the thermal velocity for both species in turn and converting the result into our dimensionless parameters for the isotropic temperature case, we obtain the key relations

(5.5a,b)\begin{equation} \frac{v_\parallel}{v_{ti}} = \frac{\omega}{k_\parallel v_{ti}} = \bar{\omega}\beta_i^{{-}1/2},\quad \frac{v_\parallel}{v_{te}} = \frac{\omega}{k_\parallel v_{te}} = \bar{\omega} \left(\frac{\tau}{\beta_i \mu}\right)^{1/2}, \end{equation}

where we define the dimensionless wave frequency normalized by the Alfvén wave frequency in the MHD limit

(5.6)\begin{equation} \bar{\omega} \equiv \frac{\omega}{k_\parallel v_A}. \end{equation}

For Alfvén waves, this normalized frequency $\bar {\omega }$ is well approximated by (Howes et al. Reference Howes, Klein and TenBarge2014)

(5.7)\begin{equation} \bar{\omega}(k_\perp\rho_i,\beta_i, \tau) = \sqrt{1 + \frac{(k_\perp\rho_i)^2}{\beta_i + 2/(1+1/\tau)}}, \end{equation}

an expression that is valid in the limits that the ion plasma beta is not very low $\beta _i \gg \mu ^{-1}$ and that the wave frequency remains smaller than the ion cyclotron frequency $\omega /\varOmega _i \ll 1$, equivalent to the condition $k_\parallel d_i \ll 1$, or $k_\parallel \rho _i \ll \beta _i^{1/2}$.

Whether damping by the Landau-resonant mechanisms is weak or strong primarily depends on where the parallel wave phase velocity falls within the equilibrium particle velocity distribution, as quantified by the expressions in (5.5a,b). The application of the weak-growth-rate approximation to the kinetic theory for Landau damping (Krall & Trivelpiece Reference Krall and Trivelpiece1973) indicates that the linear damping rate is proportional to the slope of the distribution function at the resonant velocity, as depicted in figure 6. For phase velocities that fall deep in the core of the velocity distribution $\omega /k_\parallel \ll v_{ts}$ (blue dotted) where the slope of $f_s(v_\parallel )$ is relatively flat, the damping is weak because nearly equal numbers of particles with $v_\parallel < \omega /k_\parallel$ are accelerated by the wave as particles with $v_\parallel > \omega /k_\parallel$ are decelerated by the wave, leading to little net energy exchange. Similarly, for phase velocities far out in the tail of the velocity distribution $\omega /k_\parallel \gg v_{ts}$ (black dashed), there are few particles to resonate with the waves and the slope of $f_s(v_\parallel )$ is also quite flat, so the damping is weak. Only if the parallel phase velocity falls in the steep region of the velocity distribution at $\omega /k_\parallel \sim v_{ts}$ (red solid) is there a significant net energy transfer to the particles: more particles with $v_\parallel < \omega /k_\parallel$ are accelerated by the wave than those particles with $v_\parallel > \omega /k_\parallel$ are decelerated by the wave, so that particles gain net energy at the expense of the wave. Note that this qualitative picture of collisionless damping by the Landau resonance being controlled by where the parallel wave phase velocity falls within the distribution function is independent of the form of the equilibrium velocity distribution; thus, one can apply the same reasoning to estimate the strength of the Landau-resonant damping for other non-Maxwellian forms of the equilibrium velocity distribution, such as the kappa distributions often applied to model in situ measurements in space plasmas (Vasyliunas Reference Vasyliunas1968; Abraham-Shrauner & Feldman Reference Abraham-Shrauner and Feldman1977; Gosling et al. Reference Gosling, Asbridge, Bame, Feldman, Zwickl, Paschmann, Sckopke and Hynds1981; Lui & Krimigis Reference Lui and Krimigis1981; Armstrong et al. Reference Armstrong, Paonessa, Bell and Krimigis II1983; Summers & Thorne Reference Summers and Thorne1991; Thorne & Summers Reference Thorne and Summers1991; Summers, Xue & Thorne Reference Summers, Xue and Thorne1994; Livadiotis & McComas Reference Livadiotis and McComas2013; Livadiotis, Desai & Wilson Reference Livadiotis, Desai and Wilson III2018).

Figure 6. Diagram demonstrating how Landau-resonant damping qualitatively depends on where the parallel phase velocity $\omega /k_\parallel$ of the wave falls within the velocity distribution of particle species $s$, with weak damping for $\omega /k_\parallel \ll v_{ts}$ or $\omega /k_\parallel \gg v_{ts}$, and strong damping for $\omega /k_\parallel \sim v_{ts}$.

The ultimate goal is to characterize the turbulent plasma heating – specifically $Q_i/Q_e$, $Q_{\perp,i}/Q_{i}$ and $Q_{\perp,e}/Q_{e}$ – due to Landau damping and transit-time damping as a function of the fundamental dimensionless parameters of turbulence. To accomplish this goal, it is necessary to integrate the damping rates over the full range of the turbulent cascade, as discussed in § 5.2 and illustrated in figures 4 and 5. This integration over all scales effectively eliminates the dependence of the damping mechanisms on the components of the wavevector $k_\perp \rho _i$ and $k_\parallel \rho _i$, leaving only the desired dependence on the plasma and turbulence parameters enumerated in tables 3 and 4.

With the scale dependence on $k_\perp \rho _i$ eliminated, the ion Landau damping (iLD) and ion transit-time damping (iTTD) rates depend primarily on $\beta _i$, and the electron Landau damping (eLD) and electron transit-time damping (eTTD) rates depend on $\beta _i$, $\tau$ and $\mu$, as shown by (5.5a,b). Note that the dependence of eLD and eTTD on $\tau$ and $\mu$ primarily arises from those two parameters controlling the scale separation of the ion and electron Larmor-radius scales, $\rho _i/\rho _e = (\tau \mu )^{1/2}$.

It is worthwhile pointing out that the physics of collisionless damping via the Landau resonance is fully captured in the gyrokinetic approximation, with the linear wave frequencies and damping rates both linearly proportional to $k_\parallel$, as shown in (5.4). Thus, the ratio of the damping-to-cascade rateFootnote ¹³ , $\gamma /\omega _{nl} \simeq \gamma /\omega$, is effectively independent of $k_\parallel$, leaving the expected parameter dependence on the plasma parameters $\beta _i$ and $\tau$, along with the mass ratio $\mu$, once the cascade is integrated over the perpendicular scales of the cascade.

In addition to the fact that the dependence of the linear gyrokinetic dispersion relation on $(k_\perp \rho _i, \beta _i,\tau )$ in (5.4) is consistent with the results of our analysis of the parameter dependence of the Landau-resonant damping mechanisms, the linear gyrokinetic dispersion relation also determines the amplitudes and phases of the key electromagnetic fields $E_\parallel$ and $\delta B_\parallel$ that mediate the particle energization by Landau damping and transit-time damping. As done in the H10 turbulent heating model (see § 4), the linear dispersion relation can simply be used to calculate directly the linear collisionless damping rates as a function of $(k_\perp \rho _i, \beta _i,\tau )$. It is possible, of course, that these linear damping rates may be modified under strongly turbulent plasma conditions, but modest quantitative changes to the effective damping rates are typically absorbed into the adjustable constant typically used in cascade models to specify the balance between the nonlinear cascade rate and collisionless damping rate (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a, Reference Howes, Tenbarge and Dorland2011a). However, field-particle correlation analyses of numerical simulations (Klein et al. Reference Klein, Howes and TenBarge2017, Reference Klein, Howes, TenBarge and Valentini2020; Horvath, Howes & McCubbin Reference Horvath, Howes and McCubbin2020; Conley, Howes & McCubbin Reference Conley, Howes and McCubbin2023; Huang, Howes & McCubbin Reference Huang, Howes and McCubbin2024) and spacecraft observations (Chen et al. Reference Chen, Klein and Howes2019; Afshari et al. Reference Afshari, Howes, Kletzing, Hartley and Boardsen2021) of plasma turbulence have clearly shown that Landau damping and transit-time damping do indeed play a role in the damping of strong plasma turbulence.

One final step is to estimate how the ion-to-electron heating ratio $Q_i/Q_e$ depends on the remaining key parameters $\beta _i$, $\tau$ and $\mu$. As illustrated in figure 5(c) and 7, the ion damping rate $\gamma _i$ peaks at scales $k_\perp \rho _i \sim 1$, and numerical solutions of the Vlasov–Maxwell dispersion relation (not presented here) show that this Landau-resonant damping onto the ions is essentially independent of the ion-to-electron temperature ratio. The electron damping rate $\gamma _e$, on the other hand, increases monotonically as the perpendicular wavevector increases to the electron Larmor-radius scale $k_\perp \rho _e \rightarrow 1$, eventually yielding a normalized collisionless damping rate that is believed to be sufficiently strong to terminate the turbulent cascade with $\gamma _e/\omega \rightarrow 1$, as shown in figure 7. Thus, it is expected that any energy in the turbulent cascade that passes beyond the ion kinetic scales at $k_\perp \rho _i \sim 1$ will ultimately be transferred to the electrons. Since the ion damping rate depends dominantly on $\beta _i$ with very weak dependence on $\tau$, yielding $Q_i(\beta _i)$ to lowest order, this means that the ratio $Q_i/Q_e$ is effectively determined by the ion plasma beta alone. This simplification arises because the total turbulent damping rate must eventually balance the total cascade rate, $Q =Q_i+ Q_e \sim \varepsilon _0$, so it can be shown that $Q_e/Q_i= \varepsilon _0/Q_i(\beta _i) - 1$. This dominant dependence of the Landau-resonant collisionless damping mechanisms on $\beta _i$ is consistent with existing turbulent heating models (Howes Reference Howes2010) and kinetic numerical simulations (Kawazura et al. Reference Kawazura, Barnes and Schekochihin2019) of Alfvénic turbulence.

Figure 7. Plot of the normalized damping rate $\gamma /\omega$ (black dashed) for the Alfvén wave and KAW over the perpendicular wavenumber range $10^{-2} \le k_\perp \rho _i \le 10^2$ for a plasma with $\beta _i=3$, $T_i/T_e=1$, and $m_i/m_e=1836$, where the complex eigenfrequency is given by $\omega _c=\omega -i\gamma$, so that $\gamma >0$ corresponds to wave damping and $\gamma <0$ to wave growth. The contributions from the ion damping rate $\gamma _i$ (thin blue) and electron damping rate $\gamma _e$ (thin magenta) are separately plotted, with the decomposition by species failing at strong damping rates with $\gamma /\omega \gtrsim 0.5$, occurring for $k_\perp \rho _i \gtrsim 40$. Ion damping due to Landau damping $\gamma _{i,{\rm LD}}>0$ (red dashed) and transit-time damping $\gamma _{i,{\rm TTD}}>0$ (green long dashed) are separated, showing a region of net energy transfer from ions to the wave by the Landau resonance with the ions yielding $\gamma _{i,{\rm LD}}<0$ (red dotted) at $k_\perp \rho _i \le 2$.

One final unexpected physical behaviour of Landau damping and transit-time damping is also illustrated in figure 7. For isotropic Maxwellian velocity distributions of the ions and electrons, the imaginary component of the complex frequency, $\omega _c \equiv \omega -{\rm i} \gamma$, determined by the linear dispersion relation is always negative (solid black curve), corresponding to a damping of the wave with damping rate $\gamma >0$. However, for the particular plasma parameter choices in figure 7 of $\beta _i=3$, $T_i/T_e=1$, and $m_i/m_e=1836$, the individual energy transfer rates to the ions through Landau dampingFootnote ¹⁴ mediated by $E_\parallel$ are actually negative over the range $k_\perp \rho _i \lesssim 2$ (red dotted), indicating a net transfer of energy from the ions to the wave fields when averaged over the full period (or, alternatively, over $2{\rm \pi}$ in phase) of the Alfvén wave. But the contribution from transit-time damping mediated by $E_\perp$ (Howes et al. Reference Howes, Huang and Felix2024) is positive and much larger (green long dashed), leading to a net damping of the wave ($\gamma _i >0$, thin blue) by the sum of the Landau-resonant ($n=0$) interactions of ions with the Alfvén wave. Whether the wave-period-averaged energy transfer to ions is positive or negative depends on the relative phases of $E_\parallel$ and $E_\perp$ to the self-consistently determined components of the ion current density $\boldsymbol {j}_i$. For the case of $\beta _i=3$ pictured in figure 7, at scales $k_\perp \rho _i \lesssim 2$ (red dotted), this relative phase $\delta \phi$ between $E_\parallel$ and $j_{\parallel,i}$ must fall in the range ${\rm \pi} /2 \le \delta \phi \le 3 {\rm \pi}/2$ (corresponding to wave growth at the expense of ion energy), while the relative phase between $E_\perp$ and $j_{\perp,i}$ must fall in the range $-{\rm \pi} /2 \le \delta \phi \le {\rm \pi}/2$ (corresponding to wave damping and ion energy gain). Because each individual mechanism can possibly lead to negative energy transfer while the sum always leads to damping for an isotropic Maxwellian distribution, it suggests that perhaps a separation of the particle energization rates by Landau damping and transit-time damping is not particularly important in unravelling the turbulent heating of the plasma (especially since both of these mechanisms ultimately lead to energization of the parallel degree of freedom of the particles), so that only the summed effect of the Landau-resonant collisionless damping rates is important.

Summarizing our results for Landau-resonant collisionless damping of the turbulent cascade, the resulting ion-to-electron heating is predicted to be dominantly a function of the ion plasma beta, $Q_i/Q_e(\beta _i)$, with a weak dependence on the ion-to-electron temperature ratio $\tau$ and mass ratio $\mu$. Both Landau damping, mediated by the parallel component of the electric field $E_\parallel$, and transit-time damping, mediated by the perpendicular component of the electric field $E_\perp$, lead to energization of the degrees of freedom parallel to the magnetic field, so $Q_{\perp,i}/Q_{i}=0$ and $Q_{\perp,e}/Q_{e}=0$.

5.5. Cyclotron damping

For the fundamental ($n=1$) cyclotron resonance, the resonance condition in (5.1) simplifies to $(\omega -k_\parallel v_\parallel )/\varOmega _s =1$. The numerator involves the fluctuation frequency modified by the Doppler shift of the particle's motion parallel to the magnetic field. When that Doppler-shifted frequency matches the ion cyclotron frequency, the perpendicular component of the electric field $E_\perp$ can resonantly interact with the particles to energize the particles in the directions perpendicular to the magnetic field, yielding $Q_{\perp,s}/Q_s=1$. To estimate the conditions required for cyclotron resonant damping, it is often sufficient to require the condition $\omega /\varOmega _s \rightarrow 1$.

In typical space and astrophysical plasmas, the turbulent cascade can reach the ion cyclotron frequency $\omega /\varOmega _i \rightarrow 1$ under the conditions derived below, but since the electron cyclotron frequency is higher by the mass ratio, $\varOmega _e/\varOmega _i = m_i/m_e$, there is typically little energy remaining in the cascade at the very high frequencies needed to transfer a significant amount of energy to the electrons through their cyclotron resonance. Thus, here, we consider only the physics of the ion cyclotron resonance as a possible physical mechanism for the dissipation of plasma turbulence.

Evaluating the conditions needed to obtain $\omega /\varOmega _i \rightarrow 1$ in terms of the dimensionless parameters of turbulence in the isotropic temperature case, we find

(5.8)\begin{equation} \frac{\omega}{\varOmega_i} = \bar{\omega} k_\parallel d_i = \bar{\omega} \frac{k_\parallel \rho_i}{\beta_i^{1/2}}, \end{equation}

where $\bar {\omega }$ is defined by (5.6). Note that the middle expression indicates that the most natural unit to normalize the length scale parallel to the magnetic field is the ion inertial length $d_i=c/\omega _{pi}=v_A/\varOmega _i$, but to maintain consistency with the length normalizations in tables 3 and 4, we substitute $d_i= \rho _i/ \beta _i^{1/2}$. Utilizing the general scaling of the parallel wavenumber in terms of the isotropic driving wavenumber

(5.9)\begin{equation} k_\parallel \rho_i = (k_0 \rho_i)^{(1+\alpha)/(3+\alpha)}(k_\perp \rho_i)^{2/(3+\alpha)}, \end{equation}

to determine the frequency achieved when the cascade reaches $k_\perp \rho _i=1$, we obtain the condition

(5.10)\begin{equation} \frac{\omega}{\varOmega_i} = \frac{\bar{\omega}} {\beta_i^{1/2}} (k_0 \rho_i)^{(1+\alpha)/(3+\alpha)}. \end{equation}

This indicates clearly that ion cyclotron damping is more likely to occur for small ion plasma beta $\beta _i \ll 1$ or a smaller turbulent inertial range (meaning a value of $k_0 \rho _i$ that is not too small). Note that, for anisotropic driving with $k_{\parallel 0}/k_{\perp 0} < 1$, such as may occur in the solar corona at low values of $\beta _i$, the effective isotropic driving wavenumber $k_0 \rho _i$, given by (3.11), will actually be smaller than the perpendicular driving scale $k_{\perp 0} \rho _i$, leading to yet smaller values of $\omega /\varOmega _i$. For the GS95 scaling with $\alpha =0$, we obtain a condition for the onset of ion cyclotron damping of $k_0 \rho _i =\beta _i^{3/2}/\bar {\omega }^3$; for the B06 scaling with $\alpha =1$, we obtain $k_0 \rho _i=\beta _i/\bar {\omega }^2$. For simplicity, note that at the small-scale end of the inertial range with $k_\perp \rho _i \lesssim 1$, the normalized wave frequency $\bar {\omega }$ for Alfvén waves can be approximated in the limit $k_\parallel d_i \rightarrow 1$ by the expression (Di Mare & Howes Reference Di Mare and Howes2024)

(5.11)\begin{equation} \bar{\omega} \simeq \frac{1}{\sqrt{1+(k_\parallel d_i)^2}}. \end{equation}

In summary, ion cyclotron damping is predicted to play a role in the damping of plasma turbulence for conditions with a smaller inertial range and low plasma beta values $\beta _i \ll 1$. The nature of the resulting plasma heating is strictly ion energization $Q_i/Q=1$ in the perpendicular degrees of freedomFootnote ¹⁵ $Q_{\perp,i}/Q_i=1$.

Note that recent theoretical and numerical studies of turbulence with a sufficiently large amount of imbalance $Z_0^+/Z_0^- \gg 1$, such as may be encountered in the inner heliosphere, has uncovered a new phenomenon denoted the helicity barrier, which can lead to the generation of sufficiently high-frequency fluctuations, specifically ion cyclotron waves, that may cause energization of the ions by the cyclotron resonance (Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021; Squire et al. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2022, Reference Squire, Meyrand and Kunz2023).

5.6. Stochastic heating

The stochastic heating of ions has long been proposed as a means to energize ions in turbulent plasmas (McChesney et al. Reference McChesney, Stern and Bellan1987; Chen et al. Reference Chen, Lin and White2001; Johnson & Cheng Reference Johnson and Cheng2001; White et al. Reference White, Chen and Lin2002; Voitenko & Goossens Reference Voitenko and Goossens2004; Bourouaine et al. Reference Bourouaine, Marsch and Vocks2008), with more recent work providing significant theoretical and numerical evidence for ion stochastic heating as a significant mechanism for turbulent damping in space and astrophysical plasmas (Chandran Reference Chandran2010; Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010, Reference Chandran, Dennis, Quataert and Bale2011; Bourouaine & Chandran Reference Bourouaine and Chandran2013; Chandran et al. Reference Chandran, Verscharen, Quataert, Kasper, Isenberg and Bourouaine2013; Klein & Chandran Reference Klein and Chandran2016; Vech, Klein & Kasper Reference Vech, Klein and Kasper2017; Arzamasskiy et al. Reference Arzamasskiy, Kunz, Chand Ran and Quataert2019; Martinović et al. Reference Martinović, Klein, Kasper, Case, Korreck, Larson, Livi, Stevens, Whittlesey and Chandran2020; Cerri et al. Reference Cerri, Arzamasskiy and Kunz2021).

Chandran et al. (Reference Chandran, Li, Rogers, Quataert and Germaschewski2010) used test particle simulations of randomly phased Alfvén and kinetic Alfvén waves to develop an empirical formula for the rate of ion stochastic heating by turbulent fluctuations with frequencies much lower than the ion cyclotron frequency $\omega /\varOmega _i \ll 1$. Following this, Chandran et al. (Reference Chandran, Dennis, Quataert and Bale2011) developed a refined implementation to predict the partitioning of turbulent energy among parallel ion heating, perpendicular ion heating, and electron heating (the C11 model in § 4), with a formula for the damping rate by ion stochastic heating given by

(5.12)\begin{equation} \gamma_{{\rm SH}} = c_1 \left(\frac{\delta v_\perp^{(\rho_i)}}{v_{ti}}\right) \varOmega_i \exp({-c_2 v_{ti}/\delta v_\perp^{(\rho_i)}}), \end{equation}

where $\delta v_\perp ^{(\rho _i)}$ is the root-mean-square amplitude of the perpendicular velocity fluctuations at the perpendicular scale of the ion Larmor radius, $k_\perp \rho _i=1$. Here, $c_1$ and $c_2$ are fitting constants with suggested values $c_1=0.18$ and $0.15 \lesssim c_2 \lesssim 0.34$ (Chandran Reference Chandran2010; Chandran et al. Reference Chandran, Dennis, Quataert and Bale2011; Bourouaine & Chandran Reference Bourouaine and Chandran2013; Xia et al. Reference Xia, Perez, Chandran and Quataert2013). The key dimensionless parameter controlling the onset of stochastic heating in this formulation is the ratio of the perpendicular velocity fluctuation to the thermal velocity, $\epsilon \equiv \delta v_\perp ^{(\rho _i)}/v_{ti}$; these studies found a critical threshold $\epsilon _{{\rm crit}} \simeq 0.19$ that determines the minimum fluctuation amplitude required for the onset of stochastic heating (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010).

Using the relation for Alfvén waves to connect the perpendicular fluid velocity fluctuations to the perpendicular magnetic field fluctuations $\delta v_\perp /v_A =\pm \bar {\omega } \delta B_\perp /B_0$ (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a) and the MHD scaling of the turbulent perpendicular magnetic field fluctuation amplitude given by (3.4) with the B06 scaling with $\alpha =1$, we obtain the expression

(5.13)\begin{equation} \frac{\delta v_\perp}{v_{ti}} = \bar{\omega} \beta_i^{{-}1/2} \left(\frac{k_0\rho_i}{k_{{\perp}} \rho_i}\right)^{1/4}. \end{equation}

At $k_\perp \rho _i \sim 1$, we obtain the normalized frequency $\bar {\omega } \sim 1$ and find the key relation for the amplitude in terms of our fundamental turbulence parameters

(5.14)\begin{equation} \epsilon \equiv \frac{\delta v_\perp^{(\rho_i)}}{v_{ti}} = \beta_i^{{-}1/2} (k_0\rho_i)^{1/4}. \end{equation}

To estimate the importance of damping by ion stochastic heating relative to the turbulent cascade rate, we can calculate $\gamma _{{\rm SH}}/\omega _{nl} \sim \gamma _{{\rm SH}}/\omega$ to obtain the final result

(5.15)\begin{equation} \frac{\gamma_{{\rm SH}}}{\omega} = c_1(k_0\rho_i)^{{-}1/4} \exp({-c_2 \beta_i^{1/2} (k_0\rho_i )^{{-}1/4}}), \end{equation}

where we have used the substitution $\varOmega _i/\omega = \beta _i^{1/2}/k_\parallel \rho _i$.

In summary, the proposed turbulent damping mechanism of ion stochastic heating depends primarily on the fundamental parameters $\beta _i$ and $k_0\rho _i$, where lower ion plasma beta and larger isotropic driving wavenumber (which corresponds to a smaller inertial range) both conspire to yield larger turbulent amplitude at the ion Larmor-radius scale $k_\perp \rho _i \sim 1$, enhancing stochastic heating of the ions. The exponential effectively provides a threshold turbulent amplitude at the ion scales below which stochastic ion heating is ineffective. Note that the exponential dependence on the constant $c_2$ makes its accurate determination critical to assess the importance of ion stochastic heating in any particular case; its value is likely to depend on the detailed nature of the turbulent fluctuations at the perpendicular scale of the ion Larmor radius (Chandran et al. Reference Chandran, Dennis, Quataert and Bale2011). This heating mechanism is predicted to lead to strictly ion energization $Q_i/Q=1$ in the perpendicular degrees of freedom $Q_{\perp,i}/Q_i=1$.

Note that, because the amplitude of the turbulent fluctuations is typically found to decrease monotonically as the turbulent cascade progresses to smaller scales, it is unlikely that electrons will experience significant stochastic heating, so we do not explore electron stochastic heating here. Further research will be necessary to evaluate thoroughly the potential for electron stochastic heating to damp turbulence at the perpendicular scale of the electron Larmor radius, $k_\perp \rho _e \sim 1$. Furthermore, a consideration of the intermittency – which can enhance the amplitude of individual fluctuations – of the turbulent fluctuations at the ion Larmor radius scale suggests that the rate of ion stochastic heating may be dramatically increased in the presence of significant intermittency (Mallet et al. Reference Mallet, Klein, Chand Ran, Grošelj, Hoppock, Bowen, Salem and Bale2019).

5.7. Magnetic pumping

Magnetic pumping (Spitzer & Witten Reference Spitzer and Witten1953; Berger et al. Reference Berger, Newcomb, Dawson, Frieman, Kulsrud and Lenard1958; Lichko et al. Reference Lichko, Egedal, Daughton and Kasper2017; Lichko & Egedal Reference Lichko and Egedal2020; Montag & Howes Reference Montag and Howes2022) is a particle energization mechanism in which variations of magnetic field magnitude at low frequency $\omega$ lead to an oscillating transfer of energy from parallel to perpendicular degrees of freedom and back, obeying the double-adiabatic (Chew–Goldberger–Low, CGL) equations of state (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956). In the absence of any collisions, the net change of particle energy over a full oscillation is zero. But, in the presence small but finite collisionalityFootnote ¹⁶ , there can be a non-zero transfer of energy from the oscillating electromagnetic fields to the particles over the course of a full oscillation.

As a proposed mechanism for the damping of turbulence, magnetic pumping necessarily requires turbulent fluctuations that involve changes of the magnetic field magnitude from the equilibrium field strength $B_0$, or $\delta B \equiv |\boldsymbol {B}|-B_0 \ne 0$. For Alfvénic turbulence, which is dominantly incompressibleFootnote ¹⁷ in the MHD regime at $k_\perp \rho _i \ll 1$, significant fluctuations of the magnetic field magnitude arise only within the KAW regime at $k_\perp \rho _i \gtrsim 1$ under conditions of ion plasma beta $\beta _i > 1$. Alternatively, if the turbulence has a significant component of compressible magnetosonic fluctuations with $E_{{\rm comp}}/E_{{\rm inc}}>0$, because the turbulent amplitudes typically decrease monotonically with decreasing scale lengths, the largest-amplitude compressible fluctuations occur at the driving scale and may lead to magnetic pumping, directly energizing particles from these large-scale fluctuations.

We determine first the parameters controlling magnetic pumping in Alfvénic turbulence, following a similar procedure subsequently to analyse the magnetic pumping arising from compressible turbulence with $E_{{\rm comp}}/E_{{\rm inc}}>0$. The foundation of magnetic pumping is the double-adiabatic evolution of the weakly collisional plasma which conserves the particle magnetic moment $\mu _m \equiv m v_\perp ^2/2B$ and the action integral of the parallel bounce motion $J_\parallel \equiv \oint v_\parallel dl$ (Lichko et al. Reference Lichko, Egedal, Daughton and Kasper2017; Montag Reference Montag2018). Conservation of these adiabatic invariants requires a scale separation of the perpendicular wavenumber relative to the particle species Larmor radius $k_\perp \rho _s \ll 1$ and of the wave parallel phase velocity relative to the particle thermal velocity $\omega /k_\parallel v_{ts} \ll 1$.

Kinetic Alfvén waves give rise to non-negligible magnetic field magnitude variations over the perpendicular wavenumber range $1 \lesssim k_\perp \rho _i \lesssim (\tau \mu )^{1/2}$ for values of $\beta _i \gtrsim 1$. This range of $\delta B$ fluctuations is incompatible with the necessary condition $k_\perp \rho _i \ll 1$ for the ions, so we expect that magnetic pumping by KAWs cannot energize the ions. For the electrons, on the other hand, the requisite condition $k_\perp \rho _e \ll 1$ can be rewritten as $k_\perp \rho _i \ll (\tau \mu )^{1/2}$, so we expect that KAWs may be able to damp turbulent fluctuations and energize the electrons under ion plasma beta conditions with $\beta _i \gtrsim 1$. The phase velocity condition for electrons can be expressed as $\bar {\omega } \ll (\beta _i \mu /\tau )^{1/2}$. Under the ion plasma beta condition $\beta _i \gtrsim 1$, the equation for the dimensionless frequency of kinetic Alfvén waves given by (5.7) can be approximated by $\bar {\omega } \simeq k_\perp \rho _i/\beta _i^{1/2}$, so this requirement can be converted to a condition on the perpendicular wavenumber given by $k_\perp \rho _i \ll \beta _i (\mu /\tau )^{1/2}$, which can certainly be satisfied for the perpendicular wavenumber range of KAWs with significant magnetic field magnitude fluctuations. Therefore, the necessary conditions for the magnetic pumping of electrons to play a role in the dissipation of Alfvénic plasma turbulence with $\beta _i \gtrsim 1$ appear to be satisfied.

Previous studies have shown that the ratio of the damping rate due to magnetic pumping $\gamma _{{\rm MP}}$ to the pump wave frequency $\omega$ is given by (Lichko et al. Reference Lichko, Egedal, Daughton and Kasper2017; Montag & Howes Reference Montag and Howes2022)

(5.16)\begin{equation} \frac{\gamma_{{\rm MP}}}{\omega} = \frac{\chi_{c,s}}{1+\chi_{c,s}^2}\left(\frac{\delta B}{B_0}\right)^2, \end{equation}

where $\chi _{c,s} \equiv 6 \nu _s/\omega$ is a dimensionless measure of the collision frequency $\nu _s$ for species $s$ and $\delta B$ is the amplitude of the magnetic field magnitude fluctuations.

To determine the dependence of the magnetic pumping of electrons on the plasma and turbulence parameters in tables 3 and 4, we assume isotropic driving $k_{\parallel 0}/k_{\perp 0}=1$ for simplicity, write the electron–electron collision frequency as $\nu _e = v_{te}/\lambda _{{\rm mfp},e}$, and use the general parallel wavenumber scaling for arbitrary $\alpha$ given by (3.1) to obtain

(5.17)\begin{equation} \chi_{c,e} = \frac{6}{\bar{\omega} (k_{{\parallel} 0} \lambda_{{\rm mfp},e})} \left(\frac{\beta_i \mu}{\tau}\right)^{1/2} \left(\frac{k_0 \rho_i}{k_\perp \rho_i}\right)^{2/(3+\alpha)}. \end{equation}

For a sufficiently large turbulent inertial range with $k_0\rho _i \ll 1$, KAWs over the perpendicular wavenumber range $1 \lesssim k_\perp \rho _i \lesssim (\tau \mu )^{1/2}$ will have relatively small turbulent amplitudes $\delta B_\perp /B_0 \ll 1$, so the magnetic field magnitude fluctuations can be expressed to lowest order as variations of the parallel component of the perturbed magnetic field $\delta B \simeq \delta B_\parallel$ (Huang et al. Reference Huang, Howes and McCubbin2024). The turbulent magnetic field magnitude variations may then be expressed as

(5.18)\begin{equation} \frac{\delta B}{B_0} = \left(\frac{\delta B_\parallel}{\delta B_\perp}\right) \left( \frac{k_0 \rho_i}{k_\perp \rho_i} \right)^{1/(3+\alpha)}, \end{equation}

where the ratio $\delta B_\parallel /\delta B_\perp$ for KAWs over the relevant range from the linear dispersion relation is strictly a function of $\beta _i$, or $\delta B_\parallel /\delta B_\perp =f(\beta _i)$. Therefore, electron magnetic pumping by Alfvénic turbulent fluctuations at scales $k_\perp \rho _i \gtrsim 1$ depends on the parameters $\beta _i$, $\tau$, $\mu$, $k_0 \rho _i$ and $k_{\parallel 0}\lambda _{{\rm mfp},e}$.

Next we consider that case that the turbulence is driven compressibly with $E_{{\rm comp}}/E_{{\rm inc}}>0$. In this situation, the magnetic field magnitude fluctuations at the driving scale, which typically have the largest amplitude of $\delta B/B_0$, may contribute to magnetic pumping of both the ions and electrons. In this case of turbulent damping at the driving scale, it is necessary to use the three driving parameters $(k_{\perp 0} \rho _i,k_{\parallel 0}/k_{\perp 0},\chi _0)$ in table 3, rather than the simplified parameterization through the isotropic driving wavenumber $k_0 \rho _i$, to determine the parameter dependence of the turbulent damping via magnetic pumping.

First, the normalized fast and slow magnetosonic wave phase velocities can be expressed by the linear dispersion relation (Klein et al. Reference Klein, Howes, TenBarge, Bale, Chen and Salem2012)

(5.19)\begin{equation} \tilde{\omega}^4 - \tilde{\omega}^2(1 + \beta) + \beta \cos^2 \varTheta=0, \end{equation}

which has four solutions given by

(5.20)\begin{equation} \tilde{\omega}^2 = \tfrac{1}{2}[1+ \beta \pm \sqrt{(1 + \beta)^2 - 4 \beta \cos^2 \varTheta}], \end{equation}

where the plus (minus) sign corresponds to the two fast (slow) magnetosonic wave modes. Here, we define the normalized wave frequency $\tilde {\omega } \equiv \omega /k v_A$, the total (MHD) plasma beta $\beta = \beta _i (1 + 1/\tau )$, and the wavevector components $(k_\parallel,k_\perp )$ are alternatively parameterized by $(k,\varTheta )$, where $k_\parallel = k \cos \varTheta$, $k_\perp = k \sin \varTheta$, $k=(k_\parallel ^2+k_\perp ^2)^{1/2}$ and $\varTheta =\tan ^{-1}(k_\perp /k_\parallel )$. The dimensionless form of the solutions in (5.20) makes clear that the fast and slow magnetosonic wave phase velocities depend only on two dimensionless parameters, $\tilde {\omega }( \beta,\varTheta )$.

We may determine the parameter dependence of the normalized ion collision frequency by using the relation for the ion–ion collision rate in terms of the electron-electron collision rate $\nu _i = \mu ^{-1/2} \tau ^{-3/2} \nu _e$ and writing the pumping wave frequency at the driving scale as $\omega = \tilde {\omega }_0 k_0 v_A$ to obtain

(5.21)\begin{equation} \chi_{c,i} = \frac{6}{ \tilde{\omega}_0 (k_{{\parallel} 0} \lambda_{{\rm mfp},e}) \tau^2} \beta_i^{1/2} \cos \varTheta_0. \end{equation}

Note that the ion collisional mean free path can be expressed in terms of the electron collisional mean free path by the relation $\lambda _{{\rm mfp},i}= \lambda _{{\rm mfp},e} \tau ^2$, so the denominator is actually independent of the ion-to-electron temperature ratio $\tau$, depending rather only on the normalized ion mean free path, $k_{\parallel 0} \lambda _{{\rm mfp},i}$.

To evaluate the turbulent magnetic field magnitude variations for ion magnetic pumping by compressible fluctuations at the driving scale, we first note that the phenomenon of dynamic alignment, parameterized by the alignment angle $\theta _0$, will not play a role: dynamic alignment develops through the nonlinear interactions as the turbulence cascades to smaller scales, but does not come into play at the driving scale, so we set $\theta _0 =1$, thus reducing the nonlinearity parameter at the driving scale to $\chi _0 = k_{\perp 0} \delta B_{\perp 0}/(k_{\parallel 0} B_0)$. If the magnetic perturbations at the driving scale have $|\delta \boldsymbol {B}|/B_0 \sim 1$, the linear theory for waves breaks down, and the nonlinearity must be taken into account in the lowest-order description of the wave dynamics; in this case, one may directly take the measured magnetic field magnitude fluctuations $\delta B /B_0$ as the parameter to estimate the damping rate by magnetic pumping in (5.16). However, if the magnetic perturbations at the driving scale have $\epsilon \equiv |\delta \boldsymbol {B}|/B_0 <1$, the binomial expansion may be used to obtain the expression

(5.22)\begin{equation} \delta B /B_0 \simeq \delta B_\parallel{/}B_0 + (\delta B_\parallel{/}B_0)^2/2 + (\delta B_\perp{/}B_0)^2/2, \end{equation}

where the first term is of order $\epsilon$ and the next two terms are of order $\epsilon ^2$. Thus, to lowest order, we may estimate $\delta B \simeq \delta B_\parallel$, where the error in this estimation scales as $\epsilon ^2/2$, yielding a fractional error of 25 % for $|\delta \boldsymbol {B}|/B_0 =1/2$ and of 6 % for $|\delta \boldsymbol {B}|/B_0 =1/4$.

Using this approximation, we can derive an expression for the normalized fluctuations of the magnetic field magnitude by $\delta B /B_0 \simeq \delta B_\parallel /B_0 = (\delta B_\perp /B_0) (\delta B_\parallel /\delta B_\perp )$. If we estimate $\delta B_\parallel /\delta B_\perp \sim (E_{{\rm comp}}/E_{{\rm inc}})^{1/2}$, we can then obtain the expression for the normalized amplitude of the magnetic field magnitude fluctuations

(5.23)\begin{equation} \frac{\delta B}{B_0} = \frac{k_{{\parallel} 0}}{k_{{\perp} 0}} \chi_0 \left(\frac{E_{{\rm comp}}}{E_{{\rm inc}}}\right)^{1/2}. \end{equation}

Combining the dependencies of $\chi _{c,i}$ in (5.21) and of $\delta B /B_0$ in (5.23), ion magnetic pumping by compressible turbulent fluctuations at the driving scale depends on the parameters $\beta _i$, $k_{\parallel 0}\lambda _{{\rm mfp},i}= (k_{\parallel 0}\lambda _{{\rm mfp},e}) \tau ^2$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$ and $E_{{\rm comp}}/E_{{\rm inc}}$.

For magnetic pumping of electrons by compressible turbulent fluctuations at the driving scale, we can follow an analogous procedure to obtain the normalized electron collision frequency

(5.24)\begin{equation} \chi_{c,e} = \frac{6}{\tilde{\omega}_0 (k_{{\parallel} 0} \lambda_{{\rm mfp},e})} \left(\frac{\beta_i \mu}{\tau}\right)^{1/2} \cos \varTheta_0, \end{equation}

and the same expression in (5.23) for the normalized fluctuations of the magnetic field magnitude by $\delta B /B_0$, which is independent of the species properties. Thus, electron magnetic pumping by compressible turbulent fluctuations at the driving scale depends on the parameters $\beta _i$, $\tau$, $\mu$, $k_{\parallel 0}\lambda _{{\rm mfp},e}$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$ and $E_{{\rm comp}}/E_{{\rm inc}}$.

In summary, magnetic pumping by Alfvénic turbulent fluctuations can damp turbulence and energize electrons through turbulent fluctuations at perpendicular scales $1 \lesssim k_\perp \rho _i \lesssim (\tau \mu )^{1/2}$ for plasma conditions with $\beta _i \gtrsim 1$, but ions are not expected to be energized by Alfvénic fluctuations. The resulting magnetic pumping of electrons by Alfvénic turbulent fluctuations depends on the parameters $\beta _i$, $\tau$, $\mu$, $k_0 \rho _i$, and $k_{\parallel 0}\lambda _{{\rm mfp},e}$. If the turbulence is driven compressibly with $E_{{\rm comp}}/E_{{\rm inc}}>0$, we expect that the turbulent fluctuations at the driving scales, which have the largest amplitude of $\delta B/B_0$, may lead to damping of the turbulence through the magnetic pumping of both the ions and electrons. The resulting magnetic pumping of ions by compressible turbulence at the driving scales depends on the parameters $\beta _i$, $k_{\parallel 0}\lambda _{{\rm mfp},i}= (k_{\parallel 0}\lambda _{{\rm mfp},e}) \tau ^2$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$, and $E_{{\rm comp}}/E_{{\rm inc}}$. Analogously, the magnetic pumping of electrons by compressible turbulence depends on the parameters $\beta _i$, $\tau$, $\mu$, $k_{\parallel 0}\lambda _{{\rm mfp},e}$, $k_{\parallel 0}/k_{\perp 0}$, $\chi _0$, and $E_{{\rm comp}}/E_{{\rm inc}}$. These parameter dependencies for magnetic pumping of ions and electrons are summarized in table 5. Note that magnetic pumping energizes particles in energy over all pitch angles (Montag & Howes Reference Montag and Howes2022), so that both the parallel and perpendicular degrees of freedom gain energy, yielding for either ion or electron species $Q_{\perp,s}/Q_{\parallel,s} \sim 1$.

5.8. Kinetic viscous heating

As discussed above for the case of magnetic pumping, low-frequency turbulent fluctuations that generate variations in the magnetic field magnitude will lead to a double-adiabatic (CGL) (Chew et al. Reference Chew, Goldberger and Low1956) evolution of weakly collisional plasmas that preserves to lowest order the particle magnetic moment $\mu _m$ and the action integral of the parallel bounce motion $J_\parallel$ (Lichko et al. Reference Lichko, Egedal, Daughton and Kasper2017; Montag Reference Montag2018; Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023). This dynamics leads to the development of anisotropies in the velocity distribution of the plasma particles, yielding temperature and pressure anisotropies with $p_\parallel \ne p_\perp$. These deviations from local thermodynamic equilibrium can trigger rapidly growing kinetic instabilities that can disrupt the linear wave dynamics (Squire et al. Reference Squire, Quataert and Schekochihin2016, Reference Squire, Kunz, Quataert and Schekochihin2017a, Reference Squire, Schekochihin and Quataertb) and mediate the non-local transfer of turbulent energy directly from the large driving scales to kinetic length scales (Kunz et al. Reference Kunz, Abel, Klein and Schekochihin2018; Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023), as illustrated by the red arrow in figure 1(b). The effect of these kinetic instabilities on the turbulent dynamics has recently been modelled as an effective ‘viscosity’ (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023).

Here we use the term ‘kinetic viscous heating’ to describe the transfer of energy by this effective viscosity mediated by temperature anisotropy instabilities (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023) to avoid confusion with the usual viscous heating arising from microscopic Coulomb collisions. We note that, although this effect can be grossly modelled as an effective viscosity that leads to diffusion of the bulk fluid velocity field, its properties are entirely distinct from the standard viscosity in fluid dynamics mediated by microscopic collisions between neutral particles (Chapman & Cowling Reference Chapman and Cowling1970) or by microscopic Coulomb collisions between charged particles in a plasma under strongly collisional conditions (Spitzer Reference Spitzer1962; Grad Reference Grad1963; Braginskii Reference Braginskii1965). Standard viscosity is a linear mechanism that may be modelled mathematically by a Laplacian differential operator in the fluid momentum equation, and thus its effect increases monotonically with $k^2$ in Fourier space; this kinetic viscosity is a nonlinear mechanism operating only under weakly collisional conditions that strongly depends on the value of the ion plasma beta $\beta _i$, and its effect diminishes at higher wavenumber $k$ with the monotonically decreasing amplitudes of the turbulent fluctuations. Determining how this kinetic viscosity may tap some fraction of the energy from the turbulent fluctuations at the large, driving scales and non-locally transfer that energy directly to fluctuations at the small, ion kinetic length scales remains an open line of investigation.

To properly describe the kinetic physics of turbulence in temperature anisotropic plasmas, it is necessary to employ the full set of plasma and turbulence parameters in table 3. The relevant kinetic instabilities driven by anisotropies in the ion velocity distribution are the four ion temperature anisotropy instabilities: (i) the parallel (or whistler) firehose instability (Kennel Reference Kennel1966; Gary et al. Reference Gary, Montgomery, Feldman and Forslund1976) and (ii) the Alfvén (or oblique) firehose instability (Hellinger & Matsumoto Reference Hellinger and Matsumoto2000), which are both relevant to plasmas with an ion temperature anisotropy sufficiently less than unity $A_i=T_{\perp,i}/T_{\parallel,i}<1$ and with a parallel ion plasma beta $\beta _{\parallel,i} >1$; and (iii) the mirror instability (Vedenov & Sagdeev Reference Vedenov and Sagdeev1958; Tajiri Reference Tajiri1967; Southwood & Kivelson Reference Southwood and Kivelson1993) and (iv) the ion cyclotron instability (Gary et al. Reference Gary, Montgomery, Feldman and Forslund1976), which are both relevant to plasmas with an ion temperature anisotropy sufficiently greater than unity $A_i=T_{\perp,i}/T_{\parallel,i}>1$ for any value of the parallel ion plasma beta $\beta _{\parallel,i}$. When the plasma exceeds a threshold value of the ion temperature anisotropy (Boris & Manheimer Reference Boris and Manheimer1977; Matteini et al. Reference Matteini, Landi, Hellinger and Velli2006; Hellinger & Trávníček Reference Hellinger and Trávníček2008; Klein & Howes Reference Klein and Howes2015), these instabilities can tap the free energy associated with the anisotropic ion temperature, driving electromagnetic fluctuations and ultimately reducing the temperature anisotropy, thereby moving the plasma back toward a state of marginal stability.

Hellinger et al. (Reference Hellinger, Trávníček, Kasper and Lazarus2006) has compiled values for the marginal stability boundaries for these four ion temperature anisotropy instabilities over the $(\beta _{\parallel,i}, T_{\perp,i}/T_{\parallel,i})$ plane. Using a linear Vlasov-Maxwell dispersion relation solver with bi-Maxwellian equilibrium ion velocity distributions in a fully ionized, hydrogenic (proton and electron) plasma, the marginal stability boundary is determined by calculating the complex eigenfrequency $\omega (\boldsymbol {k})$ for a fixed $\beta _{\parallel,i}$ over all possible wavevectors $\boldsymbol {k}$. The ion temperature anisotropy $T_{\perp,i}/T_{\parallel,i}$ is then varied until the most unstable wavevector has a growth rate of $|\gamma |/\varOmega _i=10^{-3}$, thus establishing the instability criterion. For the four ion temperature anisotropy instabilities, the corresponding instability criteria on the $(\beta _{\parallel,i}, T_{\perp,i}/T_{\parallel,i})$ plane are generally well fit by an expression of the form

(5.25)\begin{equation} \frac{T_{{\perp},i}}{T_{{\parallel},i}}=1+\frac{a}{(\beta_{{\parallel},i}-\beta_0)^b}, \end{equation}

where $a, b,$ and $\beta _0$ are unique values for each of the four instabilities, computed by Hellinger et al. (Reference Hellinger, Trávníček, Kasper and Lazarus2006) and found here in table 6. These instability thresholds are plotted in figure 8, and observational studies of measured intervals in the solar wind generally show that the measurements are more or less constrained by these threshold limits on the $(\beta _{\parallel,i}, T_{\perp,i}/T_{\parallel,i})$ plane (Hellinger et al. Reference Hellinger, Trávníček, Kasper and Lazarus2006; Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009; Maruca, Kasper & Bale Reference Maruca, Kasper and Bale2011).

Figure 8. Marginal linear stability thresholds for the four ion temperature anisotropy instabilities using (5.25) and the parameters in table 6.

Table 6. Instability threshold parameters for a maximum growth rate $|\gamma |/\varOmega _i=10^{-3}$ from Hellinger et al. (Reference Hellinger, Trávníček, Kasper and Lazarus2006).

Recent theoretical and numerical studies have highlighted the fact that, at sufficiently high parallel ion plasma beta $\beta _{\parallel,i} \gg 1$, the marginal stability boundary for these kinetic ion temperature anisotropy instabilities asymptotes towards isotropy, as shown in figure 8, so the double-adiabatic evolution of turbulent fluctuations can easily exceed these boundaries and trigger these rapidly growing instabilities. Large-scale Alfvén waves in the MHD limit at $k_\perp \rho _i \ll 1$ yield only perpendicular magnetic field perturbations $\delta B_\perp \ne 0$ with $\delta B_\parallel =0$, so the variation in the magnetic field magnitude $\delta B$ is negligible in the limit of small-amplitude waves $\delta B_\perp /B_0 \ll 1$. However, at sufficiently large amplitudes $\delta B_\perp /B_0 \rightarrow 1$, these fluctuations can lead to significant variation in the magnetic field magnitude, as made clear by (5.22). Thus, at sufficiently high values of parallel ion plasma beta $\beta _{\parallel,i} \gg 1$, the physics of linearly polarizedFootnote ¹⁸ Alfvén waves can be ‘interrupted’ by these kinetic ion temperature anisotropy instabilities (Squire et al. Reference Squire, Quataert and Schekochihin2016, Reference Squire, Kunz, Quataert and Schekochihin2017a, Reference Squire, Schekochihin and Quataertb).

Further examination of the effect of kinetic instabilities on the turbulent cascade under conditions of high parallel ion plasma beta, $\beta _{\parallel,i} \gg 1$, has found that these instabilities can non-locally transfer energy directly to kinetic length scales (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023), as illustrated by the red arrow in figure 1(b). This new channel of energy transfer can compete with the traditional local transfer of energy to smaller scales of the turbulent cascade that is mediated by nonlinear interactions, as illustrated by the black arrows in figure 1(a). The effect of these kinetic instabilities on the large-scale turbulent dynamics has been quantified as an effective viscosity by Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023). Here, we denote this kinetic-instability-mediated viscosity as the kinetic viscosityFootnote ¹⁹ to distinguish it from the usual definition of viscosity that arises through Coulomb collisions among the charged plasma particles.

These ion kinetic instabilities act to limit the ion pressure anisotropy in the plasma to values corresponding to marginal stability, given by $\Delta p_i/p_i \lesssim 1/\beta _{\parallel,i}$, where $\Delta p_i \equiv p_{\perp i} - p_{\parallel i}$ and $p_i \equiv (2 p_{\perp i} + p_{\parallel i})/3$ (Kunz, Schekochihin & Stone Reference Kunz, Schekochihin and Stone2014; Melville, Schekochihin & Kunz Reference Melville, Schekochihin and Kunz2016; Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023). Using this marginal stability scaling with $\beta _{\parallel,i}$, one can estimate an effective collision frequency (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023)

(5.26)\begin{equation} \nu_{{\rm eff}} \sim \beta_{{\parallel},i} \omega (\delta B_\perp{/}B_0)^2, \end{equation}

where $\omega$ is the linear wave frequency at the scale where the dynamics generates the largest pressure anisotropy, thus leading to the maximum instability growth rate. Since the turbulent fluctuation amplitudes generally decrease monotonically with scale, the kinetic temperature anisotropy instabilities will be driven most strongly by turbulent fluctuations at the driving scale. As with the case for magnetic pumping by compressible fluctuations that is dominated by the driving-scale dynamics, here, we must use the three driving parameters $(k_{\perp 0} \rho _i,k_{\parallel 0}/k_{\perp 0},\chi _0)$ in table 3, rather than $k_0 \rho _i$, to estimate the turbulent damping rate by the kinetic viscosity.

Using the effective collision frequency $\nu _{{\rm eff}}$ above as an estimate of the damping rate by this kinetic viscous heating $\gamma _{{\rm VH}}$, we can calculate the ratio of the damping rate arising from kinetic viscous heating to the wave frequency at the driving scale, obtaining

(5.27)\begin{equation} \frac{\gamma_{{\rm VH}}}{\omega}\sim \beta_{{\parallel},i} \left(\chi_0 \frac{k_{{\parallel} 0}}{k_{{\perp} 0}}\right)^2, \end{equation}

where the definition of the nonlinearity parameter $\chi _0 = (k_{\perp 0}/k_{\parallel 0})(\delta B_{\perp 0}/B_0) \theta _0$ with $\theta _0=1$ at the driving scale is used to replace the amplitude of the driving-scale turbulent fluctuations $\delta B_{\perp 0}/B_0$.

In summary, the kinetic viscous heating due to ion temperature anisotropy instabilities arising from the non-thermal velocity distributions driven by the large-scale turbulent motions occurs in weakly collisional turbulent plasmas with parallel ion plasma beta $\beta _{\parallel,i} \gg 1$. The kinetic viscous damping rate depends on the parallel ion plasma beta $\beta _{\parallel,i}$, the nonlinearity parameter $\chi _0$, and the anisotropy of the turbulent driving $k_{\perp 0}/k_{\parallel 0}$. It is important to keep in mind that the scaling used to derive the form of the effective collision frequency $\nu _{{\rm eff}}$ depends on the scaling of the marginal stability boundaries that are generally calculated in quiescent plasmas in the absence of turbulence. The scaling of the marginal stability boundaries – such as those quantified by (5.25) and the coefficients in table 6 – may differ in the presence of strong plasma turbulence, although evidence from spacecraft observations (Hellinger et al. Reference Hellinger, Trávníček, Kasper and Lazarus2006; Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009; Maruca et al. Reference Maruca, Kasper and Bale2011) and from hybrid kinetic ion and fluid electron simulations of temperature anisotropic plasma turbulence (Kunz et al. Reference Kunz, Schekochihin and Stone2014; Bott et al. Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021; Arzamasskiy et al. Reference Arzamasskiy, Kunz, Squire, Quataert and Schekochihin2023) appear to suggest that any modifications of these boundaries due to the presence of turbulence are relatively small.

5.9. Collisionless magnetic reconnection

Numerical simulations of plasma turbulence demonstrate the ubiquitous development of current sheets at small scales (Matthaeus & Montgomery Reference Matthaeus and Montgomery1980; Meneguzzi et al. Reference Meneguzzi, Frisch and Pouquet1981; Biskamp & Welter Reference Biskamp and Welter1989; Spangler Reference Spangler1998, Reference Spangler1999; Biskamp & Müller Reference Biskamp and Müller2000; Maron & Goldreich Reference Maron and Goldreich2001; Merrifield et al. Reference Merrifield, Müller, Chapman and Dendy2005; Greco et al. Reference Greco, Chuychai, Matthaeus, Servidio and Dmitruk2008), and it has been found that the dissipation of the turbulence is largely concentrated in the vicinity of these current sheets (Uritsky et al. Reference Uritsky, Pouquet, Rosenberg, Mininni and Donovan2010; Wan et al. Reference Wan, Matthaeus, Karimabadi, Roytershteyn, Shay, Wu, Daughton, Loring and Chapman2012; Karimabadi et al. Reference Karimabadi, Roytershteyn, Wan, Matthaeus, Daughton, Wu, Shay, Loring, Borovsky and Leonardis2013; TenBarge & Howes Reference TenBarge and Howes2013; Wu et al. Reference Wu, Perri, Osman, Wan, Matthaeus, Shay, Goldstein, Karimabadi and Chapman2013; Zhdankin et al. Reference Zhdankin, Uzdensky, Perez and Boldyrev2013), giving rise to the proposal that magnetic reconnection may play a role in the damping of plasma turbulence. Motivated by these findings, statistical analyses of observations of turbulence in the solar wind have sought evidence for such spatially localized heating (Borovsky & Denton Reference Borovsky and Denton2011; Osman et al. Reference Osman, Matthaeus, Greco and Servidio2011, Reference Osman, Matthaeus, Hnat and Chapman2012a, Reference Osman, Matthaeus, Wan and Rappazzob; Perri et al. Reference Perri, Goldstein, Dorelli and Sahraoui2012; Wang et al. Reference Wang, Tu, He, Marsch and Wang2013; Wu et al. Reference Wu, Perri, Osman, Wan, Matthaeus, Shay, Goldstein, Karimabadi and Chapman2013; Osman et al. Reference Osman, Matthaeus, Gosling, Greco, Servidio, Hnat, Chapman and Phan2014b). Yet the physical mechanisms that lead to the development of these current sheets in plasma turbulence is yet to be fully understood.

Under the B06 theory for the scaling of MHD turbulence, the phenomenon of dynamic alignment (Boldyrev Reference Boldyrev2006) is predicted to lead to the development of current sheets with a width $w\sim 1/k_i$ to thickness $a\sim 1/k_\perp$ ratio in the plane perpendicular to the magnetic field given by $w/a = k_\perp /k_i = [(k_\perp \rho _i)/(k_0\rho _i)]^{1/4}$, according to the scaling in (3.2) for $\alpha =1$. The geometry of the current sheets expected to develop in plasma turbulence is illustrated in figure 9, where for sufficiently small perpendicular scales relative to the driving scale, $k_\perp \rho _i \gg k_0 \rho _i$, the current sheets have a parallel extent along the mean magnetic field $\boldsymbol {B}_0$, given by $l\sim 1/k_\parallel$, and a scale separation of thickness to width to length obeying $a \ll w \ll l$; the corresponding three components of the wavevector have the complementary scaling $k_\parallel \ll k_i \ll k_\perp$. In figure 9, the wavelength of a mode that is unstable to the collisionless tearing instability is illustrated by $\lambda \sim 1/k$ (red). Note that the minimum value of this tearing unstable wavenumber $k$ is equal to the intermediate wavenumber $k_i$ corresponding to the width of the current sheet, $w \sim 1/k_i$. The maximum value of $k$ assumed in calculations of the linear collisionless tearing instability (Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a; Mallet et al. Reference Mallet, Schekochihin and Chandran2017b) must be much smaller than the perpendicular wavenumber $k_\perp$ characterizing the current sheet thickness ($a \sim 1/k_\perp$), thereby yielding $k/k_\perp = k a \ll 1$. Thus, the range of the unstable wavenumbers $k$ is given by

(5.28)\begin{equation} k_i \rho_i = (k_0 \rho_i)^{1/4}(k_\perp \rho_i)^{3/4} \lesssim k \rho_i \ll k_\perp \rho_i \lesssim 1, \end{equation}

where the final constraint $k_\perp \rho _i \lesssim 1$ is needed for the tearing instability to arise for a current sheet thickness within the MHD regime, otherwise different scalings for electron-only reconnection (Phan et al. Reference Phan, Eastwood, Shay, Drake, Sonnerup, Fujimoto, Cassak, Øieroset, Burch and Torbert2018) in the regime $k_\perp \rho _i \gtrsim 1$ would need to be incorporated.

Figure 9. Diagram of the current sheet geometry generated self-consistently in plasma turbulence in the MHD regime at $k_\perp \rho _i \lesssim 1$, with the scalings of $(k_\perp,k_i,k_\parallel )$ given by the B06 scaling. The current sheet has thickness $a\sim 1/k_\perp$ and width $w\sim 1/k_i$ in the perpendicular plane and length $l \sim 1/k_\parallel$ along the equilibrium magnetic field $\boldsymbol {B}_0$. The wavelength of the mode unstable to the collisionless tearing instability is given by $\lambda \sim 1/k$ (red), and must fall within the bounds determined by the current sheet scaling, yielding an ordering $k_i \lesssim k \ll k_\perp$.

Note that, in turbulence obeying the B06 scaling, the self-consistent generation of current sheets arises only at sufficiently small perpendicular scales relative to the driving scale, $k_\perp \rho _i \gg k_0 \rho _i$. Therefore, the ratio of the perpendicular magnetic field perturbation to the mean magnetic field will be small, $\delta B_\perp /B_0 = (k_0 \rho _i/k_\perp \rho _i)^{1/4} \ll 1$, so the reconnection that arises in plasma turbulence necessarily falls in the limit of strong-guide-field magnetic reconnection, as illustrated in figure 9 (green). An implication of this fact is that, if numerical simulations are employed to assess the importance of magnetic reconnection in turbulence, it is critical that those simulations be modelled in three spatial dimensions. Previous studies of two-dimensional vs. three-dimensional plasma turbulence have shown that three-dimensional simulations exhibit a significantly faster turbulent cascade rate (Li et al. Reference Li, Howes, Klein and TenBarge2016), so enabling full three-dimensional evolution is likely to be essential to determine the competition between the turbulent cascade and magnetic reconnection.

The tearing instability growth rate typically increases with an increasing perpendicular aspect ratio $w/a$ ratio of the current sheets. If this ratio reaches a sufficiently large value, it was recently recognized that the growth rate of the tearing instability could exceed the turbulent cascade rate at that scale, $\gamma /\omega _{nl}\gtrsim 1$, triggering a reconnection flow that would disrupt the turbulent dynamics, deviating from the scaling relations summarized in § 3. Such changes in the scaling of the turbulent fluctuations may potentially lead to differences in which physical mechanisms dominate the damping of the turbulent cascade at small scales (Boldyrev & Loureiro Reference Boldyrev and Loureiro2017; Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a, Reference Loureiro and Boldyrevb; Mallet et al. Reference Mallet, Schekochihin and Chandran2017a, Reference Mallet, Schekochihin and Chandranb; Walker et al. Reference Walker, Boldyrev and Loureiro2018).

It is worth noting here that an alternative explanation for the development of current sheets in plasma turbulence due to the nonlinear interactions between counterpropagating Alfvén waves – a phenomenon often denoted Alfvén wave collisions (Howes & Nielson Reference Howes and Nielson2013) – was proposed by Howes (Reference Howes2016), but no scaling for the resulting aspect ratio of the current sheets was provided. If the scaling of this alternative mechanism differs from that of the B06 theory, than the scaling predictions for the onset of reconnection would also change.

Here, we adopt the B06 scaling for the development of current sheets, with the goal to determine the dependence of magnetic reconnection as a turbulent damping mechanism on the fundamental plasma and turbulence parameters for the isotropic temperature case in table 4. This calculation follows previous investigations that assess the conditions under which collisionless magnetic reconnection may disrupt the turbulent cascade and thereby impact how the turbulent energy is channelled into particle energy (Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a; Mallet et al. Reference Mallet, Schekochihin and Chandran2017b).

The growth rate of the collisionless tearing instability, which is the instability that initiates the process of magnetic reconnection under the weakly collisional conditions typical of space and astrophysical plasmas, has been derived in the limit of very low electron plasma beta $\beta _e \sim \mu ^{-1}$ (Zocco & Schekochihin Reference Zocco and Schekochihin2011), where we convert $\beta _e = \beta _i/\tau$. The physics of the collisionless tearing instability depends on three length scales: (i) the thickness of the current sheet $a \sim 1/k_\perp$, assumed here to occur on MHD scales with $k_\perp \rho _i \ll 1$; (ii) the ion sound Larmor radius $\rho _s \equiv \rho _i/(2 \tau )^{1/2}$ where two fluid effects begin to lead to differences between the ion and electron dynamics; and (iii) the inner boundary layer scale $\delta _{in}$ where the magnetic flux can be unfrozen from the electron fluid flow. We will restrict our analysis to $\beta _e \gg \mu ^{-1}$, such that the flux unfreezing arises from electron inertia (Zocco & Schekochihin Reference Zocco and Schekochihin2011; Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a; Mallet et al. Reference Mallet, Schekochihin and Chandran2017b), occurring approximately at the electron inertial length scale $d_e \equiv c/\omega _{pe} = \rho _i/(\beta _i \mu )^{1/2}$. These scales are ordered by $\delta _{in} \sim d_e \ll \rho _s \ll a$. Together, these restrictions limit us to moderately small plasma beta

(5.29)\begin{equation} \mu^{{-}1} \ll\beta_i/\tau \ll 1.\end{equation}

The tearing instability depends on the tearing mode instability parameter $\varDelta '$, determined by the outer-region solution (at MHD scales $k_\perp \rho _i \ll 1$) (Zocco & Schekochihin Reference Zocco and Schekochihin2011). This parameter can be expressed over two limiting regimes in a dimensionless form (Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a) as

(5.30)\begin{equation} \Delta' a = (ka)^{{-}n} \quad \mbox{for}\ 1 < n \le 2, \end{equation}

where $k$ is the unstable wavenumber along the width of the current layer governing the reconnection flow, as illustrated in figure 9 by the red sinusoidal curve. Here $n=1$ corresponds to a Harris current sheet (a hyperbolic tangent profile of the reconnecting magnetic field across the current sheet thickness, Harris Reference Harris1962), while $n=2$ corresponds to a sinusoidal variation of the reconnection magnetic field (Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a), which may be the more likely case for the onset of reconnection at current sheets that self-consistently arise in a turbulent plasma. The tearing instability growth rates are solved in two limits of $\varDelta '$ (Zocco & Schekochihin Reference Zocco and Schekochihin2011; Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a; Mallet et al. Reference Mallet, Schekochihin and Chandran2017b): (i) $\Delta ' \delta _{in} \ll 1$

(5.31)\begin{equation} \gamma \sim k v_{A \perp} \frac{d_e \rho_s \varDelta'}{a} \quad \mbox{with}\ \delta_{in} \sim d_e \rho_s^{1/2} \varDelta'^{1/2}, \end{equation}

and (ii) $\Delta ' \delta _{in} \gtrsim 1$

(5.32)\begin{equation} \gamma \sim k v_{A \perp} \frac{d_e^{1/3} \rho_s^{2/3}}{a} \quad \mbox{with}\ \delta_{in} \sim d_e^{2/3} \rho_s^{1/3}, \end{equation}

where $v_{A \perp } \equiv \delta B_\perp /(4 {\rm \pi}n_i m_i)^{1/2}$ is the Alfvén speed based on the magnitude of only the reconnecting component of the magnetic field $\delta B_\perp$. Since the growth rates above are calculated in the limit $k_\perp \rho _i \ll 1$, we will focus our assessment on whether magnetic reconnection will disrupt the turbulent cascade at MHD scales $k_\perp \rho _i \lesssim 1$, which is the same approach as taken in previous studies (Loureiro & Boldyrev Reference Loureiro and Boldyrev2017a; Mallet et al. Reference Mallet, Schekochihin and Chandran2017b).

Using (5.30) and either (5.31) or (5.32) for $\delta _{in}$, the normalized tearing mode instability parameter $\Delta ' \delta _{in}$ can be expressed in terms of the dimensionless quantities in table 4 by

(5.33)\begin{equation} \Delta' \delta_{in} = \frac{k_\perp \rho_i}{(2 \tau \beta_i^2 \mu^2)^{1/6}} \left(\frac{k_\perp \rho_i}{k \rho_i}\right)^n. \end{equation}

This equation shows that $\Delta ' \delta _{in}$ increases as the perpendicular wavenumber $k_\perp \rho _i$ increases, corresponding to a thinner current sheet, or as the unstable wavenumber $k\rho _i$ decreases, requiring to a wider current sheet. Since the width $w\sim 1/k_i$ of the current sheets in B06 scaling theory governed by (3.2) bounds the unstable wavenumber $k$ from below by the intermediate wavenumber, $k>k_i$, the maximum value of instability parameter $\Delta ' \delta _{in}$ increases with the perpendicular aspect ratio of the current sheets $w/a = k_\perp /k_i \gg 1$. Note also that $\Delta ' \delta _{in} \propto \beta _i^{-1/3}$, so that lower values of ion plasma beta also lead to a larger instability parameter.

Note that the value of the unstable wavenumber $k \rho _i$ at the transition between the limiting regimes of the tearing mode instability parameter, where $\Delta ' \delta _{in}=1$, is given by

(5.34)\begin{equation} k \rho_i = k_t \rho_i \equiv \frac{(k_\perp \rho_i)^{(n+1)/n}}{(2 \tau \beta_i^2 \mu^2)^{1/6n}}, \end{equation}

which provides a bound on the value of $k \rho _i$ for the validity of the tearing growth rate formulas in (5.31) and (5.32).

To assess whether the collisionless tearing instability grows fast enough to enable the onset of magnetic reconnection to disrupt the turbulent cascade, we will assess $\gamma _{{\rm RXN}}/\omega _{nl} \sim \gamma _{{\rm RXN}}/\omega$, assuming critical balance of linear and nonlinear time scales $\omega _{nl}\sim \omega$ (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Mallet et al. Reference Mallet, Schekochihin and Chandran2015). Taking the linear wave frequency for MHD Alfvén waves $\omega = k_\parallel v_A = (k_\parallel \rho _i)^{1/2} (k_\perp \rho _i)^{1/2} v_A/\rho _i$ using (3.1) with $\alpha =1$ for the B06 scaling of $k_\parallel$, and substituting $v_{A \perp } = v_A (\delta B_\perp /B_0) = v_A [(k_0 \rho _i)/(k_\perp \rho _i)]^{1/4}$ using (3.4), we find that the growth rate in the $\Delta ' \delta _{in} \ll 1$ limit, which corresponds to the large unstable wavenumber $k$ limit by inspection of (5.33), can be expressed as

(5.35)\begin{equation} \frac{\gamma_{{\rm RXN}}}{\omega} \sim \frac{(k_\perp \rho_i)^{(5+4n)/4}}{(k \rho_i)^{n-1} (k_0 \rho_i)^{1/4}(2 \tau \beta_i \mu)^{1/2}}, \end{equation}

where this expression is formally valid in the MHD limit with $k_\perp \rho _i \ll 1$. Combining the limitations on the possible unstable wavenumber from (5.28) and the requirement $\Delta ' \delta _{in} \lesssim 1$ using (5.34), the formula in (5.35) is valid over the range of normalized unstable wavenumbers

(5.36)\begin{equation} \max\left[(k_0 \rho_i)^{1/4}(k_\perp \rho_i)^{3/4}, \frac{(k_\perp \rho_i)^{(n+1)/n}}{(2 \tau \beta_i^2 \mu^2)^{1/6n}} \right] \lesssim k \rho_i \ll k_\perp \rho_i \lesssim 1. \end{equation}

Note that this valid unstable wavenumber range can be used to show that no range with $\Delta ' \delta _{in} \lesssim 1$ exists if $\beta _i \lesssim (2 \tau \mu ^2)^{-1/2}$. For typical values of $\tau \sim 1$, however, the very low values of $\beta _i$ required to violate this condition would already have violated the assumed moderate beta limit given by (5.29), so this limitation is of no practical concern.

The tearing growth rate in (5.35) demonstrates that, in the large $k$ limit, the normalized growth rate is independent of $k$ with $\gamma _{{\rm RXN}}/\omega \propto k^0$ for Harris-type current sheets with $n=1$, and is inversely proportional to $k$ with $\gamma _{{\rm RXN}}/\omega \propto k^{-1}$ for the $n=2$ case of sinusoidal variations of the perpendicular magnetic field component that may be expected to occur in plasma turbulence.

The maximum growth rate $\gamma _{{\rm RXN}}/\omega$ from (5.35) increases with $k_\perp \rho _i$, but also requires the limit $k_\perp \rho _i \ll 1$ for the validity of this linear tearing growth-rate calculation. Therefore, to estimate quantitatively the maximum normalized tearing growth rate, we take $k_\perp \rho _i \rightarrow \mathcal {F}$, where $\mathcal {F} \ll 1$ is the maximum value of the normalized perpendicular wavenumber for which (5.35) remains valid; numerical simulations of collisionless magnetic reconnection will be invaluable to provide estimates of the maximum value of $\mathcal {F}$. The maximum growth rate $\gamma _{{\rm RXN}}/\omega$ also occurs for the minimum unstable wavenumber $k\rho _i$ allowed by the range in (5.36). Thus, taking these two limits and requiring a collisionless tearing instability growth rate sufficient to disrupt the turbulent cascade, $\gamma _{{\rm RXN}}/\omega \gtrsim 1$, yields the constraint

(5.37)\begin{equation} k_0 \rho_i \lesssim \frac{\mathcal{F}^{(4+5n)/n}}{(2 \tau )^{(4n+2)/3n} (\beta_i \mu)^{(2n+4)/3n}}, \end{equation}

if the minimum $k\rho _i$ is constrained by the condition $\Delta ' \delta _{in} \lesssim 1$. Alternatively, if the minimum $k\rho _i$ is constrained by the width of the current sheet given by (5.28), then the condition $\gamma _{{\rm RXN}}/\omega \gtrsim 1$ yields the constraint

(5.38)\begin{equation} k_0 \rho_i \lesssim \frac{\mathcal{F}^{(5+4n)/n}}{ (2 \tau \beta_i \mu)^{2/n}}. \end{equation}

The normalized growth rates in the $\Delta ' \delta _{in} \lesssim 1$ limit vs. normalized unstable wavenumber $k \rho _i$ are illustrated in figure 10(a) for $n=1$ (thin blue) and $n=2$ (thick blue) for turbulent parameters $\beta _i=0.01$, $\tau =1$, $\mu =1836$ and $k_0\rho _i=10^{-4}$ and for three values of the maximum of $k_\perp \rho _i$ given by $\mathcal {F}=1.0$ (solid), $0.3$ (dashed) and $0.1$ (dotted). The corresponding values of $\Delta ' \delta _{in}$ vs. $k \rho _i$ are presented in figure 10(b) for each of the $n$ and $k_\perp \rho _i$ cases. For the turbulence parameters specified in this example, we find $\gamma _{{\rm RXN}}/\omega \gtrsim 1$ only for the $\mathcal {F} =k_\perp \rho _i=1.0$ case over the unstable wavenumber range $0.3 \lesssim k\rho _i \lesssim 1$ for the $n=1$ case (thin solid blue) and over $0.6 \lesssim k\rho _i \lesssim 1$ for the $n=2$ case (thick solid blue).

Figure 10. (a) Scaling the normalized collisionless tearing instability growth rate $\gamma _{{\rm RXN}}/\omega$ for a turbulent plasma with $\beta _i=0.01$, $\tau =1$, $\mu =1836$ and $k_0\rho _i=10^{-4}$ as a function of the normalized unstable wavenumber $k \rho _i$. Colours indicate the solutions in the low $k$ limit (corresponding to $\Delta ' \delta _{in} \gtrsim 1$) (red) and in the high $k$ limit (corresponding to $\Delta ' \delta _{in} \ll 1$) (blue). Line thickness indicates the either a sinusoidal current profile with $n=2$ (thick) or a Harris-like hyperbolic tangent current profile with $n=1$ (thin). Line styles indicate different values of maximum perpendicular wavenumber $\mathcal {F} = 1.0$ (solid), $0.3$ (dashed) and $0.1$ (dotted). (b) The corresponding normalized tearing mode instability parameter $\Delta ' \delta _{in}$ vs. $k \rho _i$ for each case.

Next, we assess the initiation of collisionless tearing in the $\Delta ' \delta _{in} \gtrsim 1$ limit with the growth rate given by (5.32). Using similar substitutions for $\omega$ and $v_{A \perp }$ (as detailed above for the $\Delta ' \delta _{in} \ll 1$ limit), we find that the growth rate in the $\Delta ' \delta _{in} \gtrsim 1$ limit, which is the small unstable wavenumber $k$ limit as shown by (5.33), can be expressed as

(5.39)\begin{equation} \frac{\gamma_{{\rm RXN}}}{\omega} \sim \frac{(k \rho_i) (k_\perp \rho_i)^{1/4}} { (k_0 \rho_i)^{1/4}(4 \tau^2 \beta_i \mu)^{1/6}}, \end{equation}

where this expression is formally valid in the MHD limit with $k_\perp \rho _i \ll 1$. Combining the limitations on the possible unstable wavenumber from (5.28) and the requirement $\Delta ' \delta _{in} \gtrsim 1$ using (5.34), the formula in (5.39) is valid over the range of normalized unstable wavenumbers

(5.40)\begin{equation} (k_0 \rho_i)^{1/4}(k_\perp \rho_i)^{3/4} \lesssim k \rho_i \ll \max\left[k_\perp \rho_i,\frac{(k_\perp \rho_i)^{(n+1)/n}}{(2 \tau \beta_i^2 \mu^2)^{1/6n}} \right]\lesssim 1. \end{equation}

This valid unstable wavenumber range can be used to show that no range with $\Delta ' \delta _{in} \gtrsim 1$ exists if

(5.41)\begin{equation} k_\perp \rho_i \lesssim (k_0 \rho_i)^{n/(n+4)}(2 \tau \beta_i^2 \mu^2)^{2/3(n+4)}. \end{equation}

This finding is consistent with the fact that the value of $\Delta ' \delta _{in}$ decreases with decreasing $k_\perp \rho _i$ in (5.33), so by decreasing $k_\perp \rho _i$ enough, one will always reach a point where there is no range of $k\rho _i$ with $\Delta ' \delta _{in} \gtrsim 1$.

The tearing growth rate in (5.39) shows that, in the small $k$ limit, the normalized growth rate increases linearly with the unstable wavenumber $k$ and increases weakly with the perpendicular wavenumber $k_\perp \rho _i$, independent of whether the turbulently generated current sheets are Harris-like or sinusoidal (thus, independent of $n$).

Once again the maximum growth rate $\gamma _{{\rm RXN}}/\omega$ from (5.39) increases with $k_\perp \rho _i$ (although weakly with the $1/4$ power), so once again we take the limit $k_\perp \rho _i \rightarrow \mathcal {F}$, where $\mathcal {F} \ll 1$ is the maximum value of the normalized perpendicular wavenumber for which (5.39) remains valid. The maximum growth rate will also occur for the maximum unstable wavenumber $k\rho _i$ given by (5.40). Thus, taking these two limits and requiring a collisionless tearing instability growth rate sufficient to disrupt the turbulent cascade, $\gamma _{{\rm RXN}}/\omega \gtrsim 1$, yields the constraint

(5.42)\begin{equation} k_0 \rho_i \lesssim \frac{\mathcal{F}^{(4+5n)/n}} {(2 \tau )^{(4n+2)/3n} (\beta_i \mu)^{(2n+4)/3n}},\end{equation}

if the maximum $k\rho _i$ is constrained by the condition $\Delta ' \delta _{in} \gtrsim 1$. Note that this condition is the same as (5.37), as it must be since this maximum growth rate occurs at $\Delta ' \delta _{in} = 1$, where the growth rates for the $\Delta ' \delta _{in} \lesssim 1$ and $\Delta ' \delta _{in} \gtrsim 1$ limits cross. Alternatively, if the maximum $k\rho _i$ is constrained by the maximum $k_\perp \rho _i \rightarrow \mathcal {F}$, then the condition $\gamma _{{\rm RXN}}/\omega \gtrsim 1$ yields the constraint

(5.43)\begin{equation} k_0 \rho_i \lesssim \frac{\mathcal{F}^{(5+4n)/n}}{(4 \tau^2 \beta_i \mu)^{2/3}}. \end{equation}

The normalized growth rates in the $\Delta ' \delta _{in} \gtrsim 1$ limit vs. normalized unstable wavenumber $k \rho _i$ are illustrated in figure 10(a) (thick red) for turbulent parameters $\beta _i=0.01$, $\tau =1$, $\mu =1836$ and $k_0\rho _i=10^{-4}$ and for three values of the maximum of $k_\perp \rho _i$ given by $\mathcal {F}=1.0$ (solid), $0.3$ (dashed) and $0.1$ (dotted). The corresponding values of $\Delta ' \delta _{in}$ vs. $k \rho _i$ are presented in figure 10(b) for each of the $n$ and $k_\perp \rho _i$ cases. For the turbulence parameters specified in this example, we find $\gamma _{{\rm RXN}}/\omega \gtrsim 1$ only for the $k_\perp \rho _i=1.0$ case over the unstable wavenumber range $0.2 \lesssim k\rho _i \lesssim 0.3$ for the $n=1$ case (thick solid red) and over $0.2 \lesssim k\rho _i \lesssim 0.6$ for the $n=2$ case (thick solid red).

It is necessary to highlight here a few important caveats in using the linear collisionless tearing mode growth rates and MHD turbulence scaling theory to predict when magnetic reconnection may disrupt the dynamics of the turbulent cascade. First, the linear tearing mode growth rates used here are rigorously derived in asymptotic limits, but we employ these same rates up to the extreme boundary of those limits. For example, the linear growth rates formally require $k_\perp \rho _i \ll 1$, but we take $k_\perp \rho _i \rightarrow 1$ to estimate the maximum growth rate; also, the $\Delta ' \delta _{in} \ll 1$ limit is used to estimate the tearing mode growth rates as $\Delta ' \delta _{in} \rightarrow 1$. Such rates derived in asymptotic limits, however, may not remain quantitatively accurate as we extend beyond their formal regime of validity. We adopt the viewpoint that the growth rates calculated in these limits are at least reasonable lowest-order estimates of the rates even beyond their formally limited range. Furthermore, the growth rates for the collisionless tearing mode are determined under smooth conditions, whereas to determine whether reconnection arises in the current sheets that naturally arise in plasma turbulence, one needs to determine the collisionless tearing growth rates in the presence of turbulence; once again, the estimated tearing mode growth rates used here are taken simply as the lowest-order estimation. Computing the tearing mode growth rates in a turbulent plasma and determining the competition of the tearing mode with the nonlinear cascade of energy due to turbulence are both necessary for more accurate calculations, which likely will require detailed kinetic numerical simulations. Finally, alternative proposed mechanisms for the development of current sheets in plasma turbulence – such as the nonlinear interaction of counterpropagating Alfvén wave collisions (Howes Reference Howes2016), as opposed to dynamic alignment (Boldyrev Reference Boldyrev2006) – may yield a different scaling for the width and thickness of current sheets; the resulting predictions for the onset of magnetic reconnection in a turbulent plasma would thereby also likely change.