2.1 Parent gyrofluid model

As starting point for the derivation of the model we consider the following gyrofluid system (Scott Reference Scott2010), in normalized units, consisting of the evolution equations

(2.1) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}N_{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},N_{\unicode[STIX]{x1D6FC}}]-[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert },U_{\unicode[STIX]{x1D6FC}}]+[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},T_{\bot _{\unicode[STIX]{x1D6FC}}}]-[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}A_{\Vert },Q_{\bot _{\unicode[STIX]{x1D6FC}}}]+{\displaystyle \frac{\unicode[STIX]{x2202}U_{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x2202}z}}=0,\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}\left({\displaystyle \frac{m_{\unicode[STIX]{x1D6FC}}}{m_{i}}}U_{\unicode[STIX]{x1D6FC}}+\operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert }\right)+\left[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},{\displaystyle \frac{m_{\unicode[STIX]{x1D6FC}}}{m_{i}}}U_{\unicode[STIX]{x1D6FC}}+\operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert }\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}A_{\Vert }]+{\displaystyle \frac{m_{\unicode[STIX]{x1D6FC}}}{m_{i}}}[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},Q_{\bot _{\unicode[STIX]{x1D6FC}}}]-{\displaystyle \frac{T_{\unicode[STIX]{x1D6FC}}}{T_{e}}}[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert },N_{\unicode[STIX]{x1D6FC}}]-{\displaystyle \frac{T_{\unicode[STIX]{x1D6FC}}}{T_{e}}}[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert },T_{\Vert _{\unicode[STIX]{x1D6FC}}}]\nonumber\\ \displaystyle & & \displaystyle \quad -\,{\displaystyle \frac{T_{\unicode[STIX]{x1D6FC}}}{T_{e}}}[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}A_{\Vert },T_{\bot _{\unicode[STIX]{x1D6FC}}}]+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}}\left(\operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711}+{\displaystyle \frac{T_{\unicode[STIX]{x1D6FC}}}{T_{e}}}N_{\unicode[STIX]{x1D6FC}}+{\displaystyle \frac{T_{\unicode[STIX]{x1D6FC}}}{T_{e}}}T_{\Vert _{\unicode[STIX]{x1D6FC}}}\right)=0,\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\Vert _{\unicode[STIX]{x1D6FC}}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},T_{\Vert _{\unicode[STIX]{x1D6FC}}}]-2[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert },U_{\unicode[STIX]{x1D6FC}}]-[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert },Q_{\Vert _{\unicode[STIX]{x1D6FC}}}]\nonumber\\ \displaystyle & & \displaystyle \quad -\,2[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}A_{\Vert },Q_{\bot _{\unicode[STIX]{x1D6FC}}}]+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}}(2U_{\unicode[STIX]{x1D6FC}}+Q_{\Vert _{\unicode[STIX]{x1D6FC}}})=0,\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\bot _{\unicode[STIX]{x1D6FC}}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},T_{\bot _{\unicode[STIX]{x1D6FC}}}]+[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},N_{\unicode[STIX]{x1D6FC}}]+2[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},T_{\bot _{\unicode[STIX]{x1D6FC}}}]-[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert },Q_{\bot _{\unicode[STIX]{x1D6FC}}}]-[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}A_{\Vert },U_{\unicode[STIX]{x1D6FC}}]\nonumber\\ \displaystyle & & \displaystyle \quad -\,2[\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}A_{\Vert },Q_{\bot _{\unicode[STIX]{x1D6FC}}}]+{\displaystyle \frac{\unicode[STIX]{x2202}Q_{\bot _{\unicode[STIX]{x1D6FC}}}}{\unicode[STIX]{x2202}z}}=0,\end{eqnarray}$$

completed by the quasi-neutrality relation

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FC}}\left(\operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})(\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}N_{\unicode[STIX]{x1D6FC}}+\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}T_{\bot _{\unicode[STIX]{x1D6FC}}})+{\displaystyle \frac{T_{e}}{T_{\unicode[STIX]{x1D6FC}}}}(\unicode[STIX]{x1D6E4}_{0_{\unicode[STIX]{x1D6FC}}}-1)\unicode[STIX]{x1D711}\right)=0, & \displaystyle\end{eqnarray}$$

and by the parallel Ampère’s law

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }=-{\displaystyle \frac{\unicode[STIX]{x1D6FD}_{e}}{2}}\mathop{\sum }_{\unicode[STIX]{x1D6FC}}\operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})(\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}U_{\unicode[STIX]{x1D6FC}}+\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}Q_{\bot _{\unicode[STIX]{x1D6FC}}}). & \displaystyle\end{eqnarray}$$

The above system is formulated on a domain ${\mathcal{D}}=\{(x,y,z):-L_{x}\leqslant x\leqslant L_{x},-L_{y}\leqslant y\leqslant L_{y},-L_{z}\leqslant z\leqslant L_{z}\}$ with $L_{x},L_{y},L_{z}$ corresponding to normalized lengths. We assume the presence of a strong and uniform component of the magnetic field along $\hat{\boldsymbol{z}}$ . The transverse component of the magnetic field is indicated with $\unicode[STIX]{x1D735}A_{\Vert }\times \hat{\boldsymbol{z}}$ , with $A_{\Vert }$ denoting the fluctuations of the magnetic flux function, which are assumed to be small compared to the amplitude of the guide field. Fluctuations of the magnetic field along $\hat{\boldsymbol{z}}$ , moreover, are neglected. The symbol $\unicode[STIX]{x1D6FB}_{\bot }^{2}\equiv \unicode[STIX]{x2202}_{xx}+\unicode[STIX]{x2202}_{yy}$ , also, indicates the perpendicular Laplacian operator, where derivatives are taken only with respect to coordinates describing the plane perpendicular to the guide field. Electrostatic potential fluctuations are indicated with $\unicode[STIX]{x1D711}$ . The index $\unicode[STIX]{x1D6FC}$ denotes the particle species, with $\unicode[STIX]{x1D6FC}=e$ and $\unicode[STIX]{x1D6FC}=i$ referring to electrons and ions, respectively. The fields $N_{\unicode[STIX]{x1D6FC}}$ , $U_{\unicode[STIX]{x1D6FC}}$ , $T_{\Vert _{\unicode[STIX]{x1D6FC}}}$ , $T_{\bot _{\unicode[STIX]{x1D6FC}}}$ , $Q_{\Vert _{\unicode[STIX]{x1D6FC}}}$ , $Q_{\bot _{\unicode[STIX]{x1D6FC}}}$ are all defined on ${\mathcal{D}}$ and denote the fluctuations of the gyrocentre density, fluid velocity along the magnetic field, parallel and perpendicular temperature, parallel and perpendicular heat flux, respectively. Note that, with respect to Scott (Reference Scott2010), we adopted a different definition for the parallel heat flux fluctuations $Q_{\Vert _{\unicode[STIX]{x1D6FC}}}$ which corresponds to the heat flux fluctuations of Scott (Reference Scott2010) multiplied by a factor $2$ . The constants $q_{\unicode[STIX]{x1D6FC}}$ , $m_{\unicode[STIX]{x1D6FC}}$ and $T_{\unicode[STIX]{x1D6FC}}$ , on the other hand, indicate the charge, the mass of the particle and the equilibrium temperature of the species $\unicode[STIX]{x1D6FC}$ , respectively. We indicated with $\unicode[STIX]{x1D6E4}_{0_{\unicode[STIX]{x1D6FC}}}$ , $\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}$ and $\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}$ the gyroaverage operators, relative to the species $\unicode[STIX]{x1D6FC}$ , which commonly appear in gyrofluid models (Brizard Reference Brizard1992; Dorland & Hammett Reference Dorland and Hammett1993; Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993; Beer & Hammett Reference Beer and Hammett1996; Scott Reference Scott2010). In Fourier space, assuming periodic boundary conditions along the $y$ and $z$ directions, these are defined, in terms of dimensional quantities, as

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}_{0_{\unicode[STIX]{x1D6FC}}}=I_{0}\left(b_{\unicode[STIX]{x1D6FC}}\right)\exp \left(-b_{\unicode[STIX]{x1D6FC}}\right), & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{1}{n_{0}}}\displaystyle \int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,{\mathcal{F}}_{eq_{\unicode[STIX]{x1D6FC}}}J_{0}\left(\hat{k}_{\bot }\sqrt{2\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}m_{\unicode[STIX]{x1D6FC}}c^{2}/(|q_{\unicode[STIX]{x1D6FC}}|^{2}B_{0})}\right), & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{1}{n_{0}}}\displaystyle \int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,\left({\displaystyle \frac{\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}B_{0}}{T_{\unicode[STIX]{x1D6FC}}}}-1\right){\mathcal{F}}_{eq_{\unicode[STIX]{x1D6FC}}}J_{0}\left(\hat{k}_{\bot }\sqrt{2\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}m_{\unicode[STIX]{x1D6FC}}c^{2}/(|q_{\unicode[STIX]{x1D6FC}}|^{2}B_{0})}\right), & \displaystyle\end{eqnarray}$$

where $b_{\unicode[STIX]{x1D6FC}}=\hat{k}_{\bot }^{2}\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}^{2}$ , with $\hat{k}_{\bot }=(\hat{k}_{x}^{2}+\hat{k}_{y}^{2})^{1/2}$ corresponding to the dimensional wavenumber in the direction perpendicular to the guide field. The wavenumbers $\hat{k}_{x}$ and $\hat{k}_{y}$ are defined as $\hat{k}_{x}=2\unicode[STIX]{x03C0}m/(2\hat{L}_{x})$ and $\hat{k}_{y}=2\unicode[STIX]{x03C0}n/(2\hat{L}_{y})$ , respectively, with $m\in \mathbb{Z}$ , $n\in \mathbb{Z}$ and with $\hat{L}_{x}=\unicode[STIX]{x1D70C}_{s}L_{x}$ , $\hat{L}_{y}=\unicode[STIX]{x1D70C}_{s}L_{y}$ corresponding to the lengths $L_{x}$ and $L_{y}$ without normalization. The sonic Larmor radius $\unicode[STIX]{x1D70C}_{s}$ is defined in (2.13). The symbol $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}$ is defined by $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}=(T_{\unicode[STIX]{x1D6FC}}/m_{\unicode[STIX]{x1D6FC}})^{1/2}/(|q_{\unicode[STIX]{x1D6FC}}|B_{0}/(m_{\unicode[STIX]{x1D6FC}}c))$ and corresponds to the thermal gyroradius characteristic of the species $\unicode[STIX]{x1D6FC}$ , with $c$ equal to the speed of light. In equations (2.7)–(2.9), $I_{0}$ and $J_{0}$ indicate the zeroth-order modified Bessel function and Bessel function of the first kind, respectively. We also denoted by $n_{0}$ a uniform equilibrium gyrocentre density, with $B_{0}$ the amplitude of the magnetic guide field, with $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}$ the magnetic moment for particles of species $\unicode[STIX]{x1D6FC}$ and with $\text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}=(2\unicode[STIX]{x03C0}B_{0}/m_{\unicode[STIX]{x1D6FC}})\,\text{d}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}\,\text{d}v$ the volume element in velocity space, with $v$ indicating the velocity coordinate along the magnetic guide field. The physical meaning of the constant parameter $\unicode[STIX]{x1D6FD}_{e}$ was already anticipated in § 1 and its explicit expression is given by $\unicode[STIX]{x1D6FD}_{e}=8\unicode[STIX]{x03C0}n_{0}T_{e}/B_{0}^{2}$ . As is common in the so-called $\unicode[STIX]{x1D6FF}f$ approach, the gyrocentre distribution functions are assumed to be close to equilibrium Maxwellian functions ${\mathcal{F}}_{eq_{\unicode[STIX]{x1D6FC}}}$ defined as

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{F}}_{eq_{\unicode[STIX]{x1D6FC}}}(v,\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}})=n_{0}\left({\displaystyle \frac{m_{\unicode[STIX]{x1D6FC}}}{2\unicode[STIX]{x03C0}T_{\unicode[STIX]{x1D6FC}}}}\right)^{3/2}\exp \left(-{\displaystyle \frac{m_{\unicode[STIX]{x1D6FC}}v^{2}}{2T_{\unicode[STIX]{x1D6FC}}}}-{\displaystyle \frac{\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}B_{0}}{T_{\unicode[STIX]{x1D6FC}}}}\right). & \displaystyle\end{eqnarray}$$

As indicated in Dorland & Hammett (Reference Dorland and Hammett1993), Beer & Hammett (Reference Beer and Hammett1996) and Scott (Reference Scott2010), approximated expressions for the gyroaverage operators $\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}$ and $\unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}$ are given by

(2.11a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}=\unicode[STIX]{x1D6E4}_{0_{\unicode[STIX]{x1D6FC}}}^{1/2},\quad \unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}=b_{\unicode[STIX]{x1D6FC}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}}{\unicode[STIX]{x2202}b_{\unicode[STIX]{x1D6FC}}}}. & & \displaystyle\end{eqnarray}$$

The system (2.1)–(2.4) is formulated making use of normalized variables, which are defined by

(2.12) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}t=\unicode[STIX]{x1D6FA}_{i}\hat{t},\quad x={\displaystyle \frac{\hat{x}}{\unicode[STIX]{x1D70C}_{s}}},\quad y={\displaystyle \frac{{\hat{y}}}{\unicode[STIX]{x1D70C}_{s}}},\quad z={\displaystyle \frac{\hat{z}}{\unicode[STIX]{x1D70C}_{s}}},\\[12.0pt] N_{\unicode[STIX]{x1D6FC}}={\displaystyle \frac{\hat{N}_{\unicode[STIX]{x1D6FC}}}{n_{0}}},\quad U_{\unicode[STIX]{x1D6FC}}={\displaystyle \frac{\hat{U} _{\unicode[STIX]{x1D6FC}}}{c_{s}}},\quad T_{\Vert _{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{\hat{T}_{\Vert _{\unicode[STIX]{x1D6FC}}}}{n_{0}T_{\unicode[STIX]{x1D6FC}}}},\quad T_{\bot _{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{\hat{T}_{\bot _{\unicode[STIX]{x1D6FC}}}}{n_{0}T_{\unicode[STIX]{x1D6FC}}}},\\[12.0pt] Q_{\Vert _{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{\hat{Q}_{\Vert _{\unicode[STIX]{x1D6FC}}}}{n_{0}T_{\unicode[STIX]{x1D6FC}}c_{s}}},\quad Q_{\bot _{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{\hat{Q}_{\bot _{\unicode[STIX]{x1D6FC}}}}{n_{0}T_{\unicode[STIX]{x1D6FC}}c_{s}}},\\[12.0pt] \unicode[STIX]{x1D711}={\displaystyle \frac{e\hat{\unicode[STIX]{x1D711}}}{T_{e}}},\quad A_{\Vert }={\displaystyle \frac{\hat{A}_{\Vert }}{B_{0}\unicode[STIX]{x1D70C}_{s}}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

In (2.12) we introduced the proton charge $e=|q_{e}|$ , we denoted with a caret the dimensional quantities and we introduced the following constants:

(2.13a-c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{i}={\displaystyle \frac{eB_{0}}{m_{i}c}},\quad c_{s}=\sqrt{{\displaystyle \frac{T_{e}}{m_{i}}}},\quad \unicode[STIX]{x1D70C}_{s}={\displaystyle \frac{c_{s}}{\unicode[STIX]{x1D6FA}_{i}}}, & \displaystyle\end{eqnarray}$$

which correspond to the ion cyclotron frequency (referred to the guide field), the sound speed and the sonic Larmor radius, respectively.

In order to make contact also with the underlying $\unicode[STIX]{x1D6FF}f$ gyrokinetic theory, we specify that the dimensional moment fluctuations are defined as

(2.14a,b ) $$\begin{eqnarray}\displaystyle \hat{N}_{\unicode[STIX]{x1D6FC}}=\int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,\hat{f}_{\unicode[STIX]{x1D6FC}},\quad \hat{U} _{\unicode[STIX]{x1D6FC}}={\displaystyle \frac{1}{n_{0}}}\int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,v\hat{f}_{\unicode[STIX]{x1D6FC}}, & & \displaystyle\end{eqnarray}$$

(2.15a,b ) $$\begin{eqnarray}\displaystyle \hat{T}_{\Vert _{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{1}{n_{0}}}\int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,(m_{\unicode[STIX]{x1D6FC}}v^{2}-T_{\unicode[STIX]{x1D6FC}})\hat{f}_{\unicode[STIX]{x1D6FC}},\quad \hat{T}_{\bot _{\unicode[STIX]{x1D6FC}}}={\displaystyle \frac{1}{n_{0}}}\int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,(\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}B_{0}-T_{\unicode[STIX]{x1D6FC}})\hat{f}_{\unicode[STIX]{x1D6FC}}, & & \displaystyle\end{eqnarray}$$

(2.16a,b ) $$\begin{eqnarray}\displaystyle \hat{Q}_{\Vert _{\unicode[STIX]{x1D6FC}}}=\int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,(m_{\unicode[STIX]{x1D6FC}}v^{2}-3T_{\unicode[STIX]{x1D6FC}})v\hat{f}_{\unicode[STIX]{x1D6FC}},\quad \hat{Q}_{\bot _{\unicode[STIX]{x1D6FC}}}=\int \text{d}{\mathcal{W}}_{\unicode[STIX]{x1D6FC}}\,(\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}B_{0}-T_{\unicode[STIX]{x1D6FC}})v\hat{f}_{\unicode[STIX]{x1D6FC}}, & & \displaystyle\end{eqnarray}$$

where $\hat{f}_{\unicode[STIX]{x1D6FC}}=\hat{f}_{\unicode[STIX]{x1D6FC}}(x,y,z,v,\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}},t)$ is the perturbation of the gyrocentre distribution function related to the species $\unicode[STIX]{x1D6FC}$ .

The parent model (2.1)–(2.6) can be seen as a simplified version of the gyrofluid model (although written with a different normalization) of Scott (Reference Scott2010), taken in the non-dissipative limit and neglecting background inhomogeneities as well as evolution equations for the heat fluxes. The system (2.1)–(2.6), in particular, is evidently not closed, for closure relations for the heat flux fluctuations were not provided at this stage. Whereas the model of Scott (Reference Scott2010) includes evolution equations for the heat flux fluctuations, the system (2.1)–(2.6) does not. In order to obtain a simple model although retaining the effects of Landau damping, we decided to truncate the gyrofluid hierarchy at the lowest level where it is of interest to introduce Landau closures, that is at the level of temperature evolution equations.

As above anticipated, the choice of a simplified version of the model of Scott (Reference Scott2010) as parent model is motivated by the fact that its gyrofluid closures are determined in such a way that the model also conserves energy in the absence of heat flux fluctuations. Also, the parent model retains electron inertia, which we require for driving magnetic reconnection. In order to facilitate the comparison with Scott (Reference Scott2010) we specify that the parent model consisting of (2.1)–(2.4), (2.5) and (2.6) can be obtained from equations (163)–(166), (132) and (133) of Scott (Reference Scott2010) (where dissipation and background inhomogeneities are already omitted and when also the two-dimensional approximation for the fluctuations is assumed) by making the following formal replacements

(2.17a-c ) $$\begin{eqnarray}\displaystyle n_{z}\rightarrow N_{\unicode[STIX]{x1D6FC}},\quad u_{z\Vert }\rightarrow U_{\unicode[STIX]{x1D6FC}},\quad T_{z\Vert }\rightarrow T_{\Vert _{\unicode[STIX]{x1D6FC}}}, & & \displaystyle\end{eqnarray}$$

(2.18a-c ) $$\begin{eqnarray}\displaystyle T_{z\bot }\rightarrow T_{\bot _{\unicode[STIX]{x1D6FC}}},\quad q_{z\Vert }\rightarrow (1/2)Q_{\Vert _{\unicode[STIX]{x1D6FC}}},\quad q_{z\bot }\rightarrow Q_{\bot _{\unicode[STIX]{x1D6FC}}}, & & \displaystyle\end{eqnarray}$$

(2.19a-c ) $$\begin{eqnarray}\displaystyle A_{\Vert }\rightarrow A_{\Vert }/\unicode[STIX]{x1D6FD}_{e},\quad \unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D711},\quad \unicode[STIX]{x1D6FD}_{e}A_{G}\rightarrow \unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert }, & & \displaystyle\end{eqnarray}$$

(2.20a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{G}\rightarrow \unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711},\quad \unicode[STIX]{x1D6FD}_{e}\unicode[STIX]{x1D712}_{G}\rightarrow \unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}A_{\Vert },\quad \unicode[STIX]{x1D6FA}_{G}\rightarrow \unicode[STIX]{x1D6E4}_{2_{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D711}, & & \displaystyle\end{eqnarray}$$

(2.21a-c ) $$\begin{eqnarray}\displaystyle p_{z\Vert }\rightarrow N_{\unicode[STIX]{x1D6FC}}+T_{\Vert _{\unicode[STIX]{x1D6FC}}},\quad \unicode[STIX]{x1D6FD}_{e}\rightarrow \unicode[STIX]{x1D6FD}_{e}/2,\quad \unicode[STIX]{x1D707}_{z}\rightarrow \operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})m_{\unicode[STIX]{x1D6FC}}/m_{i}, & & \displaystyle\end{eqnarray}$$

(2.22a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{z}\rightarrow \operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}})T_{\unicode[STIX]{x1D6FC}}/T_{e},\quad a_{z}\rightarrow \operatorname{sgn}(q_{\unicode[STIX]{x1D6FC}}),\quad \unicode[STIX]{x1D6E4}_{i}\rightarrow \unicode[STIX]{x1D6E4}_{i_{\unicode[STIX]{x1D6FC}}}\quad \text{for}~i=0,1,2. & & \displaystyle\end{eqnarray}$$

Also, because equations (163)–(166), (132) and (133) of Scott (Reference Scott2010) refer to a two-dimensional limit, in order to retrieve equations (2.1)–(2.4) of the present paper one needs to add terms accounting for the dependence on the $z$ coordinate. This is accomplished by replacing the operator $\unicode[STIX]{x1D6FD}_{e}[A_{G},\,.]$ with the operator $[\unicode[STIX]{x1D6E4}_{1_{\unicode[STIX]{x1D6FC}}}A_{\Vert },\,.]-\unicode[STIX]{x2202}_{z}$ .

2.2 Derivation of the model

We consider a plasma consisting of electrons and of one single ionized ion species. The derivation of the reduced Landau-gyrofluid model proceeds essentially by an asymptotic expansion based on a prescribed scaling, and on imposing Landau-fluid closures on the heat fluxes. The asymptotic expansion is based on two small, independent parameters $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D716}$ defined by

(2.23a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}^{2}={\displaystyle \frac{m_{e}}{m_{i}}},\quad \unicode[STIX]{x1D716}=\frac{1}{\unicode[STIX]{x1D6FA}_{i}{\mathcal{T}}}, & \displaystyle\end{eqnarray}$$

with ${\mathcal{T}}$ indicating the characteristic time scale of variation of the fluctuating fields $N_{e,i}$ , $U_{e,i}$ , $T_{\Vert _{e,i}}$ and $T_{\bot _{e,i}}$ (which is assumed to be the same for all such fields). We also introduce the parameter

(2.24) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}={\displaystyle \frac{T_{i}}{T_{e}}}, & \displaystyle\end{eqnarray}$$

expressing the equilibrium temperature ratio.

The following ordering is assumed:

(2.25a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}^{2}\ll 1,\quad \unicode[STIX]{x1D70F}\sim k_{\bot }=O(1), & \displaystyle\end{eqnarray}$$

(2.26) $$\begin{eqnarray}\displaystyle & \displaystyle U_{e}\sim Q_{\Vert _{e}}\sim Q_{\bot _{e}}=O\left({\displaystyle \frac{\unicode[STIX]{x1D716}}{\unicode[STIX]{x1D6FF}}}\right), & \displaystyle\end{eqnarray}$$

(2.27) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}\sim N_{e}\sim T_{\Vert _{e}}\sim T_{\bot _{e}}\sim N_{i}\sim T_{\Vert _{i}}\sim T_{\bot _{i}}\sim \unicode[STIX]{x1D711}=O(\unicode[STIX]{x1D716})\ll 1, & \displaystyle\end{eqnarray}$$

(2.28) $$\begin{eqnarray}\displaystyle & \displaystyle k_{z}\sim U_{i}\sim Q_{\Vert _{i}}\sim Q_{\bot _{i}}\sim A_{\Vert }=O(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D716}), & \displaystyle\end{eqnarray}$$

(2.29) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}_{e}=O(\unicode[STIX]{x1D6FF}^{2}). & \displaystyle\end{eqnarray}$$

Assuming $\unicode[STIX]{x1D6FF}^{2}\ll 1$ implies $b_{e}\ll 1$ , and consequently $\unicode[STIX]{x1D6E4}_{0e}\simeq 1,\unicode[STIX]{x1D6E4}_{1e}\simeq 1,\unicode[STIX]{x1D6E4}_{2e}\simeq 0$ , which allows us to neglect finite electron Larmor radius effects due to the gyroaverage operators. On the other hand, the ordering on $A_{\Vert },\unicode[STIX]{x1D6FD}_{e},k_{\bot }$ and $U_{e}$ retains finite electron inertia in the electron momentum equation, thus allowing for reconnection, and balances parallel current density with the parallel electron velocity in Ampère’s law. In particular, as anticipated in § 1, together with the requirement of neglecting electron Larmor radius effects, this leads to $\unicode[STIX]{x1D6FD}_{e}=O(\unicode[STIX]{x1D6FF}^{2})\ll 1$  (Zocco & Schekochihin Reference Zocco and Schekochihin2011). The ordering on $N_{e,i}$ , $k_{\bot }$ , $\unicode[STIX]{x1D711}$ , and $\unicode[STIX]{x1D714}$ reflects the ordering of the original parent model. In particular, it is consistent with the requirement that the $\boldsymbol{E}\times \boldsymbol{B}$ drift be the dominant drift on the characteristic time scale ( $\unicode[STIX]{x2202}_{t}\sim [\unicode[STIX]{x1D711},\,.]$ ) and that small amplitude electron and ion gyrocentre density fluctuations be of the same order as polarization terms at the characteristic perpendicular scale ( $N_{i}\sim N_{e}\sim (\unicode[STIX]{x1D6E4}_{0i}-1)(\unicode[STIX]{x1D711}/\unicode[STIX]{x1D70F})$ ) in the quasi-neutrality relation. We also remark that, as a consequence of the ordering, the dimensional electromagnetic potentials satisfy $\hat{\unicode[STIX]{x1D711}}\sim (v_{\text{the}}/c)\hat{A}_{\Vert }$ , with $v_{\text{the}}$ indicating the electron thermal speed. The ordering on $k_{z}$ comes from the strong anisotropy assumption ( $k_{z}\ll k_{\bot }$ ) and also from the requirement that the two contributions to the gradient operator along the magnetic field, be of the same order ( $[A_{\Vert },\,.]\sim k_{z}$ ). Investigations of magnetic reconnection (Comisso et al. Reference Comisso, Grasso, Tassi and Waelbroeck2012) performed adopting a model (Waelbroeck & Tassi Reference Waelbroeck and Tassi2012) compatible with such ordering, on the other hand, indicated that, if $U_{i}=0$ at $t=0$ , then the role of the parallel ion gyrocentre velocity is marginal on the reconnection dynamics. Based on this result, we ordered $U_{i}\ll U_{e}$ . In the present model we allow for parallel and perpendicular temperature fluctuations. Because we are primarily interested in the effect of Landau closures on the electron fluid, we order the heat fluxes in such a way that they enter the evolution equations for the electron temperatures associated with leading-order terms. We allow for finite ion temperature effects in the model, thus we set $\unicode[STIX]{x1D70F}\sim 1$ , although the cold-ion limit will be considered in the simulations.

Applying the above ordering to (2.1)–(2.6) and retaining only the lowest-order terms in each equation yields the following reduced system

(2.30) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}N_{e}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},N_{e}]+\unicode[STIX]{x1D735}_{\Vert }U_{e}=0, & \displaystyle\end{eqnarray}$$

(2.31) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}(A_{\Vert }-\unicode[STIX]{x1D6FF}^{2}U_{e})}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},A_{\Vert }-\unicode[STIX]{x1D6FF}^{2}U_{e}]-\unicode[STIX]{x1D735}_{\Vert }(N_{e}+T_{\Vert _{e}})+{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}z}}=0, & \displaystyle\end{eqnarray}$$

(2.32) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\Vert _{e}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},T_{\Vert _{e}}]+\unicode[STIX]{x1D735}_{\Vert }(2U_{e}+Q_{\Vert _{e}})=0, & \displaystyle\end{eqnarray}$$

(2.33) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\bot _{e}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},T_{\bot _{e}}]+\unicode[STIX]{x1D735}_{\Vert }Q_{\bot _{e}}=0, & \displaystyle\end{eqnarray}$$

(2.34) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}N_{i}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D6E4}_{1_{i}}\unicode[STIX]{x1D711},N_{i}]+[\unicode[STIX]{x1D6E4}_{2_{i}}\unicode[STIX]{x1D711},T_{\bot _{i}}]=0, & \displaystyle\end{eqnarray}$$

(2.35) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\bot _{i}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D6E4}_{1_{i}}\unicode[STIX]{x1D711},T_{\bot _{i}}]+2[\unicode[STIX]{x1D6E4}_{2_{i}}\unicode[STIX]{x1D711},T_{\bot _{i}}]+[\unicode[STIX]{x1D6E4}_{2_{i}}\unicode[STIX]{x1D711},N_{i}]=0, & \displaystyle\end{eqnarray}$$

(2.36) $$\begin{eqnarray}\displaystyle & \displaystyle N_{e}-\unicode[STIX]{x1D6E4}_{1_{i}}N_{i}-{\displaystyle \frac{\unicode[STIX]{x1D6E4}_{0_{i}}-1}{\unicode[STIX]{x1D70F}}}\unicode[STIX]{x1D711}-\unicode[STIX]{x1D6E4}_{2_{i}}T_{\bot _{i}}=0, & \displaystyle\end{eqnarray}$$

(2.37) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }={\displaystyle \frac{\unicode[STIX]{x1D6FD}_{e}}{2}}U_{e}, & \displaystyle\end{eqnarray}$$

where we introduced the parallel gradient operator

(2.38) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\Vert }\,f=-[A_{\Vert },f]+{\displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}z}}, & \displaystyle\end{eqnarray}$$

for a function $f$ . Equations (2.30)–(2.33) correspond to the evolution equations for the electron density, parallel momentum, parallel and perpendicular temperature, respectively. Equations (2.34) and (2.35) account for the evolution of ion gyrocentre density and perpendicular temperature, whereas equations (2.36) and (2.37) are the quasi-neutrality relation and parallel Ampère’s law.

It can be noticed that the system does not contain the evolution equations for the parallel ion gyrocentre velocity and for the parallel temperature. As a consequence of the ordering, such equations turn out to be decoupled from the system and thus we disregard them. Should one be interested in determining $U_{i}$ and $T_{\Vert _{i}}$ , the required equations read

(2.39) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}(U_{i}+\unicode[STIX]{x1D6E4}_{1_{i}}A_{\Vert })+[\unicode[STIX]{x1D6E4}_{1_{i}}\unicode[STIX]{x1D711},U_{i}+\unicode[STIX]{x1D6E4}_{1_{i}}A_{\Vert }]+[\unicode[STIX]{x1D6E4}_{2_{i}}\unicode[STIX]{x1D711},\unicode[STIX]{x1D6E4}_{2_{i}}A_{\Vert }]+[\unicode[STIX]{x1D6E4}_{2_{i}}\unicode[STIX]{x1D711},Q_{\bot _{i}}]-\unicode[STIX]{x1D70F}[\unicode[STIX]{x1D6E4}_{1_{i}}A_{\Vert },N_{i}] & \displaystyle \nonumber\\ \displaystyle & \hspace{-24.0pt}-\,\unicode[STIX]{x1D70F}[\unicode[STIX]{x1D6E4}_{1_{i}}A_{\Vert },T_{\Vert _{i}}]-\unicode[STIX]{x1D70F}[\unicode[STIX]{x1D6E4}_{2_{i}}A_{\Vert },T_{\bot _{i}}]+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}}(\unicode[STIX]{x1D6E4}_{1_{i}}\unicode[STIX]{x1D711}+\unicode[STIX]{x1D70F}N_{i}+\unicode[STIX]{x1D70F}T_{\Vert _{i}})=0, & \displaystyle\end{eqnarray}$$

(2.40) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\Vert _{i}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D6E4}_{1_{i}}\unicode[STIX]{x1D711},T_{\Vert _{i}}]=0, & \displaystyle\end{eqnarray}$$

where a closure on $Q_{\bot _{i}}$ has to be determined.

We restrict ourselves, however, to the system (2.30)–(2.37) which we close by imposing the Landau-fluid closures (Hammett & Perkins Reference Hammett and Perkins1990; Hammett et al. Reference Hammett, Dorland and Perkins1992; Snyder, Hammett & Dorland Reference Snyder, Hammett and Dorland1997; Passot & Sulem Reference Passot and Sulem2004; Sulem & Passot Reference Sulem and Passot2015)

(2.41a,b ) $$\begin{eqnarray}\displaystyle Q_{\Vert _{e}}=-2\unicode[STIX]{x1D6FC}{\mathcal{L}}T_{\Vert _{e}},\quad Q_{\bot _{e}}=-\unicode[STIX]{x1D6FC}{\mathcal{L}}T_{\bot _{e}}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{L}}$ holds for the Landau damping operator. Its modelling is discussed below. In (2.41), we furthermore introduced the constant $\unicode[STIX]{x1D6FC}=(1/\unicode[STIX]{x1D6FF})(2/\unicode[STIX]{x03C0})^{1/2}$ , originating from the linear kinetic theory. We point out that, as will be argued in § 4.1, for the electron fluid with our ordering, temperature fluctuations and heat fluxes of the particles and of the gyrocentres coincide. Consequently, closures (2.41) keep the same form when expressed in terms of particle moments.

As mentioned in Snyder et al. (Reference Snyder, Hammett and Dorland1997), a proper calculation of Landau damping would involve at each point of the computational domain and every time step, the evaluation of a convolution integral along the magnetic field line attached to this point, a task which exceeds the capability of the present day computers. A drastic simplification consists in linearizing the operator near a uniform ambient magnetic field, and thus replacing the local magnetic field lines by straight magnetic lines. The Landau damping operator then reduces to the negative of the Hilbert transform in the direction of the ambient magnetic field (defined by the unit vector $\boldsymbol{b}_{0}$ ) and, in this case, corresponds in Fourier space to the multiplicator $\text{i}\,\boldsymbol{b}_{0}\boldsymbol{\cdot }\boldsymbol{k}/|\boldsymbol{b}_{0}\boldsymbol{\cdot }\boldsymbol{k}|$ . In one space dimension, this description is exact in the weakly nonlinear asymptotics (Ott & Sudan Reference Ott and Sudan1969), and remains satisfactory in the fully nonlinear regime (Laveder et al. Reference Laveder, Marradi, Passot and Sulem2011). It however turns out to be inappropriate in several dimensions, as it may lead to spurious instabilities (Passot et al. Reference Passot, Henri, Laveder and Sulem2014). In this case, indeed, the distortion of the magnetic field lines should be taken into account. Various models have been suggested for this purpose. Sharma et al. (Reference Sharma, Hammett, Quataert and Stone2006) expressed the Landau damping operator in physical space in the form $(1/k_{L})\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}=(1/k_{L})\unicode[STIX]{x1D735}_{\Vert }$ where $\boldsymbol{b}=\boldsymbol{B}/|\boldsymbol{B}|$ is the unit vector in the direction of the local magnetic field and $k_{L}$ a positive characteristic wavenumber associated with the scale of collisionless damping. This representation however does not preserve the zeroth order of the Landau operator and its non-locality in physical space, which are replaced by a normalized first-order derivative along the corresponding magnetic field line. This choice greatly alters the linear behaviour (except at the chosen wavenumber $k_{L}$ ), effectively replacing Landau damping by a usual collisional dissipation. Differently, in Scott (Reference Scott1997), the degree of the Landau operator is preserved by reducing the operator to a multiplicative constant chosen to fit the large-scale dynamics. A Laplacian-type damping term was added in Scott (Reference Scott2002) and in subsequent publications by this author in order to provide a better matching at smaller scales. In this description, the non-locality of the Landau operator is lost. With the aim of preserving the degree of the damping operator, while retaining its space non-locality, a different model, first introduced in Passot et al. (Reference Passot, Henri, Laveder and Sulem2014), was used in Sulem et al. (Reference Sulem, Passot, Laveder and Borgogno2016) and Kobayashi et al. (Reference Kobayashi, Sahraoui, Passot, Laveder, Sulem, Huang, Henri and Smets2017). This is the model we will implement here. Indicating with $\check{f}_{\boldsymbol{k}}$ the $\boldsymbol{k}$ -Fourier coefficient of $f$ , where $\boldsymbol{k}=(k_{x},k_{y},k_{z})$ , the Landau damping operator is phenomenologically modelled in Fourier space as

(2.42)

where $\bar{\unicode[STIX]{x1D70F}}=\langle \boldsymbol{b}\otimes \boldsymbol{b}\rangle$ with $\langle \,\rangle$ the spatial average over the domain ${\mathcal{D}}$ . In particular, this form reduces to the correct formula in the linear case (i.e. when magnetic field lines are straight). We point out that, in particular for the application to the magnetic reconnection problem treated in § 4, the field $A_{\Vert }$ is further decomposed as $A_{\Vert }(x,y,z,t)=A_{\Vert eq}(x)+\unicode[STIX]{x1D716}\tilde{A}_{\Vert }(x,y,z,t)$ , which has to be taken into account when expanding the operator ${\mathcal{L}}$ in terms of the small parameter $\unicode[STIX]{x1D716}$ . Considering the ordering (2.25)–(2.29), one can then see that, according to the closures (2.41), the expansion performed on the operator ${\mathcal{L}}$ has to be of order $O(1)$ , so that $\unicode[STIX]{x1D735}_{\Vert }Q_{\Vert _{e}}\sim \unicode[STIX]{x1D735}_{\Vert }Q_{\bot _{e}}\sim \unicode[STIX]{x1D716}^{2}$ , consistently with the order of the other terms present in (2.32)–(2.33).

As a consequence of the ordering, the evolution of the ion gyrocentre moments $N_{i}$ and $T_{\bot _{i}}$ is essentially electrostatic and neglects variations along the guide field direction. Moreover, choosing $N_{i}=T_{\bot _{i}}=0$ as initial conditions, one has $N_{i}=T_{\bot _{i}}=0$ as solutions. In this limit, the model yields a fluid reduction, with electron temperature anisotropy, of the hybrid kinetic-fluid model of Zocco & Schekochihin (Reference Zocco and Schekochihin2011). Neglecting temperature fluctuations, one retrieves (up to the normalization) the two-field model derived in Schep et al. (Reference Schep, Pegoraro and Kuvshinov1994) and investigated in Borgogno et al. (Reference Borgogno, Grasso, Porcelli, Califano, Pegoraro and Farina2005).

4.1 Cold-ion limit

We consider a simplified version of the model (2.30)–(2.37), where the ion response is determined only by the quasi-neutrality relation in the cold-ion limit. More specifically, we take the solutions $N_{i}=T_{\bot _{i}}=0$ for the ion gyrocentre density and perpendicular temperature fluctuations, and consider the limit $\unicode[STIX]{x1D70F}\ll 1$ . With these assumptions, the model reduces to

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}N_{e}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},N_{e}]+\unicode[STIX]{x1D735}_{\Vert }U_{e}=0, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}(A_{\Vert }-\unicode[STIX]{x1D6FF}^{2}U_{e})-[\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FF}^{2}U_{e}]-\unicode[STIX]{x1D735}_{\Vert }(N_{e}+T_{\Vert _{e}}-\unicode[STIX]{x1D711})=0, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\Vert _{e}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},T_{\Vert _{e}}]+\unicode[STIX]{x1D735}_{\Vert }(Q_{\Vert _{e}}+2U_{e})=0, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T_{\bot _{e}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},T_{\bot _{e}}]+\unicode[STIX]{x1D735}_{\Vert }Q_{\bot _{e}}=0, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}=N_{e}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }={\displaystyle \frac{\unicode[STIX]{x1D6FD}_{e}}{2}}U_{e}. & \displaystyle\end{eqnarray}$$

In such a cold-ion limit, it is also straightforward to rewrite the model in terms of particle moments instead of gyrocentre moments. Actually, due to the assumptions $\unicode[STIX]{x1D6FD}_{e}\sim \unicode[STIX]{x1D6FF}^{2}\ll 1$ , the electron gyrofluid moments are identical to the corresponding particle moments (Brizard Reference Brizard1992), i.e. $N_{e}=n_{e}$ , $U_{e}=u_{e}$ , $T_{\Vert _{e}}=t_{\Vert _{e}}$ , $T_{\bot _{e}}=t_{\bot _{e}}$ , $Q_{\Vert _{e}}=q_{\Vert _{e}}$ and $Q_{\bot _{e}}=q_{\bot _{e}}$ where $n_{e}$ , $u_{e}$ , $t_{\Vert _{e}}$ , $t_{\bot _{e}}$ , $q_{\Vert _{e}}$ and $q_{\bot _{e}}$ indicate the particle density, parallel velocity, parallel and perpendicular temperature and heat flux fluctuations, respectively, for the electron fluid. Ion dynamics enters only through the quasi-neutrality relation, which, in terms of particle moments, yields $n_{i}=\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}=n_{e}$ , where $n_{i}$ denotes the ion density fluctuations (which differ from the ion gyrocentre fluctuations $N_{i}$ ).

In terms of particle moments, the model (4.1)–(4.6) can then be reformulated as

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}]+\frac{2}{\unicode[STIX]{x1D6FD}_{e}}\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }=0, & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}\left(A_{\Vert }-{\displaystyle \frac{2}{\unicode[STIX]{x1D6FD}_{e}}}\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }\right)-{\displaystyle \frac{2}{\unicode[STIX]{x1D6FD}_{e}}}\unicode[STIX]{x1D6FF}^{2}[\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }]-\unicode[STIX]{x1D735}_{\Vert }\left(\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}+t_{\Vert _{e}}-\unicode[STIX]{x1D711}\right)=0, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}t_{\Vert _{e}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},t_{\Vert _{e}}]+\unicode[STIX]{x1D735}_{\Vert }\left(q_{\Vert _{e}}+{\displaystyle \frac{4}{\unicode[STIX]{x1D6FD}_{e}}}\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }\right)=0, & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}t_{\bot _{e}}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},t_{\bot _{e}}]+\unicode[STIX]{x1D735}_{\Vert }q_{\bot _{e}}=0, & \displaystyle\end{eqnarray}$$

with

(4.11a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle q_{\Vert _{e}}=-2\unicode[STIX]{x1D6FC}{\mathcal{L}}t_{\Vert _{e}},\quad q_{\bot _{e}}=-\unicode[STIX]{x1D6FC}{\mathcal{L}}t_{\bot _{e}}. & \displaystyle\end{eqnarray}$$

In the absence of temperature fluctuations (i.e. $t_{\Vert _{e}}=t_{\bot _{e}}=0$ ), the model corresponds, up to the normalization, to the two-field model investigated in Borgogno et al. (Reference Borgogno, Grasso, Porcelli, Califano, Pegoraro and Farina2005) and Grasso, Borgogno & Pegoraro (Reference Grasso, Borgogno and Pegoraro2007).

We remark that, when neglecting dissipation (i.e. $q_{\Vert _{e}}=q_{\bot _{e}}=0$ ), the above model still possesses a Hamiltonian structure. Indeed, the set of functionals of $N_{e}$ , $A_{e}$ , $T_{\Vert _{e}}$ and $T_{\bot _{e}}$ is a sub-algebra with respect to the Poisson bracket (3.4), and in particular with respect to the Poisson bracket $\{\,,\}_{e}$ of (3.5). The model (4.1)–(4.6), for $Q_{\Vert _{e}}=Q_{\bot _{e}}=0$ is then a Hamiltonian system with Poisson bracket given by $\{\,,\}_{e}$ and Hamiltonian functional given by

(4.12) $$\begin{eqnarray}\displaystyle H={\displaystyle \frac{1}{2}}\int d^{3}x\,\left(N_{e}^{2}+\unicode[STIX]{x1D6FF}^{2}U_{e}^{2}+{\displaystyle \frac{2}{\unicode[STIX]{x1D6FD}}}|\unicode[STIX]{x1D735}_{\bot }A_{\Vert }|^{2}+{\displaystyle \frac{T_{\Vert _{e}}^{2}}{2}}+T_{\bot _{e}}^{2}+|\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D711}|^{2}\right). & & \displaystyle\end{eqnarray}$$

The formulation of the Hamiltonian structure in terms of particle moments follows trivially because, as seen above, in this limit, electron gyrofluid moments identify with electron particle moments.

4.2 Linear regime

Equations (4.7)–(4.10) are integrated over the domain $\{-L_{x}\leqslant x\leqslant L_{x},-L_{y}\leqslant y\leqslant L_{y},-L_{z}\leqslant z\leqslant L_{z}\}$ . For the investigation of the linear growth rates, we chose $L_{x}=46,L_{y}=24\unicode[STIX]{x03C0},L_{z}=96\unicode[STIX]{x03C0}$ . The initial condition is given by

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle A_{\Vert }(x,y,z,0)=-\unicode[STIX]{x1D706}\ln \cosh \left({\displaystyle \frac{x}{\unicode[STIX]{x1D706}}}\right)+A_{0}\left(\frac{x-x_{s}}{\unicode[STIX]{x1D706}}\right)\exp [\text{i}(k_{y}y+k_{z}z)]+\text{c.c.}, & \displaystyle\end{eqnarray}$$

(4.14a-c ) $$\begin{eqnarray}\displaystyle t_{\Vert _{e}}(x,y,z,0)=0,\quad t_{\bot _{e}}(x,y,z,0)=0,\quad \unicode[STIX]{x1D711}(x,y,z,0)=0, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}$ is a positive constant whose value, for these numerical simulations, has been set equal to $3$ . The values of the parameters $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D6FD}_{e}$ , on the other hand, have been set as $\unicode[STIX]{x1D6FF}=0.1$ and $\unicode[STIX]{x1D6FD}_{e}=0.05$ , respectively. The initial condition (4.13)–(4.14) corresponds to perturbing a magnetic equilibrium $\boldsymbol{B}_{eq}(x)=\unicode[STIX]{x1D735}A_{\Vert eq}(x)\times \hat{\boldsymbol{z}}+\hat{\boldsymbol{z}}=-\text{tanh}(x/\unicode[STIX]{x1D706})\hat{\mathbf{y}}+\hat{\boldsymbol{z}}$ with a single-helicity disturbance, associated with two wavenumbers $k_{y}=m\unicode[STIX]{x03C0}/L_{y}$ and $k_{z}=n\unicode[STIX]{x03C0}/L_{z}$ , for non-negative integers $m$ and $n$ , and where the helicity is determined by the ratio $k_{z}/k_{y}$ . The amplitude of the perturbation is given by a function $A_{0}$ centred at the resonant surface $x=x_{s}$ , the latter identified by the relation $-\text{tanh}(x_{s}/\unicode[STIX]{x1D706})=k_{z}/k_{y}$ . Because of the choice of the initial conditions, one can infer, from (4.10), that $t_{\bot _{e}}=0$ at all times, so that perpendicular temperature plays no role in these simulations.

For the numerical solution of (4.7)–(4.10) we adopt further developments of the code described in Grasso et al. (Reference Grasso, Borgogno and Pegoraro2007). In particular the code has been parallelized along both the periodic directions $y$ and $z$ and the operator ${\mathcal{L}}$ of (2.41) has been introduced, exploiting its Fourier representation given in (2.42). The new code has been benchmarked in the non-dissipative isothermal limit against previous codes solving the equations described in Borgogno et al. (Reference Borgogno, Grasso, Porcelli, Califano, Pegoraro and Farina2005) and Grasso et al. (Reference Grasso, Borgogno and Pegoraro2007).

In the absence of parallel Landau damping ( $q_{\Vert _{e}}=0$ ), given the initial conditions, the parallel temperature follows an adiabatic behaviour dictated by $t_{\Vert _{e}}=2n_{e}=2\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}$ . In this limit the model can be rewritten as

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}t}}+[\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}]+{\displaystyle \frac{2}{\unicode[STIX]{x1D6FD}_{e}}}\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }=0, & \displaystyle\end{eqnarray}$$

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}\left(A_{\Vert }-{\displaystyle \frac{2}{\unicode[STIX]{x1D6FD}_{e}}}\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }\right)-{\displaystyle \frac{2}{\unicode[STIX]{x1D6FD}_{e}}}\unicode[STIX]{x1D6FF}^{2}[\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }]-\unicode[STIX]{x1D735}_{\Vert }\left(3\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}-\unicode[STIX]{x1D711}\right)=0. & \displaystyle\end{eqnarray}$$

For purely two-dimensional disturbances ( $n=0$ ), an analytical prediction of the linear growth rate of the perturbation for the system (4.15)–(4.16), with initial conditions for $A_{\Vert }$ and $\unicode[STIX]{x1D711}$ as given in (4.13)–(4.14), can be obtained by properly adapting the result of Porcelli (Reference Porcelli1991) valid for isothermal electrons. Note, in fact, that the model equations in the adiabatic limit differ from those in the isothermal limit just by the factor $3$ in the last but one term of (4.16), which becomes a factor $1$ in the isothermal limit. This amounts to replace, in the expression for the growth rate of Porcelli (Reference Porcelli1991), the parameter $\hat{\unicode[STIX]{x1D70C}}_{\unicode[STIX]{x1D70F}}$ by $\sqrt{3}\hat{\unicode[STIX]{x1D70C}}_{\unicode[STIX]{x1D70F}}$ . Assuming that in the linear regime the perturbations for $\unicode[STIX]{x1D711}$ and $A_{\Vert }$ are of the form $\bar{f}(x/\unicode[STIX]{x1D706})\exp [\text{i}(k_{y}y+k_{z}z)-\unicode[STIX]{x1D6FE}t]+\text{c.c.}$ , where $\bar{f}$ is a two-dimensional vector $(\bar{\unicode[STIX]{x1D711}},\bar{A_{\Vert }})$ , the expression for the growth rate $\unicode[STIX]{x1D6FE}$ in the adiabatic limit then becomes

(4.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{\text{adiab}}={\displaystyle \frac{2\sqrt{3}}{\unicode[STIX]{x03C0}}}{\displaystyle \frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FD}_{e}}}k_{y}\unicode[STIX]{x1D6E5}^{\prime }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E5}^{\prime }=(2/\unicode[STIX]{x1D706})((k_{y}\unicode[STIX]{x1D706})^{-1}-k_{y}\unicode[STIX]{x1D706})$ is the stability parameter for collisionless tearing modes for the equilibrium under consideration. Expression (4.17) is valid in the limit of small values of $\unicode[STIX]{x1D6E5}^{\prime }$ . In the absence of parallel temperature fluctuations (isothermal limit), the expression for the growth rate, according to Porcelli (Reference Porcelli1991), reads

(4.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{\text{isoth}}=\unicode[STIX]{x1D6FE}_{\text{adiab}}/\sqrt{3}. & & \displaystyle\end{eqnarray}$$

Therefore, temperature fluctuations tend to increase the growth rate. Although the expressions for the growth rates considered here are valid for the purely two-dimensional case, we expect them to provide sufficiently good approximations, for the sake of the argument that follows, also in the case of oblique modes with $|k_{z}|\ll |k_{y}|$ . This assumption is supported by numerical evidence (Borgogno, Perona & Grasso Reference Borgogno, Perona and Grasso2017).

If the Landau damping term $q_{\Vert _{e}}$ is included, for the magnetic equilibrium under consideration, the action of the operator ${\mathcal{L}}$ in Fourier space corresponds to

(4.19) $$\begin{eqnarray}\displaystyle & \displaystyle \check{{\mathcal{L}}f_{k_{y},k_{z}}}={\displaystyle \frac{1}{\sqrt{(1-(\unicode[STIX]{x1D706}/L_{x})\tanh (L_{x}/\unicode[STIX]{x1D706}))k_{y}^{2}+k_{z}^{2}}}}\left(-{\displaystyle \frac{\unicode[STIX]{x2202}A_{\Vert eq}}{\unicode[STIX]{x2202}x}}\text{i}k_{y}+\text{i}k_{z}\right)\check{f}_{k_{y},k_{z}}, & \displaystyle\end{eqnarray}$$

where $A_{\Vert eq}(x)=-\unicode[STIX]{x1D706}\ln \cosh (x/\unicode[STIX]{x1D706})$ is the equilibrium magnetic flux function. Note that, if compared with the expression for ${\mathcal{L}}$ in Fourier space obtained from (2.42), no Fourier transform has been carried out along the $x$ coordinate due to the choice of a non-periodic equilibrium magnetic field.

The linearized electron parallel temperature (4.9) can then be written, in Fourier space, as

(4.20) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}t_{\Vert e_{k_{y},k_{z}}}}{\unicode[STIX]{x2202}t}} & = & \displaystyle -{\displaystyle \frac{2\unicode[STIX]{x1D6FC}}{\sqrt{(1-(\unicode[STIX]{x1D706}/L_{x})\tanh (L_{x}/\unicode[STIX]{x1D706}))k_{y}^{2}+k_{z}^{2}}}}\left(k_{y}{\displaystyle \frac{\unicode[STIX]{x2202}A_{\Vert eq}}{\unicode[STIX]{x2202}x}}-k_{z}\right)^{2}t_{\Vert e_{k_{y},k_{z}}}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{4}{\unicode[STIX]{x1D6FD}_{e}}}\text{i}\left(k_{y}{\displaystyle \frac{\unicode[STIX]{x2202}A_{\Vert eq}}{\unicode[STIX]{x2202}x}}\left(-{\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}}+k_{y}^{2}\right)+k_{y}{\displaystyle \frac{\unicode[STIX]{x2202}^{3}A_{\Vert eq}}{\unicode[STIX]{x2202}x^{3}}}-k_{z}\left(-{\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}}+k_{y}^{2}\right)\right)A_{\Vert k_{y},k_{z}},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and $t_{\Vert e_{k_{y},k_{z}}},A_{\Vert k_{y},k_{z}}$ are the Fourier coefficients of the perturbations of the magnetic flux function $\tilde{A}$ and of the parallel electron temperature $\tilde{t}_{\Vert _{e}}$ , so that

(4.21) $$\begin{eqnarray}\displaystyle \tilde{A}_{\Vert }(x,y,z,t)=\mathop{\sum }_{ky,kz=-\infty }^{+\infty }A_{\Vert k_{y},k_{z}}(x,t)\exp [\text{i}(k_{y}y+k_{z}z)] & & \displaystyle\end{eqnarray}$$

and

(4.22) $$\begin{eqnarray}\displaystyle \tilde{t}_{\Vert _{e}}(x,y,z,t)=\mathop{\sum }_{ky,kz=-\infty }^{+\infty }t_{\Vert e_{k_{y},k_{z}}}(x,t)\exp [\text{i}(k_{y}y+k_{z}z)]. & & \displaystyle\end{eqnarray}$$

From (4.20) it follows that the first term on the right-hand side, which corresponds to the Landau damping term, vanishes at the resonant surface, which identifies the reconnection region. On the other hand, such term possesses a negative definite coefficient and tends to decrease the amplitude of the squared modulus of the Fourier coefficients $|t_{\Vert e_{k_{y},k_{z}}}|^{2}$ , thus leading to damping of electron temperature fluctuations. On the basis of the above reasoning and of relation (4.18), we expect the Landau damping term to decrease the reconnecting growth rate. This statement is confirmed by numerical simulations and is illustrated in figure 1, where linear growth rates as function of $k_{z}$ , for different values of the poloidal mode number $m$ , are shown for cases including and excluding the Landau damping term. The comparison between simulations with or without the Landau damping term clearly shows that the latter reduces the linear growth rate for all modes considered. By reducing the temperature fluctuations, the Landau damping term inhibits reconnecting modes leading, for $m\geqslant 4$ , to complete stabilization at sufficiently large values of $k_{z}$ , where the corresponding modes in the absence of Landau damping are unstable.

Figure 1.

Linear growth rates of reconnecting modes for simulations in the presence (dashed lines) and in the absence (solid lines) of the Landau damping term. For all modes under consideration, the Landau damping term acts to reduce the growth rate.

We remark that the relations (4.17) and (4.18), obtained from a fluid description, can also be obtained from the formula (B71) of Zocco & Schekochihin (Reference Zocco and Schekochihin2011), which accounts for the full electron kinetic response, considering the specific limits of adiabatic and isothermal electrons, respectively. Such more general formula also pertains to the small $\unicode[STIX]{x1D6E5}^{\prime }$ regime and to modes with $k_{z}=0$ . Denoting with $\unicode[STIX]{x1D6FE}_{\text{kin}}$ the growth rate calculated with the formula (B71) of Zocco & Schekochihin (Reference Zocco and Schekochihin2011), we obtain $\unicode[STIX]{x1D6FE}_{\text{isoth}}<\unicode[STIX]{x1D6FE}_{\text{kin}}<\unicode[STIX]{x1D6FE}_{\text{adiab}}$ , and in particular $\unicode[STIX]{x1D6FE}_{\text{kin}}=0.93\unicode[STIX]{x1D6FE}_{\text{adiab}}$ .

Figure 2.

Comparison between growth rates obtained numerically from the Landau-gyrofluid model (circles) and analytical growth rates (solid line) obtained from the formula (B70) of Zocco & Schekochihin (Reference Zocco and Schekochihin2011). The growth rates refer to perturbations with $k_{z}=0$ .

In figure 2 we show a comparison between numerical growth rates obtained from our Landau-gyrofluid model and analytical growth rates calculated using the formula (B70) of Zocco & Schekochihin (Reference Zocco and Schekochihin2011), valid for arbitrary values of $\unicode[STIX]{x1D6E5}^{\prime }$ . Indeed, the comparison refers to growth rates of modes with $k_{z}=0$ and values of $k_{y}$ ranging in an interval that covers both the small and the large $\unicode[STIX]{x1D6E5}^{\prime }$ regimes (note that $\unicode[STIX]{x1D6E5}^{\prime }\rightarrow +\infty$ as $k_{y}\rightarrow 0^{+}$ ). The growth rates obtained from the Landau-gyrofluid model agree well, over the whole range of $k_{y}$ considered, with the analytical growth rates derived from a formula which takes into account the full kinetic response of the electrons. We remark that, on the basis of the analysis of figure 1, for $k_{z}=0$ , the impact of Landau damping is less pronounced with respect to modes with higher $k_{z}$ , so that for purely two-dimensional disturbances, the adiabatic model could approximate the kinetic results with errors of comparable amplitude, with respect to the Landau-gyrofluid model. The advantages of including Landau damping are then expected to become more evident for modes with larger values of $k_{z}$ .

Figure 3.

Time evolution of the magnetic island area (black) and the $x$ coordinate of the X-point position (cyan) for the cases with (circle) and without (solid line) the Landau damping contribution. Each set of data is plotted versus the normalized time $\unicode[STIX]{x1D6FE}t$ , where $\unicode[STIX]{x1D6FE}$ is the corresponding linear growth rate. Because $\unicode[STIX]{x1D6FE}$ is smaller when Landau damping is taken into account, the growth of the magnetic island with respect to the time $t$ is slower in the case with Landau damping than without.

4.3 Nonlinear regime

In order to investigate the influence of Landau damping on the nonlinear phase of magnetic reconnection, we analysed the plasma dynamics induced by a magnetic perturbation of wavenumbers $k_{y}=1/12$ , $k_{z}=1/48$ in a slab with $L_{x}=46$ , $L_{y}=12\unicode[STIX]{x03C0}$ , $L_{z}=48\unicode[STIX]{x03C0}$ . Many works in the literature showed that the numerical integration of equations (4.7)–(4.10) in the absence of Landau damping effects is non-trivial, due to the co-existence of a broad range of time and length scales (Del Sarto, Califano & Pegoraro Reference Del Sarto, Califano and Pegoraro2003; Grasso et al. Reference Grasso, Borgogno and Pegoraro2007; Chacon & Stanier Reference Chacon and Stanier2016). In the present work, we verified that the numerical solution is even stiffer when the heat fluxes are retained. Within our third order Adams–Bashforth explicit algorithm, the time step in the nonlinear phase is indeed to be decreased, when the $q_{\Vert _{e}}$ field contribution is dominant, by up to two orders of magnitude with respect to the time step used in the linear phase. This makes numerical simulations very expensive in three-dimensional geometry, even when using an efficiently scalable parallel algorithm. In order to reduce the computational time, while preserving the accuracy in the spatial discretization, we adopted an alternative approach valid for single helicity problems. As described in Borgogno et al. (Reference Borgogno, Grasso, Porcelli, Califano, Pegoraro and Farina2005), it consists in the solution of an equivalent two-dimensional problem in the new coordinates $(x,y_{\ast })$ , where $y_{\ast }=y+(k_{z}/k_{y})z$ , with a shifted equilibrium magnetic flux function $A_{\Vert eq}^{\ast }(x)=A_{\Vert eq}(x)-(k_{z}/k_{y})x$ .

Figure 4.

Time variations of the different components of the energy with respect to its initial value, divided by the total energy, for the dissipationless (a) and the dissipative case (b). The curves are plotted versus the normalized time $\unicode[STIX]{x1D6FE}t$ in order to simplify the comparison between the two regimes and refer to the nonlinear phase of the reconnection process. The range of the vertical axis is the same for both (a,b).

Figure 3 shows the evolution of characteristic properties of the magnetic island with and without the Landau term versus the time normalized with respect to the corresponding growth rate. In particular, we plotted in black the logarithm of the magnetic island area, reproducing the reconnected flux at the island X-point, and in cyan the X-point position along the $x$ direction. Note that for the non-symmetric tearing modes the radial position of the magnetic null point varies in time because it does not coincide with a stagnation point, in contrast with the symmetric configurations. Note also that in figure 3 the horizontal axis indicates the time normalized by the growth rate corresponding to each simulation. Since we know from the linear analysis shown in § 4.2 that Landau damping reduces the reconnection growth rate, the observation that the plots for the cases with and without Landau damping essentially overlap, in fact indicates that Landau damping slows down, with respect to the time $t$ , also the nonlinear phase of the reconnection process. On the other hand, figure 3 shows that the Landau term does not affect the saturation level of the magnetic island area. We remark that two-dimensional simulations based on a hybrid model showed agreement, in terms of saturated island width, with the results of a reduced fluid model, as long as the value of the normalized $\unicode[STIX]{x1D6E5}^{\prime }$ parameter did not become too small (Loureiro et al. Reference Loureiro, Schekochihin and Zocco2013). Moreover, we can state that the saturated amplitudes measured in our simulations are comparable to the ones reported in Loureiro et al. (Reference Loureiro, Schekochihin and Zocco2013).

In figure 4 we consider the variation of the different terms contributing to the Hamiltonian (4.12) with respect to their initial value and normalized to the total energy at time $t=0$ , for the two cases under considerations. In particular, the variations are defined as $(1/2)\int d^{3}x(\unicode[STIX]{x1D709}(x,y,z,t)-\unicode[STIX]{x1D709}(x,y,z,0))/H(0)$ , where the function $\unicode[STIX]{x1D709}$ can be replaced by the different contributions to the Hamiltonian $H$ . So one has $E_{pe}$ for $\unicode[STIX]{x1D709}=n_{e}^{2}$ (internal energy due to density fluctuations), $E_{T\Vert }$ for $\unicode[STIX]{x1D709}=t_{\Vert _{e}}^{2}/2$ (parallel thermal energy), $E_{ke}$ for $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D6FF}^{2}u_{e}^{2}$ (parallel kinetic energy), $E_{kp}$ for $\unicode[STIX]{x1D709}=|\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D711}|^{2}$ (perpendicular plasma kinetic energy), $E_{\text{mag}}$ for $\unicode[STIX]{x1D709}=(2/\unicode[STIX]{x1D6FD}_{e})|\unicode[STIX]{x1D735}_{\bot }A_{\Vert }|^{2}$ (magnetic energy). We also defined $E_{\text{landau}}=1-H(t)/H(0)=(1/(2H(0)))\int _{0}^{t}\,\text{d}w\int d^{3}x\,t_{\Vert _{e}}\unicode[STIX]{x1D735}_{\Vert }q_{\Vert _{e}}$ , which corresponds to the energy dissipated via Landau damping. Note that the perpendicular thermal energy $E_{T\bot }$ , corresponding to $\unicode[STIX]{x1D709}=t_{\bot _{e}}^{2}$ , is neglected because of the choice of the initial conditions $\unicode[STIX]{x1D711}(x,y,z,0)=0$ , $t_{\bot _{e}}(x,y,z,0)=0$ , that does not allow any time variation of $t_{\bot _{e}}$ , according to (4.10). We remark that, in the non-dissipative case, the percentage of total energy dissipated numerically when saturation is reached (at approximately $t=10.5/\unicode[STIX]{x1D6FE}$ ) is $|1-H(t=10.5/\unicode[STIX]{x1D6FE})/H(0)|\approx 0.002$ , thus just approximately one fifth of the variation of magnetic energy at the same time, given by $|E_{\text{mag}}(t=10.5/\unicode[STIX]{x1D6FE})|$ . By extrapolating this result to the case with Landau damping, one finds that the energy dissipated numerically is thus a small (although non-negligible) fraction of that dissipated by the Landau damping term.

Before the peak of the reconnection rate, which occurs for $\unicode[STIX]{x1D6FE}t\sim 10$ as shown in figure 5 by the parallel electric field at the island X-point, the energy balance is similar in the two cases. In this time interval, the decrease of the magnetic energy is the same, as well as the growth of the perpendicular plasma kinetic energy that represents the dominant term of the energy balance. The main differences appear in the advanced nonlinear phase, when the process moves towards the saturation. In the Hamiltonian case, due to the formation of small scales in the fields $n_{e}$ , $t_{\Vert _{e}}$ and $u_{e}$ , the contributions of the terms $E_{pe}$ , $E_{ke}$ and $E_{T\Vert }$ become comparable or even higher than $E_{kp}$ . In particular, the internal energy due to electron temperature fluctuations represents the main sink where the magnetic energy is transferred in this phase. On the other hand, when Landau damping is retained, the contribution of all these terms is significantly lower than before and the initial energy is almost exclusively dissipated into electron heating, as observed also in Loureiro et al. (Reference Loureiro, Schekochihin and Zocco2013). We remark that no time lag is observed in our dissipative simulation between the maximum of the reconnection and dissipation rates, according to the fact that this lag is a typical feature of the weakly dissipative systems (Loureiro et al. Reference Loureiro, Schekochihin and Zocco2013).

Figure 5.

Time evolution in the nonlinear regime of the parallel electric field in the dissipative (circle) and dissipationless case (solid line) versus the corresponding normalized time.

The spatial distributions of the fields $\unicode[STIX]{x1D711}$ , $n_{e}$ , $J$ and $T_{\Vert _{e}}$ in the dissipationless and in the dissipative regime we considered are shown in figures 6 and 7, respectively. No significant differences appear in the electrostatic potential, that exhibits in both cases large-scale structures along the separatrices of the magnetic islands. The major discrepancies concern the other fields. The thin layers where $n_{e}$ , $J$ and $T_{\Vert _{e}}$ are located, characteristic of collisionless reconnection (Grasso et al. Reference Grasso, Califano, Pegoraro and Porcelli2001), are smoothed out by Landau damping. The major differences clearly involve the parallel electron temperature, whose peaks are almost twice smaller in the dissipative regime than in the dissipationless case.

Figure 6.

Contour plot of the fields $\unicode[STIX]{x1D711}$ , $n_{e}$ , $u_{e}$ , $T_{\Vert _{e}}$ at the peak of the reconnection rate in the dissipationless regime. The superimposed black lines map the separatrices of the corresponding magnetic island. The range of the vertical axis is the same for both (a,b).

Figure 7.

Contour plot of the fields $\unicode[STIX]{x1D711}$ , $n_{e}$ , $u_{e}$ , $T_{\Vert _{e}}$ at the peak of the reconnection rate in the dissipative regime. The superimposed black lines map the separatrices of the corresponding magnetic island. The range of the vertical axis is the same for both (a,b).