© The Authors 2024

1. Introduction

Nonthermal 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⁴ − 10²⁰ 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 the Milky Way, and beyond our Galaxy, respectively. SEPs in a typical energy range of 10 s KeV nucleon⁻¹–100 s MeV nucleon⁻¹ are of particular interest because their flux near Earth is sufficiently high to have implications for space weather phenomena.

Solar energetic particle events can be classified into impulsive events (short duration ≤ 1 day; less intense, smaller particle fluxes with typical energies < 10 MeV nucleon⁻¹; and numerous, about 1000 per year) and gradual events (long duration, several days; orders of magnitude more intense, large particle fluxes with energies > 10 MeV nucleon⁻¹; less frequent ∼ 10 per year) (Reames 1999, 2021b). Gradual events are attributed to the acceleration of charged particles by coronal mass ejection (CME) driven shocks, while impulsive events are attributed to the acceleration in solar flares by magnetic reconnection and associated processes (Kahler et al. 2001; Mason 2007; Drake et al. 2009; Lin 2011; Reames 2013, 2022; Murphy et al. 2016; Desai & Giacalone 2016; Bučík 2020). Impulsive events show abundance enhancements of heavier ions relative to their abundances in the solar corona at energies well above their average thermal energy (∼100 eV) (Reames 2021a). In contrast, the average abundances of elements in gradual events are similar to those in solar corona (e.g., Reames 2018). 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 an ion’s charge (Q) to mass (M) ratio as (Q/M)^−α for M ≥ 4 amu. Variations in the value of α have been found from event to event (Reames et al. 2014b) with mean values α = 3.26 (at 0.375 MeV nucleon⁻¹; Mason et al. 2004), α = 3.64 ± 0.15 (at 3–10 MeV nucleon⁻¹; Reames et al. 2014a), α = 3.53 (at 160–226 KeV nucleon⁻¹), and α = 3.31 (at 320–453 KeV nucleon⁻¹; Bučík et al. 2021). The enhancement factor for the ³He isotope, however, does not obey the power law and can have very large values (up to 10⁴) 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. However, the enhancement of heavy ions is not correlated with the extraordinary enhancement of ³He (Mason et al. 1986; Reames 1999).

The abundance enhancements of heavy ions with low values of Q/M and the extraordinary abundance enhancement of ³He in impulsive SEP events are not understood yet. Since these enhancements are uncorrelated (Mason et al. 1986; Reames 1999), attempts have been made to explain the abundance enhancement of heavy ions and ³He by separate mechanisms. Abundance enhancement of heavy ions is thought to be due to their preferential acceleration by magnetic reconnection and the associated processes in solar flares from thermal to SEP energies. Models of the preferential acceleration of heavy ions consider resonant interaction of the ions with the waves in Alfvénic plasma turbulence (Miller 1998; Eichler 2014; Kumar et al. 2017; Fu et al. 2020; Shi et al. 2022) that may be generated during magnetic reconnection, for example the development of turbulent outflows, as shown by fully kinetic 3D simulations of magnetic reconnection (Lapenta et al. 2020). Turbulent energy decays with wave frequency, which means that higher frequency waves have lower energy. Heavy ions with lower values of Q/M have lower values of cyclotron frequency (∝Q/M), and thus would resonate with lower frequency but higher power waves in turbulence, favoring the preferential acceleration of heavy ions. In addition, turbulence amplitudes above a threshold can make ion orbits chaotic, thus heating ions stochastically. Phenomenological arguments and simulations of test particles interacting with turbulent spectra of Alfvén and kinetic Alfvén waves at wavelengths comparable to ion gyro-radius in plasmas with β < 1 found that, at comparable temperatures, alpha particles and heavier ions are stochastically heated more efficiently than the protons (Chandran et al. 2010). Three-dimensional hybrid-kinetic simulations of continuously driven Alfvénic turbulence at low plasma beta show that the Hall and thermoelectric effects in a generalized Ohm’s law contribute significantly to the stochastic ion heating, which is enhanced by the intermittency in the turbulence (Cerri et al. 2021). 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, for example by electrons streaming along open magnetic field lines in solar corona, or via coupling with low-frequency Alfvén waves. 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 by nonresonant mechanisms, which have also been considered to explain the preferential acceleration of heavy ions (Drake et al. 2009; Bárta et al. 2011b,a; Zhou et al. 2015, 2016). In 2.5D magnetohydrodynamic (MHD) and test particle simulations (Kramoliš et al. 2022), heavy ions were found to be preferentially accelerated by first-order Fermi processes 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 2D fully kinetic simulations of magnetic reconnection, ions entering the reconnection exhaust can behave like pickup particles and can be accelerated if their Q/M is below a threshold value, resulting in the 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 only for one ion species (⁴He²⁺), did not verify the power-law dependence. Therefore, it is not clear whether the acceleration of heavy ions by pickup mechanisms can provide the observed scaling of the abundance enhancement with Q/M. Magnetic reconnection in turbulence generated current sheets, which form in turbulence as a result of the magnetic energy cascade from large to kinetic scales where magnetic fluctuations have intermittent statistics (Cerri et al. 2019), have also been suggested to contribute to the energization of charged particles (Rueda et al. 2022; Grošelj et al. 2017).

Magnetic reconnection in collisionless plasmas is a multi-scale process that occurs at electron kinetic scales and then couples to ion scales and even larger macro-scales. An ideal simulation of magnetic reconnection requires the 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 the particle-in-cell, PIC, method) of magnetic reconnection in electron-proton plasma are computationally very demanding. The self-consistent inclusion of heavier ion species make the simulations even more demanding because now it is necessary to use a larger simulation box to accommodate the larger gyro-radii of the heavier species, while at the same time resolving the electron scales. For this reason a variety of reduced kinetic models have been used for numerical studies of magnetic reconnection, plasma turbulence, and the associated charge particle heating. Plasma turbulence at ion scale was studied using hybrid-kinetic simulations (with and without electron inertia) along with fully kinetic simulations (Cerri et al. 2019). Two-dimensional fully kinetic, gyro-kinetic, and hybrid-kinetic simulations of plasma turbulence suggest that magnetic reconnection can heat electrons and ions in turbulence generated current structures (Grošelj et al. 2017). Three-dimensional hybrid-kinetic simulations of driven Alfvén-wave turbulence show a flattened core in the perpendicular-velocity distribution and non-Maxwellian wings in the parallel-velocity distribution of ions resulting from the dissipation of kinetic Alfvén waves at sub-ion scales (Arzamasskiy et al. 2019). Two-dimensional simulations of guide field magnetic reconnection by means of a hybrid-Vlasov-Landau-fluid (HVLF) model reproduces the main features of magnetic reconnection and anisotropic heating of electrons (Finelli et al. 2021). In hybrid-kinetic simulations (with massless electron fluid) of plasma turbulence, preferential heating of alpha particles with respect to protons was observed near current sheets in spatial regions of enhanced ion vorticity (Valentini et al. 2016).

In this paper, we employ 2D hybrid-kinetic simulations (kinetic ions and inertial electron fluid) to study the acceleration of heavy ions and its consequences for their abundance enhancements by magnetic reconnection, an established acceleration source for impulsive SEP events in which heavy-ion enhancement is observed (Kahler et al. 2001; Mason 2007; Drake et al. 2009; Lin 2011; Murphy et al. 2016; Reames 2022). Treatment of electrons as an inertial fluid allows inclusion of the physics at electron inertial scale, but relaxes the numerical requirement of resolving the Debye length in the PIC method, and therefore allows larger simulation domains for the same number of grid points. The computational feasibility, however, comes at the cost of electron kinetic physics, which is acceptable as long as the study of the ion acceleration is the primary objective. 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 higher than the proton’s mass. The paper is organized as follows. Section 2 presents the simulation setup. The simulation results are presented in Sect. 3, and are discussed and concluded in Sect. 4.

2. Simulation setup

We carried out 2D hybrid-kinetic simulations of magnetic reconnection with electron inertia using a hybrid-PIC code called “Code Hybrid with Inertial Electron Fluid (CHIEF)”, which is a 3D parallelized code based on message passing interface (MPI) for high-performance computing (Muñoz et al. 2018, 2023; Jain et al. 2022). 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, 2023; Jain et al. 2022).

We initialize the 2D simulations with two Harris current sheets, $\mathbf{B}=\widehat{z}{B}_{z0}[tanh\{(\mathit{y}+L_{\mathit{y}}/4)/L\}-tanh\{(\mathit{y}-L_{\mathit{y}}/4)/L\}-1]+\widehat{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 the 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₀[sech²{(y + L_y/4)/L}+sech²{(y − L_y/4)/L}]. Harris sheets are embedded in a uniform background plasma of density 0.2 n₀ and therefore peak density at their centers is 1.2 n₀. 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₀, n_bi = 0.01 n₀ (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: $T_{\text{e}}=T_{\text{p}}=T_i=0.25m_{\text{p}}{V}_{\text{Ap}}^2$, where m_p is proton mass and $V_{\mathrm{Ap}}=B_{z0}/\sqrt{{\mu }_0{n}_0{m}_{\mathrm{p}}}$ is the proton Alfvén velocity based on B_z0 and n₀.

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 the proton-to-electron mass ratio as m_p/m_e = 25. The simulation domain L_y × L_z is 51.2 d_p × 102.4 d_p resolved by a grid spacing of 0.1 d_p in each direction. The time step is $\mathrm{\Delta }t=0.0025{\omega }_{\mathrm{cp}}^{-1}$, where ω_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

Figure 1 shows the evolution of the out-of-plane current density J_x from the initial to the late phase. In the initial phase (ω_cpt = 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 ω_cpt = 39.19, the lower and the upper current sheets 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 sheets grow in size with time by merging and/or pushing among themselves and simultaneously develop turbulence inside them. At ω_cpt = 58.79 the magnetic islands of the lower and upper current sheets grow to a big enough size that the particles accelerated in the upper (lower) current sheet may cross over to the lower (upper) current sheet.

Fig. 1. Out-of-plane current density Jx/(n0eVAp) at four different times.

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 the z-direction after a time, $\sim 50{\omega }_{\mathrm{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 ω_cpt = 48.99 and in the spatial region −23 ≤ y/d_p ≤ −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 ω_cpt ≈ 30 after which they gain energy at much faster rates. We note that magnetic islands have significantly developed by ω_cpt = 30. This can be seen in Fig. 2b, which shows that ${\langle B_y^2(y=-12.5d_{\text{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 ω_cpt = 20, and by ω_cpt = 30 has grown to ∼5% of the asymptotic magnetic energy in the anti-parallel magnetic field. The development of magnetic islands, therefore, seems 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 ω_cpt ≈ 30.

Fig. 2. Time development of the fractional change in average kinetic energy per particle, ⟨K.E.(t)⟩ptcl, from its initial value (a) and average of the square of the normal component of the magnetic field By along the central line of the lower current sheet (b) for different ion species.

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 < ω_cpt < 30). This suggests that acceleration mechanisms, different from the mechanisms after ω_cpt = 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 ω_cpt = 30 (the time by which ∼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 ω_cpt = 19.6 and 48.99 in two simulations, one with ⁴He²⁺ (Fig. 3) and the other with ¹⁶O⁷⁺ (Fig. 4). At ω_cpt = 19.6, ¹⁶O⁷⁺ ions are more energized than ⁴He²⁺ ions and the locations of the energized ¹⁶O⁷⁺ ions are in the outflow regions. However, the energization is 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}_{\mathrm{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. (2009), where m_i and Z_i = Q_i/e are the mass and charge state of ions and β_p is the proton plasma beta based on asymptotic value of the anti-parallel magnetic field. Although this condition, which becomes m_i/Z_im_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_im_p = 2 and 1.5 were not energized during ω_cpt < 30.

Fig. 3. Positions of 4He2+ ions represented by dots color-coded according to the ions’ kinetic energy at ωcpt = 19.6 (left column) and ωcpt = 48.99 (right column). The lines with arrows represent magnetic field lines. Only the lower half (y < 0) of the simulation domain is shown.

It should be noted 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 ω_cpt = 19.6 is still developing and the ion inflow velocity has not yet reached its maximum value. Using |u_iy|≈0.05 v_Ap, the value of inflow velocity at ω_cpt = 19.6 in our simulations, the threshold condition becomes m_i/Z_im_p > 2.24, which allows non-adiabatic behavior for ⁵⁶Fe¹⁴⁺ (m_i/Z_im_p = 4) and only marginally for ¹⁶O⁷⁺ (m_i/Z_im_p = 2.28), but not for ⁴He²⁺ (m_i/Z_im_p = 2) and ³He²⁺ (m_i/Z_im_p = 1.5).

The localization of energized ¹⁶O⁷⁺ ions near the upper (lower) separatrix in the left (right) outflow region at ω_cpt = 19.6 is due to the asymmetries in the lower and upper separatrices in guide field magnetic reconnection (Pritchett & Coroniti 2004; Li & Zhang 2020). 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 the energization of ions entering the outflow regions from the thicker separatrix.

The increased rate of energization for all the heavy ion species after ω_cpt = 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 nonsteady 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 pileup 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 fields in the X-point regions, acceleration by motional electric fields induced by Alfvénic outflow in the outflow regions, Fermi-like acceleration in contracting magnetic islands, magnetic curvature and gradient drifts aligned with inductive or motional electric fields in the flux pileup regions, and betatron acceleration by time dependent magnetic fields in the flux compression regions.

Figures 3 and 4 show that, at ω_cpt = 48.99, energized ions are located in the X-point regions, inside magnetic islands and outflow exhaust regions. Most energized particles (black dots in Figs. 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⁷⁺, the merging of two magnetic islands and the presence of energized ions at z ≈ 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 requires a full history of the particles’ motion and therefore the tracking of their trajectories, and will be the subject of our future studies.

Fig. 4. Positions of 16O7+ ions represented by dots color-coded according to the ions’ kinetic energy at ωcpt = 19.6 (left column) and ωcpt = 48.99 (right column). The lines with arrows represent magnetic field lines. Only the lower half (y < 0) of the simulation domain is shown.

Figure 5 shows the evolution of the energy spectra for different ions species. The energy spectra for all the ion species broaden with time. The energy spectra for protons did not develop any noticeable nonthermal features up to the final simulation time ω_cpt = 48.99. Thus, the energy transferred to protons only heats them. A nonthermal 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 ω_cpt = 39.19 and ω_cpt = 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 ω_cpt = 19.6, the energy of the accelerated ¹⁶O⁷⁺ (see Fig. 4) is in the range 1–10 $m_{\text{p}}{v}_{\text{Ap}}^2$ at which ¹⁶O⁷⁺ spectra develop a shoulder. The spectral shoulder rises with time and extends to higher energies as more and more low-energy heavy ions are accelerated to higher energies.

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

Figure 6a compares the energy-per-nucleon spectra of different ion species at ω_cpt = 48.99, the last time up to which our simulations are valid. We note that, for the purpose of comparing the spectra of different ion species independent of their total mass, the spectra in Fig. 6a are shown as a function of energy-per-nucleon unlike the spectra in Fig. 5 which are shown as a function of energy. The injection energy $\sim 0.3m_{\text{p}}{v}_{\text{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 higher for the lighter ion species. In Fig. 6b we plot the same spectra (excluding the proton spectra) as in Fig. 6a, but artificially shifted on the y-axis for better visibility of the spectral shoulder and power-law fits, E^γ, of the shoulder region of the spectra, where E is the energy-per-nucleon. The spectral indices of the fits, γ = −1.59, −1.53, and −1.68 respectively for ⁴He²⁺, ³He²⁺, and ¹⁶O⁷⁺ are quite similar. The spectra for ⁵⁶Fe¹⁴⁺, on the other hand, are somewhat softer with γ = −3.0. The dependence of γ on Q/M, shown in Fig. 7a, can be approximately fitted as γ ∝ (Q/M)^−0.64.

Fig. 6. Energy-per-nucleon spectra at ωcpt = 48.99. The original spectra of different ion species (a) and spectra of heavy ions shifted for visibility on the y-axis by a factor of 100, 101, 102, and 103 for 4He2+, 3He2+, 16O7+, and 56Fe14+, respectively (b). The dashed lines in the right panel are the power-law fit (Eγ) to the shoulder region of the spectra, where E is the energy-per-nucleon.

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. Figure 7b shows this proxy of the abundance enhancement factor as a function of Q/M at ω_cpt = 48.99 for three energy thresholds, E_t = 10 KE_th, 25 KE_th, and 50 KE_th, well above the initial thermal energy $K{E}_{\text{th}}=0.25m_{\text{p}}{v}_{\text{Ap}}^2$, which is the same for all ion species in our simulations. The value of E_t was chosen such 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 ω_cpt = 48.99. Although the value $E_t=50K{E}_{\text{th}}=12.5m_{\text{p}}{v}_{\text{Ap}}^2$ is close to the cutoff energy of the shoulder region in the energy spectra of ⁴He²⁺ at ω_cpt = 48.99 (see Fig. 5b), the abundances are dominated by the abundances in the shoulder region as the energy spectra falls off very rapidly beyond the shoulder region. We note that the threshold $E_{\text{t}}=12.5m_{\text{p}}{v}_{\text{Ap}}^2$ is in the rapidly falling region of the proton’s energy spectra at ω_cpt = 48.99. It is acceptable, however, as the proton’s energy spectra does not develop a shoulder.

Fig. 7. Spectral index γ of the energy-per-nucleon spectra in the shoulder region vs. Q/M and the fit γ ∝ (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 KEth (b). The straight lines are the fit, abundance ∝(Q/M)−α, for the data and are plotted with the color corresponding to the data.

Figure 7b also shows the power-law fit F ∝ (Q/M)^−α. The enhancement factor drops with increasing Q/M (i.e., heavier ions are preferentially accelerated). The drop in the enhancement factor becomes steeper (the value of α 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 toward the abundances decreases. This decrease in the contribution is more prominent for higher values of Q/M as the shoulder in their energy spectra extends to lower energies. This results in the increase in the value of α 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 energies many times higher 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 pileup regions near the opening of the outflow exhausts. The most efficient acceleration takes place in the flux pileup regions. Heavy ions, depending upon the low value of Q/M which allows them to be non-adiabatic while crossing from inflow to outflow regions, can also be accelerated by pickup mechanisms in outflow regions even before any magnetic flux is piled up. As a result of acceleration, heavy ions develop a shoulder, a nonthermal 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. The abundance enhancement factor, defined as the number of particles above a threshold energy normalized to the total number of particles, scales as (Q/M)^−α, where α 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 ∼1 MeV nucleon⁻¹ has been in situ observed in space (Mason et al. 2002; Bučík et al. 2018). Our simulations show the break in the energy spectra in the energy range 2–6 $m_{\text{p}}{v}_{\text{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_{\text{p}}{v}_{\text{Ap}}^2\approx 0.25$ MeV, and therefore the break in the simulations occurs in the range 0.5–1.5 MeV nucleon⁻¹, consistent with the observations. We note that the energy range of the break depends on the value of v_Ap, estimates of which in the solar corona vary significantly (Brooks et al. 2021). In the observations, the spectral index γ 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 γ before the break is in the same range in our simulations as well. The values of γ after the break are, however, much larger in our simulation than in 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 similarities and differences in the behavior of the two. In our simulations, the power-law index α of the fit F ∝ (Q/M)^−α increases with the threshold energy. The observations also show variations from event to event in the value of α in the range −5 < α < 12 with a mean value of α = 3.26 at 385 keV nucleon⁻¹ (Mason et al. 2004) and 2 < α < 8 with a mean value of α = 3.64 at 3–5 MeV nucleon⁻¹ (Reames et al. 2014b,a). The mean values in the observations are similar, even though their energies are an order of magnitude apart. The mean of the three values of α obtained from power-law fits in Fig. 7b is approximately 3.66. In observations, steeper energy spectra tend to have a steeper fall in the abundance enhancement with Q/M (i.e., a large value of α; Reames et al. 2014b). This effect is similar to the increase in the value of α 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.

Our simulations consider only a single species of heavy ions at a time. In a physical situation, multiple heavy ion species will be present in reconnection regions at the same time and accelerated simultaneously by magnetic reconnection, partitioning the available magnetic energy among the species. This can affect the maximum energy gained by a given species, but is unlikely to affect the preferential acceleration of heavy ions, at least as long as the heavy ions are minor species (in the solar corona heavy ions with Z > 2 compose only 0.1% of all ions while helium is 8.9%; Lodders 2019) so that the electric and magnetic field in the reconnection region are mainly due to the electron and proton dynamics with only a small contribution from the heavy ion species. However, simulations of magnetic reconnection with multiple heavy ion species are required to understand if the minority species can influence the accelerating fields and, if so, to what extent.

It should be noted 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 acceleration alone, but also of the transport of particles from the acceleration site to the detection site. In fact, 3D test particle modeling of the inter-planetary transport of relativistic protons from the source region on the Sun shows that the spectra of the particles is highly observer dependent and do not necessarily reflect the source spectra (Dalla et al. 2020). These transport effects are possible for the spectra of heavy ions as well. However, observations of impulsive events (their intensity time profiles, velocity dispersion, and anisotropies) indicate their nearly scatter-free propagation with scattering mean free paths > 1 au (e.g., Reames 2020, 2022) and therefore a direct comparison with acceleration models is meaningful (Mason et al. 2002 and references therein). The authors point out that while the effects of interplanetary propagation on these events are not large, a definitive resolution requires observations close to the Sun.

We conclude that 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. We note that we did not find any extraordinary 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 and in two dimensions. Magnetic reconnection in three dimensions is host to a number of additional three-dimensional instabilities and allows more complex interactions in the reconnection region. For example, flux ropes formed by the growth of a 3D oblique plasmoid instability in a low-beta reconnecting current layer undergo complex interactions due to the advection and rotation by the reconnection outflow jets, leading to a turbulent state with a stochastic magnetic field and the development of secondary flux ropes in thin current layers that form between the primary flux-ropes (Stanier et al. 2019). A finite guide magnetic field on the order of the reconnecting component and axial structure of flux ropes in 3D magnetic reconnection enable efficient electron acceleration by allowing the particles to leak out and return to acceleration regions, as shown by 3D kinetic simulations (Dahlin et al. 2017). Fully kinetic simulations show that electron acceleration by parallel electric field and curvature driven mechanisms can be enhanced by the filamentation of electric and magnetic fields caused by streaming instabilities in the nonlinear stage of 3D guide-field magnetic reconnection (Muñoz & Büchner 2018). Strong drift-wave fluctuations in the lower-hybrid frequency range, developed in fully kinetic 3D simulations of asymmetric magnetic reconnection (modeling MMS magnetopause diffusion region crossings), give rise to cross-field electron transport and significantly enhanced electron parallel heating in comparison to 2D simulations (Le et al. 2018). It seems likely that 3D plasma instabilities in magnetic reconnection will interact with ³He²⁺ and other heavy ions as well and contribute to their acceleration, as they do for electron acceleration. The details of these complex interactions are unknown and are the subjects for further studies.

Acknowledgments

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. R.B. acknowledges support by NASA grants 80NSSC21K1316 and 80NSSC22K0757.

References

All Figures

	Fig. 1. Out-of-plane current density Jx/(n0eVAp) at four different times.
In the text

	Fig. 2. Time development of the fractional change in average kinetic energy per particle, ⟨K.E.(t)⟩ptcl, from its initial value (a) and average of the square of the normal component of the magnetic field By along the central line of the lower current sheet (b) for different ion species.
In the text

	Fig. 3. Positions of 4He2+ ions represented by dots color-coded according to the ions’ kinetic energy at ωcpt = 19.6 (left column) and ωcpt = 48.99 (right column). The lines with arrows represent magnetic field lines. Only the lower half (y < 0) of the simulation domain is shown.
In the text

	Fig. 4. Positions of 16O7+ ions represented by dots color-coded according to the ions’ kinetic energy at ωcpt = 19.6 (left column) and ωcpt = 48.99 (right column). The lines with arrows represent magnetic field lines. Only the lower half (y < 0) of the simulation domain is shown.
In the text

	Fig. 5. Energy spectra of protons 1H1+ (a), 4He2+ (b), and 16O7+ (c) at different times.
In the text

	Fig. 6. Energy-per-nucleon spectra at ωcpt = 48.99. The original spectra of different ion species (a) and spectra of heavy ions shifted for visibility on the y-axis by a factor of 100, 101, 102, and 103 for 4He2+, 3He2+, 16O7+, and 56Fe14+, respectively (b). The dashed lines in the right panel are the power-law fit (Eγ) to the shoulder region of the spectra, where E is the energy-per-nucleon.
In the text

Fig. 7. Spectral index γ of the energy-per-nucleon spectra in the shoulder region vs. Q/M and the fit γ ∝ (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 KEth (b). The straight lines are the fit, abundance ∝(Q/M)−α, for the data and are plotted with the color corresponding to the data.

In the text