REDUCTION OF COMPRESSIBILITY AND PARALLEL TRANSFER BY LANDAU DAMPING IN TURBULENT MAGNETIZED PLASMAS

Hydrodynamics and magnetohydrodynamics (MHD) are the central descriptions used to study turbulence in solar wind and in a wide range of natural systems. Specifically for solar wind, MHD description yielded great success in our understanding of observational data (e.g., see reviews by Goldstein et al. 1995; Tu & Marsch 1995; Bruno & Carbone 2005; Horbury et al. 2005; Marsch 2006; Ofman 2010). Observational studies show that solar wind is typically found to be only weakly compressible (e.g., see Matthaeus et al. 1991; Bavassano & Bruno 1995, and the reviews cited above) and the usual turbulence models which predict energy spectra are derived in the framework of an incompressible MHD description (Iroshnikov 1963; Kraichnan 1965; Goldreich & Sridhar 1995, 1997; Galtier et al. 2000, 2002; Boldyrev 2005, 2006; Lithwick et al. 2007; Chandran 2008; Perez & Boldyrev 2009; Podesta & Bhattacharjee 2010, see also reviews by Cho et al. 2003; Zhou et al. 2004; Galtier 2009; Sridhar 2010). Theoretical models which describe the radial evolution of spatially averaged solar wind quantities are also usually developed in the framework of incompressible MHD formulated in Elsässer variables (e.g., Zhou & Matthaeus 1989, 1990a, 1990b; Marsch & Tu 1989; Zank et al. 1996; Smith et al. 2001; Matthaeus et al. 2004; Breech et al. 2008; Zank et al. 2011). It is however well known that solar wind is not completely incompressible and many phenomena which require compressibility are observed in solar wind, such as the evolution of density fluctuations (e.g., Armstrong et al. 1990; Coles et al. 1991; Grall et al. 1997; Woo & Habbal 1997; Bellamy et al. 2005; Wicks et al. 2009; Telloni et al. 2009), magnetic holes, solitons, and mirror mode structures (e.g., Winterhalter et al. 1994; Fränz et al. 2000; Stasiewicz et al. 2003; Stevens & Kasper 2007), or strong temperature anisotropies which trigger, and are limited by, micro-instabilities such as mirror and fire hose (e.g., Gary et al. 2001; Kasper et al. 2002; Hellinger et al. 2006; Matteini et al. 2007; Bale et al. 2009). The importance of compressibility was stressed by Carbone et al. (2009), who compared the observational data with the energy flux scaling laws of Politano & Pouquet (1998), which are exact relations of incompressible MHD. Carbone et al. determined that the scaling relations can better fit the data if the relations are phenomenologically modified to account for compressibility. The incorporation of weakly compressional density fluctuations was partially addressed by so-called nearly incompressible models, which expand the compressible equations with respect to small sonic Mach number (Matthaeus & Brown 1988; Zank & Matthaeus 1990, 1991, 1993) and which were recently formulated in the presence of a static large-scale inhomogeneous background (Hunana & Zank 2010; Hunana et al. 2006, 2008, see also Bhattacharjee et al. 1998). These models however specifically assume, and cannot explain, why solar wind is only weakly compressible. The theoretical compressible MHD models developed in the wave turbulence formalism (e.g., Kuznetsov 2001; Chandran 2005) also cannot address this issue.

Describing solar wind with fully compressible MHD or Hall-MHD formalisms yields several problems. Most importantly, these compressible descriptions introduce sound waves and slow magnetosonic waves. As elaborated by Howes (2009), slow magnetosonic waves are strongly damped by Landau damping in the kinetic Maxwell–Vlasov description. The presence of fast and slow magnetosonic waves naturally implies higher level of compressibility and overestimates the parallel energy transfer. Also, numerical simulations of compressible Hall-MHD performed by Servidio et al. (2007) in the context of the magnetopause boundary layer showed that compared to the usual turbulence in compressible MHD simulations, which consists of Alfvén waves, the Hall term is responsible for decoupling of magnetic and velocity field fluctuations. In the compressible Hall-MHD regime, Servidio et al. observed spontaneous generation of magnetosonic waves which transform to a regime of quasi perpendicular "magnetosonic turbulence." This is in contrast with observational studies which typically show that turbulence in solar wind predominantly consists of quasi perpendicular (kinetic) Alfvén waves (e.g., Bale et al. 2005; Sahraoui et al. 2009, 2010). Finally, all usual MHD or Hall-MHD models cannot address how the energy is actually dissipated at small scales. It is well known that solar wind plasma is almost collisionless and therefore the classical von Karman picture of energy being dissipated via viscosity is not applicable to solar wind. It is evident that new and more realistic models have to be introduced which can overcome some of these drawbacks.

The most realistic approach is of course the fully kinetic Vlasov–Maxwell description. It is however analytically intractable, and even the biggest kinetic simulations cannot resolve the large-scale turbulence dynamics. Two leading approaches appear to be promising to substitute MHD in describing solar wind turbulence: gyrokinetics and Landau fluids. Gyrokinetics (e.g., Schekochihin et al. 2009; Howes et al. 2006) was originally developed for simulations of fusion in tokamaks. It is a kinetic-like description, which averages out the gyrorotation of particles around a mean magnetic field and therefore makes the kinetic description more tractable, mostly by eliminating fast timescales. Derived directly from the kinetic theory, gyrokinetics has a crucial advantage of being asymptotically correct. Landau fluid description, on the other hand, is a fluid-like extension of compressible Hall-MHD, in which wave dissipation is incorporated kinetically by the modeling of linear Landau damping, thus retaining a realistic sink of energy. Other linear kinetic effects such as finite Larmor radius (FLR) corrections are also incorporated.

The simplest Landau fluid closure was considered by Hammett & Perkins (1990) and the associated dispersion relations were numerically explored by Jayanti et al. (1998). The Landau fluid model was further advanced by Snyder et al. (1997), who considered the largest MHD scales, starting from the guiding center kinetic equation. The approach was reconsidered and refined to its present form with incorporation of Hall term and FLR corrections by Passot & Sulem (2007, 2003), Goswami et al. (2005), and Bugnon et al. (2004). This Landau fluid approach starts with the Vlasov–Maxwell equations and derives nonlinear evolution equations for density, velocity, and gyrotropic pressures. In the simplest formulation the model is closed at the level of heat fluxes by matching with the linear kinetic theory in the low-frequency limit. Kinetic expressions usually contain the plasma dispersion function which is not suitable for fluid-like simulations. Landau fluid closure is therefore performed in a way as to minimize occurrences of this plasma dispersion function and, where not possible, this function is replaced by a Padé approximant. This eliminates the time non-locality and also results in the presence of a Hilbert transform with respect to the longitudinal coordinate (in the direction of the ambient magnetic field) which in the fluid formalism is associated with linear Landau damping. Further details about the development of Landau fluid models are thoroughly discussed in the papers cited above. The Landau fluid approach should however be contrasted with the more classical gyrofluid models (e.g., Dorland & Hammett 1993; Brizard 1992; Scott 2010), which are derived by taking fluid moments of the gyrokinetic equation and for which a similar closure scheme is applied afterward.

The Landau fluid approach has the following advantages. In contrast with Hall-MHD, it contains separate equations for parallel and perpendicular pressures and heat fluxes. It therefore allows for the development of temperature anisotropy, which is observed in solar wind. Noticeably, in contrast with gyrokinetics or gyrofluid models, the Landau fluid model does not average out the fast waves. Compared with these approaches, it also has an advantage in that the final equations including the FLR corrections are written for the usual quantities measured in the laboratory frame. Existing spectral MHD and Hall-MHD codes can therefore be modified to Landau fluid description relatively easily. Also importantly, even though gyrokinetics is a reduced kinetic description, it is still five-dimensional and therefore naturally quite difficult to compute. While current largest numerical simulations of gyrokinetics (Howes et al. 2008) require thousands CPU-cores for a fluid-like 128 × 64² resolution, the FLR-Landau fluid model requires computational power only slightly larger than the usual Hall-MHD simulations. The results presented here, which employ a resolution of N = 128 grid points in all three directions, were calculated using 32 CPU-cores.

Landau fluid models can be developed with several levels of complexity. For these first Landau fluid simulations of three-dimensional turbulence, we use a simplified version of the most general Landau fluid model (Passot & Sulem 2007), where we constrain ourselves to isothermal electrons and leading order corrections in terms of the ratio of the ion Larmor radius to the considered scales. A similar model was used by Borgogno et al. (2009), who studied the dynamics of parallel propagating Alfvén waves in a medium with an inhomogeneous density profile. They numerically showed that the observed Alfvén wave filamentation and later transition to the regime of dispersive phase mixing is consistent with particle-in-cell simulations. In this paper, we concentrate on freely decaying turbulence and compare these to simulations of compressible Hall-MHD. Numerical integration of the full Landau fluid model in one space dimension was presented by Borgogno et al. (2007), who investigated the dynamics very close to the mirror instability threshold and showed the presence of magnetic holes. Results considering quasi-transverse one-dimensional propagation in the full Landau fluid model are presented in Laveder et al. (2011a, 2011b).

Considering a neutral bi-fluid consisting of protons (ions) and isothermal electrons, the Landau fluid model consists of evolution equations for proton density ρ = m_pn (where m_p is the proton mass and n is the number density), proton velocity $\boldsymbol{u}_p$, proton parallel and perpendicular pressures p_∥p, p_⊥p, and heat fluxes q_∥p, q_⊥p, together with the induction equation for magnetic field $\boldsymbol{b}$. The equations are normalized and density, magnetic field and proton velocity are measured in units of equilibrium density ρ₀, ambient magnetic field B₀, and Alfvén speed $V_A=B_0/\sqrt{4\pi \rho _0}$, respectively. Pressures are measured in units of initial proton parallel pressure p⁽⁰⁾_∥p and heat fluxes in units of p⁽⁰⁾_∥pV_A. The total pressure in the momentum equation has the form of a tensor. Defining $\hat{\boldsymbol{b}}=\boldsymbol{b}/|\boldsymbol{b}|$ as a unit vector in the direction of local magnetic field, the proton pressure can be cast in the form $\boldsymbol{p}_p=p_{\perp p} \boldsymbol{n} + p_{\parallel p} \boldsymbol{\tau } + \boldsymbol{\Pi }$, where $\boldsymbol{\tau }=\hat{\boldsymbol{b}}\otimes \hat{\boldsymbol{b}}$, and $\boldsymbol{n}=\boldsymbol{I}-\hat{\boldsymbol{b}}\otimes \hat{\boldsymbol{b}}$, with $\boldsymbol{I}$ being the unit tensor. FLR corrections to the gyrotropic pressures are represented by $\boldsymbol{\Pi }$. Operator ⊗ represents the usual tensor product and in the index notation, for example, $\tau _{ij}=\hat{b}_{i}\hat{b}_{j}$. Electrons are assumed to be isothermal with the scalar pressure p_e = nT⁽⁰⁾_e, where T⁽⁰⁾_e is the electron temperature. Parameter R_i, whose inverse multiplies the Hall-term and also the FLR corrections, is defined as R_i = L/d_i, where d_i is the ion inertial length and L is the unit length. The proton plasma beta is defined with respect to parallel pressure and β_∥ = 8πp⁽⁰⁾_∥p/B²₀. The density, momentum, and induction equations of the FLR-Landau fluid model can then be expressed as

$$\begin{eqnarray} && \frac{\partial \rho }{\partial t} + \boldsymbol{\nabla }\cdot (\rho \boldsymbol{u}_p) = 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} && \frac{\partial \boldsymbol{u}_p}{\partial t} + \boldsymbol{u}_p\cdot \boldsymbol{\nabla }\boldsymbol{u}_p+ \frac{\beta _\parallel }{2\rho } \boldsymbol{\nabla }\cdot (p_{\perp p} \boldsymbol{n} + p_{\parallel p} \boldsymbol{\tau } + \boldsymbol{\Pi } + p_e \boldsymbol{I}) \nonumber \\ \ms && \quad -\frac{1}{\rho } (\boldsymbol{\nabla }\times \boldsymbol{b})\times \boldsymbol{b}= 0, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} && \frac{\partial \boldsymbol{b}}{\partial t} = \boldsymbol{\nabla }\times (\boldsymbol{u}_p\times \boldsymbol{b}) - \frac{1}{R_i}\boldsymbol{\nabla }\times \left[ \frac{1}{\rho } (\boldsymbol{\nabla }\times \boldsymbol{b})\times \boldsymbol{b}\right]. \end{eqnarray} \tag{ 3 }$$

Dropping, for simplicity, indices p for proton velocity $\boldsymbol{u}_p$ and proton pressures p_⊥p, p_∥p, the evolution equations for perpendicular and parallel pressures read (neglecting the work done by the FLR stress forces)

$$\begin{eqnarray} && \frac{\partial p_\perp }{\partial t} + \boldsymbol{\nabla }\cdot (p_\perp \boldsymbol{u}) +p_\perp \boldsymbol{\nabla }\cdot \boldsymbol{u}-p_\perp \hat{\boldsymbol{b}}\cdot (\boldsymbol{\nabla }\boldsymbol{u})\cdot \hat{\boldsymbol{b}}\nonumber \\ && \quad + \boldsymbol{\nabla }\cdot (q_\perp \hat{\boldsymbol{b}}) + q_\perp \boldsymbol{\nabla }\cdot \hat{\boldsymbol{b}}= 0, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} && \frac{\partial p_\parallel }{\partial t} + \boldsymbol{\nabla }\cdot (p_\parallel \boldsymbol{u}) + 2 p_\parallel \hat{\boldsymbol{b}}\cdot (\boldsymbol{\nabla }\boldsymbol{u})\cdot \hat{\boldsymbol{b}}\nonumber\\ &&\quad +\boldsymbol{\nabla }\cdot (q_\parallel \hat{\boldsymbol{b}}) -2q_\perp \boldsymbol{\nabla }\cdot \hat{\boldsymbol{b}}= 0. \end{eqnarray} \tag{ 5 }$$

Assuming an ambient magnetic field of amplitude B₀ in the positive z-direction, a semi-linear description of the FLR corrections in the pressure tensor neglecting heat flux contributions can be expressed as

$$\begin{eqnarray} && \Pi _{xx} = -\Pi _{yy} = -\frac{\langle p_\perp \rangle }{2R_i} (\partial _y u_x +\partial _x u_y), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} && \Pi _{xy} = \Pi _{yx} = -\frac{\langle p_\perp \rangle }{2R_i}(\partial _y u_y - \partial _x u_x), \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} && \Pi _{yz} = \Pi _{zy} = \frac{1}{R_i} [ 2\langle p_\parallel \rangle \partial _z u_x + \langle p_\perp \rangle (\partial _x u_z -\partial _z u_x) ], \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} && \Pi _{xz} = \Pi _{zx} = -\frac{1}{R_i}[ 2\langle p_\parallel \rangle \partial _z u_y + \langle p_\perp \rangle (\partial _y u_z -\partial _z u_y) ], \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} && \Pi _{zz} = 0, \end{eqnarray} \tag{ 10 }$$

where 〈p_⊥〉 and 〈p_∥〉 represent the instantaneous averaged ion pressures over the entire domain, whose time variation is aimed to take into account the evolution of the global properties of the plasma. Finally, parallel and perpendicular heat fluxes q_∥, q_⊥ evolve according to

$$\begin{eqnarray} && \left(\frac{d}{dt} + \frac{\sqrt{\pi \beta _\parallel }}{4(1-\frac{3\pi }{8})} \mathcal {H} \partial _z \right) q_\parallel = \frac{1}{1-\frac{3\pi }{8}} \frac{\beta _\parallel }{2} \partial _z (p_\parallel -\rho), \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} && \left(\frac{d}{dt} - \frac{\sqrt{\pi \beta _\parallel }}{2} \mathcal {H} \partial _z \right) q_\perp =\frac{\beta _\parallel }{2} \frac{T_{\perp p}^{(0)}}{T_{\parallel p}^{(0)}} \partial _z \nonumber \\ && \quad \times \left[ \left(1-\frac{T_{\perp p}^{(0)}}{T_{\parallel p}^{(0)}} \right) |b| - \left(\frac{T_{\parallel p}^{(0)}}{T_{\perp p}^{(0)}} p_\perp - \rho \right) \right], \end{eqnarray} \tag{ 12 }$$

where T⁽⁰⁾_⊥p, T_∥p⁽⁰⁾ are the initial perpendicular and parallel proton temperatures and d/dt is the convective derivative. The operator $\mathcal {H}$, which is defined as

$$\begin{equation} \mathcal {H} f(z) = - \frac{1}{\pi } VP \int _{-\infty }^{+\infty } \frac{f(z^{\prime })}{z-z^{\prime }} dz^{\prime }, \end{equation} \tag{ 13 }$$

reduces in the Fourier space to a simple multiplication by ik_z/|k_z| and is the signature of linear Landau damping.

To gain physical insight into a quite complicated model (Equations (1)–(12)), it is useful to momentarily consider just the largest scales by putting 1/R_i → 0, which eliminates the nongyrotropic contributions to the pressure tensor and which also eliminates the Hall term. The resulting set of equations still contains linear Landau damping and with the exception of isothermal electrons, it is analogous to the model of Snyder et al. (1997). If this model is further simplified by elimination of Equations (11) and (12) and by instead prescribing q_∥ = q_⊥ = 0, the resulting model collapses to the double adiabatic model (Chew et al. 1956, see also Kulsrud 1983) and Landau damping disappears. The presence of ion Landau damping in the system (Equations (1)–(12)) is therefore a result of closure Equations (11) and (12) for the heat fluxes, which contain the operator ${\mathcal H}$. At least in the static limit, this closure can be viewed as a modified Fick's law where the gradient operator that usually relates the heat flux to the temperature fluctuations is here replaced by a Hilbert transform that is a reminiscence of the plasma dispersion function arising in the linear kinetic theory, and is a signature of Landau (zero-frequency) wave-particle resonance. The effect of Landau damping in the system (Equations (1)–(12)) might be better understood by solving the linearized set of equations. In the Appendix, we consider waves which propagate parallel to the ambient magnetic field. It is shown that except the Alfvén waves, linear waves have a frequency with a negative imaginary part and are therefore damped. This corresponds to linear Landau damping.

The model used for the simulations presented here has several limitations. First of all, it contains electrons which are assumed isothermal, a regime in fact often assumed in hybrid simulations which provide a kinetic description of the ions and a fluid description of the electrons. In the future, more realistic simulations will be performed with inclusion of independent evolution equations for parallel and perpendicular electron pressures and heat fluxes. This will also result in the presence of electron Landau damping, which is absent in the model presented here and which seems to play an important role in solar wind turbulence. Another main limitation appears to be the form of FLR corrections, which are derived as a large-scale limit of FLR corrections of the full Landau fluid model and which are therefore significantly simplified. The FLR corrections (6)–(10) are sufficient for the simulations of freely decaying turbulence, which do not lead to significant temperature anisotropies. However, our preliminary simulations which employ forcing and lead to strong temperature anisotropies show that the FLR corrections (6)–(10) are overly simplified and might lead to artificial numerical instabilities. For simulations with strong temperature anisotropies, a more refined description of FLR corrections is required (see Passot & Sulem 2007; Borgogno et al. 2007).

To explore the behavior of the FLR-Landau fluid model, we performed simulations of freely decaying turbulence. The code we used is based on a pseudo-spectral discretization method, where spatial derivatives are evaluated in Fourier space. The time stepping is performed in real space with a third-order Runge–Kutta scheme. Spatial resolution is N³ = 128³ and the size of the simulation domain is L = 16 × (2π) in each direction. The Hall parameter is R_i = 1, implying that the lengths are measured in the units of ion inertial length d_i. Velocity and magnetic field fluctuations of root mean square amplitude $\langle \boldsymbol{u}^2 \rangle ^{1/2} = \langle \boldsymbol{b}^2 \rangle ^{1/2} = 1/8$ are initialized in Fourier space on the first four modes kd_i = m/16, where m ∈ [1, 4], with flat spectra and random phases and are constructed to be divergence free. Constant pressures p_⊥ = p_∥ = 1, density ρ = 1, and heat fluxes q_⊥ = q_∥ = 0 are initialized in the entire domain. The temperature of the electrons is T⁽⁰⁾_e = 1 and is therefore equal to proton temperatures, which can be defined as T_⊥ = 〈p_⊥/ρ〉 and T_∥ = 〈p_∥/ρ〉. In the simulations presented here, both T_⊥ and T_∥ stay rather close to their initial value. The compressible Hall-MHD model consists of Equations (1)–(3), where the divergence of the pressure tensor in Equation (2) is substituted with the usual gradient of scalar pressure and, assuming the adiabatic law p = ρ^γ, in the normalized units it is equal to $\beta _0/\gamma \boldsymbol{\nabla }\rho ^\gamma$, where β₀ is the usual plasma beta defined as β₀ = c²_s/V²_A and the sound speed c²_s = γp/ρ. Adiabatic index γ = 1.66 is used.

As shown, for example, by Hirose et al. (2004) and Howes (2009), it is not straightforward to compare Hall-MHD and Vlasov–Maxwell kinetic theory because the associated dispersion relations are quite different if in the kinetic description the proton (ion) temperature T_p is not negligible with respect to the electron temperature T_e. We choose to follow Howes (2009), who, in order to compare Hall-MHD with the kinetic theory (which is represented here by the FLR-Landau fluid model), defined the necessary relation between β₀ and β_∥ as β₀ = β_∥(1 + T_p/T_e)/(2T_p/T_e). For equal proton and electron temperatures T_p = T_e this relation yields β₀ = β_∥. Landau damping alone is not sufficient to run the code and some kind of artificial dissipation is needed to terminate the cascade. We here resorted to use a filtering in the form of $\frac{1}{2}\lbrace 1-\tanh [ (m-0.8N/4)/3]\rbrace$, where m is the mode index, applied each time step on all fields for both Hall-MHD and Landau fluid regimes. Because of the filtering, it is crucial to have identical time steps dt = 0.128 in both models. In most of the simulations, we used β₀ = β_∥ = 0.8. In normalized units, the Alfvén speed is equal to unity, the usual MHD sound speed $c_s=\sqrt{\beta _0}=0.894$ and the turbulent sonic Mach number $M_s=\langle \boldsymbol{u}^2 \rangle ^{1/2}/c_s=0.14$. In Section 4, we also consider simulations with β₀ = β_∥ = 0.25 (M_s = 0.25) and β₀ = β_∥ = 0.1 (M_s = 0.4).

A wave analysis procedure was implemented in the code, which consists in choosing few modes in spatial Fourier space for each field and recording their value every 20 time steps. After the run, the time Fourier transform is performed and frequency-power spectra are obtained for each mode. This procedure makes it possible to identify at which frequency there is maximum power for each mode, and by comparing these with frequencies obtained from theoretical dispersion relations it allows to uniquely identify which waves are present in the system. This method was previously used for incompressible MHD by Dmitruk & Matthaeus (2009), therefore detecting Alfvén waves. Considering the dynamics associated with waves propagating in the direction parallel to the ambient magnetic field (propagation angle θ = 0°), Figure 1 shows frequency-power spectra of four modes with wavenumbers k_x = 0, k_y = 0, k_zd_i = m/16, with m = 1, 2, 4, 8, recorded from the component b_x for Hall-MHD (left) and for FLR-Landau fluid (right). For Hall-MHD, these correspond to polarized Alfvén waves which obey the dispersion relation $\omega = \pm k^2/(2R_i) + k\sqrt{1+(k/2R_i)^2}$. The theoretical frequency values which are expected from this dispersion relation for given k are plotted on the top axis of the figure as black vertical lines. Figure 1 shows that the resolution of 128³ is sufficient to clearly distinguish between left and right polarized Alfvén waves. In contrast with Hall-MHD, for which it is possible to write the general dispersion relation for frequency ω as a relatively simple polynomial of sixth order, for FLR-Landau fluids the general dispersion relation would be uneconomically large to write down. In general, it is necessary to numerically solve the determinant obtained from linearized Equations (1)–(12) for a given wavenumber k after assuming linear waves. FLR-Landau fluid (Equations (1)–(12)) consists of 11 evolution equations in 11 variables and, together with the divergence free constraint for the magnetic field, therefore yields general dispersion relation for frequency ω in the form of a complex polynomial of 10th order. This represents five forward and five backward propagating waves, with some solutions having highly negative imaginary part and which are therefore strongly damped. A similar situation is encountered in the Vlasov–Maxwell kinetic theory which essentially yields an infinite number of strongly damped solutions. For the propagation angle θ = 0° and additional constraint T_⊥ = T_∥ = 1, it is however possible to obtain an analytic solution for the circularly polarized Alfvén waves as

$$\begin{equation} \omega = \pm \frac{k^2}{2R_i}\left(1+\frac{\beta _\parallel }{2} \right) +k\sqrt{1+\left(\frac{k}{2R_i}\right)^2 \left(1-\frac{\beta _\parallel }{2}\right)^2}, \end{equation} \tag{ 14 }$$

with two other solutions obtained by substituting ω with −ω. The dispersion relation is quite similar to that of Hall-MHD with additional terms proportional to β_∥ and resulting from the FLR corrections. Expected theoretical frequencies obtained from this analytic solution are plotted on the top axis of Figure 1 for the Landau fluid regime (right). They again match quite precisely. Note also the moderately strong damping of parallel Alfvén waves in Landau fluid regime, which can be seen for the last mode m = 8. Landau damping does not act directly on linear constant amplitude Alfvén waves obeying relation (14), which are exact solutions of linearized FLR-Landau fluid model. However, nonlinear parallel Alfvén waves, which are of course present in the full model, cause production of density (sound) fluctuations. Sound waves in the FLR-Landau fluid model, as well as in the kinetic theory, are heavily damped by Landau damping as shown below and this process therefore also results in damping of Alfvén waves. Mjølhus & Wyller (1988) studied the kinetic derivative nonlinear Schrödinger equation (KDLNS) for parallel propagating long-wavelength Alfvén waves where they refer to this effect as nonlinear Landau damping, because it is acting on Alfvén waves in a nonlinear way.

Zoom In Zoom Out Reset image size

Further exploring propagation angle of θ = 0°, Figure 2 shows frequency-power spectra obtained from component u_z and which should therefore predominantly display sound waves (almost identical spectra can be obtained from component ρ). The same modes with k_zd_i = m/16, where m = 1, 2, 4, 8 are shown as in Figure 1. For Hall-MHD, parallel sound waves obey the dispersion relation ω = kc_s, where the sound speed $c_s=\sqrt{\beta _0}$. Sound waves (S) are clearly presented in Figure 2 for Hall-MHD (left) with quite sharp peaks which match the theoretical dispersion values shown on the top axis. Weak presence of polarized Alfvén waves (A) for modes m = 1, 2 is also visible in this component generated by nonlinear coupling. For the Landau fluid regime (Figure 2, right), it is not possible to obtain simple analytic dispersion relation which corresponds to sound waves and correct values must be obtained numerically as explained above (see also the Appendix). This yields five frequencies for forward propagating waves, with two solutions corresponding to the polarized Alfvén waves (14) and three solutions which are highly damped. The sound wave was chosen as the least damped solution. Dispersion relation $\omega =k\sqrt{\vphantom{A^R}\smash{\hbox{${(3+T_e^{(0)})\beta _\parallel /2}$}}}$ obtained from the double adiabatic model, to which Landau fluid description collapses after assumption of zero heat fluxes and large scales, can also be used as a heuristic guide to determine which out of the three damped solutions represents the sound wave frequency. Obtained frequency values are again plotted on the top axis. Figure 2 shows that the sound waves (S) are heavily damped for the Landau fluid regime (Figure 2, right) with the nonlinearly generated Alfvén waves (A) completely overpowering the spectra. Inhibition of sound waves is consistent with the kinetic theory. As elaborated by Howes (2009), sound waves are overestimated by MHD and Hall-MHD descriptions as they represent solutions which are strongly damped by the Landau resonance. Note that the calculated sound wave frequencies in the FLR-Landau fluid regime are actually higher than the corresponding Alfvén wave frequencies. A similar result was obtained by Howes (2009) who numerically compared dispersion relations for Hall-MHD and kinetic theory and who, for the nearly parallel propagation with T_i = T_e and β₀ = β_∥ = 1, noted that the kinetic solution corresponding to the slow wave has a higher phase speed than the kinetic solution corresponding to the Alfvén wave.

Zoom In Zoom Out Reset image size

Considering the propagation angle of θ = 45°, Figure 3 shows frequency-power spectra in component b_x for wavenumbers k_y = 0, k_xd_i = k_zd_i = m/16, where m = 1, 2, 4, 6 and which therefore predominantly shows slow (S) and fast (F) magnetosonic waves. Spectra also show weaker presence of Alfvén waves (A). Again for this angle of propagation, the slow waves are strongly damped. In both Hall-MHD and FLR-Landau fluid simulations, the associated dispersion relations had to be solved numerically and the predicted frequencies for spectral peaks are shown on the top axis. Slow waves (S) in the FLR-Landau fluid regime were again identified as the least damped frequency out of the three highly damped solutions. For identical wavenumbers, the presence of Alfvén waves can be better explored in the component b_y and corresponding spectra are shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, considering purely perpendicular propagation with θ = 90°, Figure 5 shows frequency-power spectra obtained from the density field ρ and wavenumbers k_y = 0, k_z = 0, k_xd_i = m/16, where m = 1, 2, 4, 8. For Hall-MHD regime (Figure 5, left), the spectral peaks correspond to magnetosonic waves with the usual dispersion relation $\omega =k\sqrt{1+\beta _0}$. In the Landau fluid regime, linearized set of Equations (1)–(12) with the additional constraint T_⊥ = T_∥ = 1 can be shown to yield the dispersion relation for perpendicular magnetosonic waves in the form

$$\begin{equation} \omega = k\sqrt{1+\beta _\parallel \left(1+\frac{T_e^{(0)}}{2}\right) + \left(\frac{\beta _\parallel k}{4R_i}\right)^2}. \end{equation} \tag{ 15 }$$

The dispersion relation clearly shows the effect of inclusion of isothermal electrons (for simulations presented here T⁽⁰⁾_e = 1) and also the effect of FLR corrections which are represented by the last quadratic term. Theoretical predictions from Hall-MHD and FLR-Landau fluid dispersion relations are again shown on the top axis. To clearly show the shift of the peaks between the two regimes, we also added the theoretical Hall-MHD frequencies to the top axis of FLR-Landau fluid regime (Figure 5, right) and represent them with the small magenta lines.

Zoom In Zoom Out Reset image size

Compressibility of the flow can be evaluated by decomposing the velocity field into its solenoidal and non-solenoidal components and by calculating the associated energies according to

$$\begin{equation} \sum _{\boldsymbol{k}} |\boldsymbol{u_k}|^2 = \sum _{\boldsymbol{k}} \frac{|\boldsymbol{k}\times \boldsymbol{u_k}|^2}{|\boldsymbol{k}|^2} + \sum _{\boldsymbol{k}} \frac{|\boldsymbol{k}\cdot \boldsymbol{u_k}|^2}{|\boldsymbol{k}|^2}, \end{equation} \tag{ 16 }$$

where the left-hand side corresponds to the total energy E^U in velocity field, the first term in the right-hand side corresponds to the energy E_in in the solenoidal component and the second term E_c originates from the compressible one. Relation (16) can be therefore expressed as E^U = E_in + E_c and the compressibility of the flow can be evaluated as a ratio of compressible and total energy E_c/E^U. Time evolution of E_c/E^U for Hall-MHD and FLR-Landau fluid regime with β₀ = β_∥ = 0.8 is presented in Figure 6. The figure shows that the ratio E_c/E^U which represents compressibility is significantly lower in Landau fluid regime and is therefore a result of presence of Landau damping. The question then arises of the influence of the sonic Mach number in the compressibility evolution. The time evolution of E_c/E^U for β₀ = β_∥ = 0.25 (which corresponds to M_s = 0.25) is displayed in Figure 7 (left), whereas the time evolution of E_c/E^U for β₀ = β_∥ = 0.1 (which corresponds to M_s = 0.40) is shown in Figure 7 (right). The simulations were performed with the same time step dt = 0.128 and filtering, for approximately half the total integration time of the previous regime with β₀ = β_∥ = 0.8 (M_s = 0.14). The figure shows that while the compressibility is still clearly reduced in Landau fluid simulation with β₀ = β_∥ = 0.25, this reduction is almost insignificant for simulations with β₀ = β_∥ = 0.1. This is an expected effect, as the strength of Landau damping is proportional to β_∥. We note that the turbulent sonic Mach number in solar wind is typically small and, for example, analysis of observational data performed by Bavassano & Bruno (1995) (from 0.3 to 1.0 AU) showed that the most probable value is between M_s = 0.1–0.2 with the distribution having an extended tail to lower values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The question also arises of the compared evolution of total energy for Hall-MHD and FLR-Landau fluid simulations. The total energy E_tot can be evaluated as the sum of the kinetic energy E_kin, the magnetic energy E_mag, and the internal energy E_int. In Hall-MHD and FLR-Landau fluid model, the definitions of kinetic and magnetic energy are identical and equal to

$$\begin{equation} E_{{\rm kin}}=\frac{1}{2}\int \rho |\boldsymbol{u}|^2 dx^3, \quad E_{{\rm mag}}=\frac{1}{2}\int |\boldsymbol{b}|^2 dx^3. \end{equation} \tag{ 17 }$$

However, the definition of internal energy is naturally different in each model. In the Hall-MHD model, the internal energy is defined as

$$\begin{equation} \hbox{HMHD:} \quad E_{{\rm int}}=\frac{\beta _o}{\gamma (\gamma -1)}\int \rho ^\gamma dx^3, \end{equation} \tag{ 18 }$$

whereas in the FLR-Landau fluid model with isothermal electrons, the internal energy is given by

$$\begin{equation} \hbox{Landau:} \quad E_{{\rm int}} = \frac{\beta _\parallel }{2}\int \left(p_\perp +\frac{p_\parallel }{2} +T_e^{(0)}\rho \ln \rho \right) dx^3. \end{equation} \tag{ 19 }$$

It is emphasized that because of the filtering, the total energy E_tot is not exactly preserved in the Hall-MHD and Landau fluid simulations. During the simulations with β₀ = β_∥ = 0.8, the total energy in Hall-MHD decreased by approximately 1.2%, whereas in FLR-Landau fluid the decrease was approximately 0.5%. Filtering indeed dissipates kinetic and magnetic energies without turning them into heat. Therefore, in contrast with Landau fluid models, Hall-MHD simulations do not have heating at all and the internal energy could increase only through development of density fluctuations. Time evolution of the sum of kinetic and magnetic energies (normalized to their initial values) is displayed in Figure 8 (left), where the energy contained in the ambient magnetic field was subtracted. This figure shows that this sum decays faster for FLR-Landau fluid simulation. Time evolution of internal energy E_int is shown in Figure 8 (right), where energies were again normalized to their initial value. The initial jump observed in both simulations reflects a rapid adjustment from initial conditions that are not close to an equilibrium state. Later on, the internal energy of Hall-MHD simulation decreases, whereas the internal energy of FLR-Landau fluid simulation more or less smoothly increases until around time t = 1500. During this time (which corresponds to over 10⁴ time steps) Landau damping acts strongly and, by mainly damping slow waves, converts the mechanical energy E_kin + E_mag into the internal energy E_int, which represents heating of the plasma. The question also arises what fraction of mechanical energy is dissipated directly by Landau damping and what fraction is dissipated by the filtering process. Unfortunately, we are unaware of any technique how to address this question.

Zoom In Zoom Out Reset image size

We note that for simulations of freely decaying turbulence the heating is quite weak, implying that driving the system is necessary to produce significant temperature anisotropies. However, the absence of forcing, which yields only a small amount of heating, makes it easier to precisely identify various waves in the system as was presented in Section 3.

The presence of Landau damping can also be seen in the usual wavenumber velocity and magnetic field spectra. Considering first simulations with β₀ = β_∥ = 0.8, Figure 9 shows the velocity spectra with respect to perpendicular and parallel wavenumbers k_⊥, k_∥ which are defined as E^U = ∫E^u(k_⊥)dk_⊥ = ∫E^u(k_∥)dk_∥. With respect to k_⊥, the spectra for Hall-MHD and FLR-Landau fluid are almost identical (Figure 9, left) whereas with respect to k_∥ (Figure 9, right), the spectra of Landau fluid are much steeper. Landau damping therefore significantly inhibits the parallel transfer. Even though low resolution does not allow to precisely identify the slopes of the spectra, three straight lines were added to figures and correspond to power-law solutions k^s, where s = −3/2, −5/3, and −7/3. For E^u(k_⊥), the closest spectral index value appears to be −5/3, the spectral range being however quite limited. The same conclusion with the inhibition of parallel transfer is also obtained for the magnetic field spectra, which are almost identical to the velocity spectra and are shown in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast, for simulations with β₀ = β_∥ = 0.1 the parallel spectrum E(k_∥) displays a similar behavior for Hall-MHD and Landau fluids. The velocity spectra are shown in Figure 11 and the magnetic field spectra are shown in Figure 12. This is consistent with results presented in the previous section where it was shown that Landau damping is responsible for significant reduction of compressibility for simulations with β₀ = β_∥ = 0.8, whereas for simulations with β₀ = β_∥ = 0.1, Landau damping was much weaker and the reduction of compressibility almost negligible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is useful to compare our compressible simulations to incompressible simulations. Naturally, our compressible Hall-MHD code cannot be run in an incompressible regime, for which the turbulent sonic Mach number M_s → 0 and the sound speed c_s → ∞. Nevertheless, incompressible MHD simulations of decaying turbulence were performed, for example, by Bigot et al. (2008). These simulations showed that the combined velocity and magnetic field spectra (they used Elsässer variable $\boldsymbol{z}^+$) are much steeper with respect to k_∥ than with respect to k_⊥, if the ambient magnetic field is sufficiently strong. Our compressible simulations presented here have initially $\langle \boldsymbol{b}^2\rangle ^{1/2}/B_0=1/8$ and can therefore be compared to incompressible MHD simulations of Bigot et al. (2008) (their Figures 17(c) and (d)), where this ratio is 1/5 and 1/15, respectively. Considering compressible Hall-MHD, our simulations showed that spectra with respect to k_∥ (e.g., Figure 9, right, red line) are steeper than spectra with respect to k_⊥ (Figure 9, left, red line). However, these parallel spectra are nowhere near as steep as the parallel spectra of Bigot et al. (2008). Interestingly, their parallel spectra more resemble our parallel spectra for Landau fluid model (Figure 9, right, blue line), where Landau damping was strong (β_∥ = 0.8). These results show that magnetosonic waves (which are not present in an incompressible MHD description) play an important role in regulating the parallel energy cascade. Steeper parallel spectra in Landau fluid model are therefore a consequence of damping of slow magnetosonic waves by Landau damping.

We have presented the first three-dimensional fluid simulations of decaying turbulence in a collisionless plasma in conditions close to solar wind. For this purpose, we used the FLR-Landau fluid model that extends compressible Hall-MHD by incorporating low-frequency kinetic effects such as Landau damping and FLR corrections. It was shown that in spite of the turbulent regime, it is possible to precisely identify linear waves present in the system. Comparisons between compressible Hall-MHD and FLR-Landau fluid model showed that when beta is not too small, linear Landau damping yields strong damping of slow magnetosonic waves in Landau fluid simulations. These waves are indeed damped in kinetic theory described by the Vlasov–Maxwell equations but not in compressible MHD and Hall-MHD descriptions, which overestimate compressibility and parallel transfer in modeling weakly collisional plasmas. The FLR-Landau fluid model can therefore be useful for simulating solar wind, which is typically found to be only weakly compressible.

The support of INSU-CNRS "Programme Soleil-Terre" is acknowledged. Computations were performed on the Mesocentre SIGAMM machine hosted by the Observatoire de la Côte d'Azur (OCA) and on the JADE cluster of the CINES computational facilities. P.H. was supported by an OCA Poincaré fellowship and by the NASA Postdoctoral Program. The work of D.B. was supported by the European Community under the contract of Association between EURATOM and ENEA. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

To clearly understand how Landau damping acts in the present system (Equations (1)–(12)), it is useful to solve dispersion relations for linear waves propagating in parallel direction to the ambient magnetic field. A detailed analysis of linear waves for various propagation angles was elaborated by Passot & Sulem (2004). To simplify the analytic expressions, we define the proton temperature anisotropy as T⁽⁰⁾_⊥p/T_∥p⁽⁰⁾ ≡ a_p, and the normalized electron temperature as T⁽⁰⁾_e/T_∥p⁽⁰⁾ ≡ τ. It can be shown that for parallel propagation angle, the Landau fluid model contains four dispersive Alfvén waves. Two waves obey the dispersion relation which can be expressed as

$$\begin{eqnarray} \!\!\!\!\!\!\!\omega &=& \pm \frac{k^2}{2R_i} \left[ 1+\beta _\parallel \left(1-\frac{a_p}{2}\right) \right]\quad\nonumber\\ &&+ k\sqrt{ 1+ \frac{\beta _\parallel }{2}(a_p-1) + \left(\frac{k}{2R_i} \right)^2 \left[ 1-\beta _\parallel \left(1- \frac{a_p}{2}\right) \right]^2 },\nonumber\\ \end{eqnarray} \tag{ A1 }$$

with another two solutions obtained by substituting ω with −ω. Obviously, these Alfvén waves are independent of the electron temperature τ, which is a consequence of electrons being modeled as isothermal. For a more general Landau fluid model which contains evolution equation for electron pressures and heat fluxes, the electron temperature τ enters the dispersion relation for Alfvén waves. Solutions (A1) can become imaginary, if the expression under the square root becomes negative. At large scales (when 1/R_i → 0) the condition 1 + β_∥(a_p − 1)/2 < 0 represents the well-known criterion for fire-hose instability (see, for example, Ferrière & André 2002). The Hall term and FLR corrections modify the instability criterion. For isotropic temperatures (a_p = 1), solution (A1) naturally collapses to solution (14). The four Alfvén waves (A1) can be eliminated from the general dispersion relation and this yields a complex polynomial of sixth order in frequency ω. Solutions of this polynomial represent three forward and three backward propagating waves which have a negative imaginary part and are therefore damped. Importantly, it is possible to eliminate the dependence on β_∥ and wavenumber k and, after applying a substitution $\Omega =\omega /(k\sqrt{\beta _\parallel })$, the polynomial of sixth order can be simplified to

$$\begin{eqnarray} \Omega ^6 &+& \Omega ^5 i \sqrt{\pi }\frac{-4+3\pi }{-16+6\pi } +\Omega ^4 \frac{-\pi (14+3\tau)+24+8\tau }{-16+6\pi }\nonumber\\ \ms &+& \Omega ^3 i \sqrt{\pi }\frac{\tau -\pi (9+3\tau)/4}{-8+3\pi }+\Omega ^2 \frac{\pi (15+5\tau)/4-2}{-8+3\pi } \nonumber \\ \ms &+& \Omega ^1 i \sqrt{\pi }\frac{(5+3\tau)/2}{-8+3\pi } - \frac{1+\tau }{-8+3\pi } =0. \end{eqnarray} \tag{ A2 }$$

Because this polynomial in Ω does not depend on β_∥ or k, the substitution implies that all six waves are linear with k and $\sqrt{\beta _\parallel }$. The polynomial (A2) has to be solved numerically for a given value of τ. The simulations presented here use τ = 1 and numerically solving polynomial (A2) yields Ω = ±1.48106 − i0.36117, Ω = ±0.65467 − i0.88285, and Ω = ±0.55104 − i0.44311. These six waves therefore satisfy

$$\begin{eqnarray} && \omega = k\sqrt{\beta _\parallel }(\pm 1.48106 -i 0.36117), \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} && \omega = k\sqrt{\beta _\parallel }(\pm 0.65467 -i 0.88285), \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} && \omega = k\sqrt{\beta _\parallel }(\pm 0.55104 -i 0.44311). \end{eqnarray} \tag{ A5 }$$

The least damped solution (A3) represents the sound wave. The solutions of Landau fluid model (Equations (1)–(12)) for parallel propagation angle are therefore four Alfvén waves (A1), two sound waves (A3), and four waves (A4), (A5), which are highly damped. These four waves (A4), (A5) do not have an analogy in Hall-MHD description and correspond to solutions of kinetic Maxwell–Vlasov description, which contains an infinite number of highly damped solutions. Interestingly, the last solution (A5) is not dependent on the value of τ and it can be expressed analytically as $\Omega =\pm \sqrt{8-\pi }/4 -i\sqrt{\pi }/4$. After eliminating these waves from Equation (A2), the polynomial which contains the sound waves (A3) and solutions (A4) is now of fourth order in Ω and expressed as

$$\begin{eqnarray} \Omega ^4 &+&i\frac{2\sqrt{\pi }}{-8+3\pi }\Omega ^3 -\frac{1}{2}\frac{9\pi -16 +(-8+3\pi)\tau }{-8+3\pi }\Omega ^2\nonumber\\ \ms &&-i\frac{\sqrt{\pi }(3+\tau)}{-8+3\pi }\Omega +\frac{2(1+\tau)}{-8+3\pi }=0. \end{eqnarray} \tag{ A6 }$$

It is of course possible to use Ferrari–Cardano's relations to solve this polynomial analytically, the final result is however too complicated and it is still more convenient to solve Equation (A6) numerically for a given value of τ.

If this wave analysis is repeated with the model (Equations (1)–(10)) with heat flux equations q_∥ = 0, q_⊥ = 0, the same dispersion relation (A1) for Alfvén waves is obtained. However, the only other solutions present in the parallel direction are

$$\begin{equation} \omega = \pm k\sqrt{\frac{\beta _\parallel }{2}(3+\tau)}. \end{equation} \tag{ A7 }$$

These waves have frequencies which are purely real and correspond to undamped sound waves of double adiabatic model with isothermal electrons.

For completeness, considering perpendicular propagation, it can be shown that the heat fluxes vanish and the solutions of the Landau fluid model with isothermal electrons are undamped magnetosonic waves with the dispersion relation expressed as

$$\begin{equation} \omega =\pm k \sqrt{1+\beta _\parallel \left(a_p+\frac{\tau }{2}\right) + \left(\frac{a_p\beta _\parallel k}{4 R_i} \right)^2 }. \end{equation} \tag{ A8 }$$