Preferential Acceleration of Heavy Ions in a Spontaneously Fragmenting Flare Current Sheet

We study the ion acceleration in a mesoscale, spontaneously fragmenting flare current sheet (SFCS) characterized by the presence of a plasmoid cascade. The main subject of our investigation is to determine whether and how plasmoid cascades at intermediate scales in a fragmented current sheet of a solar flare can impact the (preferential) acceleration of specific ions. The time evolution of the SFCS is obtained from high-resolution 2.5D MHD simulations. The ion trajectories (in the background fields resulting from the MHD model), energies, and pitch angles are calculated using a relativistic test-particle code based on the half-acceleration–rotation–half-acceleration method. For light ions, the main acceleration effects of electromagnetic fields within the SFCS are analyzed using the guiding center approximation. We identify regions with the most-efficient ion acceleration within the SFCS, the accelerator efficiency, and spectra of the accelerated ions. The influence of the charge-to-mass ratio on ion behavior is also studied and resulting ion abundances are compared with observational data. The main ion acceleration takes place in the regions with a strong polarization term, which is part of the first-order Fermi acceleration. Because the term is mass dependent, heavier ions undergo preferential acceleration. The ion energy spectra, abundance-enhancement factors, and differential fluxes, obtained from the model, exhibit power-law profiles, in agreement with observed solar energetic particle events. Nonetheless, the obtained slopes for the abundance-enhancement factor do not exactly match the observed data. The computed slopes and profiles are not sensitive to changes in the initial plasma temperature.


Introduction
The somewhat mysterious enhancement of the abundances of heavy ions, i.e., those with a small charge-to-mass ratio Q/A (specific charge), in populations of energetic particles w.r.t. the standard composition of the solar corona (or background solar wind) is nowadays a well-confirmed observational fact. Such heavy-ion enrichments are commonly observed in situ during 3 He-rich solar energetic particle (SEP) events (e.g., Mason 2007;Bučík 2020, and references therein). The enhancement factor tends to increase as a power law with the ion Q/A (Mason et al. 2004;Reames & Ng 2004)-the heavy (Ne-Fe) and ultraheavy elements (mass >70 amu) are enhanced by a factor of 3-10 and >100, respectively. However, there is a notable exception to this power-law pattern: the light 3 He isotope, which shows enormous enhancements by factors of up to 10 4 above the coronal abundance (see, e.g., Kocharov & Kocharov 1984;Mason 2007).
There have been several scenarios proposed to explain the mechanisms for the observed anomalous enhancements. First, in general, there are strong observational indications (Reames 1999(Reames , 2013(Reames , 2017) that this acceleration should be attributed to internal processes in solar flares rather than to those related to coronal mass ejection (CME)-driven shocks. However, the flare-related scenario itself offers more possible mechanisms, which have been studied: e.g., acceleration by turbulent reconnection within a plasmoid-rich kinetic current sheet (CS) applied to heavy ions (Drake et al. 2009) or cyclotron resonance acceleration (Roth & Temerin 1997;Miller 1998;Liu et al. 2006;Kumar et al. 2017) involving various wave modes to explain anomalous abundances for both heavy ions and 3 He. In this paper, we attempt to follow the reconnection scenario but at much larger scales, mapping a substantial part of the flare volume and-consequently-a much broader range of anticipated plasmoid cascades with our model.
Magnetic reconnection is-as many observations show (see, e.g., Su et al. 2013, and references therein)-a ubiquitous process in solar flares. In the course of reconnection, free magnetic energy accumulated in the corona is converted to plasma bulk motion and accelerated particles and subsequently also to heat and radiation. It seems that there are more classes of processes exhibiting flaring activity; they likely differ in the geometries of the reconnecting current layers. Here we focus on eruptive solar flares, the current layer of which is expected to form behind the filament ejecta, further creating the CME phenomenon. This idea basically started with the "standard" CSHKP scenario (Carmichael 1964;Sturrock 1966;Hirayama 1974;Kopp & Pneuman 1976); however, it was later modified in order to explain many features observed during flares that are hard to fit into the CSHKP model. The modifications aim basically at sanitizing two main issues of the "standard" eruptive-flare model: (1) its extension from 2D to the full 3D geometry and a proper description of the observed 3D effects (see, e.g., Török & Kliem 2005;Aulanier et al. 2012;Mei et al. 2017), and (2) an explanation of the quite-short timescales of the flares, related to the fast reconnection, signatures of fragmented energy release, and the scale gap between the global flare dimensions (namely, the thickness of the flare current layer) and expected dissipation scale for the fast, collisionless kinetic reconnection. The latter sort of issue was first addressed by the pioneering paper of Shibata & Tanuma (2001), who suggested the conjecture of the so-called "fractal reconnection." It assumes a cascade of tearing-mode instabilities taking place in consecutively smaller CSs formed between separating magnetic islands (also known as plasmoids) in the stretched global flare current layer. This idea was motivated by the somewhat fast appearance of plasmoids in even low-resistive MHD simulations, contrary to theoretical expectations for the tearing mode (Furth et al. 1963). This discrepancy has been theoretically solved by Loureiro et al. (2007), who raised the concept of "plasmoid instability." Later, high-resolution MHD simulations confirmed (see, e.g., Bhattacharjee et al. 2009;Samtaney et al. 2009;Huang & Bhattacharjee 2010;Bárta et al. 2011b;Shen et al. 2011;Ni & Lukin 2018;Zhao & Keppens 2020, and references therein) the tearing cascade suggested by Shibata & Tanuma (2001), and furthermore, they revealed another current-density fragmentation between colliding plasmoids (Bárta et al. 2011b) running concurrently with the plasmoid instability in the extending main current layer. Spontaneous current-layer fragmentation and filamentation seen in high-resolution MHD simulations and its application to the physics of solar flares represent an attractive extension of the standard CSHKP model for multiple reasons: In addition to the reconnection speed-up, it also naturally forms an environment characterized by the presence of multiple small-scale CSs, where fragmented energy release (as indicated by observations) and multisite particle acceleration can easily take place. Therefore, the question of whether the fragmented, plasmoid-rich current layers really appear in the solar flares arises. And, indeed, there are several indications for that: First, the plasma parameters as estimated for the flare current layer fall well into the plasmoidmediated regime in the phase diagram of the magnetic reconnection (Ji & Daughton 2011), and, second, there is much observational evidence for the presence of plasmoids and their cascades in the reconnecting current layers in solar flares (see, e.g., Nishizuka et al. 2009;Karlický et al. 2010;Bárta et al. 2011a;Nishizuka et al. 2015;Shimojo et al. 2017;Karlický et al. 2018).
As cascading fragmentation and plasmoid-rich reconnection become an important part of our picture of flare physics, in this paper we study whether and how it can address the preferential acceleration of heavy ions, too. We are interested in the MHD scales and the impact of the plasmoid-size distribution (expressed as the spectral density of magnetic energy) on the acceleration process and the resulting characteristics of the populations of energetic particles. Namely, we are motivated by the questions of (i) whether the self-similar distribution of plasmoids in a cascade can act as an efficient accelerator of ions up to the observed energies and how this acceleration works, (ii) whether the breaks in magnetic-energy spectra (i.e., the plasmoid size distribution)-the large-scale one, related to the global size (thickness) of the flare current layer, and one on the small-scale end of the cascade, connected with the break toward the dissipation scale-can lead to preferential acceleration, either systematically, with increasing mass-to-charge ratio or selectively, for a specific ion species, and (iii) what the relation is between the spectral index in the magnetic-energy cascade and the energy distributions of accelerated particles. To a large extent, our work can be considered as an extension of the paper by Drake et al. (2009), where a similar topic is approached at very small scales via kinetic simulations, to the intermediate and large scales. The main advantage of our combined approach is rooted in the fact that the high-resolution MHD simulation can cover a larger range of spatial and temporal scales and map a substantial part of the plasmoid-instability energy cascade. Clearly, the MHD simulations themselves are insensitive to the kinetics of particles, from which the plasma is actually composed. On the other hand, full-kinetic plasma models, namely the particle-in-cell (PIC) and Vlasov simulations, can, for obvious technical reasons, model just very small segments of the CS and only for very short time periods. In order to bridge the scale gap and investigate particle acceleration and transport across the large part of the flare volume, the test-particle (TP) approach can be conveniently used. In general, this approach is used quite often, calculating the particle dynamics in the background electromagnetic fields resulting either from the underlying MHD simulations (Karlický & Kosugi 2004;Karlický & Bárta 2006;Gordovskyy et al. 2010aGordovskyy et al. , 2010bGordovskyy et al. , 2014Guo et al. 2010;Gordovskyy & Browning 2011;Zhou et al. 2015Zhou et al. , 2016 or described analytically (e.g., reconnectionlike geometries) (Mori et al. 1998;Dalla & Browning 2006;Grady et al. 2012). This method, when applied to a big ensemble of TPs injected into the flare volume, yields the time evolution of distribution functions of injected particles as its result, which can be used for a thorough investigation of the process of acceleration and a broad range of model parameters.
The TP approach allows us to study the behavior of different species of particles during current-layer evolution in detail: It is not only electron (as has been done in the papers cited above) but also heavy-ion propagation that can be modeled within the CS. One just needs to account for different sizes of the Larmor radius/ periods (w.r.t. the spatial/temporal scales of magnetic structures in the underlying MHD model) and use the appropriate set of algorithms for calculating their trajectories. Hence, the TP method is feasible for investigating the energization of ions with different masses and charges, too, and, by means of statistical processing of their calculated kinetics, eventually also for studies of the possible enrichment of energetic particles by specific species.
In this paper, we study the acceleration of ions of various masses and charges in a spontaneously fragmenting (plasmoidrich) flare CS. A model of background electromagnetic fields is taken from the high-resolution 2.5D MHD simulation by Bárta et al. (2011b). We use a smaller but still substantial subvolume of the original MHD simulation domain that covers basically the entire flare. Individual ions are tracked using Vay's relativistic algorithm (Vay 2008). Their full trajectories are then statistically processed in order to obtain substantial collective characteristics of the ion populations like energy spectra or pitch-angle distributions.
This paper is organized as follows: First, we briefly describe the underlying MHD model (Section 2.1) and algorithm for ion-trajectory calculation (Section 2.2). In Section 3, we analyze the model results and impact of various parameters on the ion distribution functions. The time evolution of the energy spectra and pitch-angle distributions for ions is studied in Section 3.1 in order to see how efficient and fast the plasmoid cascade can be as an accelerator. Our main goal, investigating the dependence of acceleration efficiency on the specific charge, is presented in Section 3.2. We attempt to understand the obtained results and identify the main physical processes and regions, where the processes take place within the fragmenting CS-such analysis is presented in Sections 3.3 and 3.4. In Section 4 we aim to compare characteristics of accelerated ion populations inferred from our modeling with those found by observations. Because there is no one-to-one correspondence between the quantities diagnosed in the acceleration source (our model) and in-situ measurements in the solar wind (i.e., far from the acceleration site), we try to find suitable proxies for this comparison. Finally, the results of our modeling are discussed w.r.t. findings of the previous works on this topic and from the point of view of comparison with the observed data.

Model Description
Our model is composed of two parts: (i) an MHD simulation, which describes the system evolution on large and intermediate scales and provides the background electromagnetic fields, and (ii) calculations of particle trajectories in those fields and their processing in the statistical sense. We now describe both parts separately.

MHD Model of a Spontaneously Fragmenting CS
In order to study the possible impact of a spontaneously fragmenting, plasmoid-rich, flare CS on the (specific) ion acceleration, we use the results of the MHD model described in Bárta et al. (2011b), who address exactly this topic. Here we summarize just the main features of the model and emphasize its aspects that are important for the subsequent TP tracking. We refer to the above paper for details.
The MHD simulation starts from a vertical (modified) Harris-type CS in the fully ionized hydrogen solar atmosphere with a gravitationally stratified plasma density (see Figure 3 in Bárta et al. 2011b). The main initial magnetic field component is in the z direction, oriented antisymmetrically along both sides of the CS. Its asymptotic value far from the CS (x → ∞) is B z0 . The current-density vector thus points in the y direction. The model also considers a weaker guide field in the same direction B y0 = 0.2B z0 . The importance of the guide field for particle acceleration in the CS was discussed by Hamilton et al. (2003Hamilton et al. ( , 2005. They showed that the guiding field tends to keep accelerated particles in the CS plane and prevents particles from drifting out of regions of high-resistive electric fields. Moreover, Bárta et al. (2011b) have shown that the guide-field magnitude also influences the rate of spontaneous fragmentation of the CS.
The evolution of magnetized plasmas has been modeled by a set of compressible one-fluid MHD equations (Priest 1984, p. 91) using the finite volume method. The standard collisional (Spitzer) resistivity is set to zero. On the other hand, the authors attempted to mimic an emergence of the other terms in the generalized Ohm's law (that arise from the Hall, two-fluid, and kinetic physics) at small (subgrid) scales by means of generalized (effective) resistivity η (Kliem et al. 2000) given by the prescription Here, v d = j/en e is the current-carrier drift velocity, v cr is the appropriately chosen value of the critical drift velocity (characterizing the high amount of current-density filamentation), and C is a proper scaling factor. The MHD equations are solved in dimensionless units with C = 0.003 and v cr = 15.0. As already written, this setup ensures the onset of generalized dissipation (E · j ≠ 0) at small scales; the scale at which this happens is controlled by the value of v cr .
In order to reach high resolution and in this way cover a broad range of scales, the simulation is restricted to 2.5D, meaning 2D geometry (with the invariant direction along the yaxis) and 3D vectors of plasma velocities and fields. The MHD model is essentially scale free, which is reflected by the dimensionless units used in the simulation. However, in order to apply it to the vertical current layer of the solar flare, the authors have used the following parameters (based on values inferred from flare observations): The plasma density (equal to the electron density n e ) at the bottom of the simulation domain is n e = 2.4 × 10 10 cm −3 (for x → ∞), the temperature T = 2 MK, and the gravity scale height H n = 120 Mm. The asymptotic value of the vertical magnetic field is B z0 = 40 G. This setup corresponds to a low-β plasma (β = 0.1), so the evolution of the CS is controlled by the magnetic field. The initial width of the CS is 1.2 Mm. The total size of the computational domain is 43.17 Mm along the x-axis and 172.8 Mm in the vertical z direction. The 2D domain is covered by 1600 × 6400 grid points forming a coarse square grid with a spatial resolution of 27 km. During the MHD simulation, the grid has been adaptively refined by a specific AMR (adaptive mesh refinement) technique: Whenever it was appropriate (because of continuing current-density filamentation), an embedded smaller MHD subsystem with a finer grid (and its own finer time-stepping) was allocated and independently evolved, interacting with the global mesh via dynamic (interpolated in space and time) boundary conditions. This setup somewhat anticipated (though with a far more simple implementation) the modern DISPATCH approach (Nordlund et al. 2018), with its convenient local time-stepping. The embedded subsystem cells are 10 times smaller than the coarse ones, giving an effective spatial resolution of 2.7 km. The entire simulation covered ≈300 s in real time, corresponding well to the duration of the impulsive phase of the flare.
Because of the very large difference in dynamic scales between the flare evolution and the motions of the energetic particles, we selected only a small subset of the original MHD results, namely a short time interval 320-330τ A at its later evolution stage, i.e., 235.8-243.2 s. This time interval is equivalent to ≈2 × 10 6 proton gyroperiods. The MHD model used as a background for our subsequent TP simulations, shown in Figure 1, thus represents already quite a developed fragmented CS, which contains a large number of plasmoids (at many scales) and several resistive regions (dissipative CSs). We use the entire MHD domain for our subsequent particle trajectory calculations. Nevertheless, for the later analysis of acceleration processes (Section 3.3), we restricted ourselves also in space-just to the region marked by the white box in the right panel of Figure 1. The underlying MHD simulation is 2.5D, with an invariant direction along the y-axis.
The plasmoid cascade, which represents an essential attribute of the spontaneous current-layer fragmentation, can be characterized by several means. As an alternative to the plasmoidsize distribution P(L) (e.g., Drake et al. 2009), Bárta et al. (2011b) calculated the spectral energy density of the magnetic field in the plasmoid by means of a Fourier transform of the B x ( i.e., perpendicular to the CS) component of the magnetic field along the CS. It is displayed in Figure 2 (see also Figure 6(c) in Bárta et al. 2011b). Quite a broad inertial range characterized by the power-law spectrum is visible, finished by a break in dissipative scale. This behavior is in qualitative agreement with the concept of fractal or-better yet-cascading reconnection (Shibata & Tanuma 2001), i.e., the fragmented CS manifests self-similar properties in the inertial range of scales. We will return to the relation between the plasmoid-size distribution and the spectral density of magnetic energy in plasmoids in the Discussion and Conclusions section (Section 5).
As already mentioned, the MHD simulations are to a large extent scale free. However, their application to real systems bounds the model with actual dimensions. At the large-scale end, it is the size of the simulated system, i.e., the flare dimension in our case, represented here by the thickness of the initial current layer (plus its reasonably large vicinity). However, there is also a small-scale threshold related to dissipation. In order to keep the dissipation scale resolved, Bárta et al. (2011b) chose the v cr parameter, which controls the onset of the effective resistivity in such a way that the thickness of the dissipative CSs does not drop below a couple of model cells. The finite thickness of the dissipative CS consequently prescribes the size of the smallest plasmoids in the cascade. Actually, this fact is exactly the reason for the small-scale knee of the plasmoid spectra displayed in Figure 2. The MHD simulations show that typically the size of the small plasmoids along the current layer is an order of magnitude larger than the thickness of the CSs that interleave the plasmoids. In Figure 2, the plasmoid scaling breaks at ≈200-300 km, corresponding to the dissipative-CS width of ≈20-30 km (i.e., indeed, a couple of MHD-model mesh cells of size ≈3 km). As adduced, this was accomplished with the appropriate choice of v cr . However, in the real world, the dissipative-CS width is limited by the onset of fast kinetic reconnection at small scales, and its typical value is around the ion inertial length d i ≡ c/ω pi (Büchner 2006). For the parameters of the flaring corona, this means (by order of magnitude) a dissipative-CS thickness of ≈10 m. To sum up: In order to match the dissipation scale (scaling break for the plasmoid distribution in Figure 2) in the MHD model with its real counterpart, we should shift the entire scaling The red rectangle in the left panel represents the area of particle generation. The white rectangle in the right panel represents the area of interest for ion-trajectory analysis. Note that only the most dynamic subvolume of the entire MHD simulation is displayed; the scale at the axes is highly anisotropic for this reason (actually, the plasmoid lengths are typically an order of magnitude larger than their widths). law in Figure 2 to the right by a factor of ≈3000. To that end, we apply a downscaling factor of 1:3000 to the spatial coordinates (meshes) of all fields resulting from the MHD model. After applying this downscaling, the actual dimensions of our box are 14.4 × 57.6 km. The underlying MHD simulation is 2.5D, with an invariant direction along the y-axis. For the subsequent full 3D TP calculations, we thus replicated the 2D fields in the third dimension (y) to form a box with a finite depth of 57.6 km (i.e., the same value as its height).
The question of how this shift toward small scales impacts the typical values and structure of the resulting fields arises: To what extent do they still correspond to the unresolved physics at the MHD subgrid scales? In the case of density and temperature, the downscaling should not change their typical values; however, it can somewhat affect their structures. Namely, the spatial variations in the maps of plasma density and temperature-if simply rescaled-do not correspond pretty well to reality. First, the gravity height scale no longer corresponds to the actual stratification of the (downscaled) solar atmosphere. Moreover, the density and temperature structure, which at large scales can be reasonably considered to be the result of adiabatic compression/ expansion in the course of MHD evolution, is likely smeared out by heat conduction and similar effects at the downscaled small dimensions. Because the density and temperature affect just the initial statistical distribution of the TPs in our subsequent TP calculations, we use uniform temperature and density, instead of those resulting directly from the downscaled MHD model, for the sake of simplicity. In our model, temperature and density match the temperature and density of the quiet corona-T = 1 MK, n e = 10 15 m −3 . This simpler setup also allows us to evaluate the impact of different initial temperatures on the resulting characteristics of accelerated ion populations -see the discussion in Section 5.
For the magnetic field, the situation is more complicated. Bárta et al. (2011b) have shown that plasmoid coalescence leads to the formation of thin transverse CSs between colliding plasmoids, which become subsequently prone to next-level fragmentation and reconnection. Because the plasmoids are mutually attracted by Lorentz force, the reconnection between them is essentially driven and hence very efficient. It is likely that reconnection in these secondary CSs is thus the main contributor to particle energization. Such a process is therefore in our highest interest. The formation of those transverse CSs is connected to a significant flux pileup, which continues until the threshold for plasmoid instability (CS aspect ratio 100; see Loureiro et al. 2007) is reached. This flux pileup causes a local enhancement of the magnetic field strength, we estimate the factor of its increase in those small-scale CSs to be ≈10. Hence, investigating those small-scale structures we should expect an order of magnitude higher magnetic field than the ambient field of 40 G used for the original large-scale MHD simulation. Therefore, we estimate the maximum strength of the magnetic fields at small scales to be 400 G and its mean level-for further estimations of typical gyroradii (see Table 4 in Section 3.3)-to be 200 G, i.e., twice that in Drake et al. (2009). Considering this local amplification of the magnetic field together with basically the preserved typical value of density, the Alfvén velocity is upscaled with B. So, the combined scaling, i.e., 1:30,000 applies for the timescales. The evolution time of our box after complete rescaling is thus 2.43 × 10 −4 s. Last but not least, because of large gyroradii and low mobility of the heavy ions, the probability of their significant interaction with the DC electric field inside the resistive regions is negligible. For this reason, we ignore the effects of resistive electric DC fields on ions and concentrate only on the "drift" acceleration, similarly to Zhou et al. (2015).

Particle Tracking
The gyroradii of most heavy-ion species studied in our simulation are comparable to or even larger than the smallest magnetic structures (plasmoids) in our (downscaled) background MHD model (see the discussion in Section 3.3). It is, therefore, generally speaking, inappropriate to use the computation-saving (and especially in electron TP simulations, popular) guiding center approximation (GCA), and we have to follow the dynamics of ions by direct integration of their equation of motion. Only the trajectories of light ions are to some extent eligible to be calculated by the GCA method, and we will take advantage of that fact in our further analysis of acceleration mechanisms in Section 3.3.
In view of what has been written above, the evolution of our TPs in the complex electromagnetic fields of a fragmented CS is calculated using Vay's relativistic TP code based on a variant of the half-acceleration-rotation-half-acceleration (HARHA) algorithm (Vay 2008): where m and q are the rest mass and charge of the particle, r and v are the particle position and velocity with u = γv, g = -= + v c u c 1 1 1 2 2 2 2 is the relativistic Lorentz factor, B and E are the magnetic induction and electric In order to get the field values at the positions and times of occurrence of individual particles, we use linear interpolation in space and time (Gordovskyy et al. 2010b). During the interpolation, the convective electric field is kept perpendicular to magnetic induction and the resistive electric field is set to zero as we are ignoring the DC field acceleration effects for the reasons explained in Section 2.1.
The integration time step is set to a fraction of a particle gyroperiod, and it is limited also by the size of the MHD cell in which the particular particle is located (on the border between fine and coarse cells, the size of the finer cell is used).
The TP simulation starts by uniformly distributing approximately 5 × 10 6 TPs (ions) in a belt centered on the z-axis of the CS with thickness 1 km, covering the regions of plasmoid occurrence-see the red rectangle in Figure 1. The directions of the initial particle velocity vectors in space are randomized in both pitch angle θ and azimuthal angle f of spherical coordinates. In this way we obtain an isotropic distribution of velocities in 3D space, i.e., a homogeneous distribution of unit velocity vector endpoints on a unit sphere. This initial distribution transforms to a uniform distribution in directional cosines (see the bottom-left panel in Figure 3). The initial energy distribution corresponds to a shifted Maxwellian with the temperature 1 MK and a displacement corresponding to the bulk motion of surrounding plasma (v = v thermal + v plasma ). The shift ensures that the velocity averaged over the gyroperiodv has its average magnitude corresponding to the Maxwellian velocity. This is apparent from the definition of the magnetic moment μ (Northrop 1963), 2 is a relativistic correction term. The equation above expresses the conservation of the particle magnetic momentum in the electromagnetic field. Because the velocity of the particle changes during a single orbit, we will use the averaged counterpart of the velocityv to obtain both initial and final energy spectra of particles.
The particle gyroperiods T o and gyroradii R o are defined using the relations The simulations are performed within the downscaled domain described in Section 2.1. To study how the mass and charge affect the behavior of the individual ions in our CS, we have selected in total six different masses and four charges of ions for our simulations. The particular values are summarized in Table 1.

Results and Their Analysis
Using Vay's full-orbit integration, we calculate the trajectories, energies, and pitch angles of several ensembles of ions in the temporally evolving electromagnetic field of a spontaneously fragmenting CS. The ensembles differ in masses A and specific charges Q/A of ions (see Table 1), while particles within each ensemble are identical. Then, by means of statistical methods, we process and analyze the simulation results to derive the time evolution of energy spectra, final energy spectra, and distribution functions for ensembles of ions with various specific charges and trajectories of ions in fields of the CS. In order to obtain insight into the processes related to ion energization and preferential acceleration of ions of various kinds, the GCA method is used.

Time Evolution
First, we study the efficiency and rate of ion acceleration in the plasmoid-rich CS. This simulation is done only for ions with charge Q = 1 (in units of the elementary charge e) and mass A = 1 (in the atomic mass units-amu). We capture the time evolution of the ion ensemble in several snapshots of the energy spectrum, pitch-angle distribution, and the distribution of energy in pitch angles. The latter distribution allows us to assess the efficiency of ion acceleration in parallel and perpendicular directions. The results are summarized in Figure 3.
The time evolution of the energy spectrum is shown in the top-left panel in Figure 3. It demonstrates the gradual energization of particles in the ensemble and the formation of high-energy tails that can be fitted with a power law, where E is the ion kinetic energy, N is a number of particles, and δ the power-law index.
For electrons, such behavior is well known from observations of the hard-X-ray nonthermal bremsstrahlung generated during impulsive phases of solar flares. The typical power-law indices for electrons accelerated in solar flares obtained from hard-X-ray observations and a thick-target model (Brown 1971) are δ ∼ 2.5-3 (Aschwanden 2005, p. 603). Observations of 3 He-rich SEPs show that half of the energetic events have either rounded spectral profiles or single or double power-law profiles (Mason et al. 2002). Miller (1998) has shown that rounded profiles can be well fitted using the cascading MHD turbulence model. However, we postpone a more detailed discussion on the relations of our results to observations to Section 4 below because the observed energetic spectra are slightly different from the particle energy distribution functions used in this section.
We performed a least-squares power-law fit to the tail of the energy spectra (top-left panel in Figure 3) and found that the power-law indices δ E evolve with time from δ E = 1.2 during the first quarter of simulation, to δ E = 0.9 at the end of the simulation. The values for the selected snapshots are summarized in Table 2. The change of slope with time is relatively small, but the whole profile gradually shifts toward higher energies. Another feature that can be seen in the spectra is the presence of a hump. We speculate that such a result could mean a beam of accelerated particles has formed in our simulation. Because traditionally derived quantities are used to describe SEP events, we postpone comparison with observations to the next section.
The gradual energization of the particle ensemble is also visible in the top-right panel of Figure 3, where the average energy gain per particle as a function of time is plotted. Here, it is also apparent that the energy gain at the end of the simulation is not saturated, so the particles could gain even more energy if the time evolution was longer.
The bottom-left panel in Figure 3 shows that the initial isotropic distribution of particle pitch-angles changes with time, and two gradually growing populations of particles oriented in the parallel and antiparallel directions relative to the magnetic field lines arise. The bottom-right panel in Figure 3 shows that the two populations are accelerated very efficiently. The alignment of pitch angles with the magnetic field also enables particles to escape the magnetic confinement and propagate to outer space, thus allowing the formation of a particle beam.
Due to diagnostic limitations of the HARHA integrator, which is unable to separate the motion of particles into more meaningful components, we postpone the identification of the main acceleration effects of the complex electromagnetic fields in the fragmenting CS until Section 3.3. In that section, we analyze the mean forces leading to the acceleration of particles using the GCA method. In Section 3.4 we justify the use of such a GCA approach by analyzing its results in detailed comparison with the full-orbit particle calculations using HARHA for light-and intermediate-mass ions.

Role of Specific Charge
In this section we address the main topic of our study-the not yet fully understood phenomenon of preferential acceleration of heavy ions during 3 He-rich SEP events. To study how the specific charge Q/A influences the behavior of the   Figure 3 t (s) δ E 5.0 × 10 −5 1.2 1.0 × 10 −4 1.0 1.5 × 10 −4 0.9 2.0 × 10 −4 0.9 individual ion ensembles in the CS, we use similar distribution functions to those in the previous section, but all the results are taken at the end of the simulation (t = 2 × 10 −4 s). No time evolution is thus captured, and the results for ensembles with various specific charges at identical moments of evolution are compared to each other. By keeping the charge Q = 1 and varying the mass only, we get six selected specific charges 1, 1/5, 1/10, 1/15, 1/50, and 1/200 covering a wide range of ion atomic masses with A = 1-200 (see the first row of Table 1). The heavier ions that have a smaller specific charge show progressively larger energy gain per particle, which leads to preferential acceleration-see the top-right panel of Figure 4. The average energy gain per specific charge forms a power-law profile with δ G = 1.2. The effect is visible also in the spectra displayed in the top-left panel of Figure 4, which are more energetic for heavier ions.
Because the specific charge can be changed not only as a consequence of varying mass but also charge, we have also studied combinations with Q > 1. As an illustration, the energy spectra for various mass-charge combinations are shown in Figures 4 and 5, and the power-law indices of their high-energy tails are summarized in Table 3 with respect to the particular Q and A. For spectra, where possible (low-mass ions), the powerlaw fits of the high-energy tails give practically identical values of δ E between 0.8 and 0.9. Again, the heavier ions show more pronounced high-energy tails, with similar spectra for ions with comparable specific charges.
The pitch-angle distribution changes in a similar fashion for all combinations of Q and A, with the exception of mass A = 200 (the bottom-left panel of Figure 4). The ensembles corresponding to low-mass ions show an increase in the number of particles oriented in the parallel or antiparallel direction to the magnetic field lines. Similar behavior can be seen in the bottom-left panel in Figure 3. On the other hand, the ultraheavy ion A = 200 shows a complete reversal, with a strong population of particles oriented in the direction perpendicular to the magnetic field.
With the growing mass of particles, the pitch-angle distribution of energy also changes and exhibits a growing asymmetry relative to θ = π/2 ( q = cos 0)-see the bottom-right panel of Figure 4. We speculate that these asymmetries are closely linked with increasing ion Larmor radii (see Table 4) and the larger presence of inhomogeneities (local changes) in the field during a single particle orbit, which can result in slightly enhanced acceleration or deceleration of particles due to nonadiabatic behavior. The probable reason for pitch-angle asymmetries is investigated further in Section 3.4.
Generally, for A 50 the particles are more influenced by local changes in the field than by magnetic drifts, which can be attributed to the very large Larmor radii and gyroperiods of these heavy ions (see Table 4) with respect to the dimensions of our simulation domain and evolution time. In this sense, we are reaching the physical limits of the simulation. The final states of their distribution functions are basically random, and utterly different profiles can be obtained for just slightly different background-field configurations. Therefore, the results for heavy ions with A 50 are not statistically significant because their proper tracking requires the inclusion of larger spatial scales and longer simulation time to provide sufficiently averaged data over a gyroperiod.

Identification of Accelerating Mechanisms by GCA Analysis
Although full-orbit integration is suitable for the precise tracking of ions in general electromagnetic fields, the GCA allows us to separate the various field characteristics (e.g., field curvature, polarization, etc. described by the individual GCA terms) and their effect on the particle orbit and energy (see Equation (8)). In this way, we can get a clue about the significance of individual effects playing a role in particle energization.
The GCA decomposes the parallel and perpendicular motions of the particle guiding center relative to the magnetic field lines and assumes the conservation of particle magnetic momentum μ (Equation (4)). For future reference, we present only the parallel component of the guiding center equations (Northrop 1963;Gordovskyy et al. 2010b), where b is the unit vector corresponding to B and (( ) 2 the relativistic Lorentz factor. The first term on the right-hand side is acceleration by a resistive field, the second term is the magnetic mirror effect, and the third term the inertial/polarization effect. The GCA approach enables us to study solely the influence of magnetic drifts on the particles and estimate their importance for acceleration efficiency.
The GCA approach is valid only for light ions, where the gyroperiod is faster than significant changes in the underlying magnetic field and the gyroradius smaller than typical structures (i.e., plasmoids) in the field. Our data set contains a wide range of scales: the smallest scale consists of plasmoids with widthL 10 min m, whereas the largest plasmoids have widthL 3 max km. Because the smallest plasmoids do not occupy a significant part of the simulation volume (see Figure 1) and the chance to be hit by ions is thus small, we estimate that GCA will be-at least to some extent-applicable for particles with Larmor radii R o < L mean ∼ 100 m and gyroperiod T o < τ mean ∼ 5 × 10 −5 s-see Table 4, which has an overview of the ion parameters used in the simulation. The same domain of GCA applicability can be guessed from the top-right panel of Figure 4, where the curves obtained by (exact) HARHA and perturbation-theory (GCA) start diverging significantly for ions with A > 15. Note, however, that the power-law profile is still present for the GCA approach with δ G = 0.8.
The energy gain for GCA can be computed using (North- The first term is the electrical acceleration, the second term is the induction effect (betatron-type acceleration), the third term is the mirror/gradient B effect, the fourth term is the inertial effect (which includes centripetal acceleration), and the last is the polarization effect. Northrop (1963) has also shown that the induction and mirror effects are responsible for the first-order type A Fermi acceleration-an acceleration by a moving magnetic mirror along a straight line of force (hereafter Fermi Ia)-whereas the inertial and polarization effects are responsible for type B Fermi I acceleration-an acceleration within a curved magnetic field (hereafter Fermi Ib), with types A and B in the notation of the original paper of Fermi (1949).
We would like to remind here that we set E || = 0 and that we use identical initial energy distributions for all ions-terms witĥ mv 2 and || mv 2 thus have comparable sizes for ions with different masses. This is not the case for magnetic moments scaled by mass mμ, where we need to account for its field component v E . Because the E-cross-B drift is not dependent on mass and charge, and thus, it is identical for all ions, and because its mean magnitude in plasmoids (v E ∼ 4 × 10 6 m s −1 -see Figure 1) is larger than typical ion thermal velocities (v thermal 1.3 × 10 5 m s −1 for 1 MK -typical proton thermal velocity), terms with 2-5 in Equation (8) will be stronger for heavier ions, thus leading to mass-dependent acceleration within a strong inhomogeneous EM field.
To distinguish which component plays the most important role during acceleration in our particular CS, we monitored individual components from the gain equation as the time average within GCA simulations. From Figure 6, we see that the most important term during acceleration is the polarization effect followed by the inertial effect, both of which are part of the Fermi Ib acceleration, which occurs in curved magnetic fields. From Equation (7) and the relation it is also apparent that both effects are responsible for the increase in the number of ions with parallel or antiparallel orientation to the magnetic field over time and their efficient acceleration (see the graphs for light ions in the bottom-left and right panels of Figure 4). As we can see from the top-right panel of the same figure, the GCA solution diverges from the HARHA solution with increasing ion mass; therefore, for heavier ions, we have to account also for nonadiabatic acceleration during a single gyroorbit.
If we compare this behavior to electrons, the main mechanism for electron acceleration in a resistive-free electromagnetic field is the curvature/centripetal acceleration term, which is part of the inertial term, and the mirror/gradient B effect (Zhou et al. 2015); therefore, electrons are accelerated by both type A and B Fermi I acceleration.

Trajectory Analysis
To confirm our findings, we visualized trajectories of randomly selected individual particles for ions with specific charge Q/ A = 1 and Q = 1, calculated by the full-orbit HARHA approach, and compared regions where particle acceleration takes place with the polarization term.
The left panel of Figure 7 shows particle propagation within the CS-see Figure 1 for the precise location in the CS marked by the white rectangle. The trajectories are projected onto the plane, disregarding any out-of-plane motion in the process. The locations where the particles are accelerated or decelerated indeed match the areas with the strongest polarization term. The disproportion in area of regions with a different sign of polarization term is also likely responsible for the asymmetry in the energy distribution per pitch angle (see Figure 4). Arguably, the inertial term can also take a significant role in such asymmetry because it is scaled by the parallel velocity of the guiding center, and we also cannot entirely eliminate the possibility that asymmetry is relevant only to the current phase of MHD evolution or that it can be reverted in ongoing flare evolution. The concrete mechanism is at this moment unknown.
Notice also in the middle panel of Figure 7 that the E-cross-B drift and subsequently the magnetic drifts are pushing particles inward, toward the plasmoids. The lengths of the trajectories can also be taken as a proxy of the acceleration rate. The longest lines are concentrated within the plasmoids, where the acceleration is most efficient.
The right panel of Figure 7 shows the particle energy gain or loss gathered on their trajectories. It can be clearly seen that the strongest acceleration occurs at the maxima of the polarization term.

Comparison with Observational Data
To compare results with the available observational data, we have to first match qualitatively our produced data to measured quantities.

Differential Fluxes
The populations of energetic ions in in-situ measurements are characterized by different quantities from the energy distribution functions such as those shown, e.g., in Figure 4. A typical characteristic is the differential flux, the particle flux per energy and solid angle. The ion events have shown power-law or double power-law spectra profiles in energy per nucleon with slopes δ F ≈ 1-3 for energies below 1 MeV nucleon −1 (Mason et al. 2002). Spectra below 1 MeV nucleon −1 are harder (show a smaller spectral index; e.g., Mason et al. 2002;Bučík et al. 2018). Furthermore, in some events, 3 He and Fe show rounded spectra to lower energies (Mason et al. 2002).
Differential flux is defined as the differential of particle flux (first-order moment of velocity distribution function) per energy and solid angle (Wüest et al. 2007, p. 183): , with energy expressed in energy per nucleon (E/A), which is a measure of ion velocity. Both energy and differential flux are expressed in a logarithmic grid. For logarithmic bins and energy measured in eV/nucleon, the energy E is defined as where x is between 0 and 1. Particle flux can be then rewritten from where we used the equality where n is plasma density. To obtain the differential flux, dx is expressed as dE:

Q [e]
A (amu) Note. The values were computed using B = 2 × 10 −2 T, T = 1 MK, v E = 4 × 10 6 m s −1 , and¯= . The resulting Figure 8 and Table 5 show the ion differential fluxes derived from our model. They exhibit power-law profiles, likewise the observational data, with slopes ranging from δ J = 0.8 to δ J = 1.4, with the average δ J = 1.2. The slopes are calculated from the tail of the flux profile in the energy range 3 × 10 3 −2 × 10 5 eV nucleon −1 . The values of the power-law indices obtained from the simulations exhibit a slight scatter, but the δ J produced by the model belongs in the range of typically observed values δ J ≈ 1-3.

Abundance-enhancement Factor
Another experimental quantity, measured in situ, is the enhancement factor defined for ions with a given specific charge Q/A. The experimental Q/A values cover the range from 1 to ∼0.08, where 1 corresponds to H and 0.08 to partially ionized ultraheavy elements. Q/A = 2/3 (1/2) corresponds to fully ionized 3 He ( 4 He) to name a few examples. The source temperature used for the calculation of charge Q for a specific charge Q/A is higher than that in our simulations, T ∼ 3.2 Mk (see also the discussion on initial temperature influence in Section 5). Its magnitude has been determined from in-situ elemental abundance ( Bučík et al. 2021). The energy of ions is in the experiment measured within the bands defined usually as D = E E 2 0 , which we use also in our study (see Table 6 for a detailed band definition).
The enhancement factor K in observational data is usually expressed as an abundance of element X relative to oxygen O (X/O) in the in-situ SEP measurement, divided by the X/O in the corona: For 3 He-rich flares, a power-law dependency of the enhancement factor on the specific charge has been found with typical slopes δ K = 3.26 (Mason et al. 2004) or δ K = 3.64 (Reames et al. 2014b).
We use two proxies of the enhancement factor: (A) the nonnormalized number of accelerated particles within a given energy band as defined in Table 6, and (B) the ratio of the number of accelerated particles above the selected energy threshold to the total number of particles. When plotting those proxies as a function of specific charges Q/A, we use the same approach as in the observed data for easier comparison; our values of Q/A correspond to individual elements detected by experimental apparatuses.
The first proxy (A) is based on the fact that our simulation uses an identical number of particles for a given element; therefore, the initial abundances are identical for each ion N (X)/N(O) = 1. As we can see from Figure 9, the first proxy for the enhancement factor versus specific charge changes with  Figure 1). Trajectories are calculated using the HARHA approach. Left panel: ion trajectories color-coded by the energy gain (red) or loss (blue) during the whole simulation. The underlying field is color-coded by the magnitude of the polarization acceleration term (see the scale at the top and Equation (8)). Middle panel: ion trajectories color-coded by distance (see the color scale at the top). As the trajectories are 3D projections onto the plane, their lengths in the figure do not correspond to their real lengths in space. Right panel: ion trajectories color-coded by absolute energy gain (difference between final particle energy and initial particle energy).
increasing energy in a particular band from almost a flat linear curve (energy band I) to a curve with a partial power-law like dependency (energy band IV) with a final slope δ K(A) = 0.8 and a gap that separates heavier and lighter ions with Q/A > 0.8.
Our second proxy for the enhancement factor (B) uses the ratio of the number of accelerated particles above the selected energy threshold E t to the total number of particles, where E t = 3.0 × 10 4 eV is selected from the tail of the initial Maxwell distribution with T = 1 MK and set to match the definition of the lowest energy band. The threshold is in this case not dependent on the number of nucleons as all ions have identical initial energy distributions. This ratio thus gives the acceleration efficiency of the coronal accelerator for a given element.
From the top-left panel of Figure 10 we can see that ions are subjected to preferential acceleration, with ions with smaller specific charges (heavier ions) accelerated more effectively than ions with larger specific charges (lighter ions). Proxy B again has a power-law profile with mean slope δ K(B) = 0.5 ± 0.1, which is relatively close to that of proxy A, but contrary to proxy A, no separation into two distinct parts is observed, which appeared in energy band IV for proxy A.
In addition to proxies A and B, which should mimic the observed enhancement factor, we have performed further statistical analyses with the aim to illustrate the dependence of acceleration efficiency on the specific charge; they are presented in Figure 10. The average energy gain per particle ΔE for various specific charges is presented in the top-right panel, which also shows a power-law profile with mean slope δ G = 1.1 ± 0.1. The bottom-left panel of the figure shows the average energy gain from the previous situation but scaled by the number of nucleons. Here no power-law dependency is observed.
If we compare the slopes of our proxies to the enhancement factor (A and B) with the value observed in impulsive events, where the slope is around 3.26 (Mason et al. 2004) or 3.64 (Reames et al. 2014b), the slopes do not match. The disagreement of slopes can have several causes, which we discuss in the Discussion and Conclusions section (Section 5).  Table 5 for the power-law indices. Table 5 Power-law Indices δ J Corresponding to Ion Differential Fluxes for Selected Charges Q and A at the End of the Simulation (see Figure 8)

Discussion and Conclusions
We have studied the impact of MHD plasmoid cascades at mesoscales on the preferential acceleration of ions. Such cascades are, for good reasons, anticipated in fragmenting current layers in solar flares, from which the energetic ions are expected to originate. The effect of the plasmoid pool on acceleration efficiency was studied earlier at kinetic scales by Drake et al. (2009). We extended that study to larger dimensions, covering a much broader range of scales, making it possible to map the cascade and its multiscale dynamics.
First, we aimed to answer the question of whether the actual presence of a multiscale pool of plasmoids (obeying a power-law scaling) in the simulation of a fragmented CS in the larger computational domain can lead to super-Alfvénic particles exhibiting power-law spectra. In some sense, it can be considered an extension of previous multiplasmoid scenario PIC studies targeting either a different context of the reconnection exhaust (Wu et al. 1998;Drake et al. 2009) and/or lacking the truly multiscale distribution of the plasmoids present in the simulation (Drake et al. 2006(Drake et al. , 2010. Here, the answer is positive: The simulations produce a weak population of ions with slightly super-Alfvénic velocities (≈1.3v A ), reaching energies of ≈0.5 MeV nucleon −1 and exhibiting a powerlaw spectrum. Taking into account the dimension and number of plasmoids in our computational domain in contrast to the total volume and number of plasmoids in CSs of real solar flares, we could expect more numerous energetic particles and even higher energies owing to multiple interactions of particles with plasmoids on their path throughout the CS.
We have analyzed the mechanisms for this quite efficient acceleration and found that it is basically a Fermi I type of process (or its generalization for heavy ions, where the GCA approach is not strictly applicable; see Section 3.3) acting in the stochastic environment of the multiscale pool of moving plasmoids, which leads to the fast energization. In our case, the ions tend to be accelerated while drifting through a single plasmoid (or visiting two); see the trajectories in Figure 7. The stochastization of the particle dynamics is thus not reached by the repeating back-andforth bouncing motions in plasmoids (as in, e.g., Drake et al. 2006Drake et al. , 2010. Nevertheless, the stochasticity (manifested by the power-law energetic spectra) is still present due to various plasmoid sizes (obeying a power-law scaling) and randomized initial particle positions and velocities. Because of the finite time spent in the acceleration region, the dynamics of ions in such an environment is strongly dependent on the initial conditions. The plasmoid-scaling law, together with the finite acceleration time (limited by the particle presence in the region), thus results in the formation of power-law spectra. The relation between the powerlaw scaling of plasmoids and the power-law energy spectra of ions is thus not related exclusively to the ion-pickup mechanism described in Drake et al. (2009, namely relation (9)), but it likely represents a more generic behavior. In our context, the Fermi I acceleration is responsible for the ion acceleration. This is the reason why the spectral indices of the plasmoid-size distribution (δ P ) and the energy spectra (δ K ) are apparently not bound by the relation δ K = δ P − 3 (see Equation (9) in Drake et al. 2009 and our discussion below).
The trajectories of individual particles and the pitch-angle dependencies of energy distributions, together with the detailed analysis of the energy gain using the GCA framework, have shown that the main contribution to the acceleration of particles represented by the guiding center is the Fermi Ib acceleration, which occurs in curved magnetic fields. This acceleration consists of a polarization effect and an inertial effect, in which the polarization effect is the more dominant component compared to the less strong inertial effect. Nevertheless, the additional influence of nonadiabatic acceleration caused by field inhomogeneities within ion gyro-orbits can be observed with increasing mass. Such behavior, which is out of the scope of adiabatic GCA analysis, can also be the main reason for the observed asymmetries in the distributions of ion energies in pitch angles. The acceleration efficiency increases with the ion massto-charge ratio, in line with the observed preferential acceleration of heavy ions measured during the SEP events. We have also shown that the preferential acceleration for heavy ions occurs when the E-cross-B drift is larger than the thermal velocity.
The broad extent of scales in our model, which covers a substantial part of the plasmoid cascade, allows us to study other relevant questions regarding preferential ion acceleration. In this respect, we have been interested whether the position of break points in the plasmoid-scaling law can be related to the preferential acceleration of some concrete ions with a favorable value of its specific charge (like the observed 3 He excess). Nevertheless, no such behavior was found, at least for our space of parameters. Another question related to the impact of plasmoidcascade scaling to the ion acceleration, which we attempt to answer, is the relation of the power-law index of the plasmoid distribution to the characteristics of populations of energized ions, namely, ion energy spectra and pitch-angle distributions, differential fluxes, or enhancement factors of abundances of energized ions w.r.t. the background chemical composition. We discuss our findings in this direction now, step by step: The high-energy tails of the ion energy distribution function obtained from simulations 3.0 × 10 4 -7.2 × 10 4 7.2 × 10 4 -1.7 × 10 5 1.7 × 10 5 -4.2 × 10 5 4.2 × 10 5 -1.0 × 10 6 Figure 9. Ion abundance-enhancement factor proxy (A) at the end of the simulation for energy bands I-IV. Energy bands are measured in energy per nucleon; the abundance-enhancement factor is computed as an average of all combinations of specific charge. The dashed line is a power-law fit.
exhibit single power-law profiles with slopes ranging from δ E = 0.8 to δ E = 0.9. Our results further show that the simulated differential fluxes obey a power law, in line with observed 3 He-rich SEP events (Mason et al. 2002). In this case, even the quantitative agreement is quite satisfactory: The spectral slopes (δ J = 0.8-1.4) are very close to those found in observations. Moreover, the rounding of spectra observed for low ion energies is in good match with observational data.
As the last characteristic of the energized ion distribution, we investigated the abundance-enhancement factor via two possible proxies. The first proxy for the enhancement factor increases progressively with Q/A, thus exhibiting a systematic preferential acceleration of the heavy ions. At the highest energy band, the heavy-ion portion can be fitted with a power law. The second proxy increases with Q/A as a power law, again in line with the observations. Nevertheless, the agreement reached in this case is just at the qualitative level-the slopes found in our model for both abundance-enhancement factor proxies (δ K = 0.4-0.8) do not match observed data.
This discrepancy could be caused by multiple reasons: (a) different CS evolution compared to naturally occurring fragmentation-discrepancy in the size distribution of plasmoids, (b) an insufficient number of simulated ions, (c) short simulation time, (d) initial ion temperatures not matching, (e) use of masses that do not exactly correspond to measured ones, (f) transport effect (adiabatic deceleration, diffusion, scattering) in the interplanetary space that is not included in our simulations and (g) the turbulent nature of the CS not the most significant factor in observed phenomena. Among them, option (a) is the most serious candidate for the explanation, and we will discuss it in more detail separately below. To eliminate the possible culprit from the rest, we have simulated an identical situation with a different temperature T = 3 MK. The results exhibit almost identical curves and slope values, caused by the still small thermal velocities compared to the E-cross B drift; therefore, variant (d) can be eliminated. Because the results remained qualitatively the same, we can also reject variant (b), which would produce different profiles each time the simulation is relaunched. Option (c) is related to the problem with heavy ions, which we have discussed in Section 3.2. Such ions require larger spatial and timescales that match their large gyro-orbits. Consequently, statistics obtained from the individual trajectories of those ultraheavy ions is likely being biased by the concrete geometry of the fields in the MHD box selected for the particle tracking.
Options (e) and (f) are out of the scope of this study and their influence was not taken into account during the simulations, thus no proper estimate can be made at this moment. Of course, we cannot preclude (option (g)) that the acceleration in the plasmoid cascade is the most relevant mechanism as an alternative explanation, e.g., a resonant wave acceleration, exists, too (Roth & Temerin 1997;Miller 1998;Liu et al. 2006;Kumar et al. 2017). Because studying the relation between the plasmoid turbulent cascade and the observed scaling of ion acceleration efficiency with the specific charge Q/A is the main highlight of this paper, we would like to discuss the option (a)-a possible mismatch between the MHD-simulated plasmoid-cascade scaling law and that present in real flares more extensively.
First, we would like to emphasize again that the power-law scaling of plasmoid sizes (as actually present in our MHD fields, in contrast to the previous studies) results naturally from the power-law statistics of the accelerated particles in our simulation. Hence, at least at the qualitative level, there is an agreement. Let us now focus on the quantitative relation between the plasmoid-scaling power-law index and the corresponding indices in TP statistics.
The relation between the plasmoid-cascade scaling and the acceleration efficiency increasing systematically with Q/A has been quantitatively discussed in Drake et al. (2009). Their analysis was performed for the "pickup" acceleration process (i.e., the "magnetization" of ions at the very beginning of the reconnection exhaust), which is a quite different context from the Fermi I process identified as the main acceleration effect in our study. Hence, we can expect a different quantitative relation between the scaling law for plasmoids and spectral indices for energetic particle statistics. Nevertheless, because Drake et al. (2009) is-to our knowledge-the only paper that presents such a quantitative relation (namely Equation (9) in the cited paper), let us spend a few lines here to compare.
The authors in Drake et al. (2009) estimate that in a pool of multiscale plasmoids, whose number scales with their size L as ( ) µ d -P L L P (in the original paper, the notation δ P ≡ α), the efficiency (expressed as a production rate of energetic ions) scales with their mass and charge as ( For the slope δ K = 3.26 found in Mason et al. (2004) this would mean δ P = 6.26 to reproduce the observed scaling of the enhancement factor.
In order to check how our scaling law for plasmoids, expressed as the spectral density of magnetic energy ( )  k M (Fourier power spectrum) displayed in Figure 2, is related to the slope of the enhancement factor (or its proxies), we have to first find its relation with the plasmoid-size distribution P(L) for the number of plasmoids of size L. To that end, we first estimate the magnetic energy in a single plasmoid (SP) of size L: The simplest assumption for B x (x = 0, y, z) as a function of distance from the plasmoid center z o (where it is zero-this place is an O-type magnetic "null point," where the quotation marks express that there is still a guide-field component present) toward its edge, where it reaches a typical value of the ambient field B 0 , is linear: 2 1 x 0 Let us consider a cylindrical-shaped (i.e., circular cross section) plasmoid and ignore the guide field B y for simplicity. Then, the magnetic energy contained in it reads where L y is the length of the plasmoid (flux rope) in the third invariant direction, assumed here as a constant factor independent of the plasmoid cross-section size. Because the magnetic energy at scale L holds: , we finally arrive at the power spectrum for magnetic-energy density in the form That is, the scaling law for (spectral) magnetic-energy density is a power law ( ) µ - k k M s , again, with the spectral index s = 4 − δ P . The spectral index s for the inertial domain in Figure 2 has been fitted to s = 2.14 (Bárta et al. 2011b), which corresponds to a plasmoid-size distribution slope δ P = 1.86. Let us note that this value is quite close to that found in the analytical (MHD) work on the so-called "chain plasmoid instability" by Uzdensky et al. (2010), where the scaling law assumes δ P = 2.
Hence, if we blindly apply the relation between power-law scalings of the acceleration efficiency as a function of the specific charge and the plasmoid-size distribution as anticipated in Drake et al. (2009), we arrive at the expected scaling of the abundance-enhancement factor d d = -= -3 1.14 K P expect . This result would indicate that we should actually expect easier acceleration for light ions. This, however, contradicts our results, which show preferential acceleration of heavier ions (albeit with a milder slope than that in observations).
In order to verify the preferential acceleration of heavier ions in our model, we attempted to characterize the acceleration efficiency by other quantities more pertinent to the analysis of the model, which, however, do not have direct observational counterparts: the energy gain dependence on Q/A and the gain normalized per nucleon. Both characteristics show a (albeit mild) preference for heavier ion acceleration, confirming thus the earlier results.
To sum up, the relation δ K = δ P − 3, devised for the ionpickup mechanism by the small islands (plasmoids) at the kinetic scale, does not hold for acceleration at larger scales where we have identified the Fermi I acceleration process as a dominant effect, as our analysis of individual ion trajectories has shown. It is also in line with the super-Alvénic particles reported above. Instead of that, we can present a (numeric) experimental relation between those two indices: The abundance-enhancement scaling in our simulation ranges in the interval δ K = 0.5-0.8 (depending whether proxy A or B is used), while the plasmoid-size scaling (under the assumptions above) can be estimated as δ P = 1.86. Hence, (the difference ranging in the interval 1.06-1.36). Applying it to the observational scaling of the abundance factor δ K = 3.26, one can find the required plasmoid-size scaling δ P ≈ 4.47, which is still quite steep but milder than the value of 6.26 found in Drake et al. (2009) for the pickup mechanism.
Concerning the plasmoid-scaling law, there is, of course, another important question regarding the applicability of our rescaling of the computational domain done in order to match the scaling break toward the dissipation scale. The scaling found by MHD modeling is a single-fluid MHD scaling. Its break was reached in a controlled way, so that the dissipative CSs remained resolved by the computational mesh. It was achieved by switching on the effective resistivity (representing actually the other terms in the generalized Ohm's law) when the current-density filamentation reached a critical scale. In reality, the terms of the generalized Ohm's law have different scales at which they become consecutively important (e.g., the Hall term starts to be essential in earlier stages of currentdensity filamentation compared to the off-diagonal pressure terms, etc.). It is therefore very likely that the scaling law has several break points (knees) separating scale domains with different (likely consecutively steeper) scaling laws. This may change the acceleration physics significantly.
In light of one of the questions that motivated our studywhether the break points at the plasmoid-scaling law can be related to the acceleration efficiency of a specific ion with favorable Q/A (having 3 He excess in mind)-one detail of our analysis may be of interest: We have observed a gap in Figure 9, which shows the enhancement factor (proxy), in the high-energy channel for a specific narrow range of Q/A. Should this finding be statistically significant, it represents a kind of "antieffect" to that observed for 3 He excess. Nevertheless, it could indicate that the positions of the knees at the plasmoid-scaling law may be relevant for the excessive acceleration efficiency of specific ions.
To conclude, the acceleration of ions in a fragmenting plasmoid-rich macroscopic current layer is efficient, capable of producing super-Alvénic particles. The mechanism contributing most to the energy gain is the first-order type B Fermi acceleration. Qualitatively, we found that there is a systematic preference for the acceleration of heavier ions; however, the scaling of acceleration efficiency with the specific charge does not quantitatively match the observed data. The exact form of the plasmoid-scaling law, namely the positions of its breaking points, may even be related to excessive acceleration efficiency for specific ions.