# The Impact of Multi-fluid Effects in the Solar Chromosphere on the
Ponderomotive Force Under LTE and NEQ Ionization Conditions

Juan Martínez-Sykora Lockheed Martin Solar & Astrophysics Laboratory, 3251
Hanover Street, Palo Alto, CA 94304, USA  Bay Area Environmental Research
Institute, NASA Research Park, Moffett Field, CA 94035, USA.  Rosseland Center
for Solar Physics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo,
Norway  Institute of Theoretical Astrophysics, University of Oslo, P.O. Box
1029 Blindern, NO-0315 Oslo, Norway  Bart De Pontieu Lockheed Martin Solar &
Astrophysics Laboratory, 3251 Hanover Street, Palo Alto, CA 94304, USA
Rosseland Center for Solar Physics, University of Oslo, P.O. Box 1029
Blindern, NO-0315 Oslo, Norway  Institute of Theoretical Astrophysics,
University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway  Viggo H.
Hansteen Lockheed Martin Solar & Astrophysics Laboratory, 3251 Hanover Street,
Palo Alto, CA 94304, USA  Bay Area Environmental Research Institute, NASA
Research Park, Moffett Field, CA 94035, USA.  Rosseland Center for Solar
Physics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway
Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029
Blindern, NO-0315 Oslo, Norway  Paola Testa Harvard-Smithsonian Center for
Astrophysics, 60 Garden Street, Cambridge, MA 02193, USA  Q. M. Wargnier
Lockheed Martin Solar & Astrophysics Laboratory, 3251 Hanover Street, Palo
Alto, CA 94304, USA  Bay Area Environmental Research Institute, NASA Research
Park, Moffett Field, CA 94035, USA.  Mikolaj Szydlarski Rosseland Center for
Solar Physics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo,
Norway  Institute of Theoretical Astrophysics, University of Oslo, P.O. Box
1029 Blindern, NO-0315 Oslo, Norway

###### Abstract

The ponderomotive force is suggested to be the main mechanism to produce the
so-called first ionization potential (FIP) effect — the enrichment of low FIP
elements observed in the outer solar atmosphere. It is well known that the
ionization of these elements occurs within the chromosphere. Therefore, this
phenomenon is intimately tied to the plasma state in the chromosphere and the
corona. In addition, the chromosphere is a highly complex region with a large
variation in the ion-neutral collision frequencies, and hydrogen and helium
ionization is largely out of equilibrium. For this study, we combine IRIS
observations, a single fluid 2.5D radiative magnetohydrodynamics (MHD) model
of the solar atmosphere, including ion-neutral interaction effects and non-
equilibrium (NEQ) ionization effects, and a novel multi-fluid multi-species
(MFMS) numerical code (Ebysus). Nonthermal velocities of Si IV measured from
IRIS spectra can provide an upper limit for the strength of any high-frequency
Alfvén waves. With the single-fluid model, we investigate the possible impact
of NEQ ionization within the region where the FIP may occur and the plasma
properties in those regions. These models suggest that regions with strong
enhanced network and type II spicules are associated with the presence of
large ponderomotive forces. We use the plasma properties from the single-fluid
MHD model and the IRIS observations to initialize our multi-fluid models to
investigate the multi-fluid effects on the ponderomotive force associated with
Alfvén waves. Our multi-fluid analysis reveals that collisions and NEQ
ionization effects dramatically impact the behavior of the ponderomotive force
in the chromosphere, and existing theories may need to be revisited.

Magnetohydrodynamics (MHD) —Methods: numerical — Radiative transfer — Sun:
atmosphere — Sun: Corona

##  1 Introduction

The chromosphere is where high FIP elements remain neutral while low FIP
elements are ionized. It is here that chemical fractionation in the solar
atmosphere must begin. Observations suggest that neutral elements remain well
mixed with hydrogen, while ionized elements are preferentially pushed upwards
(Meyer, 1985; Feldman, 1992; Testa, 2010; Testa et al., 2015). Other stars
also exhibit fractionation effects in their coronal composition compared with
photospheric abundances, and their coronal chemical composition shows some
dependence on the stellar activity level (see reviews by Testa 2010; Testa et
al. 2015, and references therein). To further understand this process, one has
to keep in mind that the chromosphere is highly dynamic because of shocks and
magnetic field dynamics, and the magnitude of the magnetic field and density
rapidly decreases with height (e.g., Carlsson & Stein, 1992). In such a
dynamic and highly structured chromosphere, ionization shows strong gradients
(of many orders of magnitude) and is not in equilibrium. Instead, ionization
depends on the history of the plasma due to the non-equilibrium effects
(Leenaarts et al., 2007; Golding et al., 2016; Martínez-Sykora et al., 2020).
This means that ion-neutral collision frequencies vary by many orders of
magnitude on small scales within the solar atmosphere (e.g., Martínez-Gómez et
al., 2015; Martínez-Sykora et al., 2020; Nóbrega-Siverio et al., 2020;
Khomenko et al., 2021; Rempel & Przybylski, 2021; Przybylski et al., 2022;
Wargnier et al., 2022).

The chemical composition of the solar wind is also an indicator of the source
region on the Sun (Geiss et al., 1995) and is critical to establishing the
magnetic connectivity from the wind to the surface. A well-known example
corresponds to in-situ measurements of helium abundances at 1AU that depend on
the solar wind speed and phase of the activity cycle (Kasper et al., 2012).
The smallest helium abundances are observed during solar minimum. These
dependencies point to mechanisms that affect the second most abundant
constituent of the solar corona and whose effectiveness changes depending on
the levels of activity on the Sun. Similarly, Landi & Testa (2015) have found
the coronal Ne/O abundance ratio to significantly vary ($\gtrsim\times 2$,
i.e., more than a factor of two) during the solar cycle and peaking at the
cycle minimum. The new NASA mission Parker Solar Probe (PSP, Fox et al., 2016)
and the Solar Orbiter (SO, Müller et al., 2020) approach our star and measure
the fields and particles much closer than 1 AU, where the plasma is less
evolved and still preserves many of the original properties (composition,
charge, kinetic properties, etc.). However, establishing direct links to the
surface (observed with remote sensing instruments) remains a major challenge.
One of the most striking recent discoveries by PSP is the Alfvénic magnetic
field reversals known as switchbacks in the inner heliosphere (Bale et al.,
2021). The structures appear to be separated on spatial scales associated with
the supergranular scales. Therefore, it is thought that these structures
originate at the base of the lower solar atmosphere. A possible mechanism to
drive the switchbacks is interchange reconnection in the lower corona (Bale et
al., 2022). However, the switchbacks are also linked with the enrichment of
alpha particles, depletion of electron temperature, and slow solar wind
(Fargette et al., 2021; Woolley et al., 2021). The chemical fractionation,
hence the origin of the enhancement of alpha particles, should occur in or
close to the chromosphere, where plasma is partially ionized, suggesting the
driving mechanism of these phenomena may be connected to this region. Studies
such as these highlight the challenges and opportunities of establishing the
connectivity between the solar surface and the solar wind. This is important
also because understanding the connectivity between stellar surfaces and their
winds, the connection to the drivers, and their impact on the heliosphere and
astrospheres is key to predicting space weather.

The ponderomotive force associated with strong Alfvén waves is currently
regarded as the most likely mechanism to explain the FIP effect as modeled for
instance by Laming (2004, 2015). Nonetheless, this model uses a simplified
approach based on 1D semi-empirical, static atmospheres from Avrett (2007).
These atmospheric models estimate the ionization state from averaged
observations or assuming statistical equilibrium without coupling to the
overlying wind. Such semi-empirical models do not capture the complex dynamics
of the chromosphere and, as such, are not a good test bed for wave propagation
and fractionation studies. Under those assumptions, and in order to reproduce
the observables in most solar targets, the model requires Alfvén waves that
are generated in the corona and that propagate (from above) into the
chromosphere and that are there reflected. This theory postulates that
nanoflares and reconnection in the coronal volume drive the downward
propagating Alfvén waves. Observational evidence for low-frequency Alfvén
waves, however, shows a predominance of upward propagating Alfvén waves that
originate in the chromosphere (e.g., De Pontieu et al., 2007; Okamoto & De
Pontieu, 2011). Similarly, observations in the corona indicate a predominance
of upward propagating Alfvén waves (Tomczyk et al., 2007; McIntosh et al.,
2011), or a mix of upward and downward waves in which the upward flux exceeds
the downward flux by a significant factor (McIntosh & De Pontieu, 2009). It is
of course possible that high frequency Alfvén waves, which are much more
difficult to detect due to observational limitations, somehow show
significantly different behavior.

Dahlburg et al. (2016) computed the ponderomotive acceleration from a 3D MHD
model with optically thin radiation and thermal conduction. This ponderomotive
acceleration in the model occurs at the footpoints of coronal loops as a
byproduct of wave-driven coronal heating. The first part of this present work
broadly extends this study on the ponderomotive force in realistic radiative
MHD simulations, which includes non-equilibrium ionization effects and ion-
neutral interaction effects (Section 4).

In the following, the ponderomotive force theory is briefly described in
Section 2. In Section 3 we describe how we use IRIS observations to estimate
an upper limit to the amplitude of high-frequency Alfvén waves (Section 3). As
mentioned, the results regarding the characterization of the ponderomotive
force in the radiative MHD model under LTE and NEQ conditions are described in
Section 4. We present the multi-fluid models where we investigate the
collisional effects on the ponderomotive acceleration for Alfvén waves
(Section 5). Finally, the manuscript ends with a conclusion and discussion
(Section 6).

##  2 Theoretical considerations regarding the ponderomotive force

The Lorentz force is a key force for understanding the dynamics of plasmas. It
is a straightforward term that acts on charged particles. However, it is not
always practical to use this term in complex scenarios such as when
electromagnetic waves or complex environments are involved. For such
conditions, the ponderomotive force is often used to simplify understanding of
the plasma dynamics. The ponderomotive force is based on the time average of
nonlinear forces that act in the presence of oscillating electromagnetic
fields. The calculation of this force is however rather cumbersome as it
includes non-linear aspects, and it usually involves quite limiting
assumptions, rendering the results approximate. There are a variety of
formulations of the ponderomotive force, which differ in the underlying
assumptions. All variations of the ponderomotive force are derived by a
nonlinear dependence on the amplitude of the electric field oscillations. This
force can be treated as the acceleration in a electron-ion plasma, or in ions,
as nicely described in Whitelam et al. (2002). They described the acceleration
along the magnetic field associated with the ponderomotive force resulting
from nonlinear propagation of circularly polarized electro-magnetic (CPEM)
waves in an electron-ion plasma as follows:

$F=-\left(\frac{\partial}{\partial
z}+\frac{k\Omega_{cj}}{\omega(\omega+\Omega_{cj})}\frac{\partial}{\partial
t}\right)\frac{q_{j}^{2}|\vec{E}|^{2}}{m\omega(\omega+\Omega_{cj})}~{}$ (1)

where $\vec{E}$ is the electric field; $\omega$, and $k$ are the wave
frequency and number, respectively; and $\Omega_{cj}$, $m_{j}$ and $q_{j}$ are
the gyrofrequency, mass, and charge for species $j$ ($j$ can be any ionized
species or electrons). Whitelam et al. (2002) also considered a case in which
the ion acceleration is produced by electro-magnetic electron-cyclotron (EMEC)
or ion-cyclotron Alfvén (EMICA) waves in which they assumed that the
ponderomotive force dominates over the background Lorentz force and pressure
gradients:

$F=-\frac{q_{i}^{2}}{m_{i}m_{e}\omega(\omega-|\Omega_{ce}|)}\left(\frac{\partial}{\partial
z}-\frac{k\Omega_{ce}}{\omega(\omega-|\Omega_{ce}|)}\frac{\partial}{\partial
t}\right)|\vec{E}|^{2}~{}$ (2)

For the case of a multi-ions species and the ponderomotive force from EMICA
waves, the $\alpha$ species is governed by:

$\displaystyle F$ $\displaystyle=$
$\displaystyle-\frac{q_{i\alpha}^{2}}{m_{i\alpha}^{2}\Omega_{ci\alpha}}\left[\frac{1}{\omega-\Omega_{ci\alpha}}\frac{\partial}{\partial
z}+\right.$
$\displaystyle\left.\frac{k}{\omega^{2}}\left(1-\frac{\Omega_{ci\alpha}^{2}}{(\omega-\Omega_{ci\alpha})^{2}}\right)\frac{\partial}{\partial
t}\right]|\vec{E}|^{2}~{}$

Note that the spatial and temporal components are considered in the literature
are sometimes considered as different forces. A nice summary of the various
calculations of different forces can be found in Lundin & Guglielmi (2006). In
short, and to list a few from that paper:

  * •

Lundin-Hultqvist force or magnetic moment pumping (MMP) is given by:

$F=-\frac{mc^{2}}{2}\frac{E^{2}}{B^{2}}\nabla_{||}\ln B$ (4)

(Lundin & Hultqvist, 1989). Here $\vec{B}$ is the magnetic field. This
ponderomotive force assumes a low frequency regime, that the electric field
oscillation is perpendicular to the magnetic field and a very slow variability
in $B$. It is always in the direction of decreasing B (independent of the wave
propagation direction).

  * •

Miller force or gradient ponderomotive force

$F=\frac{mc^{2}}{4B^{2}}\frac{\partial E^{2}}{\partial z}$ (5)

(Miller, 1958). This force is related to the spatial inhomogeneity of the wave
field. The derivation of the force assumes a static magnetic field with a
perpendicular oscillation between the electric field and magnetic field, where
the magnetic field perturbations are negligible compared to the ambient
magnetic field. The formula provided is for Alfvén waves with low frequency.
This force is positive (in direction of wave propagation) if the gradient of
the wave field is positive (i.e., increasing in the direction of wave
propagation).

  * •

Abraham force is

$F=\pm\frac{mc^{2}}{2v_{A}B^{2}}\frac{\partial E^{2}}{\partial t}.$ (6)

This description of the ponderomotive force neglects standing waves and
applies to waves traveling in the direction of the magnetic field. The
positive (negative) sign corresponds to a wave propagating in the direction
(opposite) to the magnetic field.

The chromosphere is partially ionized and none of those studies considered the
ponderomotive effects of damping of Alfvén waves from ion-neutral
interactions. Haerendel (1992) takes into account the presence of collisions
and the ponderomotive force has a very different dependence:

$F=\rho_{n}u_{x}^{2}\frac{\omega^{2}}{2v_{A}\nu_{ni}}~{}$ (7)

where $\rho_{n}$, and $\nu_{ni}$ are the neutral density, and collision
frequency between ions and neutrals, respectively. $v_{A}$ is the Alfvén
speed. Note that this force is a factor of two smaller than originally
reported by Haerendel (1992). The corrected formula is provided in de Pontieu
& Haerendel (1998). This formula is based on the WKB approximation of a weakly
damped wave and is averaged over a wave period. A derivation for strongly
damped waves is provided by (Song & Vasyliūnas, 2011). The nature of this
force is such that it is proportional to the damping of the wave field along
the direction of propagation, i.e., $F\sim-\frac{\partial E^{2}}{\partial z}$.
Note that the frequency dependence is different from Eq. 1, which decreases
with increasing frequency. The ion-neutral damping of the waves leads to a
damping length $L$ for the wave field:

$L=\frac{2v_{A}\nu_{ni}\rho_{t}}{\omega^{2}\rho_{n}}~{}$ (8)

where $\rho_{t}$ is the total mass density, and $\rho_{n}$ is the neutral mass
density.

As we will show with our multi-fluids models (Section 5.3), the many
assumptions underlying some of these equations are not necessarily easy to
justify for solar conditions or even the simplified scenarios we consider here
(no density stratification, dynamics, and only Alfvén waves in the low-
frequency regime).

##  3 IRIS Observations: non-thermal velocity

We have selected a decaying AR with plage and a large presence of spicules
near the limb to characterize the high-frequency Alfvén waves traveling along
the spicules, which will be used to constrain our models in the following
sections. The exposure time of the selected IRIS observation is such as can
add a constrain to the high-frequencies, i.e., as short as possible (see below
for details).

###  3.1 Target and data calibration

For this study, we use IRIS observations acquired from 2014-12-11 22:39:14 to
23:32:274 in a plage region near the limb $(x,y)=(-910,-233)$″. IRIS observed
this active region (AR) with a very large sparse 64-step raster, with no roll,
exposures of 2 s, a cadence of 3.2 s, raster-scans of 199 s, and 2x2 binning
(i.e., a spatial plate scale of 0.33″and spectral sampling of 5.4 km s-1),
using OBSID 3820105390. The observation includes fourteen full rasters. This
study focuses only on the 13${}^{\text{th}}$ raster, which is free of
radiation hits from passage through the South Atlantic Anomaly (SAA). The slit
jaw image (SJI) data are obtained in the 1330 and 2796 passbands. We use IRIS
calibrated level 2 data, which have been processed for dark current, flat
field, and geometrical corrections (De Pontieu et al., 2014).

###  3.2 Analysis of the IRIS observations

We fit Si IV 1339 Å profiles with a single Gaussian to extract the Doppler
shifts and non-thermal widths (see Fig 1). The non-thermal velocity is
determined from the 1/e width of the single Gaussian fit by taking the square
root of the square of the 1/e width after subtracting the squares of the
instrumental and thermal contributions (see, Testa et al., 2016). Since the
observation is near the limb, the non-thermal width contains mostly unresolved
velocities that are expected to be almost perpendicular to the magnetic field.
One can recognize spicules and moss in the SJI and Si IV line intensity map
(left two panels). In plage and spicules, most non-thermal velocities range
from $10$ to $25$ km s-1. Note that the exposure time is 2 s, so the
assumption that the unresolved velocity comes from Alfvén waves suggests that
they may have high frequencies (1 s or smaller) with amplitudes within
$[15,36]$ km s-1. We note that there is also plentiful observational evidence
for Alfvén waves at lower frequencies (De Pontieu et al., 2007; Okamoto & De
Pontieu, 2011; McIntosh et al., 2011; Tomczyk et al., 2007).

Figure 1:  From left to right, IRIS SJI 1330 at 23:28:32, the total intensity,
Doppler shift, and non-thermal width from Si IV 1339 Å from fitting a single
Gaussian. The contours correspond to the SJI 1330 at 100 DN. The vertical red
lines in the left panel limit the region of the raster scan.

##  4 Characterizing the ponderomotive acceleration and chemical fractionation
in rMHD simulation: LTE vs NEQ

This section briefly describes the single fluid radiative MHD numerical model
and plasma properties in the simulated chromosphere, including the
ponderomotive force (calculated after the simulation was run), to understand
the conditions under which a FIP effect could possibly occur.

###  4.1 Radiative MHD numerical model including NEQ ionization and ambipolar
diffusion

For the forward analysis, we used an already analyzed radiative MHD model
computed with Bifrost (Gudiksen et al., 2011), which includes non-equilibrium
ionization for hydrogen and helium (Leenaarts et al., 2007; Golding et al.,
2016) as well as ambipolar diffusion (Nóbrega-Siverio et al., 2020). In short,
this 2.5D numerical model mimics two plage regions connected by $\sim 40$ Mm
long loops and reveals features resembling type I and II spicules, low-lying
loops, and other features. The simulation spans a vertical domain stretching
from $\sim 3$ Mm below the photosphere to $40$ Mm above into the corona with a
non-uniform vertical grid size of $12$ km in the photosphere and chromosphere
and $14$ km grid size in the horizontal axis. For further details about the
numerical model, we refer to Martínez-Sykora et al. (2020). We have chosen
this model because it includes NEQ effects and produces type II spicules and
the associated Alfvén waves (Martínez-Sykora et al., 2017).

###  4.2 Characterizing the simulated chromospheric properties for the
chemical fractionation

Figure 2:  Ionization fractions are shown for various species with low and
high FIP (hydrogen in panels a and c, helium in panels b and d, iron in panel
e, silicon in panel f, and calcium in panel g). Panels a and b assume LTE, and
panels c and d are in NEQ for hydrogen and helium. Note the color bar is in
log-scale. Labels in the panels show the FIP energy in eV.

In this numerical model, we computed the ionization fraction for various
species with both low (calcium, iron, silicon) and high (hydrogen and helium)
FIP (Figure 2). We derive from the model the ionization fraction for hydrogen
and helium in LTE (panels a and b) as well as in NEQ ionization (panels c and
d). We assumed photospheric abundances (Asplund et al., 2009), and statistical
and thermal equilibrium for the ionization and density populations for
neutrals and ions in the LTE cases. Low FIP species, e.g., calcium, iron, and
silicon, are highly ionized in a large fraction of the chromosphere in LTE,
whereas high FIP species, including hydrogen and helium, are mostly neutral in
the chromosphere in LTE. However, NEQ ionization effects significantly
increase the ionization fraction for hydrogen (panel c) and helium (panel d)
and lead to a larger height range where these species are partially ionized
(Leenaarts et al., 2007; Golding et al., 2016).

Figure 3:  The ratio of the ionization fraction of hydrogen with other species
(helium, iron, silicon, and calcium from top to bottom) reveals where the
relative chemical fractionation may happen. Red means hydrogen is ionized,
whereas the other is neutral and blue vice-versa. The left column is in LTE,
and the right column NEQ for hydrogen and helium.

A comparison between hydrogen ionization fraction and any other species
ionization fraction allows us to visualize the region of interest where the
first ionization occurs for the various species (for low and high FIP) and
chemical fractionation may occur. Figure 3 shows the ratio of hydrogen
ionization fraction with the ionization fraction of any other species. We
considered this comparison in LTE (left column) and NEQ ionization for
hydrogen and helium (right column). The red or blue regions are regions of
interest in terms of fractionation since in these regions the ionization
degree between low and high FIP elements is significantly different (either
higher or lower, respectively). If a physical process like the ponderomotive
force were to act in these regions, it would potentially lead to relative
enrichment (or depletion) of elements with different FIP. In addition, in NEQ,
the ionization degree may be inverted in extended areas (red color), i.e., the
hydrogen is ionized when the other species are neutral, e.g., panels f-i).

Figure 4:  Collision frequency between the neutral hydrogen and the other
ionized species (helium, iron, silicon, and calcium from top to bottom) maps
are shown in LTE (left column) and NEQ ionization (right column).

However, we note that to have chemical fractionation, a different degree of
ionization between the various species is not sufficient since: 1. a physical
process needs to act on the various species (e.g., ponderomotive force) to
accelerate the particles upward, 2. the dynamics of the various species need
to be somewhat decoupled from the timescales involved in the acceleration
(i.e., collisions between high and low FIP elements cannot dominate since
these can couple the dynamics of the various elements). We compute the
collisional rates between hydrogen, which is the most abundant fluid and will
provide a good estimate of coupling, and other ionized species, again,
considering LTE (left) and NEQ (right) in Figure 4. The upper chromosphere and
transition region are weakly collisional. Weakly collisional regions is where
collisional timescales are a significant fraction of the wave period.
Therefore, we consider weakly collisional is for collision frequencies lower
than 10 Hz, and the driver of the chemical fractionation may occur in those
locations. Note that in the simulation, the collision frequencies for ionized
iron, silicon, and calcium are the same in NEQ and LTE because they do not
depend on hydrogen number density. Their contribution to the momentum exchange
nevertheless depends on the hydrogen number density (e.g., see Wargnier et
al., 2022).

A key aspect is where in the model the assumed physical process (ponderomotive
force) is likely to be important. Dahlburg et al. (2016) estimated the
ponderomotive acceleration in their 3D model considering the electric field
perturbations along the magnetic field. Similarly, we compute the electric
field perturbation along the magnetic field as follows:

$\displaystyle\delta_{s}E=\frac{1}{B^{2}}\frac{\partial E^{2}}{\partial
s},~{}$ (9)

where $s$ is along the loop. We mask regions with high collision rates ($\geq
10$ Hz) for each element as well as regions where they have the same
ionization state as hydrogen (Figure 5). By masking out these regions, we get
the region of interest (the non-masked region), where the action of
$|\delta_{s}E|$ may lead to chemical fractionation for each species. Different
structures show different $|\delta_{s}E|$, with type II spicules experiencing
the largest values. Note that the region of interest changes depending on
which species are considered. Those masks reveal the complexity of this
problem due to the highly structured and dynamic atmosphere.

Figure 5:  Maps of the electric field variations along the loop (Eq. 9). In
these maps, we have masked out regions with high collision rates ($>10$ Hz)
for each element (Figure 4) as well as regions where they have the same
ionization state as hydrogen (Figure 3). The region of interest (non-masked
region) is where the ponderomotive force may produce chemical fractionation.
The distribution of the various panels is the same as in Figure 4.

Since the models in Dahlburg et al. (2016) were driven by waves, expression 9
estimates the ponderomotive acceleration. This is not necessarily true in
self-consistent radiative MHD since many other physical processes change the
electric field and not only waves. To obtain a broader view, we show in Figure
6, in addition to $\delta_{s}E$, the Alfvén wave power (by computing the
velocity projections onto three characteristic directions, e.g., Khomenko &
Cally, 2011); $|\vec{J}|/|\vec{B}|$, which reveals regions with strong tension
or reconnection; the Poynting flux; the Lorentz force; and the advective part
of the electric field $|\vec{u}\times\vec{B}|$. One can appreciate that there
is not a single physical process that correlates or connects directly with
$\delta_{s}E$ or Alfvén waves. Some locations with high $\delta_{s}E$ are
associated with several other physical processes. Still, regions with spicules
show either the largest values of $|\vec{J}|/|\vec{B}|$ or Lorentz force
($x=20$, 35 or 40 Mm) as well as $\delta_{s}E$. On top of this, those regions
generate the strongest Alfvén wave power that penetrates the corona. We find
that regions with large $|\vec{J}|/|\vec{B}|$ are associated with spicules or
low-lying loops interacting with small and large-scale fields. The field lines
connected to large $|\vec{J}|/|\vec{B}|$, and hence large $\delta_{s}E$ in the
chromosphere, shows large Alfvén wave amplitude in the transition region and
corona. Note that these regions are connected to regions with enhanced
magnetic field networks or plage regions. These then are the regions that seem
to harbor conditions that are favorable for the ponderomotive force to act.
Observationally these are often regions that show enhanced non-thermal line
broadening, as we show in Section 3.2. In principle, high-frequency Alfvén
waves can produce large non-thermal velocities of transition region lines.
This is because these waves are often unresolved at the temporal and spatial
resolution of current instruments (see Section 3.2).

Figure 6:  Various physical processes play a different role in different
features. $\delta_{s}E$ following Eq 9, the Alfvén wave power, $|J|/|B|$,
Poynting flux, Lorentz force and advection term of the electric field are
shown from left to right and top to bottom, respectively. The solid line is at
temperature $10^{5}$ K.

In order to initialize our MFMS simulations, we compute within the regions of
interest (the non-masked regions in Figure 5), the populations of each species
(four left panels of Figure 7), magnetic field (panels e and f) and
temperature (panels g and h). The top row is under LTE conditions, and the
bottom row is for NEQ. Panels a, b, e-h are within the regions of interest for
each species, and panels c and d are for the region of interest of H-Ca
(panels d and h in Figure 5). Note that the regions of interest are where the
chemical fractionation may happen since it is relatively weakly collisional
and the ionization degree differs between high and low FIP elements. To
inspire the initial conditions of our MFMS simulations (see below), we use the
region of interest of H-Ca, because it covers the largest region where
chemical fractionation may happen. From the models it appears that the field
strength ranges from $\sim 10$ G to $\sim 120$ G and the median is around 30 G
in those regions where chemical fractionation may occur. For the temperature,
there is a clear difference between LTE and NEQ scenarios; the former has a
narrower temperature range ($\log(T(K))=[3.4,4]$) than in NEQ
($\log(T(K))=[3.2,4.3]$) and both roughly peaking at similar temperature
($\log(T(K))\sim 3.7$), but for LTE this varies more for each species.

Figure 7:  Histogram of the density population for neutrals (solid, panels
a-d) and ions (dashed, panels a-d), magnetic field (third column) and
temperature (right column) within the regions of interest (non-masked out
regions) in Figure 5, for each species (panels a and b, and e-h). Panels c and
d are for the region of interest of H-Ca. LTE case is in the top row and NEQ
in the bottom row.

For the following study (Section 5) and to limit the parameter range of the
multi-fluid simulations (which are computationally expensive), we restrict the
variation of the magnetic field along the loop, and initialize the densities
as constant in space (Section 2). This means that we consider 1D loops that
expand or constrict with height (and/or direction of wave propagation). One
would expect that on the Sun the degree of expansion with height of the field
would depend on the type of region, e.g., coronal holes, quiet sun, plage,
active regions, etc. We compute the magnetic field strength and Alfvén speed
in the region of interest shown in panel h of Figure 5 from the 2.5D numerical
model and find that in our plage simulation the configuration is much more
complex than a simple expansion with height, and this would need to be taken
into account in the future (Figure 8). For now, inspired by the model field
topology, we will parameterize the expansion of the magnetic field as listed
in Table 2. The complex magnetic field expansion or canopies have roughly a
variation of the magnetic field of $\Delta B\sim\pm 20$ G.

Figure 8:  Magnetic field strength (top) and Alfvén speed (bottom) maps, where
we applied the same mask as in panel h in Figure 5, reveal the complexity of
the expansion or canopies of the magnetic field topology within the regions of
interest.

##  5 Multi-fluid and multi-species simulations of the ponderomotive force

In this second part of this manuscript, we introduce the multi-fluid multi-
species (MFMS) numerical code Ebysus, including a short description of the
physics included and the numerical scheme (Section 5.1). As a continuation, we
describe the setup of the parametric study of the various numerical models
(Section 5.2) and, finally, the results of the MFMS models (Section 5.3).

###  5.1 The MFMS numerical Ebysus code

The Ebysus code is the first of its kind to treat each excited/ionized level
for each desired species as a separate fluid. The multi-fluid equations in SI
implemented in the code are as follows:

$\displaystyle\frac{\partial\rho_{\mathrm{aIE}}}{\partial
t}+\nabla\cdot\left(\rho_{\mathrm{aIE}}\,{\vec{u}_{\mathrm{aIE}}}\right)$
$\displaystyle=$
$\displaystyle\sum_{I^{\prime}E^{\prime},a}m_{\mathrm{aIE}}(n_{\mathrm{aI^{\prime}E^{\prime}}}\Gamma^{ion}_{\mathrm{aI^{\prime}E^{\prime}IE}}-n_{\mathrm{aIE}}\Gamma^{rec}_{\mathrm{aIEI^{\prime}E^{\prime}}})$
(10)
$\displaystyle\frac{\partial(\rho_{\mathrm{a\hat{I}E}}{\vec{u}_{\mathrm{a\hat{I}E}}})}{\partial
t}+\nabla\cdot(\rho_{\mathrm{a\hat{I}E}}{\vec{u}_{\mathrm{a\hat{I}E}}}\otimes{\vec{u}_{\mathrm{a\hat{I}E}}}-\hat{\tau}_{\mathrm{a\hat{I}E}}+P_{\mathrm{a\hat{I}E}}\mathbf{\mathbb{I}})$
$\displaystyle=$
$\displaystyle\rho_{\mathrm{a\hat{I}E}}{\vec{g}}+n_{\mathrm{a\hat{I}E}}q_{\mathrm{a\hat{I}}}\left({\vec{E}}+{\vec{u}_{\mathrm{a\hat{I}E}}}\times{\vec{B}}\right)+$
(11)
$\displaystyle\sum_{\hat{I^{\prime}}E^{\prime},a}(\Gamma_{\mathrm{a\hat{I^{\prime}}E^{\prime}\hat{I}E}}^{ion}m_{\mathrm{a\hat{I}E}}\vec{u}_{\mathrm{a0E^{\prime}}}-\Gamma_{\mathrm{a\hat{I}EI^{\prime}E^{\prime}}}^{rec}m_{\mathrm{a\hat{I}E}}\vec{u}_{\mathrm{a\hat{I}E}})+\sum_{\mathrm{a^{\prime}I^{\prime}E^{\prime}}}\vec{R}_{\mathrm{a\hat{I}E}}^{\mathrm{a\hat{I}Ea^{\prime}I^{\prime}E^{\prime}}}$
$\displaystyle\frac{\partial(\rho_{\mathrm{a0E}}{\vec{u}}_{\mathrm{a0E}})}{\partial
t}+\nabla\cdot(\rho{\vec{u}}_{\mathrm{a0E}}\otimes{\vec{u}}_{\mathrm{a0E}}-\hat{\tau}_{\mathrm{a0E}}+P_{\mathrm{a0E}}\mathbf{\mathbb{I}})$
$\displaystyle=$
$\displaystyle\rho_{\mathrm{a0E}}{\vec{g}}+\sum_{\mathrm{a^{\prime}I^{\prime}E^{\prime}}}\vec{R}_{\mathrm{a0E}}^{\mathrm{a0Ea^{\prime}I^{\prime}E^{\prime}}}-$
(12)
$\displaystyle\sum_{\mathrm{I^{\prime}E^{\prime},a}}\left(\Gamma_{\mathrm{aI^{\prime}E^{\prime}}}^{ion}m_{\mathrm{a0E}}\vec{u}_{\mathrm{a0E}}+\Gamma_{\mathrm{a0E}}^{rec}m_{\mathrm{aI^{\prime}E^{\prime}}}\vec{u}_{\mathrm{aI^{\prime}E^{\prime}}}\right)$
$\displaystyle\frac{\partial{e_{\mathrm{aIE}}}}{\partial{t}}+\nabla\cdot\left(e_{\mathrm{aIE}}\vec{u}_{\mathrm{aIE}}\right)+P_{\mathrm{aIE}}\nabla\cdot\vec{u}_{\mathrm{aIE}}=$
$\displaystyle
Q^{visc}_{\mathrm{aIE}}+Q^{ir}_{\mathrm{aIE}}+Q^{\mathrm{aIEe}}_{\mathrm{aIE}}+\sum_{\mathrm{a^{\prime}I^{\prime}E^{\prime}}}^{\mathrm{a^{\prime}I^{\prime}E^{\prime}}\neq\mathrm{aIE}}Q_{\mathrm{aIE}}^{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}$
(13)
$\displaystyle\frac{\partial{e_{e}}}{\partial{t}}+\nabla\cdot\left(e_{e}\vec{u_{e}}\right)+P_{e}\nabla\cdot\vec{u_{e}}=Q^{visc}_{e}$
$\displaystyle+$ $\displaystyle
Q^{ir}_{e}+\sum_{\mathrm{a^{\prime}I^{\prime}E^{\prime}}}Q_{e}^{e\mathrm{a^{\prime}I^{\prime}E^{\prime}}}+Q_{e}^{\mathrm{spitz}}$
(14)

where the ionization states are referred as $\mathrm{I}$, i.e., $\mathrm{I}=0$
denotes neutrals and $\hat{\mathrm{I}}=\mathrm{I}\geq 1$ ions. The excited
levels are marked with $\mathrm{E}$ and the identity of the chemical species
(or molecules) is indicated by $\mathrm{a}$. Consequently, each set of
particles in a given micro-state will be described with $\mathrm{aIE}$. For
electrons, the notation $\mathrm{aIE}$ is reduced to just $\mathrm{e}$. For
simplicity, $\sum_{\mathrm{a^{\prime}}}$ is the sum over all the species
$\mathrm{a^{\prime}}$, $\sum_{\mathrm{I^{\prime},a}}$ is the sum over all
ionization levels, including neutrals, for a given species $\mathrm{a}$ and
$\sum_{\mathrm{E}^{\prime},\mathrm{aI}}$ is the sum over all the excited
levels for a given ionized species $\mathrm{aI}$. For clarity, we define
$\sum_{\mathrm{a^{\prime}I^{\prime}E^{\prime}}}=\sum_{\mathrm{a^{\prime}}}\sum_{\mathrm{I^{\prime},a^{\prime}}}\sum_{\mathrm{E^{\prime},a^{\prime}I^{\prime}}}$,
and
$\sum_{I^{\prime}E^{\prime},a}=\sum_{\mathrm{I^{\prime},a}}\sum_{\mathrm{E^{\prime},aI^{\prime}}}$.
The mass density is $\rho_{\mathrm{aIE}}=m_{\mathrm{aIE}}\,n_{\mathrm{aIE}}$,
where ${\vec{u}_{\mathrm{aIE}}}$, $n_{\mathrm{aIE}}$ and $m_{\mathrm{aIE}}$
are the velocity, number density and particle mass for a given micro-state.
$\Gamma^{rec}_{\mathrm{aIEI^{\prime}E^{\prime}}}$, and
$\Gamma^{ion}_{\mathrm{aI^{\prime}E^{\prime}IE}}$ are the transition rate
coefficients between levels $I^{\prime}E^{\prime}$ and $IE$ due to
recombination or de-excitation, and ionization or excitation, respectively.
$q_{\mathrm{aI^{\prime}E^{\prime}IE}}$,
$P_{\mathrm{aI^{\prime}E^{\prime}IE}}$, and
$\hat{\tau_{\mathrm{aI^{\prime}E^{\prime}IE}}}$ are the ion charge, gas
pressure, and viscous tensor for a specific species. ${\vec{g}}$, ${\vec{E}}$,
and ${\vec{B}}$ are gravity acceleration, and electric and magnetic field,
respectively.
$\vec{R}_{\mathrm{aIE}}^{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}$ is the
momentum exchange where
$\mathrm{aIE}\neq\mathrm{a^{\prime}I^{\prime}E^{\prime}}$. For the collision
integrals, we consider neutral-neutral, ion-neutral including charge exchange,
or Maxwell molecular collisions, and Coulomb collisions between ions following
Bruno et al. (2010) (see Wargnier et al., 2022, for further details on the
collisions and cross-sections). The momentum exchange can then be expressed as
follows:

$\displaystyle\vec{R}_{\mathrm{aIE}}^{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}=m_{\mathrm{aIE}}\,n_{\alpha}\,\nu_{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}(\vec{u}_{\mathrm{a^{\prime}I^{\prime}E^{\prime}}}-\vec{u}_{\mathrm{aIEa}}),$
(15)

where $\nu_{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}$ is the collision
frequency. Note that, to guarantee the conservation of the total momentum,
$\vec{R}_{\mathrm{aIE}}^{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}=-\vec{R}_{\mathrm{a^{\prime}I^{\prime}E^{\prime}}}^{\mathrm{a^{\prime}I^{\prime}E^{\prime}aIE}}$.
$Q^{visc}$, $Q^{ir}_{\mathrm{aIE}}$,
$Q_{\mathrm{aIE}}^{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}$, and
$Q_{\mathrm{spitz}}$ are the viscous heating due to hyper-diffusion,
ionization and recombination, momentum exchange, and Spitzer thermal
conduction heating/cooling terms, respectively.
$Q_{\mathrm{aIE}}^{\mathrm{aIEa^{\prime}I^{\prime}E^{\prime}}}$ includes the
thermalization between the various fluids.

We ignore electron inertia and its time variation. So, in Ebysus, we consider
the magnetic induction equation:

$\displaystyle\frac{\partial{\vec{B}}}{\partial{t}}$ $\displaystyle=$
$\displaystyle-\nabla\times\vec{E}$ $\displaystyle=$
$\displaystyle\nabla\times\left(\vec{u_{e}}\times\vec{B}-\frac{\nabla
P_{e}}{n_{\mathrm{e}}q_{\mathrm{e}}}+\frac{\sum_{\mathrm{aIE}}\vec{R}_{\mathrm{e}}^{\mathrm{eaIE}}}{n_{\mathrm{e}}q_{\mathrm{e}}}\right)$

the second and third terms are the Bierman battery and the ohmic diffusion,
respectively, and we ignored the ionization recombination contribution to the
electric field. Since electrons move so fast and their mass is negligible we
assume quasi-neutrality:
$n_{\mathrm{e}}=\sum_{\mathrm{aIE}}n_{\mathrm{aIE}}\,Z_{\mathrm{aI}}$ where
$Z_{\mathrm{aI}}$ is the ionized state, and we neglect their inertia. Since
several ionized species are considered, the electron velocity, thanks to the
assumption of quasi-neutrality, reads as follows:

$\displaystyle\vec{u}_{e}=\left(\sum_{\mathrm{a\hat{I}E}}\frac{n_{\mathrm{a\hat{I}E}}q_{\mathrm{a\hat{I}E}}\vec{u}_{\mathrm{a\hat{I}E}}}{n_{\mathrm{e}}\,q_{\mathrm{e}}}\right)-\frac{\vec{J}}{q_{\mathrm{e}}\,n_{\mathrm{e}}},~{}$
(17)

and $\vec{J}=(\nabla\times\vec{B})/\mu_{0}$. We refer to the total ion
velocity as ${\vec{u}}_{c}=\sum_{\mathrm{a\hat{I}E}}{\vec{u}}_{\mathrm{aIE}}$
and for neutrals as
${\vec{u}}_{n}=\sum_{\mathrm{a0E}}{\vec{u}}_{\mathrm{a0E}}$.

Ebysus has inherited the numerical methods from Bifrost (Gudiksen et al.,
2011), and the first results of this code can be found in Martínez-Sykora et
al. (2020). In short, the numerical mesh is defined in a staggered Cartesian
box. The spatial derivatives and interpolation of the variables are 6th and
5th-order polynomials, respectively. To advance in time, we can select which
source terms are advanced explicitly or implicitly using either Lie or
Additive operator splittings. For the explicit part, we stepped forward in
time using the modified explicit third-order predictor-corrector Hyman method
(Hyman, 1979), or we can also choose a 3rd-order Runge-Kutta method. The
implicit part can be solved using different methods, some of which are
mentioned below and further details will be described in a dedicated
manuscript for the code. Finally, the numerical noise is suppressed using a
high-order artificial diffusion.

Due to the numerical stiffness, the ionization, recombination, excitation and
de-excitation, momentum exchange, and the ohmic terms have been added and are
solved by applying a Newton-Raphson implicit method, or the fifth-order
implicit Runge Kutta method known as Radau IIA method (Hairer & Wanner, 1999,
1996).

For this research, we neglect ionization and recombination, gravity, and the
heating or cooling sources from radiation, ohmic diffusion, and thermal
conduction. These assumptions allow us to “isolate” and limit the parameter
range of the considered dependencies of the ponderomotive force on Alfvén
waves.

###  5.2 Initial and boundary conditions

We investigate the effects of MFMS interactions and LTE vs NEQ ionization on
the ponderomotive force due to Alfvén waves in a parametric study for
chromospheric conditions. The density and temperature are taken from the
median of the histograms in Figure 7 assuming the region of interest from
panels d (LTE) or h (NEQ) in Figure 5 and listed in Table 1. We consider one
LTE case. For computational reasons, we had to select higher ion densities to
be able to run this LTE scenario for hydrogen and helium. This is because the
very different range of values of number density for the various fluids leads
to computational complications. The temperature is constant with $1.6\times
10^{4}$ K. Finally, the magnetic field follows a quarter of a sinusoidal
profile which varies within the values listed in the sixth column of Table 2
where the first number is at the bottom and the last is at the top. The range
for the Alfvén speed is listed in the sixth column. The values and range of
the magnetic field and Alfvén speed are similar to the parameter range within
the region of interest shown in Figure 8 except near the transition region
where the Alfvén speed variations with height can become much larger than the
range covered here. Because we assume constant densities, the variation of the
Alfvén speed in the box has the opposite sign to that of the magnetic field.
This is contrary to what would be expected on the Sun since the density
stratification leads to an increase in Alfvén speed with height, despite the
decrease in magnetic field strength expected from an expanding flux tube. This
should be considered when we refer to the magnetic field variation along the
loop as the “expansion” or “constriction” of field lines in loops. The
variation of the Alfv’en speed along the loop is expected to play a
significant role in ponderomotive force and is thus not fully considered here,
given our constant density assumption.

The numerical domain in all 1.5D numerical experiments is along the loop,
which for the MFMS models will be the $z$-axis, and covers the range $z=[0,3]$
Mm with a uniform grid of 1000 points, which is roughly the length where the
fractionation may occur in panel j) of Figure 5. The grid spacing has been
chosen to resolve spatially and temporally the Alfvén waves for any
simulation. Note that $z=0$ is not the photosphere as in the previous section
but the range $z=[0,3]$ Mm covers a subregion within the chromosphere. In all
of our simulations, waves are launched from z=0 Mm, and propagate towards
positive values of z. As the wave propagates it will either encounter an
increase of magnetic field strength (i.e., a constricting flux tube), or a
decrease of magnetic field strength (i.e., an expanding flux tube).

Table 1: Density number of the various fluids. From left to right: the species, and the number density for the neutral and ionized fluid in cm-3 Species | Neutrals (NEQ, LTE) | Ions (NEQ, LTE)
---|---|---
H |  $10^{9.8}$, $10^{11.1}$ |  $10^{10.7}$, $10^{8.5}$
He |  $10^{9.7}$, $10^{10.6}$ |  $10^{8.3}$, $10^{4.8}$
Fe |  $10^{4.6}$, $10^{5.25}$ |  $10^{5.6}$, $10^{6.4}$
Si |  $10^{4.7}$, $10^{5.7}$ |  $10^{5.4}$, $10^{6.2}$
Ca |  $10^{2.7}$, $10^{2.5}$ |  $10^{5.3}$, $10^{5.9}$

The boundaries are open at the “top” of the domain. At the “bottom,” we drive
an Alfvén wave along the component $x$ of the magnetic field with different
frequencies and amplitudes listed in the two rightmost columns of Table 2. We
put quotation marks around the words “top” and “bottom” because our numerical
experiments do not include gravity. The bottom to top direction is really the
direction of wave propagation. Note that the frequencies are inspired from the
observational constraints in Section 3 and the selected amplitudes are a
fraction of the derived velocity amplitudes from the non-thermal velocity
assuming that not all of the broadening is caused by Alfvén waves. Another
reason to not use large amplitudes is to avoid non-linearity effects in this
initial study. Columns 2-5 list which physical processes are included in each
of the simulations.

Table 2: List of numerical simulations. From left to right: list the name, physical processes included or not, magnetic field configuration, Alfvén speed range, and frequency and amplitude of the Alfvén wave driver. In the fifth column, the first value of the $B$ range is the bottom boundary, and the second is the top, so $[70,50]$ G is a loop that its magnetic field strength decreases with $z$, i.e., direction of wave propagation. Name | Ohm | Coll | Hall |  $B$ range (G) |  $v_{A}$ (km s-1) | Fr. (Hz) |  $B$ Amp. wave (G)
---|---|---|---|---|---|---|---
NH_B7050_F1_D1 | No | No | No | $[70,50]$ | $[290,210]$ | 1 | 0.1
NC_B7050_F1_D1 | No | No | Yes | $[70,50]$ | $[290,210]$ | 1 | 0.1
NC_B5070_F1_D1 | No | No | Yes | $[50,70]$ | $[210,290]$ | 1 | 0.1
C_B7050_F1_D1 | No | Yes | Yes | $[70,50]$ | $[290,210]$ | 1 | 0.1
C_B5070_F1_D1 | No | Yes | Yes | $[50,70]$ | $[210,290]$ | 1 | 0.1
A_B8040_F1_D1 | Yes | Yes | Yes | $[80,40]$ | $[330,170]$ | 1 | 0.1
A_B6555_F1_D1 | Yes | Yes | Yes | $[60,55]$ | $[270,220]$ | 1 | 0.1
A_B5070_F1_D1 | Yes | Yes | Yes | $[50,70]$ | $[210,290]$ | 1 | 0.1
A_B7050_F10_D1 | Yes | Yes | Yes | $[70,50]$ | $[290,210]$ | 10 | 0.1
A_B7050_F2_D1 | Yes | Yes | Yes | $[70,50]$ | $[290,210]$ | 2 | 0.1
A_B7050_F01_D1 | Yes | Yes | Yes | $[70,50]$ | $[290,210]$ | 0.2 | 0.1
A_B7050_F1_D5 | Yes | Yes | Yes | $[70,50]$ | $[290,210]$ | 1 | 0.5
A_B7050_F1_D1 | Yes | Yes | Yes | $[70,50]$ | $[290,210]$ | 1 | 0.1
A_B7050_F1_D05 | Yes | Yes | Yes | $[70,50]$ | $[290,210]$ | 1 | 0.05
A_B7050_F1_D01 | Yes | Yes | Yes | $[70,50]$ | $[290,210]$ | 1 | 0.01
LTE_B5070_F1_D1 | Yes | Yes | Yes | $[70,50]$ | $[225,160]$ | 1 | 0.1

###  5.3 Results of the MFMS models

In Section 2, we emphasize the challenging work of selecting the proper
assumptions to choose the appropriate ponderomotive force. In contrast, the
Ebysus code self-consistently solves the multi-fluid equations for various
species without considering any of those particular assumptions. As we will
show, the assumptions required for calculating the ponderomotive forces can be
very limiting, which can lead to an overestimation or misjudgement of their
impact on the plasma dynamics.

The ponderomotive force depends on the variation of the electric field along
the loop, while the electric field associated with Alfvén waves depends on the
wave energy. Therefore, first, we investigate the velocity and magnetic field
variation as a function of space. This will help to understand the variation
of the electric field with the distance along the loop (i.e., magnetic field
lines). These gradients are ultimately responsible for the acceleration along
the loop.

Let us focus first on a case without collisions, e.g., NC_B7050_F1_D1. Figure
9 shows various physical parameters as a function of space once the wave has
propagated through most of the domain. Since there are no collisions, neutral
fluids do not experience any variation in space and time as a result of an
Alfvén wave. However, it is interesting to see a rotation of the ion
velocities (second row) and magnetic field (bottom row in red) perpendicular
to the guided field, as shown with the two perpendicular components ($x$, and
$y$ in the first and second columns). The rotation is due to the ion coupling
between the various ionized fluids as reported by (Martínez-Sykora et al.,
2020) and not because of the Hall term. The equivalent simulation without the
Hall term (NH_B7050_F1_D1) also has this rotation (not shown here). The ion
coupling leads to a velocity drift between different ions (fourth row).
Furthermore, due to the expansion of the magnetic field, i.e., a magnetic
field strength decrease with $z$ and similarly with the Alfvén speed, there is
an upflow in the ion velocity (right column), which can produce a chemical
fractionation depending on the ionization fraction of each species. It is also
interesting to see the amplitude of the ion velocities (panels e-h) and the
perpendicular component of magnetic field (panels q-s) increase. The
wavelength decreases with $z$, i.e., with the decrease of the magnetic field
(panel t, red), This is expected from a propagating Alfvén wave in a medium
where the Alfvén speed decreases with $z$. Note that the Alfvén speed follows
the magnetic field strength since the density and temperature is constant.

Figure 9:  Spatial variation at t=5 s of velocity for all the neutral species
(top row), ion species (second row), velocity drift for each neutral fluid
with respect the total velocity of all fluids (third row), velocity drift for
individual ion-fluids with respect to the combined fluid of ions (fourth row).
In the bottom row, we show the magnetic field (red) and single fluid velocity,
i.e., the total velocity of all fluids (black). For all rows, we show the $x$
component, the $y$ component, the absolute value of the component
perpendicular to the guide field, and the $z$ component, from left to right
columns, respectively. This figure is for simulation NC_B7050_F1_D1, i.e., it
excludes collisions. In this case the field strength decreases in the
direction of wave propagation (e.g., for a wave propagating upwards in a loop
that expands with height).

This picture changes drastically when collisions are included, e.g.,
A_B7050_F1_D1 (Figure 10). Neutrals are being dragged by the ions (top row),
and the velocity patterns for each neutral fluid are shifted spatially because
of a phase shift, which depends on the collision frequency with other ions.
Still, neutrals are not completely coupled with ions, and a small velocity
drift is present (third row). As a result, the wave power, i.e., the velocity
perpendicular to the magnetic field, drops while it propagates along the
$z$-axis, i.e., dissipation occurs. Interestingly, vertical flows are much
larger than the collisionless scenario shown above. Like for the amplitudes,
the vertical velocities drop drastically with $z$ in contrast to the previous
case. This is because of the presence of collisions. Neutrals show relatively
large velocity drifts, which result from a vertical velocity variation with
$z$ due to collisional coupling, which does not occur in the collisionless
scenario shown previously. Consequently, the multi-fluid interactions lead to
a mix of wave phases between the various fluids, producing dissipation and a
large contribution to the ponderomotive force (as shown below) compared to the
collisionless case shown above.

Figure 10:  Same layout as Figure 9 for simulation A_B7050_F1_D1, i.e., with
collisions and for a decreasing field strength in the direction of the wave
propagation (i.e., for a wave propagating upwards in a loop that expands with
height).

We now study the role of the parameters listed in Table 2 on the variation of
the wave properties as a function of $z$ (i.e., the direction of wave
propagation, not necessarily height). To this end we calculate (and show in
Fig. 11), the gradient with $z$ of the peaks of the oscillations for the
various types of velocities: perpendicular to the field (first and third rows,
for neutrals and ions respectively), along the field (second and fourth rows,
for neutrals and ions respectively), and the total velocity of all fluids
combined (fifth row). As we can see, most of the cases we studied show a
negative gradient with $z$, i.e., a decrease with increasing $z$, except for a
few cases without collisions. An increase in the amplitude (first column) or
frequency (second column) of the driver leads to stronger gradients of
velocities with $z$ and greater drifts for neutrals (top panels) in both
perpendicular and longitudinal motions. The slopes of the various variables as
a function of the background magnetic field variation with $z$ provide a
fascinating result. For all the collisional cases, the slopes are negative,
regardless of whether the loop is expanding ($\Delta B<0$) or constricting
($\Delta B>0$) in the direction of wave propagation. The collisionless cases
(shown with dots) clearly have a different sign of the slopes (panels m and
n), i.e., it is positive for $\Delta B<0$ (expanding) and negative for $\Delta
B>0$ (constricting). For the collisional cases, the gradient with $z$ for all
variables increases with increasing magnetic field expansion (decreasing
$\Delta B$ for all negative $\Delta B$ cases). We will go into further details
below. Finally, one can appreciate that the important physical processes are
the presence of collisions and, in more minor roles, the ohmic diffusion and
the Hall term (right column). In the right column, we added the impact of the
LTE initial setup. Due to the differences in the populations of the various
species, the interactions between species and the gradients of the electric
field, hence the ponderomotive force, change dramatically compared to the NEQ
case. In LTE, the neutral density is higher than the ion density which results
in a stronger dissipation of the Alfvén wave and coupling between ion and
neutral fluids. Here again, we added dots for the non-collisional cases for
different magnetic field configurations for comparison.

Figure 11:  The median of the spatial gradient with $z$ of the wave-associated
velocities that are perpendicular to the field (first and third rows) and
longitudinal (second and fourth rows) for neutrals (top two rows) and ions
(third and fourth rows). We plot these as a function of amplitude and
frequency of the Alfvén wave driver, variation of the guide magnetic field,
and inclusion of various physical approaches, respectively from left to right.
The last row shows the spatial gradient with $z$ for wave-associated
velocities (for all fluids combined) that are perpendicular to the field
(black) and longitudinal (blue). We also show in the bottom row the spatial
variation of the perpendicular component of the magnetic field (red). The
collisionless cases have been added with dots in the two right columns.

In order to study the variation of the acceleration for different fluids as a
function of the physical parameters, first, we visualize the same variables as
in Figure 10 (which are as a function of space for a given time), but now as a
function of time at $z=0.18$ Mm. Note that there is no variation in the wave
amplitude with time (Figure 12). Therefore, the electric field has no
dependence on time. Consequently, the second term in the parenthesis of Eq. 1
is expected to be negligible. Neutral fluids, as mentioned above, experience a
different phase speed due to the differences in collisional frequencies. We
note that neutral fluids experience more significant differences in vertical
velocities than ions because the latter are magnetized and move with the
magnetic field. A similar effect can be seen for the gradients with $z$, as
shown in Figure 10. Again, one can appreciate that there is a large phase
mixing between the different types of fluids, which leads to the ponderomotive
force, as we show below.

Figure 12:  Temporal variation at $z=0.18$ Mm of velocity for all the neutral
species (top row), ion species (second row), velocity drift for each neutral
fluid with respect to the total velocity of all fluids (third row), velocity
drift for individual ion-fluids with respect to the combined ion-fluid (fourth
row). In the bottom row, we show the magnetic field (red) and single fluid
velocity, i.e., the total velocity of all fluids (black). For all rows, we
show the $x$ component, the $y$ component, the absolute value of the component
perpendicular to the guide field, and the $z$ component, from left to right
columns, respectively. This figure is for simulation A_B7050_F1_D1, i.e., with
collisions and a decreasing magnetic field strength with $z$.

We now calculate in Figure 13 the gradients with time of the local maximum
velocities in the $z$ direction (right column) of neutrals (top) and ions
(bottom) shown with red dashed lines, similarly as for Figure 11 but in time.
This figure shows the estimates of the acceleration (in the direction of wave
propagation) of the various fluids. The comparison of these accelerations, for
each physical process, initial plasma properties, or imposed driver, allows us
to further understand the action of the ponderomotive force associated with
Alfvén waves in the chromosphere for various fluids (Figure 13). The
acceleration increases with the amplitude and frequency of the driver. We
remind the reader that the selected frequencies are inspired by IRIS
observations (Section 3). It is very intriguing that the acceleration
increases with frequency whereas Eqs. 1-2 are proportional to $1/\omega^{2}$.
As shown in the third column, in the presence of collisions, there is always
an acceleration in the direction of the wave propagation, no matter if the
flux tube is expanding or constricting. There is some decrease in the
acceleration with an increase in the expansion of the magnetic field, i.e.,
values of decreasing $\Delta B$ for the cases where $\Delta B$ is negative. If
we assume that a loop expands with height, and assume a collisional case, the
plasma is accelerated in the same direction as the propagating waves. However,
in the collisionless cases, if the loop constricts ($\Delta B>0$) plasma is
accelerated in the direction of wave propagation, while, if the loop expands
($\Delta B<0$), plasma is accelerated in the direction opposite of the wave
propagation. If we assume that a loop expands with height, and assume no
collisions, this would mean that plasma is accelerated downward for upward
propagating waves, and downward for downward propagating waves. It appears
that, for the collisionless case, the sign disagrees with equations 2-2, but
agrees with the sign of equation 1. In any case, the presence of collisions
leads to a major qualitative and quantitative difference.

Finally, the large accelerations result from collisions which lead to a vital
phase mixing between the fluids, as shown in the right column in Figure 13.
Note also the large difference between the LTE and NEQ cases. In LTE, the
number density of neutrals is very large and collisions are high.
Consequently, waves are damped within a very short traveled distance and the
upflows and acceleration is highly localized closer to the source of the wave
($z\sim 0$).

We notice that for higher wave frequencies or amplitude, collisions dissipate
the wave faster. Therefore, the accelerations in those cases are highly
localized closer to the driver. Further away from the driver (which is at
$z=0$), the spatial gradients of the variables and acceleration are close to
zero. This is similar, as mentioned above, to the LTE case. In the latter case
this due to the strong coupling with neutrals.

Figure 13:  The median of the vertical acceleration at $z=0.18$ Mm for
neutrals top(row) and ions (bottom row) as a function of Alfvén wave amplitude
and frequency, the guided magnetic field variation, and physics are shown from
left to right.

In order to understand the cause of the vertical acceleration as a function of
the various parameters, we compute (and show in Fig. 14), as a function of
$z$, the electric field perpendicular to the magnetic field (top row),
$-dE^{2}_{\bot}/dz$ (again calculated based on peaks associated with the wave)
as it quantifies the ponderomotive acceleration (middle row), and the Alfvén
speed (bottom row). Note that in our simulations, we find a negligible value
for $dE^{2}/dt$, so we only consider the spatial gradient for $E^{2}$ in our
proxy for the ponderomotive force. The dependence of the acceleration on the
parameters shown in Figure 13 can be understood by investigating the panel
that shows $-dE^{2}_{\bot}/dz$ in the middle row of Figure 14:

Figure 14:  Electric field perpendicular to the magnetic field (top row),
$-dE^{2}_{\bot}/dz$ (middle row) and Alfvén speed (bottom row) as a function
of $z$-axis. We show this for a variety of simulations with varying amplitude,
and varying frequency of the Alfvén wave driver, the variation of the guide
magnetic field, and different physical assumptions, as shown, respectively,
from left to right.

  * •

The greater the amplitude of the Alfvén wave (panels a and b in Figure 13),
the larger the values in the estimated ponderomotive force
($-dE^{2}_{\bot}/dz$) (panel e of Figure 14).

  * •

The higher the frequency of the Alfvén wave (panels c and d in Figure 13), the
larger the values in the estimated ponderomotive force (panel f of Figure 14).
Note that the force drops faster with $z$, when the frequency is higher. This
is because of the strong dissipation.

  * •

Due to collisions, $E_{\bot}$ decreases with $z$ for any background field
configuration considered here ($\Delta B$). This results in a positive
ponderomotive force (i.e., in the direction of wave propagation) for
constricting and expanding loops. There is a decrease in the acceleration and
ponderomotive force when increasing the expansion of the field lines, i.e.,
with decreasing $\Delta B$. Still, the presence of collisions has a stronger
impact on the ponderomotive force and acceleration than the variation of the
field configuration and Alfvén speed (panel k).

  * •

The collisionless cases reveal the large role of collisions. Indeed, for an
expanding or constricting loop, the sign of the acceleration changes, being
positive for constricting loops ($\Delta B>0$), and negative for expanding
loops ($\Delta B<0$) (dots in panels e-h in Figure 13). This is in agreement
with the $E_{\bot}$ dependence with $z$, i.e., $E_{\bot}$ decreases in the
direction of the propagating wave (with $z$) for a constricting loop (red and
purple lines in panels d, h, and l in Figure 14) (since the Alfvén speed
increases) and as a result the estimated ponderomotive force is positive.
Whereas for an expanding loop, $E_{\bot}$ increases in the direction of the
propagating wave (brown lines panels d, h and l in Figure 14) (since the
Alfvén speed decreases) and, as a result, the estimated ponderomotive force is
negative. Note that the sign of this acceleration follows Eq. 1.

  * •

The initial setup of the density population (NEQ vs. LTE) greatly impacts the
electric field changes with $z$, hence the estimated ponderomotive force and
acceleration of the various fluids, as shown previously.

  * •

The motions in the perpendicular plane (and electric field) are very similar
between the NC_B7050_F1_D1 case and NH_B7050_F1_D1, i.e., with and without the
Hall term.

It is clear that collisions play a critical role in the ponderomotive force
within the parameter range considered here (Tables 1 and 2), which is inspired
by chromospheric values. This is shown in the right columns in Figures 11, 13
and 14). The sign of the ponderomotive acceleration agrees with that found in
Eq. 1. However, as shown in Figure 15, the ratio of the acceleration with
(-$dE^{2}_{\bot}/dz$) does not follow $1/\omega^{2}$, but follows
$\omega^{2}$. This dependence is, however, in agreement with the predictions
from damping of Alfvén waves from ion-neutral interactions (Haerendel, 1992;
De Pontieu & Haerendel, 1998). The predicted damping length $L$ for the wave
field is given in Eq. 8 and can be compared with the distance over which the
wave field $E_{\bot}^{2}$ decays in the top row of Fig. 14. The damping length
in our simulations is similar to what is predicted and shows only a
discrepancy of $\sim 2$.

Figure 15:  The ratio of the acceleration with $dE^{2}_{\bot}/dz$ for ionized
helium as a function of frequency. The blue line follows a dependence of
$\omega^{2}$ and matches what is expected from the ponderomotive force
associated with damping of Alfvén waves from ion-neutral damping (see text for
details).

Note that the setup of the simulations allows considering scenarios where the
Alfvén wave is driven from “above”, i.e., corona, or “below”, the photosphere.
Under those highly simplified scenarios, i.e., no density stratification nor
reflections and in the collisionless case, for a wave coming from the corona,
an expanding loop with height, i.e., a constricting loop in the direction of
the wave propagation, will produce an acceleration opposite to the wave
propagation direction, i.e., the FIP effect. In contrast, a constricting loop
with height (which may be quite rare, as it is not a condition thought to
commonly occur), i.e., an expanding loop in the direction of a wave coming
from the corona, will produce an acceleration in the direction of the wave,
which could conceivably lead to an inverse FIP effect. However, collisions
will lead to a ponderomotive force in the direction of the wave propagation
independently of the field configuration, producing accelerations in the
direction of the wave propagation. This could imply an inverse FIP effect for
a wave coming from the corona.

##  6 Discussions and Conclusions

One of the most generally accepted models to explain the FIP effect is based
on the effects of the ponderomotive force, using semi-empirical models (e.g.,
Laming, 2004). In this study, by combining IRIS observations with single and
multi-fluid models, we go beyond the limitations of those assumptions. We
investigate the possible role of the ponderomotive force on accelerating the
plasma by considering multi-fluid effects and the ionization under LTE and NEQ
conditions.

Being aware of complexity of the mixture of waves in the solar atmosphere and
LOS effects (Okamoto & De Pontieu, 2011; De Pontieu et al., 2004), we used
non-thermal velocities from a transition region line, i.e., Si IV from IRIS
observations to provide constraints on the properties of unresolved high-
frequency waves. Assuming that the non-thermal broadening can be ascribed to
waves, one can estimate their lowest frequencies and highest amplitudes. We
found that plage and network field regions experience the largest non-thermal
velocities.

Using the radiative-MHD numerical model described in Martínez-Sykora et al.
(2020), we investigate the properties of the regions of interest where low and
high FIP elements first ionize in the solar atmosphere for a plage region. Our
analysis reveals that, in the chromosphere, the largest perturbations of the
electric field component perpendicular to the magnetic field occur in areas of
enhanced network field or plage regions, associated with type II spicules and
large changes in the magnetic field topology. Chromospheric spicules occur
most frequently around the network field that is concentrated at the edges of
supergranular cells, i.e., where non-thermal broadening in the low TR is
enhanced (see also Pereira et al., 2014). In simulations of spicules, they are
associated with strong currents and drive Alfvén waves (Martínez-Sykora et
al., 2017), suggesting that, if the ponderomotive force indeed drives the FIP
effect, this could lead to chemical fractionation. Our results reveal that the
NEQ effects largely impact the region where the chemical fractionation may
occur within the chromosphere and needs to be considered when studying the FIP
effect. We note that the recent work by Wargnier et al. (2022) provides
another possible chemical fractionation mechanism that depends on the history
of the multi-fluid evolution for reconnection events.

In the last part of this manuscript, inspired by our observational and single
fluid analysis, we performed multi-fluid numerical models of Alfvén waves in a
non-uniform magnetic field configuration. One has to keep in mind that the
multi-fluid models used here are highly simplified: they do not include
density stratification or gravity, nor do they capture the chromospheric
dynamics, ionization/recombination, thermal conduction, or reflections of
waves. Further work is obviously needed to investigate further the various
effects that were ignored in this initial experiment.

However, in our results, the presence of different fluids (ionized and neutral
fluids for various species) leads to both collisional and electric coupling.
This, in turn, causes a (wave) phase offset between the various fluids and
results in damping of the Alfvén waves due to collisions while propagating
through the atmosphere. This coupling and damping leads to a positive
ponderomotive acceleration, i.e., in the direction of the propagating wave. If
the wave propagates upward, the ponderomotive force is upward, and for
downward propagating waves, the force is downward. The role of collisions is
much stronger than the variations of the magnetic field (and thus the Alfvén
speed) along the loop, which are meant to mimic a loop that expands or
constricts with height. Our results suggest that the presence of ion-neutral
damping in the chromosphere has a large potential to dominate any other type
of ponderomotive forcing, and may invalidate previous approaches that are
based on other forms of the ponderomotive force. For example, the direction of
the ponderomotive force in our simulations is different from what has been
proposed in the recent literature (see below). Our results suggest that if the
ponderomotive force associated with Alfvén waves were the dominant mechanism
driving the FIP effect, the Alfvén waves should propagate from the
chromosphere to corona. This is contrast to what has been assumed so far in
current theoretical models of the ponderomotive force (Laming, 2017). However,
it brings that theory in better alignment with the current observational
evidence that indicates predominant upward propagation from the chromosphere
into the corona (Okamoto & De Pontieu, 2011) and a dominance of upward
propagating waves within the coronal volume at lower frequencies (Tomczyk et
al., 2007).

We find that in our models, the ponderomotive force increases with increasing
wave amplitude or frequency of the Alfvén wave. We also find an increase in
ponderomotive force in an environment in which the magnetic field strength
increases in the direction of the wave propagation. We run cases assuming
different approaches to ionizations with both LTE populations and NEQ
populations. The two scenarios are very different, leading to different
ponderomotive accelerations and damping mechanisms. The LTE case tends to damp
the waves faster due to the larger presence of neutrals.

Our results lead us to speculate about what this means for sunspots, which
have much lower densities than other solar regions or phenomena. We speculate
that this would lead to a reduced collisional coupling, so that our
collisionless case might apply. In that case, this could potentially explain
the inverse FIP effect that is sometimes seen in sunspots (Doschek & Warren,
2016). However, this speculation needs further investigation within the
parameter range of a sunspot stratification.

We finish with the caveat that this first exploration of the relative motions
and accelerations of ions with different first ionization potential is highly
simplified and does not yet include various effects, including stratification.
The latter has the potential to lead to larger spatial variations in the
Alfvén speed than we considered here, which were only based on the varying
guide field (not the density). At a minimum, this would mean that the
conditions that we have called “expanding” or “constricting” loops (purely
based on guide field changes) may well change meaning when stratification is
included, given the expected change in sign of $dv_{A}/ds$. This is important
to note since the ponderomotive force also depends on changes in $v_{A}$. In
any case, large changes in Alfvén speed have the potential to significantly
change the ponderomotive force, including through reflections from the
transition region. Future numerical multi-fluid models that include more
complexity are needed to move beyond the simplified analytical approaches of
the current state-of-the-art FIP models and better understand how the
ponderomotive force acts in the low solar atmosphere and can lead to
fractionation.

It is clear from our results however that the presence of collisions in the
partially ionized chromosphere introduces a ponderomotive force that, for the
upper chromospheric thermodynamic and magnetic field conditions we have
considered here, may well dominate other types of ponderomotive force.

We gratefully acknowledge support by NASA grants 80NSSC18K1285, 80NSSC20K1272,
80NSSC21K0737, 80NSSC21K1684, and contract NNG09FA40C (IRIS). 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 Pleiades through the
computing project s1061, s2601, and s8305. This research is also supported by
the Research Council of Norway through its Centres of Excellence scheme,
project number 262622, and through grants of computing time from the Programme
for Supercomputing. JMS is very grateful to his family, and especially to his
wife, for their support over the many years that allowed him to progress in
his career, develop Ebysus with the incredible help from the Ebysus team and
culminate with this work.

## References

  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Avrett (2007) Avrett, E. H. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 368, The Physics of Chromospheric Plasmas, ed. P. Heinzel, I. Dorotovič, & R. J. Rutten, 81
  * Bale et al. (2021) Bale, S. D., Horbury, T. S., Velli, M., et al. 2021, ApJ, 923, 174, doi: 10.3847/1538-4357/ac2d8c
  * Bale et al. (2022) Bale, S. D., Drake, J. F., McManus, M. D., et al. 2022, arXiv e-prints, arXiv:2208.07932.  https://arxiv.org/abs/2208.07932
  * Bruno et al. (2010) Bruno, D., Catalfamo, C., Capitelli, M., et al. 2010, Physics of Plasmas, 17, 112315, doi: 10.1063/1.3495980
  * Carlsson & Stein (1992) Carlsson, M., & Stein, R. F. 1992, ApJ, 397, L59, doi: 10.1086/186544
  * Dahlburg et al. (2016) Dahlburg, R. B., Laming, J. M., Taylor, B. D., & Obenschain, K. 2016, ApJ, 831, 160, doi: 10.3847/0004-637X/831/2/160
  * De Pontieu et al. (2004) De Pontieu, B., Erdélyi, R., & James, S. P. 2004, Nature, 430, 536, doi: 10.1038/nature02749
  * de Pontieu & Haerendel (1998) de Pontieu, B., & Haerendel, G. 1998, A&A, 338, 729
  * De Pontieu & Haerendel (1998) De Pontieu, B., & Haerendel, G. 1998, A&A, 338, 729
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2007, Science, 318, 1574, doi: 10.1126/science.1151747
  * De Pontieu et al. (2014) De Pontieu, B., Rouppe van der Voort, L., McIntosh, S. W., et al. 2014, Science, 346, D315, doi: 10.1126/science.1255732
  * Doschek & Warren (2016) Doschek, G. A., & Warren, H. P. 2016, ApJ, 825, 36, doi: 10.3847/0004-637X/825/1/36
  * Fargette et al. (2021) Fargette, N., Lavraud, B., Rouillard, A. P., et al. 2021, ApJ, 919, 96, doi: 10.3847/1538-4357/ac1112
  * Feldman (1992) Feldman, U. 1992, Phys. Scr, 46, 202, doi: 10.1088/0031-8949/46/3/002
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7, doi: 10.1007/s11214-015-0211-6
  * Geiss et al. (1995) Geiss, J., Gloeckler, G., & von Steiger, R. 1995, Space Sci. Rev., 72, 49, doi: 10.1007/BF00768753
  * Golding et al. (2016) Golding, T. P., Leenaarts, J., & Carlsson, M. 2016, ApJ, 817, 125, doi: 10.3847/0004-637X/817/2/125
  * Gudiksen et al. (2011) Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154+, doi: 10.1051/0004-6361/201116520
  * Haerendel (1992) Haerendel, G. 1992, Nature, 360, 241, doi: 10.1038/360241a0
  * Hairer & Wanner (1996) Hairer, E., & Wanner, G. 1996, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Vol. 14, doi: 10.1007/978-3-662-09947-6
  * Hairer & Wanner (1999) —. 1999, Journal of Computational and Applied Mathematics, 111, 93, doi: https://doi.org/10.1016/S0377-0427(99)00134-X
  * Hyman (1979) Hyman, J. M. 1979, in Advances in Computer Methods for Partial Differential Equations - III, 313–321
  * Kasper et al. (2012) Kasper, J. C., Stevens, M. L., Korreck, K. E., et al. 2012, ApJ, 745, 162, doi: 10.1088/0004-637X/745/2/162
  * Khomenko & Cally (2011) Khomenko, E., & Cally, P. S. 2011, Journal of Physics: Conference Series, 271, 012042, doi: 10.1088/1742-6596/271/1/012042
  * Khomenko et al. (2021) Khomenko, E., Collados, M., Vitas, N., & González-Morales, P. A. 2021, Philosophical Transactions of the Royal Society of London Series A, 379, 20200176, doi: 10.1098/rsta.2020.0176
  * Laming (2004) Laming, J. M. 2004, ApJ, 614, 1063, doi: 10.1086/423780
  * Laming (2015) —. 2015, Living Reviews in Solar Physics, 12, 2, doi: 10.1007/lrsp-2015-2
  * Laming (2017) —. 2017, ApJ, 844, 153, doi: 10.3847/1538-4357/aa7cf1
  * Landi & Testa (2015) Landi, E., & Testa, P. 2015, ApJ, 800, 110, doi: 10.1088/0004-637X/800/2/110
  * Leenaarts et al. (2007) Leenaarts, J., Carlsson, M., Hansteen, V., & Rutten, R. J. 2007, A&A, 473, 625, doi: 10.1051/0004-6361:20078161
  * Lundin & Guglielmi (2006) Lundin, R., & Guglielmi, A. 2006, Space Sci. Rev., 127, 1, doi: 10.1007/s11214-006-8314-8
  * Lundin & Hultqvist (1989) Lundin, R., & Hultqvist, B. 1989, J. Geophys. Res., 94, 6665, doi: 10.1029/JA094iA06p06665
  * Martínez-Gómez et al. (2015) Martínez-Gómez, D., Soler, R., & Terradas, J. 2015, A&A, 578, A104, doi: 10.1051/0004-6361/201525785
  * Martínez-Sykora et al. (2017) Martínez-Sykora, J., De Pontieu, B., Hansteen, V. H., et al. 2017, Science, 356, 1269, doi: 10.1126/science.aah5412
  * Martínez-Sykora et al. (2020) Martínez-Sykora, J., Leenaarts, J., De Pontieu, B., et al. 2020, ApJ, 889, 95, doi: 10.3847/1538-4357/ab643f
  * McIntosh & De Pontieu (2009) McIntosh, S. W., & De Pontieu, B. 2009, ApJ, 707, 524, doi: 10.1088/0004-637X/707/1/524
  * McIntosh et al. (2011) McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477, doi: 10.1038/nature10235
  * Meyer (1985) Meyer, J.-P. 1985, ApJS, 57, 151, doi: 10.1086/191000
  * Miller (1958) Miller, M. A. 1958, Russia
  * Müller et al. (2020) Müller, D., St. Cyr, O. C., Zouganelis, I., et al. 2020, A&A, 642, A1, doi: 10.1051/0004-6361/202038467
  * Nóbrega-Siverio et al. (2020) Nóbrega-Siverio, D., Moreno-Insertis, F., Martínez-Sykora, J., Carlsson, M., & Szydlarski, M. 2020, A&A, 633, A66, doi: 10.1051/0004-6361/201936944
  * Okamoto & De Pontieu (2011) Okamoto, T. J., & De Pontieu, B. 2011, ApJ, 736, L24, doi: 10.1088/2041-8205/736/2/L24
  * Pereira et al. (2014) Pereira, T. M. D., De Pontieu, B., Carlsson, M., et al. 2014, ApJ, 792, L15, doi: 10.1088/2041-8205/792/1/L15
  * Przybylski et al. (2022) Przybylski, D., Cameron, R., Solanki, S. K., et al. 2022, A&A, 664, A91, doi: 10.1051/0004-6361/202141230
  * Rempel & Przybylski (2021) Rempel, M., & Przybylski, D. 2021, ApJ, 923, 79, doi: 10.3847/1538-4357/ac2c6d
  * Song & Vasyliūnas (2011) Song, P., & Vasyliūnas, V. M. 2011, Journal of Geophysical Research (Space Physics), 116, 9104, doi: 10.1029/2011JA016679
  * Testa (2010) Testa, P. 2010, Space Sci. Rev., 157, 37, doi: 10.1007/s11214-010-9714-3
  * Testa et al. (2016) Testa, P., De Pontieu, B., & Hansteen, V. 2016, ApJ, 827, 99, doi: 10.3847/0004-637X/827/2/99
  * Testa et al. (2015) Testa, P., Saar, S. H., & Drake, J. J. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140259, doi: 10.1098/rsta.2014.0259
  * Tomczyk et al. (2007) Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192, doi: 10.1126/science.1143304
  * Wargnier et al. (2022) Wargnier, Q. M., Martínez-Sykora, J., Hansteen, V. H., & De Pontieu, B. 2022, ApJ, 933, 205, doi: 10.3847/1538-4357/ac6e62
  * Wargnier et al. (2022) Wargnier, Q. M., Martinez-Sykora, J., Hansteen, V. H., & De Pontieu, B. 2022, Multi-Fluid Simulations of Upper Chromospheric Magnetic Reconnection with Helium-Hydrogen mixture, arXiv, doi: 10.48550/ARXIV.2211.02157
  * Whitelam et al. (2002) Whitelam, S., Ashbourn, J. M. A., Bingham, R., Shukla, P. K., & Spicer, D. S. 2002, Sol. Phys., 211, 199, doi: 10.1023/A:1022408206824
  * Woolley et al. (2021) Woolley, T., Matteini, L., McManus, M. D., et al. 2021, MNRAS, 508, 236, doi: 10.1093/mnras/stab2281