# Multi-Fluid Simulations of Upper Chromospheric Magnetic Reconnection with
Helium-Hydrogen mixture

Q. M. Wargnier  Lockheed Martin Solar & Astrophysics Laboratory, 3251 Hanover
St, Palo Alto, CA 94304, USA  Bay Area Environmental Research Institute, NASA
Research Park, Moffett Field, CA 94035, USA  J. Martínez-Sykora  Lockheed
Martin Solar & Astrophysics Laboratory, 3251 Hanover St, Palo Alto, CA 94304,
USA  Bay Area Environmental Research Institute, NASA Research Park, Moffett
Field, CA 94035, USA  Rosseland Centre for Solar Physics, University of Oslo,
P.O. Box 1029 Blindern, NO0315 Oslo, Norway  V. H. Hansteen  Lockheed Martin
Solar & Astrophysics Laboratory, 3251 Hanover St, Palo Alto, CA 94304, USA
Bay Area Environmental Research Institute, NASA Research Park, Moffett Field,
CA 94035, USA  Rosseland Centre for Solar Physics, University of Oslo, P.O.
Box 1029 Blindern, NO0315 Oslo, Norway  Institute of Theoretical Astrophysics,
University of Oslo, P.O. Box 1029 Blindern, N-0315 Oslo, Norway  B. De Pontieu
Lockheed Martin Solar & Astrophysics Laboratory, 3251 Hanover St, Palo Alto,
CA 94304, USA  Rosseland Centre for Solar Physics, University of Oslo, P.O.
Box 1029 Blindern, NO0315 Oslo, Norway  Institute of Theoretical Astrophysics,
University of Oslo, P.O. Box 1029 Blindern, N-0315 Oslo, Norway

###### Abstract

Our understanding of magnetic reconnection (MR) under chromospheric conditions
remains limited. Recent observations have demonstrated the important role of
ion-neutral interactions in the dynamics of the chromosphere. Furthermore, the
comparison between spectral profiles and synthetic observations of
reconnection events suggest that current MHD approaches appear to be
inconsistent with observations. First, collisions and multi-thermal aspects of
the plasma play a role in these regions. Second, hydrogen and helium
ionization effects are relevant to the energy balance of the chromosphere.

This work investigates multi-fluid multi-species (MFMS) effects on MR in
conditions representative of the upper chromosphere using the multi-fluid
Ebysus code. We compare an MFMS approach based on a helium-hydrogen mixture
with a two-fluid MHD model based on hydrogen only. The simulations of MRs are
performed in a Lundquist number regime high enough to develop plasmoids and
instabilities. We study the evolution of the MR and compare the two approaches
including the structure of the current sheet and plasmoids, the decoupling of
the particles, the evolution of the heating mechanisms, and the composition.

The presence of helium species leads to more efficient heating mechanisms than
the two-fluid case. This scenario, which is out of reach of the two-fluid or
single-fluid models, can reach transition region temperatures starting from
upper chromospheric thermodynamic conditions, representative of a quiet Sun
scenario. The different dynamics between helium and hydrogen species could
lead to chemical fractionation and, under certain conditions, enrichment of
helium in the strongest outflows. This could be of significance for recent
observations of helium enrichment in the solar wind in switchbacks and CMEs.

Magnetohydrodynamics (MHD) —numerical methods — magnetic reconnection — Sun:
atmosphere — Sun: chromosphere

††journal: ApJ

##  1 Introduction

Several plasma state transitions are found in the chromosphere. Indeed, the
plasma can be partially or fully ionized, weakly or strongly magnetized,
weakly or strongly collisional, all depending on the region considered (see
Vernazza & Mason, 1978; Martínez-Sykora et al., 2020; Ni et al., 2020).
Furthermore, the solar atmosphere is a complex environment composed of many
species with a large range of ionization levels. In partially ionized plasmas,
charged species interact with other charged and neutral particles through
collisions. Recent studies have suggested that ion-neutral interaction effects
play an important role in the dynamics of the low solar atmosphere and its
energy balance (see Martínez-Sykora et al., 2015; Ballester et al., 2018;
Soler & Ballester, 2022, and references therein). Recent observations show
that many dynamic and heating events, including through episodic magnetic
reconnection, occur at chromospheric and transition region heights De Pontieu
et al. (2009, 2014a); Peter et al. (2014); Hansteen et al. (2017); Guo et al.
(2020).

Magnetic reconnection (MR) is a relaxation of the magnetic field’s topology
created by a local change in the magnetic field connectivity as described by
Parker (1955). This process is due to relevant dissipative effects in
localized regions, called current sheets, where the magnetic field lines are
reconnecting. This phenomenon is continuously present in the solar
chromosphere (Gontikakis, C. et al., 2013; Bharti, L. et al., 2013; Cargill,
2013; Klimchuk, 2015) as revealed by a wide range of observations, including
the Interface Region Imaging Spectrograph (IRIS, De Pontieu et al., 2014b) and
the Swedish 1-m Solar Telescope (SST) telescopes (e.g., Peter et al., 2014;
Guo et al., 2020; Rouppe van der Voort et al., 2017) which provide unique
observations of the dynamic outer solar atmosphere. In particular, recent
observations with IRIS and SST (Rouppe van der Voort et al., 2017; Young,
2018) of strong brightenings in Si iv emission lines, “UV bursts”, reveal
evidence for the formation of plasmoids in the presence of MR events.
Observations with SST reveal the presence of blob-like features at Alfvénic
speed in the Ca ii K line. However, the comparison between spectral profiles
of observed plasmoids and synthetic observations of reconnection suggests that
1- the number of plasmoids is large and 2- the observed motions within the
reconnection site are much more complicated than those predicted by a
classical single-fluid MHD approach, (Nóbrega-Siverio et al., 2017). UV bursts
may span over extended regions from the upper photosphere, the chromosphere to
the transition region (TR), thus, collisional and multi-thermal aspects of the
plasma may play a role and need to be investigated (Hansteen et al., 2019).

On the one hand, several studies try to approximate collisions and multi-fluid
aspects of the plasma by including the generalized Ohm’s Law with the
ambipolar diffusion, see for instance Leake & Arber (2006); Khomenko et al.
(2018); Martínez-Sykora et al. (2020); Rempel & Przybylski (2021). This
approach is generally used to simulate the solar atmosphere and accounts for
ion-neutral interactions while retaining the simplicity of a single-fluid
approach. On the other hand, several two-fluid MHD models (e.g., Alvarez
Laguna et al., 2016; Ni & Lukin, 2018; Popescu Braileanu, B. et al., 2019;
Wójcik, D. et al., 2020; Ni et al., 2020; Pelekhata, M. et al., 2021;
Niedziela, R. et al., 2021; Leake & Linton, 2013) have been constructed to
investigate collisions and multi-thermal aspects of the plasma in the solar
chromosphere. In that framework, the two fluids can interact through
collisions and reactions, while the ions and electrons can interact with the
electromagnetic field. Alternatively, the multicomponent model (e.g., Graille
et al., 2009; Wargnier et al., 2018; Wargnier, 2019) has been derived from
kinetic theory and developed from a multi-scale analysis of the non-
dimensional Boltzmann equations for electrons and heavy species using a
generalized Chapman-Enskog expansion. Note that we refer to heavy species or
fluids to any species except electrons. This model is a sound alternative to
two-fluid models and is suitable for both partially- and fully-ionized
multicomponent plasmas. Moreover, it accounts for the thermal non-equilibrium
between the heavy particles (any ions and neutrals) and the electrons with a
rigorous and detailed description of the dissipative effects based on kinetic
theory of gases.

In Leake et al. (2012); Leake & Linton (2013), a study of MR with a two-fluid
simulation of MR for a weakly ionized reacting plasma, focusing on the solar
chromosphere by considering collisions, thermal non-equilibrium between ions
and neutrals, ionization/recombination, as well as optically thin radiative
losses, was performed. The studies cited above focus on a two-fluid approach
based on a single species with one ionized level of hydrogen. However, the Sun
has additional species (such as helium and ionized metals) all with several
excited and ionized levels with different degree of magnetization and
collisional cross-sections that potentially could change the thermodynamics of
the MR and its topology. Additionally, recent studies show evidence of the
importance of helium ionization in the energy balance of the solar
chromosphere, as shown in Hansteen et al. (1997); Golding et al. (2014, 2016);
Golding, T. P. et al. (2017); Martínez-Sykora et al. (2020). It has been
demonstrated that the helium ionization state is often far from Local
Thermodynamic Equilibrium (LTE) and that models require that both the temporal
state and photo- and collisional-ionization effects are considered in order to
be consistent with observations. Therefore, the most recent multi-fluid
simulations of reconnection lack the self-consistent treatment of the relevant
physical processes involved in the solar atmosphere, such as decoupling
between various types of species, like hydrogen or helium species, including
their respective ionized and excited levels.

In solar chromospheric conditions, and also in other astrophysical
environments (Gosling, 2007; Pucci et al., 2020b), MR processes mostly evolve
in a high Lundquist number regimes. In this regime, current sheets become
rapidly unstable to resistive tearing instabilities which produces plasmoids
(e.g., Bhattacharjee et al., 2009; Bhattacharjee et al., 2010; Pucci et al.,
2020b). Plasmoids are believed to play a major role in speeding up
reconnection by having an influence on the variation of the current sheet
size, as shown by Murtas et al. (2021). Under the formation of plasmoids, the
current sheet breaks into fragments or sections (see Leake et al., 2020) and
the resulting high current densities in each of these sections facilitate a
high reconnection rate, as shown by Leake et al. (2012); Leake & Linton
(2013). Plasmoid formation due to the instability of Sweet-Parker current
sheets has been extensively examined through numerical studies with single-
fluids (Guo et al., 2020; Loureiro & Uzdensky, 2015; Steinolfson & Hoven,
1984; Ugai, 1995; Loureiro et al., 2005; Leake et al., 2020). Recent works
have proved that in fully ionized plasmas it is possible to trigger plasmoids
for a critical Lunquist number $S_{c}$ ranging from $10^{3}$ to $10^{4}$, as
described by Bhattacharjee et al. (2009); Cassak et al. (2009); Huang &
Bhattacharjee (2010); Ni et al. (2010). However, many studies have also
demonstrated that this critical Lundquist number may vary by several orders of
magnitude depending on the initial setup, such as the current sheet
configuration, the amplitudes of viscosity and perturbations, and the plasma
beta coefficient $\beta_{p}$, (see Ni et al., 2012, 2013; Huang et al., 2017).
Additionally, the role of partial ionization on the onset of MR and
development of the resistive tearing instability has been demonstrated Zweibel
et al. (2011); Pucci et al. (2020a); Murtas et al. (2021). However, most of
these studies consider only hydrogen species and do not consider additional
species such as helium. Therefore, we aim to expand the studies of plasmoid
formation and instabilities by including several species.

Recent measurements from SWEAP and FIELD on board of Parker Solar Probe
(Kasper et al., 2015) revealed that Alfvénic magnetic field reversals called
switchbacks (SBs) are associated with a significant increase in alpha particle
abundances at the upper layers of the solar atmosphere Bale et al. (2021).
Even though the SBs origin is not clear and could be associated with different
possible mechanisms (magnetic reconnections, nanoflares, wave dissipation,
footpoint shearing, or turbulence), observed SBs are characterized by strong
periodic spatial modulations of the amplitude of the magnetic fields combined
with large increase of alpha particles at coronal height. This spatial
characteristic scale is consistent with solar super-granulation scales
suggesting that the driver is connected to the super-granular network magnetic
field (Bale et al., 2021; Fargette et al., 2021). These are also distinct
features of slow solar wind Woolley et al. (2021). The chemical fractionation
involving helium species should occur in the chromosphere, where plasma is
partially ionized, suggesting the driver mechanism of these phenomena may be
connected to this region. Addressing the chemical fractionation could help us
understand the switchbacks’ drivers and slow solar wind, and may provide
insight into the enhancement of alpha particles in the upper layers of the
solar atmosphere.

In this work, a magnetic reconnection event in conditions typical of the upper
chromosphere is studied. To do so, we consider a Lundquist number environment
greatly above the critical Lundquist number $S_{c}\approx 4\times 10^{3}$ that
leads to an instability regime characterized by the production of plasmoids,
see Murtas et al. (2021) and references therein. The multi-fluid multi-species
(MFMS) model accounts for multiple species with different ionized and excited
levels, as described by Martínez-Sykora et al. (2020). We aim to study many
fundamental physical processes and assess their impact on the reconnection
process. By focusing on one scenario of a magnetic reconnection event, we
perform two simulations based on two different mixtures: a two-fluid (TF)
hydrogen mixture based on H, H+ and electrons and a multi-fluid multi-species
(MFMS) hydrogen helium mixture based on H, H+, He, He+ and electrons. This
allows highlighting the impact of helium species with different ionized levels
on the reconnection process. Additionally, we assess the impact of the MFMS
and two-fluid approach on the structure of plasmoids and their formation.

The structure of this paper is as follows. In Section 2, the MFMS and TF model
and their assumptions are described. Section 3 briefly describes the numerical
methods considered in Ebysus. The initial conditions of the MRs simulations
are inspired by upper-chromospheric conditions (Section 4). In Section 5, we
analyze the results of the numerical simulations based on MFMS and TF models
at different levels: the structure of the current sheet and plasmoids, the
velocity fields and decoupling between species, the evolution of the heating
processes and composition of the plasma. Finally, Section 6 contains the
discussion and the summary.

##  2 Multi-fluid multi-species and two-fluid governing equations

We briefly summarize the two-fluid (TF) and multi-fluid multi-species (MFMS)
MHD governing equations and their assumptions. Note that all equations are
written in SI units. For MFMS, we focus on a hydrogen-helium mixture
$\mathcal{M}$ such as $\mathcal{M}=\mathcal{H}\cup\\{\mathfrak{e}\\}$, where
$\mathcal{H}=\left\\{\text{H},\text{H}^{+},\text{He},\text{He}^{+}\right\\}$
and $\mathfrak{e}$ denotes the set of heavy particles and electrons,
respectively. Whereas TF has a pure hydrogen mixture
$\mathcal{M}_{\text{TF}}=\mathcal{H}_{\text{TF}}\cup\\{\mathfrak{e}\\}$ where
$\mathcal{H}_{\text{TF}}=\left\\{\text{H},\text{H}^{+}\right\\}$. No excited
levels will be considered in this work. $\mathfrak{n}_{\mathrm{\mathcal{M}}}$
and $\mathfrak{n}_{\mathrm{\mathcal{H}}}$ denote the number of species in
$\mathcal{M}$ and $\mathcal{H}$, respectively.

For clarity and consistency, we use a similar nomenclature used by Ballester
et al. (2018); Khomenko et al. (2014) with minor adjustments. The ionization
states are referred as $\mathcal{I}$, i.e., $\mathcal{I}=0$ denotes neutrals
and $\mathcal{\hat{I}}=\mathcal{I}\geq 1$ ions. The identity of the chemical
species (here hydrogen or helium) is indicated by $\mathfrak{a}$.
Consequently, each set of particles in a given micro-state is described by
$\mathfrak{a}\mathcal{I}\in\mathcal{M}$ or
$\mathfrak{a}\mathcal{I}\in\mathcal{M}_{\text{TF}}$. For electrons, the
notation $\mathfrak{a}\mathcal{I}$ is reduced to just $\\{\mathfrak{e}\\}$. In
both TF and MFMS models, we neglect the electron inertia and assume
quasineutrality. Additionally, we neglect the viscous effects associated with
electrons. Each species has its own temperature $T_{\mathfrak{a}\mathcal{I}}$
and no internal degrees of freedom. Finally, we neglect the effect of gravity.

In the following, the governing equations will be written for the general
case, i.e., for MFMS based on $\mathcal{M}$. The set of equations for TF based
on $\mathcal{M}_{\text{TF}}$ will be presented as a simplification of MFMS.

###  2.1 Continuity Equations for Heavy Particles

The mass density for each type of fluid in a given micro-state
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$ is governed by the continuity equation
in this generic form:

$\displaystyle\partial_{t}\rho_{\mathrm{\mathfrak{a}\mathcal{I}}}+\nabla\cdot\left(\rho_{\mathrm{\mathfrak{a}\mathcal{I}}}\,{\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}}\right)=$
$\displaystyle m_{\mathrm{\mathfrak{a}\mathcal{I}}}$
$\displaystyle(\Gamma^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+\Gamma^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}),\quad\mathfrak{a}\mathcal{I}\in\mathcal{H},$
(1)

where
$\rho_{\mathrm{\mathfrak{a}\mathcal{I}}}=m_{\mathrm{\mathfrak{a}\mathcal{I}}}\,n_{\mathrm{\mathfrak{a}\mathcal{I}}}$
and $m_{\mathrm{\mathfrak{a}\mathcal{I}}}$ are the mass density and particle
mass respectively, $\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}$ is the
velocity, $\Gamma^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}$, and
$\Gamma^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}$ are the mass
transition rate due to recombination or de-excitation, and ionization or
excitation, respectively. For any $\mathfrak{a}\mathcal{I}\in\mathcal{H}$,
these rates are defined as

$\Gamma^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}=\sum_{\mathfrak{a}\mathcal{I}^{\prime}\in\mathcal{H}\neq\mathfrak{a}\mathcal{I}}\begin{cases}n_{\mathrm{\mathfrak{a}\mathcal{I}^{\prime}}}\nu_{\mathrm{\mathfrak{a}\mathcal{I},\mathcal{I}^{\prime}}}^{\text{ion}},&\text{if}\,\mathcal{I}>\mathcal{I}^{\prime}\\\
-n_{\mathrm{\mathfrak{a}\mathcal{I}}}\nu_{\mathrm{\mathfrak{a}\mathcal{I},\mathcal{I}^{\prime}}}^{\text{ion}},&\text{otherwise,}\end{cases}$
(2)

and

$\Gamma^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}=\sum_{\mathfrak{a}\mathcal{I}^{\prime}\in\mathcal{H}\neq\mathfrak{a}\mathcal{I}}\begin{cases}-n_{\mathrm{\mathfrak{a}\mathcal{I}}}\nu_{\mathrm{\mathfrak{a}\mathcal{I},\mathcal{I}^{\prime}}}^{\text{rec}},&\text{if}\,\mathcal{I}>\mathcal{I}^{\prime}\\\
n_{\mathrm{\mathfrak{a}\mathcal{I}^{\prime}}}\nu_{\mathrm{\mathfrak{a}\mathcal{I},\mathcal{I}^{\prime}}}^{\text{rec}},&\text{otherwise,}\end{cases}$
(3)

where
$\nu_{\mathrm{\mathfrak{a}\mathcal{I},\mathcal{I}^{\prime}}}^{\text{ion}}$ and
$\nu_{\mathrm{\mathfrak{a}\mathcal{I},\mathcal{I}^{\prime}}}^{\text{rec}}$ are
the ionization and recombination frequencies between particles with micro-
state $\mathrm{\mathfrak{a}\mathcal{I}}$ and
$\mathrm{\mathfrak{a}\mathcal{I}^{\prime}}\in\mathcal{H}$. In both MFMS and
TF, the frequencies associated with the transition rates $\nu^{rec}$ and
$\nu^{ion}$ include radiative and dielectronic recombination where the rates
have been taken from Aldrovandi & Pequignot (1973). The ionization and three-
body recombination rates have been taken from Janev et al. (1989).
Photoionization is roughly modeled by simply assuming constant rates taken
from Vernazza et al. (1981) for the H-H+ and He-He+ transitions. Note that for
TF the continuity equations presented in Eq. 1 are identical to MFMS by simply
changing the set of heavy particles $\mathcal{H}$ by
$\mathcal{H}_{\text{TF}}$.

###  2.2 Momentum Equations for Heavy Particles

In MFMS, the momentum equation associated with a given species defined by its
microstate $\mathfrak{a}\mathcal{I}\in\mathcal{H}$ is defined as follows:

$\displaystyle\partial_{t}(\rho_{\mathrm{\mathfrak{a}\mathcal{I}}}{\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}})$
$\displaystyle+\nabla\cdot\left(\rho_{\mathrm{\mathfrak{a}\mathcal{I}}}{\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}}\otimes{\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}}+P_{\mathrm{\mathfrak{a}\mathcal{I}}}\mathbf{\mathbb{I}}-\boldsymbol{\tau}_{\mathfrak{a}\mathcal{I}}\right)=$
(4) $\displaystyle
n_{\mathrm{\mathfrak{a}\mathcal{I}}}q_{\mathrm{\mathfrak{a}\mathcal{I}}}\left(\mathbb{E}+{\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}}\wedge\mathbb{B}\right)$
$\displaystyle+\mathbb{R}^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+\mathbb{R}^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+\mathbb{R}_{\mathrm{\mathfrak{a}\mathcal{I}}}^{\text{col}},\quad\mathfrak{a}\mathcal{I}\in\mathcal{H}$

where $q_{\mathfrak{a}\mathcal{I}}$ and $P_{\mathfrak{a}\mathcal{I}}$ are the
ion charge and gas pressure of a given species
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$,
$\boldsymbol{\tau}_{\mathfrak{a}\mathcal{I}}$ is the viscous stress tensor
associated with species $\mathfrak{a}\mathcal{I}\in\mathcal{H}$, $\mathbb{E}$
and $\mathbb{B}$ are the electric and magnetic field,
$\mathbb{R}^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}$ and
$\mathbb{R}^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}$ are the changes
in momentum for species $\mathfrak{a}\mathcal{I}$ due to ionization and
recombination processes respectively,
$\mathbb{R}_{\mathfrak{a}\mathcal{I}}^{\text{col}}$ is the sum of all the
momentum exchange between particles with a given micro state
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$ with any other particles with another
micro state $\mathfrak{a}^{\prime}\mathcal{I}^{\prime}\in\mathcal{M}$.

In TF, the momentum equation written in Eq. 4 slightly differs. We solve for
1- a momentum equation for neutrals (here H) and 2- a momentum equation for
all ionized species that is the sum of the momentum equations for ions (here
$\text{H}^{+}$) and electrons. Therefore, the Lorentz force, corresponding to
the first right-hand side term of Eq. 4, reads $\mathbb{J}\wedge\mathbb{B}$
where $\mathbb{J}$ is the total current density. Additionally, in the momentum
equations for ions, the pressure term in the convective fluxes is simply
replaced by the total pressure of charged particles
$\nabla\cdot\left(\left[P_{\text{H}^{+}}+P_{\mathfrak{e}}\right]\mathbf{\mathbb{I}}\right)$.
Note that in MFMS the term associated with the gradient of pressure of the
electrons is involved in the definition of $\mathbb{E}$.

In this study, we consider plasmas evolving in a high Reynolds number regime
leading to characteristic length scales associated with viscous effects that
are generally much smaller than any chosen grid size. In this context, an
artificial subgrid-scale model of these viscous terms
$\boldsymbol{\tau}_{\mathfrak{a}\mathcal{I}},\,\mathfrak{a}\mathcal{I}\in\mathcal{M}$
or $\mathcal{M}_{\text{TF}}$ has been considered in order to represent shocks
and turbulence. This approach is inspired from Nordlund (1982); Gudiksen et
al. (2011) and has been extended for MFMS and TF models.

In MFMS, the sum of all the momentum exchange for a fluid of a specific
particle with a given micro state $\mathfrak{a}\mathcal{I}\in\mathcal{H}$ is
defined as:

$\displaystyle\mathbb{R}_{\mathfrak{a}\mathcal{I}}^{\text{col}}$
$\displaystyle=\sum_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}\in\mathcal{M}}\mathbb{R}^{\text{col}}_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}\,$
(5)
$\displaystyle=\sum_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}\in\mathcal{M}}m_{\mathfrak{a}\mathcal{I}}\,n_{\mathfrak{a}\mathcal{I}}\nu^{\text{col}}_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}(\mathbb{u}_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}-\mathbb{u}_{\mathfrak{a}\mathcal{I}}),\;$

where
$\nu^{\text{col}}_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}$
is the collision frequency between particle $\mathfrak{a}\mathcal{I}$ and
$\mathfrak{a}^{\prime}\mathcal{I}^{\prime}$. The momentum exchange operators
are symmetric, i.e.,
$\forall\,(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{H}\times\mathcal{M},\;\mathbb{R}^{\text{col}}_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}=-\mathbb{R}^{\text{col}}_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime},\mathfrak{a}\mathcal{I}}$.
Note that Eq. 5 does not include the term associated with the difference in
heat fluxes between species, as described by Hansteen et al. (1997) or in the
13N moment model of Zhdanov (2002). The collision frequencies
$\nu^{\text{col}}_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}$
are described in detail in Wargnier et al. (2022). In this context, the
collisional frequencies are consistent with the definition given by Zhdanov
(2002) in the context of the 13-N moment model for multicomponent plasmas,
calculated from generalized Chapman-Cowling collision integrals. The collision
integrals data for $\mathcal{M}$ and $\mathcal{M}_{\text{TF}}$ have been taken
from Bruno et al. (2010) based on a review of experimental measurements. These
collision integrals include charge-exchange interactions for H-H+ and He-He+
interactions.

###  2.3 Thermal Energy Equations

In both MFMS and TF models, we account for thermal decoupling between all
particles, including electrons. Therefore, we solve for the equations of
thermal energy for each particle $\mathfrak{a}\mathcal{I}\in\mathcal{M}$ or
$\mathcal{M}_{\text{TF}}$.

In MFMS, the energy equations for each heavy particle and electrons read as
follows,

$\displaystyle\partial_{t}e_{\mathrm{\mathfrak{a}\mathcal{I}}}+\nabla\cdot$
$\displaystyle\left(e_{\mathrm{\mathfrak{a}\mathcal{I}}}\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}\right)+P_{\mathrm{\mathfrak{a}\mathcal{I}}}\nabla\cdot\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}=$
(6) $\displaystyle
Q^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+Q^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+Q_{\mathrm{\mathfrak{a}\mathcal{I}}}^{\text{col}}+Q_{\mathrm{\mathfrak{a}\mathcal{I}}}^{\boldsymbol{\tau}},\quad\mathfrak{a}\mathcal{I}\in\mathcal{H}$

$\displaystyle\partial_{t}e_{\mathrm{\mathfrak{e}}}+\nabla\cdot$
$\displaystyle\left(e_{\mathfrak{e}}\mathbb{u}_{\mathfrak{e}}\right)+P_{\mathfrak{e}}\nabla\cdot\mathbb{u}_{\mathfrak{e}}=$
(7) $\displaystyle
Q^{\text{ion}}_{\mathfrak{e}}+Q^{\text{rec}}_{\mathfrak{e}}+Q_{\mathrm{\mathfrak{e}}}^{\text{col}}+Q_{\mathrm{\mathfrak{e}}}^{\mathrm{spitz}},$

where $e_{\mathfrak{a}\mathcal{I}}$ and ${e_{\mathfrak{e}}}$ are the thermal
energies of any heavy particles $\mathfrak{a}\mathcal{I}\in\mathcal{H}$ and
electrons respectively.
$Q^{\text{ion}}_{\mathfrak{a}\mathcal{I}}+Q^{\text{rec}}_{\mathfrak{a}\mathcal{I}}$
and $Q^{\text{ion}}_{\mathfrak{e}}+Q^{\text{rec}}_{\mathfrak{e}}$ are the
heating/cooling term due to the ionization and recombination processes
associated with particles $\mathfrak{a}\mathcal{I}\in\mathcal{H}$ and
electrons respectively, $Q_{\mathfrak{e}}^{\mathrm{spitz}}$ is the heating
term due to the electron thermal conduction along the magnetic field,
$Q_{\mathrm{\mathfrak{a}\mathcal{I}}}^{\boldsymbol{\tau}}$ is the heating term
due to viscous effects associated with particle
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$,
$Q_{\mathfrak{a}\mathcal{I}}^{\text{col}}$ and $Q_{\mathfrak{e}}^{\text{col}}$
are the heating terms due to collisions associated with heavy particle
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$ and electrons respectively. In TF, the
energy equations Eqs. 6 and 7 are obtained by substituting
$\mathfrak{a}\mathcal{I}\in\mathcal{H}_{\text{TF}}$ to
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$.

Note that $Q^{\text{ion}}_{\mathrm{\mathfrak{e}}}$ corresponds to the
optically thin radiative losses term which involve hydrogen and helium
ionization potential coefficients, as described by Leake et al. (2012); Leake
& Linton (2013). The Spitzer term is anisotropic with respect to the magnetic
field and is identical to the definition from Spitzer (1956); Gudiksen et al.
(2011). The heating term due to viscous effects
$Q_{\mathrm{\mathfrak{a}\mathcal{I}}}^{\boldsymbol{\tau}}$ is derived from the
viscous terms in momentum equations
$\boldsymbol{\tau}_{\mathfrak{a}\mathcal{I}}$ leading to a similar expression
introduced by Nordlund (1982); Gudiksen et al. (2011) for a single-fluid MHD
model.

The collisional terms associated with any particle
$\mathfrak{a}\mathcal{I}\in\mathcal{M}$ are defined as

$Q_{\mathfrak{a}\mathcal{I}}^{\text{col}}=\sum_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}\in\mathcal{M}}Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\text{col}}=\sum_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}\in\mathcal{M}}Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\textbf{u},\text{col}}+Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{T,\text{col}}.$
(8)

The first component of Eq. 8 corresponds to the heating produced by the drift
velocity between species
$(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
due to collisions. The second component corresponds to a thermalization
process where temperature $T_{\mathfrak{a}\mathcal{I}}$ relaxes towards
$T_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}$ at a collisional rate
$\nu^{\text{col}}_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}$.
The heating produced by the drift velocity and thermalization terms between
species $\mathfrak{a}\mathcal{I}$ and
$\mathfrak{a}^{\prime}\mathcal{I}^{\prime}$ are defined as

$\begin{split}Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\textbf{u},\text{col}}&=\nu^{\text{col}}_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}\,m_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}n_{\mathfrak{a}\mathcal{I}}|\mathbb{u}_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}-\mathbb{u}_{\mathfrak{a}\mathcal{I}}|^{2},\\\
Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{T,\text{col}}&=\nu^{\text{col}}_{\mathfrak{a}\mathcal{I}\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}\left(\tfrac{3m_{\mathfrak{a}\mathcal{I}}}{m_{\mathfrak{a}\mathcal{I}}+m_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}}\right)n_{\mathfrak{a}\mathcal{I}}k_{B}\left(T_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}-T_{\mathfrak{a}\mathcal{I}}\right).\end{split}$
(9)

Note that for any
$(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
or $\mathcal{M}^{2}_{\text{TF}}$, thermalization terms
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{T,\text{col}}$
do not contribute to the total thermal energy (sum of all the thermal energy
equations Eqs. 6 and 7) of the multi-fluid plasma since these terms are
symmetric with respect to species. Similarly, for any
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$ or $\mathcal{H}_{\text{TF}}$,
ionization/recombination terms of heavy species
$Q^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+Q^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}$
do not contribute to the total thermal energy. However, the heating produced
by drift velocities, optically thin radiative losses and Spitzer terms
contribute to the total thermal energy.

####  2.3.1 Equation of States

In both MFMS and TF, we assume equation of states for electrons and heavy
species following classical ideal gases laws such as

$P_{\mathrm{\mathfrak{a}\mathcal{I}}}=n_{\mathrm{\mathfrak{a}\mathcal{I}}}k_{B}T_{\mathrm{\mathfrak{a}\mathcal{I}}}=(\gamma-1)e_{\mathrm{\mathfrak{a}\mathcal{I}}},\quad\mathfrak{a}\mathcal{I}\in\mathcal{H}\,\,\,\text{or}\,\,\,\mathcal{H}_{\text{TF}}$
(10)

and

$P_{\mathrm{\mathfrak{e}}}=n_{\mathrm{\mathfrak{e}}}k_{B}T_{\mathrm{\mathfrak{e}}}=(\gamma-1)e_{\mathrm{\mathfrak{e}}},$
(11)

where $\gamma=5/3$ is the adiabatic constant taken as identical for all
species.

###  2.4 Generalized Ohm’s law and magnetic Induction Equation

In both MFMS and TF, we compute the electric field from a simplified electron
momentum equation in the asymptotic limit where we assume that 1- the electron
inertia and its time variation is neglected, 2- the mass of electrons is much
smaller than the mass of any heavy particles and 3- the effects associated
with ionization and recombination are neglected. In this framework, in both
MFMS and TF, the Ohm’s law reads

$\displaystyle\mathbb{E}=-\mathbb{u}_{\mathfrak{e}}\wedge\mathbb{B}+\frac{1}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\left[\nabla
P_{\mathfrak{e}}-\mathbb{R}_{\mathfrak{e}}^{\text{col}}\right],$ (12)

where $\mathbb{R}_{\mathfrak{e}}^{\text{col}}$ has the same definition as Eq.
5 by replacing a given heavy species $\mathfrak{a}\mathcal{I}\in\mathcal{H}$
by $\mathfrak{a}\mathcal{I}=\\{\mathfrak{e}\\}$. The electromagnetic field is
not considered as an external force but their evolution is governed by the set
of Maxwell equations. Therefore, combining Eq. 12 with Maxwell-Faraday
equation,

$\partial_{t}\mathbb{B}=-\nabla\wedge\mathbb{E},$ (13)

we obtain the magnetic induction equation.

In MFMS, since several ionized species are considered and the quasi-neutrality
approximation is assumed, the electron velocity can be expressed as a function
of the hydrodynamic velocity of each ion and the total current as follows

$\mathbb{u}_{\mathfrak{e}}=\sum_{\mathfrak{a}\mathcal{I}}\frac{n_{\mathfrak{a}\mathcal{I}}q_{\mathfrak{a}\mathcal{I}}\mathbb{u}_{\mathfrak{a}\mathcal{I}}}{n_{\mathfrak{e}}q_{\mathfrak{e}}}-\frac{\mathbb{J}}{n_{\mathfrak{e}}q_{\mathfrak{e}}}.$
(14)

Note that the contribution of the total current $\mathbb{J}$ of the electron
velocity in the equation of the electric field in Eq. 12, also known as the
Hall term, has not been considered for this study (however, we will
investigate the effect of this term in future studies). In this study, since
the characteristic speed of each particle is much smaller than the speed of
light in vacuum, the Maxwell-Ampere’s law can be simplified by neglecting the
displacement current term, which leads to the following relation of the total
current

$\mathbb{J}=(\nabla\wedge\mathbb{B})/\mu_{0},$ (15)

where $\mu_{0}$ is the vacuum permeability. Note that from Eq. 12 it is
possible to identify one component of the resistivity of the MFMS model as
function of collision integrals, see Section A for further details.

The set of equations associated with MFMS and TF described in Eqs. 1, 4, 6, 7
and 13 is a closed system of equations which can be solved and integrated with
a specific numerical strategy.

##  3 Numerical methods

The two sets of governing equations for MFMS and TF models presented in Eqs.
1, 4, 6, 7, and 13 have been implemented in a code named Ebysus (Martínez-
Sykora et al., 2020), which is an extension of the single-fluid radiative MHD
code Bifrost (Gudiksen et al., 2011). Note that the governing equations can be
rewritten in a more compact form which allows us to distinguish the nature of
the different terms involved, as shown in Appendix B.

Concerning the spatial discretization, a sixth-order differential operator has
been considered for the convective fluxes, as in Bifrost. The conservative
quantities are not co-located in space and are considered on a staggered grid,
thus, a fifth-order interpolation is used to relocate variables as needed.
Since a high-order finite difference scheme has been considered, it is
convenient to numerically stabilize shocks and discontinuities with artificial
numerical terms. In the literature, several type of artificial terms have been
introduced, as described by Wilkins (1980); Cook & Cabot (2005) and references
therein. By inheritance of the Bifrost code, the structure of the artificial
terms are the hyperdiffusive terms described by Nordlund (1982); Gudiksen et
al. (2011). These terms are considered in order to capture shocks and
discontinuities while treating correctly turbulence in a high Reynolds number
regime. In Ebysus the hyperdiffusive terms have been expanded for the multi-
fluid case from those considered in Bifrost since they depend on the velocity
and speeds of sound for each species and the Alfven speed of all ions. The
constant coefficients introduced in these terms have been chosen based on fits
to a wide variety of shock and turbulence test problems.

The integration of the two presented systems requires some specific temporal
integration methods. Indeed, these systems can be seen as convective-diffusive
systems with stiff source terms (see Abdulle & Vilmart, 2013; Duarte et al.,
2013; Descombes & Massot, 2004), as presented in Appendix B. Some of these
source terms (in particular those associated with collisions and
ionization/recombination) are numerically stiff since they are associated with
physical processes where timescales often are much smaller than the convective
timescales. Therefore, these terms are subject to strong CFL constraints
leading to small timesteps and drastically increased computational costs. To
avoid this, our strategy is to consider a first order Lie splitting approach
in order to split the temporal integration on different subsystems. First, the
convective terms have been integrated following a third-order predictor-
corrector Hyman (Hyman et al., 1979). Alternately, it is possible to consider
a third-order explicit Runge-Kutta method. Then, the diffusive term (i.e., the
Spitzer term) has been integrated by solving a non-Fickian (hyperbolic)
diffusion equation as described by (see Rempel, 2016). Note that, in Ebysus,
it is possible to integrate this term following an implicit multi-grid method,
as described in Gudiksen et al. (2011). Finally, the stiff source terms (see
$S^{\text{stiff}}$ in Appendix B) are integrated following an Ordinary
Differential Equation (ODE) solver based on a fifth-order implicit Runge Kutta
method also known as Radau IIA method (see Hairer & Wanner, 1999, 1996). In
order to maintain the solenoidal constraint on the magnetic field, a Hodge
projection method has been used (see Brackbill & Barnes, 1980).

Ebysus allows the introduction of as many species as desired using DIPER
format atom files and the CHIANTI database (see Judge, 2007; Del Zanna, G. et
al., 2015). The number of ionized and excited levels, transitions, and rates
are detailed in these atom files. The code allows one to easily perform
numerical simulations with any chosen mixtures following identical numerical
strategies focusing on the system of equations presented in the previous
section, as long as the atom files associated with all species considered are
provided. The Ebysus code also has the possibility of running a single-fluid
MHD option allowing a comparison with multi-fluid MHD models using the same
numerical code.

##  4 Methodology and initial conditions

In this study, we focus on a magnetic reconnection configuration under
specific thermodynamic conditions that aim to represent the upper solar
chromosphere. For identical initial thermodynamic conditions, we perform a
comparison between TF based on $\mathcal{M}_{\text{TF}}$ (similar to what was
performed by Alvarez Laguna et al., 2016; Leake et al., 2012; Leake & Linton,
2013) and MFMS based on $\mathcal{M}$. Our comparison highlights the impact of
helium on the reconnection process.

###  4.1 Initial reconnection setup

The reconnection setup is inspired from Leake et al. (2012); Leake & Linton
(2013), i.e., a magnetic field configuration forming a Harris current sheet
placed in a constant thermodynamic and plasma beta environment. In this work,
we consider an initial temperature of $T_{0}=16\,000$ K and a total number
density $n_{0}=7.47\times 10^{9}$ $\text{cm}^{-3}$. Then, the population
associated with each species involved in $\mathcal{M}$ and
$\mathcal{M}_{\text{TF}}$ is calculated. In both MFMS and TF, the plasma beta
coefficient is calculated from the total pressure of the plasma
$P_{0}=n_{0}k_{B}T_{0}$ and is defined as $\beta_{p}=0.2$. We have considered
a 2.5D domain based on a uniform grid, with $800\times 800$ grid points where
each grid cell size is$dy=dz=5.333$ km. In this configuration, the whole
computational domain represents a 2.5D box of size $L_{0}=4.2664$ Mm. Open
boundary conditions have been considered in all directions for all the
variables and fields with linear extrapolations while minimizing strong
gradients. In order to mitigate reflections, the last $0.1332$ Mm grid points
inside the domain are dedicated for extrapolating the fields and damping waves
that could eventually disturb the reconnection site.

The Harris current sheet, as described by Leake et al. (2012); Leake & Linton
(2013), is defined as follows

$\mathbf{A_{x}}=-\frac{B_{0}}{2}\lambda_{\psi}\log\left(\cosh(\frac{y}{\lambda_{\psi}})\right)\mathbf{e_{x}},\;\mathbf{B}=\left[\nabla\wedge\mathbf{A_{x}}\right],$
(16)

where $\lambda_{\psi}$ is the initial thickness of the current sheet defined
as $\lambda_{\psi}=0.1L_{0}$, $B_{0}$ is the amplitude of the magnetic field
(in balance with the total pressure
$P_{0}=\beta_{p}\frac{B_{0}^{2}}{2\mu_{0}}$) equal to $1.42$ G,
$\mathbf{A_{x}}$ is the initial magnetic potential in the x direction, and
$\mathbf{e_{x}}$ is the elementary vector in the x direction.

The pressure of the ionized species $P^{\text{ions}}$ is locally increased in
the current sheet in order to balance the decrease of the Lorentz force in the
current sheet, as follows:

$P^{\text{ions}}(y)=P^{\text{ions}}_{0}+\frac{P_{0}^{\text{ions}}}{\cosh^{2}(\frac{y}{\lambda_{\psi}})},$
(17)

where $P^{\text{ions}}_{0}$ is the total pressure of the ions computed from
the initial thermodynamic conditions ($T_{0}$ and $\rho_{0}$). In order to
trigger the reconnection we have considered a local-centered perturbation on
the magnetic field (while guaranteeing the divergence free of the magnetic
field) as follows:

$\mathbf{B^{\alpha}_{z}}(y,z^{*})=-\alpha(y-\frac{L_{0}}{2})f(y,z^{*})\mathbf{e_{z}},\;(y,z^{*})\in\left[0,L_{0}\right]^{2}$
(18)
$\mathbf{B^{\alpha}_{y}}(y,z^{*})=\alpha(z^{*}-\frac{L_{0}}{2})f(y,z^{*})\mathbf{e_{y}},\;(y,z^{*})\in\left[0,L_{0}\right]^{2}$
(19)

where

$f(y,z^{*})=\exp(-\left[\frac{z^{*}-L_{0}/2}{\lambda_{\psi}}\right]^{2})\exp(-\left[\frac{y-L_{0}/2}{\lambda_{\psi}}\right]^{2}),$
(20)

where $\alpha$ is the amplitude of the perturbation which has been chosen in
such a way that the maximum amplitude of
$|\mathbf{B^{\alpha}}|=\sqrt{|\mathbf{B^{\alpha}_{y}}|^{2}+|\mathbf{B^{\alpha}_{z}}|^{2}}$
corresponds to approximately $10\%$ of the amplitude of the magnetic field of
the Harris current sheet $B_{0}$. This amplitude is chosen in accordance with
Alvarez Laguna et al. (2016); Leake et al. (2012).

For both MFMS and TF, the characteristic initial Alfven speed is calculated as
$V^{\text{Alf}}_{0}=B_{0}/\sqrt{\mu_{0}\rho^{\text{ions}}_{0}}=60$ km s-1,
where $\rho^{\text{ions}}_{0}$ is the total initial mass density of ions. The
initial thermal speed is $V^{\text{therm}}_{0}=\sqrt{8k_{B}T_{0}/\pi
m_{0}}=18.3$ km s-1, where $m_{0}$ is the characteristic mass chosen here as
the mass of protons. No initial velocity field is considered, rather the
velocity field is naturally triggered during the reconnection process. The
numerical simulations based on MFMS and TF have been run until
$t_{\text{end}}=454$ s and $325$ s, respectively. For the rest of this work,
we will focus on a restricted domain in the $z$ direction where $z\in[0,4]$ Mm
and considering the change of variable $z=z^{*}-\delta_{z}$ with
$\delta_{z}=0.1332$ Mm.

###  4.2 Initial composition

| H | H+ | He | He+
---|---|---|---|---
$\mathcal{M}$ | $\mathcal{M}_{\text{TF}}$ | $\mathcal{M}$ | $\mathcal{M}_{\text{TF}}$
$n$ $\left[\%\right]$ | 34.55 | 37.5 | 57.55 | 62.5 | 7.9 | $3.17\times 10^{-6}$
$\rho$ $\left[\%\right]$ | 28 | 37.5 | 46.6 | 62.5 | 25.4 | $10^{-5}$
Table 1: Initial number and mass densities in percentage for MFMS and TF. The
initial number and mass densities percentage are calculated as
$\alpha_{\mathfrak{a}\mathcal{I}}\left[\%\right]=100\times\alpha_{\mathfrak{a}\mathcal{I}}/\sum_{\mathfrak{a}\mathcal{I}\in\beta}\alpha_{\mathfrak{a}\mathcal{I}},\;\alpha\in\\{n,\rho\\}$
and $\beta\in\\{\mathcal{H},\,\mathcal{H}_{\text{TF}}\\}$.

In both TF and MFMS, the initial number and mass densities for all species are
calculated from the total number density $n_{0}$ and temperature $T_{0}$
assuming statistical and thermal equilibrium. The initial total number density
of the plasma is the same for MFMS and TF but the total mass density differs
due to the presence of neutral helium species in MFMS. The initial composition
in number and mass densities for mixtures $\mathcal{M}$ and
$\mathcal{M}_{\text{TF}}$ are given in Table 1. Note that, in MFMS, the
initial abundance of helium species, i.e., with respect to the total hydrogen
number density, is approximately $8.5\%$, i.e., photospheric abundances
(Asplund et al., 2009).

###  4.3 Initial multi-fluid parameters

| H-H+ | H-He | H+-He | He+-H | He+-H+ | He+-He |  $\mathfrak{e}$-H |  $\mathfrak{e}$-He
---|---|---|---|---|---|---|---|---
$\nu^{\text{col}}$ [s-1] | $\approx 10\text{ -- }100$ | $\approx 0.1\text{ -- }1$ | $\approx 0.1\text{ -- }1$ | $\approx 0.1\text{ -- }1$ | $\approx 100\text{ -- }1000$ | $\approx 1\text{ -- }5$ | $100\text{ -- }500$ | $\approx 20\text{ -- }50$
$\lambda_{p}$ [km] | $\approx 0.1\text{ -- }1$ | $\approx 10\text{ -- }100$ | $\approx 10\text{ -- }100$ | $\approx 10\text{ -- }100$ | $\approx 0.01\text{ -- }0.1$ | $\approx 3\text{ -- }7$ | $\approx 2\text{ -- }3$ | $\approx 15\text{ -- }30$
Table 2: Initial range of values (minimum/maximum respectively) of collision
frequencies in s-1 and mean free path in km for MFMS in orders of magnitude.
All the other interactions that are not shown in this table (for example
$\mathfrak{e}$-H+) have a mean free path largely below 0.01 km.

In our initial setup, the rates associated with ionization or recombination
processes have orders of magnitude ranging from $10^{-4}$ to $10^{-1.5}$ s-1
for helium and hydrogen respectively. For MFMS, the initial values of
collisional frequencies and mean free paths (taken as the thermal speed of the
colliding species times collisional frequencies between the two species) of
the different interactions in orders of magnitude are presented in Table 2.
These values highlight the large disparities in the coupling between fluids
due to collisions in MFMS. Note that in the TF model, the collisional rate
between H and H+ is slightly higher than in Table 2 for MFMS due to the
initial mass density of H+, which is higher in TF than in MFMS, as shown in
Table 1.

These preliminary calculations demonstrate that the different fluids involved
in MFMS are expected to be strongly or weakly coupled due to collisions at the
initial stage of the reconnection. For example, interactions involving He
appears to be less coupled due to collisions than interaction H-H+ by two
orders of magnitude. However, interactions associated with He+ and H+ are more
strongly coupled due to collisions than the H-H+ interaction. Note that
timescales associated with ionization and recombination are much slower than
collisional timescales, thus, collisional processes are much faster than
ionization and recombination processes.

One can notice that some of the characteristic mean free paths are not
spatially resolved by the grid size $dy=dz$, in particular those associated
with electrons. Therefore, the coupling due to collisions associated with
electrons appears to be not fully resolved by the current numerical setup and
would require a better spatial resolution. However, the chosen grid size
numerically resolves the mean free paths associated with heavy particles
within the simulation. Therefore, the following simulations resolve the
coupling and decoupling of the various fluids except for electrons.

##  5 Results

In this section, we analyze the results of our two simulations based on MFMS
and TF models, keeping in mind that the mixture associated with MFMS is
$\mathcal{M}=\left\\{\text{H},\text{H}^{+},\text{He},\text{He}^{+},\mathfrak{e}\right\\}$,
and the mixture associated with TF is
$\mathcal{M}_{\text{TF}}=\left\\{\text{H},\text{H}^{+},\mathfrak{e}\right\\}$.
The analysis is performed at different levels. First, we focus on the
structure of the current sheet and its evolution (Section 5.1), and, then, the
velocity fields and highlight the decoupling between fluids (Section 5.2). We
continue with Section 5.3 in which we identify the main heating processes in
both mixtures. Finally, we analyze the evolution of the ionization level in
the current sheet in MFMS and TF and the chemical fractionation between helium
and hydrogen species in MFMS during the reconnection process (Section 5.4).

###  5.1 Current sheet and plasmoids

In both MFMS and TF, our calculations show that two phases of the reconnection
process can be identified. Figure 1 and associated animation illustrates the
distribution of the orthogonal current $\mathbb{J}_{x}$ and total velocity
fields in both $y$ and $z$ directions for both MFMS and TF, at different
times. The total velocity fields are calculated as
$\mathbb{u}_{\text{tot}}=\left(\rho\mathbb{u}\right)_{\text{tot}}/\rho_{\text{tot}}$.
Note that further details about each individual velocity are given in the next
Sections 5.2.

First, the simulation has a laminar or a Sweet-Parker reconnection phase,
where only one elongated current sheet is formed and develops inflow and
outflow velocity fields corresponding to the $y$ and $z$ direction,
respectively (see Parker, 1955, 1963; Sweet, 1958) due to the initial
structure of the magnetic field. This can be clearly identified for both MFMS
and TF in the first column of Fig. 1, as well as in the animation
corresponding to this figure. After this first phase, we have a plasmoid
instability phase where many plasmoids are formed and are advected from the
middle of the reconnection region towards the outflow region along the $z$
axis (the last column of Fig. 1). During this phase, the presence of plasmoids
breaks the current sheet into multiple thinner current sheets, which reconnect
faster than in the Sweet-Parker regime (see Murtas et al., 2021). This is also
in agreement with previous 2D (e.g., Huang & Bhattacharjee, 2010) and 3D
(e.g., Huang & Bhattacharjee, 2016) simulations.

At the very beginning of the laminar phase ($t\approx 54$ s), the magnitude of
the current inside the reconnection region is slightly higher in TF than in
MFMS, as shown in the animation of Fig. 1. Additionally, the magnitude of the
total velocity fields close to the current sheet and the magnitude of the
total velocity in the $z$ and $y$ directions are slightly higher in TF than in
MFMS. The current sheet is thinned sooner in TF than in MFMS, as shown in the
animation of Fig. 1 and in the top left panel of Fig. 2. At the end of the
laminar phase and during the plasmoid instability phase, it remains difficult
to compare the evolution of the current sheet and total velocity fields as it
is unstable and fields change rapidly on short length scales. Therefore,
further analysis has been performed and is presented in Fig. 2.

Figure 1: Distribution of the component $x$ of the current density
$\mathbb{J}_{x}$ in Ampère and total velocity in km s-1 in the $y$ and $z$
directions corresponding to the inflow and outflow directions of the
reconnection region respectively, for the MFMS case (first three rows) and TF
case (last three rows) in a restricted domain defined as
$(y,z)\in[1.9,2.5]\times[0,4]$ Mm. From left to right: evolution of these
distributions at $t=54,\,100$ and $250$ s. The vertical line corresponds to
$z=2.45$ Mm. The associated animation of this figure shows the time evolution.

Figure 2 shows the evolution of the thickness of the current sheet at a fixed
$z=2.45$ Mm (top left), the mean thickness of the whole current sheet (top
right), the length of the current sheet (bottom left), and the maximum total
velocity in the $z$ direction $u_{\text{tot},z}$ (bottom right) taken at the
main reconnection site, for both MFMS and TF. To calculate the length at a
fixed $z$, first, we locate the maximum of the magnitude of current for any
$y\in[2.1,\,2.3]$ Mm. Then, we consider its full-width half maximum (FWHM) in
the $z$ direction to estimate the length. In our simulations, the maximum
current is generally located at $y=L_{0}/2$. Similarly, the $z$ component of
the maximum total velocity at the reconnection site has been taken within the
FWHM along the $z$ axis of the maximum current. Concerning the thickness at a
fixed $z=2.45$ Mm, a similar procedure as the estimation of the length of the
current sheet has been performed. First, we locate the maximum magnitude of
the current for any $y$ in the domain at a fixed $z=2.45$ Mm. We choose the
FWHM as it accurately represents the termination of the reconnection site for
both the laminar and plasmoid instability phases. We extend these calculations
to any $z\in[0.8,3.2]$ Mm in order to calculate the evolution of the mean
thickness and deviation of the current sheet.

Focusing on the left column of Fig. 2, from 0 to approximately $25$ s, for
both MFMS and TF, the thickness at $z=2.45$ Mm and length of the current sheet
decreases drastically. It corresponds to a fast transition from the initial
conditions to the beginning of the laminar phase leading to the formation of
the current sheet. The initial thickness of the current sheet is
$\lambda_{\phi}=0.1L_{0}$ and the length $L_{0}=426$ km. In both MFMS and TF,
from $0$ to $\sim 100$ s, the thickness at $z=2.45$ Mm of the current sheet
decreases to $\sim 25$ km corresponding to the numerical length scale.
Additionally, if we focus on the mean thickness of the current sheet (top
right of Fig. 2), it is clear that the thickness is larger in MFMS than in TF
during the laminar phase. Consequently, TF reduces the thickness of the
current sheet slightly faster than in MFMS, leading to higher velocities
during the laminar regime, as shown in the bottom right of Fig. 2. Note that
there are no large disparities between MFMS and TF in the evolution of the
length of the current sheet.

At $t\approx 100~{}s$, the plasmoid instability phase of the reconnection
process starts. This phase is characterized by strong variation in the
thickness of the current sheet as shown in top left plot of Fig. 2. In MFMS,
we can clearly see that the thickness at $z=2.45$ Mm increases around at
$t=200$ s from 25 to $80$ km. Similarly in TF, a small variation of the
thickness at $z=2.45$ Mm is obtained around $t=140$ s where the thickness
increases from 25 to approximately $37$ km. In general, at any time, the
average size of the plasmoids appears to be slightly larger in MFMS than in TF
meaning that the dilatation of the plasmoids during their formation is more
important in MFMS than in TF, as shown in the top right panel of Fig. 2.

The plasmoid instability phase impacts the length of the current sheet that is
associated with the maximum current $\mathbb{J}_{x}$ in both MFMS and TF
(bottom left plot of Fig. 2). Indeed, after $t\approx 100$ s, in both MFMS and
TF, strong variations of the length of the current sheet are obtained, but on
average decrease. This is because the presence of plasmoids breaks the initial
current sheet formed during the laminar phase into multiple thinner current
sheets. In summary, even during the plasmoid instability phase, on average,
the mean frequency at which the plasmoids are formed and length scales between
them is similar between MFMS and TF.

At any time of the reconnection process, in both MFMS and TF, the maximum
total velocity in the $z$ direction $u_{\text{tot},z}$ is increasing, as
illustrated in the bottom right of Fig. 2. However, it is clear that the
velocity is higher in TF than in MFMS, meaning that stronger flows are
obtained in TF than in MFMS, in both the laminar and plasmoid instability
phase. Note that a similar result has been obtained in the $y$ direction.

Figure 2: Evolution in time of the thickness at a fixed $z=2.45$ Mm (top
left), mean thickness for $y\in[1,\,3.5]$ Mm (top right) and length of the
current sheet (bottom left), and maximum total velocity in the $z$ direction
estimated at the ouflow of the main reconnection site (bottom right), by
taking the full with half maximum (FWMH) of the maximum current density
$\mathbb{J}_{x}$. The dashed and full line corresponds to TF and MFMS
respectively. Note that $z=2.45$ Mm is illustrated in Fig. 1 by the black
vertical line.

In summary, the average thickness of the current sheet and the size of the
plasmoids appears to be slightly larger in MFMS than in TF, and, reconnection
with the TF model develop stronger flows than the MFMS model. In this context,
the presence of helium species appears to only slightly impact the evolution
of the current sheet as well as the flows. The length of the current sheet and
the frequency at which the plasmoids are formed are similar between the two
models. In the next sections, we will see that, even though these properties
only slightly differ between the two models, the decoupling of the particles
and composition of the plasma in the current sheet show large disparities,
leading to different evolution of the heating processes and temperatures.

###  5.2 Decoupling of the particles and velocity fields

In this section, the decoupling of the particles and velocity fields are
investigated. To do so, we focus on the $z$ and $y$ components of the velocity
for each fluid separately. These correspond to the outflow and inflow
directions of the magnetic reconnection, respectively.

The first four rows (starting from the top) of Figures 3 (Fig. 4) represent
the distribution of the $y$ ($z$) component of the velocity of protons H+ and
the drift velocities with all the other particles in MFMS. The last two rows
of the figures correspond to the velocity of protons H+ and drift velocity
between H and H+ in TF. From left to right columns, it shows the time
evolution at times $t=89$ s, $t=282$ s and $t=327.5$ s, respectively.

Figure 3: Distribution of the $y$ component of the velocity of H+ and drift
velocities with other species in km s-1 at different times for MFMS and TF.
From first to fourth row: distribution of the velocity $u_{y,\text{H}^{+}}$,
the drift velocities $u_{y,\text{H}^{+}}-u_{y,\text{H}}$,
$u_{y,\text{H}^{+}}-u_{y,\text{He}^{+}}$ and
$u_{y,\text{H}^{+}}-u_{y,\text{He}}$ in MFMS respectively. The last two rows
corresponds to the distribution of the velocity $u_{y,\text{H}^{+}}$ and the
drift velocity $u_{y,\text{H}^{+}}-u_{y,\text{H}}$ in TF. From left to right
column: different times $t=136$, $t=278.5$ and $327.5$ s are presented. An
animation of this figure is available.

Figure 4: Similar layout as Fig. 3 for the $z$ component of the velocities.
Unlike Fig. 3, the drift velocity between H+ and He+ is not shown since its
magnitude is negligible compared to the other quantities. An animation of this
figure is available.

During the laminar phase, the velocities of H+ increases up to a maximum
achieved at approximately $1$ km s-1 and 20 km s-1 $y$ and $z$ components in
both simulations. Note that TF seems to show slightly higher velocities, as
described in Section 5.1. This maximum is reached at the end of the laminar
phase which corresponds to $t\approx 100$ s for both MFMS and TF. This
corresponds to the time where the maximum length of the current sheet is
achieved, as shown in Fig. 2. In the first column of Figs. 3, and 4, the drift
velocities between particles are small compared to $u_{y,\text{H}^{+}}$ and
$u_{z,\text{H}^{+}}$ which implies that there is a large collisional coupling
between species during the laminar phase. The decoupling between H+ and He is
highest, as shown in panel D of Fig. 3 and F of Fig. 4. Neutral helium species
moves slightly slower than H+ by approximately $0.25$ and 0.15 km s-1 in the y
and z directions respectively. Concerning TF, no significant decoupling
between H and H+ is seen.

In order to understand the evolution of the velocity fields, we have also
calculated and analyzed the distributions of all the forces acting on each
fluid involved in MFMS or TF, at different times. In the following, for the
sake of simplicity, we briefly summarize them.

During the laminar phase, the initial structure of the magnetic field allows
the development of a Lorentz force and leads to the formation of outflow and
inflow velocities. The Lorentz force and gradient of pressure associated with
H+ are the dominant forces of the plasma in both MFMS and TF. Therefore,
during the laminar phase, the dynamics of the reconnection process itself in
the current sheet are dominated by these forces. This is also consistent with
the fact that H+ is the most dominant species in number and mass density in
both MFMS and TF at that time. In MFMS, no drift velocities associated with
H+-H and H+-He+ occur since the species H+, H and He+ are highly coupled due
to collisions. In other words, in each individual momentum equation associated
with H and He+, the magnitude of the momentum exchange terms due to collisions
with H+ is large enough to highly couple the momentum equations associated
with H, H+ and He+ altogether. Similarly in TF, the fluids H and H+ are highly
coupled due to collisions during the laminar phase.

In contrast, in MFMS, the neutral helium species is decoupled from the other
fluids. This drift at the current sheet is due to the fact that neutral He
decouple from H+ when the thickness of the current sheet thins down to the
order of the mean free path between He and H+ during the laminar phase, and,
neutral He are not dragged by collisions into the current sheet as quickly as
H+ or H. The decoupling between He and H+ due to the lack of collisions is
clearly illustrated in panel D of Fig. 3 and in panel C of Fig. 4. The fact
that the decoupling between He and H+ is stronger than other interactions is
attributed to the collisional frequencies that involves neutral helium species
which are much smaller by more than two orders of magnitude than all the other
collisional frequencies, as shown in Table 2 (for the initial conditions).
These differences are attributed to the magnitude of the collision integrals
at the kinetic level that leads to smaller collisional frequencies, as
described with further details in Wargnier et al. (2022). This effect leads to
a slower velocities for neutral helium species with respect to H+, as
illustrated by panel D of Fig. 3 and panel C of Fig. 4.

In the second and last column of Figs. 3 and 4 corresponding to the plasmoid
instability phase, we note that drift velocities are increasing during the
formation of plasmoids for both MFMS and TF. This is also illustrated in the
animations. When plasmoids are produced, local velocity fields in opposite
directions are generated due to a local reconnection of the magnetic field
lines. It should be noted that the structure and evolution of the velocity of
H+ is similar for both MFMS and TF, however, the drift velocities differ
between the two. In the case of MFMS, the drift between H and H+ is stronger
at the boundary of plasmoids and inside the current sheet, reaching
approximately 1 and 0.1 km s-1 for the y and z components respectively, as we
can see in panels H and N of Fig. 3 and panels G and L of Fig. 4. In TF, the
drift velocity between H and H+ is smaller than in MFMS by a factor of three
and is showing some oscillating pattern structure parallel to the current
sheet, as shown in panel R of Fig. 3. Similar to the decoupling between He and
H+, the drift between H and H+ develops during the plasmoid instability phase,
as the current sheet thins down to the order of the mean free path (a few
kilometers corresponding to a few grid points) between H and H+ and neutral H
does not collide enough with H+ to be dragged into the current sheet as
quickly as H+. This decoupling is consistent with the results obtained by
Leake et al. (2012); Leake & Linton (2013).

The drift between H and H+ is higher in MFMS than in TF because the
collisional frequency associated with H-H+ interactions is smaller in MFMS
than in TF. As shown in Wargnier et al. (2022), the collisional frequencies
are calculated from the products of 1- the collision integral associated with
H-H+ interaction, 2- H and H+ number densities and 3- thermal speed associated
with H-H+ interaction. In our chosen conditions, the number densities of H and
H+ in MFMS are smaller than in TF due to the presence of helium species.
Additionally, the temperature inside the current sheet is higher in MFMS than
in TF, as shown later in Section 5.3, leading to a decrease of the collision
integral associated with H-H+ interactions, as shown in Wargnier et al.
(2022). These two effects contribute to decrease the collisional frequency
associated with H-H+ interactions in MFMS compared with TF. Note that if the
temperature increases, the thermal speed contributes to increase the
collisional frequency, as shown in the definition from Wargnier et al. (2022).
However, in the chosen conditions and for H-H+ interactions, this contribution
is smaller than the contribution of the H-H+ collision integral and number
densities.

In panels L and R of Figure 3 and panels J and O of Figure 4, the oscillating
pattern of the drift velocity between H and H+ in TF is attributed to the
spatial distribution momentum exchange term associated with H-H+ interactions.
These collisions compensate for the localized unbalance of the forces i.e.,
the Lorentz force and pressure gradients, when the current sheet thins down to
the mean free path between H and H+. A similar effect has been obtained in the
distribution of drift velocity between H+ and He+ in panel O of Fig. 3 where
the characteristic wavelength and magnitude of this pattern are much smaller
than for the drift between H and H+. However, this effect is attributed to two
phenomena: the momentum exchange term associated with H+-He+ interactions in
MFMS tries to compensate for the localized unbalance of forces (pressure
gradients of H, H+, He and He+ and Lorentz forces acting on H+ and He+
separately) at much smaller spatio-temporal scales than in TF. In addition,
the ions experience a synchronized cyclotron-type motion along the magnetic
field lines, also known as the ion-coupling effect, leading to a sinusoidal
ion drift velocity distribution perpendicular to the magnetic field line, as
described with further details in Martínez-Sykora et al. (2020). Note that the
$z$ component of the drift velocity associated with H+-He+ interactions is
negligible. Consequently, the distribution of the latter has not been shown
here.

As in the laminar phase, the decoupling between He and H+ remains important
due to the lack of collisions between the two species. This drift velocity is
high both inside the current sheet and plasmoids. The maximum drift velocity
between He and H+ reaches more than 1 km s-1 and 0.3 km s-1 in the $y$ and $z$
directions respectively. We can also notice some internal structure of the
drift velocity inside the plasmoids (panels J and P of Fig. 3 and panels H and
M of Fig. 4).

###  5.3 Multi-fluid heating processes

####  5.3.1 Heating processes due to departure in velocities

In order to analyze the heating processes and identify the type of collisions
mainly responsible for the increase in temperature, we focus on the dominant
heating terms due to drift velocities
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\textbf{u},\text{col}}$
in both models MFMS and TF. As explained in Section 2.3, the heating terms
associated with drift velocities are not adiabatic terms and can increase the
total thermal energy of the plasma. Therefore, the heating processes for both
MFMS and TF could be strongly connected to the structure of the drift
velocities which have been investigated in previous Sections 5.2.

Under the chosen thermodynamic conditions presented in Section 4, the electron
thermal conduction and the terms associated with ionization or recombination
processes $Q_{\mathfrak{a}\mathcal{I}}^{\text{ion}}$ or
$Q_{\mathfrak{a}\mathcal{I}}^{\text{rec}}$ have negligible impact on the
heating processes, in both MFMS and TF. Indeed, as described in Section 2, the
heating terms associated with ionization or recombination processes for any
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$ do not contribute to the total thermal
energy of the plasma (apart from radiative losses term). Thus, they do not
contribute to the current sheet’s temperature increase. In all thermal energy
equations, the magnitude of these terms is negligible compared to the
thermalization terms due to collisions since ionization and recombination
timescales are bigger than collisional timescales, as described in Section 4.2
and also in Section 5.4.1.

Figure 5 represents the distribution of the rate of the heating terms
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\textbf{u},\text{col}}$
in a 1D cut of the current sheet located at $z=2.45$ Mm for all the possible
interactions involved in both MFMS (top row) and TF (middle row) at different
times $t=11.5$ s and $t=28.5$ s from left to right column. Note that these
times correspond to a different stage of the laminar phase. The rates are
calculated by dividing the heating terms, produced by drift velocities, by the
total thermal energy
$e^{\alpha}_{\text{tot}}=\sum_{\mathfrak{a}\mathcal{I}\in\alpha}e_{\mathfrak{a}\mathcal{I}},\,\alpha\in\\{\mathcal{M},\mathcal{M}_{\text{TF}}\\}$.
We also show the distribution of the temperatures of each species in MFMS and
TF (bottom row) at $z=2.45$ Mm. Note that an animation of this figure is
available.

At the beginning of the laminar phase in MFMS, it is clear that the heating
processes produced by the drifts between H+ and He, H and He dominate in the
current sheet, as observed in panel A of Fig. 5. This effect is consistent
with the previous Sections 5.2 showing that the drift velocities between
neutral helium and hydrogen species are much higher than any other
interactions during the laminar phase. Additionally, the heating due to
collisions between electrons and all heavy species (ohmic collisions) is at
the same order of magnitude as H+-He and H-He interactions but localized
around the current sheet. The heating due to the drift between H+ and H is not
negligible in the current sheet. Still, it is much lower than the other
interactions (panel A of Fig. 5), which is also consistent with the magnitude
of the drift velocity in Sections 5.2. In the case of TF, the heating produced
by the drift between electrons and heavy species (H and H+) is dominating, as
shown in panel B. Additionally, the latter is slightly higher in magnitude as
the heating is produced by the drift between electrons and heavy species in
panel A (black dashed line in panel A). However, if we sum all the heating
terms of panel A and compare it with the sum of all the heating terms in panel
B, it is clear that the total heating produced by velocity drifts is higher,
by approximately a factor two, in MFMS than TF due to the heating terms
$Q_{\text{H},\text{He}}^{\textbf{u},\text{col}}+Q_{\text{He},\text{H}}^{\textbf{u},\text{col}}$
and
$Q_{\text{H}^{+},\text{He}}^{\textbf{u},\text{col}}+Q_{\text{He},\text{H}^{+}}^{\textbf{u},\text{col}}$.

At a later time of the laminar phase, corresponding to the second column of
Fig. 5, the temperatures of all species in MFMS increase and are higher than
in TF. The temperature achieves its maximum at the center of the current sheet
at $27\,000$ K in MFMS. In contrast, in TF, the temperature is approximately
$19\,000$ K at the center of the current sheet (panel F). By comparing panels
A and D, the magnitude of the heating terms due to collisions associated with
interactions H-He, H+-He, and electrons-heavy species increase in time inside
the current sheet. This increase is attributed to drift velocities increasing
with time during the laminar phase. Similarly, if we sum all the heating terms
in panels A and B, the sum of all the heating terms is higher in MFMS than in
TF. Consequently, temperatures inside the current sheet are higher in MFMS
than in TF. Note that temperatures of all species are equal during the laminar
phase in both MFMS and TF.

Figure 5: Distribution of the rates of the heating terms
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\mathbb{u},\text{col}}+Q_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime},\mathfrak{a}\mathcal{I}}^{\mathbb{u},\text{col}},\,(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
or $\mathcal{M}^{2}_{\text{TF}}$ associated with drift velocities divided by
the total thermal energy
$e_{\text{tot}}=\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{M}}e_{\mathfrak{a}\mathcal{I}}$
in s-1 (top row for MFMS and middle row for TF) and temperatures of each
particles (bottom row) along the $z$ axis in the current sheet at position
$y=L_{0}/2$ in both MFMS and TF, respectively. From left to right column:
distributions at different times $t=11.5$ s and $t=28.5$ s corresponding to
the laminar phase. Note that the vertical axis ranges from $0$ to $1.4\times
10^{-3}$ s-1. An animation of this figure is available.

Figure 6: From top to bottom, for any $(y,z)\in[1.9,2.5]\times[0,4]$ Mm:
distribution of the decimal logarithm of the rates of the heating terms
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\mathbb{u},\text{col},S},\,(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
divided by the total thermal energy
$e_{\text{tot}}=\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{M}}e_{\mathfrak{a}\mathcal{I}}$
in s-1 in MFMS, the temperature of H+ in MFMS, the decimal logarithm of the
heating terms
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\text{col}},\,(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}_{\text{TF}}^{2}$
divided by the total thermal energy
$e_{\text{tot}}=\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{M}_{\text{TF}}}e_{\mathfrak{a}\mathcal{I}}$
in s-1 in TF and the temperature of H+ in TF. The other heating terms have not
been shown here since their magnitude are negligible compared to the other
terms. From left to right: distributions at different times $t=210$, $250$ and
$282.5$ s corresponding to different stages of the plasmoid instability phase.
An animation of this figure is available.

Figure 6 shows the heating rates due to the drift velocities
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\mathbb{u},\text{col}}+Q_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime},\mathfrak{a}\mathcal{I}}^{\mathbb{u},\text{col}},\,(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
or $\mathcal{M}_{\text{TF}}^{2}$ divided by the total thermal energy and the
temperature of protons $T_{H^{+}}$, for MFMS (first five rows) and TF (last
two rows), at different times of the simulation corresponding to the plasmoid
instability phase. From left to right columns, it corresponds to times
$t=210$, $250$ and $282.5$ s respectively. Higher temperatures are reached
inside the current sheet and plasmoids than during the laminar phase, as shown
in panels E, L, S, G, N, and U from Fig. 6 by comparing with temperatures
presented in Fig. 5.

Additionally, higher temperatures are reached in MFMS than in TF, as
illustrated in panels E, L, and S in comparison with panels G, N, and U of
Fig. 6. Indeed, in TF, the temperatures achieve a maximum of $35\,000$ K,
whereas in MFMS the temperatures can reach up to approximately $50\,000$ K. In
TF, this temperature is reached with the H and H+ drift heating (panels F, M,
and T). These panels show that the rate of the heating terms ranges from
$10^{-2.8}$ to $10^{-3}$ s-1 inside the current sheet and plasmoids,
increasing the temperature $T_{H^{+}}$ from $22\,000$ K (panel G) to $30\,000$
K. In the plasmoid located at $z\approx 1.5$ Mm, the temperature achieves a
maximum of $35\,000$ K due to the H-H+ drift (panel U).

In MFMS, since many species and interactions are involved, the heating
produced by drift velocities is more complex to interpret than TF and requires
the analysis of several components of the heating terms during the plasmoid
instability phase.

First, if we compare panels A, H, and O with F, M, and T, the rate associated
with the heating terms due to the drift between H and H+ is much higher in
MFMS than in TF. The disparity can reach up to one order of magnitude in some
regions of the current sheet. For example, at $z\approx 1.6$ or $z\approx 3.2$
Mm of panel O, the rate achieves a maximum at $10^{-2}$ s-1. Therefore, the
temperature $T_{H^{+}}$ increases in these regions from $30\,000$ K to $\sim
50000$ K between $t=210$ and $t=282.5$ s, (panels E, L, and S). Even though
most of the heating comes from the H and H+ collision, the heating due to the
H+ and He+ collision is not negligible at $z\approx 1.6$ or $z\approx 3.2$ Mm
(panel R). This last heating rate is approximately $10^{-3.1}$ s-1(panels D, K
and R).

Focusing on the formation of the plasmoid located at $z\approx 1.9$ Mm in
panels E, L, and S of Fig. 6, one can notice that the thickness of this
plasmoid increases. This dilation is the consequence of the balance of the
thermal pressure of the plasma inside the current sheet that increases with
respect to the magnetic pressure, which is constant. Therefore, during this
dilation, the temperature inside the plasmoid is constant, and $\sim 50\,000$
K, while its thickness increases from $\sim 0.08$ Mm to $\sim 0.2$ Mm.
Therefore, the dilation process of the thermodynamic system defined as the
volume of plasma located inside the plasmoid at $z\approx 1.9$ Mm between
$t=210$ and $t=282.5$ s is not adiabatic. However, if we look at panels A, H,
and O or D, K, and R, it is clear that the heating terms due to the drift
between H and H+ or H+ and He+ are negligible inside the plasmoid. Therefore,
one has to focus on other heating sources inside the plasmoids to explain the
non-adiabaticity of the dilation.

Indeed, the heating due to the friction between He and He+ is not negligible
inside the plasmoid, as shown in panels C, J, and Q of Fig. 6. In panel J, at
$z\approx 1.9$ Mm, the heating rate associated with He-He+ collisions is
approximately equal to $10^{-3.4}$ s-1 which allows maintaining the plasma
temperature inside the plasmoid to $50\,000$ K during the dilatation. The
cooling rate associated with the adiabatic expansion of the plasmoids is
negligible compared to this heating rate due to He-He+ collisions ($\approx
10^{-4.2}$ s-1). The temperature of the second plasmoid located at $z\approx
3$ Mm is slightly lower than the other plasmoid and ranges from $35\,000$ K to
$40\,000$ K. This is because the initial plasma temperature in the current
sheet where this plasmoid is initially formed is lower than the initial
plasmoid temperature located at $z\approx 1.9$ Mm. Indeed, during the plasmoid
instability phase, the temperature distribution is not uniform along the $z$
axis leading to some variations in the temperature of the formed plasmoids.

Furthermore, the heating rate due to He and He+ collisions in the plasmoid at
$z\approx 3$ Mm is lower than at $z\approx 1.9$ Mm, and approximately equal to
$10^{-4}$ s-1, as we can see in the panel J of Fig. 6. The heating due to the
friction between He and H+ slightly contributes to the heating of the
plasmoids, as we can see in panels B, I, and P. Note that this contribution is
not included in the TF model. Thus, during the expansion of the plasmoids, the
temperature inside the plasmoids decreases and diffuses with time.

In summary, the temperature is higher in MFMS than in TF both during the
laminar (with temperatures of $30\,000$ K for MFMS and $20\,000$ K for TF) and
plasmoid instability phases ($65\,000$ K for MFMS and $35\,000$ K for TF).
Consequently, in MFMS and for the chosen initial thermodynamic conditions, the
plasma reaches a sufficient temperature inside the current sheet to eventually
ionize helium species rapidly. During the plasmoid instability phase, the main
heating process is due to the collisions between H and H+. It allows
increasing the temperature inside the current sheet from $30\,000$ K to
$50\,000$ K. Then, plasmoids are produced and expanded following a non-
adiabatic process. In this context, neutral helium is rapidly ionized (as
described with further details in Section 5.4) and collides with all the other
particles, in particular He, inside the plasmoids. It allows the maintenance
of the temperature of the plasmoids at $50\,000$ K (for the plasmoid located
at $z=1.8$ Mm). A maximum temperature of $65\,000$ K is achieved inside the
plasmoids at a later time of the plasmoid instability phase.

####  5.3.2 Decoupling of the temperatures in MFMS

In this section, we focus only on MFMS at a later time of the plasmoid
instability phase. Indeed, these effects appear to be out of reach in the TF
model. In MFMS, our previous results from Section 5.3.1 have shown that most
of the temperatures are in equilibrium during the laminar and plasmoid
instability phase. This is because the thermalization terms for any possible
interactions
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{T,\text{col}},\,(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
have larger magnitude than any other heating terms such as
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\textbf{u},\text{col}},\,(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
or $Q^{\text{spitz}}_{\mathfrak{e}}$. In other words, the characteristic
timescale of the thermalization terms range from 10 s to just a few seconds
whereas the heating terms associated with drift velocities range from 100 s to
higher, as shown in Fig. 6.

Departures in temperature occur for a couple of species
$(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$
if the rate associated with all the heating processes of particles
$\mathfrak{a}\mathcal{I}$ becomes much higher than the rate of the
thermalization term
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{T,\text{col}}/e_{\mathfrak{a}\mathcal{I}}$.
In other words, for any
$(\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime})\in\mathcal{M}^{2}$,
neglecting the electron thermal conduction and radiative losses, this
condition can be written as

$\frac{Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{T,\text{col}}}{e_{\mathfrak{a}\mathcal{I}}}\ll\frac{1}{e_{\mathfrak{a}\mathcal{I}}}\left(\sum_{\mathfrak{a}^{*}\mathcal{I}^{*}\in\mathcal{M}}Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{*}\mathcal{I}^{*}}^{\textbf{u},\text{col}}\right)=\frac{Q_{\mathfrak{a}\mathcal{I}}^{\textbf{u},\text{col}}}{e_{\mathfrak{a}\mathcal{I}}},$
(21)

Figure 7: From top to bottom row: distribution along the $z$ axis at
$y=L_{0}/2$ of the rates of the total heating term
$Q_{\text{He}}^{\textbf{u},\text{col}}/e_{\text{He}}$ and all the
thermalization terms
$Q_{\text{He},\mathfrak{a}\mathcal{I}}^{T,\text{col}}/e_{\text{He}}$, with
$\mathfrak{a}\mathcal{I}\in\mathcal{M}$, distribution of all the temperatures
and distribution of temperature differences $\Delta T$ with respect to
$T_{\text{He}}$ in the current sheet at position $y=L_{0}/2$ for MFMS. From
left to right column: distributions at different times $t=307.5$ s and
$t=314.5$ s.

In order to highlight the thermal decoupling effects occurring during the
reconnection process, we focus only on neutral helium. Neutral helium shows
strong thermal decoupling effects due to the high decoupling in collisions
with other species, as shown in the previous sections. Figure 7 shows the rate
$Q_{\text{He}}^{\textbf{u},\text{col}}/e_{\text{He}}$ distribution compared
with all the rates associated with the thermalization terms
$Q_{\text{He},\mathfrak{a}\mathcal{I}}^{T,\text{col}}/e_{\text{He}},\,\mathfrak{a}\mathcal{I}\in\mathcal{M}$
at times $t=307.5$ s, and $t=314.5$ s in the current sheet at $y=L_{0}/2$. On
top of that, we also show the distributions of all the temperatures (panels B
and E) as well as the differences in temperature $\Delta
T=T_{\text{He}}-T_{\mathfrak{a}\mathcal{I}},\,\mathfrak{a}\mathcal{I}\in\mathcal{H}$
with respect to the neutral helium temperature $T_{\text{He}}$ (panels C and
F).

First, if we focus on the plasmoid located at $z\approx 2$ Mm at $t=307.5$ s,
in panel A, the rate of the total heating term
$Q_{\text{He}}^{\textbf{u},\text{col}}/e_{\text{He}}$ is slightly higher than
the rates $Q_{\text{He},\text{H}^{+}}^{T,\text{col}}/e_{\text{He}}$ and
comparable to $Q_{\text{He},\text{He}^{+}}^{T,\text{col}}/e_{\text{He}}$.
After a few seconds, corresponding to the panel D of Fig. 7, the rate of the
total heating term $Q_{\text{He}}^{\textbf{u},\text{col}}/e_{\text{He}}$
increases drastically (500% in 7 seconds) and becomes much larger than the
rates of thermalization terms
$Q_{\text{He},\text{H}^{+}}^{T,\text{col}}/e_{\text{He}}$ and
$Q_{\text{He},\text{H}}^{T,\text{col}}/e_{\text{He}}$ and stay comparable to
$Q_{\text{He},\text{He}^{+}}^{T,\text{col}}$. In this framework, the condition
in Eq. 21 is satisfied for the couples of species $(\text{He},\text{H}^{+})$
and $(\text{He},\text{H})$. This leads to a fast and strong thermal decoupling
where the departures in temperatures $T_{\text{He}}-T_{\text{H}^{+}}$ and
$T_{\text{He}}-T_{\text{H}}$ achieve a maximum at $\approx 3000$ K, occurring
in just a few seconds, at the center of the plasmoid.

During the reconnection process, thermal decoupling effects can be stronger
leading to a neutral helium temperature that can become much larger and reach
up to $90\,000$ K where other species temperature stay at $70\,000$ K. Note
that this effect occurs at later time of the plasmoid instability phase.
Additionally, a small temperature decoupling is seen in neutral hydrogen, as
illustrated with the blue line in panels C and F. However, this effect remains
smaller than temperature decoupling effects associated with neutral helium.

###  5.4 Evolution of the composition

####  5.4.1 Evolution of the ionization levels in MFMS and TF

Figures 9 (first two rows) and 8 show the distribution of the total ionization
fraction of MFMS and TF. Figure 8 is a cut in the middle of the current sheet
at $y=L_{0}/2$ during the laminar phase (on the left) and plasmoid instability
phase (on the right). In Fig. 8, the distributions of the ionization fraction
of hydrogen and helium species are shown. Figure 10 has the same layout as
Fig. 9 but for the distribution of the ionization and recombination
frequencies $\nu^{\text{ion}}_{\text{H},\text{H}^{+}}$,
$\nu^{\text{rec}}_{\text{H}^{+},\text{H}}$ for MFMS and TF, and
$\nu^{\text{ion}}_{\text{He},\text{He}^{+}}$ for MFMS. These frequencies
allows the assesment of the evolution of ionization and recombination effects
and the composition of the plasma is altered during the reconnection process.
The recombination frequency $\nu^{\text{rec}}_{\text{He}^{+},\text{He}}$ is
negligible and much smaller than the other ionization and recombination
frequencies, thus, it is not shown here.

Figure 8: Distribution of the total ionization fraction for both MFMS and TF,
ionization of hydrogen and helium species in MFMS in the middle of the current
sheet at position $y=L_{0}/2$ at times $t=113.5$ s (left) and $t=282$ s
(right) respectively. The full black, dashed black and dot black lines
corresponds to the total ionization fraction, the ionization of hydrogen
species $x_{\text{H}^{+}}=n_{\text{H}^{+}}/(n_{\text{H}^{+}}+n_{\text{H}})$
and the ionization of helium species
$x_{\text{He}^{+}}=n_{\text{He}^{+}}/(n_{\text{He}^{+}}+n_{\text{He}})$ for
MFMS. The full red line correspond to the total ionization fraction in TF. An
animation of this figure is available.

Figure 9: Distribution of the total ionization fraction in TF (first row) and
MFMS (second row), molar fraction of helium species (third row)
$(n_{\text{He}}+n_{\text{He}^{+}})/(\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{M}}n_{\mathfrak{a}\mathcal{I}})$
in MFMS, at different times $t=113.5$ s, $t=282$ s and $327.5$ s for any
$(y,z)\in[1.8,2.6]\times[0,L_{0}]$ Mm (from left to right column)
respectively. The vertical line corresponds to $z=0.825$ Mm. An animation of
this figure is available.

Figure 10: Distribution of the decimal logarithm of the ionization and
recombination frequencies in s-1 of MFMS (first three rows) and TF (last two
rows) at different times $t=113.5$ s, $t=282$ s and $327.5$ s for any
$(y,z)\in[1.8,2.6]\times[0.4,4]$ Mm (from left to right column) respectively.
From top to bottom row: distribution of
$\nu^{\text{ion}}_{\text{H},\text{H}^{+}}$,
$\nu^{\text{rec}}_{\text{H}^{+},\text{H}}$ and
$\nu^{\text{ion}}_{\text{He},\text{He}^{+}}$ in MFMS, and distribution of
$\nu^{\text{ion}}_{\text{H},\text{H}^{+}}$,
$\nu^{\text{rec}}_{\text{H}^{+},\text{H}}$ in TF. Note that the distribution
of $\nu^{\text{rec}}_{\text{He}^{+},\text{He}}$ is not shown since the
magnitude is negligible and much smaller than $10^{-4}$ s-1. An animation of
this figure is available.

During the laminar phase, represented in the first column of Figures 9 (panels
A and B) and 10 (panels A to E), and left plot of Fig. 8, the ionization level
is increasing in the current sheet for both MFMS and TF due to the ionization
of hydrogen species. Indeed, the ionization frequency
$\nu^{\text{ion}}_{\text{H},\text{H}^{+}}$ in the current sheet, in both MFMS
and TF, is higher by more than one order of magnitude than the recombination
frequency $\nu^{\text{rec}}_{\text{H}^{+},\text{H}}$, as shown in Fig. 10.
Furthermore, the recombination frequency
$\nu^{\text{rec}}_{\text{H}^{+},\text{H}}$ in MFMS (panel B) is lower than in
TF (panel E). This leads to an ionization level that is much higher in MFMS
(ranging from $85\%$ to $95\%$) than in TF (ranging from $70\%$ to $80\%$)
during the laminar phase, as shown on the left of Fig. 8. These disparities
obtained in the ionization and recombination frequencies between MFMS and TF
are attributed to the difference in temperature of the plasma in the current
sheet during the laminar phase. Indeed, as described in the previous Section
5.3, the presence of helium species that collide with hydrogen species in MFMS
leads to higher temperature in the current sheet than in TF. This increase of
temperature gives rise to a larger hydrogen ionization in MFMS than in TF
during the laminar phase.

The ionization levels in plasmoids are drastically increasing for both MFMS
(right panel in Fig. 8) and TF compared to the laminar phase (left panel). In
TF, we can distinguish three plasmoids located at $y=1.5$ Mm, 3.1 Mm and 3.7
Mm, and, in MFMS, two plasmoids located at $y=1.9$ Mm and 3 Mm. In TF, at the
middle of the plasmoid located at $y=1.5$ Mm, the ionization level increases
up to more than 98%. This is attributed to a local increase of the temperature
inside plasmoids, giving rise to an increase of H+ due to ionization process,
as illustrated in panel I of Fig. 10.

In MFMS, the friction of hydrogen species with neutral helium allows a higher
temperature level in the plasmoids than in TF during the plasmoid instability
phase, as described in the previous Section 5.3. This level of temperature is
high enough to ionize helium, as illustrated in Fig. 8 characterized by a
large increase of the black dot line. This is also illustrated in panels H and
M of Fig. 10 where we can clearly see that
$\nu^{\text{ion}}_{\text{He}^{+},\text{He}}$ increases drastically, especially
in the core of plasmoids. At this moment, a self-sustaining mechanism occurs.
Indeed, when He+ species are produced, they collide with both H+ and He, as
described in Section 5.3. In this context, an additional heating process
occurs due to collisions involving He+. This effect will slightly increase the
temperature at the center of the plasmoid, leading to even more He+ due to the
ionization process. This is characterized by the increase of
$\nu^{\text{ion}}_{\text{He}^{+},\text{He}}$ in the core of plasmoids, as
illustrated in panels H and M of Fig. 10.

####  5.4.2 Chemical fractionation in MFMS

In this section, we focus only on MFMS and chemical fractionation between
helium and hydrogen species during the reconnection process. In Fig. 9 (last
row), the distributions of the helium species molar fraction at different
times, corresponding to the laminar and plasmoid instability phase, is shown.
The molar fraction is defined as
$(n_{\text{He}}+n_{\text{He}^{+}})/(\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{M}}n_{\mathfrak{a}\mathcal{I}})$.

During the laminar phase, in MFMS, the helium species molar fraction (panel C
of Fig. 9) decreases in the current sheet (from 7.9% to 6%). Additionally, we
also notice a small increase of the helium species molar fraction at the
boundary of the current sheet close to the outflow regions. Indeed, since
helium is not ionized, these effects are attributed to the accumulation of
neutral helium species which are decoupled from hydrogen species nearby the
current sheet, as described in Section 5.2.

During the plasmoid instability phase, the helium species molar fraction
increases (from 6% to 10%) in the core of plasmoids, as shown in panel F of
Fig. 9). This separation process, also called chemical fractionation, occurs
inside the plasmoids and remains even when the latter are advected along the
$z$ axis. This separation process occurs due to the ionization of helium
species during the instability phase.

At a later time of the plasmoid instability phase, the chemical fractionation
process gets stronger, leading to more helium species inside the current sheet
that are advected towards outflow regions. This effect is illustrated in panel
I of Fig. 9. Indeed, an increase of helium species, up to 12%, in the middle
of the current sheet at $y=L_{0}/2$ and $z\in[0.4,1.5]$ Mm is shown. Note that
the fusion of plasmoids also leads to a large increase of helium abundance as
we can see in the same figure, in the current sheet at $z=2.9$ Mm.

To characterize the chemical fractionation in the outflows, the fluxes in
number densities of helium and hydrogen species have been calculated and
compared at the exhaust of the reconnection event until $t_{\text{end}}=454$
s. These calculations have been performed at a fixed $z=0.825~{}$ Mm. For any
species $\mathfrak{a}\mathcal{I}\in\mathcal{M}$, its corresponding flux in
number density is calculated as

$\phi\left(\mathfrak{a}\mathcal{I}\right)=dt\,dx\int_{l_{y}}\left[u_{z,\mathfrak{a}\mathcal{I}}n_{\mathfrak{a}\mathcal{I}}\right]\left(y,\,z=0.825~{}Mm\right)dy,$
(22)

where $dx$ is the grid size in the $x$ direction and set to one, $l_{y}$ is
the length in the $y$ direction limited to the outflow region. We define the
outflow regions as the FWHM of the maximum outflow at $z=0.825~{}$ Mm.
Similarly as the thickness of the current sheet in Section 5.1, this length
has been calculated by locating the maximum outflow velocity at $z=0.825~{}$
Mm. Then, we consider the termination of the exhaust at approximately $90\%$
of this maximum in the $y$ direction to estimate $l_{y}$.

Figure 11 shows the time evolution of the maximum and mean total outflow
velocity (left) and the ratio of flux of helium species with respect to the
flux of all species
$\phi(\text{He}+\text{He}^{+})/(\phi(\text{He}+\text{He}^{+})+\phi(\text{H}+\text{H}^{+}))$
(right). In Fig. 11, if the ratio is above 0.079 (corresponding to the initial
total number density of helium species 7.9% shown in Table 1) means helium
enrichment and, otherwise, hydrogen enrichment.

Our calculations show that during most of the reconnection process, the flux
in the number density of helium species represents from $6$ to $8\%$ of the
total flux of all species, as shown in the bottom left panel of Fig. 11.
Helium is enhanced at the exhaust when the helium ionization starts at
$t\approx 250$ s during the plasmoid instability regime. It increases the
helium flux up to 14% of the total flux, as shown in the bottom left panel of
Fig. 11. These calculations are consistent with our figures presented in Fig.
9. The helium enrichment in the outflows during the strongest flows and late-
formed plasmoids also depends on the amount of helium species collected during
the formation of plasmoids and the current sheet. In some localized events on
the right-hand side of the reconnection, we found the opposite behavior when
He has been previously deployed, and hydrogen is left in the latest plasmoids
leading to the opposite, i.e., enrichment of hydrogen.

Note that, from $t\approx 250$ s to $t_{\text{end}}=454$ s, strong variations
of the maximum and mean total outflow velocity at the exhaust appear, as shown
in the top left of Fig. 11 ranging between 15 km s-1 to 24 km s-1 for the
maximum outflow velocity. These peaks coincide with peaks in the flux of
helium species with respect to the total flux, as shown in the bottom left of
Fig. 11. We note that helium species are mostly ionized at the exhaust, as
illustrated in Fig. 8 on the right.

Figure 11: Distribution of the evolution in time of the maximum and mean total
outflow velocity (left) and the ratio of fluxes
$\phi(\text{He}+\text{He}^{+})/(\phi(\text{He}+\text{He}^{+})+\phi(\text{H}+\text{H}^{+}))$
(right) for MFMS at the exhaust. On the left figure, the full and dashed black
lines correspond to the maximum and mean total outflow velocity.

Taking into account the complexity of the chemical evolution, chromospheric
reconnection may play a role in helium enrichment. In short, at the laminar
phase, the ionization fraction for hydrogen is higher, and it is slowly
expelled. At that time, outflows are relatively small. However, during the
second phase, with stronger flows and plasmoids, the number of helium atoms is
larger and helium is ionized due to the temperature increase, resulting in an
enrichment of helium at the outflows. However, as mentioned above, in some
locations of the right-hand side of the reconnection, stronger flows also
could produce hydrogen enrichment depending on the history of the plasma
evolution inside the current sheet, so care must be taken when interpreting
these results.

##  6 Conclusions and Discussion

In order to study a magnetic reconnection event under upper chromospheric
conditions, we have focused on a multi-fluid multi-species (MFMS) MHD model
and compared it with a two-fluid (TF) MHD model using Ebysus code. The MFMS
model is an extension of the TF approach, generalized to any possible species
with different ionization levels, considered as separate fluids and coupled
altogether with collisions and ionization or recombination processes. Both
models include anisotropic electron thermal conductivity and radiative losses.
In this work, we have focused on two different mixtures: a helium-hydrogen
mixture (MFMS model) $\mathcal{M}$ and a hydrogen or two-fluid (TF model)
mixture $\mathcal{M}_{\text{TF}}$. By doing so, our strategy allows the
comparison of the classical two-fluids approach (see Leake et al., 2012, 2013;
Murtas et al., 2021) with a MFMS approach and highlight the role of the helium
species in a chromospheric reconnection event. Unlike previous studies of
chromospheric magnetic reconnection from the literature (see Leake et al.,
2012, 2013; Murtas et al., 2021), the collisional frequencies are not assumed
constant but are calculated accurately based on a formalism with collision
integrals consistent at the kinetic level with the 13N moment model, as
described with further details in Wargnier et al. (2022).

For the two approaches TF and MFMS, we have considered a perturbed Harris
current sheet in upper chromospheric conditions, inspired by Leake et al.
(2012); Leake & Linton (2013). We have compared and analyzed results in both
cases at different levels: the evolution of the total current sheet, the
velocity fields and decoupling of the particles, the composition and evolution
of the ionization levels in MFMS and TF, and the main heating processes.

For both cases, we have two different phases: 1- a laminar phase where an
elongated current sheet is formed and 2- a plasmoid instability phase where
many plasmoids are formed and are advected from the middle of the reconnection
region towards the outflow regions. Our calculations have shown that the
evolution of the current sheet and the structure of the total velocity fields
slightly differs between MFMS and TF. Indeed, in MFMS, the presence of helium
species slightly slows the flow, leading to slightly smaller velocities in
MFMS than in TF. Additionally, the average thickness of the current sheet and
plasmoids is larger in MFMS than in TF while the average length of the current
sheet is similar between the two.

During the laminar phase, outflow and inflow velocities are formed due to the
Lorentz force generated by the initial structure of the magnetic field. In
both MFMS and TF, velocities reach $20$ km s-1 and $3.5$ km s-1 in the $z$ and
$y$ directions, respectively. In both directions, in MFMS, our calculations
have shown that the drift velocity between H+ and He is much higher than any
other drift velocities. This decoupling is due to the fact that neutrals He
decouples from H and H+ when the thickness of the current sheet thins down to
the order of the mean free path between He and H+ or H, and, neutral He atoms
are not dragged into the current sheet as quickly as H and H+. The fact that
He decouples from H and H+ and not H with H+ is attributed to smaller
collisional rate for He-H+ or He-H interactions than for H-H+ interactions in
upper chromospheric conditions. We refer to Wargnier et al. (2022) for further
details about the magnitude of these collisional rates. The relatively large
He and H+ drift velocity in MFMS does not occur in the TF model since the
latter considers only hydrogen species, i.e., H and H+.

During the plasmoid instability phase, the amplitude of all the drift
velocities is increasing in both MFMS and TF. The amplitude of the drift
velocity associated with He-H+ interactions is the highest inside plasmoids.
In MFMS, the drift between H and H+ is stronger at the boundary of plasmoids
and inside the current sheet, while it is approximately three times smaller
and shows an oscillating pattern parallel to the current sheet in TF. In both
MFMS and TF, the decoupling between H and H+ occurs because the thickness of
the current sheet continues to thin down to the order of the mean free path
between H and H+, and, neutral H does not collide enough with H+ during the
plasmoid instability phase, similarly to He-H+ interactions in MFMS during the
laminar phase. The amplitude of the drift between H and H+ is smaller in TF
than in MFMS - meaning that the collisional coupling is stronger in TF than in
MFMS due to the presence of helium species, as described with further details
in Section 5.2.

In both models, the main heating processes are due to the drift velocities
between any possible particles, also written as
$Q_{\mathfrak{a}\mathcal{I},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}^{\textbf{u},\text{col}}$.
During the laminar phase, our calculations have shown that adding helium in
the model allows to considerably increase the temperature of the plasma in the
current sheet compared to TF. Inside the current sheet, the temperature is
approximately 28000 K in MFMS and 19000 K in TF. In this context, a sufficient
level of temperature is reached due to collisions associated with H-H+ in the
current sheet to ionize helium. Finally, collisions associated with He-He+ and
He-H+ produce a small additional heating in the core of plasmoids and allows
to maintain the temperature level to 60000 K whereas only 30000 K is reached
with TF. Note that since the electrons are strongly coupled by collisions with
H+ species, the maximum electron temperature reached is also at 60000 K.

In the MFMS approach, due to the weak coupling between neutral helium and
other species, strong thermal decoupling effects have been obtained at a later
time of the plasmoid instability phase. Our results have demonstrated that the
temperature of neutral helium inside plasmoids can increase drastically up to
$90\,000$ K while the temperatures of other species reach a maximum of
$70\,000$ K. Note that a potential improvement of these models would be to
include photoionization of helium species, as highlighted by Golding et al.
(2014, 2016).

By combining a model where the collisional rates have been calculated
accurately with an efficient numerical strategy, we have demonstrated that in
upper chromospheric magnetic reconnection the inclusion of helium species in
the system allows a considerable increase temperature of the plasma while the
flows are weaker than a model where only hydrogen species are considered. This
large increase of temperature compared to a classical two-fluid case is
associated with a weak coupling between neutral helium and hydrogen that leads
to a high energy release due to large drift velocities in the current sheet
and inside plasmoids. In addition, during the plasmoid instability phase,
inside the plasmoids, two efficient mechanisms can occur and allow an
additional increase of the temperature: the production of He+ due to
ionization and recombination processes and collisions with other species, and
thermal decoupling effects which can increase the temperature of neutral
helium species locally. These mechanisms are able to easily increase the
temperature of the plasma during the plasmoid instability regime of a magnetic
reconnection event from upper chromospheric temperatures ($16\,000$ K) to
transition region temperatures ($90\,000$ K) on typical spatial and temporal
scales of Megameters and minutes respectively.

During the laminar phase, the total ionization level is much higher in MFMS
(from 85% to 95%) than in TF (from 70%) to 80%) due to more efficient heating
mechanisms in MFMS, as described previously, leading to higher temperature in
the current sheet giving rise to a large increase of the total ionization
level, mostly populated by H+. During the plasmoid instability phase, in both
MFMS and TF, the total ionization level increases drastically in the current
sheet and plasmoids. In MFMS, both neutral hydrogen and helium species are
ionizing, thus, the plasmoids and current sheet are mostly populated by H+ and
He+. In TF, the total ionization level inside the plasmoids, mostly due to the
presence of H+ species, is slightly smaller than in MFMS and is ranging from
80 to 98%.

Note at a later time of the plasmoid instability phase, chemical fractionation
effects have been captured in MFMS. Indeed, as shown in Section 5.4, helium
and hydrogen species are separated inside the current sheet, characterized by
an increase of helium species (from 5 to 12%) and a decrease of hydrogen
species (from 98% to 89%), which are advected in the outflow directions, along
the $z$ axis. During the plasmoid instability phase, the temperature increase
inside the current sheet and results in the ionization of helium species
leading to helium species enrichment at the exhaust of the reconnection.

Our study suggest that reconnection could lead to an enrichment of helium
during the strong outflow regimes under certain conditions. This mechanism,
out of reach of single-fluid or two-fluid MHD approaches, might explain the
enhancement of helium species observed at upper layers of the solar atmosphere
associated to the switchbacks (Fargette et al., 2021; Bale et al., 2021) from
the upper chromosphere. Similarly, intriguing observations of enhanced Helium
abundance in CMEs have been linked to reconnection events in the chromosphere
(Fu et al., 2020). Our results suggest that in the upper layer of the solar
atmosphere, any event associated with a large variation or reconfiguration of
the magnetic field could in principle release a sufficient level of energy to
produce this chemical fractionation between hydrogen and helium species while
increasing the temperature of the plasma.However, a simulation of a full
stratified atmosphere with a helium-hydrogen mixture, including alpha
particles, while guaranteeing that most of the mean free paths in that mixture
are spatially resolved, would be necessary to further investigate this issue.

##  7 Acknowledgments

We gratefully acknowledge support by NASA contracts, NNX17AD33G,
80NSSC18K1285, 80NSSC20K1272, and NNG09FA40C (IRIS), NSF contract AST1714955.
Resources supporting this work were provided by the NASA High-End Computing
(HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames
Research Center. The simulations have been run on clusters from the Notur
project, and the Pleiades cluster through the computing project s1061, s2053
and s8305. We thankfully acknowledge the support of the Research Council of
Norway through contract 230938/F50 and through grants of computing time from
the Programme for Supercomputing.

##  Appendix A Definition of the MFMS resistivity in small drift velocities
asymptotic limit

By considering the definition of the electric field and electron velocity from
Eqs. 12 and 14 it is possible to provide a new definition of the resistivity
of the MFMS model in terms of collision integrals
$\Omega^{1,1}_{\mathfrak{e},\mathfrak{a}\mathcal{I}}(\mu_{\mathfrak{e},\mathfrak{a}\mathcal{I}}),\,\mathfrak{a}\mathcal{I}\in\mathcal{H}$
following Wargnier et al. (2022). Indeed, assuming that the mass of the
electrons is much smaller than the mass of any heavy particles, the
collisional frequencies between $\mathfrak{e}$ and any
$\mathfrak{a}\mathcal{I}\in\mathcal{H}$ reads

$\nu^{\text{col}}_{\mathfrak{e},\mathfrak{a}\mathcal{I}}=\tfrac{4}{3}\tfrac{m_{\mathfrak{e},\mathfrak{a}\mathcal{I}}}{m_{\mathfrak{e}}}n_{\mathfrak{a}\mathcal{I}}|\mathbb{u}^{\text{th}}_{\mathfrak{e},\mathfrak{a}\mathcal{I}}|\,\Omega^{1,1}_{\mathfrak{e},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}(\mu_{\mathfrak{e},\mathfrak{a}\mathcal{I}})\approx\tfrac{4}{3}n_{\mathfrak{a}\mathcal{I}}|\mathbb{u}^{\text{th}}_{\mathfrak{e}}|\,\Omega^{1,1}_{\mathfrak{e},\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}(\mu_{\mathfrak{e}}),$
(A1)

where the thermal speed is simply defined as
$|\mathbb{u}^{\text{th}}_{\mathfrak{e}}|=\sqrt{8/\pi\mu_{\mathfrak{e}}}$. If
we consider the definition of Eq. A1 and the third term of Eq. 12 we obtain

$-\frac{1}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\mathbb{R}^{\text{col}}_{\mathfrak{e}}=\tfrac{4}{3}\tfrac{m_{\mathfrak{e}}|\mathbb{u}^{\text{th}}_{\mathfrak{e}}|}{q^{2}_{\mathfrak{e}}\,x_{\mathcal{\hat{I}}}}\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{H}}x_{\mathfrak{a}\mathcal{I}}\Omega^{1,1}_{\mathfrak{e},\mathfrak{a}\mathcal{I}}(\mu_{\mathfrak{e}})\left[\sum_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}\in\mathcal{H}}n_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}q_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}\mathbb{u}_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}}-\mathbb{J}-n_{\mathfrak{e}}q_{\mathfrak{e}}\mathbb{u}_{\mathfrak{a}\mathcal{I}}\right].$
(A2)

where
$x_{\mathcal{\hat{I}}}=n_{\mathfrak{e}}/(n_{\mathfrak{e}}+\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{H}}n_{\mathfrak{a}\mathcal{I}})$
is the ionization fraction of the whole plasma and
$x_{\mathfrak{a}\mathcal{I}}=n_{\mathfrak{a}\mathcal{I}}/(n_{\mathfrak{e}}+\sum_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}\in\mathcal{H}}n_{\mathfrak{a}^{\prime}\mathcal{I}^{\prime}})$
is the molar fraction of particle $\mathfrak{a}\mathcal{I}\in\mathcal{H}$. Eq.
A2 can be interpreted as the resistive term of the MFMS model. In the limit
where drift velocities are smaller than the hydrodynamic velocity of the
plasma or if the plasma is fully ionized or if the plasma is composed only of
one type of ionized species Eq. A2 simplifies to

$\frac{1}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\mathbb{R}^{\text{col}}_{\mathfrak{e}}\approx\left[\tfrac{4}{3}\tfrac{m_{\mathfrak{e}}|\mathbb{u}^{\text{th}}_{\mathfrak{e}}|}{q^{2}_{\mathfrak{e}}\,x_{\mathcal{\hat{I}}}}\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{H}}x_{\mathfrak{a}\mathcal{I}}\Omega^{1,1}_{\mathfrak{e},\mathfrak{a}\mathcal{I}}(\mu_{\mathfrak{e}})\right]\mathbb{J}=\eta\mathbb{J},$
(A3)

where the resistivity is defined as
$\eta=\tfrac{4}{3}\tfrac{m_{\mathfrak{e}}|\mathbb{u}^{\text{th}}_{\mathfrak{e}}|}{q^{2}_{\mathfrak{e}}\,x_{\mathcal{\hat{I}}}}\sum_{\mathfrak{a}\mathcal{I}\in\mathcal{H}}x_{\mathfrak{a}\mathcal{I}}\Omega^{1,1}_{\mathfrak{e},\mathfrak{a}\mathcal{I}}(\mu_{\mathfrak{e}})$.
Therefore, we have been able to define the resistivity of the MFMS model in
terms of collision integrals for any possible mixture valid in the asymptotic
limit where where the drift velocities are assumed to be much smaller than the
hydrodynamic velocities of the MFMS plasma. Note that this resistivity
coefficient is always true for the TF case based on
$\mathcal{M}_{2}{\text{TF}}$.

##  Appendix B Compact form of the MFMS governing equations

The presented system 1,4,6,7 and 13 of the MFMS model can be rewritten in a
compact form. By doing so, we can easily distinguish the nature of the
different terms involved in the system. The system can be rewritten as

$\partial_{t}\mathrm{U}+\nabla\cdot\mathrm{F}\left(\mathrm{U}\right)+\nabla\cdot\mathrm{D}\left(\mathrm{U},\nabla\mathrm{U}\right)=\mathrm{S}^{\text{conv}}\left(\mathrm{U},\nabla\mathrm{U}\right)+\mathrm{S}^{\text{stiff}}\left(\mathrm{U}\right),$
(B1)

where $\mathrm{U}\in\mathbb{R}^{5\mathfrak{n}_{\mathrm{\mathcal{H}}}+4}$,
$\mathrm{F}\left(\mathrm{U}\right)\in\mathbb{R}^{3\times(5\mathfrak{n}_{\mathrm{\mathcal{H}}}+4)}$,
$\mathrm{D}\left(\mathrm{U},\nabla\mathrm{U}\right)\in\mathbb{R}^{3\times(5\mathfrak{n}_{\mathrm{\mathcal{H}}}+4)}$,
$\mathrm{S}^{\text{conv}}\left(\mathrm{U},\nabla\mathrm{U}\right)\in\mathbb{R}^{5\mathfrak{n}_{\mathrm{\mathcal{H}}}+4}$
and
$\mathrm{S}^{\text{stiff}}\left(\mathrm{U}\right)\in\mathbb{R}^{5\mathfrak{n}_{\mathrm{\mathcal{H}}}+4}$.
These terms read

$\displaystyle\mathrm{U}$ $\displaystyle=$
$\displaystyle\left((\rho_{\mathfrak{a}\mathcal{I}})_{\mathfrak{a}\mathcal{I}\in\mathcal{H}},\;(\rho_{\mathfrak{a}\mathcal{I}}\mathbb{u}_{\mathfrak{a}\mathcal{I}}^{T})_{{\mathfrak{a}\mathcal{I}}\in\mathcal{H}},\;(e_{\mathfrak{a}\mathcal{I}})_{{\mathfrak{a}\mathcal{I}}\in\mathcal{H}},\;e_{\mathfrak{e}},\;\mathbb{B}^{T}\right)^{T},$
$\displaystyle\mathrm{F}\left(\mathrm{U}\right)$ $\displaystyle=$
$\displaystyle\left((\rho_{\mathfrak{a}\mathcal{I}}\mathbb{u}_{\mathfrak{a}\mathcal{I}})_{{\mathfrak{a}\mathcal{I}}\in\mathcal{H}},\;(\rho_{\mathfrak{a}\mathcal{I}}\mathbb{u}_{\mathfrak{a}\mathcal{I}}\otimes\mathbb{u}_{\mathfrak{a}\mathcal{I}}+P_{\mathfrak{a}\mathcal{I}}\mathbf{\mathbb{I}})_{{\mathfrak{a}\mathcal{I}}\in\mathcal{H}},\;(e_{\mathfrak{a}\mathcal{I}}\mathbb{u}_{\mathfrak{a}\mathcal{I}})_{{\mathfrak{a}\mathcal{I}}\in\mathcal{H}},\;e_{\mathfrak{e}}\mathbb{u}_{\mathfrak{e}},\;\mathbf{\mathbb{I}}\wedge\left(\mathbb{u}_{\mathfrak{e}}\wedge\mathbb{B}\right)\right)^{T},$
$\displaystyle\mathrm{D}\left(\mathrm{U},\nabla\mathrm{U}\right)$
$\displaystyle=$
$\displaystyle\left(0_{\mathfrak{n}_{\mathrm{\mathcal{H}}}},\;0_{3\times\mathfrak{n}_{\mathrm{\mathcal{H}}}},\;0_{\mathfrak{n}_{\mathrm{\mathcal{H}}}},\;\mkern
1.5mu\overline{\mkern-1.5mu\mkern
1.5mu\overline{\mkern-1.5mu\kappa_{\mathfrak{e}}\mkern-1.5mu}\mkern
1.5mu\mkern-1.5mu}\mkern 1.5mu\cdot\nabla
k_{B}T_{\mathfrak{e}},\;0_{3}\right)$
$\displaystyle\mathrm{S}^{\text{conv}}\left(\mathrm{U},\nabla\mathrm{U}\right)$
$\displaystyle=$
$\displaystyle\Biggl{(}0_{\mathfrak{n}_{\mathrm{\mathcal{H}}}},\;\left(n_{\mathrm{\mathfrak{a}\mathcal{I}}}q_{\mathrm{\mathfrak{a}\mathcal{I}}}\left(\left[\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}+\mathbb{u}_{\mathfrak{e}}\right]\wedge\mathbb{B}-\frac{\nabla
P_{\mathfrak{e}}}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\right)\right)_{\mathfrak{a}\mathcal{I}\in\mathcal{H}},$
$\displaystyle\,\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad-\left(P_{\mathrm{\mathfrak{a}\mathcal{I}}}\nabla\cdot\mathbb{u}_{\mathrm{\mathfrak{a}\mathcal{I}}}\right)_{\mathfrak{a}\mathcal{I}\in\mathcal{H}},\;-P_{\mathfrak{e}}\nabla\cdot\mathbb{u}_{\mathfrak{e}},\;-\nabla\wedge\left(\frac{1}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\nabla
P_{\mathfrak{e}}\right)\biggr{)}$
$\displaystyle\mathrm{S}^{\text{stiff}}\left(\mathrm{U}\right)$
$\displaystyle=$
$\displaystyle\biggl{(}\left(m_{\mathfrak{a}\mathcal{I}}\Gamma^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+m_{\mathfrak{a}\mathcal{I}}\Gamma^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}\right)_{\mathfrak{a}\mathcal{I}\in\mathcal{H}},\;\left(\mathbb{R}^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+\mathbb{R}^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+\mathbb{R}_{\mathrm{\mathfrak{a}\mathcal{I}}}^{\text{col}}+\left[\frac{n_{\mathfrak{a}\mathcal{I}}q_{\mathfrak{a}\mathcal{I}}}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\right]\mathbb{R}_{\mathrm{\mathfrak{e}}}^{\text{col}}\right)_{\mathfrak{a}\mathcal{I}\in\mathcal{H}},\;$
$\displaystyle\,\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\left(Q^{\text{ion}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+Q^{\text{rec}}_{\mathrm{\mathfrak{a}\mathcal{I}}}+Q_{\mathrm{\mathfrak{a}\mathcal{I}}}^{\text{col}}\right)_{\mathfrak{a}\mathcal{I}\in\mathcal{H}},\;Q^{\text{ion}}_{\mathrm{\mathfrak{e}}}+Q^{\text{rec}}_{\mathrm{\mathfrak{e}}}+Q_{\mathrm{\mathfrak{e}}}^{\text{col}},\;\nabla\wedge\left(\frac{\mathbb{R}_{\mathrm{\mathfrak{e}}}^{\text{col}}}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\right)\biggr{)}$

where $\mathrm{U}$ are the conservative variables of the system, $\mathrm{F}$
the convective fluxes, $\mathrm{D}$ the diffusive fluxes, $\mkern
1.5mu\overline{\mkern-1.5mu\mkern
1.5mu\overline{\mkern-1.5mu\kappa_{\mathfrak{e}}\mkern-1.5mu}\mkern
1.5mu\mkern-1.5mu}\mkern 1.5mu$ the anisotropic electron thermal conductivity,
$\mathrm{S}^{\text{conv}}$ the source terms as part of the convective system
and $\mathrm{S}^{\text{stiff}}$ the stiff source terms associated with
collisions and ionization/recombination processes,
$0_{\mathfrak{n}_{\mathrm{\mathcal{H}}}}$ is a null vector of size
$\mathfrak{n}_{\mathrm{\mathcal{H}}}$ and $\mathbf{\mathbb{I}}$ is the
identity matrix.

Note that in this system the only diffusive flux is the one associated with
the Spitzer term in the equation of thermal energy of the electrons. The
source terms $\mathrm{S}^{\text{conv}}$ have characteristic timescales which
typically can be much larger than $\mathrm{S}^{\text{stiff}}$. Indeed,
$\mathrm{S}^{\text{stiff}}$ have timescales at the order of the collisional,
ionization and recombination rates between all the particles
$\mathfrak{a}\mathcal{I}\in\mathcal{M}$. One can notice that the Lorentz force
terms in momentum equations have been split into two parts, one in the
$\mathrm{S}^{\text{conv}}$ part and the other one in
$\mathrm{S}^{\text{stiff}}$. This is because the part of the Lorentz force as
$\left[\frac{n_{\mathfrak{a}\mathcal{I}}q_{\mathfrak{a}\mathcal{I}}}{n_{\mathfrak{e}}q_{\mathfrak{e}}}\right]\mathbb{R}_{\mathrm{\mathfrak{e}}}^{\text{col}}$
has a characteristic timescale which is at the order of the collisional time
between electrons and heavy-particles, which can be much smaller than the
other terms of the Lorentz force.

## References

  * Abdulle & Vilmart (2013) Abdulle A., Vilmart G., 2013, Journal of Computational Physics, 242, 869, PIROCK: A swiss-knife partitioned implicit–explicit orthogonal Runge–Kutta Chebyshev integrator for stiff diffusion–advection–reaction problems with or without noise
  * Aldrovandi & Pequignot (1973) Aldrovandi S. M. V., Pequignot D., 1973, A&A, 25, 137, Radiative and Dielectronic Recombination Coefficients for Complex Ions
  * Alvarez Laguna et al. (2016) Alvarez Laguna A., Lani A., Deconinck H., Mansour N. N., Poedts S., 2016, Journal of Computational Physics, 318, 252, A fully-implicit finite-volume method for multi-fluid reactive and collisional magnetized plasmas on unstructured meshes
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481, The Chemical Composition of the Sun
  * Bale et al. (2021) Bale S. D., Horbury T. S., Velli M., et al., 2021, The Astrophysical Journal, 923, 174, A Solar Source of Alfvénic Magnetic Field Switchbacks: In Situ Remnants of Magnetic Funnels on Supergranulation Scales
  * Ballester et al. (2018) Ballester J. L., Alexeev I., Collados M., et al., 2018, Space Sci. Rev., 214, 58, Partially Ionized Plasmas in Astrophysics
  * Bharti, L. et al. (2013) Bharti, L. , Hirzberger, J. , Solanki, S. K. , 2013, A&A, 552, L1, Fine structures in the atmosphere above a sunspot umbra
  * Bhattacharjee et al. (2009) Bhattacharjee A., Huang Y.-M., Yang H., Rogers B., 2009, Physics of Plasmas, 16, 112102, Fast reconnection in high-Lundquist-number plasmas due to the plasmoid Instability
  * Bhattacharjee et al. (2010) Bhattacharjee A., Sullivan B., Huang Y.-M., 2010, Current sheet formation and the plasmoid instability in large, hyperresistive Hall MHD systems, en APS Division of Plasma Physics Meeting Abstracts, Vol. 52, p. CP9
  * Brackbill & Barnes (1980) Brackbill J. U., Barnes D. C., 1980, Journal of Computational Physics, 35, 426, The Effect of Nonzero div B on the numerical solution of the magnetohydrodynamic equations
  * Bruno et al. (2010) Bruno D., Catalfamo C., Capitelli M., et al., 2010, Physics of Plasmas, 17, 112315, Transport properties of high-temperature Jupiter atmosphere components
  * Cargill (2013) Cargill P., 2013, Astronomy and Geophysics, 54, 3.16, From flares to nanoflares: magnetic reconnection on the Sun
  * Cassak et al. (2009) Cassak P. A., Shay M. A., Drake J. F., 2009, Physics of Plasmas, 16, 120702, Scaling of Sweet–Parker reconnection with secondary islands
  * Cook & Cabot (2005) Cook A. W., Cabot W. H., 2005, Journal of Computational Physics, 203, 379, Hyperviscosity for shock-turbulence interactions
  * De Pontieu et al. (2009) De Pontieu B., McIntosh S. W., Hansteen V. H., Schrijver C. J., 2009, ApJ, 701, L1, Observing the Roots of Solar Coronal Heating in the Chromosphere
  * De Pontieu et al. (2014a) De Pontieu B., Rouppe van der Voort L., McIntosh S. W., et al., 2014a, Science, 346, D315, On the prevalence of small-scale twist in the solar chromosphere and transition region
  * De Pontieu et al. (2014b) De Pontieu B., Title A. M., Lemen J. R., et al., 2014b, Sol. Phys., 289, 2733, The Interface Region Imaging Spectrograph (IRIS)
  * Del Zanna, G. et al. (2015) Del Zanna, G. , Dere, K. P. , Young, P. R. , Landi, E. , Mason, H. E. , 2015, A&A, 582, A56, CHIANTI - An atomic database for emission lines. Version 8
  * Descombes & Massot (2004) Descombes S., Massot M., 2004, Numerische Mathematik, 97, 667, Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: Singular perturbation and order reduction
  * Duarte et al. (2013) Duarte M., Descombes S., Tenaud C., Candel S., Massot M., 2013, Combustion and Flame, 160, 1083, Time-space adaptive numerical methods for the simulation of combustion fronts
  * Fargette et al. (2021) Fargette N., Lavraud B., Rouillard A. P., et al., 2021, The Astrophysical Journal, 919, 96, Characteristic Scales of Magnetic Switchback Patches Near the Sun and Their Possible Association With Solar Supergranulation and Granulation
  * Fu et al. (2020) Fu H., Harrison R. A., Davies J. A., et al., 2020, ApJ, 900, L18, The High Helium Abundance and Charge States of the Interplanetary CME and Its Material Source on the Sun
  * Golding et al. (2014) Golding T. P., Carlsson M., Leenaarts J., 2014, The Astrophysical Journal, 784, 30, Detailed and simplified nonequilibrium helium ionization in the solar atmosphere
  * Golding et al. (2016) Golding T. P., Leenaarts J., Carlsson M., 2016, The Astrophysical Journal, 817, 125, Non-equilibrium helium ionization in an mhd simulation of the solar atmosphere
  * Golding, T. P. et al. (2017) Golding, T. P. , Leenaarts, J. , Carlsson, M. , 2017, A&A, 597, A102, Formation of the helium extreme-UV resonance lines
  * Gontikakis, C. et al. (2013) Gontikakis, C. , Winebarger, A. R. , Patsourakos, S. , 2013, A&A, 550, A16, Spectral diagnostic of a microflare. Evidences of resonant scattering in C IV 1548 Å, 1550 Å lines
  * Gosling (2007) Gosling J., 2007, The Astrophysical Journal, 671, L73, Observations of magnetic reconnection in the turbulent high-speed solar wind
  * Graille et al. (2009) Graille B., Magin T. E., Massot M., 2009, Mathematical Models and Methods in Applied Sciences, 19, 527, Kinetic theory of plasmas: translational energy
  * Gudiksen et al. (2011) Gudiksen B. V., Carlsson M., Hansteen V. H., et al., 2011, A&A, 531, A154, The stellar atmosphere simulation code Bifrost. Code description and validation
  * Guo et al. (2020) Guo L.-J., Pontieu B. D., Huang Y.-M., Peter H., Bhattacharjee A., 2020, The Astrophysical Journal, 901, 148, Observations and Modeling of the Onset of Fast Reconnection in the Solar Transition Region
  * Hairer & Wanner (1996) Hairer E., Wanner G., 1996, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Vol. 14
  * Hairer & Wanner (1999) Hairer E., Wanner G., 1999, Journal of Computational and Applied Mathematics, 111, 93, Stiff differential equations solved by Radau methods
  * Hansteen et al. (2019) Hansteen V., Ortiz A., Archontis V., et al., 2019, A&A, 626, A33, Ellerman bombs and UV bursts: transient events in chromospheric current sheets
  * Hansteen et al. (2017) Hansteen V. H., Archontis V., Pereira T. M. D., et al., 2017, ApJ, 839, 22, Bombs and Flares at the Surface and Lower Atmosphere of the Sun
  * Hansteen et al. (1997) Hansteen V. H., Leer E., Holzer T. E., 1997, ApJ, 482, 498, The Role of Helium in the Outer Solar Atmosphere
  * Huang & Bhattacharjee (2010) Huang Y.-M., Bhattacharjee A., 2010, Physics of Plasmas, 17, 062104, Scaling laws of resistive magnetohydrodynamic reconnection in the high-Lundquist-number, plasmoid-unstable regime
  * Huang & Bhattacharjee (2010) Huang Y.-M., Bhattacharjee A., 2010, Physics of Plasmas, 17, 062104, Scaling laws of resistive magnetohydrodynamic reconnection in the high-Lundquist-number, plasmoid-unstable regime
  * Huang & Bhattacharjee (2016) Huang Y.-M., Bhattacharjee A., 2016, ApJ, 818, 20, Turbulent Magnetohydrodynamic Reconnection Mediated by the Plasmoid Instability
  * Huang et al. (2017) Huang Y.-M., Comisso L., Bhattacharjee A., 2017, The Astrophysical Journal, 849, 75, Plasmoid Instability in Evolving Current Sheets and Onset of Fast Reconnection
  * Hyman et al. (1979) Hyman J., Vichnevtsky R., Stepleman R., 1979, Adv. in Comp. Meth, PDE’s-III, 313
  * Janev et al. (1989) Janev R., Harrison M., Drawin H., 1989, Nuclear Fusion, 29, 109, Atomic and molecular database for fusion plasma edge studies
  * Judge (2007) Judge P., 2007, The HAO Spectral Diagnostic Package for Emitted Radiation (HAOS-DIPER) Reference Guide (Version 1.0)
  * Kasper et al. (2015) Kasper J. C., Abiad R., Austin G. W., et al., 2015, Space Science Reviews, 204, 131, Solar Wind Electrons Alphas and Protons (SWEAP) Investigation: Design of the Solar Wind and Coronal Plasma Instrument Suite for Solar Probe Plus
  * Khomenko et al. (2014) Khomenko E., Collados M., Díaz A., Vitas N., 2014, Physics of Plasmas, 21, 092901, Fluid description of multi-component solar partially ionized plasma