# [

###### Abstract

We present a theoretical framework for describing electromagnetic kinetic turbulence in a multi-species, magnetized, pressure-anisotropic plasma. The turbulent fluctuations are assumed to be small compared to the mean field, to be spatially anisotropic with respect to it, and to have frequencies small compared to the ion cyclotron frequency. At scales above the ion Larmor radius, the theory reduces to the pressure-anisotropic generalization of kinetic reduced magnetohydrodynamics (KRMHD) formulated by Kunz et al., 2015, J. Plasma Phys., vol. 81, 325810501. At scales at and below the ion Larmor radius, three main objectives are achieved. First, we analyse the linear response of the pressure-anisotropic gyrokinetic system, and show it to be a generalisation of previously explored limits. The effects of pressure anisotropy on the stability and collisionless damping of Alfvénic and compressive fluctuations are highlighted, with attention paid to the spectral location and width of the frequency jump that occurs as Alfvén waves transition into kinetic Alfvén waves. Secondly, we derive and discuss a very general gyrokinetic free-energy conservation law, which captures both the KRMHD free-energy conservation at long wavelengths and dual cascades of kinetic Alfvén waves and ion entropy at sub-ion-Larmor scales. We show that non-Maxwellian features in the distribution function change the amount of phase mixing and the efficiency of magnetic stresses, and thus influence the partitioning of free energy amongst the cascade channels. Thirdly, a simple model is used to show that pressure anisotropy, even within the bounds imposed on it by firehose and mirror instabilities, can cause order-of-magnitude variations in the ion-to-electron heating ratio due to the dissipation of Alfvénic turbulence. Our theory provides a foundation for determining how pressure anisotropy affects the turbulent fluctuation spectra, the differential heating of particle species, and the ratio of parallel and perpendicular phase mixing in space and astrophysical plasmas.

Kinetic turbulence in pressure-anisotropic plasmas]Astrophysical gyrokinetics:

Turbulence in pressure-anisotropic plasmas at ion scales and beyond
M. W. Kunz and others]M. W. Kunz^{†}^{†}thanks: Email address for correspondence: mkunz@princeton.edu,
I. G. Abel,
K. G. Klein,
and A. A. Schekochihin

## 1 Introduction

In a previous paper (Kunz et al. 2015, hereafter Paper I), we presented a theoretical framework for low-frequency electromagnetic (drift-)kinetic turbulence valid at scales larger than the particles’ Larmor radii (“long” wavelengths) in a collisionless, multi-species plasma. That result generalised reduced magnetohydrodynamics (RMHD; Kadomtsev & Pogutse 1974; Strauss 1976, 1977; Zank & Matthaeus 1992) and kinetic RMHD (Schekochihin et al. 2009, hereafter S09) to the case where the mean distribution function of the plasma is pressure-anisotropic and different ion species are allowed to drift with respect to each other – a situation routinely encountered in the solar wind (e.g. Hundhausen et al. 1967; Feldman et al. 1973; Marsch et al. 1982a, b; Marsch 2006) and presumably ubiquitous in hot dilute astrophysical plasmas such as the intracluster medium of galaxy clusters (e.g. Schekochihin et al. 2005; Schekochihin & Cowley 2006). This framework was obtained via two routes: one starting from Kulsrud’s formulation of kinetic MHD (Kulsrud 1964, 1983) and one starting from applying the nonlinear gyrokinetic reduction (e.g. Frieman & Chen 1982; Howes et al. 2006) of the Vlasov-Maxwell set of equations. The latter approach also enables the study of fluctuations at and below the ion Larmor scale, the subject of this Paper.

Before embarking on any quantitative analysis or even qualitative discussion of what the gyrokinetic framework entails, we catalogue the principal theoretical achievements and implications of Paper I. First, we showed that the main physical feature of low-frequency, long-wavelength plasma turbulence survives the generalisation to non-Maxwellian equilibrium distribution functions: Alfvénic and compressive fluctuations are energetically decoupled, with the latter passively advected by the former. The Alfvénic cascade is fluid, satisfying RMHD equations (with the Alfvén speed modified by pressure anisotropy and interspecies drifts), whereas the compressive cascade is kinetic and subject to collisionless damping. For a bi-Maxwellian plasma, the kinetic cascade splits into three independent collisionless cascades. Secondly, the organising principle of this long-wavelength turbulence was elucidated in the form of a conservation law for the appropriately generalised kinetic free energy. Using this alongside linear theory, we showed that certain non-Maxwellian features in the distribution function can reduce the rate of collisionless damping and the efficacy of magnetic stresses, and that these changes influence the partitioning of free energy amongst the various cascade channels. As the firehose or mirror instability thresholds are approached, the dynamics of the plasma are modified so as to reduce the energetic cost of bending magnetic-field lines or of compressing them.

In this paper, we concentrate on the sub-Larmor-scale “dissipation” range. We investigate the linear properties of kinetic Alfvén waves (KAWs) and their nonlinear phase-space cascade in a plasma whose mean particle distribution functions exhibit pressure anisotropy and interspecies drifts. We find that the stability conditions imposed on the KAWs by the anisotropy of the distribution functions are a combination of those experienced by Alfvén waves and compressive fluctuations at long wavelengths. We further show that, similar to the dual Alfvénic-kinetic cascade of free energy in the inertial range (Paper I), there are two sub-ion-Larmor-scale kinetic cascades: one of KAWs, which is governed by a set of fluid-like electron reduced magnetohydrodynamic equations, and a passive phase-space cascade of ion-entropy fluctuations. These cascades have already been considered for a single-ion-species plasma whose equilibrium distribution function is isotropic and Maxwellian (S09). Here we focus on whether and to what extent their results carry over to the more general case. Special attention is paid to the transition from the inertial range across the ion-Larmor scale to the kinetic range, and to the effect of pressure anisotropy on the spectral location of this transition and on the amount of Landau-damped energy that ultimately makes its way to collisional scales.

All sections and equations in Paper I are referenced using the prefix “I–”; e.g. (I–C1) refers to equation (C1) of Paper I and §I–4 refers to section 4 of Paper I.

## 2 Prerequisites

### 2.1 Basic equations and notation

For completeness, we provide here the basic equations derived in Paper I from which this paper’s results follow, as well as the notation introduced in Paper I by which this paper’s results may be understood.^{1}^{1}1A glossary of frequently used symbols can be found in appendix E of Paper I. This recapitulation starts with the Vlasov-Landau equation,

(2.0) |

governing the space-time evolution of the particle distribution function of species , , where is the velocity-space variable and is the real-space variable. The charge and mass of species are denoted and , respectively; is the speed of light. The electric field and magnetic field are expressed in terms of scalar and vector potentials:

(2.0) |

where is the guide magnetic field, taken to lie along the axis, and (the Coulomb gauge). These fields satisfy the plasma quasineutrality constraint,

(2.0) |

and the pre-Maxwell version of Ampère’s law,

(2.0) |

where and are the number density and mean velocity of species and is the current density.

The term on the right-hand side of (2.1) represents the effect of collisions on the distribution function; in this paper, collisions are assumed to be sub-dominant and thus its specific form will not be required (precisely what ‘sub-dominant’ means will be stated in short order). The assumption of weak collisionality gives the pressure tensor

(2.0) |

the freedom to be anisotropic, even in the mean (zeroth-order) background. An example of such a pressure tensor is that describing a gyrotropic plasma (see §2.3),

(2.0) |

where is the unit dyadic, is the unit vector in the direction of the magnetic field, the subscript () denotes the component perpendicular (parallel) to , and

(2.0) | |||

(2.0) |

are the parallel and perpendicular pressures, respectively, of species . An oft-employed distribution function that exhibits such pressure anisotropy is the bi-Maxwellian

(2.0) |

where

(2.0) |

are the parallel and perpendicular thermal speeds of species . Pressure anisotropy is caused in a weakly collisional plasma by adiabatic invariance: conservation of the magnetic moment implies that a slow change in magnetic-field strength must be accompanied by a proportional change in the perpendicular temperature of species (Chew et al. 1956). While such velocity-space anisotropy is generically exhibited by the gyrokinetic fluctuations regardless of whether the mean distribution function is proved (or assumed) to be isotropic and Maxwellian, in what follows we also allow for the possibility of a background pressure anisotropy.

### 2.2 Gyrokinetic ordering

Our aim is to reduce (2.1)–(2.1) so that they describe only those fields whose fluctuating parts are small compared to the mean field, are spatially anisotropic with respect to it, have frequencies small compared to the Larmor frequency , and have parallel length scales large compared to the Larmor radius . While such specifications may appear to be quite restrictive, modern theories (e.g. Goldreich & Sridhar 1995) and numerical simulations (e.g. Shebalin et al. 1983; Oughton et al. 1994; Cho & Vishniac 2000; Maron & Goldreich 2001) of magnetized turbulence provide a strong foundation for expecting such anisotropic low-frequency fluctuations to comprise much of the energy in the turbulent cascade. Such spatial anisotropy is also now routinely measured in the solar wind (e.g. Bieber et al. 1996; Horbury et al. 2008; Podesta 2009; Wicks et al. 2010; Chen et al. 2011; Chen 2016) and suggested by observations of turbulent density fluctuations in the interstellar medium (e.g. Armstrong et al. 1990; Rickett et al. 2002).

The reduction is carried out in detail in appendix C of Paper I; here we describe its primary ingredients and principal consequences. The fields are split into their mean parts (denoted with a subscript ‘0’) and fluctuating parts (denoted with ), the former characterized by spatial homogeneity. The latter are taken to satisfy the asymptotic ordering

(2.0) |

where we have expanded the distribution function in powers of :

(2.0) |

Note that the fluctuations are permitted to have perpendicular scales on the order of the Larmor radius. We further assume that the collision frequency , thereby allowing non-Maxwellian (cf. §A2.2 of Howes et al. 2006).

The gyrokinetic ordering guarantees that (to lowest order) all species drift perpendicularly to the magnetic field with identical velocities, . It then follows that the mean drift of any species relative to the centre-of-mass velocity must be in the parallel direction, viz., , with

(2.0) |

Our collisionless ordering permits parallel interspecies drifts (denoted by ) in the background state, and we formally order for all species . We further assume that the Alfvén speed

(2.0) |

where is the mean mass density of the plasma. This implies that the parallel and perpendicular plasma beta parameters,

(2.0) |

respectively, are considered to be of order unity in the gyrokinetic expansion. The other dimensionless parameters in the system – namely, the electron-ion mass ratio , the charge ratio , the parallel and perpendicular temperature ratios

(2.0) |

and the temperature anisotropy of species – are all considered to be of order unity as well. Subsidiary expansions with respect to these parameters can (and will) be made after the gyrokinetic expansion is performed.

Since we have , fast magnetosonic fluctuations are ordered out of our equations. Such fast-wave fluctuations are rarely seen in the solar wind (Howes et al. 2012). Observations of turbulence in the solar wind confirm that it is primarily Alfvénic (e.g. Belcher & Davis 1971; Chen 2016) and that its compressive component is approximately pressure-balanced (Burlaga et al. 1990; Roberts 1990; Marsch & Tu 1993; McComas et al. 1995; Bavassano et al. 2004; Bruno & Carbone 2005). A more serious limitation of our analysis is perhaps the exclusion of cyclotron resonances, which have been traditionally considered necessary to explain the strong perpendicular heating observed in the solar wind (Leamon et al. 1998; Isenberg 2001; Kasper et al. 2013). Larmor-scale fluctuations whose amplitudes are large enough to break adiabatic invariance and thus drive chaotic gyromotion and stochastic particle heating (Chandran et al. 2010, 2013) are also precluded. That being said, the gyrokinetic framework does capture much of the physics governing both the inertial and dissipative ranges of kinetic turbulence, and so it is a sensible step to incorporate realistic background distribution functions into the gyrokinetic description of weakly collisional astrophysical plasmas. It is with that goal in mind that we commence with a presentation of the gyrokinetic theory.

### 2.3 Gyrokinetic reduction

#### 2.3.1 Gyrotropy of the background distribution function

Under the ordering (2.2), the largest term in the Vlasov-Landau equation (2.1) corresponds to Larmor motion of the mean distribution about the uniform guide field:

(2.0) |

This directional bias allows us to set up a local Cartesian coordinate system and decompose the particle velocity in terms of the parallel velocity , the perpendicular velocity , and the gyrophase angle ,

(2.0) |

Equation (2.3.1) then takes on the simple form

(2.0) |

which states that the mean distribution function is gyrotropic (independent of the gyrophase):

(2.0) |

All velocity-space derivatives of that enter (2.1) are thus with respect to and , viz.

(2.0) |

where

(2.0) |

are dimensionless derivatives of a species’ mean distribution function with respect to the square of the parallel velocity (peculiar to the species drift velocity) and the perpendicular velocity, respectively. Their weighted difference,

(2.0) |

measures the velocity-space anisotropy of the mean distribution function. For a bi-Maxwellian distribution (2.1), and

(2.0) |

where is the temperature (or, equivalently, pressure) anisotropy of the mean distribution function of species .

#### 2.3.2 Boltzmann response

At , we learn from (2.1) that the first-order distribution function may be split into two parts. The first of these is the so-called adiabatic (or ‘Boltzmann’) response,

(2.0) |

where

(2.0) |

is the fluctuating electrostatic potential in the frame of the parallel-drifting species . This part of represents the (leading-order) evolution of under the influence of the perturbed electromagnetic fields. To see this, we first introduce the total particle energy in the parallel-drifting frame,

(2.0) |

and the (gyrophase-dependent part of the) first adiabatic invariant,

(2.0) |

both written out to first order in the fluctuation amplitudes (e.g. Kruskal 1958; Hastie et al. 1967; Taylor 1967; Catto et al. 1981; Parra 2013). It is then straightforward to show by using

(2.0) |

(see §I–C.4) that the sum of the mean distribution function and the Boltzmann response is simply

(2.0) |

In other words, the Boltzmann response does not change the form of the mean distribution function if the latter is written as a function of the constants of the motion (calculated sufficiently accurately).

#### 2.3.3 Gyrokinetic response

The second part of , which we denote by , represents the response of rings of charge to the fluctuating fields, and is thus referred to as the gyrokinetic response. It satisfies

(2.0) |

where we have transformed the derivative taken at constant position to one taken at constant guiding centre

(2.0) |

Thus, is independent of the gyrophase angle at constant guiding centre (but not at constant position ):

(2.0) |

#### 2.3.4 Gyrokinetic equation

At , we find from (2.1) that the gyrokinetic response evolves via the gyrokinetic equation

(2.0) |

where

(2.0) |

is the gyrokinetic potential and

(2.0) |

denotes the ring average of at fixed guiding centre . The Poisson bracket

(2.0) |

represents the nonlinear interaction between the gyrocentre rings and the ring-averaged electromagnetic fields.

The gyrokinetic equation (2.3.4) can also be written in the following, perhaps more physically illuminating, form:

(2.0) |

where

(2.0) |

is the ring velocity,

(2.0) |

is the ring-averaged rate of change of the particle energy (2.3.2), and

(2.0) |

is the ring-averaged rate of change of the (gyrophase-dependent part of the) first adiabatic invariant (2.3.2). The right-hand side of (2.0) represents the effect of collisionless work done on the rings by the fields (the wave-ring interaction). Written in this way, (2.3.4) is simply the ring-averaged Vlasov equation,

(2.0) |

to lowest order in .

It is a manifestly good idea in much of what follows to absorb the final term of (2.3.4) (and, likewise, of (2.0)), into by writing the latter in terms of the velocity-space coordinates , where

(2.0) |

is the full adiabatic invariant, viz., . At long wavelengths satisfying ,

(2.0) |

which is simply the magnetic moment of a particle in a magnetic field of strength drifting across said field at the velocity,

(2.0) |

Then, introducing^{2}^{2}2Our is equivalent to of Frieman & Chen (1982) – see their equation (42), with their being our .

(2.0a) | ||||

(2.0b) |

the gyrokinetic equation reads

(2.0) |

This form of the gyrokinetic equation is particularly well suited for deriving the gyrokinetic invariants (§4). Its right-hand side represents the collisionless work done on the rings by the fields in a frame comoving with the parallel drift velocity of species .

It will also prove useful in what follows to modify the energy variable to obtain

(2.0) |

which is the kinetic energy of the particle as measured in the frame moving with the and drifts (e.g. Parra 2013); indeed,

(2.0) |

at long wavelengths. If the mean distribution function is expressed in terms of these new velocity-space variables, viz. , then the perturbed distribution function becomes (see (I–C52))

(2.0a) | ||||

(2.0b) |

This particular form of the perturbed distribution function is quite useful; it is the generalisation of the perturbed distribution function that prominently features in the generalised free energy of KRMHD (§I–5.1), and thus is anticipated to appear in the generalised free energy of the gyrokinetic theory. The latter is derived in §4.

#### 2.3.5 Field equations

The equations governing the electromagnetic potentials are most easily obtained by substituting the decomposition

(2.0) |

into the leading-order expansions of the quasineutrality constraint (2.1) and Ampère’s law (2.1). The result is (see §I–C.3)

(2.0) |

(2.0) |

(2.0) |

where

(2.0) |

is the temperature anisotropy of species augmented by the parallel ram pressure from background parallel drifts, are parallel moments of the perpendicular-differentiated mean distribution function (all of which equate to unity for a drifting bi-Maxwellian distribution; see appendix A), and

(2.0) |

denotes the gyro-average of at fixed . Together with the gyrokinetic equation (2.3.4), the field equations (2.0)–(2.0) constitute a closed system that describes the evolution of a gyrokinetic plasma with non-Maxwellian and parallel interspecies drifts.

This completes our abbreviated review of the material derived in Paper I on the gyrokinetic framework for homogeneous, non-Maxwellian plasmas. We now proceed to analyse the linear and nonlinear behaviour of the perturbations governed by this system of equations.

## 3 Linear gyrokinetic theory

### 3.1 From rings to gyrocentres

The most straightforward way of making contact with the results of Paper I, while facilitating the extension of the theoretical framework into the kinetic range, is via the linear gyrokinetic theory (pioneered by Rutherford & Frieman 1968; Taylor & Hastie 1968; Catto 1978; Antonsen & Lane 1980; Catto et al. 1981). This is obtained most easily by shifting the description of the plasma from one composed of extended rings of charge that move in a vacuum to one of a gas of point-particle-like gyrocentres moving in a polarizable medium. This transformation is enacted by working with the gyrocentre distribution function

(3.0a) | ||||

(3.0b) | ||||

(3.0c) |

This new function not only helps simplify the algebra involved in deriving the linear theory, but also makes a good deal of physical sense. In the electrostatic limit, the use of (which, in this limit, equals ) aids in the interpretation of polarization effects within gyrokinetics (Krommes 2012), places the gyrokinetic equation in a numerically convenient characteristic form (Lee 1983), and arises naturally from the Hamiltonian formulation of gyrokinetics (Dubin et al. 1983; Brizard & Hahm 2007). In the electromagnetic case, introducing takes advantage of the fact that the Alfvénic fluctuations have a gyrokinetic response that is approximately cancelled at long wavelengths by the Boltzmann response (see §I–C.4), i.e., for long-wavelength Alfvénic fluctuations.

Using (3.1) to replace in the gyrokinetic equation (2.0), we find that evolves according to

(3.0) |

where

(3.0) |

is the spatial derivative along the perturbed magnetic field and

(3.0) |

We have used compact notation in writing out the nonlinear terms: , where the first Poisson bracket involves derivatives with respect to and the second with respect to . We now develop the linear theory.

### 3.2 Linear gyrokinetic equation

We begin by linearizing the gyrokinetic equation (3.1) in the fluctuations’ amplitudes:

(3.0) |

Decomposing the perturbed distribution function and the fluctuating electromagnetic potentials and into plane-wave solutions,

and substituting these expressions into (3.0), we find that

(3.0) |

where and are, respectively, the zeroth- and first-order Bessel functions of (cf. equation I–B1). The Bessel functions arise from performing the ring averages in the Fourier space (see §I–C6 for details).

### 3.3 Gyrokinetic field equations

Next, we insert (3.0) into the field equations (2.0)–(2.0). This procedure involves computing several -, -, and Bessel-function–weighted Landau-like integrals over the mean distribution function. These integrals (denoted and for integer and ) are defined in appendix A and evaluated to leading order in . Using these definitions, the quasineutrality constraint (2.0) and the parallel (2.0) and perpendicular (2.0) components of Ampère’s law may be written, respectively, as

(3.0) |