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

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 ³He 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 ³He-rich events. The enhancement of heavy ions is, however, not correlated with the extra-ordinary enhancement of ³He (Mason et al., 1986; Reames, 1999). Abundance enhancement of ³He 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 ³He 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 ³He 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 ³He consider absorption of some wave energy by cyclotron resonance of ³He 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 ³He 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 (⁴He²⁺), 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.

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 : ⁴He²⁺, ³He²⁺, ¹⁶O⁷⁺ and ⁵⁶Fe¹⁴⁺. 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.

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$.

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$.

Heavier ions (¹⁶O⁷⁺ and ⁵⁶O¹⁴⁺ 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 ⁴He²⁺ (Fig. 3) and the other with ¹⁶O⁷⁺ (Fig. 4). At $\omega_{cp}t=19.6$, ¹⁶O⁷⁺ ions are more energized than ⁴He²⁺ ions and the locations of the energized ¹⁶O⁷⁺ 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 ⁴He²⁺ and ³He²⁺ 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 ⁵⁶Fe¹⁴⁺ ($m_{i}/Z_{i}m_{p}=4$) and only marginally for ¹⁶O⁷⁺ ($m_{i}/Z_{i}m_{p}=2.28$) but not for ⁴He²⁺ ($m_{i}/Z_{i}m_{p}=2$) and ³He²⁺ ($m_{i}/Z_{i}m_{p}=1.5$).

The localization of energized ¹⁶O⁷⁺ 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 ¹⁶O⁷⁺, 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.

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 ⁴He²⁺ and ¹⁶O⁷⁺ at $\omega_{cp}t=39.19$ and $\omega_{cp}t=19.6$, respectively. The early development of the spectral shoulder for heavier ions, ¹⁶O⁷⁺ and ⁵⁶Fe¹⁴⁺ (Figure not shown), in comparison to the lighter ions, ⁴He²⁺ and ⁴He³⁺ (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 ¹⁶O⁷⁺ (see Fig. 4) is in the range 1-10 $m_{p}v_{Ap}^{2}$ in which ¹⁶O⁷⁺ 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.

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 ⁴He²⁺, ³He²⁺ and ¹⁶O⁷⁺ are quite similar. The spectra for ⁵⁶Fe¹⁴⁺, 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}$.

In ³He-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 ⁴He²⁺ 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}$.

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 ³He²⁺ in our simulations. This, however, does not rule out the role of magnetic reconnection in the abundance enhancement of ³He²⁺ as our simulations are carried out only for limited range of parameters.