TOWARD A THEORY OF ASTROPHYSICAL PLASMA TURBULENCE AT SUBPROTON SCALES

The simplest description of plasma dynamics on large scales (compared to microscales such as ion gyroradius, skin depth, etc.) is provided by one-fluid magnetohydrodynamics (MHD). In the presence of a background magnetic field, ideal MHD possesses exact nonlinear solutions—Alfvén waves propagating up and down along the magnetic field lines (e.g., Biskamp 2003; Kulsrud 2005). Interactions of counter-propagating waves redistribute energy over scales and, when the dissipation is not significant, give rise to a turbulent cascade by which energy is transferred to progressively smaller scales. It has been established that such an energy cascade is intrinsically anisotropic, in that it predominantly supplies energy to the modes with mostly field-perpendicular wavenumbers, $k_\perp \gg k_\Vert$ (e.g., Shebalin et al. 1983; Goldreich & Sridhar 1995; Galtier et al. 2000; Cho & Vishniac 2000; Boldyrev & Perez 2009; Perez & Boldyrev 2010; Perez et al. 2012).

When the energy cascade reaches particle gyroscales, the one-fluid approximation breaks down and the character of the turbulence changes qualitatively. The resulting microscale plasma turbulence plays a role in laboratory experiments with strongly magnetized plasmas and it has been an important ingredient in theories of plasma confinement and magnetic reconnection (Biskamp et al. 1999; Cho & Lazarian 2004; Gürcan et al. 2009; Chandran et al. 2010; Kletzing et al. 2010). Recent high-resolution in situ observations of the solar wind also revealed the presence of significant magnetic, density, and electric fluctuations at subproton scales, often referred to as dispersion or dissipation scales, indicating a greater importance of subrange plasma dynamics in natural systems (Bale et al. 2005; Smith et al. 2006; Alexandrova et al. 2009, 2011; Kiyani et al. 2009; Sahraoui et al. 2009; Chen et al. 2010, 2012a, 2012b, 2013b; Salem et al. 2012; Šafránková et al. 2013; Podesta 2013; Kiyani et al. 2013). Studies of subproton turbulence shed light on some fundamental problems of plasma astrophysics, such as energy dissipation and particle heating by magnetic turbulence (e.g., Howes et al. 2008a; Schekochihin et al. 2009; Howes 2010; Chandran et al. 2010; Chandran 2010; Cranmer & van Ballegooijen 2012; Boldyrev & Perez 2012; Perri et al. 2012).

Microscale plasma turbulence has been studied to a much lesser extent compared with its Alfvénic counterpart. Both the qualitative nature of microscale turbulence, e.g., the plasma modes involved, the relative contributions of electrons and ions to the dynamics, and its quantitative characteristics (e.g., spectra of turbulence, degree of anisotropy, rate of energy dissipation) attract considerable interest.

Various approaches have been employed in analytical and numerical studies of plasma turbulence; these approaches have advantages or disadvantages depending on the specific goals of the study. While a full kinetic description in terms of the Vlasov–Maxwell equations is the most accurate approach, it is not practical in nonlinear numerical simulations. The gyrokinetic description provides a significant advantage over the full kinetic description in terms of computational cost, however, it is limited to frequencies below the ion-cyclotron frequency and it is not physically transparent. Fluid models are physically the most transparent and are simpler to implement in numerical simulations, however, their limits of validity are often not well established.

The two electromagnetic modes of particular interest at subproton scales are the whistler mode and the kinetic-Alfvén mode. These modes correspond to disparate regions in frequency space and they are often studied within different frameworks. In particular, the higher-frequency electron whistler modes have been discussed in the framework of the cold plasma model (e.g., Aleksandrov et al. 1984; Stix 1992) and, in the nonlinear case, in terms of electron MHD, two-fluid and Hall MHD, and shell models (e.g., Gordeev et al. 1994; Ghosh et al. 1996; Biskamp et al. 1999; Cho & Lazarian 2004; Galtier & Buchlin 2007; Shaikh & Shukla 2009; Saito et al. 2010). The lower-frequency kinetic-Alfvén modes have been discussed in kinetic, gyrokinetic, and fluid frameworks (e.g., Hollweg 1999; Howes et al. 2006; Gary & Smith 2009; Schekochihin et al. 2009; Sahraoui et al. 2012).

The subproton whistler and kinetic-Alfvén modes have characteristics in common as well as essential physical differences, which should be taken into account in phenomenological models and in the interpretation of observational data. The goal of the present work is to provide an analysis of subproton plasma turbulence, specifically concentrating on the role of kinetic-Alfvén and whistler modes. In both linear and nonlinear cases, we employ unifying analytic approaches that allow us to treat kinetic-Alfvén and whistler waves on the same footing and to effectively compare their dynamics. As subproton electromagnetic fluctuations are rarely covered in textbooks, we pay special attention to methodical aspects of the work: where possible, derivations are carried out rigorously and in detail and necessary assumptions and simplifications are clearly explained. We start with the derivation of the properties of the linear waves and then proceed with the derivation of the corresponding nonlinear equations and the analysis of the results of the numerical simulations. Our work complements previous fluid dynamics studies, gyrokinetic simulations, and known numerical solutions of the linearized Vlasov–Maxwell equations.

In the first part of this work (Sections 2–5), we present a rigorous derivation of subproton whistler and kinetic-Alfvén waves starting from the collisionless kinetic equation. We emphasize similarities and differences between these modes by establishing the sectors of the phase space where they exist, the role of the electrons and the ions in their dynamics and in their dissipation, and by analyzing the associated density fluctuations and electron and magnetic compressibilities. In the second part of the paper (Sections 6–9), we extend the analysis to the nonlinear kinetic-Alfvén and whistler dynamics. Starting from the picture of electron drift motion, in Sections 6 and 7 we derive in the equations for nonlinear kinetic-Alfvén waves in the lower-frequency region. In the higher-frequency region, we derive the equations for the whistler waves, that is, the electron MHD (EMHD) equations. We demonstrate that in the limit of strong background magnetic fields and oblique propagation, the EMHD equations simplify to the reduced EMHD (REMHD) equations; we compare and contrast the REMHD equations with the kinetic-Alfvén equations. In Sections 8 and 9, we discuss the phenomenological models of kinetic-Alfvén and whistler turbulence and present the results of numerical simulations. Finally, in Section 10, we summarize our results and discuss their relation to some other recent studies of subproton plasma turbulence.

In this section, we define our notation and outline the method that we use to analyze electromagnetic subproton plasma fluctuations. We start with the collisionless kinetic equation for the distribution function of particles of sort α:

$$\begin{eqnarray} &&\frac{\partial }{\partial t}f_{\alpha }+{\bf v}\frac{\partial }{\partial {\bf r}}f_{\alpha }+q_{\alpha }\left({\bf E} + \frac{1}{c}\left[{\bf v}\times {\bf B} \right]\right) \frac{\partial }{\partial {\bf p}}f_{\alpha }=0. \end{eqnarray} \tag{ 1 }$$

We then Fourier transform this equation in time and space and assume that the system has a uniform magnetic field B₀ applied in the z-direction. We furthermore assume that the wave vector k is in the x−z plane, k = (k_⊥, 0, k_z). For small perturbations, the particle distribution function can be represented as f_α = f_0α + δf_α, where f_0α is a Maxwellian and δf_0α is a small perturbation. The kinetic equation (1) can then be solved for the perturbation, which in cylindrical coordinates, v = (v_⊥cos ϕ, v_⊥sin ϕ, v_z), has the form (e.g., Aleksandrov et al. 1984; Stix 1992)

$$\begin{eqnarray} \delta f_{\alpha }&=&\frac{iq_{\alpha }}{T_{\alpha }} f_{0\alpha } \sum \limits _{n=-\infty }^{+\infty } \frac{\exp \left\lbrace -in\phi +i\frac{k_\perp v_\perp }{\Omega _{\alpha }}\sin \phi \right\rbrace }{\omega -k_z v_z-n\Omega _{\alpha }} \nonumber \\ &&\times \left[ \frac{n\Omega _{\alpha }}{k_\perp } J_n \left(\frac{k_\perp v_\perp }{\Omega _{\alpha }}\right)E_x + iv_\perp J^{\prime }_n \left(\frac{k_\perp v_\perp }{\Omega _{\alpha }}\right)E_y \right. \nonumber \\ &&+ \left. v_z J_n\left(\frac{k_\perp v_\perp }{\Omega _{\alpha }} \right)E_z \right], \end{eqnarray} \tag{ 2 }$$

where J_n is a Bessel function of the order of n. In this expression, Ω_α = q_αB₀/(m_αc) is the gyrofrequency of the particles of sort α, which can be either positive or negative depending on the sign of the charge. Substituting Equation (2) into the Maxwell equations, one derives the dielectric tensor _ij(ω, k) in a plasma with a uniform magnetic field, whose components are given by (Aleksandrov et al. 1984)

$$\begin{eqnarray} \ \, \epsilon _{xx} &=& 1- \sum \limits _{\alpha, n} \frac{n^2 \omega _{p\alpha }^2}{\omega (\omega -n\Omega _\alpha)} \frac{A_n(z_\alpha)}{z_\alpha } J_+(\beta _{n\alpha }),\qquad\qquad \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \epsilon _{yy}&=&\epsilon _{xx}+2 \sum \limits _{\alpha, n} \frac{\omega _{p\alpha }^2 z_\alpha }{\omega (\omega - n \Omega _\alpha)}A^\prime _n (z_\alpha)J_+(\beta _{n\alpha }),\qquad\quad\ \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \epsilon _{xy}&=&-\epsilon _{yx}=-i \sum \limits _{\alpha, n} \frac{n\omega _{p\alpha }^2}{\omega (\omega - n \Omega _\alpha)}A^\prime _n(z_\alpha) J_+(\beta _{n\alpha }),\quad\ \ \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \epsilon _{xz}&=& \epsilon _{zx} = \sum \limits _{\alpha, n} \frac{n \omega _{p\alpha }^2 k_\perp }{\omega \Omega _\alpha k_z} \frac{A_n(z_\alpha)}{z_\alpha } \left[1-J_+(\beta _{n\alpha }) \right],\qquad\ \ \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \epsilon _{yz}&=& -\epsilon _{zy}=-i \sum \limits _{\alpha, n} \frac{\omega _{p\alpha }^2 k_\perp }{\omega \Omega _\alpha k_z} A^\prime _n(z_\alpha)\left[1-J_+(\beta _{n\alpha })\right],\quad \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \epsilon _{zz}&=&1+ \sum \limits _{\alpha, n}\frac{\omega _{p\alpha }^2 (\omega -n\Omega _\alpha)}{\omega k_z^2 v_{T\alpha }^2} A_n(z_\alpha)\left[1-J_+(\beta _{n\alpha }) \right].\quad\ \ \end{eqnarray} \tag{ 8 }$$

In these expressions, we used the short-hand notation

$$\begin{eqnarray} &&\ \ \ \sum \limits _{\alpha, n}\equiv \sum \limits _{\alpha }\sum \limits _{n=-\infty }^{+\infty } \end{eqnarray} \tag{ 9 }$$

and we denoted

$$\begin{eqnarray} z_\alpha &=&\frac{k_\perp ^2 v_{T\alpha }^2}{\Omega _\alpha ^2}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \ \ \ \beta _{n\alpha }&=&\frac{\omega -n\Omega _\alpha }{|k_z|v_{T\alpha }}, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \quad\, A_n&=&e^{-z_\alpha }I_n(z_\alpha), \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \qquad\quad\ \ J_+(x)&=& x e^{-x^2/2}\int ^x_{i\infty }e^{\xi ^2/2}\, d\xi, \end{eqnarray} \tag{ 13 }$$

with the limiting forms

$$\begin{eqnarray} J_+(x) &\approx & 1+\frac{1}{x^2}+\cdots \ -i\sqrt{\frac{\pi }{2}}x e^{-{x^2}/{2}},\quad |x|\gg 1, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} J_+(x)&\approx & -i\sqrt{\frac{\pi }{2}}x, \quad |x|\ll 1. \qquad\qquad\qquad \end{eqnarray} \tag{ 15 }$$

In what follows, we need to know the limiting forms of the few first functions A_n for z_α ≪ 1, which are, to first order

$$\begin{eqnarray} &{ }& A_0\approx 1,\, A_0^{\prime }\approx -1, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &{ }& A_1= A_{-1} \approx z_\alpha /2. \end{eqnarray} \tag{ 17 }$$

The following notation is adopted here and throughout the paper: $\omega _{p\alpha }=\sqrt{4\pi n_{0\alpha }q^2_{\alpha }/m_\alpha }$ is the plasma frequency and $v_{T\alpha }=\sqrt{T_\alpha /m_\alpha }$ is the thermal velocity associated with the particles of sort α. For a plasma consisting of electrons and ions, the so-called ion-acoustic velocity can be defined $v_{s}=\sqrt{T_e/m_i}$. It is also convenient to introduce the Alfvén speed $v_A=B_0/\sqrt{4 \pi n_0 m_i}$ and the plasma beta, which is the ratio of the thermal energy of the particles to the magnetic energy of the plasma. The plasma beta can be different for ions and electrons if the plasma is non-isothermal, $\beta _i=2v_{Ti}^2/v_A^2$, $\beta _e=2v_{s}^2/v_{A}^2$. In the above expressions, I_n is a modified Bessel function of the order of n and the function J₊ is related to the plasma dispersion function (it should not be confused with the Bessel function).

Our discussion is mostly motivated by astrophysical plasmas (e.g., the solar wind at 1 AU, the warm interstellar medium), where the thermal energies of particles and the magnetic field energies can be on the same order, that is, the plasma beta is of the order of 1: β_α = 8πn_αT_α/B² ∼ 1. Under this assumption, the particle gyroscales and the inertial scales are on the same order, that is, $\rho _\alpha =v_{T\alpha }/\Omega _\alpha \sim c/\omega _{p_\alpha }$. Regarding the wavelengths under consideration, we are interested in subproton but superelectron scales, that is, 1/ρ_i ≪ k_⊥ ≪ 1/ρ_e. We also restrict ourselves to the frequencies below the electron cyclotron frequency, ω ≪ Ω_e.

Even within these assumptions, there are different regions that should be considered separately. To identify these regions, we first examine, in wavenumber space, the rays defined by k_z = θk_⊥, where θ is some constant. Then, for each such ray, we can divide the two-dimensional frequency–wavenumber space ω − k_⊥ into sectors that are separated by the lines ω = k_zv_Te, ω = k_⊥v_Ti, ω = Ω_i, and k_⊥ = 1/ρ_i. One can check that all the possibilities for subproton oscillations are covered by considering the regions numbered I–IV in Figures 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each of these regions, the wave modes can be analytically derived, if we stay sufficiently far from the boundaries. We will call those regions asymptotic regions. We will show that in asymptotic region I, the obtained dispersion relation corresponds to the kinetic-Alfvén wave. In this region, the ions are dynamically important; they are not spatially magnetized and they quickly adjust to the fluctuating electric potential. In Section 6, we demonstrate that the nonlinear dynamics in this regime can be captured by a two-field kinetic-Alfvén model, which can be used to describe the kinetic-Alfvén turbulence discussed in Section 8.

The solution in asymptotic region II corresponds to the whistler mode. The principal difference between region II and region I is the dynamics of the ions. While in the low-frequency region (region I), the ions move to adjust to the fluctuating electric potential, in the high-frequency region (region II), the ions are immobile and dynamically irrelevant. As a consequence of quasi-neutrality, the whistler waves are then not associated with plasma density fluctuations. The electron nature of waves in region II allows one to derive the nonlinear equations for the whistler waves, the so-called EMHD equations. In Section 7, we derive the EMHD equations from the electron drift picture, similar to the derivation of the kinetic-Alfvén model. We then demonstrate how the EMHD equations are simplified in the case of a strong background magnetic field and oblique wave propagation, thus leading to the REMHD model. We then compare and contrast this model with the kinetic-Alfvén model derived in Section 6. In Section 9, we discuss both weak and strong turbulence of whistlers. Finally, it can be shown that there are no solutions in regions III and IV.

In order to find the linear perturbations (waves) that can propagate in a plasma, one needs to solve the linear system of equations

$$\begin{eqnarray} &&D_{lm}\, E_m(\omega, {\bf k})=0, \end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray} &&D_{lm}\equiv k^2\delta _{lm} -k_lk_m -\frac{\omega ^2}{c^2}\epsilon _{lm}(\omega, {\bf k}). \end{eqnarray} \tag{ 19 }$$

E(ω, k) is the polarization of the electric field and summation over repeated indices is assumed. In a general case, the solution of this equation cannot be obtained analytically. However, the analytic solutions can be found for the asymptotic regions I–IV in Figures 1 and 2. These solutions are derived below.

Here, we consider the solution of Equation (18) in region I, that is, we apply the assumptions ω ≪ k_⊥v_Ti and ω ≪ k_zv_Te. The general expressions (3)–(8) can be simplified in this limit in the following way

$$\begin{eqnarray} \epsilon _{xx}&\simeq & 1+\frac{\omega _{pe}^2}{\Omega _e^2} + \frac{ \omega _{pi}^2}{k_\perp ^2 v_{Ti}^2} \simeq \frac{ \omega _{pi}^2}{k_\perp ^2 v_{Ti}^2},\quad\qquad\quad\ \ \ \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} \ \epsilon _{yy}&\simeq & 1+\frac{\omega _{pe}^2}{\Omega _e^2} \simeq \frac{\omega _{pe}^2}{\Omega _e^2},\qquad\qquad\qquad\qquad \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \epsilon _{xy}&= &-\epsilon _{yx}\simeq i\frac{\omega _{pe}^2}{\omega \Omega _e},\qquad\qquad\qquad\qquad \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \epsilon _{xz}&= &\epsilon _{zx}= O\left(k_z\rho _i \frac{\omega _{pi}^2}{k_\perp ^2 v_{Ti}^2}\right),\qquad\qquad\quad\ \ \ \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \epsilon _{yz}&= &-\epsilon _{zy} \simeq i\frac{\omega _{pe}^2 k_\perp }{\omega \Omega _e k_z},\qquad\qquad\qquad\quad\ \ \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} \epsilon _{zz}&\simeq & 1+\frac{\omega _{pe}^2}{k_z^2 v_{Te}^2} \simeq \frac{\omega _{pe}^2}{k_z^2 v_{Te}^2}.\qquad\qquad\qquad\ \ \end{eqnarray} \tag{ 25 }$$

The derivation requires an explanation. While the electron contribution is derived straightforwardly from Equations (3)–(8) by keeping only the leading terms in these expressions, the ion contribution is not as easy to derive from Equations (3)–(8) since one needs to keep many similar terms in the series. For example, in the electron contribution to _xx in Equation (3), one needs to retain only the lowest-order terms n = ±1. In contrast, in the ion contribution one needs to keep the terms up to $n \sim \sqrt{z_i}\gg 1$, as it follows from the asymptotic form of the modified Bessel function I_n(z) for $1 \ll n \sim \sqrt{z}$

$$\begin{eqnarray} &&I_n(z)\sim \frac{e^{z}e^{-n^2/(2z)}}{\sqrt{2\pi z}}. \end{eqnarray} \tag{ 26 }$$

See, e.g., Abramowitz & Stegun (1972, Chapter 9) and Olver (1997). The summation over n in Equation (3) can be done with the aid of the result

$$\begin{eqnarray} &&\sum \limits _{n=-\infty }^{+\infty }I_n(z)=e^z, \end{eqnarray} \tag{ 27 }$$

which leads to

$$\begin{eqnarray} &&\epsilon ^{(i)}_{xx}=\frac{\omega _{pi}^2}{k_\perp ^2 v_{Ti}^2}. \end{eqnarray} \tag{ 28 }$$

In this derivation, we also assumed that k_z ≪ k_⊥, which can be verified a posteriori (cf. the discussion following Equation (40)). In a similar fashion, one can find the ion contributions to the other components of the dielectric tensor, where, together with Equation (27), one needs to use the formula

$$\begin{eqnarray} &&\sum \limits _{n=-\infty }^{+\infty }{n^2}{I_n(z)}= z e^z \end{eqnarray} \tag{ 29 }$$

and the asymptotic result

$$\begin{eqnarray} &&\sum \limits _{n=1}^{+\infty }\frac{I_n(z)}{n^2}\rightarrow \frac{\pi ^2}{6}\frac{e^z}{\sqrt{2\pi z}}, \quad z\rightarrow \infty. \end{eqnarray} \tag{ 30 }$$

The derivation of the ion contribution from the series in Equations (3)–(8) is however not physically transparent and it is not very convenient. A more efficient way is to note that in the limit k_z ≪ k_⊥, the dominant ion contribution comes through the field-perpendicular components of the ion dielectric tensor (Equations (20)–(22)). For those components, the limit of large z_i corresponds to the limit of small magnetic field, in which case the ion contribution $\epsilon ^{(i)}_{lm}$ to the dielectric tensor $\epsilon _{lm}=\delta _{lm}+\epsilon ^{(e)}_{lm} +\epsilon ^{(i)}_{lm}$ is the same as in the non-magnetized case. Taking into account ω ≪ kv_Ti and k_z ≪ k_⊥, we have, to leading order

$$\begin{eqnarray} &&\epsilon ^{(i)}_{lm}=\frac{k_l k_m}{k_\perp ^2}\frac{\omega _{pi}^2}{k_\perp ^2 v_{Ti}^2}. \end{eqnarray} \tag{ 31 }$$

This expression for the ion contribution was used in the derivation of Equations (20)–(22). We can now write the expressions for the tensor D_lm in Equation (19)

$$\begin{eqnarray} D_{xx}&\simeq & k_z^2-\frac{\omega ^2}{c^2} \frac{\omega _{pi}^2}{k_\perp ^2 v_{Ti}^2}, \qquad\, \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} D_{yy}&\simeq & k^2 \simeq k_\perp ^2,\qquad\qquad\ \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} D_{xy}&=& -D_{yx}\simeq -i \frac{\omega \omega _{pe}^2}{\Omega _e c^2},\quad \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} \ \ D_{xz}&=& D_{zx}\simeq -k_\perp k_z,\qquad\quad \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} D_{yz}&=& -D_{zy}\simeq -i\frac{\omega _{pe}^2\omega k_\perp }{\Omega _e c^2 k_z}, \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} \ \ D_{zz}&\simeq & k_\perp ^2-\frac{\omega ^2 \omega _{pe}^2}{k_z^2 v_{Te}^2 c^2}.\qquad\quad \end{eqnarray} \tag{ 37 }$$

The dispersion equation (18) now allows us to write the equation for the wave frequency, $\det (D_{lm})\equiv A\omega ^4-B\omega ^2=0$, where

$$\begin{eqnarray} A&=&\frac{\omega _{pi}^2\omega _{pe}^2}{c^4v_{Te}^2v_{Ti}^2k_z^2}+\frac{\omega _{pe}^6}{\Omega _e^2c^6 v_{Te}^2k_z^2}+\frac{\omega _{pe}^4\omega _{pi}^2}{c^6 v_{Ti}^2\Omega _e^2 k_z^2}, \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} B&=& \frac{k_\perp ^2\omega _{pi}^2}{c^2 v_{Ti}^2}+\frac{k_\perp ^2\omega _{pe}^2}{c^2 v_{Te}^2},\qquad\qquad\qquad\quad\ \ \end{eqnarray} \tag{ 39 }$$

with the solution

$$\begin{eqnarray} &&\omega ^2 =\frac{\left[ \frac{\omega _{pi}^2}{v_{Ti}^2} +\frac{\omega _{pe}^2}{v_{Te}^2}\right]}{\left[ \frac{\omega _{pe}^2\omega _{pi}^2}{c^2 v_{Ti}^2 v_{Te}^2}+\frac{\omega _{pe}^6}{\Omega _e^2 c^4 v_{Te}^2}+\frac{\omega _{pi}^2 \omega _{pe}^4}{\Omega _e^2 c^4 v_{Ti}^2}\right]}k_z^2 k_\perp ^2. \end{eqnarray} \tag{ 40 }$$

In the above expressions, we have approximated k ≈ k_⊥; now, we can explain this assumption. By substituting the obtained frequency (Equation (40)) into the condition ω ≪ kv_Ti, we derive k_zρ_i ≪ 1, which, together with our main assumption k_⊥ρ_i ≫ 1, ensures that the wave propagation must be oblique. This dispersion relation can be re-written in a more familiar form by assuming singly charged ions (n_i = n_e) and by introducing the ion-acoustic scale ρ_s = v_s/Ω_i, where $v^2_s=T_e/m_i$ is the ion-acoustic speed

$$\begin{eqnarray} &&\omega ^2 = \frac{v_A^2\rho _s^2\left(1+\frac{T_i}{T_e}\right)}{1+\frac{v_s^2}{v_A^2}\left(1+\frac{T_i}{T_e}\right)}k_z^2k_\perp ^2. \end{eqnarray} \tag{ 41 }$$

This expression is the dispersion relation for the kinetic-Alfvén waves (e.g., Camargo et al. 1996; Terry et al. 2001; Howes et al. 2006; Schekochihin et al. 2009).

For completeness, we present the polarization of the electric field, which can be found from Equation (18)

$$\begin{eqnarray} &&E_z=-\frac{T_e}{T_i}\frac{k_z}{k_\perp } E_x, \qquad\qquad\qquad \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} &&E_y=i \sqrt{\frac{\left(1+\frac{T_e}{T_i}\right)\frac{v_{Ti}^2}{v_A^2}}{1+\frac{v_{Ti}^2}{v_A^2}\left(1+\frac{T_e}{T_i} \right)}} \left(1+\frac{T_e}{T_i} \right)\frac{|k_z|}{k_\perp }E_x. \end{eqnarray} \tag{ 43 }$$

For large angles of propagation, k_z ≪ k_⊥, the x-component of the electric field dominates over the other two components, implying that the electric field is almost potential (e.g., Howes et al. 2006; Schekochihin et al. 2009; Sahraoui et al. 2012). We note, however, that from Equation (42) the small E_z component is not given by the potential part of the electric field (in which case we would have E_z/k_z = E_x/k_⊥); rather, it has comparable contributions from both the potential and the solenoidal parts. Kinetic-Alfvén waves are essentially electromagnetic, not electrostatic, waves.

We also see that the oblique kinetic-Alfvén waves are right-hand polarized, that is, their electric field rotates around the $\hat{z}$ axis similarly to the electrons. This situation is in contrast with the kinetic-Alfvén waves propagating along the guide field, which are left-hand polarized. This phenomenon, observed numerically by Gary (1986), explains why oblique kinetic-Alfvén waves do not resonate with the ions and why they can exist both below and above the ion-cyclotron frequency.

Once the polarization of the electric field is obtained, the magnetic field polarization can be straightforwardly found from the induction equation, b_k = (c/ω)[k × E_k]. In particular, we obtain the following expression for the amplitude of the magnetic fluctuations in kinetic-Alfvén waves

$$\begin{eqnarray} &&b_k^2 = \frac{c^2 \Omega _i^2}{ v_A^2 v_{Ti}^2}\left[ 2\frac{v_{Ti}^2}{v_A^2} \left(1+\frac{T_e}{T_i} \right)^2 + \left(1+\frac{T_e}{T_i} \right) \right]\frac{E_x^2}{k_\perp ^2}. \quad \end{eqnarray} \tag{ 44 }$$

The result (Equations (42) and (43)) also allows one to analyze the electron and ion-density perturbations in the kinetic-Alfvén wave. For that, we integrate δf_e in Equation (2) over velocity and obtain

$$\begin{eqnarray} &&\delta n_e= i\frac{n_0 e}{T_e k_z} E_z, \end{eqnarray} \tag{ 45 }$$

where e is the modulus of the electron charge. We see that it is the z-component of the electric field that couples to the electron density, reflecting the fact that the electron motion is non-magnetized along the direction of the magnetic field.

Further insight into the physics of density fluctuations is gained if we use Equation (42) to express the electron-density perturbations through the field-perpendicular component of the electric field, E_x, which is almost potential

$$\begin{eqnarray} &&\delta n_e=-i\frac{n_0 e}{T_i k_\perp } E_x=-\frac{n_0 q_i}{T_i}\phi _k. \end{eqnarray} \tag{ 46 }$$

Noting that the electron and ion perturbations are the same due to quasi-neutrality, δn_i = δn_e, we observe that the ion-density perturbations follow the electric potential. This fact is not surprising if one notes that, as we discussed above, the ions are spatially non-magnetized and, therefore, their low-frequency response can be obtained from the non-magnetic collisionless kinetic equation, which leads to the result

$$\begin{eqnarray} &&\frac{\delta n_i}{n_0}=-\frac{q_i \phi _k}{T_i}\left[ 1-J_+\left(\frac{\omega }{kv_{Ti}} \right) \right]\approx -\frac{q_i\phi _k}{T_i}, \end{eqnarray} \tag{ 47 }$$

where ϕ_k = ik · E/k² describes the potential part of the electric field and the asymptotic form of the function J₊ is given in Equation (15). This expression is the same as Equation (46).

It may seem puzzling that the ions and the electrons follow different potentials, ϕ_k and ${\phi _k^{\prime }}=iE_z/k_z$, respectively. The explanation is that the electrons are moving with a velocity v_e that can be found from the integral n_ev_e = ∫vδf_ed³v and the magnetic field lines are advected by the same velocity. The electrons thus follow the potential ϕ' existing in the moving frame, to which they quickly adjust along the magnetic field lines. The ions, on the other hand, are not spatially magnetized and they adjust to the "laboratory-frame" potential ϕ, which is different from ϕ', since the electric field is different in different frames.

This effect reveals important differences in the electron and ion motion, which we now study in more detail. The short-cut method for calculating the electron velocity is to use the electron conductivity tensor, which is expressed through the electron part of the dielectric tensor

$$\begin{eqnarray} \sigma ^{(e)}_{lm}&=&\frac{\omega }{4\pi i}\epsilon _{lm}^{(e)} \equiv -\frac{c^2}{4 \pi i\omega }D_{lm} -\frac{\omega }{4\pi i}\epsilon _{lm}^{(i)} \nonumber \\ &&+ \frac{c^2}{4 \pi i \omega }\left(k^2\delta _{lm}-k_lk_m-\frac{\omega ^2}{c^2}\delta _{lm}\right), \end{eqnarray} \tag{ 48 }$$

where D_lm is given by Equation (18) and $\epsilon ^{(i)}_{lm}$ is given by Equation (31). The electron velocity is related to the electron current density,

$$\begin{eqnarray} &&{J}_{e, l}=-en_0v_{e, l}=\sigma ^{(e)}_{lm} E_m, \end{eqnarray} \tag{ 49 }$$

which is easy to calculate from Equation (48) since the wave electric field satisfies D_lmE_m = 0. As a result, we obtain the following equations for the electron velocity field

$$\begin{eqnarray} v_{e, x}&=& c\frac{E_y}{B_0},\qquad\qquad\ \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} v_{e, y}&=&-c\frac{E_x}{B_0}-c\frac{T_e}{T_i}\frac{E_x}{B_0}, \end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray} v_{e, z}&=&-\frac{ic}{4 \pi n_0 e}k_\perp b_y. \ \ \ \end{eqnarray} \tag{ 52 }$$

This result has the following physical interpretation. The electron-fluid velocity across the magnetic field is just the "E × B" drift plus the diamagnetic drift due to the gradient of the (isothermal) electron pressure, p_e = n_eT_e. This fact implies that the magnetic field lines are "frozen" into the electron fluid. The electron velocity along the magnetic field lines is related to the total parallel current J_z = (ic/4π)[k × b]_z, since, as can be verified from Equation (48), the ion contribution to the parallel current is negligible.

The ion velocity field can be similarly calculated from the ion part of the conductivity tensor,

$$\begin{eqnarray} &&J_{i, l}=\frac{\omega }{4\pi i}\epsilon _{lm}^{(i)}E_m=\frac{\omega \omega _{pi}^2}{4\pi i k_\perp ^4 v_{Ti}^2}({\bf k} \cdot {\bf E})k_l, \end{eqnarray} \tag{ 53 }$$

which implies that the ions, being spatially non-magnetized, adjust to the electric potential. It can be demonstrated that Equation (53) can also be obtained from the continuity equation for the ions, where the ion density is expressed through the electric potential according to Equation (47). The ions have significant velocity (comparable to that of the electrons) only in the x-direction, while their velocities in the y- and z-directions can be neglected compared with the corresponding electron velocities.

The obtained expressions also demonstrate that in the electron co-moving frame, the electric field becomes purely potential with the potential ϕ'. The electron and ion behavior that we observed in the linear case will be preserved in the case of strong kinetic-Alfvén turbulence. We will return to the fluid-like picture when we discuss the nonlinear equations for the kinetic-Alfvén waves in Section 6.

For practical applications, it may be useful to know the so-called electron compressibility, that is, the ratio of the normalized electron-density fluctuations to the normalized magnetic fluctuations (e.g., Gary & Smith 2009). From Equations (44), (46), and (42), we obtain

$$\begin{eqnarray} &&C_e^{ka}\equiv \frac{(\delta n_e/n_0)^2}{(b/B_0)^2} =\frac{1}{{ 2\frac{v_{Ti}^4}{v_A^4} \left(1+\frac{T_e}{T_i} \right)^2 + \frac{v^2_{Ti}}{v^2_A}\left(1+\frac{T_e}{T_i} \right) } }.\quad \end{eqnarray} \tag{ 54 }$$

We note that this ratio is independent of wavenumber and that it strongly depends on the plasma beta $\beta _i=2v_{Ti}^2/v_A^2$. In a similar fashion, one can obtain the so-called magnetic compressibility, that is, the ratio of the z-component of the magnetic fluctuations to the total amplitude of the fluctuations (e.g., Gary & Smith 2009). By evaluating b_z from the induction equation and using Equation (44), we obtain

$$\begin{eqnarray} &&C_\Vert ^{ka}\equiv \frac{b_z^2}{b^2}=\frac{\frac{v_{Ti}^2}{v_A^2}\left(1+\frac{T_e}{T_i} \right)}{1+2\frac{v_{Ti}^2}{v_A^2}\left(1+\frac{T_e}{T_i} \right)}. \end{eqnarray} \tag{ 55 }$$

The magnetic compressibility is also independent of the wavenumber and it strongly depends on the plasma beta.

We point out that the ions are dynamically important in kinetic-Alfvén waves, which distinguishes the physics of the kinetic-Alfvén modes from the physics of the whistler modes where the ions are immobile and dynamically irrelevant. The whistler mode is discussed in the next section.

We now consider the asymptotic region II defined by kv_Ti ≪ ω ≪ k_zv_Te. In this case, the electron contribution to the dielectric function is the same as in region I and it can be inferred from Equations (20)–(25). The ion contribution is however different, since the ions can now be treated as "cold." The ion contribution is given by

$$\begin{eqnarray} &&\ \ \ \ \epsilon _{lm}^{(i)}=-\delta _{lm}\frac{\omega _{pi}^2}{\omega ^2}. \qquad\ \end{eqnarray} \tag{ 56 }$$

Keeping the dominant terms in the tensor D_lm, we obtain

$$\begin{eqnarray} \ D_{xx}&\simeq & k_z^2+\frac{\omega _{pi}^2}{c^2},\qquad \end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray} \ D_{yy}&\simeq & k^2,\qquad\qquad\ \end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray} \qquad D_{xy}&=& -D_{yx}\simeq -i \frac{\omega \omega _{pe}^2}{\Omega _e c^2},\ \ \end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray} D_{xz}&=& D_{zx}\simeq -k_\perp k_z,\ \end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray} \qquad\ \ D_{yz}&=& -D_{zy}\simeq -i\frac{\omega _{pe}^2\omega k_\perp }{\Omega _e c^2 k_z}, \end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray} D_{zz}&\simeq & k_\perp ^2-\frac{\omega ^2 \omega _{pe}^2}{k_z^2 v_{Te}^2 c^2}. \end{eqnarray} \tag{ 62 }$$

It is worth noting here that the only difference between the wave equations in the whistler and kinetic-Alfvén cases is the ion contributions in the D_xx terms in Equations (32) and (57).

The equation for the wave frequency follows from Equation (18) and it has a biquadratic form: $\det (D_{lm})\equiv A \omega ^4-B \omega ^2+C=0$, where

$$\begin{eqnarray} A&=&\frac{\omega _{pe}^6}{\Omega _e^2c^6k_z^2v_{Te}^2},\qquad\qquad\qquad \end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray} \ \ B&=&\frac{\omega _{pe}^2k^2}{v_{Te}^2c^2}+\frac{k^2 \omega _{pi}^2 \omega _{pe}^2}{k_z^2 v_{Te}^2 c^4}+\frac{k_\perp ^2 \omega _{pe}^4 \omega _{pi}^2}{k_z^2 \Omega _e^2 c^6}, \end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray} C&=&\frac{\omega _{pi}^2k_\perp ^2 k^2}{c^2}.\qquad\qquad\qquad\ \ \ \ \end{eqnarray} \tag{ 65 }$$

It can be checked that only one solution of this equation can be consistent with the assumption kv_Ti ≪ ω. For that, one has to require that k_zρ_i ≫ 1, which implies that for the considered plasma parameters, the whistlers cannot be too oblique. In this case, the first term in Equation (64) dominates over the other two terms, C can be neglected, and a solution is found:

$$\begin{eqnarray} &&\omega ^2=\frac{\Omega _e^2 c^4}{\omega _{pe}^4}k_z^2 k^2 = \frac{v_A^4 \rho _s^2 }{v_s^2}k_z^2 k^2. \end{eqnarray} \tag{ 66 }$$

This relation is the dispersion relation for the subproton whistler waves.³$^,$⁴

The polarization of the electric field can be found from Equation (18) in this approximation

$$\begin{eqnarray} \ \ \ E_y&=&i\frac{|k_z|}{k}E_x, \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} E_z&=&0. \quad\ \ \ \end{eqnarray} \tag{ 68 }$$

This result, together with the fact that the electron-density fluctuations couple to the E_z-component of the electric field (see Equation (45)), implies that density fluctuations are negligible in whistler waves. This fact can be demonstrated directly from the non-magnetized collisionless kinetic equation for the ions (due to quasi-neutrality, n_e = n_i),

$$\begin{eqnarray} &&\frac{\delta n_e}{n_0}=\frac{\delta n_i}{n_0}=-\frac{q_i \phi _k}{T_i}\left[ 1-J_+\left(\frac{\omega }{kv_{Ti}} \right) \right], \end{eqnarray} \tag{ 69 }$$

where ϕ_k = ik · E/k² describes the potential part of the electric field. For high-frequency fluctuations, kv_Ti/ω ≪ 1, the asymptotics of the J₊ function are given by Equation (14), which shows that the density variations are indeed negligibly small,

$$\begin{eqnarray} &&\left| \frac{\delta n_e}{n_0} \right| =\frac{q_i \phi _k}{T_i}\frac{k^2 v_{Ti}^2}{\omega ^2} \ll \frac{q_i \phi _k}{T_i}. \end{eqnarray} \tag{ 70 }$$

In particular, one can compare the density fluctuations to the magnetic fluctuations by calculating the electron compressibility, similarly to our discussion in Equation (54). From the magnetic induction equation and the electric field polarization (Equations (67) and (68)), we obtain the magnetic field fluctuations in the whistler wave,

$$\begin{eqnarray} &&b^2=\frac{2c^2}{\omega ^2}k_z^2 E_x^2, \end{eqnarray} \tag{ 71 }$$

and the electron compressibility for the whistler wave is found to be

$$\begin{eqnarray} &&C_e^w\equiv \frac{(\delta n_e/n_0)^2}{(b/B_0)^2}=\frac{1}{2}\frac{\omega _{pi}^4}{k_z^4 c^4}. \end{eqnarray} \tag{ 72 }$$

Therefore, the whistler electron compressibility is small, as whistlers exist only for k_z > 1/ρ_i ∼ ω_pi/c. Our analytic result explains the behavior previously observed numerically by Gary & Smith (2009); moreover, it turns out that Equation (72) formally derived in the asymptotic limit (ω_pi/(k_zc) ≪ 1) holds quite well already for ω_pi/(k_zc) ≲ 1.

In addition, the properties of the whistler waves can be characterized by the magnetic compressibility, which is defined as

$$\begin{eqnarray} &&C_\Vert ^w\equiv \frac{b_z^2}{b^2}=\frac{k_\perp ^2}{2k^2}. \end{eqnarray} \tag{ 73 }$$

The whistler magnetic compressibility is independent of the plasma beta, and it is practically independent of wavenumber for k_z ≪ k_⊥.

The obtained results reveal essential differences between the kinetic-Alfvén modes and the whistler modes. The kinetic-Alfvén modes are compressible modes where the ion dynamics play an essential role. For the considered plasma parameters, they exist at the subproton scales in the frequency range ω ≪ kv_Ti, k_zv_Te (depicted as region I in Figures 1 and 2) and in the wavenumber range k_zρ_i ≪ 1. In contrast, the whistler waves are electron waves, where the ion dynamics are not essential. The ions provide a uniform stationary background, which, together with the quasi-neutrality condition, explains the incompressible character of whistler waves. We note that it is the high frequency of the whistler waves, not the large plasma beta, that ensures the immobility of the ions and the resulting incompressibility of the whistler modes. For the considered plasma parameters, the whistler waves exist at subproton scales in the frequency range kv_Ti ≪ ω ≪ k_zv_Te (depicted as region II in Figure 1) and in the wavenumber range k_zρ_i ≫ 1.

Therefore, the kinetic-Alfvén modes and the whistler modes correspond to different frequency and wavenumber ranges and they are governed by different physical mechanisms. This fact is reflected in, e.g., the different normalizations of the dispersion relations, (cf. Equations (41) and (66)) and the different variation of electron and magnetic compressibilities with plasma parameters (cf. Equations (54) and (55) with Equations (72) and (73)). Finally, it can be checked that under the considered plasma parameters, there are no wave modes in regions III and IV.

So far, we considered the subproton plasma fluctuations under the assumption that dissipation can be neglected. In this section, we evaluate the damping of subproton fluctuations due to weak collisionless dissipation. The terms responsible for the dissipation in the dielectric tensor (Equations (3)–(8)) arise from the small imaginary parts of the function J₊ in Equations (14) and (15), which should now be retained together with the real parts. As a result, the components of the dielectric tensor will acquire small corrections, $\epsilon _{lm}= \epsilon _{lm}^0+\delta \epsilon _{lm}$, and the corresponding small corrections will appear in the tensor $D_{lm}= D_{lm}^0+\delta D_{lm}$, where $\delta D_{lm}=-\frac{\omega ^2}{c^2}\delta \epsilon _{lm}$ and the superscript 0 denotes the components where the dissipation parts are neglected.

The solution of the new dispersion equation $\det (D_{lm})=0$ will now take the form ω = ω₀ + iγ, where ω₀ is the frequency of a non-dissipative wave, and γ is the small dissipation rate. Owing to the smallness of δ_lm and γ, the dispersion equation can be re-written in a simplified form. Let us denote $D=\det (D_{lm})$ and M_lm as the cofactor corresponding to the matrix element D_lm. We then have

$$\begin{eqnarray} &&\delta D= M_{lm}^0 \delta D_{lm}(\omega _0)+\frac{\partial D^0}{\partial \omega _0}i\gamma =0, \end{eqnarray} \tag{ 74 }$$

where we used the identity ∂D/∂D_lm = M_lm and summation over repeated indices is assumed. We therefore derive the general expression for the dissipation rate

$$\begin{eqnarray} &&\gamma =i\frac{M_{lm}^0 \delta D_{lm}(\omega _0)}{(\partial D^0/\partial \omega _0)}. \end{eqnarray} \tag{ 75 }$$

Further simplification comes from the observation that in both the kinetic-Alfvén case and the whistler case, the determinant has a simple general form,

$$\begin{eqnarray} &&\det \left(D^0_{lm}\right)=A\omega _0^4-B\omega _0^2, \end{eqnarray} \tag{ 76 }$$

which leads to ∂D⁰/∂ω₀ = 2Bω₀. The coefficient B for the kinetic-Alfvén case is found in Equation (39); for the whistler case, B is found in Equation (64).

We now consider the dissipation separately for the kinetic-Alfvén waves and the whistler waves. In the case of β_i ∼ 1, the dissipation of oblique kinetic-Alfvén waves is dominated by the electron contribution, see, e.g., the numerical solutions of the Vlasov–Maxwell equations in Howes et al. (2008a). We will therefore consider only electron dissipation in this case. In contrast, in the case of the whistler waves, ion dissipation becomes significant and it may quickly overcome the electron dissipation as the obliquity increases. In the whistler case, we therefore derive both electron and ion contributions to the dissipation.

In the kinetic-Alfvén case, the dominant dissipation terms in the tensor D_lm are

$$\begin{eqnarray} \ \delta D^{(e)}_{yy}&=&-i\sqrt{2\pi } \,\frac{\omega \omega _{pe}^2k_\perp ^2 v_{Te}}{c^2\Omega _e^2|k_z|}, \end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray} \qquad\quad\ \delta D^{(e)}_{yz}&=&-\delta D_{zy}^{(e)}=\sqrt{\frac{\pi }{2}} \frac{\omega ^2 \omega _{pe}^2 k_\perp }{\Omega _e v_{Te}c^2 |k_z|k_z}, \end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray} \delta D^{(e)}_{zz}&=&-i\sqrt{\frac{\pi }{2}}\frac{\omega ^3\omega ^3_{pe}}{v_{Te}^3c^2|k_z|k_z^2}. \ \end{eqnarray} \tag{ 79 }$$

A straightforward but somewhat lengthy calculation using Equation (75) then gives

$$\begin{eqnarray} &&\frac{\gamma }{\omega _0}=-F_{{\rm kaw}}^{(e)}k_\perp \rho _e, \qquad\ \end{eqnarray} \tag{ 80 }$$

where

$$\begin{eqnarray} &&\qquad\ F_{{\rm kaw}}^{(e)}=\sqrt{\frac{\pi }{2}}\sqrt{\frac{{\beta _i}}{{(1+a)(2+\beta _i(1+a))}}} \nonumber \\ &&\times \left[1+a+\frac{2}{\beta _i(2+\beta _i(1+a))} \right], \end{eqnarray} \tag{ 81 }$$

where we denoted a = T_e/T_i. Formulae (80) and (81) describe the collisionless damping of the subproton kinetic-Alfvén waves; they agree with previous findings (e.g., Howes et al. 2006, 2008b). For the case of β_i = 1 and a = 1, we estimate $F_{{\rm kaw}}^{(e)}\approx 1.1$.

We now turn to the discussion of collisionless damping of whistler waves. The dissipative contributions to the D_lm tensor will now be given by both the electrons and the ions. The electron contribution is given by the same expressions (77)–(79) as in the kinetic-Alfvén case. The ion contribution is obtained from the dielectric tensor for non-magnetized ions. Taking into account that whistlers exist in the region ω ≫ kv_Ti, we obtain

$$\begin{eqnarray} \delta D^{(i)}_{xx}&=&-i\left[1+\frac{k_\perp ^2}{k^2} \right]G, \end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray} \delta D^{(i)}_{yy}&=&-iG,\qquad\quad\ \ \end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray} \quad \delta D^{(i)}_{xz}&=&\delta D^{(i)}_{zx}= -i\frac{k_zk_\perp }{k^2}G,\, \end{eqnarray} \tag{ 84 }$$

$$\begin{eqnarray} \delta D^{(i)}_{zz}&=& -i\left[1+\frac{k_z^2}{k^2} \right]G, \end{eqnarray} \tag{ 85 }$$

where we have denoted

$$\begin{eqnarray} &&G=\sqrt{\frac{\pi }{2}}\frac{\omega ^3}{kv_{Ti}c^2}\frac{\omega _{pi}^3}{k^2 v_{Ti}^2}\exp {\left(-\frac{\omega ^2}{2k^2 v_{Ti}^2}\right)}. \end{eqnarray} \tag{ 86 }$$

A straightforward computation taking into account both electron and ion contributions leads to the dissipation rate of the subproton whistler waves:

$$\begin{eqnarray} &&\frac{\gamma }{\omega _0}=-\sqrt{\frac{\pi }{2}}\frac{k_\perp }{k}k_\perp \rho _e - \frac{2\sqrt{\pi }|k_z| c}{{\beta _i}^{3/2}\omega _{pi}}\exp \left(-\frac{k_z^2 c^2}{\beta _i\omega _{pi}^2}\right).\quad \end{eqnarray} \tag{ 87 }$$

We see that the electron dissipation of the whistler waves (the first term in Equation (87)) is comparable to that of the kinetic-Alfvén waves, while the ion contribution to the dissipation (the second term in Equation (87)) may become significantly larger when the propagation is oblique. For example, for β_i = 1, the ion dissipation rate is about 1.6% for k_zc/ω_pi = 2.5 and it reaches 14% for k_zc/ω_pi = 2. Thus, the ion dissipation rate increases quickly as k_zc/ω_pi decreases so that the whistlers become strongly damped by the ions at k_zc/ω_pi ≲ 2.

From Equation (87), it follows that the electron damping of whistlers limits k_⊥ from above while the ion damping limits k_z from below, imposing restrictions on the obliquity of propagation. For a given k_⊥, as k_z decreases and the propagation becomes more oblique, the dissipation becomes stronger and the whistler wave ceases to exist at k_z ∼ ω_pi/c. It is easy to see that this boundary corresponds to the line ω = k_⊥v_Ti in Figure 1, separating the whistler region II from the kinetic-Alfvén region I in phase space. For even more oblique propagation (below the boundary ω = k_⊥v_Ti), a new mode appears in region I, the kinetic-Alfvén wave, whose dispersion relation is different from that of a whistler wave.

The discontinuity between the dispersion relations of the kinetic-Alfvén and whistler waves, Equations (41) and (66), reflects the differences in the physics governing the two modes. These differences must be taken into account in order to derive the nonlinear equations for the whistler and kinetic-Alfvén waves. This analysis is presented in the following sections. We demonstrate that the whistler waves are described by the EMHD equations or, in the case of oblique propagation, by the REMHD equations. In contrast, the kinetic-Alfvén waves are described by the kinetic-Alfvén system of equations, which is different from the EMHD equations.

As we discussed in Section 3, kinetic-Alfvén waves exist in the wavenumber region k_z ≪ k_⊥ for the plasma parameters we consider, where k_z and k_⊥ are the typical wavenumbers of fluctuations in the directions parallel and perpendicular to the magnetic field, respectively. As will be explained in more detail below, the nonlinear effects in this case become significant for relatively small fluctuations of the magnetic field, satisfying k_zB₀ ∼ k_⊥b, where B₀ is the non-fluctuating component and b is the fluctuating part of the magnetic field. This fact implies that even in the case of strong turbulence, fluctuations of the distribution function remain small and they can be obtained through perturbation theory.

To demonstrate this finding, we consider the electron kinetic equation (Equation (1)) with a small electric field. This equation has a solution of the form f_e(x, v, t) ≃ f_0e(v) + δf_e, where f_0e is Maxwellian⁵ and the small perturbation δf_e is given by Equation (2). From this expression, we estimate δf_e ∼ (eϕ'/T)f_0e, where we denote ϕ' = iE_z/k_z. The applicability of the perturbation procedure thus requires that eϕ'/T ≪ 1. This condition does not necessarily imply that the corresponding turbulence is weak, however. Indeed, the condition for strong turbulence, k_⊥b ∼ k_zB₀, implies that eϕ'/T ∼ k_z/k_⊥ ≪ 1. Thus, the perturbation expansion holds even in the case of strong turbulence if the wave propagation is oblique (e.g., Howes et al. 2008a). The same conclusion could also be reached starting from the ion kinetic equation, with the quasi-neutrality condition ensuring that q_iϕ/T_i ∼ −eϕ'/T_e ≪ 1.

In order to derive fluid-like equations for nonlinear kinetic-Alfvén waves, we now consider the first two moments of the electron kinetic equation (Equation (1)). First, we integrate this equation over velocity and obtain the electron continuity equation

$$\begin{eqnarray} &&\partial n_e/\partial t+\nabla \cdot (n_e {\bf v}_e)=0, \end{eqnarray} \tag{ 88 }$$

where, by definition, n_e = ∫f_ed³v and n_ev_e = ∫f_evd³v. Second, we multiply the equation by v and integrate over velocity to obtain the momentum equation

$$\begin{eqnarray} &&-\nabla _m \Pi ^{lm} - \frac{e n_e}{m_e} E^l-\frac{e n_e}{m_e}[{\bf v}_e\times {\bf B}]^l =0, \end{eqnarray} \tag{ 89 }$$

where Π^lm = ∫v^lv^mf_ed³v, summation over repeated indices is assumed, and we neglected the time derivative (or the electron inertia term) since we are interested in low-frequency fluctuations, ω ≪ k_zv_Te.

To calculate the term Π^lm in the resulting force-balance equation, we need to know the electron distribution function. In principle, the electron distribution function can be calculated from the kinetic equation perturbatively assuming a small parameter eϕ'/T_e. However, the form of the distribution function can be established in a simpler way if we use the conditions ω ≪ k_zv_Te and k_⊥ρ_e ≪ 1. The first condition tells us that we can neglect the term ∂f_e/∂t in the kinetic equation. The second condition allows us to look for the solution of the form (Aleksandrov et al. 1984)

$$\begin{eqnarray} &&f_e=\frac{n_0 {m_e}^{3/2}}{(2\pi T_e)^{3/2}}\exp \left(\frac{e \phi ^{\prime }}{T_e}-\frac{m_e({\bf v}-{\bf v}_e)^2}{2T_e} \right), \end{eqnarray} \tag{ 90 }$$

where T_e = const and the electron drift velocity v_e is much smaller than the electron thermal velocity. The functions ϕ'(x, t) and v_e(x, t) are related to the electric and magnetic fields through

$$\begin{eqnarray} &&-\nabla \phi ^{\prime } ={\bf E}+\frac{1}{c}\left[{\bf v}_e\times {\bf B}\right]. \end{eqnarray} \tag{ 91 }$$

The solution (Equation (90)) approaches the Maxwellian distribution for eϕ'/T_e → 0.⁶ Equation (91) explains the fact already observed in the linear case that the electrons adjust to the potential existing in the co-moving frame. We now obtain from Equation (90) the pressure tensor in the form:

$$\begin{eqnarray} &&\Pi ^{lm}=\delta ^{lm}\frac{T_e}{m_e}n_e, \end{eqnarray} \tag{ 92 }$$

where

$$\begin{eqnarray} &&n_e=n_0\exp (e\phi ^{\prime }/T_e). \end{eqnarray} \tag{ 93 }$$

With the aid of Equations (88) and (92), the force balance in the electron momentum equation can be re-written in a fluid-like form

$$\begin{eqnarray} &&-\frac{1}{n_e m_e}\nabla p_e - \frac{e}{m_e} {\bf E}-\frac{e}{m_ec }[{\bf v}_e\times {\bf B}]=0, \end{eqnarray} \tag{ 94 }$$

where p_e = n_eT_e is the electron pressure. The force-balance equation thus agrees with that of an isothermal electron fluid.

Analogously, in the low-frequency regime, ω ≪ k_⊥v_Ti, the solution for the ion kinetic equation is expected to have the form

$$\begin{eqnarray} &&f_i=\frac{n_0 {m_i}^{3/2}}{(2\pi T_i)^{3/2}}\exp \left(-\frac{q_i \phi }{T_i}-\frac{m_i({\bf v}-{\bf v}_i)^2}{2T_i} \right), \end{eqnarray} \tag{ 95 }$$

where T_i = const and v_i is much smaller than the ion thermal speed. Similarly to the electron case, the potential ϕ and the ion drift velocity satisfy the equation

$$\begin{eqnarray} &&-\nabla _\perp \phi ={\bf E}_\perp +\frac{1}{c}\left[{\bf v}_i\times {\bf B}\right], \end{eqnarray} \tag{ 96 }$$

although in the ion case this equation holds only for the component of the electric field perpendicular to the magnetic field, since the ion motion along the field lines may not satisfy the low-frequency condition $\omega \ll k_\Vert v_{Ti}$. It can be checked that due to k_⊥ρ_i ≫ 1, the last term in Equation (96) must be small⁷ and we can approximate −∇_⊥ϕ ≈ E_⊥. Noting that in the case of $k_\perp \gg k_\Vert$ the perpendicular component of the electric field dominates, $E_\perp \gg E_\Vert$, we see that ϕ represents the potential part of the electric field. Similarly to the derivation of Equation (94), we now integrate the ion kinetic equation with v and derive the part perpendicular to the magnetic field

$$\begin{eqnarray} &&-\frac{1}{n_i m_i}\nabla _\perp p_i + \frac{q_i}{m_i} {\bf E}_\perp +\frac{q_i}{m_i c}[{\bf v}_i\times {\bf B}]=0, \end{eqnarray} \tag{ 97 }$$

where

$$\begin{eqnarray} &&n_i=n_0\exp (-q_i\phi /T_i), \end{eqnarray} \tag{ 98 }$$

and p_i = n_iT_i is the ion pressure. Equation (97) can be used to remove the electric field from the electron force-balance equation (Equation (94)), which allows us to re-write its part perpendicular to the magnetic field in a more familiar fluid-like form

$$\begin{eqnarray} &&-\nabla _\perp p_e -\nabla _\perp p_i +(1/c){\bf J}\times {\bf B}=0, \end{eqnarray} \tag{ 99 }$$

Equations (88), (94), and (97), along with the magnetic induction equation, form a complete set of equations for kinetic-Alfvén waves.

We represent the magnetic field in terms of the mean and fluctuating parts as follows ${\bf B}=(B_0+{b}_z){\hat{z}} + {\bf b}_\perp$. The fact that the z-component of the magnetic fluctuations should be retained, together with B₀, follows from the force-balance equation (Equation (99)). Taking into account $k_\Vert \ll k_\perp$, we derive from this relation

$$\begin{eqnarray} &&-\frac{1}{4\pi }B_0\nabla _\perp b_z - \nabla _\perp p_e - \nabla _\perp p_i=0, \end{eqnarray} \tag{ 100 }$$

which gives an estimate b_z/B₀ ∼ β(δn_e/n₀). For β ≪ 1, the fluctuations of b_z can be neglected, while for β ∼ 1 (the case we consider in this work) they should be retained.

From the force balance in the electron momentum equation (Equation (94)), it follows that the electrons are advected across the magnetic field by the "E × B" drift velocity and the diamagnetic drift velocity, respectively,

$$\begin{eqnarray} &&{\bf v}_{e\perp }= \frac{c}{B^2}{\bf E}\times ({\bf B}_0+b_z{\hat{z}}) +\frac{c}{n_eeB^2}{\nabla p_e}\times ({\bf B}_0+b_z{\hat{z}}),\quad \, \end{eqnarray} \tag{ 101 }$$

where an approximate expression $B^2\approx B_0^2+2B_0b_z$ should be substituted in the denominators of Equation (101). As for the motion of the electrons parallel to the magnetic field it can be expressed through the current $J_\Vert = - e n_e v_{e \Vert }$, since the velocity of ions parallel to the field is relatively small. The perpendicular component of the fluctuating magnetic field can be expressed through the flux function ${\bf b}_\perp = {\hat{z}} \times \nabla \psi$, such that

$$\begin{eqnarray} &&J_\Vert \approx J_z =(c/4\pi)\nabla _\perp \times {\bf b}_\perp = (c/4\pi)\nabla _\perp ^2 \psi. \end{eqnarray} \tag{ 102 }$$

The flux function is the (negative) component of the vector potential parallel to the field, ψ = −A_z. In the considered case of strong anisotropy, $k_\Vert \ll k_\perp$, the component of the fluctuating magnetic field parallel to the field is related to the (non-potential) component of the electric field perpendicular to the field through the magnetic induction equation,

$$\begin{eqnarray} &&-\frac{1}{c}\frac{\partial b_z}{\partial t} =[\nabla \times {\bf E}_\perp ]_z. \end{eqnarray} \tag{ 103 }$$

To derive the equations for nonlinear kinetic-Alfvén waves, we use the electron continuity equation (Equation (88)) and the magnetic induction equation. The force balance parallel to the magnetic field in the electron momentum equation (Equation (94)) gives

$$\begin{eqnarray} &&- \nabla _\Vert p_e -n_e e{\bf E}_\Vert = 0, \end{eqnarray} \tag{ 104 }$$

where the electric field is E = −∇ϕ − (1/c)∂_tA. Supplementing this equation with the electron continuity equation,

$$\begin{eqnarray} &&\partial n_e/\partial t+\nabla \cdot (n_e {\bf v}_{e\perp }) -(1/e) \nabla _{\Vert } J_{\Vert }=0, \end{eqnarray} \tag{ 105 }$$

and using the induction equation (Equation (103)) to express the non-potential part of the electric field, one obtains the following system of equations for the fluctuating parts of the magnetic and density fields:

$$\begin{eqnarray} &&\frac{1}{c}\frac{\partial }{\partial t} {\psi }-\nabla _\Vert {\phi } +\frac{1}{n_e e}\nabla _\Vert {p_e}= 0, \end{eqnarray} \tag{ 106 }$$

$$\begin{eqnarray} &&\frac{\partial }{\partial t} \left[\frac{ n_e}{n_0}-\frac{b_z}{B_0} \right] - \frac{c}{B_0}\nabla {\phi }\times {\hat{z}}\cdot \nabla \left[ \frac{n_e}{n_0}- \frac{b_z}{B_0}\right] \nonumber \\ &&\quad\quad-\frac{c}{eB_0}\nabla \left(\frac{p_e}{n_0}\right)\times {\hat{z}}\cdot \nabla \left(\frac{b_z}{B_0}\right)- \frac{1}{en_0}\nabla _\Vert {J}_\Vert =0. \quad \quad \qquad \end{eqnarray} \tag{ 107 }$$

These equations are written in general form, where the pressure term for the electrons can be arbitrary as long as the fluid description of Equations (88) and (94) holds. Such equations were derived and studied previously in various regimes, (e.g., Hazeltine 1983; Scott et al. 1985; Camargo et al. 1996; Terry et al. 2001; Schekochihin et al. 2009; Smith & Terry 2011; Boldyrev & Perez 2012). In the case relevant to our consideration, the equations of state are isothermal.

In the nonlinear case when the magnetic fluctuations obey k_zB₀ ∼ k_⊥b, we must distinguish the gradients along the guide field B₀, denoted as ∇_z, from the gradients along the local field B = B₀ + b, denoted as $\nabla _\Vert$. The gradients in Equations (106) and (107) parallel to the field are thus defined as

$$\begin{eqnarray} &&\nabla _\Vert =\nabla _z + \frac{1}{B_0}{\hat{z}} \times \nabla {\psi } \cdot \nabla \,. \end{eqnarray} \tag{ 108 }$$

In our discussion of strong kinetic-Alfvén turbulence, we will assume that the fluctuations are anisotropic with respect to the magnetic field in such a way that the linear and nonlinear terms in Equations (107) and (110) are of the same order, that is,

$$\begin{eqnarray} &&B_0 {\nabla }_z \sim {\hat{z}} \times \nabla {\psi } \cdot \nabla. \end{eqnarray} \tag{ 109 }$$

This ordering is analogous to the so-called critical balance condition, k_zB₀ ∼ k_⊥b, of Goldreich & Sridhar (1995); see also Cho & Lazarian (2004), Howes et al. (2011), and TenBarge & Howes (2012). This ordering shows that for k_z ≪ k_⊥, nonlinearity becomes important already for small perturbations. When condition (109) is satisfied, Equations (106) and (107) are essentially nonlinear and three-dimensional.

For completeness, we also note that in the limit of small plasma beta, we have b_z → 0 and the electron continuity equation (107) turns into

$$\begin{eqnarray} &&\frac{\partial }{\partial t} { n_e} - \frac{c}{B_0}\nabla {\phi }\times {\hat{z}}\cdot \nabla {n_e}- \frac{1}{e}\nabla _\Vert {J}_\Vert =0. \end{eqnarray} \tag{ 110 }$$

This equation demonstrates that in the limit of small beta the electrons are advected across the field lines only by the "E × B" drift. Indeed, when the magnetic field strength does not change, that is, b_z = 0, the diamagnetic drift does not advect the electron density. Physically, this situation occurs because the guide centers of particles do not move when a diamagnetic current is present. That is why the diamagnetic drift does not enter Equation (110) even in the case of a general equation of state. However, if the magnetic field strength changes, the magnetic curvature effects do affect the density advection and terms with derivatives of b_z do not cancel out.

The system of equations involving three fields, n_e, ψ, and ϕ ((106) and (110)) or the system of equations involving four fields n_e, b_z, ψ, and ϕ ((106) and (107)) are incomplete, as they have more independent fields than equations. The uniqueness is restored when the systems are supplemented by the equations for the ions. The situation here depends on the scales considered. Briefly deviating from our main discussion, we mention that above the ion-cyclotron scale, ρ_i = v_Ti/Ω_i, a fluid description can be justified for the ions if the ions are cold, which is essentially the limit of low beta. (This case is applicable to most laboratory plasmas, the context in which the corresponding equations were originally derived.) In this case, the ions move across the magnetic field due to the "E × B" drift and the so-called polarization drift. One can write the charge conservation law, $\partial \rho /{\partial t} +\nabla _\perp {\bf J}_\perp +\nabla _\Vert J_\Vert =0$, where the parallel current is given by the electrons, while the perpendicular current is due to the polarization drift of the ions (the "E × B" drifts are the same for ions and electrons; they do not lead to charge separation and do not contribute to the current). It can be shown that quasi-neutrality is preserved for ω_pα ≫ Ω_α, therefore, the charge continuity equation can be reduced to $\nabla _\perp {\bf J}_\perp +\nabla _\Vert J_\Vert =0$, which can be re-written as (e.g., Terry et al. 2001)

$$\begin{eqnarray} &&\frac{n_im_ic^2}{B_0^2}\left[\frac{\partial }{\partial t}\nabla ^2\phi -\frac{c}{B_0}\nabla \phi \times {\hat{z}}\cdot \nabla \nabla ^2\phi \right]=\nabla _{\Vert } J_{\Vert }. \end{eqnarray} \tag{ 111 }$$

Equations (106), (110), and (111) provide the closed three-field system for evolution of electric, magnetic, and density fields in the case of low plasma beta.

If the plasma beta is not small (the case of interest in this work), the ions require kinetic description and a simple fluid model is not justified. We, however, are interested in the subproton, dispersive kinetic-Alfvén waves, that is, we consider scales smaller than the ion gyroscale k_⊥ρ_i ≫ 1. At such scales, the ions are (spatially) not magnetized. Moreover, since we are interested in frequencies smaller than kv_Ti, we have the simple response for the ion-density fluctuations (Equations (98) and (69)). The quasi-neutrality condition n_i = n_e then relates the electric potential to the electron density,

$$\begin{eqnarray} &&\phi =-(T_i/e)\ln (n_e/n_0). \end{eqnarray} \tag{ 112 }$$

Similarly, in the three-field system, the b_z field can be removed from Equation (107) with the aid of the force-balance equation perpendicular to the magnetic field (Equation (100)):

$$\begin{eqnarray} &&b_z=-4\pi (T_i+T_e)\delta n_e/B_0. \end{eqnarray} \tag{ 113 }$$

These expressions for the electric potential ϕ and the magnetic field b_z leave us with only two independent fields, the electron density and the magnetic flux function. Let us normalize these expressions as follows

$$\begin{eqnarray} {\tilde{n}}&=&\frac{v_s}{v_A} \left(1+\frac{T_i}{T_e}\right)^{1/2}{\left[1+\frac{v_s^2}{v_A^2} \left(1+\frac{T_i}{T_e}\right)\right]^{1/2}} \frac{n_e}{n_0}, \quad \end{eqnarray} \tag{ 114 }$$

$$\begin{eqnarray} {\tilde{\psi }} &=& \frac{v_s e}{c T_e} \psi,\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ \end{eqnarray} \tag{ 115 }$$

and normalize the time and the length according to

$$\begin{eqnarray} {\tilde{t}}&=&\frac{\left(1+{T_i}/{T_e}\right)^{1/2}}{\left({\rho _s}/{v_A}\right)\left[ 1+({v_s}/{v_A})^2\left(1+{T_i}/{T_e} \right) \right]^{1/2}}\,t, \end{eqnarray} \tag{ 116 }$$

$$\begin{eqnarray} {\tilde{\bf x}}&=&\frac{\bf x}{\rho _s}. \qquad\qquad\qquad\qquad\qquad\quad\ \ \ \end{eqnarray} \tag{ 117 }$$

In the rest of this section, we will use only the normalized variables and we will omit tilde symbol. The system (Equations (106) and (107)) then takes the form

$$\begin{eqnarray} &\ \partial _t {\psi } + \nabla _\Vert {n}= 0, \end{eqnarray} \tag{ 118 }$$

$$\begin{eqnarray} & \partial _t {n} -\nabla _\Vert \nabla _\perp ^2 \psi =0,\quad\ \end{eqnarray} \tag{ 119 }$$

where $\nabla _\Vert =\nabla _z+ {\hat{z}} \times \nabla {\psi } \cdot \nabla _\perp$. Notice that the nonlinearity enters this system only through the magnetic-line bending effects in $\nabla _\Vert$.

The presented ideal system conserves total energy E and the cross-correlation H:

$$\begin{eqnarray} & E^{ka}=\int \left(|\nabla _\perp \psi |^2 + n^2 \right)d^3 x, \end{eqnarray} \tag{ 120 }$$

$$\begin{eqnarray} & H^{ka}=\int \psi n d^3 x.\qquad\qquad \end{eqnarray} \tag{ 121 }$$

The system (Equations (118) and (119)) possesses linear waves, n_k∝ψ_k∝exp (− iωt + ikx). The linearization is achieved by neglecting the second term in the right-hand side of Equation (108), that is, by replacing ∇_|| → ∇_z, which gives the (dimensionless) dispersion relation for the kinetic-Alfvén waves:

$$\begin{eqnarray} &&\omega ^2= k_z^2 k^2_\perp. \end{eqnarray} \tag{ 122 }$$

The linear modes are characterized by the equipartition between the density and magnetic fluctuations, n_k = ±k_⊥ψ_k.

It is easy to check that the linearized equations produce the same dispersion relation and wave polarization given by the kinetic theory, that is, Equations (41)–(43). Also, note that we do not require the frequencies to be smaller than the ion gyrofrequency, as is implied, e.g., in gyrokinetic treatments, therefore, our results are applicable to region I in Figures 1 and 2 both below and above Ω_i.

As we demonstrated in Section 4, for the considered plasma parameters the whistler waves exist only in the wavenumber region k_z/k_⊥ ≫ 1/(k_⊥ρ_i), which implies that k_z cannot be too small compared with k_⊥. According to the critical balance condition (Equation (109)), the nonlinear effects become important for the whistler waves for magnetic field fluctuations that are significantly larger compared with those in the kinetic-Alfvén case. It is therefore reasonable to expect that whistler turbulence, if it is present in those systems at subproton scales, should predominantly be weak. In this case, the kinetic treatment presented in Section 4 should suffice.

If, however, nonlinearities need to be taken into account, it is convenient to use a nonlinear system for whistler waves. This system is derived if one notes that the frequencies of whistler modes are high, ω ≫ kv_Ti, so that the ions can be considered cold. Moreover, since we are studying subproton scales, k_⊥ρ_i ≫ 1, the ions are non-magnetized both spatially and temporally, ω ≫ Ω_i. The fluid equations for "cold" non-magnetized ions then agree with the single-particle motion equations,

$$\begin{eqnarray} &{ }& d{\bf v}_i/dt=(q_i/m_i){\bf E}, \end{eqnarray} \tag{ 123 }$$

$$\begin{eqnarray} &{ }&\partial n_i/\partial t+\nabla \cdot (n_i{\bf v}_i)=0. \end{eqnarray} \tag{ 124 }$$

For the density fluctuations, these equations lead to exactly the same expression as Equation (70), which demonstrates that for high frequencies, ω ≫ kv_Ti, the ion-density fluctuations can be neglected. The ion fluctuation velocity is estimated from these equations as v_i ∼ (qϕ/T_i)(kv_Ti/ω)v_Ti, while the electron velocity is estimated from the current to be v_e ∼ J/(ne) ∼ ckB/(ne). We therefore obtain that the ion velocity is negligible compared with the electron velocity,

$$\begin{eqnarray} &&\frac{v_i}{v_e}\sim \frac{q\phi _k}{T_i}(k \rho _i) \left(\frac{\omega _{pi}^2}{k^2 c^2} \right)\left(\frac{kv_{Ti}}{\omega } \right) \ll 1, \end{eqnarray} \tag{ 125 }$$

where we used kρ_i ∼ kc/ω_pi ≫ 1 and kv_Ti/ω ≪ 1.

Similar to the kinetic-Alfvén case, we now use the force balance in the electron momentum equation (Equation (94)) and the magnetic induction equation. The force balance is easy to write if one neglects the pressure-gradient term (indeed, the electrons are isothermal and their density does not change)

$$\begin{eqnarray} &&{\bf E}+\frac{1}{c}{\bf v}_e \times {\bf B}=0, \end{eqnarray} \tag{ 126 }$$

where v_e = −J/(n₀e). The induction equation then takes the form

$$\begin{eqnarray} &&\frac{\partial }{\partial t}{\bf B}=-\nabla \times \frac{c}{4\pi n_0 e}\left[ \left(\nabla \times {\bf B} \right)\times {\bf B}\right]. \end{eqnarray} \tag{ 127 }$$

This equation is the electron MHD (EMHD) equation (e.g., Gordeev et al. 1994). A linearization of Equation (127) around a uniform background field B₀ leads to the dispersion relation of whistler waves (Equation (66)).

We now demonstrate how the EMHD equations can be simplified in the presence of a uniform magnetic field B₀ and anisotropic fluctuations, k_z/k_⊥ ≪ 1. In the case of oblique propagation, the system becomes strongly nonlinear for a relatively small fluctuating field, satisfying the critical balance condition b/B₀ ∼ k_z/k_⊥. The EMHD equation (Equation (127)) can then be simplified by keeping only the leading terms in this small parameter (c.f., e.g., Schekochihin et al. 2009). As a result, one obtains what we call the reduced electron MHD (REMHD) equations, which can be conveniently written in terms of the two variables, the z-component of the magnetic field, b_z, and the z-component of the vector potential, ψ = −A_z:

$$\begin{eqnarray} \partial _t \psi &=&-\frac{cB_0}{4\pi n_0 e}\nabla _\Vert b_z, \end{eqnarray} \tag{ 128 }$$

$$\begin{eqnarray} \partial _t b_z &=& \frac{cB_0}{4\pi n_0 e} \nabla _\Vert \nabla ^2_\perp \psi, \end{eqnarray} \tag{ 129 }$$

with the same definition of $\nabla _\Vert$ as in Equation (108).

To elucidate the similarities and differences between the REMHD equations (Equations (128) and (129)) and the previously derived kinetic-Alfvén equations (Equations (118) and (119)), it is instructive to also derive the REMHD equations directly from the electron drift picture, similar to our derivation of the kinetic-Alfvén system (Equations (106) and (107)). For that, we first analyze the force balance in the electron momentum equation. When the electron-density fluctuations are neglected, this equation takes the form of Equation (126). For the considered case of oblique propagation satisfying the critical balance condition, we obtain the following equation for the component perpendicular to the magnetic field

$$\begin{eqnarray} &&{\bf E}_\perp =\nabla _\perp \frac{B_0}{4\pi n_0 e}b_z, \end{eqnarray} \tag{ 130 }$$

which demonstrates that the perpendicular component of the electric field is almost potential, that is, E_⊥ = −∇_⊥ϕ, where

$$\begin{eqnarray} &&\phi =-\frac{B_0}{4\pi n_0 e}b_z. \end{eqnarray} \tag{ 131 }$$

This result is also consistent with the linear case, where the whistler electric field is given by Equations (67) and (68), and it becomes almost potential for k_z ≪ k_⊥.

Substituting the expression for the electric potential (Equation (131)) into Equations (106) and (107) and neglecting the electron-density fluctuations, we rederive the REMHD equations (Equations (128) and (129)). Note the essential difference in the expressions for the electric potential and the electron density used in the derivation of the REMHD equations (Equation (131) and δn_e ≈ 0) and the corresponding expressions used in the derivation of the kinetic-Alfvén equations (Equations (112) and (113)). These differences reflect the different dynamics of the ions in the two systems.

Similarly to the kinetic-Alfvén case, one can introduce the dimensionless variables

$$\begin{eqnarray} {\tilde{\psi }} &=& \frac{v_s e}{c T_e} \psi, \quad {\tilde{b}}_z = \frac{b_z}{B_0}, \end{eqnarray} \tag{ 132 }$$

$$\begin{eqnarray} {\tilde{t}} &=& \frac{v_s^2}{v_A^2}\Omega _i t,\quad {\tilde{\bf x}}=\frac{\bf x}{\rho _s}, \end{eqnarray} \tag{ 133 }$$

which cast the REMHD system (Equations (128) and (129)) into the form

$$\begin{eqnarray} &&\partial _t {\tilde{\psi }} +\nabla _\Vert {\tilde{b}}_z =0, \end{eqnarray} \tag{ 134 }$$

$$\begin{eqnarray} &&\partial _t {\tilde{b}}_z - \nabla _\Vert \nabla ^2_\perp {\tilde{\psi }}=0,\quad \ \ \end{eqnarray} \tag{ 135 }$$

which is structurally equivalent to the kinetic-Alfvén system (Equations (118) and (119)). The corresponding conservation laws are now the conservation of energy and helicity of the magnetic fluctuations

$$\begin{eqnarray} &&E^w=\int \left(|\nabla _\perp {\tilde{\psi }}|^2 + {\tilde{b}}_z^2 \right)d^3 x, \end{eqnarray} \tag{ 136 }$$

$$\begin{eqnarray} &&H^w=\int {\tilde{\psi }} {\tilde{b}}_z d^3 x.\qquad\qquad \end{eqnarray} \tag{ 137 }$$

Although the mathematical structure of the obtained system is formally equivalent to that of the kinetic-Alfvén system, the normalization of the fields (Equations (132) and (134)) is different owing to the difference in the ion dynamics between the two systems. The mathematical similarity of the REMHD system and the kinetic-Alfvén system, however, implies similarity in the scaling relations of the corresponding turbulent states, which are discussed in the following sections.

Before we leave the discussion of the nonlinear equations, we note that the considered subproton kinetic-Alfvén and EMHD equations provide examples of systems where a fluid description can be applied to a collisionless, or weakly collisional, plasma. The applicability of the fluid approximation to collisionless plasmas in certain cases is a known fact; other examples and further discussion can be found in, e.g., Ginzburg (1970) and Aleksandrov et al. (1984).

In a turbulent state, the energy E cascades toward small scales while the cross-correlation H cascades toward large scales. A standard phenomenological approach addressing the energy cascade at small scales is discussed in, e.g., Vainshtein (1973), Biskamp et al. (1999), Ng et al. (2003), Cho & Lazarian (2004), Howes et al. (2008b) and Cho & Lazarian (2009). This approach assumes that in strong turbulence the critical balance condition, which ensures that both linear and nonlinear terms in Equation (108) are of the same order, should be satisfied at all scales. Let us denote n_λ and ψ_λ the typical (rms) fluctuations of the density and the magnetic flux function, respectively, at the field-perpendicular scale λ, and l as the corresponding scale of those fluctuations parallel to the field. By balancing the linear and nonlinear terms in Equation (108), we estimate l ∼ λ²/ψ_λ, in which case the time of nonlinear interaction is comparable to the inverse of the linear frequency (Equation (41)), τ ∼ 1/ω ∼ lλ ∼ λ³/ψ_λ. In addition, we estimate from Equations (118) and (119) that n_λ ∼ ψ_λ/λ. The energy associated with the scale λ can therefore be estimated as $E_\lambda \sim n_\lambda ^2$, and the condition of constant energy flux in the turbulent cascade leads to $n_\lambda ^2/\tau ={\rm const}$, which translates into a scaling of n_λ ∼ ψ_λ/λ ∼ λ^2/3 for turbulent fields. The resulting Fourier energy spectrum of strong kinetic-Alfvén turbulence is

$$\begin{eqnarray} &&E_{KA}(k_\perp)\,dk_\perp \propto k_\perp ^{-7/3}\,dk_\perp. \end{eqnarray} \tag{ 138 }$$

We will numerically study kinetic-Alfvén turbulence based on the dimensionless equations (118) and (119). As was previously explained, the critical balance condition requires that k_z ≪ k_⊥ in the case of strong turbulence. We note, however, that the ideal equations (118) and (119) admit a rescaling ∂/∂t → ∂/∂t, ∇_z → ∇_z, n → n, and ψ → ψ with arbitrary , which preserves the critical balance. This result reflects the fact that equations (118) and (119) lack any frequency scale, since they are derived in the limit of an infinitely large electron gyrofrequency. We may therefore always simultaneously rescale the fields in these equations and the size of the simulation box in the z-direction to satisfy k_z ∼ n ∼ ψ ∼ 1. Such a rescaling is assumed in the numerical simulations discussed below.

To study steady-state turbulence, we supplement Equations (118) and (119) with large-scale random forces that supply the energy to the system:

$$\begin{eqnarray} &&\partial _t {\psi } + \nabla _\Vert {n}= \eta \nabla _\perp ^2 \psi + f_{\psi },\ \end{eqnarray} \tag{ 139 }$$

$$\begin{eqnarray} &&\partial _t {n} - \nabla _\Vert \nabla _\perp ^2 \psi = \nu \nabla _\perp ^2 n + f_n. \qquad \end{eqnarray} \tag{ 140 }$$

The small dissipation terms serve to remove energy at small scales and they are mostly needed to stabilize the code. We solve these equations on a triply periodic cubic domain (L³, L = 1) using standard pseudo-spectral methods. The random forces f_ψ and f_n are applied in Fourier space at wavenumbers 2π/L ⩽ k_{x, y} ⩽ 2(2π/L), k_z = 2π/L. The Fourier coefficients outside the above range are zero and are Gaussian random numbers with amplitudes chosen so that |∇ψ|_rms ∼ 1 inside that range. The individual random values are refreshed independently on average every τ = 0.1L/(2π|∇ψ|_rms). We choose ν = η = 0.01.

The energy spectrum obtained numerically is shown in Figure 3. It turns out that it is steeper than the spectrum −7/3 predicted by phenomenological theories based on dimensional arguments. It is interesting that a spectrum steeper than −7/3 was also inferred from observations of subproton magnetic and density fluctuations in the solar wind (e.g., Sahraoui et al. 2006; Chen et al. 2010, 2012b; Alexandrova et al. 2011, 2012).

Zoom In Zoom Out Reset image size

Various explanations have been proposed for the steeper than −7/3 spectrum of subproton turbulence observed in the solar wind. These explanations include steepening of the spectrum by Landau damping, turbulence weakening, wave–particle interactions, etc. (e.g., Rudakov et al. 2011; Howes et al. 2011). In our analytical and numerical approach, wave–particle interactions are absent, however, the steeper spectrum persists. A possible explanation proposed in Boldyrev & Perez (2012) invoked intermittent corrections that result from two-dimensional structures formed by density and magnetic fluctuations. It was proposed that the spectrum should be close to −8/3, a value consistent with observations and numerical simulations. This fact points to an interesting possibility that the observed scaling is not an artifact of non-universal or dissipative effects; rather, it is an inherent property of the nonlinear turbulent dynamics. The spectrum may therefore be universal, analogous to the Kolmogorov spectrum of fluid turbulence. A definitive numerical study that requires higher numerical resolution will be conducted elsewhere.

We also note that the magnetic energy exceeds the kinetic energy at all scales in the inertial interval. Quite interestingly, the ratio of the density energy to the magnetic energy, the nonlinear electron compressibility, is independent of the wavenumber in the inertial interval, which is similar to the behavior observed in the linear case (cf. Equation (54)). Our additional runs (not shown here) demonstrate that the excess of the magnetic energy is not related to the forcing routine; the excess of magnetic energy in the inertial interval is also observed in simulations where only the density field is forced or where both the density and magnetic fields are forced at large scales. From Figure 3 we estimate |b_k⊥|²/|n_k|² ≈ 1.25 in dimensionless units (Equations (114)–(117)).⁸ We propose that this effect is analogous to the generation of residual energy in MHD turbulence, where magnetic energy also provides the dominant contribution (e.g., Frisch et al. 1975; Müller & Grappin 2005; Podesta et al. 2007; Tessein et al. 2009; Mininni & Pouquet 2009; Boldyrev & Perez 2009; Wang et al. 2011; Boldyrev et al. 2011, 2012; Borovsky 2012; Chen et al. 2013a). We also point out that it is the total energy rather than its density and magnetic field components that exhibits a good scaling in an interval of limited extent, in close analogy with MHD turbulence.

Finally, we discuss the spectrum of the electric field. In our model, its potential part is rigidly related to the density field (cf. Equation (112)). So, we have that E_⊥∝k_⊥n_k and the spectrum of the mostly potential perpendicular electric field should be given by

$$\begin{eqnarray} &&E_{E_\perp }(k_\perp)\,dk_\perp \propto k_\perp ^{-2/3}\, d k_\perp. \end{eqnarray} \tag{ 141 }$$

The spectrum of the parallel electric field is then found from $E_\Vert \propto k_\Vert n_k$, which gives

$$\begin{eqnarray} &&E_{E_\Vert }(k_\perp)\,d k_\perp \propto k_\perp ^{-4/3}\, d k_\perp, \end{eqnarray} \tag{ 142 }$$

where we take into account the anisotropy predicted in the phenomenological model of Boldyrev & Perez (2012), $k_\Vert \propto k_\perp ^{2/3}$. The results of the measurements are shown in Figure 4. The top panel shows the spectrum of the perpendicular component of the electric field, compensated by the predicted k^2/3. This component, by definition, follows the density spectrum in Figure 3. In contrast with the total energy spectrum, the density spectrum by itself does not exhibit a good power law in our inertial range of limited extent and neither does the spectrum of E_⊥. The bottom panel of Figure 4 shows the spectrum of $E_\Vert \equiv {\bf E}\cdot {\bf B}/B$, compensated by the best fit power law in the interval k ∈ [5, 15] (the inertial interval in Figure 3). Although no definite power law is observed, the best fit based on this entire range is given by $k_\perp ^{-2}$, with $k_\perp ^{-1.4}$ fitting better at the lower end and $k_\perp ^{-2.6}$ fitting better at the higher end of this range. Although the −1.4 scaling is close to our phenomenological estimate, the much steeper behavior in the rest of the interval is not yet understood. Note, however, that a steep $E_\Vert$ spectrum at subproton scales has also been recently observed in the solar wind (Mozer & Chen 2013).

Zoom In Zoom Out Reset image size

From a practical point of view, the electric field is not easy to measure since it depends of the frame of the observations (e.g., Bale et al. 2005; Chen et al. 2011). The electric field measured in the moving frame must be transformed to the rest frame of the plasma before its spectrum is calculated. In the case of the solar wind, this procedure requires subtraction of a large quantity, v_sw × B, from the electric field measured by a satellite, where v_sw is the solar-wind velocity relative to the satellite. In contrast, the ${E}_\Vert$ electric field is Galilean invariant and its measurements are frame-independent.

As we pointed out in Section 7, the kinetic-Alfvén equations (118) and (119) are structurally identical to the reduced EMHD equations (128) and (129). This fact ensures that the scaling properties of kinetic-Alfvén turbulence and EMHD turbulence should be identical in the limit $k_\Vert \ll k_\perp$. Strong EMHD turbulence has been numerically studied in the literature previously (e.g., Biskamp et al. 1999; Cho & Lazarian 2009), where a magnetic energy spectrum closer to −7/3 has been reported. As discussed in the previous section, the spectrum observed in our simulations is steeper than −7/3; a complete understanding of this steepening will probably require higher numerical resolution. We should, however, point out that the systems studied in previous works are not identical to our system (Equations (139) and (140)). For example, the studies of Biskamp et al. (1999) addressed three-dimensional decaying EMHD turbulence without a background field. In addition, their equations were also different from Equation (127), as they included terms associated with the finite electron inertia length d_e = c/ω_pe.

Work by Cho & Lazarian (2009) contains simulations of EMHD with a strong background magnetic field, which make their setup closer to ours. However, these simulations also concentrate on decaying turbulence, which does not reach a steady state. The initial conditions in those runs deviate from the critical balance that would be needed in order to ensure strongly nonlinear coupling. Considering their ranges of initially populated $k_\Vert$ and k_⊥, one estimates that the parameter $k_\Vert B_0/(k_\perp b)$ is somewhere between 0.7 and 5 at large scales, so turbulence may be already on the "weaker" side at the beginning of the evolution (in comparison, in our runs, this parameter is kept between 0.5 and 1). As the fluctuations decay, turbulence becomes progressively weaker since the background field does not change. In addition, both papers (Biskamp et al. 1999; Cho & Lazarian 2009) employ hyperviscosities and run at lower spatial resolutions compared with our runs.

The regime of weak whistler turbulence has also been addressed in the literature for the cases of both parallel and oblique propagation with respect to the magnetic field (Livshitz & Tsytovich 1972; Boldyrev 1995; Galtier & Bhattacharjee 2003, 2005). In the limit of oblique propagation, which is more relevant to our consideration, the spectrum of weak whistler turbulence was studied in Boldyrev (1995) and Galtier & Bhattacharjee (2003). In these works, it was found that the kinetic equation for the energy spectrum had a formal steady biscaling solution

$$\begin{eqnarray} &&E(k_\Vert, k_\perp)dk_\Vert dk_\perp \propto k_\Vert ^{-1/2} k_\perp ^{-5/2} dk_\Vert dk_\perp, \end{eqnarray} \tag{ 143 }$$

which should correspond to a constant energy flux in k-space. Galtier & Bhattacharjee (2003) also presented an elegant heuristic analysis that allowed the deviation of the anisotropic whistler spectrum (Equation (143)) from dimensional arguments. The same energy spectrum has also been derived for weak kinetic-Alfvén turbulence by Voitenko (1998), which agrees with the established similarity between the kinetic-Alfvén equations and the REMHD equations; see our discussion at the end of Section 7. It was, however, pointed out in Boldyrev (1995) that the formal solution (Equation (143)) is non-local, that is, the collision integral in the wave kinetic equation diverges logarithmically for such a spectrum at low wavenumbers. The question of whether the formal solution (Equation (143)) or a solution with the exponents close to those in Equation (143), is physically realizable therefore requires further study (e.g., Zakharov et al. 1992, Chapter 3).

We presented a systematic study of subproton turbulence in a collisionless plasma, in the case of plasma beta ∼1. Our work is based on a full kinetic derivation, which complements previously available fluid-like, gyrokinetic derivations and numerical solutions of the Vlasov–Maxwell equations (e.g., Hollweg 1999; Howes et al. 2006, 2008a; Gary & Smith 2009; Schekochihin et al. 2009; Sahraoui et al. 2012). As we demonstrated, the full kinetic derivation allows one to effectively address the dynamics of kinetic-Alfvén and whistler modes in both linear and nonlinear regimes at subproton scales.

These modes have principal differences, which can be characterized from both formal mathematical and physical points of view. The regions of phase space where the kinetic-Alfvén and the whistler modes exist are shown in Figures 1 and 2. The kinetic-Alfvén modes occupy region I in these figures, both below and above the ion-cyclotron frequency Ω_i. More specifically, they exist in the phase space region ω ≪ kv_Ti and k_zρ_i ≪ 1. The whistler modes occupy an essentially different region II in Figure 1; more specifically, the region k_zv_Te ≫ ω ≫ kv_Ti and k_zρ_i ≫ 1.

The physical difference between kinetic-Alfvén and whistler modes is in the dynamics of the ions. In kinetic-Alfvén fluctuations, the ions are dynamically important; they rapidly adjust to the fluctuating electric potential in the plasma. The kinetic-Alfvén fluctuations are essentially compressible and they are described by the kinetic-Alfvén system (Equations (118) and (119)). In the case of whistler modes, on the other hand, the ions are dynamically irrelevant; they provide a uniform stationary background, which, together with the quasi-neutrality condition, ensures effective incompressibility of the whistler modes. The whistler fluctuations are electron fluctuations; they are described by the EMHD equations (Equation (127)) or the REMHD equations (Equations (128) and (129)).

The kinetic-Alfvén and whistler modes share the same dispersion relation scaling, which may cause ambiguity when these waves are identified in numerical simulations or observations. In addition, as we discussed in Section 7, the corresponding non-dimensional fluid-like models have similar structures that may artificially blur their differences. Analysis of the ion dynamics helps to remove a possible ambiguity in such cases. In particular, analysis of the ion dynamics allows one to classify the subproton oscillations obtained by different methods in the literature. For example, according to our results, the modes obtained numerically in Sahraoui et al. (2012, Figure 3) for large propagation angles θ > 80° would be classified as kinetic-Alfvén waves, not whistlers. Indeed, the ion dynamics are essential for these modes since ω < k_⊥v_Ti.⁹ Similarly, the system of equations derived from gyrokinetics in Schekochihin et al. (2009, Equations (226) and (227)) is essentially the kinetic-Alfvén system in our classification; it agrees with our Equations (118) and (119). Although, as we have demonstrated in the present work, these equations are applicable both above and below the ion gyrofrequency, they are different from the REMHD equations (Equations (134) and (135)); they do not describe oblique whistler waves, for example.

Finally, we comment on an interesting possibility that in the case when the kinetic-Alfvén frequency is comparable to the ion-cyclotron frequency, the kinetic-Alfvén modes may couple to the ion-Bernstein modes (e.g., Howes et al. 2008b; Podesta 2012; Chen et al. 2013b). The ion-Bernstein waves have a frequency satisfying $\omega -\Omega _i\sim \Omega _i/\sqrt{k_\perp \rho _i}$. The collisionless ion damping of these modes is negligible when $k_\Vert v_{Ti}/(\omega -\Omega _i)\ll 1$. From the dispersion relation of kinetic-Alfvén waves (Equation (41)), we see that ω^ka ∼ Ω_i is satisfied when $k_\Vert v_{Ti}\sim \Omega _i/(k_\perp \rho _i)$, where we assumed β_i ∼ 1. The collisionless damping is therefore negligible when $\sqrt{k_\perp \rho _i} \gg 1$. In this limit, however, the ion-Bernstein modes occupy a narrow band in frequency space and they may not couple effectively to the kinetic-Alfvén modes. In the opposite case, when k_⊥ρ_i is not large, these modes are relatively strongly damped compared with the kinetic-Alfvén modes (e.g., Podesta 2012) and they also may not be effectively generated.

We are grateful to Christopher Chen and Vladimir Mirnov for useful comments. This work was supported by the U.S. DOE award DE-SC0003888, the DOE grant DE-SC0001794, the NSF grant PHY-0903872, and the NSF Center for Magnetic Self-organization in Laboratory and Astrophysical Plasmas at the University of Wisconsin–Madison. High Performance Computing resources were provided by the Texas Advanced Computing Center (TACC) at the University of Texas at Austin under the NSF-Teragrid Project TG-PHY080013N.