# Noncollisional plasmoid instability based on a gyrofluid and gyrokinetic
integrated approach

C. Granier 1,a, R. Numata2 D. Borgogno3 E. Tassi1 and D. Grasso3

<EMAIL_ADDRESS>
1Université Côte d’Azur, CNRS, Observatoire de la Côte d’Azur, Laboratoire J.
L. Lagrange, Boulevard de l’Observatoire, CS 34229, 06304 Nice Cedex 4, France
2Graduate School of Information Science, University of Hyogo, Kobe 650-0047,
Japan
3Istituto dei Sistemi Complessi - CNR and Dipartimento di Energia, Politecnico
di Torino, Torino 10129, Italy

Abstract

In this work, the development of two-dimensional current sheets with respect
to tearing-modes, in collisionless plasmas with a strong guide field, is
analysed. During their non-linear evolution, these thin current sheets can
become unstable to the formation of plasmoids, which allows the magnetic
reconnection process to reach high reconnection rates. We carry out a detailed
study of the impact of a finite $\beta_{e}$, which also implies finite
electron Larmor radius effects, on the collisionless plasmoid instability.
This study is conducted through a comparison of gyrofluid and gyrokinetic
simulations. The comparison shows in general a good capability of the
gyrofluid models in predicting the plasmoid instability observed with
gyrokinetic simulations. We show that the effects of $\beta_{e}$ promotes the
plasmoid growth. The impact of the closure applied during the derivation of
the gyrofluid model is also studied through the comparison of the energy
variation.

##  1 Introduction

Magnetic reconnection is a change of topology of the magnetic field lines
taking place in regions of intense localized current, referred to as current
sheets. This fundamental process ultimately converts magnetic energy into bulk
flow and particle heating, and is responsible for the explosive release of
magnetic energy in astrophysical and laboratory plasmas. The instabilities of
very elongated reconnecting current sheets leading to the formation of
secondary magnetic islands, called plasmoids, have generated a lot of
interest, as they are believed to achieve fast reconnection. Plasmoids have
been greatly studied through the most standard reconnection model based on the
Sweet-Parker (SP) theory in the resistive magnetohydrodynamics (RMHD)
framework (Biskamp (1986); Loureiro et al. (2005)). In Biskamp (1986), it has
been shown that collisional current sheets become unstable when
$S=\mu_{0}L_{cs}v_{A}/\eta>10^{4}$, where $L_{cs}$ is the length of the
current sheet, $\eta$ is the resistivity and $v_{A}$ is the Alfvén speed. Much
work has followed and allowed to identify the plasmoid regime as a function of
the Lundquist number $S$ and of the characteristic scale of a dynamic of the
ions (ion-sound Larmor radius $\rho_{s}$ or ion skin depth $d_{i}$ scales) at
which a transition to a non-collisional regime, dominated by kinetic effects,
occurs (Ji and Daughton (2011); Daughton and Roytershteyn (2012); Loureiro and
Uzdensky (2015); Uzdensky et al. (2010a, b); Bhat and Loureiro (2018)). This
extension of the resistive reconnection regime with the inclusion of the ion
dynamics enlarged the study to a broader parameter space, but also suggested
that plasmoids are fundamental features of reconnecting current sheets,
regardless the value of the Lundquist number (Ji and Daughton (2011); Daughton
and Roytershteyn (2012)).

The plasma in the magnetosphere and solar wind, which regularly undergoes
reconnection, is so dilute that collisions between particles are extremely
infrequent. In such plasmas, electron inertia becomes particularly relevant to
drive reconnection in thin current sheets. Indeed, recent observations
revealed many reconnection onsets driven by electrons, in the presence of a
strong guide field, close to the dayside magnetopause and magnetosheath Burch
et al. (2016); Phan et al. (2018) with current sheets having a thickness of
the order of the electron inertial length. Regarding experiments, a study by
Olson et al. (2016a) also gave direct experimental proof of plasmoid formation
at the electron scale in a weakly collisional regime. In these collisionless,
magnetized environments, effects of the finite electron Larmor radius (FLR) on
the reconnection process can also become non-negligible, in particular when
$\beta_{e}$, being defined as the ratio between equilibrium thermal electron
pressure and guide field magnetic pressure, is not much smaller than unity.
This motivates the study of the formation of plasmoids in non-collisional
current sheets, and in particular, the impact of the effects relevant at the
electron scales such as the electron skin depth and the electron Larmor
radius.

In Granier et al. (2022a) instability thresholds for the purely collisionless
plasmoid onset were clearly identified in the regime $\beta_{e}\rightarrow 0$.
In this article we relax the assumption of small $\beta_{e}$ and carry out a
detailed study of the impact of a finite $\beta_{e}$, on the collisionless
plasmoid instability, in the case of a strong guide field. We consider
inertial reconnection, and finite electron FLR effects arise from the
combination of electron inertia and finite $\beta_{e}$ parameters. This study
is conducted through a comparison of gyrofluid and gyrokinetic simulations.
Both approaches are assuming that the plasma is immersed in a strong guide
field directed along the $z$ direction. As a by-product of our analysis, we
also obtain a way to validate, by means of gyrokinetic simulations, part of
the results on collisionless plasmoid instability obtained by Granier et al.
(2022a) with a gyrofluid approach in the $\beta\rightarrow 0$ limit (later
referred to as fluid limit)

The gyrofluid model is the 2-field system presented in (Granier et al.
(2022b)) and assumes cold and immobile ions along the guide field direction.
Gyrofluid models, although greatly simplified with respect to the original
gyrokinetic system, are useful tools for studying collisionless reconnection,
in which the microscopic scales, such as the electron skin depth and the
electron Larmor radius, can be more important than resistivity. In addition,
the gyrofluid framework is less costly in terms of computational resources,
and physically more intuitive when compared to the kinetic or gyrokinetic
framework. So far, gyrofluid modelling allowed to gain a good understanding of
the role of collisionless effects (e.g. Sarto et al. (2011); Comisso et al.
(2013); Tassi et al. (2018); Granier et al. (2021, 2022b)).

The gyrokinetic model, adopted for the comparison, is a $\delta f$ model which
solves the electromagnetic gyrokinetic Vlasov-Maxwell system. The gyrokinetic
equations are solved by means of the AstroGK code, presented and used in
Numata et al. (2010); Numata and Loureiro (2015). One of the main advantages
of using the AstroGK code for a comparison with the gyrofluid results, is
that, in a specific limit, the gyrokinetic system solved by AstroGK reduces to
the one that was taken to derive the 2-field gyrofluid model used in this
study (Howes et al., 2006). This allows to study the impact of the closure
applied on the moments, performed during the derivation of the gyrofluid
model, on the distribution and conversion of energy during reconnection and
identify the possible limitations of the gyrofluid approach. The specific
limit in which the AstroGK code has to be used, in order to reproduce the
parent gyrokinetic model of the gyrofluid model, is that corresponding to a
straight guide field, with no density and temperature gradients and without
collisions. To be consistent with the gyrofluid approach, the ions are assumed
to be cold.

The article is organized as follow. In Sec. 2 we present the gyrofluid and
gyrokinetic systems, as well as the numerical set up. In Sec. 3 we present the
results concerning the plasmoid instability obtained from a comparison of the
two approaches. In Sec. 4 we compare the energy variations in the two
frameworks and discuss the impact of the closure hypothesis on the energy
conversion. Section 5 is devoted to conclusions.

##  2 Adopted models

###  2.1 Gyrofluid

The gyrofluid model used for our analysis is the one considered by Granier et
al. (2022b), which consists of the following evolution equations

$\displaystyle\frac{\partial N_{e}}{\partial
t}+[G_{10e}\phi-\rho_{s}^{2}2G_{20e}B_{\parallel},N_{e}]-[G_{10e}A_{\parallel},U_{e}]=0,$
(1) $\displaystyle\frac{\partial A_{e}}{\partial
t}+[G_{10e}\phi-\rho_{s}^{2}2G_{20e}B_{\parallel},A_{e}]+\rho_{s}^{2}[G_{10e}A_{\parallel},N_{e}]=0,$
(2)

complemented by the relations

$\displaystyle\left(\frac{G_{10e}^{2}-1}{\rho_{s}^{2}}+\nabla_{\perp}^{2}\right)\phi-\left(G_{10e}2G_{20e}-1\right)B_{\parallel}=G_{10e}N_{e},$
(3) $\displaystyle\nabla_{\perp}^{2}A_{\parallel}=G_{10e}U_{e},$ (4)
$\displaystyle\left(G_{10e}2G_{20e}-1\right)\frac{\phi}{\rho_{s}^{2}}-\left(\frac{2}{\beta_{e}}+4G_{20e}^{2}\right)B_{\parallel}=2G_{20e}N_{e}.$
(5)

Equation (1) corresponds to the electron gyrocenter continuity equation,
whereas Eq. (2) refers to the electron momentum conservation law, along the
guide field direction. The static relations (3), (4) and (5) descend from
quasi-neutrality and from the projections of Ampère’s law along directions
parallel and perpendicular to the guide field, respectively. The guide field
is directed along the $z$ axis of a Cartesian coordinate system $x,y,z$, and,
in the present 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 gyrocenter 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\pi
e^{2}n_{0}}}/L$ is the normalized electron skin depth and $G_{10e}$ is an
electron gyroaverage operator, defined later in Eq. (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
normalized as

$\displaystyle t=\frac{v_{A}}{L}\hat{t},\qquad x=\frac{\hat{x}}{L},\qquad
y=\frac{\hat{y}}{L},$ (6) $\displaystyle
d_{i}N_{e,i}=\frac{\hat{N}_{e,i}}{n_{0}},\qquad
d_{i}U_{e,i}=\frac{\hat{U}_{e,i}}{v_{A}},$ (7) $\displaystyle
A_{\parallel}=\frac{\hat{A}_{\parallel}}{LB_{0}},\qquad
d_{i}B_{\parallel}=\frac{\hat{B}_{\parallel}}{B_{0}},\qquad\phi=\frac{c}{v_{A}}\frac{\hat{\phi}}{LB_{0}},$
(8)

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\pi
m_{i}n_{0}}$ is the Alfvén speed, with $m_{i}$ indicating the ion mass. The
normalized ion skin depth is defined by $d_{i}=\sqrt{m_{i}c^{2}/{4\pi
e^{2}n_{0}}}/L$, where $e$ indicates the proton charge. In Eq. (7) we also
introduced the quantities $N_{i}$ and $U_{i}$, corresponding to the ion
gyrocenter density and parallel velocity fluctuations, respectively. Such
moments do not evolve in the model (1)-(5), and the assumptions on such
quantities will be discussed later in this section, as well as in Sec. 4. We
find it also useful to write explicitly the expression for the magnetic field
normalized with respect to the guide field amplitude. In the present 2D
setting, by virtue of the normalization (6)-(8), such expression is given by

$\mathbf{B}(x,y,t)=\mathbf{z}+d_{i}B_{\parallel}(x,y,t)\mathbf{z}+\nabla
A_{\parallel}(x,y,t)\times\mathbf{z},$ (9)

where $\mathbf{z}$ is the unit vector along $z$. Independent parameters in the
model are $\beta_{e}=8\pi n_{0}T_{0e}/B_{0}^{2}$,
$\rho_{s}=\sqrt{{T_{0e}m_{i}c^{2}}/{(e^{2}B_{0}^{2})}}/L$111According to a
customary notation, in the symbols $\rho_{s}$, the subscript $s$ is to
indicate a sonic quantity and not the particle species. and $d_{e}=$, where
$T_{0e}$ is the uniform equilibrium electron temperature and $m_{e}$ is the
electron mass. These three parameters correspond to the ratio between
equilibrium electron pressure and magnetic guide field pressure, to the
normalized 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, 1992)

$G_{10e}=2G_{20e}\rightarrow\mathrm{e}^{-k_{\perp}^{2}\frac{\beta_{e}}{4}d_{e}^{2}},$
(10)

where $k_{\perp}^{2}=k_{x}^{2}+k_{y}^{2}$ is the squared perpendicular wave
number and $k_{x}=m\pi/L_{x}$, $k_{y}=n\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, Eqs. (1) and (2) are expressed in terms of
gyrocenter 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 et al., 2020) 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 will be briefly
recalled in Sec. 4) the normalized density fluctuations and parallel velocity
fluctuations of the electrons, indicated with $n_{e}$ and $u_{e}$,
respectively, are related to those of the corresponding gyrocenters by

$\displaystyle
N_{e}=G_{10e}^{-1}\left(n_{e}+\left(G_{10e}^{2}-1\right)\frac{\phi}{\rho_{s}^{2}}-G_{10e}^{2}B_{\parallel}\right),$
(11) $\displaystyle U_{e}=G_{10e}^{-1}u_{e}.$ (12)

Also, in our gyrofluid model we neglect the contributions due to the density
and parallel velocity fluctuations of the ion gyrocenters, 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 (3), Ampère’s law (4)-(5), combined with
Eqs. (11)-(12), we can obtain the relations

$\displaystyle
n_{e}=\frac{2}{2+\beta_{e}}\nabla_{\perp}^{2}\phi=-\frac{2}{\beta_{e}}B_{\parallel},$
(13) $\displaystyle u_{e}=\nabla_{\perp}^{2}A_{\parallel},$ (14)

that permit to express the electron particle (as opposed to gyrocenter)
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 Ref.
Granier et al. (2022a) 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:

$G_{10e}f(x,y)=2G_{20e}f(x,y)=f(x,y)+O(\beta_{e}).$ (15)

Using this development in Eqs. (1) - (5) and neglecting the first order
corrections, we obtain the evolution equations (Schep et al. (1994))

$\frac{\partial n_{e}}{\partial t}+[\phi,n_{e}]-[A_{\parallel},u_{e}]=0,$ (16)
$\begin{split}\frac{\partial}{\partial
t}\left(A_{\parallel}-d_{e}^{2}u_{e}\right)+\left[\phi,A_{\parallel}-d_{e}^{2}u_{e}\right]-\rho_{s}^{2}[n_{e},A_{\parallel}]=0,\end{split}$
(17)

where the static relations (3) - (5) are replaced by

$\displaystyle\nabla_{\perp}^{2}\phi=N_{e}=n_{e},$ (18)
$\displaystyle\nabla_{\perp}^{2}A_{\parallel}=U_{e}=u_{e},$ (19)
$\displaystyle B_{\parallel}=0.$ (20)

In this limit, particle density and parallel velocity fluctuations coincide
with the corresponding gyrocenter counterparts. The system (16)-(17),
complemented by the static relations (18)-(19) corresponds to the model by
Cafaro et al. (1998), 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 will refer to the model (16)-(19) as to the
fluid limit of the general gyrofluid model (1)-(5).

###  2.2 Gyrokinetic

In this section, we present the electromagnetic $\delta f$ gyrokinetic model
used in this work (Howes et al., 2006; Numata et al., 2010), from which the
gyrofluid model can be derived with appropriate approximations and closure
hypotheses Tassi et al. (2020). The gyrokinetic model is formulated in terms
of the perturbation of the gyrocenter 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 center coordinates is given by

$\textbf{X}_{s}=\mathbf{x}+\frac{v_{{th}_{s}}}{\omega_{cs}}\textbf{v}\times\mathbb{z},$
(21)

where $\mathbf{x}$ is the particle position, v is the particle velocity,
$v_{{th}_{s}}=\sqrt{T_{0s}/m_{s}}$ is the thermal speed and
$\omega_{cs}=eB_{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 normalization scheme with the
gyrofluid model, the gyrokinetic system can be written in the following way,

$\displaystyle\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)+\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]=0,$
(22)
$\displaystyle\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}}_{{eq}_{s}}(1-\mathcal{J}_{0s})\right)+B_{\parallel}\sum_{s}Z_{s}\left(\frac{1}{n_{0}}\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}v_{\perp}^{2}\mathcal{J}_{0s}\mathcal{J}_{1s}\right),$
(23) $\displaystyle\sum_{s}Z_{s}\bar{u}_{s}=-\nabla_{\perp}^{2}A_{\parallel},$
(24)
$\displaystyle\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}}_{{eq}_{s}}v_{\perp}^{2}\mathcal{J}_{0s}\mathcal{J}_{1s}\right)-B_{\parallel}\left(\frac{2}{\beta_{e}}+\sum_{s}\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}(v_{\perp}^{2}\mathcal{J}_{0s}\mathcal{J}_{1s})^{2}\right)$
(25)

Eq. (22) is the gyrokinetic equation, whereas Eqs. (23), (24) and (25)
correspond to the quasi-neutrality relation and to the parallel and
perpendicular projection of Ampère’s law. We have introduced the following
additional normalizations

$\displaystyle g_{s}=\frac{\hat{g}_{s}}{\hat{\mathcal{F}}_{{eq}_{s}}},\qquad
v_{\parallel,\perp}=\frac{\hat{v}_{\parallel,\perp}}{v_{{th}_{s}}},\qquad
d_{i}\bar{p}_{s}=\frac{\hat{\bar{p}}_{s}}{n_{0}T_{0e}}$ (26)

where the Maxwellian equilibrium distribution function in the dimensional form
is

$\hat{\mathcal{F}}_{{eq}_{s}}(\hat{v}_{\parallel},\hat{v}_{\perp})=n_{0}\left(\frac{m_{s}}{{2\pi
T_{0s}}}\right)^{3/2}\mathrm{e}^{-\frac{m_{s}\hat{v}_{\parallel}^{2}}{2T_{0s}}-\frac{m_{s}\hat{v}_{\perp}^{2}}{2T_{0s}}}.$
(27)

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 Eqns. (23)-(25)
are defined by

$\displaystyle d_{i}\bar{n}_{s}$
$\displaystyle=\frac{1}{n_{0}}\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}\mathcal{J}_{0s}g_{s},$
(28) $\displaystyle d_{i}\bar{u}_{s}$
$\displaystyle=\sqrt{\frac{\tau_{s}\beta_{e}}{2\sigma_{s}}}\frac{1}{n_{0}}\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}v_{\parallel}\mathcal{J}_{0s}g_{s},$
(29) $\displaystyle d_{i}\bar{p}_{s}$
$\displaystyle=\tau_{s}\frac{1}{n_{0}}\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}v_{\perp}^{2}\mathcal{J}_{1s}g_{s},$
(30)

where the volume element $d\hat{\mathcal{W}}$ in velocity space is defined as
$d\hat{\mathcal{W}}=\pi v_{{th}_{s}}^{3}dv_{\parallel}dv_{\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 Eq. (10), in terms of Fourier multipliers in
the following way:

$\displaystyle\mathcal{J}_{0s}\rightarrow
J_{0}\left(\alpha_{s}\right),\qquad\mathcal{J}_{1s}\rightarrow\frac{J_{1}(\alpha_{s})}{\alpha_{s}},$
(31)

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_{{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

$\hat{g}_{s}(\textbf{X}_{s},v_{\parallel},v_{\perp},t)=\hat{\mathcal{F}}_{{eq}_{s}}\sum_{n,m=0}^{\infty}H_{m}(v_{\parallel})L_{n}\left(\frac{v_{\perp}^{2}}{2}\right)f_{mn_{s}}(\textbf{X}_{s},t)$
(32)

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:

$f_{{mn}_{s}}=\frac{1}{n_{0}\sqrt{m!}}\int\text{d}\hat{\mathcal{W}}\,\hat{\mathcal{F}}_{{eq}_{s}}g_{s}H_{m}\left(v_{\parallel}\right)L_{n}\left(\frac{v_{\perp}^{2}}{2}\right).$
(33)

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 gyrocenter electron parallel temperature
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.
(2020) 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

$f_{{mn}_{s}}=-\delta_{m0}\left(G_{1ns}\frac{1}{d_{i}}\frac{2Z_{s}}{\tau_{s}\beta_{e}}\phi+2G_{2ns}d_{i}B_{\parallel}\right)$
(34)

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 Eq. (34), we also introduced the gyroaverage
operators which, as Fourier multipliers, are given by

$\displaystyle
G_{1n}(b_{s})\rightarrow\frac{\mathrm{e}^{-b_{s}/2}}{n!}\left(\frac{b_{s}}{2}\right)^{n},\qquad
n\geq 0,$ (35) $\displaystyle
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),\qquad
n\geq 1,$ (36)

with
$b_{s}=k_{\perp}^{2}v_{{th}_{s}}^{2}/(L\omega_{cs})^{2}=k_{\perp}^{2}d_{i}^{2}(\sigma_{s}\tau_{s}\beta_{e}/(2Z_{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

$\frac{\hat{\omega}}{\hat{k}_{y}}\ll v_{{th}_{s}}$ (37)

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 (22) linearized about an equilibrium
${\phi}^{(0)}={B}_{\parallel}^{(0)}=0,A_{\parallel}^{(0)}=ax$, with a constant
$a$. Such linearization assumed that one can identify, also in real space, the
Fourier multiplier $-k_{\perp}^{2}$ with the operator $\nabla_{\perp}^{2}$ and
then make use of the identity
$J_{0}(\alpha_{s})=\sum_{r=0}^{\infty}(1/(r!)^{2})(v_{\perp}v_{{th}_{s}}/(L\omega_{cs}))^{2r}{\nabla_{\perp}^{2r}}$,
so that $\mathcal{J}_{0s}x\approx x$.The condition (37) 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
equations (1)-(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 (22), with $s=e$.

Concerning $N_{i}$ and $U_{i}$, as above stated, we assume the conditions

$N_{i}=0,\qquad U_{i}=0$ (38)

to hold. This assumption effectively decouples the ion gyrofluid dynamics from
the electron gyrofluid dynamics, leaving Eqs. (1)-(5) a closed system.

The assumption (38), as can be derived from the four-field model in Granier et
al. (2022b) is valid for $\beta_{e}\ll 1$.222With this regard it could be
useful to mention here a misprint in Eq. (2.2) in Granier et al. (2022b),
where $U_{i}$ should have been multiplied by the factor
$2\rho_{s_{\perp}}^{2}/\beta_{e_{\perp}}$. Therefore we expect, in the
gyrokinetic simulations, a departure from the condition (38) as $\beta_{e}$
increases.

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

$\tau\ll 1,$ (39)

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

$\phi^{(0)}=0,\quad\quad A_{\parallel}^{(0)}=A_{\parallel
0}^{eq}/\cosh^{2}\left(x\right),\quad\quad B_{\parallel}^{(0)}=0,$ (40)

with $A_{\parallel 0}^{eq}=1.299$, in order to have
$\mathrm{max}_{x}(B_{y}^{eq}(x))=1$. Due to periodicity assumption in the
simulations333In the AstroGK code, a shape function $S_{h}(x)$ is multiplied
to $A_{\parallel}^{(0)}$ to enforce periodicity (Numata et al., 2010). This
minor difference in the simulation set-up between two models practically
introduces no difference in the following results., the dimensionless
equilibrium $A_{\parallel}^{(0)}$ is replaced by

$A_{\parallel}^{(0)}=\sum_{n=-30}^{n=30}A_{\parallel
0}^{eq}a_{n}\mathrm{e}^{in2\pi x},\qquad$ (41)

where $a_{n}$ are the Fourier coefficients of the function
$1/\cosh^{2}\left(x\right)$.

Note that, in order to satisfy Eqs. (3)-(5) at equilibrium, the corresponding
equilibrium density and parallel electron velocity fluctuations of the
electron gyrocenters are fixed as

$N_{e}^{(0)}=0,\qquad\nabla_{\perp}^{2}A_{\parallel}^{(0)}=G_{10e}U_{e}^{(0)}.$
(42)

Also in the gyrokinetic simulations, in accordance with Eq. (42), the
equilibrium current density is assumed to be entirely due to the parallel
electron velocity (we recall that, as discussed in Sec. 2.3, the gyrokinetic
admits, unlike the gyrofluid model, also a finite parallel ion flow). The
tearing stability parameter (Furth et al., 1963) for the equilibrium is
Porcelli et al. (2002)

$\Delta^{\prime}=2\frac{\left(5-k_{y}^{2}\right)\left(k_{y}^{2}+3\right)}{k_{y}^{2}(k_{y}^{2}+4)^{1/2}}.$
(43)

The equilibrium (40) is tearing unstable when $\Delta^{\prime}>0$, which
corresponds to wave numbers $k_{y}=\pi m/L_{y}<\sqrt{5}$. In the development
of the tearing process, after the saturation of the unstable dominant mode
($m=1$), eventually, a thinning current sheet located at the X point of the
initial tearing configuration will form. In this evolving current sheet,
tearing-mode-like perturbations can develop to form secondary magnetic islands
denoted as plasmoids, when they enter their non linear phase.

The fluid numerical solver SCOPE3D (Solver Collisionless Plasma Equations in
3D) (Granier et al. (2022a)) 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., 2010). Although AstroGK employs
some sophisticated techniques for the treatment of linear terms, it uses
essentially the same pseudo-spectral and temporal schemes.

##  3 Results on the plasmoid onset

An extensive numerical simulation campaign, reported in the tables 1 and 2,
was carried out to study the physical conditions under which plasmoids appear.

No. | $\rho_{s}$ | $\Delta^{\prime}$ | $\beta_{e}$ | $m_{e}/m_{i}$ | $\gamma_{L}$ | $\gamma_{max}$ | Plasmoid
---|---|---|---|---|---|---|---
$1_{GF1}$ | 0.3 | 14.3 | 0.2491 | 0.01 | 0.214 | 0.285 | One small
$1_{GF2}$ | 0.3 | 14.3 | 0.06228 | 0.0025 | 0.225 | 0.322 | No
$1_{F}$ | 0.3 | 14.3 | 0 | 0 | 0.230 | 0.337 | No
$2_{GF1}$ | 0.3 | 29.09 | 0.2491 | 0.01 | 0.211 | 0.342 | One plasmoid
$2_{GF2}$ | 0.3 | 29.09 | 0.1246 | 0.005 | 0.218 | 0.367 | One plasmoid
$2_{GF3}$ | 0.3 | 29.09 | 0.06228 | 0.0025 | 0.231 | 0.378 | One plasmoid
$2_{GF4}$ | 0.3 | 29.09 | 0.00622 | 0.0005 | 0.241 | 0.385 | Several plasmoids
$2_{F}$ | 0.3 | 29.09 | 0 | 0 | 0.242 | 0.386 | Several plasmoids
$3_{GF1}$ | 0.5 | 14.3 | 0.692 | 0.01 | 0.286 | 0.334 | Several plasmoids
$3_{GF2}$ | 0.5 | 14.3 | 0.3460 | 0.005 | 0.310 | 0.383 | Several plasmoids
$3_{F}$ | 0.5 | 14.3 | 0 | 0 | 0.338 | 0.448 | Several plasmoids
$4_{F}$ | 0.06 | 14.3 | 0 | 0 | 0.081 | 0.188 | No
Table 1: Gyrofluid and fluid simulations. No. | $\rho_{s}$ | $\Delta^{\prime}$ | $\beta_{e}$ | $m_{e}/m_{i}$ | $\gamma_{L}$ | $\gamma_{max}$ | Plasmoid
---|---|---|---|---|---|---|---
$1_{GK1}$ | 0.3 | 14.3 | 0.2491 | 0.01 | 0.2245 | 0.308 | One small
$1_{GK2}$ | 0.3 | 14.3 | 0.06228 | 0.0025 | 0.2438 | 0.342 | No
$2_{GK1}$ | 0.3 | 29.09 | 0.2491 | 0.01 | 0.2165 | 0.352 | One large
$2_{GK2}$ | 0.3 | 29.09 | 0.1246 | 0.005 | 0.2267 | 0.389 | One large
$2_{GK3}$ | 0.3 | 29.09 | 0.06228 | 0.0025 | 0.2329 | 0.401 | One large
$3_{GK1}$ | 0.5 | 14.3 | 0.692 | 0.01 | 0.3040 | 0.362 | One
$3_{GK2}$ | 0.5 | 14.3 | 0.3460 | 0.005 | 0.3286 | 0.410 | One
$3_{GK3}$ | 0.5 | 14.3 | 0.1730 | 0.0025 | 0.3472 | 0.453 | One
$4_{GK1}$ | 0.06 | 14.3 | 0.009965 | 0.01 | 0.08617 | 0.207 | No
$4_{GK2}$ | 0.06 | 14.3 | 0.002491 | 0.0025 | 0.08779 | 0.209 | No
Table 2: Gyrokinetic simulations.

Each simulation is identified by a code of the form $p_{F/GF/GK\,r}$, where
$p$ and $r$ are integers and $F$, $GF$ and $GK$ indicate whether the
simulation is carried out in the fluid limit, with the gyrofluid model or with
gyrokinetic model, respectively. For all the simulations, the value of the
electron skin depth is fixed to $d_{e}=0.085$. Simulations with the same
number $p$ are characterized by the same values of $d_{e}$, $\rho_{s}$ and
$\Delta^{\prime}$. For a fixed $p$, different values of the index $r$, on the
other hand, indicate different values of $\beta_{e}$ (and, consequently, of
$m_{e}/m_{i}$), with $\beta_{e}$ decreasing as $r$ increases. Not all the
simulations of Table 1 have a corresponding simulation in Table 2 and
viceversa, although this is the case for most of the simulations. In
particular, we point out that, because gyrokinetic simulations always have a
finite value of $\beta_{e}$, strictly speaking there is no gyrokinetic
counterpart for the fluid simulations, which formally correspond to the
$\beta_{e}\rightarrow 0$ limit.

For all the gyrokinetic simulations, the temperature ratio is set to
$\tau=10^{-3}$, where the ion Larmor radius is $\sqrt{\tau}\rho_{s}$. As
mentioned before, the gyrofluid model assumes $\tau\rightarrow 0$. Therefore,
in both the gyrofluid and gyrokinetic approach, the ion Larmor radius effects
are neglected.

As a first general comment, we observe, by comparing gyrofluid and gyrokinetic
simulations with the same indices $p$ and $r$, that, in terms just of
appearance or absence of plasmoids, gyrofluid simulations agree with the
gyrokinetic ones. Therefore, in this respect, we can conclude that the quasi-
static closure for the electrons and the suppression of ion gyrocenter
fluctuations, do not affect critically the stability of the nonlinear current
sheet. However, as will be discussed in the next Sections, differences appear
in terms of the number and size of plasmoids. In particular, when more than
one plasmoid is observed, this is indicated in table 1, generically, as
’several plasmoids’. The number of plasmoids in the same simulation can indeed
vary in time, as plasmoids can form at different times and pairs of plasmoids
can merge into a single one.

These simulations are compared in the following sections. In particular, with
the set of simulations $p=1$ and $p=2$ the effects of $\beta_{e}$ on the
formation of plasmoids are investigated in Section 3.3. With the set of
simulations $p=3$ and $p=4$ we verify that the regime $\rho_{s}\gg d_{e}$
promotes the formation of plasmoids in Section 3.4.

###  3.1 Growth rates

Before discussing in detail the plasmoid instability, we briefly comment about
the linear growth rate of the tearing mode excited by the perturbation of the
initial equilibria (40). In the tables, we reported the value of the linear
and maximum growth rate of the tearing instability, evaluated measuring the
following quantity at the X-point

$\gamma=\frac{d}{dt}\log\left|A_{\parallel}^{(1)}\left(\frac{\pi}{2},0,t\right)\right|.$
(44)

The two approaches give close growth rate values for small $\beta_{e}$. One
can note a discrepancy on the value of $\gamma_{max}$ when $\beta_{e}$ is
large, although, as shown on Fig. 1, the gyrofluid and gyrokinetic simulations
yield the same dependence of $\gamma_{max}$ on the parameters, with gyrofluid
simulations slightly underestimating the growth rate, in general. This
discrepancy suggests that, for large values of $\beta_{e}$ and during the non-
linear phase, the efficiency of the gyrofluid model to reproduce the
gyrokinetic results becomes limited. As commented in Sec. 2.3, one reason for
this might be the absence of ion gyrocenter density and parallel velocity
fluctuations, which occurs in the gyrofluid model, even for large $\beta_{e}$,
due to the imposed condition (38).

Figure 1: Maximum growth rates of the collisionless tearing mode as a function
of $\beta_{e}$, for the cases $p=2$ and $p=3.$

By comparing the growth rate results, for a fixed mass ratios, of simulations
$p=1-3-4$, we note that increasing $\beta_{e}$ and $\rho_{s}$, as
$\rho_{s}\sim\sqrt{\beta_{e}/2}$, destabilizes the tearing mode. Increasing
these parameters can be seen as fixing the background density, the ion mass
and the guide field amplitude, while increasing the electron temperature. It
was shown numerically in Numata and Loureiro (2015); Granier et al. (2022b)
that, in this latter situation, the linear tearing growth rate is first ruled
by the destabilizing effect of the sonic Larmor radius. However, in cases
where the electron temperature is high enough for the effects of $\rho_{e}$ to
take over those of $\rho_{s}$, the linear growth rate is damped. Here, we find
ourselves in the first case, for which the effects of the sonic Larmor radius
are visibly dominant.

###  3.2 Remarks on the numerical resolution

It is important to anticipate the role of the resolution in this study. In the
forming nonlinear current sheet, tearing modes grow and can become unstable at
different times. The current sheet can therefore be broken by multiple
dominant modes, and the number of plasmoids is highly sensitive to the
resolution used. Given that the fluid simulations $2_{F}$ and $3_{F}$ were
those which allowed the formation of several plasmoids, we carried out
resolution tests with the gyrofluid code on these two simulations to determine
the necessary number of points along $y$, that does not prevent the growth of
large mode numbers.

Table 3 reports the number of visible plasmoids for simulation $2_{F}$ as a
function of the number of points and indicates their order of appearance. The
convergence is reached for a resolution of $2304^{2}$.

$n_{x}\times n_{y}$ | # plasmoids | Comments on the order of appearance
---|---|---
$200\times 120$ | 1 | 1 at the center
$200^{2}$ | 3 | 2 symmetrically with respect to the center then 1 at the center
$200\times 400$ | 3 | 2 symmetrically with respect to the center then 1 at the center
$2304^{2}$ | 7 | 6 symmetrically with respect to the center then 1 at the center
$3400\times 4800$ | 7 | 6 symmetrically with respect to the center then 1 at the center
Table 3: number of visible plasmoids for simulation $2_{F}$ for different
grids.

For $3_{F}$, which is close to marginal stability, a spatial discretization
smaller than $2L_{y}/n_{y}\sim 0.0078$ was needed. Unfortunately, it is not
foreseeable to perform gyrokinetic simulations with such a high resolution.
Therefore, a grid of $256\times 128$ points has been used for all the
gyrokinetic runs. We compared these gyrokinetic simulations to
fluid$/$gyrofluid simulations performed with a nearly identical resolution.
However, since the fluid code is much less demanding in computation time, we
also performed the fluid simulations with grids up to $2304^{2}$ points.

###  3.3 Effect of $\beta_{e}$ on the plasmoid onset

In this Section we present how the $\beta_{e}$ parameter changes the
characteristics of the forming current sheet and promotes the plasmoid
formation. We measure the current sheet aspect ratio using the current density
$j_{\parallel}$. The length is defined such that the current distribution from
the X point to $L_{\rm{cs}}/2$ equals a specific value $\alpha$

$\frac{1}{N}\sum^{N}_{i=1}\left(j_{\parallel}|_{X}-j_{\parallel}(0,i\Delta
y,t)\right)^{2}=\alpha j_{\parallel}|_{X},$ (45)

where $\Delta y$ is the mesh length along $y$ and $N$ indicates the number of
points from $y=0$ to $y=L_{\rm{cs}}/2$. The constant $\alpha$ is taken to
$1/3$ as it gives good measurement of the length of the region with a strong
current. Formula (45) allows to apply a single consistent method for all the
simulations, while taking into account the reduction of the current intensity
along the layer. The width of the current sheet corresponds to the distance
between the two points along $x$ where the value of $j_{\parallel}$ takes the
same value as at the point $(0,L_{\rm{cs}}/2)$.

We focus first on the comparison of the series of simulations for $p=1$,
starting with the higher $\beta_{e}$ case, for which $\beta_{e}=0.2491$. The
contour plots of the parallel electron velocity $u_{e}$ (proportional to the
parallel current density), for the gyrofluid simulation $1_{GF1}$ and of the
current density, $j_{\parallel}$, for the gyrokinetic simulations $1_{GK1}$,
are shown on Fig. 2. Isolines of the magnetic potential, showing the topology
of the magnetic field, are overplotted. Both approaches indicate the formation
of a plasmoid. For the fluid simulation, the aspect ratio is
$A_{\rm{cs_{f}}}=4.90$. In the gyrokinetic case, we measure
$A_{\rm{cs_{k}}}=4.03$. We observe a persistent difference between the value
of the gyrofluid and gyrokinetic aspect ratios, which is explained by the
difference in resolution. However, their evolution according to the parameters
are in agreement.

Figure 2: Top: Contour of the parallel electron velocity $u_{e}$ (proportional
to the parallel current density) for simulation $1_{GF1}$. Bottom: Contour of
the parallel current density $j_{\parallel}$ of simulation $1_{GK1}$. Isolines
of the magnetic potential are superimposed on all the contours.

Figure 3: Top: Contour of the electron velocity $u_{e}$ (proportional to the
parallel current density) for simulation $1_{GF2}$. Bottom: Contour of the
parallel current density $j_{\parallel}$ of simulation $1_{GK2}$. Isolines of
the magnetic potential are superimposed on all the contours.

For the lowest $\beta_{e}$ cases, for which $\beta_{e}=0.06228$, the contour
plots of the simulations $1_{GF2}$ and $1_{GK2}$, are shown on Fig. 3. The two
simulations lead to the formation of a stable current sheet having an aspect
ratio decreasing in time. The maximum aspect ratio is reached when the growth
rate has reached its maximum value and the process enters the saturation
phase. From the gyrofluid simulation we measured a maximum aspect ratio
$A_{\rm{cs_{f}}}=5.11$. In the gyrokinetic case, the aspect ratio is
$A_{\rm{cs_{k}}}=4.14$. The measured aspect ratio are very close to those
obtained for $\beta_{e}=0.2491$, and yet, no plasmoids develop. In this first
series of tests we are at the frontier between stability and instability, and
the role of $\beta_{e}$ seems crucial to switch to an unstable case.

In the series of simulations for $p=2$, the idea is to consider the same
parameters as those for $p=1$ but with a longer forming current sheet. Since
highly unstable primary reconnecting modes favour the formation of extended
secondary current sheets we consider a larger domain size along the y
direction, with $L_{y}=1.4\pi$, that corresponds to $\Delta^{\prime}=29.9$.
The other parameters are kept the same. In this case, even the
small/negligible $\beta_{e}$ simulations become plasmoid unstable. Figure 4
shows the plasmoids and the growth rate evolutions obtained for simulations
$2_{GK1}$ \- $2_{GK3}$. On Fig. 4, in order to better show the plasmoid size,
we show the contour plots of the magnetic potential $A_{\parallel}$. The
magnetic potential contour is shown as the plasmoid reaches its maximum size,
which occurs in the saturation phase of the tearing instability. It can be
seen that increasing $\beta_{e}$ results in larger sized plasmoids, although,
from the aspect ratio measurement, increasing $\beta_{e}$ reduces the aspect
ratio obtained just before the plasmoid onset. Here, increasing $\beta_{e}$
(considering therefore larger mass ratios) seems to have a similar effect to
increasing $\rho_{s}$, and allows the plasmoid instability to grow better in
current sheets whose dimensions are not particularly favourable (low aspect
ratio). A recent result obtained in a 2D, collisionless, fluid model (Grasso
and Borgogno (2022)) has shown that rhos significantly enlarge the spectrum of
the linear unstable reconnecting modes that develop in presence of a sheared
flow and magnetic field

Figure 4: Left: contours of $A_{\parallel}$ for simulations $2_{GK3}$,
$2_{GK2}$ and $2_{GK1}$. Right: growth rate evolution as function of time and
aspect ratio of the forming current sheet as a function of time for the same
simulations.  Figure 5: Growth rate evolution as function of time (top) and
aspect ratio of the forming current sheet as a function of time (bottom) for
the simulations $4_{GK1}$ and $4_{GK2}$.

In comparison, Fig. 5 shows the aspect ratio and the growth rate obtained for
the simulations $4_{GK1}$ and $4_{GK2}$. For this set of simulations, the
effects of $\beta_{e}$ are negligible and the parameters $\rho_{s}$ and
$\Delta^{\prime}$ are smaller than those of simulations $2_{GK1}$ \-
$2_{GK3}$. Nevertheless, despite a very different set of parameters, the two
set of simulations lead to the formation of current sheets whose aspect ratio
is almost identical. Yet, unlike cases $2$, cases $4$ remain stable. If we
assume, as done in the theory by Comisso et al. (2016) and Granier et al.
(2022a), that the growth of the perturbation in the secondary current sheet
does not depend on the initial tearing equilibrium (and that therefore the
initial $\Delta^{\prime}$ has no role in the plasmoid size) then plasmoid
formation is determined only by the values of $\rho_{s}$ and $\beta_{e}$ which
are greater in case $2$.

Figure 6 shows in detail the evolution of the instability for $2_{F}$,
$2_{GF4}$, $2_{GF3}$ and $2_{GF1}$ having the highest resolution. For
negligible $\beta_{e}$ we can see several plasmoids forming in a row, whereas,
for $\beta_{e}>0.06$ we see only a single central plasmoid. The effects of
$\beta_{e}$ eventually prevent the development of large modes inside the
current sheets.

Figure 6: Contour of $u_{e}$ with isolines of $A_{\parallel}$. From top to
bottom panels: : $2_{F}$ ($2304^{2}$), $2_{GF4}$, $2_{GF3}$, $2_{GF1}$
($2304\times 2400$).

###  3.4 Validation of the plasmoid regime for $\rho_{s}\gg d_{e}$

A theory and numerical study developed by Granier et al. (2022a) stated that,
for a current sheet close to marginal stability, the regime $\rho_{s}\gg
d_{e}$ promotes the plasmoid formation. In this Ref., the simulations were
carried out with the fluid model (16) - (17) which assumes a negligible mass
ratio and a negligible $\beta_{e}$. In this subsection, we present a
gyrokinetic validation of these results. In addition to observing a possible
role played by the closure, we also compare the fluid results with those
including a finite mass ratio of $m_{e}/m_{i}=0.005$, and consequently a small
$\beta_{e}$. Moreover, as already recalled, the evolution of ion quantities
such as $N_{i}$ and $U_{i}$, prevented by the gyrofluid model, but present in
the gyrokinetic simulations, might in principle also play a role. Therefore,
in this subsection we focus on the low $\beta_{e}$ regime and compare the
simulation set for $p=3$, for which $\rho_{s}\gg d_{e}$, with the simulation
set for $p=4$, for which $\rho_{s}\ll d_{e}$. These two sets of simulation
lead to the formation of a secondary current layer close to the instability
threshold.

Figure 7 shows the evolution of the instability for the simulations $3_{GK3}$
(lowest $\beta_{e}$ gyrokinetic case of this series) and $3_{F}$. For the two
approaches, the current sheet becomes plasmoid unstable. Also in this case the
resolution plays an important role. With a resolution of $1728^{2}$, three
plasmoids were visible in the simulations $3_{F}$. However, the same fluid
simulation performed with a resolution $500\times 360$ shows only one
plasmoid. Since a resolution higher than that was not foreseeable with the
gyrokinetic code, we used a grid of $256\times 128$ points that allowed to
observe one single plasmoid at the center. In the regime $\rho_{s}\gg d_{e}$,
the current aligns with the magnetic field lines, thus forming a cross shaped
current sheet (Cafaro et al., 1998). This behavior is retrieved by the
gyrokinetic simulation.

Figure 7: Top: contour of the parallel current density $u_{e}=U_{e}$ for
simulation $3_{F}$ ($1728^{2}$). Bottom: Contour of the parallel current
density $j_{\parallel}$ of simulation $3_{GK3}$. Isolines of the magnetic
potential are superimposed on all the color maps.

Figure 8 shows the evolution of the secondary current sheet for the cases
$4_{GK3}$ (lowest $\beta_{e}$ gyrokinetic case of this series) and $4_{F}$.
The current sheet formed in the two frameworks does not follow the
separatrices but remains mainly aligned along $x=0$.

Figure 8: Top: Contour of the parallel current density $u_{e}=U_{e}$ for
simulation $4_{F}$ ($1728^{2}$). Bottom: Contour of the parallel current
density $j_{\parallel}$ of simulation $4_{GK2}$. Isolines of the magnetic
potential are superimposed on all the color maps.

This comparison makes it possible to show that, by simply adding bi-fluid
effects resulting from a large sonic Larmor radius, one can switch from a
marginally stable case to a marginally unstable case. This result, shown by
Granier et al. (2022a) by means of the fluid model, is thus confirmed by
gyrokinetic simulations.

##  4 Energy partition - Similarities and differences between gyrokinetics and
gyrofluid

###  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., 2006;
Schekochihin et al., 2009), normalized by $B_{0}^{2}/(4\pi)$

$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}}_{{eq}_{s}}\mathcal{J}_{0s}\delta{\mathcal{F}}_{s}^{2}+|\nabla_{\perp}A_{\parallel}|^{2}+d_{i}^{2}|B_{\parallel}|^{2}\right)$
(46)

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

$\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}^{\prime},$
(47)

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}^{\prime}$
contains all higher moments of the perturbed distribution function. By
definition,

$\displaystyle\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}\mathcal{J}_{0s}h_{s}^{\prime}$
$\displaystyle=0,$
$\displaystyle\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}v_{\parallel}\mathcal{J}_{0s}h_{s}^{\prime}$
$\displaystyle=0.$ (48)

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

$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}}_{{eq}_{s}}\mathcal{J}_{0s}{h_{s}^{\prime}}^{2}\right)+|\nabla_{\perp}A_{\parallel}|^{2}+d_{i}^{2}|B_{\parallel}^{2}|\right).$
(49)

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 equations (1)
- (5) possesses a conserved Hamiltonian given by

$\displaystyle H_{gf}(N_{e},A_{e})=\frac{1}{2}\int
dxdy\,\left(\rho_{s}^{2}N_{e}^{2}+d_{e}^{2}U_{e}^{2}+|\nabla_{\perp}^{2}A_{\parallel}|^{2}-N_{e}(G_{10e}\phi-\rho_{s}^{2}2G_{20e}B_{\parallel})\right).$
(50)

We remark that, as will be shown in the Appendix, the form of the Hamiltonian
(50), 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 (13) and (14) we can also write the Hamiltonian in terms of
particle variables as follows:

$\displaystyle H_{p}(n_{e},A_{e})=$ $\displaystyle\frac{1}{2}\int
dxdy\,\left(\rho_{s}^{2}n_{e}G_{10e}^{-2}n_{e}+d_{e}^{2}\left(G_{10e}^{-1}u_{e}\right)^{2}+|\nabla_{\perp}^{2}A_{\parallel}|^{2}+d_{i}^{2}|B_{\parallel}|^{2}\right.$
$\displaystyle\left.+n_{e}\left(1-2G_{10e}^{-2}\right)\phi+\phi\left(G_{10e}^{-2}-1\right)\frac{\phi}{\rho_{s}^{2}}\right).$
(51)

When we consider the limit $\beta_{e}$, $m_{e}/m_{i}\rightarrow 0$ the
Hamiltonian of the gyrofluid equations is reduced to

$\displaystyle H_{p}(n_{e},A_{e})=\frac{1}{2}\int
dxdy\,\left(\rho_{s}^{2}n_{e}^{2}+d_{e}^{2}u_{e}^{2}+|\nabla_{\perp}^{2}A_{\parallel}|^{2}+|\nabla_{\perp}\phi|^{2}\right),$
(52)

which is namely the Hamiltonian of the fluid Eqs. (16) - (17). In Eq. (52),
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 $\textbf{E}\times\textbf{B}$ flow energy.

###  4.2 Negligible $\beta_{e}$: fluid vs gyrokinetic

On Fig. 9 we present the comparison between the energy variation of the fluid
case $1_{F}$ and that of the low $\beta_{e}$ gyrokinetic case $1_{GK3}$
($\beta_{e}=0.062$). The variations are defined as
$(1/2)\int\text{d}\mathbb{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.

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 et al. (2013); Numata
and Loureiro (2015)), 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 cinetic energy increases in the fluid case, whereas it is quasi-
constant or decreasing in the gyrokinetic one (Fig. 9). 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_{eq}^{2}$),
the perturbation part ($\tilde{u}_{e}^{2}$) and the cross term
($2\tilde{u}_{e}u_{eq}$). The change of each component is shown on the bottom
panel of Fig. 9. 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.

Figure 9: Top: Time evolution of the energy variations for the cases $1_{F}$
(left) resolution and $1_{GK3}$ (right). Bottom: Change of the parallel
kinetic energy for the same simulations. No plasmoid in this case.

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 Fig.
9 that the main component of the total ion energy consists of the
perpendicular kinetic energy, which is included in the ${h_{s}}^{\prime}$ part
in Eq. (49)444Since the $\textbf{v}_{\perp}$ moments are generally not
orthogonal, we cannot clearly separate each of them., given by

$\frac{1}{2}\int\text{d}x\text{d}yd_{i}^{2}\textbf{u}_{\perp,i}^{2}.$ (53)

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

$\displaystyle
d_{i}\textbf{u}_{\perp,i}=\sqrt{\frac{\tau_{i}\beta_{e}}{2\sigma_{i}}}\frac{1}{n_{0}}\int\text{d}\hat{\mathcal{W}}{\mathcal{F}}_{eq_{i}}\textbf{v}_{\perp}\delta{\mathcal{F}}_{i}$
(54)

Notice that the perpendicular flow holds the identity from the definition

$d_{i}\textbf{u}_{\perp,i}=(-\nabla\phi-\rho_{s}^{2}\nabla\cdot\texttt{p}_{\perp\perp,i})\times\textbf{z}$
(55)

where the perpendicular pressure tensor is given by

$\displaystyle
d_{i}\texttt{p}_{\perp\perp,i}=\tau_{i}\frac{1}{n_{0}}\int\text{d}\hat{\mathcal{W}}{\mathcal{F}}_{eq_{i}}\textbf{v}_{\perp}\textbf{v}_{\perp}\delta{\mathcal{F}}_{i}$
(56)

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

In spite of the closure, the evolution of the energy component (53) is very
similar to that of the $\textbf{E}\times\textbf{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 vs gyrokinetic

When $\beta_{e}$ is very small, the FLR corrections become negligible and the
particle and gyrocenter variables coincide. On the other hand, for non-
negligible $\beta_{e}$, the electron Larmor radius becomes finite and the
relations (11) and (12) allow us to relate the density and parallel velocity
of the particles to those of the gyrocenters. On Fig. 10 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 10: Time evolution of the energy variations for the cases $3_{G}F$ and
$3_{GK}$.

In the plot referring to the gyrofluid energy, we show the variation of both
the particles and gyrocenters energy. For instance, the curve referring to
”Kin∥e” corresponds to the variation of $(1/2)\int dxdyd_{e}^{2}u_{e}^{2}$,
which is comparable to the second term of the gyrokinetic energy (49). The one
referring to ”Gyrocenter Kin∥e” corresponds to the variation of $(1/2)\int
dxdyd_{e}^{2}U_{e}^{2}$. By increasing $\beta_{e}$, the difference between the
variation of the energy of the gyrocenters and that of the particles broadens.
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 already
mentioned, 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 Fig. 10
which shows that, in the gyrokinetic case, for sufficiently large $\beta_{e}$,
the ion fluid is actually accelerated along the $\mathbb{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$.

##  5 Conclusion

In this work, we have numerically investigated the plasmoid formation
employing both gyrofluid and gyrokinetic simulations, assuming a finite, but
small $\beta_{e}$. This made it possible to show that the effect of finite
$\beta_{e}$, associated to finite electron Larmor radius effects, promotes the
plasmoid growth. These results contribute to shed light on collisionless
reconnection mediated by the plasmoid instability, and in particular on the
role of the effects present at the electron scale.

This work showed the ability of the reduced gyrofluid model to achieve
relevant new insights into current-sheet stability and magnetic reconnection.
In particular, predictions on marginal stability on current sheets, obtained
by Granier et al. (2022a) in the fluid limit, were confirmed by gyrokinetic
simulations. It also indicates that the fluid and gyrofluid models make it
possible to obtain accurate results in short computational times.

The comparison between the gyrofluid and the gyrokinetic models reveals key
similarities and differences between the two frameworks, which gives insight
into the important underlying physical effects. Indeed, the adopted
gyrokinetic model is a $\delta f$ model from which the gyrofluid model can be
derived with appropriate approximations and closure hypotheses. This allowed
to directly identify possible limitations of the closures applied to the
gyrofluid moments, that distinguish the gyrofluid model from its gyrokinetic
parent model. We therefore presented the impact of the closure on the
distribution and conversion of energy during reconnection. The closure, which
does not allow for parallel temperature fluctuations, implies that the energy
must be converted into fluctuations of density and parallel velocity of the
electrons. This is not in agreement with the gyrokinetic simulations, but does
not seem to interfere with the formation of plasmoids. In particular, for
relatively small but finite $\beta_{e}$, the hypothesis of absent parallel ion
motion made in the gyrofluid framework is valid and does not affect the
plasmoid instability. The gyrokinetic perpendicular ion velocity is well
represented by the fluid $\textbf{E}\times\textbf{B}$ velocity. On the other
hand, gyrokinetic simulations show a large fraction of magnetic energy
transferred to fluctuations of higher order moments.

### Acknowledgements

The work of RN was partly supported by JSPS KAKENHI Grant Number 22K03568. The
numerical simulations were performed using the EUROfusion high performance
computer Marconi Fusion and Galileo100 hosted at CINECA (project FUA35-FKMR
and IsC86 MR-EFLRA), the computing facilities provided by the Mesocentre
SIGAMME hosted by the Observatoire de la Côte d’Azur, as well as JFRS-1
supercomputer system at the Computational Simulation Centre of the
International Fusion Energy Research Centre (IFERC-CSC) at the Rokkasho Fusion
Institute of QST (Aomori, Japan) and on the ’Plasma Simulator’ (NEC SX-Aurora
TSUBASA) of NIFS with the support and under the auspices of the NIFS
Collaboration Research program (NIFS22KISS019).

### Appendix. Comparison between the gyrofluid isothermal and the quasi-static
closure

In this Appendix we first show how the gyrofluid Hamiltonian (50) can be
obtained from the gyrokinetic Hamiltonian (46) by applying the quasi-static
closures and the assumptions described in Sec. 2.3. Subsequently, we compare
with the gyrofluid Hamiltonian obtained by applying what we refer to as
gyrofluid isothermal closure.

The conserved energy (46) of the $\delta f$ gyrokinetic model used for the
comparison in this paper can be expressed in terms of the gyrocenter perturbed
distribution function

$g_{s}=\delta{\mathcal{F}}_{s}+\frac{1}{d_{i}}\frac{2Z_{s}}{\tau_{s}\beta_{e}}\left(\phi-\left\langle\phi-
d_{i}\sqrt{\frac{\tau_{s}\beta_{e}}{2\sigma_{s}}}\textbf{v}_{\perp}\cdot\textbf{A}_{\perp}\right\rangle_{\textbf{X}_{s}}\right),$
(57)

where $\langle\ \rangle_{\textbf{X}_{s}}$ denotes the gyroaverage at constant
guiding center $\textbf{X}_{s}$ and $\textbf{A}_{\perp}$ is the magnetic
vector potential associated with parallel magnetic perturbations, so that
$\nabla\times\textbf{A}_{\perp}=B_{\parallel}\mathbf{z}$. In this way we
obtain

$\displaystyle
W(\delta{\mathcal{F}}_{e},\delta{\mathcal{F}}_{i})=H_{gy}(g_{e},g_{i})$
$\displaystyle=\frac{1}{2}\int\text{d}x\text{d}y\sum_{s}\frac{1}{n_{0}}\int\text{d}\hat{\mathcal{W}}\hat{\mathcal{F}}_{{eq}_{s}}\left(\frac{\tau_{s}\beta_{e}}{2}g_{s}^{2}\right.$
$\displaystyle\left.+\frac{Z_{s}}{d_{i}}g_{s}\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}\right)\right)$
(58)

Note that, in Eq. (58), we already took the spatial 2D limit, in order to
directly obtain the gyrofuid Hamiltonian.

The gyrocenter perturbed distribution function can be developed as a series of
its gyrocenter moments using Hermite and Laguerre polynomials. Here, we will
retain only the first two moments of the hierarchy for the two species, and
apply two different closures, namely the quasi-static closure Tassi et al.
(2020), and a gyrofluid isothermal closure, where the perpendicular and
parallel gyrocenter temperature fluctuations $T_{\perp s}=T_{\parallel s}=0$,
are set equal to zero (as all the other higher order gyrofluid moments). Such
kind of closure is applied, for instance, by Scott (2010), although in this
Reference, gyrocenter temperature, as well as heat flux fluctuations, are
retained and all the other higher order moments are set equal to zero.

For the quasi static closure, the expansion of the gyrocenter perturbed
distribution functions for the two species are given by

$\displaystyle g_{e}$
$\displaystyle=d_{i}N_{e}+\sqrt{\frac{2\sigma_{e}}{\beta_{e}}}d_{i}v_{\parallel}U_{e}-\sum_{n=1}^{\infty}L_{n}\left(\frac{v_{\perp}^{2}}{2}\right)\left(G_{1ne}\frac{1}{d_{i}}\frac{2Z_{e}}{\beta_{e}}\phi+2G_{2ne}d_{i}B_{\parallel}\right),$
(59) $\displaystyle g_{i}$
$\displaystyle=-\sum_{n=1}^{\infty}L_{n}\left(\frac{v_{\perp}^{2}}{2}\right)\left(G_{1ni}\frac{1}{d_{i}}\frac{2}{\tau_{i}\beta_{e}}\phi+2G_{2ni}d_{i}B_{\parallel}\right)$
(60)

The difference between electron and ion treatments in Eqs. (59) and (60), is
clearly due to the assumption (38). We mention that, by retaining, in Eq.
(60), also ion gyrocenter density and parallel velocity fluctuations, and
applying the same procedure described in the following, one can derive the
energy of the 4-field Hamiltonian gyrofluid model described by Granier et al.
(2022b).

In the case of the gyrofluid isothermal closure the truncated expansion simply
gives

$\displaystyle g_{e}$
$\displaystyle=d_{i}N_{e}+\sqrt{\frac{2\sigma_{e}}{\beta_{e}}}d_{i}v_{\parallel}U_{e},$
(61) $\displaystyle g_{i}$ $\displaystyle=0.$ (62)

To simplify the infinite sums in (59) and (60), we can make use of the
following relations (Szegö, 1975)

$\displaystyle J_{0}(\alpha_{s})$
$\displaystyle=e^{-b_{s}/2}\sum_{n=0}^{\infty}\frac{L_{n}\left(v_{\perp}^{2}/2\right)}{n\text{!}}\left(\frac{b_{s}}{2}\right)^{n},$
(63) $\displaystyle 2\frac{J_{1}(\alpha_{s})}{\alpha_{s}}$
$\displaystyle=e^{-b_{s}/2}\sum_{n=0}^{\infty}\frac{L_{n}^{(1)}\left(v_{\perp}^{2}/2\right)}{(n+1)\text{!}}\left(\frac{b_{s}}{2}\right)^{n},$
(64)

where $L_{n}^{(1)}$ are associated Laguerre polynomials and of the expression
of the operators (35)-(36). We therefore obtain the following equality, that
can be injected in (59) and (60),

$\displaystyle\sum_{n=1}^{\infty}L_{n}\left(\frac{v_{\perp}^{2}}{2}\right)\left(G_{1ns}\frac{1}{d_{i}}\frac{2Z_{s}}{\tau_{s}\beta_{e}}\phi+2G_{2ns}d_{i}B_{\parallel}\right)$
$\displaystyle=\left(J_{0}(\alpha_{s})-G_{10s}\right)\frac{1}{d_{i}}\frac{2Z_{s}}{\tau_{s}\beta_{e}}\phi+\left(v_{\perp}^{2}\frac{J_{1}(\alpha_{s})}{\alpha_{s}}-2G_{20s}\right)d_{i}B_{\parallel}.$
(65)

We can now reduce the gyrokinetic Hamiltonian (58) to gyrofluid ones by
injecting the two sets of truncated perturbed gyrocenter distribution
functions. In the quasi-static case, thanks to the relation (65), all the
contributions involving $G_{1n_{s}}$ and $G_{2n_{s}}$, with $n\geq 1$, cancel.
It turns out that the two resulting gyrofluid Hamiltonians can be written in
an identical form, which corresponds to the following

$\displaystyle H(N_{e},A_{e})$
$\displaystyle=\frac{1}{2}\int\text{d}x\text{d}y\,\left(\rho_{s}^{2}N_{e}^{2}-A_{e}\mathcal{L}_{U_{e}}(A_{e})-N_{e}(G_{10e}\mathcal{L}_{\phi}(N_{e})\right.$
$\displaystyle~{}~{}~{}\left.-\rho_{s}^{2}2G_{20e}\mathcal{L}_{B}(N_{e}))\right),$
(66)

where we recall that $A_{e}=G_{10e}A_{\parallel}-d_{e}^{2}U_{e}$ and the
linear operators $\mathcal{L}_{U_{e}}$, $\mathcal{L}_{\phi}$ and
$\mathcal{L}_{B}$ are given by

$\displaystyle
B_{\parallel}=\mathcal{L}_{B}(N_{e}),\qquad\phi=\mathcal{L}_{\phi}(N_{e}),$
(67) $\displaystyle U_{e}=\mathcal{L}_{U_{e}}(A_{e}).$ (68)

through the quasi-neutrality relation and the two components of Ampère’s law.
The expression (66) coincides, up to integration by parts, to the Hamiltonian
(50).

Evidently, the quasi-static quasineutrality equation and Ampère’s law will
differ from the isothermal ones. Therefore, the total conserved energy are
actually evolving differently and the operators $\mathcal{L}_{B}$,
$\mathcal{L}_{\phi}$, $\mathcal{L}_{U_{e}}$ are closure-dependent operators.

For instance, the explicit form of the quasi-neutrality relation, in the
quasi-static case, writes

$\displaystyle-
G_{10e}N_{e}+\nabla_{\perp}^{2}\phi+(G_{10e}^{2}-1)\frac{\phi}{\rho_{s}^{2}}+(1-G_{10e}2G_{20e})B_{\parallel}=0,$
(69)

while in the case of the gyrofluid isothermal closure, we have

$-G_{10e}N_{e}+\nabla_{\perp}^{2}\phi+(\Gamma_{0e}-1)\frac{\phi}{\rho_{s}^{2}}+(1-\Gamma_{0e}-\Gamma_{1e})B_{\parallel}=0,$
(70)

and where the $\Gamma_{0,1e}$ operators are defined in Fourier space in the
following way

$\Gamma_{0e}\rightarrow
I_{0}\left(k_{\perp}^{2}\frac{\beta_{e}}{2}d_{e}^{2}\right)e^{-k_{\perp}^{2}\frac{\beta_{e}}{2}d_{e}^{2}},\qquad\Gamma_{1e}\rightarrow
I_{1}\left(k_{\perp}^{2}\frac{\beta_{e}}{2}d_{e}^{2}\right)e^{-k_{\perp}^{2}\frac{\beta_{e}}{2}d_{e}^{2}},$
(71)

where $I_{n}$ are the modified Bessel functions of order $n$.

A first comment is that the two closures, and consequently the two sets of
static equations are in fact identical if we assume

$G_{1ne}=G_{2ne}=0,\quad\text{for}\quad n\geq 1,$ (72)

which gives the approximations

$\Gamma_{0e}(b_{e})^{1/2}=G_{10e},\qquad(\Gamma_{0e}(b_{e})-\Gamma_{1e}(b_{e}))^{1/2}=2G_{10e}G_{20e}.$
(73)

On the other hand, the relations (73) can be interpreted in a different way,
i.e. considering the exact expressions (71) for $\Gamma_{0e}$ and
$\Gamma_{1e}$, and assuming that the expressions for $G_{10e}$ and $G_{20e}$
can be adapted to match the quasi-neutrality relations following from the two
closures. In this way, from the first relation in Eq. (73), one retrieves, for
the case $s=e$, the approximate expression for the operators $G_{10s}$
introduced by Dorland and Hammett (1993). The advantages of this approach have
been more recently discussed also by Mandell et al. (2018). The above approach
is indeed reminiscent of the approach used by Dorland and Hammett (1993) in
order to find an expression for $G_{10s}$ yielding a better agreement of liner
gyrofluid theory with the linear gyrokinetic theory. A similar approach,
accounting also for $G_{20s}$ in a finite-$\beta$ gyrofluid model, was
followed by Despain (2011).

A second point is that, in the limit of negligible ion and electron Larmor
radius, when considering $\tau_{i}\rightarrow 0$ and $\beta_{e}\rightarrow 0$,
the two sets of static relations become identical and we obtain the same fluid
Hamiltonian

For instance, both (70) and (69) will reduce to the quasi-neutrality relation

$N_{e}=\nabla_{\perp}^{2}\phi,$ (74)

which is going to give rise to the $\textbf{E}\times\textbf{B}$ flow energy in
the fluid Hamiltonian.

## References

  * Bhat and Loureiro [2018] Pallavi Bhat and Nuno F. Loureiro.  Plasmoid instability in the semi-collisional regime.  _Journal of Plasma Physics_ , 84(6):905840607, December 2018.  doi: 10.1017/S002237781800106X.
  * Bhattacharjee et al. [2009] A. Bhattacharjee, Yi-Min Huang, H. Yang, and B. Rogers.  Fast reconnection in high-lundquist-number plasmas due to the plasmoid instability.  _Physics of Plasmas_ , 16(11):112102, 2009.  doi: 10.1063/1.3264103.  URL https://doi.org/10.1063/1.3264103.
  * Biskamp [1986] D. Biskamp.  Magnetic reconnection via current sheets.  _The Physics of Fluids_ , 29(5):1520–1531, 1986\.  doi: 10.1063/1.865670.  URL https://aip.scitation.org/doi/abs/10.1063/1.865670.
  * Brizard [1992] A. Brizard.  Nonlinear gyrofluid description of turbulent magnetized plasmas.  _Phys. Fluids B_ , 4:1213–1228, 1992.
  * Burch et al. [2016] J. L. Burch, R. B. Torbert, T. D. Phan, L.-J. Chen, T. E. Moore, R. E. Ergun, J. P. Eastwood, D. J. Gershman, P. A. Cassak, M. R. Argall, S. Wang, M. Hesse, C. J. Pollock, B. L. Giles, R. Nakamura, B. H. Mauk, S. A. Fuselier, C. T. Russell, R. J. Strangeway, J. F. Drake, M. A. Shay, Yu. V. Khotyaintsev, P.-A. Lindqvist, G. Marklund, F. D. Wilder, D. T. Young, K. Torkar, J. Goldstein, J. C. Dorelli, L. A. Avanov, M. Oka, D. N. Baker, A. N. Jaynes, K. A. Goodrich, I. J. Cohen, D. L. Turner, J. F. Fennell, J. B. Blake, J. Clemmons, M. Goldman, D. Newman, S. M. Petrinec, K. J. Trattner, B. Lavraud, P. H. Reiff, W. Baumjohann, W. Magnes, M. Steller, W. Lewis, Y. Saito, V. Coffey, and M. Chandler.  Electron-scale measurements of magnetic reconnection in space.  _Science_ , 352(6290):aaf2939, 2016.  doi: 10.1126/science.aaf2939.  URL https://www.science.org/doi/abs/10.1126/science.aaf2939.
  * Cafaro et al. [1998] E. Cafaro, D. Grasso, F. Pegoraro, F. Porcelli, and A. Saluzzi.  Invariants and geometric structures in nonlinear hamiltonian magnetic reconnection.  _Phys. Rev. Lett._ , 80:4430–4433, May 1998.  doi: 10.1103/PhysRevLett.80.4430.  URL https://link.aps.org/doi/10.1103/PhysRevLett.80.4430.
  * Comisso et al. [2013] L. Comisso, D. Grasso, F. L. Waelbroeck, and D. Borgogno.  Gyro-induced acceleration of magnetic reconnection.  _Phys. Plasmas_ , 20(9):092118, 2013.  doi: 10.1063/1.4821840.
  * Comisso et al. [2016] L. Comisso, M. Lingam, Y.-M. Huang, and A. Bhattacharjee.  General theory of the plasmoid instability.  _Physics of Plasmas_ , 23(10):100702, 2016.  doi: 10.1063/1.4964481.  URL https://doi.org/10.1063/1.4964481.
  * Daughton et al. [2009] W. Daughton, V. Roytershteyn, B. J. Albright, H. Karimabadi, L. Yin, and Kevin J. Bowers.  Transition from collisional to kinetic regimes in large-scale reconnection layers.  _Phys. Rev. Lett._ , 103:065004, Aug 2009.  doi: 10.1103/PhysRevLett.103.065004.  URL https://link.aps.org/doi/10.1103/PhysRevLett.103.065004.
  * Daughton and Roytershteyn [2012] William Daughton and Vadim Roytershteyn.  Emerging Parameter Space Map of Magnetic Reconnection in Collisional and Kinetic Regimes.  _Space Science Reviews - SPACE SCI REV_ , 172(1-4):271–282, November 2012.  doi: 10.1007/s11214-011-9766-z.
  * Del Sarto et al. [2003] D. Del Sarto, F. Califano, and F. Pegoraro.  Secondary instabilities and vortex formation in collisionless-fluid magnetic reconnection.  _Phys. Rev. Lett._ , 91:235001, Dec 2003.  doi: 10.1103/PhysRevLett.91.235001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.91.235001.
  * Despain [2011] K. M. Despain.  _Gyrofluid modeling of turbulent, kinetic physics_.  PhD thesis, University of Maryland, 2011.
  * Dorland and Hammett [1993] W. Dorland and G. W. Hammett.  Gyrofluid turbulence models with kinetic effects.  _Phys. Fluids B_ , 5:812–835, 1993.
  * Furth et al. [1963] H.P. Furth, J. Killeen, and M. N. Rosenbluth.  Finite resistivity instabilities of a sheet pinch.  _Phys. Fluids_ , 6:459, 1963.
  * Granier et al. [2021] C. Granier, E. Tassi, D. Borgogno, and D. Grasso.  Impact of electron temperature anisotropy on the collisionless tearing mode instability in the presence of a strong guide field.  _Physics of Plasmas_ , 28(2):022112, 2021.  doi: 10.1063/5.0037227.  URL https://doi.org/10.1063/5.0037227.
  * Granier et al. [2022a] C. Granier, D. Borgogno, L. Comisso, D. Grasso, E. Tassi, and R. Numata.  Marginally stable current sheets in collisionless magnetic reconnection.  _Phys. Rev. E_ , 106:L043201, Oct 2022a.  doi: 10.1103/PhysRevE.106.L043201.  URL https://link.aps.org/doi/10.1103/PhysRevE.106.L043201.
  * Granier et al. [2022b] C. Granier, D. Borgogno, D. Grasso, and E. Tassi.  Gyrofluid analysis of electron $\beta_{e}$ effects on collisionless reconnection.  _Journal of Plasma Physics_ , 88(1):905880111, 2022b.  doi: 10.1017/S0022377822000010.
  * Grasso et al. [2001] D. Grasso, F. Califano, F. Pegoraro, and F. Porcelli.  Phase mixing and island saturation in hamiltonian reconnection.  _Phys. Rev. Lett._ , 86:5051–5054, May 2001.  doi: 10.1103/PhysRevLett.86.5051.  URL https://link.aps.org/doi/10.1103/PhysRevLett.86.5051.
  * Grasso and Borgogno [2022] Daniela Grasso and Dario Borgogno.  _Fluid models for collisionless magnetic reconnection_.  2053-2563. IOP Publishing, 2022.  ISBN 978-0-7503-3559-1.  doi: 10.1088/978-0-7503-3559-1ch6.  URL https://dx.doi.org/10.1088/978-0-7503-3559-1ch6.
  * Howes et al. [2006] Gregory G. Howes, Steven C. Cowley, William Dorland, Gregory W. Hammett, Eliot Quataert, and Alexander A. Schekochihin.  Astrophysical gyrokinetics: Basic equations and linear theory.  _The Astrophysical Journal_ , 651(1):590–614, nov 2006.  doi: 10.1086/506172.  URL https://doi.org/10.1086/506172.
  * Ji and Daughton [2011] Hantao Ji and William Daughton.  Phase diagram for magnetic reconnection in heliophysical, astrophysical, and laboratory plasmas.  _Physics of Plasmas_ , 18(11):111207, 2011.  doi: 10.1063/1.3647505.  URL https://doi.org/10.1063/1.3647505.
  * Kunz et al. [2015] M. W. Kunz, A. A. Schekochihin, C. H. K. Chen, I. G. Abel, and S. C. Cowley.  Inertial-range kinetic turbulence in pressure-anisotropic astrophysical plasmas.  _Journal of Plasma Physics_ , 81(5):325810501, 2015.  doi: 10.1017/S0022377815000811.
  * Loureiro and Uzdensky [2015] N F Loureiro and D A Uzdensky.  Magnetic reconnection: from the sweet–parker model to stochastic plasmoid chains.  _Plasma Physics and Controlled Fusion_ , 58(1):014021, nov 2015.  doi: 10.1088/0741-3335/58/1/014021.  URL https://doi.org/10.1088/0741-3335/58/1/014021.
  * Loureiro et al. [2005] N. F. Loureiro, S. C. Cowley, W. D. Dorland, M. G. Haines, and A. A. Schekochihin.  $x$-point collapse and saturation in the nonlinear tearing mode reconnection.  _Phys. Rev. Lett._ , 95:235003, Nov 2005.  doi: 10.1103/PhysRevLett.95.235003.  URL https://link.aps.org/doi/10.1103/PhysRevLett.95.235003.
  * Loureiro et al. [2007] N. F. Loureiro, A. A. Schekochihin, and S. C. Cowley.  Instability of current sheets and formation of plasmoid chains.  _Physics of Plasmas_ , 14(10):100703, 2007.  doi: 10.1063/1.2783986.  URL https://doi.org/10.1063/1.2783986.
  * Loureiro et al. [2013] N. F. Loureiro, A. A. Schekochihin, and A. Zocco.  Fast collisionless reconnection and electron heating in strongly magnetized plasmas.  _Phys. Rev. Lett._ , 111:025002, Jul 2013.  doi: 10.1103/PhysRevLett.111.025002.  URL https://link.aps.org/doi/10.1103/PhysRevLett.111.025002.
  * Mandell et al. [2018] N. R. Mandell, W. Dorland, and M. Landreman.  Laguerre-Hermite pseudo-spectral velocity formulation of gyrokinetics.  _J. Plasma Phys._ , 84:905840108, 2018.
  * Matthaeus and Lamkin [1985] W. H. Matthaeus and S. L. Lamkin.  Rapid magnetic reconnection caused by finite amplitude fluctuations.  _The Physics of Fluids_ , 28(1):303–307, 1985\.  doi: 10.1063/1.865147.  URL https://aip.scitation.org/doi/abs/10.1063/1.865147.
  * Numata and Loureiro [2015] Ryusuke Numata and N. F. Loureiro.  Ion and electron heating during magnetic reconnection in weakly collisional plasmas.  _Journal of Plasma Physics_ , 81(2):305810201, 2015.  doi: 10.1017/S002237781400107X.
  * Numata et al. [2010] Ryusuke Numata, Gregory G Howes, Tomoya Tatsuno, Michael Barnes, and William Dorland.  Astrogk: Astrophysical gyrokinetics code.  _Journal of computational physics_ , 229(24):9347 – 9372, 2010.  ISSN 0021-9991.
  * Numata et al. [2011] Ryusuke Numata, William Dorland, Gregory G. Howes, Nuno F. Loureiro, Barrett N. Rogers, and Tomoya Tatsuno.  Gyrokinetic simulations of the tearing instability.  _Physics of Plasmas_ , 18(11):112106, 2011.  doi: 10.1063/1.3659035.  URL https://doi.org/10.1063/1.3659035.
  * Olson et al. [2016a] J. Olson, J. Egedal, S. Greess, R. Myers, M. Clark, D. Endrizzi, K. Flanagan, J. Milhone, E. Peterson, J. Wallace, D. Weisberg, and C. B. Forest.  Experimental demonstration of the collisionless plasmoid instability below the ion kinetic scale during magnetic reconnection.  _Phys. Rev. Lett._ , 116:255001, Jun 2016a.  doi: 10.1103/PhysRevLett.116.255001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.116.255001.
  * Olson et al. [2016b] Joseph Olson, J. Egedal, S. Greess, R. Myers, M. Clark, Doug Endrizzi, Ken Flanagan, Jason Milhone, E. Peterson, John Wallace, D. Weisberg, and C. Forest.  Experimental demonstration of the collisionless plasmoid instability below the ion kinetic scale during magnetic reconnection.  _Physical Review Letters_ , 116, 06 2016b.  doi: 10.1103/PhysRevLett.116.255001.
  * Phan et al. [2018] T. Phan, Jonathan Eastwood, Michael Shay, J. Drake, Bengt Sonnerup, M. Fujimoto, P. Cassak, Marit Oieroset, J. Burch, Roy Torbert, A. Rager, J. Dorelli, Daniel Gershman, Craig Pollock, Prayash Sharma Pyakurel, Colby Haggerty, Y. Khotyaintsev, B. Lavraud, Yoshifumi Saito, and Werner Magnes.  Electron magnetic reconnection without ion coupling in earth’s turbulent magnetosheath.  _Nature_ , 557, 05 2018.
  * Porcelli et al. [2002] F Porcelli, D Borgogno, F Califano, D Grasso, M Ottaviani, and F Pegoraro.  Recent advances in collisionless magnetic reconnection.  _Plasma Physics and Controlled Fusion_ , 44(12B):B389–B405, nov 2002.  doi: 10.1088/0741-3335/44/12b/327.  URL https://doi.org/10.1088/0741-3335/44/12b/327.
  * Sarto et al. [2011] D Del Sarto, C Marchetto, F Pegoraro, and F Califano.  Finite larmor radius effects in the nonlinear dynamics of collisionless magnetic reconnection.  _Plasma Physics and Controlled Fusion_ , 53(3):035008, jan 2011.  doi: 10.1088/0741-3335/53/3/035008.  URL https://dx.doi.org/10.1088/0741-3335/53/3/035008.
  * Schekochihin et al. [2009] A. A. Schekochihin, S. C. Cowley, W. Dorland, G. W. Hammett, G. G. Howes, E. Quataert, and T. Tatsuno.  Astrophysical Gyrokinetics: Kinetic and Fluid Turbulent Cascades in Magnetized Weakly Collisional Plasmas.  _The Astrophys. J. Suppl. Series_ , 182(1):310–377, May 2009.  doi: 10.1088/0067-0049/182/1/310.
  * Schep et al. [1994] T. J. Schep, F. Pegoraro, and B. N. Kuvshinov.  Generalized two-fluid theory of nonlinear magnetic structures.  _Phys. Plasmas_ , 1:2843–2851, 1994.
  * Scott [2010] B.D. Scott.  Derivation via free energy conservation constraints of gyrofluid equations with finite-gyroradius electromagnetic nonlinearities.  _Phys. Plasmas_ , 17:102306, 2010.
  * Shibata and Tanuma [2001] Kazunari Shibata and Syuniti Tanuma.  Plasmoid-induced-reconnection and fractal reconnection.  _Earth, Planets and Space_ , 53:473–482, June 2001.  doi: 10.1186/BF03353258.
  * Szegö [1975] G Szegö.  _Orthogonal Polynomials_.  (Providence, RI: American Mathematical Society), 1975.
  * Tassi et al. [2018] E. Tassi, D. Grasso, D. Borgogno, T. Passot, and P. L. Sulem.  A reduced landau-gyrofluid model for magnetic reconnection driven by electron inertia.  _Journal of Plasma Physics_ , 84(4):725840401, 2018.  doi: 10.1017/S002237781800051X.
  * Tassi et al. [2020] E. Tassi, T. Passot, and P.L. Sulem.  A Hamiltonian gyrofluid model based on a quasi-static closure.  _J. Plasma Phys._ , 86:835860402, 2020.
  * Uzdensky et al. [2010a] D. A. Uzdensky, N. F. Loureiro, and A. A. Schekochihin.  Fast magnetic reconnection in the plasmoid-dominated regime.  _Phys. Rev. Lett._ , 105:235002, Dec 2010a.  doi: 10.1103/PhysRevLett.105.235002.  URL https://link.aps.org/doi/10.1103/PhysRevLett.105.235002.
  * Uzdensky et al. [2010b] D. A. Uzdensky, N. F. Loureiro, and A. A. Schekochihin.  Fast magnetic reconnection in the plasmoid-dominated regime.  _Phys. Rev. Lett._ , 105:235002, Dec 2010b.  doi: 10.1103/PhysRevLett.105.235002.  URL https://link.aps.org/doi/10.1103/PhysRevLett.105.235002.