# Preferential acceleration of heavy ions in magnetic reconnection: Hybrid-
kinetic simulations with electron inertia

Neeraj Jain  Zentrum für Astronomie und Astrophysik, Technische Universität
Berlin, Hardenbergstr. 36, D-10623, Berlin, Germany  Jörg Büchner  Zentrum für
Astronomie und Astrophysik, Technische Universität Berlin, Hardenbergstr. 36,
D-10623, Berlin, Germany  Miroslav Bárta  Astronomical Institute of the
Academy of Sciences of the Czech Republic, Fric̆ova 298, Ondr̆ejov, 251 65,
Czech Republic  Radoslav Buc̆ík  Southwest Research Institute, 6220 Culebra
Road, San Antonio, TX 78238, USA

###### Abstract

Solar energetic particles (SEPs) in the energy range 10s KeV/nucleon - 100s
MeV/nucleon originate from Sun. Their high flux near Earth may damage the
space borne electronics and generate secondary radiations harmful for the life
on Earth and thus understanding their energization on Sun is important for
space weather prediction. Impulsive (or 3He-rich) SEP events are associated
with the acceleration of charge particles in solar flares by magnetic
reconnection and related processes. The preferential acceleration of heavy
ions and the extra-ordinary abundance enhancement of 3He in the impulsive SEP
events are not understood yet. In this paper, we study ion acceleration in
magnetic reconnection by two dimensional hybrid-kinetic plasma simulations
(kinetic ions and inertial electron fluid). All the ions species are treated
self-consistently in our simulations. We find that heavy ions are
preferentially accelerated to energies many times larger than their initial
thermal energies by a variety of acceleration mechanisms operating in
reconnection. Most efficient acceleration takes place in the flux pileup
regions of magnetic reconnection. Heavy ions with sufficiently small values of
charge to mass ratio ($Q/M$) can be accelerated by pickup mechanism in outflow
regions even before any magnetic flux is piled up. The energy spectra of heavy
ions develop a shoulder like region, a non-thermal feature, as a result of the
acceleration. The spectral index of the power law fit to the shoulder region
of the spectra varies approximately as $(Q/M)^{-0.64}$. Abundance enhancement
factor, defined as number of particles above a threshold energy normalized to
total number of particles, scales as $(Q/M)^{-\alpha}$ where $\alpha$
increases with the energy threshold. We discuss our simulation results in the
light of the SEP observations.

magnetic reconnection, heavy ion acceleration

††journal: The Astrophysical Journal

##  1 Introduction

Non-thermal acceleration of charged particles is a widespread process in space
and astrophysical plasmas ranging from planetary magnetospheres to clusters of
galaxies. The accelerated particles, mostly ions, are detected near and on
Earth in a wide energy range , $10^{4}-10^{20}$ eV, either directly by
particle detectors on satellites and high altitude balloons or indirectly by
detecting the secondary particles and electromagnetic radiation produced by
their interaction with Earth’s atmosphere. These particles are generally
categorized as solar energetic particles (SEPs), galactic cosmic rays and
extra-galactic cosmic rays based on their source regions in our star the Sun,
our galaxy Milky Way and beyond our galaxy, respectively. Solar energetic
particles (SEPs) in typical energy range of 10s KeV/nucleon - 100s MeV/nucleon
are of particular interest because their flux near Earth is sufficiently high
to have implications for space weather phenomena.

SEP events can be classified into impulsive (short duration $\leq$ 1 day, less
intense — relatively smaller particle fluxes with typical energies $<10$
MeV/nucleon and numerous — about 1000 per year) and gradual (long duration
$\sim$ several days, orders of magnitude more intense — large particle fluxes
with energies $>$ 10 MeV/nucleon, less frequent $\sim$ 10 per year ) events
(Reames, 1999, 2021a). Gradual events are attributed to the acceleration of
charged particles by Coronal Mass Ejection (CME) driven shocks while impulsive
events to the acceleration in solar flares by magnetic reconnection and
associated processes (Reames, 2013; Desai & Giacalone, 2016; Bučík, 2020).
Impulsive events are electron rich and associated with type-III radio bursts
while gradual events are proton rich and associated with type-II radio bursts
(Cane et al., 1986; Reames, 2013). Impulsive events show abundance
enhancements of heavier ions relative to their abundances in the solar corona
at energies well above their average thermal energy ($\sim$ 100 eV) (Reames,
2021b). The enhancement factor, defined as the ratio of the relative (to a
reference element, usually Oxygen) abundances of an element X in impulsive SEP
events and solar corona, exhibits power law dependence on ion’s charge ($Q$)
to mass ($M$) ratio as $(Q/M)^{-\alpha}$. Event-to-event variations in the
value of $\alpha$ have been found (Reames et al., 2014a) with mean values
$\alpha=3.26$ (at 0.375 MeV/nucleon) (Mason et al., 2004), $\alpha=3.64\pm
0.15$ (at 3-10 MeV/nucleon) (Reames et al., 2014b), and $\alpha=3.53$ (at
160-226 KeV/nucleon) and $\alpha=3.31$ (at 320-453 KeV/nucleon) (Bučík et al.,
2021). Enhancement factor for 3He isotope, however, does not obey the power
law and can have very large values (up to $10^{4}$) compared to those
calculated using the power law (Kocharov & Kocharov, 1984; Mason, 2007). For
this reason, impulsive SEP events are also called 3He-rich events. The
enhancement of heavy ions is, however, not correlated with the extra-ordinary
enhancement of 3He (Mason et al., 1986; Reames, 1999). Abundance enhancement
of 3He is sometimes observed in gradual events as well (Desai et al., 2016;
Buc̆ík et al., 2023). These are associated either with the acceleration of the
remnant material in the interplanetary space from earlier flares by CME-driven
shocks or with the simultaneous activity in the corona.

The preferential acceleration of heavy ions with small value of $Q/M$ from
their thermal energies to SEP energies and the extra-ordinary abundance
enhancement of 3He are not understood yet. Since these enhancements are un-
correlated (Mason et al., 1986; Reames, 1999), attempts have been made to
explain the preferential acceleration of heavy ions and 3He by separate
mechanisms. Models of the preferential acceleration of heavy ions consider
resonant interaction of the ions with the waves in Alfvénic plasma turbulence
(Eichler, 2014; Kumar et al., 2017; Fu et al., 2020; Shi et al., 2022; Miller,
1998). Turbulent energy decays with wave frequency, i.e., higher frequency
waves have lower energy. Heavy ions with lower values of $Q/M$ have lower
values of cyclotron frequency ($\propto Q/M$) and thus would resonate with
lower frequency but higher power waves in turbulence, favoring preferential
acceleration of heavy ions. The mechanisms proposed for the preferential
acceleration of 3He consider absorption of some wave energy by cyclotron
resonance of 3He with the wave (Fisk, 1978; Temerin & Roth, 1992; Liu et al.,
2006). These waves can be produced by electron beams, e.g., by electrons
streaming along open magnetic field lines in solar corona, or via coupling
with low-frequency Alfvén waves. In fact, efficient acceleration of 3He by
ion-cyclotron resonance has been observed in nuclear fusion devices, which has
implications for space plasmas as well (Kazakov et al., 2017).

Magnetic reconnection can also accelerate electrons and ions in solar flares
(Zhou et al., 2015, 2016; Bárta et al., 2011b, a). It has also been considered
to explain the preferential acceleration of heavy ions. In 2.5-D
magnetohydrodynamic (MHD) and test particle simulations (Kramoliš et al.,
2022), heavy ions were found to be preferentially accelerated by first order
Fermi process in cascading plasmoids generated by spontaneous magnetic
reconnection in a meso-scale current sheet. The ion energy spectra and
abundance enhancement factors exhibit power-law profiles. The index of the
power law for the abundance enhancement factor, however, was not in agreement
with the observations. The authors suspected that the disagreement could be
due to the limitation of the MHD model which lacks the kinetic physics
essential for magnetic reconnection in collisionless plasmas. In 2-D fully
kinetic simulations of magnetic reconnection, ions entering the reconnection
exhaust can behave like pickup particles and gets accelerated if their $Q/M$
is below a threshold value, and thus preferential acceleration of heavy ions.
An explanation for the power-law dependence of the enhancements on $Q/M$ was
proposed based on the pickup mechanism (Drake et al., 2009; Knizhnik et al.,
2011). These simulations, carried out for only for one ion species (4He2+),
did not verify the power-law dependence. It is, therefore, not clear if
acceleration of heavy ions by pickup mechanism can provide the observed
scaling of the abundance enhancement with $Q/M$.

Magnetic reconnection in collisionless plasmas is a multi-scale process which
occurs at electron kinetic scales and then couples to ion and even larger
macro-scales. An ideal simulation of magnetic reconnection requires kinetic
treatment of electrons and ions and covering scales from electron kinetic
scales all the way up to large macro-scales well above the ion kinetic scales.
Such fully kinetic simulations, e.g., using Particle-in-Cell (PIC) method, of
magnetic reconnection in electron-proton plasma are computationally very
demanding. Self-consistent inclusion of heavier ion species make the
simulations even more demanding because now one has to use larger simulation
box to accommodate the larger gyro-radii of the heavier species while at the
same time resolve the electron scales. Computationally less demanding hybrid-
kinetic simulations, which treats ions kinetically and electron as an inertial
fluid, can be employed at the cost of electron kinetic physics which is
acceptable as far as the study of the ion acceleration is the primary
objective. Treatment of electrons as an inertial fluid allows to include the
physics at electron inertial scale but relaxes the numerical requirement of
resolving the Debye length in PIC method and therefore allowing larger
simulation domains for the same number of grid points.

In this paper, we employ a hybrid-kinetic plasma model (kinetic ions and
inertial electron fluid) to simulate magnetic reconnection in a 2-D plane with
the objective of studying ion acceleration in reconnection. All the ion
species are treated self-consistently in our simulations. We use the term
“heavy ion” to mean any ion species whose mass $M$ is larger than the proton’s
mass. The paper is organized as follows. Section 2 presents the simulation
setup. Simulation results are presented in section 3 and discussed and
concluded in section 4.

##  2 Simulation setup

We carry out 2-D hybrid-kinetic simulations of magnetic reconnection with
electron inertia using the hybrid-PIC code CHIEF (Code Hybrid with Inertial
Electron Fluid) which is a 3-D code parallelized based on Message Passing
Interface (MPI) for high performance computing (Muñoz et al., 2018; Jain et
al., 2022; Muñoz et al., 2023). In the hybrid code CHIEF, electrons are
treated as an inertial and isothermal fluid whose equations are coupled to
Maxwell’s equation to obtain the electric and magnetic fields. All the ions
(protons and heavier ions) , on the other hand, are treated self-consistently
as kinetic species and their equations of motion in electric and magnetic
fields are solved using the PIC method. The code CHIEF treats the inertial
effects of the electron fluid without any of the approximations used by other
electron-inertial hybrid-kinetic code (Jain et al., 2022). The details of the
hybrid-kinetic model used in CHIEF, its numerical implementation and
parallelization are discussed in our other publications (Muñoz et al., 2018;
Jain et al., 2022; Muñoz et al., 2023).

We initialize the 2-D simulations with two Harris current sheets,
$\mathbf{B}=\hat{z}B_{z0}[\tanh\\{(y+L_{y}/4)/L\\}-\tanh\\{(y-L_{y}/4)/L\\}-1]+\hat{x}B_{x0}$,
in a y-z plane and a guide magnetic field $B_{x0}=0.2\,B_{z0}$ perpendicular
to the plane. The half-thickness $L=d_{p}$ of the current sheets is taken to
be equal to a proton inertial length ($d_{p}$) and $L_{y}$ is the length of
the simulation domain along y-direction. A small initial perturbation is added
to the Harris equilibrium to form an X-point and an O-point in the current
sheets centered at $y=-L_{y}/4$ and $y=L_{y}/4$, respectively. The plasma of
Harris current sheets consists of quasi-neutral populations of proton
particles and electron fluid with Harris equilibrium density profile
$n_{e}=n_{p}=n_{0}[\text{sech}^{2}\\{(y+L_{y}/4)/L\\}+\text{sech}^{2}\\{(y-L_{y}/4)/L\\}]$.
Harris sheets are embedded in a uniform background plasma of density
$0.2\,n_{0}$ and therefore peak density at their centers is $1.2\,n_{0}$. The
background plasma consists of electron fluid and particle populations of heavy
ions (only one species) and protons with densities $n_{be}=0.2\,n_{0}$,
$n_{bi}=0.01\,n_{0}$ (5% of the background plasma density) and
$n_{bp}=n_{be}-Z\,n_{bi}$, respectively, where $Z$ is the charge state of
heavy ions with charge of heavy ion species given by $Q=Z\,e$. The initial
velocity distribution of the protons and heavy ions is Maxwellian. The initial
temperatures of electrons, protons and heavy ions are the same, viz.,
$T_{e}=T_{p}=T_{i}=0.25\,m_{p}V_{Ap}^{2}$ where $m_{p}$ is proton mass and
$V_{Ap}=B_{z0}/\sqrt{\mu_{0}n_{0}m_{p}}$ is the proton Alfvén velocity based
on $B_{z0}$ and $n_{0}$.

Each of our simulations consists of only one species of heavy ion. We consider
the following four different species of heavy ions to be included (one at a
time) in our simulations : 4He2+, 3He2+, 16O7+ and 56Fe14+. We take proton to
electron mass ratio as $m_{p}/m_{e}=25$. The simulation domain $L_{y}\times
L_{z}$ is $51.2\,d_{p}\times 102.4\,d_{p}$ resolved by a grid spacing of
$0.1\,d_{p}$ in each direction. The time step is $\Delta
t=0.0025\,\omega_{cp}^{-1}$, where $\omega_{cp}=eB_{z0}/m_{p}$ is the proton
cyclotron frequency. We take 500 and 200 particles per cell for protons and
heavy ions, respectively. Plasma resistivity is taken to be zero and electron
inertia allows the magnetic reconnection. Boundary conditions are periodic.

##  3 Simulation results

Fig. 1 shows evolution of out-of-plane current density $J_{x}$ from initial to
late phase. In the initial phase ($\omega_{cp}t=19.6$), an X-point forms in
the lower current sheet ($y<0$) while an O-point forms in the upper current
sheet ($y>0$), as per the initialized perturbation. By $\omega_{cp}t=39.19$,
the lower (also the upper) current sheet spontaneously develop magnetic
islands as a result of magnetic reconnection at multiple sites, even though
the initialized perturbation was chosen to initiate reconnection at a single
site in the simulation domain. The magnetic islands on each of the current
sheet grow in size with time by merging and/or pushing among themselves and
simultaneously develop turbulence inside them. At $\omega_{cp}t=58.79$, the
magnetic islands of the lower and upper current sheet grow to big enough size
that the particles accelerated in the upper (lower) current sheet may cross
over to the lower (upper) current sheet.

Our objective here is to study the ion acceleration without the influence of
the periodic boundaries. Since the upper current sheet is affected by the
periodic boundary conditions along z-direction from the beginning (the initial
X-points form at the z-boundaries), we focus on the lower current sheet for
our studies. The lower current sheet is also likely to be affected by periodic
boundaries along z-direction after a time $\sim 50\omega_{cp}^{-1}$, Alfvén
waves take to cross half the simulation domain. In order to avoid particles
crossing from the upper half to the lower half of the simulation domain and
the influence of periodic boundaries, we limit our analysis of results only up
to the time $\omega_{cp}t=48.99$ and in the spatial region $-23\leq
y/d_{p}\leq-3$.

$\omega_{cp}t=19.6$$\omega_{cp}t=39.19$
$\omega_{cp}t=48.99$$\omega_{cp}t=58.79$

Figure 1: Out-of-plane current density $J_{x}/(n_{0}eV_{Ap})$ at four
different times.

Figure 2a shows the evolution of fractional change in average kinetic energy
from its initial value for different ion species. All the heavy ions gain
energy first at slow rates up to $\omega_{cp}t\approx 30$ after which they
gain energy at much faster rates. Note that magnetic islands have
significantly developed by $\omega_{cp}t=30$. This can be seen in Fig. 2b
which shows that $\langle B_{y}^{2}(y=-12.5\,d_{p},z)\rangle_{z}/B_{z0}^{2}$
(average of the magnetic energy in the normal component of the magnetic field
along the central line of the lower current sheet — a proxy for the magnetic
island development) begins to grow at $\omega_{cp}t=20$ and, by
$\omega_{cp}t=30$, has grown to $\sim$ 5% of the asymptotic magnetic energy in
the anti-parallel magnetic field. The development of magnetic islands,
therefore, seem to be linked with the efficient energization of heavy ions.
Proton energization, on the other hand, does not seem to be affected
significantly by the formation of magnetic islands as it does not enhance
significantly after $\omega_{cp}t\approx 30$.

(a)(b)

Figure 2: Time development of the fractional change in average kinetic energy
per particle, $\langle K.E.(t)\rangle_{ptcl}$, from its initial value (a) and
average of the square of the normal component of the magnetic field $B_{y}$
along the central line of the lower current sheet (b) for different ion
species.

Heavier ions (16O7+ and 56O14+ in Fig. 2) get energized significantly in
comparison to the lighter Helium ions even during early phase of acceleration
($0<\omega_{cp}t<30$). This suggests that acceleration mechanisms, different
from the mechanisms after $\omega_{cp}t=30$, operate in the early phase only
on the heavier ions, implying a threshold for the energization based on the
heaviness of the ions. In order to understand the threshold behavior and the
acceleration mechanisms before and after $\omega_{cp}t=30$ (the time by which
$\sim$5% of the magnetic energy of the asymptotic anti-parallel magnetic field
has contributed to the development of magnetic islands), we show the locations
of the energized particles in the reconnection region at $\omega_{cp}t$=19.6
and 48.99 in two simulations — one with 4He2+ (Fig. 3) and the other with
16O7+ (Fig. 4). At $\omega_{cp}t=19.6$, 16O7+ ions are more energized than
4He2+ ions and the locations of the energized 16O7+ ions are in the outflow
regions. The energization is, however, not uniformly distributed in the
outflow regions: the energized particles are located near the upper (lower)
separatrix in the left (right) outflow region. This observations combined with
the threshold behavior of energization suggests the pick-up mechanism of the
energization in which non-adiabatic ions are energized as they cross the
reconnection separatrices (Drake et al., 2009). A threshold condition,
$m_{i}/Z_{i}m_{p}>5\sqrt{2\beta_{p}}/\pi$, for non-adibaticity of ions
crossing the sepratrix in guide field reconnection was obtained by Drake et
al. (Drake et al., 2009), where $m_{i}$ and $Z_{i}=Q_{i}/e$ are the mass and
charge state of ions and $\beta_{p}$ is the proton plasma beta based on
asymptotic value of the anti-parallel magnetic field. Although this condition,
which becomes $m_{i}/Z_{i}m_{p}>1.12$ for our simulation parameters, is
satisfied by all the species of heavy ions we have considered, the 4He2+ and
3He2+ ions with $m_{i}/Z_{i}m_{p}$ = 2 and 1.5 were not energized during
$\omega_{cp}t<30$.

Note that the threshold condition was obtained assuming that the ions cross
the separatrix with a velocity equal to $0.1v_{Ap}$ which is almost the upper
bound on the inflow velocity in a fully developed steady state magnetic
reconnection (Liu et al., 2017). In our simulations, magnetic reconnection at
$\omega_{cp}t=19.6$ is still developing and the ion inflow velocity has not
yet reached its maximum value. Using $|u_{iy}|\approx 0.05\,v_{Ap}$, the value
of inflow velocity at $\omega_{cp}t=19.6$ in our simulations, the threshold
condition becomes $m_{i}/Z_{i}m_{p}>2.24$ which allows non-adiabatic behavior
for 56Fe14+ ($m_{i}/Z_{i}m_{p}=4$) and only marginally for 16O7+
($m_{i}/Z_{i}m_{p}=2.28$) but not for 4He2+ ($m_{i}/Z_{i}m_{p}=2$) and 3He2+
($m_{i}/Z_{i}m_{p}=1.5$).

The localization of energized 16O7+ ions near the upper (lower) separatrix in
the left (right) outflow region at $\omega_{cp}t=19.6$ is due to the
asymmetries in the lower and upper separatrices in guide field magnetic
reconnection (Li & Zhang, 2020; Pritchett & Coroniti, 2004). The thickness of
the upper (lower) separatrix in the left (right) outflow region is smaller
than that of the lower (upper) separatrix. This limits the non-adiabaticity
and thus energization of ions entering outflow regions from the thicker
separatrix.

The increased rate of energization for all the heavy ion species after
$\omega_{cp}t=30$ is somehow linked to the growth of the normal component
$B_{y}$ of the magnetic field in the current sheet (Fig. 2). The normal
component $B_{y}$ in current sheet can grow large by non-steady reconnection
with increasing rate and/or the compression along the current sheet of the
reconnected magnetic field lines. The compression can occur due to the pile up
of the magnetic flux reconnected at two neighboring sites on the magnetic
island between the two sites, contraction of magnetic islands and/or mutual
pushing among magnetic islands. Several acceleration mechanisms may be
associated with these scenarios: direct acceleration by inductive reconnection
electric field in the X-point regions, acceleration by motional electric field
induced by Alfvénic outflow in the outflow regions, Fermi-like acceleration in
contracting magnetic islands, magnetic curvature and gradient drifts aligned
with inductive/motional electric field in the flux pile-up regions and
betatron acceleration by time dependent magnetic field in the flux compression
regions.

Figures 3 and 4 show that, at $\omega_{cp}t=48.99$, energized ions are located
in the X-point regions, inside magnetic islands and outflow exhaust regions.
Most energized particles (black dots in Figures 3 and 4) are mostly
concentrated near the opening of the exhaust regions, where the reconnected
magnetic field lines usually pile up. In the case of simulations with 16O7+,
merging of two magnetic islands and presence of energized ions at $z\approx 0$
can also be seen (Fig. 4). Although the locations of the energetic particles
at a given time are not necessarily in the neighborhood of their acceleration
sites, it seems from Figs. 3 and 4 that a number of acceleration mechanisms
out of those discussed above are operating simultaneously in different regions
of reconnection. A given particle might experience acceleration due to
different mechanisms it encounters on its trajectory in different reconnection
regions. Disentangling these mechanisms require a fully history of particles’
motion and therefore tracking of their trajectories and would be the subject
of our future studies.

$\omega_{cp}t=19.6$$\omega_{cp}t=48.99$

Figure 3: Positions of 4He2+ ions represented by dots colored by ion’s kinetic
energy at $\omega_{cp}t=19.6$ (left column) and $\omega_{cp}t=48.99$ (right
column). Lines with arrows represent magnetic field lines. Only the lower half
($y<0$) of the simulation domain is shown.

$\omega_{cp}t=19.6$$\omega_{cp}t=48.99$

Figure 4: Positions of 16O7+ ions represented by dots colored by ion’s kinetic
energy at $\omega_{cp}t=19.6$ (left column) and $\omega_{cp}t=48.99$ (right
column). Lines with arrows represent magnetic field lines. Only the lower half
($y<0$) of the simulation domain is shown.

Fig. 5 shows the evolution of the energy spectra for different ions species.
Energy spectra for all the ion species broaden with time. Energy spectra for
protons did not develop any noticeable non-thermal feature up to the final
simulation time $\omega_{cp}t=48.99$. Thus energy transferred to protons only
heat them. Non-thermal feature, a shoulder in the spectra in the intermediate
energy range after which spectra falls off rapidly, develops in the spectra of
4He2+ and 16O7+ at $\omega_{cp}t=39.19$ and $\omega_{cp}t=19.6$, respectively.
The early development of the spectral shoulder for heavier ions, 16O7+ and
56Fe14+ (Figure not shown), in comparison to the lighter ions, 4He2+ and 4He3+
(Figure not shown), is consistent with the early energization of heavier ions
by pickup mechanism. At $\omega_{cp}t=19.6$, the energy of the accelerated
16O7+ (see Fig. 4) is in the range 1-10 $m_{p}v_{Ap}^{2}$ in which 16O7+
spectra develop a shoulder. The spectral shoulder rises with time as well as
extends to higher energies as more and more low energy heavy ions are
accelerated to higher energies.

(a)(b)(c)

Figure 5: Energy spectra of protons 1H1+ (a), 4He2+ (b) and 16O7+ (c) at
different times.

(a)(b)

Figure 6: Energy-per-nucleon spectra at $\omega_{cp}t=48.99$. Original spectra
of different ion species (a) and spectra of heavy ions shifted for visibility
on y-axis by a factor of $10^{0}$, $10^{1}$, $10^{2}$ and $10^{3}$ for 4He2+,
3He2+, 16O7+ and 56Fe14+, respectively (b). The dashed lines in the right
panel are the power law fit ($E^{\gamma}$) to the shoulder region of the
spectra where $E$ is the energy-per-nucleon.

Figure 6a compares the energy-per-nucleon spectra of different ion species at
$\omega_{cp}t=48.99$ — the last time up to which our simulations are valid.
Note that, for the purpose of comparing the spectra of different ion species
independent of their total mass, the spectra in Fig. 6a is shown as a function
of energy-per-nucleon unlike the spectra in Fig.5 which is shown as a function
of energy. The injection energy $\sim 0.3m_{p}v_{Ap}^{2}$, around which the
spectral shoulder begins, is similar for heavy ion species. The cutoff energy,
up to which the shoulder extends, and the maximum gained energy per nucleon
are larger for the lighter ion species. In Fig. 6b, we plot the same spectra
(excluding proton’s spectra) as in Fig. 6a but artificially shifted on y-axis
for better visibility of the spectral shoulder and power law fits,
$E^{\gamma}$, of the shoulder region of the spectra, where $E$ is the energy-
per-nucleon. The spectral indices of the fits, $\gamma$=-1.59, -1.53 and -1.68
respectively for 4He2+, 3He2+ and 16O7+ are quite similar. The spectra for
56Fe14+, on the other hand, is somewhat softer with $\gamma=-3.0$. The
dependence of $\gamma$ on $Q/M$, shown in Fig. 7a, can be approximately fitted
as $\gamma\propto(Q/M)^{-0.64}$.

(a)(b)

Figure 7:  Spectral index $\gamma$ of the energy-per-nucleon spectra in the
shoulder region vs. $Q/M$ and the fit $\gamma\propto(Q/M)^{-0.64}$ (a).
Scaling of the proxy of the ion abundance, defined as number of particles
above a threshold energy normalized to the total number of particles, with
$Q/M$ for threshold energy 10 (red filled circles), 25 (blue plus signs) and
50 (black filled squares) times the initial thermal energy $KE_{th}$ (b).
Straight lines are the fit, abundance $\propto(Q/M)^{-\alpha}$, for the data
and are plotted with the color corresponding to the data.

In 3He-rich SEP events, ion abundances are enhanced relative to their solar
abundances at energies well above their thermal energy. We consider a proxy
for abundance enhancement factor $F$ as number of particles above a threshold
energy normalized to the total number of particles. Fig. 7b shows this proxy
of the abundance enhancement factor as a function of $Q/M$ at
$\omega_{cp}t=48.99$ for three energy thresholds
$E_{t}=10\,\mathit{KE}_{th},25\,\mathit{KE}_{th}$ and $50\,\mathit{KE}_{th}$
well above the initial thermal energy $\mathit{KE}_{th}=0.25m_{p}v_{Ap}^{2}$
which is the same for all ion species in our simulations. The value of $E_{t}$
is so chosen that it is in the tail region of the initial energy spectra which
is the same for all the ion species and at the same time is not outside the
shoulder regions of the energy spectra of the heavy ions at
$\omega_{cp}t=48.99$. Although the value
$E_{t}=50\,\mathit{KE}_{th}=12.5\,m_{p}v_{Ap}^{2}$ is close to the cutoff
energy of the shoulder region in the energy spectra of 4He2+ at
$\omega_{cp}t=48.99$ (see Fig. 5b), the abundances will be dominated by the
abundances in the shoulder region as the energy spectra falls off very rapidly
beyond the shoulder region. Note that the threshold
$E_{t}=12.5\,m_{p}v_{Ap}^{2}$ is in the rapidly falling region of the proton’s
energy spectra at $\omega_{cp}t=48.99$. It is, however, acceptable as proton’s
energy spectra does not develop a shoulder.

Fig. 7b also shows the power law fit $F\propto(Q/M)^{-\alpha}$. The
enhancement factor drops with increasing $Q/M$, i.e., heavier ions are
preferentially accelerated. The drop of the enhancement factor becomes steeper
(the value of $\alpha$ increases) with the increasing threshold energy. As the
value of $E_{t}$ increases, the contribution of the shoulder region relative
to the contribution of the rapidly falling region of the energy spectra
towards the abundances decreases. This decrease in the contribution is more
prominent for larger values of $Q/M$ as the shoulder in their energy spectra
extends to smaller energies. This results in the increase in the value of
$\alpha$ with $E_{t}$.

##  4 Discussion and conclusion

We carried out hybrid-kinetic simulations (with electron inertia) to study the
acceleration of heavy ions ($Q/M<1$) by magnetic reconnection. We find that
heavy ions can be accelerated to high energies many times larger than their
initial energies by a variety of acceleration mechanisms. Heavier ions are
preferentially accelerated in the sense that energy gain averaged over
particles increases with decreasing $Q/M$. They are primarily accelerated in
magnetic islands and flux pile up regions near the opening of the outflow
exhausts. Most efficient acceleration takes place in the flux pileup regions.
Heavy ions, depending upon the smallness of $Q/M$ which allows them to be non-
adiabatic while crossing from inflow to outflow regions, can also be
accelerated by pickup mechanism in outflow regions even before any magnetic
flux is piled up. As a result of acceleration, heavy ions develop a shoulder,
a non-thermal feature, in their energy spectra. The spectral index obtained
from the power law fit in the shoulder region of the spectra varies
approximately as $(Q/M)^{-0.64}$. Abundance enhancement factor, defined as
number of particles above a threshold energy normalized to total number of
particles, scales as $(Q/M)^{-\alpha}$ where $\alpha$ increases with the
energy threshold.

Energy spectra with a shoulder or in other words double power law with a break
in the spectra at energy $\sim$ 1 MeV/nucleon has been in-situ observed in
space (Mason et al., 2002; Buc̆ík et al., 2018). Our simulations show the
break in the energy spectra in the energy range 2-6 $m_{p}v_{Ap}^{2}$ per
nucleon depending upon the ion species. The value of $v_{Ap}$ in the active
regions of solar corona, where acceleration takes place, has been estimated to
lie in the range 2000-9000 km/s (Brooks et al., 2021). For a typical value
$v_{Ap}$=5000 km/s, $m_{p}v_{Ap}^{2}\approx$0.25 MeV and therefore break in
the simulations occurs in the range 0.5-1.5 MeV/nucleon, consistent with the
observations. Note that the energy range of the break depends on the value of
$v_{Ap}$, estimates of which in solar corona vary significantly (Brooks et
al., 2021). In the observations, the spectral index $\gamma$ of the power law
before the break is in the range 1-3 for the ions of Helium-4, Helium-3,
Oxygen and Iron. The values of $\gamma$ before the break is in the same range
in our simulations as well. The values of $\gamma$ after the break are,
however, much larger in our simulation in comparison to the observations.

We defined a proxy $F$ for the abundance enhancement factor to study its
variation with $Q/M$. Although this proxy does not exactly correspond to the
abundance enhancement factor used in observations, we point out some
similarity and differences in the behavior of the two. In our simulations, the
power law index $\alpha$ of the fit $F\propto(Q/M)^{-\alpha}$ increases with
the threshold energy. Observation also show event to event variation in the
value of $\alpha$ in the range $-5<\alpha<12$ with a mean value of
$\alpha=$3.26 at 385 keV/nucleon (Mason et al., 2004) and $2<\alpha<8$ with a
mean value of $\alpha=3.64$ at 3-5 MeV/nucleon (Reames et al., 2014a, b). The
mean values in the observations are similar despite their energies being an
order of magnitude apart. The mean of the three values of $\alpha$ obtained
from power law fits in Fig. 7b is approximately 3.66. In observations, steeper
energy spectra tend to have steeper fall of the abundance enhancement with
$Q/M$, i.e., large value of $\alpha$ (Reames et al., 2014a). This effect is
similar to the increase in the value of $\alpha$ with $E_{t}$ in our
simulations. Larger values of $E_{t}$ increases the contribution of the
steeper part of the energy spectra in the calculations of the abundances.

Note that the SEP events are detected in space far away from their
acceleration sites on Sun. The observed scaling of abundance enhancement with
$Q/M$ and other features of these events may, therefore, not necessarily be
the effect of only acceleration but also of the transport of particles from
the acceleration site to the detection site. In fact, 3-D test particle
modeling of the inter-planetary transport of relativistic protons from the
source region on Sun shows that the spectra of the particles is highly
observer dependent and do not necessarily reflect the source spectra (Dalla et
al., 2020). Such transport effects are possible for the spectra of heavy ions
as well. Nevertheless, magnetic reconnection is a potential candidate for the
preferential acceleration of heavy ions which may provide a power law
dependence of the abundance enhancement on $Q/M$, as suggested by the results
presented here. More detailed studies on the relative roles of the different
acceleration mechanisms operating in magnetic reconnection in the abundance
enhancement will be presented in a future publication. Note that we did not
find any extra-ordinary abundance enhancement of 3He2+ in our simulations.
This, however, does not rule out the role of magnetic reconnection in the
abundance enhancement of 3He2+ as our simulations are carried out only for
limited range of parameters.

The authors thank Patricio Muñoz for fruitful discussion and his help. We
gratefully acknowledge the financial support by the German Science Foundation
(DFG), projects JA 2680-2-1 and BU 777-17-1 as well as the Czech Republic GAĈR
project 20-09922J M, personally Markus Rampp and Meisam Tabriz of the Max
Planck Computing and Data Facility (MPCDF) for their support of the
development of the hybrid-kinetic code CHIEF used for the simulation studies
presented in this paper. The simulations were carried out on the MPS
supercomputers at the MPCDF, Garching, Germany.

## References

  * Bárta et al. (2011a) Bárta, M., Büchner, J., Karlický, M., & Kotrč, P. 2011a, The Astrophysical Journal, 730, 47, doi: 10.1088/0004-637X/730/1/47
  * Bárta et al. (2011b) Bárta, M., Büchner, J., Karlický, M., & Skála, J. 2011b, The Astrophysical Journal, 737, 24, doi: 10.1088/0004-637X/737/1/24
  * Brooks et al. (2021) Brooks, D. H., Warren, H. P., & Landi, E. 2021, The Astrophysical Journal Letters, 915, L24, doi: 10.3847/2041-8213/ac0c84
  * Bučík et al. (2021) Bučík, R., Mulay, S. M., Mason, G. M., et al. 2021, The Astrophysical Journal, 908, 243, doi: 10.3847/1538-4357/abd62d
  * Buc̆ík et al. (2018) Buc̆ík, R., Wiedenbeck, M. E., Mason, G. M., et al. 2018, The Astrophysical Journal Letters, 869, L21, doi: 10.3847/2041-8213/aaf37f
  * Buc̆ík et al. (2023) Buc̆ík, R., Mason, G. M., Gómez-Herrero, R., et al. 2023, Astronomy & Astrophysics, 669, A13, doi: 10.1051/0004-6361/202245037
  * Bučík (2020) Bučík, R. 2020, Space Science Reviews, 216, 24, doi: 10.1007/s11214-020-00650-5
  * Cane et al. (1986) Cane, H. V., McGuire, R. E., & von Rosenvinge, T. T. 1986, The Astrophysical Journal, 301, 448
  * Dalla et al. (2020) Dalla, S., de Nolfo, G. A., Bruno, A., et al. 2020, A&A, 639, A105, doi: 10.1051/0004-6361/201937338
  * Desai & Giacalone (2016) Desai, M., & Giacalone, J. 2016, Living Reviews in Solar Physics, 13, 3, doi: 10.1007/s41116-016-0002-5
  * Desai et al. (2016) Desai, M. I., Mason, G. M., Dayeh, M. A., et al. 2016, The Astrophysical Journal, 816, 68, doi: 10.3847/0004-637X/816/2/68
  * Drake et al. (2009) Drake, J. F., Cassak, P. A., Shay, M. A., Swisdak, M., & Quataert, E. 2009, The Astrophysical Journal, 700, L16, doi: 10.1088/0004-637X/700/1/L16
  * Eichler (2014) Eichler, D. 2014, The Astrophysical Journal, 794, 6, doi: 10.1088/0004-637X/794/1/6
  * Fisk (1978) Fisk, L. A. 1978, Astrophysical Journal, 224, 1048, doi: 10.1086/156456
  * Fu et al. (2020) Fu, X., Guo, F., Li, H., & Li, X. 2020, The Astrophysical Journal, 890, 161, doi: 10.3847/1538-4357/ab6d68
  * Jain et al. (2022) Jain, N., Muñoz, P. A., Farzalipour Tabriz, M., Rampp, M., & Büchner,Jörg. 2022, Physics of Plasmas, 29, 053902, doi: 10.1063/5.0087103
  * Kazakov et al. (2017) Kazakov, Y. O., Ongena, J., Wright, J. C., et al. 2017, Nature Physics, 13, 973, doi: 10.1038/nphys4167
  * Knizhnik et al. (2011) Knizhnik, K., Swisdak, M., & Drake, J. F. 2011, The Astrophysical Journal Letters, 743, L35, doi: 10.1088/2041-8205/743/2/L35
  * Kocharov & Kocharov (1984) Kocharov, L. G., & Kocharov, G. E. 1984, Space Sci. Rev., 38, 89.  https://doi.org/10.1007/BF00180337
  * Kramoliš et al. (2022) Kramoliš, D., Bárta, M., Varady, M., & Bučík, R. 2022, The Astrophysical Journal, 927, 177, doi: 10.3847/1538-4357/ac4fc9
  * Kumar et al. (2017) Kumar, R., Eichler, D., Gaspari, M., & Spitkovsky, A. 2017, The Astrophysical Journal, 835, 295, doi: 10.3847/1538-4357/835/2/295
  * Li & Zhang (2020) Li, Z., & Zhang, M. 2020, The Astrophysical Journal, 888, 5, doi: 10.3847/1538-4357/ab5aea
  * Liu et al. (2006) Liu, S., Petrosian, V., & Mason, G. M. 2006, The Astrophysical Journal, 636, 462, doi: 10.1086/497883
  * Liu et al. (2017) Liu, Y.-H., Hesse, M., Guo, F., et al. 2017, Phys. Rev. Lett., 118, 085101, doi: 10.1103/PhysRevLett.118.085101
  * Mason (2007) Mason, G. M. 2007, Space Sci. Rev., 130, 231, doi: 10.1007/s11214-007-9156-8
  * Mason et al. (2004) Mason, G. M., Mazur, J. E., Dwyer, J. R., et al. 2004, The Astrophysical Journal, 606, 555, doi: 10.1086/382864
  * Mason et al. (1986) Mason, G. M., Reames, D. V., Klecker, B., Hovestadt, D., & von Rosenvinge, T. T. 1986, The Astrophysical Journal, 303, 849, doi: 10.1086/164133
  * Mason et al. (2002) Mason, G. M., Wiedenbeck, M. E., Miller, J. A., et al. 2002, The Astrophysical Journal, 574, 1039, doi: 10.1086/341112
  * Miller (1998) Miller, J. A. 1998, Space Science Reviews, 86, 79, doi: 10.1023/A:1005066209536
  * Muñoz et al. (2018) Muñoz, P., Jain, N., Kilian, P., & Büchner, J. 2018, Computer Physics Communications, 224, 245, doi: https://doi.org/10.1016/j.cpc.2017.10.012
  * Muñoz et al. (2023) Muñoz, P. A., Jain, N., Tabriz, M. F., Rampp, M., & Jörg Büchner. 2023, Electron inertia effects in 3D hybrid-kinetic collisionless plasma turbulence.  https://arxiv.org/abs/2303.01322
  * Pritchett & Coroniti (2004) Pritchett, P. L., & Coroniti, F. V. 2004, Journal of Geophysical Research: Space Physics, 109, A01220, doi: https://doi.org/10.1029/2003JA009999
  * Reames (1999) Reames, D. V. 1999, Space Sci. Rev., 90, 413.  https://doi.org/10.1023/A:1005105831781
  * Reames (2013) —. 2013, Space Science Reviews, 175, 53, doi: 10.1007/s11214-013-9958-9
  * Reames (2021a) —. 2021a, Lecture Notes in Physics, Vol. 978, Solar Energetic Particles: A Modern Primer on Understanding Sources, Acceleration and Propagation (Cham: Springer International Publishing), doi: 10.1007/978-3-030-66402-2
  * Reames (2021b) —. 2021b, Space Science Reviews, 217, 72, doi: 10.1007/s11214-021-00845-4
  * Reames et al. (2014a) Reames, D. V., Cliver, E. W., & Kahler, S. W. 2014a, Solar Physics, 289, 4675, doi: 10.1007/s11207-014-0589-4
  * Reames et al. (2014b) —. 2014b, Sol. Phys., 289, 3817.  https://doi.org/10.1007/s11207-014-0547-1
  * Shi et al. (2022) Shi, Z., Muñoz, P. A., Jörg Büchner, & Liu, S. 2022, The Astrophysical Journal, 941, 39, doi: 10.3847/1538-4357/ac9fd7
  * Temerin & Roth (1992) Temerin, M., & Roth, I. 1992, The Astrophysical Journal, 391, L105, doi: 10.1086/186408
  * Zhou et al. (2015) Zhou, X., , Büchner, J., et al. 2015, The Astrophysical Journal, 815, 6, doi: 10.1088/0004-637X/815/1/6
  * Zhou et al. (2016) Zhou, X., Büchner, J., Bárta, M., Gan, W., & Liu, S. 2016, The Astrophysical Journal, 827, 94, doi: 10.3847/0004-637X/827/2/94