2.1. Gyrofluid

The gyrofluid model used for our analysis is that considered by Granier et al. (Reference Granier, Borgogno, Grasso and Tassi2022b), which consists of the following evolution equations

(2.1)\begin{gather} \frac{\partial N_e}{\partial t}+[G_{10e} \phi - \rho_s^2 2 G_{20e} B_\parallel , N_e]- [G_{10e} A_{{\parallel}} , U_e]=0, \end{gather}

(2.2)\begin{gather}\frac{\partial A_e}{\partial t}+[G_{10e} \phi - \rho_s^2 2 G_{20e} B_\parallel , A_e]+\rho_s^2[G_{10e} A_{{\parallel}} ,N_e]=0, \end{gather}

complemented by the relations

(2.3)\begin{gather} \left( \frac{ G_{10e}^2 -1}{\rho_s^2}+\nabla_{{\perp}}^2\right) \phi-\left( G_{10e} 2 G_{20e} -1\right)B_\parallel=G_{10e} N_e, \end{gather}

(2.4)\begin{gather}\nabla_{\perp}^2 A_{{\parallel}}=G_{10e} U_e, \end{gather}

(2.5)\begin{gather}( G_{10e} 2 G_{20e}-1)\frac{\phi}{\rho_s^2} -\left(\frac{2}{\beta_e}+ 4 G_{20e}^2\right)B_\parallel=2 G_{20e} N_e. \end{gather}

Equation (2.1) corresponds to the electron gyrocentre continuity equation, whereas (2.2) refers to the electron momentum conservation law, along the guide field direction. The static relations (2.3)–(2.5) descend from quasi-neutrality and from the projections of Ampère's law along directions parallel and perpendicular to the guide field, respectively. As mentioned previously, the guide field is directed along the $z$ axis of a Cartesian coordinate system $x,y,z$, and, in the present two-dimensional (2D) version of the model, the dependent variables are functions only of $x$ and $y$, as well as of the time variable $t$. We indicated with $N_e$ and $U_e$ the fluctuations of the electron gyrocentre density and parallel velocity, respectively, whereas $\phi$ and $B_\parallel$ indicate the fluctuations of the electrostatic potential and of the magnetic field along the guide field. The variable $A_e$ is defined by $A_e=G_{10e} A_{\parallel } - d_e^2 U_e$, where $A_{\parallel }$ is the $z$-component of the magnetic vector potential, $d_e=\sqrt {m_e c^2/{4 {\rm \pi}e^2 n_0}}/L$ is the normalised electron skin depth and $G_{10e}$ is an electron gyroaverage operator, defined later in (2.10). The operator $[ \, , \, ]$ is the canonical Poisson bracket and is defined by $[f,g]=\partial _x f \partial _y g - \partial _y f \partial _x g$, for two functions $f$ and $g$. The perpendicular Laplacian operator $\nabla _{\perp }^2$ is defined by $\nabla _{\perp }^2 f=\partial _{xx}f + \partial _{yy} f$. The variables are normalised as

(2.6a–c)\begin{gather} t=\frac{v_A}{L}\hat{t}, \quad x=\frac{\hat{x}}{L}, \quad y=\frac{\hat{y}}{L}, \end{gather}

(2.7a,b)\begin{gather}d_{i} N_{e,i}= \frac{\hat{N}_{e,i}}{n_0}, \quad d_i U_{e,i}=\frac{\hat{U}_{e,i}}{v_A}, \end{gather}

(2.8a–c)\begin{gather}A_{{\parallel}}=\frac{\hat{A}_\parallel}{L B_0}, \quad d_i B_\parallel{=} \frac{\hat{B}_\parallel}{B_0}, \quad \phi=\frac{c}{v_A} \frac{\hat{\phi}}{L B_0}, \end{gather}

where the hat indicates dimensional variables, $c$ is the speed of light, $L$ is a characteristic scale length, $n_0$ is the equilibrium uniform density, $B_0$ is the amplitude of the guide field and $v_A=B_0/\sqrt {4 {\rm \pi}m_i n_0}$ is the Alfvén speed, with $m_i$ indicating the ion mass. The normalised ion skin depth is defined by $d_i=\sqrt {m_i c^2/{4 {\rm \pi}e^2 n_0}}/L$, where $e$ indicates the proton charge. In (2.7a,b) we also introduced the quantities $N_i$ and $U_i$, corresponding to the ion gyrocentre density and parallel velocity fluctuations, respectively. Such moments do not evolve in the model (2.1)–(2.5), and the assumptions on such quantities are discussed later in this section, as well as in § 4. We find it also useful to write explicitly the expression for the magnetic field normalised with respect to the guide field amplitude. In the present 2D setting, by virtue of the normalisation (2.6a–c)–(2.8a–c), such expression is given by

(2.9)\begin{equation} \boldsymbol{B}(x,y,t) = \boldsymbol{z}+ d_i B_\parallel(x,y,t)\boldsymbol{z} + \boldsymbol{\nabla} A_{{\parallel}} (x,y,t)\times \boldsymbol{z}, \end{equation}

where $\boldsymbol {z}$ is the unit vector along $z$. Independent parameters in the model are $\beta _e=8 {\rm \pi}n_0 T_{0e}/B_0^2$, $\rho _s=\sqrt {{T_{0e} m_i c^2}/{(e^2 B_{0}^2)}}/L$ and $d_e$.Footnote ¹ These three parameters correspond to the ratio between equilibrium electron pressure and magnetic guide field pressure, to the normalised sonic Larmor radius and to the electron skin depth, respectively.

The model is formulated on a domain $\{ (x,y) : -L_x \leq x \leq L_x , -L_y \leq y \leq L_y \}$, with $L_x$ and $L_y$ positive constants. Periodic boundary conditions are assumed. This allows to express gyroaverage operators in terms of the corresponding Fourier multipliers. In particular, we associate the electron gyroaverage operators $G_{10e}$ and $G_{20e}$ with corresponding Fourier multipliers in the following way (Brizard Reference Brizard1992)

(2.10)\begin{equation} G_{10e}=2G_{20e} \rightarrow \mathrm{e}^{{-}k_{{\perp}}^2 ({\beta_e}/{4})d_e^2}, \end{equation}

where $k_{\perp }^2=k_x^2 + k_y^2$ is the squared perpendicular wave number and $k_x=m {\rm \pi}/L_x$, $k_y=n {\rm \pi}/L_y$ are the $x$ and $y$ components of the wave vector, with $m$ and $n$ positive integers. As is customary with gyrofluid models, (2.1) and (2.2) are expressed in terms of gyrocentre variables. However, for the sake of the subsequent analysis, it can be useful also to express their relation with particle variables. Such relation, in particular, is affected by the quasi-static assumption, used in the derivation of the model (Tassi, Passot & Sulem Reference Tassi, Passot and Sulem2020) to obtain a closure on the infinite hierarchy of moment equations obtained from a parent gyrokinetic system. As a consequence of such quasi-static closure (which are briefly recalled in § 4) the normalised density fluctuations and parallel velocity fluctuations of the electrons, indicated with $n_e$ and $u_e$, respectively, are related to those of the corresponding gyrocentres by

(2.11)\begin{gather} N_{e}= G_{10e}^{{-}1} \left( n_e + ( G_{10e}^2 -1 ) \frac{\phi}{\rho_s^2} - G_{10e}^2 B_\parallel \right), \end{gather}

(2.12)\begin{gather}U_{e}= G_{10e}^{{-}1} u_e. \end{gather}

In addition, in our gyrofluid model we neglect the contributions due to the density and parallel velocity fluctuations of the ion gyrocentres, by imposing that $N_i=0$, $U_i=0$. Furthermore, ions are assumed to be cold, i.e. $\tau \rightarrow 0$, where $\tau =T_{0i}/T_{0e}$ is the ratio between ion and electron equilibrium temperature.

In terms of the ion particle density and parallel velocity fluctuations, denoted as $n_i$ and $u_i$, respectively, such assumptions lead to the relations $n_i= \nabla _{\perp }^2 \phi + B_\parallel = n_e$ and $u_i=0$. From the quasi-neutrality relation (2.3), Ampère's law (2.4) and (2.5), combined with (2.11) and (2.12), we can obtain the relations

(2.13)\begin{gather} n_e = \frac{2}{2 + \beta_e} \nabla_{{\perp}}^2 \phi={-} \frac{2}{\beta_e} B_\parallel, \end{gather}

(2.14)\begin{gather}u_e = \nabla_{{\perp}}^2 A_{{\parallel}}, \end{gather}

that permit to express the electron particle (as opposed to gyrocentre) density and parallel velocity fluctuations, in terms of electromagnetic perturbations such as $\phi$, $B_\parallel$ and $A_{\parallel }$.

It is also particularly relevant to consider the limit $\beta _e \rightarrow 0$ with $d_e$ and $\rho _s$ remaining finite (which implies $m_e /m_i \rightarrow 0$). This corresponds to suppressing the effects of parallel magnetic perturbations and electron FLR effects. One of the purposes of our investigation is indeed to consider possible modifications, due to kinetic effects, of the plasmoid instability scenario described by Granier et al. (Reference Granier, Borgogno, Comisso, Grasso, Tassi and Numata2022a) and which was conceived namely in the regime with $\beta _e \rightarrow 0$ and finite $d_e$ and $\rho _s$. In this limit, the gyroaverage operators can be approximated in the Fourier space in the following way:

(2.15)\begin{equation} G_{10e} f(x,y) = 2 G_{20e} f(x,y) =f(x,y) + O(\beta_e). \end{equation}

Using this development in (2.1)–(2.5) and neglecting the first-order corrections, we obtain the evolution equations (Schep, Pegoraro & Kuvshinov Reference Schep, Pegoraro and Kuvshinov1994)

(2.16)\begin{gather} \frac{\partial n_e}{\partial t} + [\phi, n_e] - [ A_{{\parallel}}, u_e] = 0, \end{gather}

(2.17)\begin{gather}\frac{\partial}{\partial t} ( A_{{\parallel}} - d_e^2 u_e) + [\phi , A_{{\parallel}} - d_e^2 u_e] - \rho_s^2[n_e, A_{{\parallel}}]=0, \end{gather}

where the static relations (2.3)–(2.5) are replaced by

(2.18)\begin{gather} \nabla_{{\perp}}^2 \phi = N_e =n_e, \end{gather}

(2.19)\begin{gather}\nabla_{\perp}^2 A_{{\parallel}}= U_e= u_e, \end{gather}

(2.20)\begin{gather}B_\parallel=0. \end{gather}

In this limit, particle density and parallel velocity fluctuations coincide with the corresponding gyrocentre counterparts. The system (2.16) and (2.17), complemented by the static relations (2.18) and (2.19) corresponds to the model by Cafaro et al. (Reference Cafaro, Grasso, Pegoraro, Porcelli and Saluzzi1998), adopted for describing collisionless reconnection in the presence of electron inertia and finite sonic Larmor radius effects. Because of the absence of FLR effects, we refer to the model (2.16)–(2.19) as to the fluid limit of the general gyrofluid model (2.1)–(2.5). This model was used extensively for studying reconnection in the past, and a dispersion relation was derived in Porcelli (Reference Porcelli1991).

We point out that it is possible to take the small FLR limit of (2.2), with parameters satisfying ($d_e \ll 1$, $\rho _s \ll 1,$ $d_e/\rho _s \ll 1$, $\beta _e = O(1))$ and $\nabla _{\perp }^2 = O(1)$, to obtain an Ohm's law given by

(2.21)\begin{align} & \frac{\partial}{\partial t}\left( A_{{\parallel}} +\left(\frac{\beta_e}{4}-1\right) d_e^2 \nabla_{{\perp}}^2 A_{{\parallel}} \right)+\left[ \phi , A_{{\parallel}} +\left(\frac{\beta_e}{4}-1\right) d_e^2 \nabla_{{\perp}}^2 A_{{\parallel}} \right] \nonumber\\ & \quad + \rho_s^2\left(\frac{\beta_e}{2 + \beta_e} -1\right)[\nabla_{{\perp}}^2 \phi , A_{{\parallel}}]=0. \end{align}

This equation retains first-order corrections proportional to $(\beta _e/4)d_e^2$ and $\beta _e\rho _s^2/(2+\beta _e)$, that arise from both electron FLR (assuming $G_{10e} \approx 1 +(\beta _e/4)d_e^2 \nabla _\perp ^2$ for $d_e \ll 1$) and finite $B_\parallel$ effects, respectively. The dispersion relation of Porcelli (Reference Porcelli1991), which was obtained by adopting boundary layer and asymptotic matching techniques for $\beta _e=0$, can be extended by identifying an effective electron skin depth $d_e'$ and effective sonic Larmor radius $\rho _s'$, given by $d_e'/d_e = \sqrt {1 - \beta _e/4}$ and $\rho _s'/\rho _s=\sqrt {2/(\beta _e +2})$, respectively. This assumes that one can perform the matching asymptotic expansion as in the case without electron FLR effects, apart from the correction embedded in $d_e '$. This excludes, in particular, the role played, in the dispersion relation, by a possible innermost boundary layer of size proportional to some positive power of $\beta _e$. However, as shown in figure 1, the comparison with numerical simulations suggests that, at least for the small values of $\beta _e$ considered, such approximation, appears to describe effectively the dependence on $\beta _e$ of the linear growth rate.

Figure 1.

Evolution of linear growth rate as a function of $k_y$ for $d_e = 0.085$, $\rho _s = 0.3$. For $\beta _e = 0$, the solution corresponds to that of Porcelli (Reference Porcelli1991), while the cases $\beta _e = 0.12$ and $\beta _e = 0.24$ represent an extension of this solution accounting for small electronic FLR effects and parallel magnetic perturbations.

2.2. Gyrokinetic

In this section, we present the electromagnetic $\delta f$ gyrokinetic model used in this work (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Numata et al. Reference Numata, Howes, Tatsuno, Barnes and Dorland2010), from which the gyrofluid model can be derived with appropriate approximations and closure hypotheses (Tassi et al. Reference Tassi, Passot and Sulem2020). The gyrokinetic model is formulated in terms of the perturbation of the gyrocentre distribution function $g_s = g_s(\textbf {X}_s, v_{\parallel }, v_\perp, t)$ where $v_{\parallel }$, $v_\perp$ are the parallel and perpendicular velocity coordinates. The guiding centre coordinates are given by

(2.22)\begin{equation} \boldsymbol{X}_s =\boldsymbol{x} + \frac{v_{{\rm th}_s}}{\omega_{cs}} \boldsymbol{v}\times \boldsymbol{z}, \end{equation}

where $\boldsymbol {x}$ is the particle position, $\boldsymbol {v}$ is the particle velocity, $v_{{\rm th}_s}=\sqrt {T_{0s}/m_s}$ is the thermal speed and $\omega _{cs}=e B_0 /(m_s c)$ is the cyclotron frequency. The index $s$ labels the particle species, with $s=e$ for electrons and $s=i$ for ions. For simplicity, we assume a uniform background plasma, and two-dimensionality ($\partial /\partial z = 0)$. By adopting the same normalisation scheme with the gyrofluid model, the gyrokinetic system can be written in the following way:

(2.23)\begin{align} & \frac{\partial }{\partial t}\left(g_s + \frac{1}{\rho_s} \frac{ Z_s}{\sqrt{\sigma_s\tau_s}} v_{{\parallel}} \mathcal{J}_{0s} A_{{\parallel}}\right) \nonumber\\ & \quad ={-} \left[ \mathcal{J}_{0s} \phi - \sqrt{\frac{\tau_s \beta_e}{2\sigma_s}} v_{{\parallel}} \mathcal{J}_{0s} A_{{\parallel}} + \frac{\tau_s}{Z_s} \rho_s^2 v_\perp^2 \mathcal{J}_{1s} B_\parallel , g_s + \frac{1}{\rho_s} \frac{Z_s}{\sqrt{\sigma_s \tau_s}} v_{{\parallel}} \mathcal{J}_{0s} A_{{\parallel}} \right], \end{align}

(2.24)\begin{gather} \sum_{s} Z_{s} \bar{n}_s = \phi \sum_{s} \frac{Z_s}{\rho_s^2 \tau_s} \left( \frac{1}{n_0} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} (1 - \mathcal{J}_{0s}) \right) \nonumber\\ +\, B_\parallel \sum_{s} Z_s \left(\frac{1}{n_0} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} v_\perp^2 \mathcal{J}_{0s} \mathcal{J}_{1s} \right), \end{gather}

(2.25)\begin{gather} \sum_{s} Z_{s} \bar{u}_s ={-} \nabla_{{\perp}}^{2} A_{{\parallel}}, \end{gather}

(2.26)\begin{gather} \sum_{s} \bar{p}_s ={-} \phi \sum_{s} \frac{1}{\rho_s^2}\left(\frac{1}{n_0} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} v_\perp^2 \mathcal{J}_{0s} \mathcal{J}_{1s} \right) \nonumber\\ - \,B_\parallel \left( \frac{2}{\beta_e} + \sum_{s} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} ( v_\perp^2 \mathcal{J}_{0s} \mathcal{J}_{1s})^2 \right). \end{gather}

Equation (2.23) is the gyrokinetic equation, whereas (2.24)–(2.26) correspond to the quasi-neutrality relation and to the parallel and perpendicular projection of Ampère's law. We have introduced the following additional normalisations

(2.27a–c)\begin{equation} g_s = \frac{\hat{g}_s}{\hat{\mathcal{F}}_{{\rm eq}_s}}, \quad v_{{\parallel},\perp} = \frac{\hat{v}_{{\parallel},\perp}}{v_{{\rm th}_s}}, \quad d_i \bar{p}_s = \frac{\hat{\bar{p}}_s}{n_0 T_{0e}}, \end{equation}

where the Maxwellian equilibrium distribution function in the dimensional form is

(2.28)\begin{equation} \hat{\mathcal{F}}_{{\rm eq}_s}(\hat{v}_\parallel,\hat{v}_\perp)= n_0 \left(\frac{m_{s}}{{2 {\rm \pi}T_{0s}}}\right)^{3/2} \,\mathrm{e}^{-{m_{s} \hat{v}_\parallel^2}/{2 T_{0s}}-{m_s \hat{v}_\perp^2}/{2 T_{0s}}}. \end{equation}

The species-dependent parameters are the mass ratio $\sigma _s = m_s/m_i$, temperature ratio $\tau _s = T_{0s}/T_{0e}$ and charge number ratio $Z_s = q_s/q_i$. Obviously, $\sigma _i=1$, $\tau _e = 1$, $Z_i=1$. We fix $Z_e=-1$, i.e. $q_i=e$, throughout this work. We may occasionally denote the non-trivial ones as $\sigma _e = \sigma$, $\tau _i = \tau$.

The velocity moments of the distribution function appear in (2.24)–(2.26) are defined by

(2.29)\begin{gather} d_i \bar{n}_{s} = \frac{1}{n_{0}} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} \mathcal{J}_{0s} g_{s} , \end{gather}

(2.30)\begin{gather}d_i \bar{u}_{s} = \sqrt{\frac{\tau_s \beta_e}{2\sigma_s}} \frac{1}{n_{0}} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} v_{{\parallel}} \mathcal{J}_{0s} g_{s}, \end{gather}

(2.31)\begin{gather}d_i \bar{p}_{s} = \tau_{s} \frac{1}{n_{0}} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} v_\perp^2 \mathcal{J}_{1s} g_{s}, \end{gather}

where the volume element ${\rm d} \hat {\mathcal {W}}$ in velocity space is defined as ${\rm d} \hat {\mathcal {W}} = {\rm \pi}v_{{\rm th}_s}^3 \, {\rm d} v_{\parallel } \, {\rm d} v_\perp ^2$. Note that these quantities are moments of $g_s$, thus are different from the actually particle density, flow and pressure, i.e. the moments of the total distribution function $\delta {\mathcal {F}}_s$ defined later on.

Finally, the gyroaverage operators $\mathcal {J}_{0s}$ and $\mathcal {J}_{1s}$ can be expressed, analogously to (2.10), in terms of Fourier multipliers in the following way:

(2.32a,b)\begin{equation} \mathcal{J}_{0s} \rightarrow J_0 (\alpha_s ), \quad \mathcal{J}_{1s} \rightarrow \frac{J_{1}(\alpha_s)}{\alpha_s}, \end{equation}

where $J_0$ and $J_1$ are the zeroth- and first-order Bessel functions, respectively, and the argument is defined by $\alpha _s = k_{\perp } v_\perp v_{{\rm th}_s}/(L \omega _{cs}) = k_{\perp } v_\perp (\rho _s \sqrt {\sigma _s \tau _s}/Z_s)$.

2.3. Connection between the gyrofluid and the gyrokinetic models

We assume the distribution function can be written as

(2.33)\begin{equation} \hat{g}_s(\boldsymbol{X}_s,v_{{\parallel}},v_\perp,t) = \hat{\mathcal{F}}_{{\rm eq}_s} \sum_{n,m=0}^{\infty} H_{m}(v_{{\parallel}}) L_{n}\left(\frac{v_\perp^2}{2}\right) f_{mn_{s}}(\boldsymbol{X}_s,t) \end{equation}

where $H_{m}$ and $L_{n}$ indicate the Hermite and Laguerre polynomials, respectively, of order $m$ and $n$, with $m$ and $n$ non-negative integers. From the orthogonality properties of the Hermite and Laguerre polynomials, the following relation holds:

(2.34)\begin{equation} f_{{mn}_s}=\frac{1}{n_0 \sqrt{m!}}\int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} g_s H_m(v_{{\parallel}}) L_n\left(\frac{v_\perp^2}{2} \right). \end{equation}

The functions $f_{mn_{s}}$ are coefficients of the expansion and are proportional to fluctuations of the gyrofluid moments. Indeed, for instance, $f_{{20}_e}$ is proportional to gyrocentre electron parallel temperature fluctuations and $f_{{00}_e}$ is proportional to gyrocentre electron density fluctuations.

In the present 2D case with an isotropic equilibrium temperature, the system is closed by a closure called ‘quasi-static’ which was derived in Tassi et al. (Reference Tassi, Passot and Sulem2020) and which implies that, with the exception of $N_{e,i}$ and $U_{e,i}$, all the other gyrofluid moments are constrained by the relations

(2.35)\begin{equation} f_{{mn}_s} ={-} \delta_{m0} \left( G_{1ns} \frac{1}{d_i}\frac{2Z_s}{\tau_s\beta_e} \phi + 2 G_{2ns} d_i B_\parallel \right), \end{equation}

where $\delta _{m0}$ is a Kronecker delta and $m$ and $n$ are non-negative integers, with $(m,n) \neq (0,0)$ and $(m,n) \neq (1,0)$, namely to exclude $N_{e,i}$ and $U_{e,i}$. In (2.35), we also introduced the gyroaverage operators which, as Fourier multipliers, are given by

(2.36)\begin{gather} G_{1n} (b_s) \rightarrow \frac{\mathrm{e}^{{-}b_s /2}}{n!}\left(\frac{b_s}{2}\right)^n, \quad n\geq 0, \end{gather}

(2.37a,b)\begin{gather}G_{20} (b_s) \rightarrow \frac{\mathrm{e}^{{-}b_s /2}}{2}, \quad G_{2n} (b_s) \rightarrow -\frac{\mathrm{e}^{{-}b_s/2}}{2} \left( \left(\frac{b_s}{2}\right)^{n-1} \frac{1}{(n-1)!}-\left(\frac{b_s}{2}\right)^{n} \frac{1}{n !}\right), \quad n \geq 1, \end{gather}

with $b_s=k_\perp ^2 v_{{\rm th}_s}^2/(L\omega _{cs})^2 = k_\perp ^2 d_i^2 (\sigma _s \tau _s \beta _e/(2 Z_s^2))=k_{\perp }^2\rho _s^2 (\sigma _s\tau _s/Z_s^2)$. Expressed in terms of particle variables, this closure implies that, with the exception of $n_{e,i}$ and $u_{e,i}$, the fluctuations of the particle moments are zero.

The (2D) quasi-static closure is valid when

(2.38)\begin{equation} \frac{\hat{\omega}}{\hat{k}_y} \ll v_{{\rm th}_s}, \end{equation}

for $s=e,i$ are satisfied, where $\hat {k}_y$ is the $y$ component of the wave vector and $\hat {\omega }$ is the frequency obtained from the dispersion relation of the gyrokinetic equation linearised about an equilibrium ${\phi }^{(0)}={B}_{\parallel }^{(0)}=0,A_{\parallel }^{(0)}=a x$, with a constant $a$ (strictly speaking, in order to perform the gyroaverage of the function $A_{\parallel }^{(0)}=a x$, the linearisation is not carried out on (2.23), which assumes periodic fields, but on its slightly more general form, in which Bessel function operators are replaced by gyroangle averages, so that, for instance, $\mathcal {J}_{0s} f$, for a function $f$, is replaced by $< f (\boldsymbol {x})>_{\boldsymbol {X}_s}=(1/(2{\rm \pi} ))\int _0^{2{\rm \pi} }\, {\rm d} \theta f(\boldsymbol {X}_s - (v_{{\rm th}_s}/\omega _{cs}) \boldsymbol {v}\times \boldsymbol {z})$, with $\theta$ indicating the gyroangle). The condition (2.38) is better fulfilled for waves with small phase velocity along $y$, which justifies the term quasi-static. With regard to the moments not fixed by the quasi-static closure, we have that the dynamics of $N_e$ and $U_e$ is governed by the evolution (2.1) and (2.2), that can then be obtained from the zeroth- and first-order moment, with respect to the parallel velocity coordinate $v_{\parallel }$, of the electron gyrokinetic equations (2.23), with $s=e$.

In the context of the tearing instability, the dispersion relations of the tearing mode usually satisfy the condition $\omega /k_y \ll v_{{\rm th}_e}$. This relation indicates that electrons have time to thermalise along the field lines while the tearing mode develops. Similar comments have already been made in the context of the MHD model with pressure anisotropies. Shi, Lee & Fu (Reference Shi, Lee and Fu1987) discussed the following two equations of state: double adiabatic and isothermal. According to Kulsrud (Reference Kulsrud1983), a double adiabatic closure requires $L/t \gg v_{{\rm th}_e}$, with $t$ the characteristic scale of time variation and $L$ the scale of the spatial variation. Shi et al. (Reference Shi, Lee and Fu1987) indicated that $L/t \gg v_{{\rm th}_e}$ cannot be satisfied by most tearing modes with $t$ and $L$ taken to be the growth time and the wavelength of the mode, respectively. On the other hand, the isothermal closure for the electrons, valid in the opposite regime ($L/t \ll v_{{\rm th}_e}$) can be a better approximation for the study of the tearing instability.

Concerning $N_i$ and $U_i$, as stated previously, we assume the conditions

(2.39a,b)\begin{equation} N_i=0, \quad U_i=0, \end{equation}

to hold. This assumption effectively decouples the ion gyrofluid dynamics from the electron gyrofluid dynamics, leaving (2.1)–(2.5) a closed system.

Whereas all the terms in (2.1)–(2.5) do not assume any ordering on $\beta _e$, the assumption (2.39a,b), as can be derived from the four-field model in Granier et al. (Reference Granier, Borgogno, Grasso and Tassi2022b), is valid for $\beta _e \ll 1$.Footnote ² Thus, in this respect, the adopted gyrofluid model is not derived from the parent gyrofluid model (or from gyrokinetics) under a consistent ordering in $\beta _e$. For this reason, our analysis will be restricted to moderate values of $\beta _e$ (we take $\beta _e \leq 0.692$), where effects of $\beta _e$ associated with electron FLR and parallel magnetic perturbations are nevertheless appreciable. Checking that, in gyrokinetic simulations, the energy associated with the perturbations $N_i$ and $U_i$ remains small, will be an empirical way to make sure that the condition (2.39a,b) is approximately fulfilled. On the other hand, we expect, in the gyrokinetic simulations, a departure from the condition (2.39a,b) as $\beta _e$ increases.

We finally mention that in both gyrofluid and gyrokinetic simulations, we consider the cold-ion case, i.e.

(2.40)\begin{equation} \tau \ll 1, \end{equation}

where we recall that $\tau =\tau _i = T_{0i}/T_{0e}$.

2.4. Numerical set-up

We assume an equilibrium in which the electromagnetic quantities are given by

(2.41a–c)\begin{equation} \phi^{(0)} = 0, \quad A_{{\parallel}}^{(0)} = A_{{\parallel} 0}^{\rm eq} /\cosh^2 ( x ), \quad B_\parallel^{(0)}=0, \end{equation}

with $A_{\parallel 0}^{\rm eq}=1.299$, in order to have $\max _x(B_{y}^{\rm eq}(x))=1$. This equilibrium corresponds to a current sheet centred at $x=0$, of dimensionless length $2L_y=2\hat {L}_y/L$, and of dimensionless width corresponding to unity. Due to periodicity assumption in the simulations,Footnote ³ the dimensionless equilibrium $A_{\parallel }^{(0)}$ is replaced by

(2.42)\begin{equation} A_{{\parallel}}^{(0)}=\sum_{n={-}30}^{n=30} A_{{\parallel} 0}^{\rm eq} a_n \mathrm{e}^{{\rm i}n 2{\rm \pi} x}, \end{equation}

where $a_n$ are the Fourier coefficients of the function $1/\cosh ^2 ( x )$.

Note that, in order to satisfy (2.3)–(2.5) at equilibrium, the corresponding equilibrium density and parallel electron velocity fluctuations of the electron gyrocentres are fixed as

(2.43a,b)\begin{equation} N_e^{(0)}=0, \quad \nabla_{{\perp}}^2 A_{{\parallel}}^{(0)}=G_{10e} U_e^{(0)}. \end{equation}

Also in the gyrokinetic simulations, in accordance with (2.43a,b), the equilibrium current density is assumed to be entirely due to the parallel electron velocity (we recall that, as discussed in § 2.3, the gyrokinetic admits, unlike the gyrofluid model, also a finite parallel ion flow). For both the gyrofluid and gyrokinetic simulations, the perturbation of the equilibrium magnetic flux function $A_{\parallel }^{(0)}$ is of the form $A_{\parallel }^{(1)} \propto \cos (k_y y)$ and is initially excited by the mode $m=1$. The stability condition is given by the tearing parameter (Furth, Killeen & Rosenbluth Reference Furth, Killeen and Rosenbluth1963), which for our equilibrium is given by

(2.44)\begin{equation} \varDelta'= 2 \frac{(5- k_y^2) (k_y^2+3)}{ k_y^2 ( k_y^2+4)^{1/2}}. \end{equation}

The equilibrium (2.41a–c) is tearing unstable when $\varDelta '(k_y) >0$, which corresponds to wave numbers $k_y = {\rm \pi}m/L_y < \sqrt {5}$. In this article, we refer to $\varDelta '$ as that associated with the initially excited mode $m=1$. However, depending on the ratio between the width and the length of the initial current sheet, several wavenumbers $k_y$ with a positive tearing parameter can result from nonlinear interactions of the mode $m=1$. After the saturation of the tearing modes, eventually, the $X$-point collapse and a secondary thinning current sheet forms. The secondary current sheet becomes thinner until reaching a minimum width, and is subject to an inflow that compresses it. Therefore, the instability threshold of this secondary current sheet is not indicated by $\varDelta '$ which is specific to the initial static current sheet. In this secondary current sheet, small perturbations can grow and cause the emergence of other islands, when they enter their nonlinear phase. This dynamics is due to the superposition of several unstable modes with different wavenumbers and growth rates, among which the fastest mode is obviously the dominant one.

The fluid numerical solver SCOPE3D (Solver Collisionless Plasma Equations in 3D) (Granier et al. Reference Granier, Borgogno, Comisso, Grasso, Tassi and Numata2022a) is pseudo-spectral and the advancement in time is done through a third-order Adams–Bashforth scheme. The numerical solver SCOPE3D has been adapted to solve the gyrofluid equations. The gyrokinetic model is solved by AstroGK (Numata et al. Reference Numata, Howes, Tatsuno, Barnes and Dorland2010). Although AstroGK employs some sophisticated techniques for the treatment of linear terms, it uses essentially the same pseudo-spectral and temporal schemes.

4.1. Energy components

As we consider here a plasma with no collisions, the gyrokinetic system solved by AstroGK conserves the total energy (Hamiltonian) (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), normalised by $B_{0}^2/(4{\rm \pi} )$

(4.1)\begin{equation} W(\delta {\mathcal{F}}_{e},\delta {\mathcal{F}}_{i}) = \frac{1}{2} \int \text{d} x \, \text{d} y \left( \frac{\tau_s\beta_e}{2} \sum_s \frac{1}{n_0} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} \delta {\mathcal{F}}_{s}^2 + |\nabla_{{\perp}} A_{{\parallel}}|^2 + d_i^2 |B_\parallel|^2 \right), \end{equation}

where $\delta {\mathcal {F}}_{s}=g_{s} + (1/di) (2Z_{s}/(\tau _s \beta _e)) (\mathcal {J}_{0s}^2-1) \phi + d_i v_\perp ^2 \mathcal {J}_{1s}\mathcal {J}_{0s} B_\parallel$ is the perturbation of the particle distribution function for the species $s$. The first term is the perturbed entropy of the species $s$, while the second term and third terms are the energy of the perpendicular and parallel perturbed magnetic field. We can extract the first two moments from the perturbed particle distribution function as

(4.2)\begin{equation} \delta {\mathcal{F}}_{s} = d_i n_s + \sqrt{\frac{2\sigma_s}{\tau_s\beta_e}} d_i v_{{\parallel}} u_s + h_s', \end{equation}

where the perturbed density and parallel velocity of the particle of species $s$ are denoted as $n_s$ and $u_s$, respectively, and $h_s'$ contains all higher moments of the perturbed distribution function. By definition,

(4.3a,b)\begin{equation} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} h_s' = 0,\quad \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} v_{{\parallel}} h_s' = 0. \end{equation}

We can therefore decompose the expression (4.1) in the following way

(4.4)\begin{align} W(\delta {\mathcal{F}}_{e}, \delta {\mathcal{F}}_{i}) & = \frac{1}{2} \int \text{d} x \, \text{d} y \left( \sum_s \left( \tau_s \rho_s^2 n_s^2 + \sigma_s d_i^2 u_s^2 + \frac{\tau_s\beta_e}{2} \frac{1}{n_0} \int \text{d} \hat{\mathcal{W}} \hat{\mathcal{F}}_{{\rm eq}_s} {h_s'}^2 \right) \right. \nonumber\\ & \quad +\left. \vphantom{\sum_s}|\nabla_{{\perp}} A_{{\parallel}}|^2 + d_i^2 |B_\parallel^2| \right). \end{align}

The first term is the energy generated by the electron density variance, the second term is the kinetic energy of the parallel electron flow, and the third term is the free electron energy.

With regard to the collisionless gyrofluid model, the system of (2.1)–(2.5) possesses a conserved Hamiltonian given by

(4.5)\begin{equation} H_{gf}(N_e , A_e )=\frac{1}{2} \int \text{d} x \, \text{d} y ( \rho_s^2 N_e^2 + d_e^2 U_e^2 + | \nabla_\perp A_{{\parallel}} |^2 - N_e ( G_{10e} \phi - \rho_s^2 2 G_{20e} B_\parallel)). \end{equation}

We remark that, as is shown in the Appendix, the form of the Hamiltonian (4.5), obtained from the quasi-static closure, is the same that one obtains by imposing what we refer to as an isothermal gyrofluid closure (the relations between $\phi, A_{\parallel }, B_\parallel$ and $N_e, U_e$ will, however, be different in the two cases).

Using the relation (2.13) and (2.14) we can also write the Hamiltonian in terms of particle variables as follows:

(4.6)\begin{align} H_{p}(n_e ,A_e )& = \frac{1}{2} \int \text{d} x \, \text{d} y \left( \rho_s^2 n_e G_{10e}^{{-}2} n_e + d_e^2 ( G_{10e}^{{-}1}u_e )^2 + | \nabla_\perp A_{{\parallel}} |^2 + d_i^2 |B_\parallel|^2 \vphantom{\frac{\phi}{\rho_s^2}}\right.\nonumber\\ & \quad + \left. n_e ( 1 - 2G_{10e}^{{-}2} )\phi + \phi ( G_{10e}^{{-}2}-1 ) \frac{\phi}{\rho_s^2} \right). \end{align}

When we consider the limit $\beta _e$, $m_e/m_i \rightarrow 0$ the Hamiltonian of the gyrofluid equations is reduced to

(4.7)\begin{equation} H_p(n_e , A_e)=\frac{1}{2} \int \text{d} x \, \text{d} y ( \rho_s^2 n_e^2 + d_e^2 u_e^2 + | \nabla_\perp A_{{\parallel}}|^2 + |\nabla_{{\perp}} \phi|^2), \end{equation}

which is namely the Hamiltonian of the fluid (2.16) and (2.17). In (4.7), the contribution from left to right are the energy generated by the electron density fluctuation, the parallel electron kinetic energy, the perpendicular magnetic energy and the perpendicular plasma kinetic energy which is essentially the $\boldsymbol {E} \times \boldsymbol {B}$ flow energy.

4.2. Negligible $\beta _e$: fluid versus gyrokinetic

On figure 11 we present the comparison between the energy variation of the fluid case $1_F$ and that of the low $\beta _e$ gyrokinetic case $1_{\rm GK3}$ ($\beta _e=0.062$). The variations are defined as $(1/2)\int \text {d} \boldsymbol {r}( \xi (x,y,t) - \xi (x,y,0)) / E(0)$ where the function $\xi$ can be replaced by the different contributions to $\hat {W}$ and $H$ (where $\hat {W}$ is also considered in the 2D limit) and $E(0)$ is the initial total energy. On the gyrokinetic plots, the four main energy channels are shown as solid lines. The solid purple line is the total ion energy variation. We also show the evolution of the variations relative to the density variance (dashed dotted), the parallel kinetic energy (densely dashed) and the perpendicular kinetic energy (loosely dashed), that are components of the total particle energy. The same channels are shown for the electrons in green.

Figure 11.

Time evolution of the energy variations for the cases $1_F$ (a) and $1_{\rm GK3}$ (b). (c) and (d) show the time evolution of the different contributions to the parallel kinetic energy for the corresponding simulations. No plasmoids in this case.

The amount of magnetic energy that is converted is identical between fluid and gyrokinetics and appears to be transferred mainly to the electrons. On the other hand, it is not identically distributed in the gyrokinetic and fluid frameworks. For the fluid simulations, the magnetic energy has no choice but to be converted into electron density fluctuations or electron parallel acceleration, whereas in the gyrokinetic case, there is little energy sent to these channels. This suggests that, in the gyrokinetic framework, the energy of the electrons increases due the fluctuations of the higher-order moments of the distribution function due to phase mixing (Loureiro, Schekochihin & Zocco Reference Loureiro, Schekochihin and Zocco2013; Numata & Loureiro Reference Numata and Loureiro2015), such as, for instance, the perpendicular and parallel electron temperature. It is likely that the magnetic energy is actually converted into thermal electron energy. Such possibility is prevented in the fluid case because, as a consequence of the closure, for $\beta _e \rightarrow 0$, no temperature fluctuations are allowed.

The striking difference between the two approaches is that the parallel electron kinetic energy increases in the fluid case, whereas it is quasi-constant or decreasing in the gyrokinetic one (figure 11). In order to investigate the origin of this difference, we performed an initial condition check and decomposed the parallel electron kinetic energy. The decomposition leads to three energy components, namely the equilibrium part ($u_{\rm eq}^2$), the perturbation part ($\tilde {u}_e^2$) and the cross term ($2\tilde {u}_eu_{\rm eq}$). The change of each component is shown on the bottom panel of figure 11. The equilibrium contribution clearly does not change in time. The quadratic perturbation part is always positive but globally the variation of parallel electron kinetic energy can decrease because of the cross term becoming negative, which is the case for the gyrokinetic simulation. For the fluid case, the perturbation term increases considerably, leading to a positive variation of the parallel kinetic energy, since the electrons are highly accelerated for conservation of the total energy.

With regard to the ions, the closure assumptions imply an even rougher approximation of the ion dynamics, in the fluid case, with respect to gyrokinetics. In the gyrokinetic case, for low $\beta _e$, we can see on figure 11 that the main component of the total ion energy consists of the perpendicular kinetic energy, which is included in the ${h_s}'$ part in (4.4),Footnote ⁴ given by

(4.8)\begin{equation} \frac{1}{2} \int \text{d} x \, \text{d} y d_i^2 \boldsymbol{u}_{{\perp},i} ^2, \end{equation}

where the perpendicular ion velocity $\boldsymbol {u}_{\perp, i}$ is calculated directly from its definition as a moment in the following way:

(4.9)\begin{equation} d_i \boldsymbol{u}_{{\perp},i} = \sqrt{\frac{\tau_i\beta_e}{2\sigma_i}}\frac{1}{n_{0}} \int \text{d} \hat{\mathcal{W}} {\mathcal{F}}_{{\rm eq}_i} \boldsymbol{v}_{{\perp}} \delta {\mathcal{F}}_{i}. \end{equation}

Note that the perpendicular flow holds the identity from the definition

(4.10)\begin{equation} d_i \boldsymbol{u}_{{\perp},i} = (-\boldsymbol{\nabla} \phi - \rho_s^2 \boldsymbol{\nabla} \boldsymbol{\cdot} \texttt{p}_{{\perp}{\perp},i}) \times \boldsymbol{z}, \end{equation}

where the perpendicular pressure tensor is given by

(4.11)\begin{equation} d_i \texttt{p}_{{\perp}{\perp},i} = \tau_i \frac{1}{n_0} \int \text{d} \hat{\mathcal{W}} {\mathcal{F}}_{{\rm eq}_i} \boldsymbol{v}_{{\perp}} \boldsymbol{v}_{{\perp}} \delta {\mathcal{F}}_{i}. \end{equation}

The perpendicular flow is given by the sum of $\boldsymbol {E}\times \boldsymbol {B}$ drift and diamagnetic drift of perturbed pressure.

In spite of the closure, the evolution of the energy component (4.8) is very similar to that of the $\boldsymbol {E}\times \boldsymbol {B}$ flow energy of the gyrofluid case. For a very small $\beta _e$, no parallel ion kinetic energy and parallel magnetic energy seems to be generated.

4.3. Finite $\beta _e$: gyrofluid versus gyrokinetic

When $\beta _e$ is very small, the FLR corrections become negligible and the particle and gyrocentre variables coincide. On the other hand, for non-negligible $\beta _e$, the electron Larmor radius becomes finite and the relations (2.11) and (2.12) allow us to relate the density and parallel velocity of the particles to those of the gyrocentres. In figure 12 we compare the gyrofluid energy variations with the gyrokinetic ones for $0 < \beta _e < 1$. For this purpose, we use the simulation set for $p=3$.

Figure 12.

Time evolution of the energy variations for the cases $3_{GF}$ and $3_{GK}$.

In the plot referring to the gyrofluid energy, we show the variation of both the particles and gyrocentres energy. For instance, the curve referring to ‘Kin$_{\parallel e}$’ corresponds to the variation of $(1/2)\int {\rm d} x \, {\rm d} y d_e^2 u_e^2$, which is comparable to the second term of the gyrokinetic energy (4.4). The one referring to ‘Gyrocentre Kin$_{\parallel e}$’ corresponds to the variation of $(1/2)\int {\rm d} x\, {\rm d} y d_e^2 U_e^2$. By increasing $\beta _e$, the difference between the variation of the energy of the gyrocentres and that of the particles increases. With finite $\beta _e$, we now note a loss of parallel kinetic energy of the electrons for the gyrofluid case, which is in better agreement with the gyrokinetic approach. Increasing $\beta _e$, will also generate more parallel magnetic energy, which is well reproduced by the gyrofluid model. On the other hand, the gyrokinetic cases indicate that a significant part of the magnetic energy is now converted into parallel ion kinetic energy. As mentioned previously, a limitation of the reduced gyrofluid model is that the ion parallel velocity has been ‘artificially’ removed by imposing $U_i=u_i=0$. The limitations of this assumption become evident, in particular, from figure 12 which shows that, in the gyrokinetic case, for sufficiently large $\beta _e$, the ion fluid is actually accelerated along the $\boldsymbol {z}$ axis. On the other hand, it seems that despite this missing element, the gyrofluid model is suitable for studying the formation of plasmoid for $0<\beta _e<1$.