Turbulent Magnetic Reconnection as an Acceleration Mechanism in Earth’s Magnetotail

Using electric and magnetic fields measured by the Magnetospheric Multiscale (MMS) mission, we construct a test-particle simulation of a turbulent magnetic reconnection region to investigate observed ion acceleration. We identify three types of energized ions: (1) ion jets, (2) Speiser-like energized ions—both of which carry significant energy but do not produce a strong energetic (>80 keV) tail in the ion distribution—and (3) a separate but sizable population of ions that are accelerated to significantly higher energies (>80 keV) by the turbulent fields. The majority of ions that undergo energization by the turbulent fields cross the magnetic null plane multiple times. By preferentially energizing these particles, the turbulence creates a separate population of ions that mostly exits in the dawn direction of the magnetotail and forms a high-energy power-law tail in the ion flux-energy distribution. We also find that the highest acceleration energies are limited by the size of the turbulent region (with respect to ion gyroradii). Since turbulence is widespread in astrophysical plasmas and has no a priori limit on scale size, the MMS observations suggest turbulence may have a significant role in particle acceleration.


Introduction
In a magnetic reconnection region of Earth's magnetotail, the Magnetospheric Multiscale (MMS) mission measured non-Gaussian distributions with accelerated ions (Ergun et al. 2020a); how these accelerated ions are energized is not well understood.Superthermal ions of energies up to 100 keV have been seen in previous magnetospheric missions (DeCoster & Frank 1979;Christon et al. 1991;Keiling et al. 2004;Grigorenko et al. 2009;Artemyev et al. 2010;Haaland et al. 2010), and theoretical literature followed proposing acceleration mechanisms.Energizing via the dawn-dusk electric field of magnetic reconnection has been considered (Speiser 1965;Zelenyi et al. 2004Zelenyi et al. , 2007) ) as well as jet acceleration (Litvinenko & Somov 1993;Zharkova & Gordovskyy 2004) and interaction with electromagnetic fluctuations (Artemyev et al. 2009;Perri et al. 2011).An acceleration mechanism that includes stochastic heating and generates magnetic holes has been shown to explain the energy distribution of electrons (Dolgonosov et al. 2013;Ergun et al. 2020b), but the case of energetic ions is more complex, due to their slower initial speeds and their partially unmagnetized orbits (e.g., serpentine orbit ;Speiser 1965;Somov 2013).When turbulence is present, observed ion distributions have a substantial nonthermal tail extending to 100 times the initial ion thermal energy.
The development of a power-law distribution generally requires an acceleration mechanism that favors energetic particles (Blandford & Eichler 1987).Fermi-like processes, of which there are several flavors, are often invoked; for example, simulations show that repeated impulses to a set of partially trapped electrons in magnetic islands can cause a power-law distribution (Drake et al. 2006).Likewise, it has been advanced that turbulence in the near-Earth magnetotail also favors energetic electrons (Ergun et al. 2020a(Ergun et al. , 2022;;Usanova & Ergun 2022).While observations indicate that turbulence plays a role in ion acceleration, it is unclear what feature of turbulence favors energetic ions.
The MMS observations of strong turbulence-with δ|B|/|B 0 | ∼ 1, where B 0 is the background magnetic field and δ|B| its fluctuations-suggest that at least three energization mechanisms are active (Ergun et al. 2018), all necessary to understand the energy balance that occurs in the reconnection region.These are stochastic turbulent acceleration, ion jets, and the dawn-dusk electric field, E y .The combined effect of the first two mechanisms-stochastic acceleration by turbulent electric fields and bulk acceleration of ions into opposing jetscarries a significant fraction of the magnetic field annihilation energy away from the x-line.Likewise, the measured reconnection electric field (E y in this article) is sufficient to energize ions to 10 times the initial thermal energy, supporting that energization, due to advancing along the E y dawn-dusk field can also be active (we refer to this energization as Speiserlike).An important question that needs to be resolved is the relative contribution of these three mechanisms and how they interact.
Test-particle simulations are often employed to elucidate the acceleration mechanism for ions.They subject ions in a magnetic field reversal region to various electromagnetic field models and allow for the analysis of their individual trajectories as well as their statistical behavior.Veltri et al. (1998), Greco et al. (2002), and Taktakishvili et al. (2003) studied the motions of ions under static magnetic fluctuations, and Greco et al. (2009) and Perri et al. (2011) included electromagnetic fluctuations as clouds of electromagnetic fields that oscillate throughout the reconnection plane; this scheme was proposed to suggest a Fermi acceleration mechanism, which is known to produce power-law tails in the most energetic parts of a distribution (Fermi 1949;Davis 1956).Dolgonosov et al. (2013) showed that a power-law index of −4.45 could indeed be produced with this mechanism in the relevant energy range (>80 keV), whereas Artemyev et al. (2009) found a negative result (no power-law tail) by adding the electromagnetic fluctuations as phase-mixed plane waves following a powerlaw spectrum.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
The test-particle simulation used in this article is based on measured turbulent magnetic (δB) and electric (δE) fields and has several advantages over many current self-consistent codes.To properly reproduce the measured δB and δE, a selfconsistent code must have an enormous range of scales, a representative large number of particles at the smallest scales, and a correct electron-to-ion mass ratio.To properly treat ion acceleration, a code must have a very large 3D simulation domain, long run times to allow ions to pass through the system, and very importantly, fully open boundaries.A testparticle code can incorporate all of these features at the expense of self-consistent behavior.Moreover, most of the parameters of our simulation are well constrained by measurements: the particle distributions (densities and temperatures) at the boundaries, as well as the size of the turbulent region, can be estimated from data (Angelopoulos et al. 1994;Ergun et al. 2018Ergun et al. , 2020a)), the only exception being the length of the x-line of magnetic reconnection.
Basing the test-particle simulation on the event described in Ergun et al. (2020a), we find that, while energization by ion jets from magnetic reconnection and Speiser-like processes have substantial contributions, the most energetic ions (>80 keV) result from turbulent fields.Dynamically, most of these energetic ions cross the magnetic null plane multiple times, which results in higher energization by the turbulent fields, and as we show, present the highest energization rate at z = 0.By preferentially energizing these serpentine ions, the turbulence creates a separate population of ions that mostly exits in the dawn direction of the magnetotail and forms a high-energy, power-law tail in the distribution.
In Section 2, we provide a brief overview of the MMS observations.In Section 3, we describe the test-particle simulation domain.In Section 4, we show how a simulation with a 2R E long x-line reproduces the MMS-measured ion flux and density depletion seen in the turbulent region.We explain the mechanics of the acceleration by following the trajectories of superthermal ions and changing the amplitude of the turbulent fields in order to expose the interplay between turbulent, Speiser-like, and ion jet energization.In Section 5, we discuss the implications of our findings for larger-scale regions of magnetic reconnection in astrophysics, and in Section 6, we end with our conclusions.

MMS Observations
A subset of MMS observations of accelerated ions in a turbulent magnetic reconnection event is displayed in Figure 1.The event is in Earthʼs magnetotail, roughly 23R E from Earth.Panel (a) displays the intensity of energetic ions from 60-600 keV measured by the Flyʼs Eye Energetic Particle Spectrometer (FEEPS; Blake et al. 2016).Given that the nominal ion omnidirectional fluxes in the plasma sheet have an average energy of 4-5 keV, these perturbed fluxes have been energized by more than a factor of 10.The concurrent magnetic fields, electric fields, and plasma density are shown in the panels below.The magnetic field signal shows that the fluctuations (δB) have nearly the same amplitude as the background field (B 0 ), suggesting strong turbulence.Intense electric fields and a strong density depletion are evident.
Figure 1(e) shows the ion intensity averaged over a 2 minute period of the turbulent region.The blue squares represent data from the Fast Plasma Investigation (FPI; Pollock et al. 2016), and the red circles represent data from FEEPS.The orange line represents a power-law fit with an index of −3.9.The appearances of the turbulent fields' amplitude and the energetic ions are strongly cotemporal; hence, the ions appear to be locally accelerated.From arguments elaborated in Ergun et al. (2020a), the high-energy ion fluxes <200 keV are predominantly protons.Measured ion fluxes with energies above 200 keV are not used in our analysis.

Test-particle Simulation
We design a test-particle simulation after the above event to better understand ion acceleration.The basic idea is to reproduce realistic background and turbulent electromagnetic fields in order to analyze the kinematic behavior of the ions under those fields by following their trajectories.Moreover, the features of the ion fluxes from the simulation can be compared to those of the measured ion fluxes.In this section, we describe the simulation domain (Section 3.1), the background fields (Section 3.2), the turbulent fields (Section 3.3), and the boundary conditions (Section 3.4).More details about the simulation can be found in the Appendix.

Simulation Domain
Figure 2 shows a sketch of the simulation domain, which contains a smaller turbulent region.Ergun et al. (2020a) estimated the physical extent of the turbulent region to be roughly (within a factor of 2) 16R E in the direction of the reconnecting magnetic fields (the x-direction), using an estimated retreat velocity of ∼100 km s −1 and the fact that the turbulence and accelerated ions are detected for roughly 16 minutes.The turbulence appears to extend in the z-direction (normal to the current sheet) to the plasma sheet boundary layers (south to north), which gives a z-extent of roughly ∼1.5R E (Baumjohann et al. 1989).This distance is in consort with the expected extent of the ion diffusion region of magnetic reconnection normal to the current sheet of several ion skin depths.
The extent along the y-direction (x-line of magnetic reconnection) cannot be established, given the close spacing of the MMS spacecraft.Observations from multispacecraft missions, such as Geotail and Cluster (Nakamura et al. 2004;Grigorenko et al. 2009;Artemyev et al. 2010;Haaland et al. 2010), indicate the x-line can be anywhere from 1-3R E .In the simulation, we treat the length of the turbulence along the xline as a free parameter, 2R E being the value that best reproduces the ion fluxes observed by MMS (Figure 1); the degree of ion energization is sensitive to this dimension.
The dimensions of the simulation domain are 16R E × 4R E × 4R E .The turbulent fields δB and δE are only present in a smaller region of dimensions 14R E × 2R E × 1.5R E .A noturbulence buffer zone of at least 1R E is introduced to allow high-energy particles to reenter the turbulent region during a gyration.

Background Fields
The event shown in Figure 1 occurs in a highly turbulent region of magnetic reconnection.Even though the turbulence visibly obscures the background fields of magnetic reconnection (B 0 and E 0 ), they can be estimated via data smoothing combined with modeling using well-established properties of magnetic reconnection.The background magnetic field in the simulation domain consists of an unperturbed, sign-reversing component in the x-direction modeled after a Harris-like current sheet.The same reversal configuration is applied in the z-direction (the north-south direction) to form a magnetic reconnection region with an 8:1 ratio.The background fields are described as where λ x = 8λ z and λ z is the half-thickness of the current sheet, for which we use an ion skin depth of w = d c i p i (where w p i is the ion plasma frequency and c is the speed of light).From the average measured plasma density, the ion skin depth is determined to be ∼1000 km (Ergun et al. 2020a).The extent of the simulation box is ∼100d i in x by ∼25d i in z.B 0x is set to the measured value of 20 nT (Ergun et al. 2020a), and B 0z is 1/8 of that value, 2.5 nT.The blue lines in Figure 2 represent the background magnetic field lines.Note that the signreversing north-south magnetic field makes the total magnetic field at z = 0 nonzero (except at points along the x-line), and hence, it plays the role of B n in Greco et al. (2002) and Perri et al. (2011).
The background electric field is with E y set to the measured value of 2.7 mV m −1 (Ergun et al. 2020a) within the turbulent region and reducing to E y = 1.35 mV m −1 in the buffer zone.

Turbulent Fields
A chief advantage of a test-particle simulation is that the turbulent electromagnetic fields, δB and δE, can be designed to closely mimic those in observations.In this effort, we employ a set of pseudo-randomized electromagnetic and electrostatic waves designed to match the measured properties of δB and δE, including amplitude, wave speeds, coherence times, coherence scales, and the measured power spectra (δB shows a Kolmogorov (1991)-like turbulent power law; see Appendix).We reconstruct δB with 200 plane waves: where A n is the amplitude of wave n, k n is its wavenumber, ω n is the angular frequency, f n is an arbitrary phase, and  r is the 3D position vector.The amplitudes and directions of A n are pseudo-random, biased to match the measured spectra, ω n and k n are set to mimic the measured speeds (see Appendix), and f n are randomly assigned (0-2π) and vary in time and space to match the measured coherency times and distances as a function of frequency.
Once δB is determined, the electromagnetic part of δE is constructed directly from Faradayʼs law.An electrostatic component of δE is also developed emphasizing the perpendicular (to B 0 ) component using the same procedure.The wavenumbers are estimated from measured speeds of δE.A small parallel contribution is added to match the δE || spectrum.A comparison of the reproduced δB and δE spectra with the measured spectra is in the Appendix.
The total fields are then given by The equation of motion for protons is mdv/dt = e(E + v × B).

Boundary Conditions
The particles that enter the plasma sheet from the lobes have a source located somewhere in the magnetosphere mantle (Russell et al. 2016), and hence historically, the temperature of the injected particles in simulations has been 0.1-1 keV (Greco et al. 2002).However, MMS measures a distribution outside of the turbulent region with a temperature of 4 keV, which we employ.The measured distribution can be described as a shifted Maxwellian with a temperature of 4 keV.
The plasma density at the simulation boundaries (outside of the domain) is based on a simple, quiescent Harris-like current sheet model, which varies only with distance from the background magnetic field reversal at z = 0. We impose a constant density and temperature as a function of x and y at the simulation boundary so that the resulting densities/temperatures inside the simulation domain can be attributed to the imposed turbulent fields.The plasma density in a Harris-like current sheet decreases to a value of n lobe at high-z: where z n0 is the thickness of the current sheet, which is set to be a skin depth of 1 ion.Observed densities when B ∼ 0 prior to and after the event range from ∼0.1 to ∼0. , so the ratio n sheet /n lobe is fixed at 4. In the simulation, ions are injected at the boundaries to mimic the prescribed density and temperature.Ions that exit the domain are removed.

Results and Analysis
In this work, most of the numerical results are obtained by averaging the results of multiple runs.A run starts by initializing 1800 particles in the simulation domain at t = 0, with the density profile described in Equation (6) and a temperature of 4 keV.During a run, one particle is randomly injected at every time step of dt = 0.01 s, following the density profile and fixed temperature (4 keV) at the boundary.Each run advances in time until it reaches a steady state at t = 100 s (no results used), after which few of the initialized ions remain.The resulting ion fluxes, temperature, and density are recorded from t = 100-300 s, during which ∼20,000 ions pass through the simulation domain.Each run has a unique set of 600 pseudorandomly generated waves (δB and δE) as described above.We distinguish a run from a simulation.A simulation is an ensemble of runs where the random variables are regenerated, including δB and δE.To compile distributions, for example, we use an ensemble average of 50 runs during which roughly 1 million particles pass through the system and 30,000 waves are generated.
The main goals of Section 4.1 are to (1) determine the extent of the turbulence in the y (dawn-dusk) direction, which is not well constrained in the observations, and (2) demonstrate that the test-particle simulation reproduces MMS observations reasonably well.In Section 4.2, we examine the acceleration mechanisms that allow ions to reach the high energies (>80 keV) seen in the data.In the acceleration study, we show that the motion near z = 0 is most critical and that the majority of the energetic particles exit through the +y (dusk) face.By analyzing individual and collective trajectories, we determine that energization favors ions crossing the magnetic reversal plane.Moreover, once energized after crossing the magnetic reversal, those particles are more likely to further gain energy, which leads to the development of a power-law tail (acceleration).To demonstrate the impact of turbulence, we compare the results of simulations with and without turbulence imposed.
We start by showing how our simulation can reproduce the MMS-measured omnidirectional fluxes, temperatures, and densities in the turbulent region.

Baseline Simulation
The baseline simulation employs parameters determined by MMS, including B, E, ion boundary temperature and density, and the extent of the turbulent region.The dawn-dusk extent of the x-line (y-direction) is a poorly constrained but critical parameter.To determine this distance, the simulation is performed with x-line lengths of 1, 2, and 4 R E .The extent in the y-direction is then chosen so that ion fluxes best agree with those observed by MMS.
Figure 3 compares the simulated fluxes with x-line lengths (extent of turbulent region in y) of 2 and 4R E to the observed fluxes in Figure 1(e).The simulated ion fluxes are plotted in physical units by properly weighting ions as they enter the domain.The plasma density (derived from electrons, Figure 1(d)) during the period of the observed distribution averages ∼0.05 cm −3 .The test-particle distribution has the same density if n sheet = 0.175 cm −3 .Comparisons are made with x-line lengths of 1, 2, and 4R E .The measured flux distribution (blue and red squares) is best matched by the simulation with an x-line length of 2R E (black).A 4R E x-line length (green) in the simulation increases the energetic (>80 keV) ion fluxes by more than a factor of 3, whereas a 1R E x-line length (not displayed) significantly decreases the energetic fluxes.As a note, the simulated fluxes are in the rest frame of the x-line, whereas the measured fluxes are in a moving frame of ∼100 km s −1 (Ergun et al. 2020a).This difference in frame does not significantly impact the simulation to observed flux comparisons, which are at energies >3 keV.
Simulated fluxes were compared with observed fluxes from two other time periods, one early in the event and one late in the event (not shown).The observed plasma densities vary by more than a factor of 2. These comparisons also support an xline length of roughly 2R E .An x-line length of 1.5R E results in the best match between observed and simulated fluxes if the simulated flux distributions are calculated from the +y side of the turbulent region (rather than over the entire turbulent region).An x-line length of 2.5R E could be justified if the simulated flux distributions are compiled from the −y side of the turbulent region.Given strongly varying observed densities and the dependence on where the simulated fluxes are compiled, we can only estimate the x-line length (turbulent region length in y) to be between 1.5 and 2.5R E .Given the quantitative comparison displayed in Figure 3, we believe that the simulation with an x-line length of 2R E reproduces the energetic acceleration well and can lead to meaningful conclusions.
Figure 4 (top) shows the ion temperature as a function of run time.The initial transient from t = 0-50 s shows rapid ion energization.One can see that the test-particle simulation achieves a steady state before 100 s.The resulting steady-state ion temperature of ∼14 keV is in good agreement with that reported in Ergun et al. (2020a).Figure 4 (bottom) shows the ion density in the center of the simulation domain as a function of x.The density inside the turbulent region is significantly depleted from its value outside of the box, which is consistent with the data shown in Figure 1(d).
As the system enters its steady state, pressure balance is maintained (within a factor of 2) because the temperature increases as the density depletion onsets.Deviations from the pressure equilibrium are to be expected, given the magnitude of the turbulent fields.It is also to be expected that the pressure balance is not exact but remains within a factor of 2 throughout the simulation, given the absence of electrons.This rough pressure balance implies that the magnitudes of the currents and fields involved in the simulation do not evolve to unrealistic values.Greco et al. (2002) refer to this as a posteriori self-consistency.Nevertheless, while the simulated ion currents remain consistent with the background fields, the lack of self-consistency still implies that ions cannot damp or grow the waves, which calls for the careful construction of the imposed turbulent fields (see Appendix).
In order to better understand how the energetic tail is generated, we ask where ions are energized the most.To do so, we record the instantaneous energization rate at 0.1 s intervals as a function of position and find that the majority of the energization (and de-energization) occurs near z = 0, where the background magnetic field is the weakest (see Figure 5).Figure 5 also shows nearly equal rates of energization and deenergization in the turbulent region.Similar results regarding the z dependence of the heating rate were reported in Greco et al. (2002) and Artemyev et al. (2009), and were mostly interpreted as acceleration due to the dawn-dusk electric field E y of ion in serpentine-type orbits (a process characterized in Speiser 1965).However, we also find a net positive energization that emerges when we sum the negative and positive changes in energy produced by the turbulent fields, which is characteristic of second-order stochastic acceleration (Miller et al. 1996;Somov 2013).Speiser-like and stochastic processes combine and interact to generate a net energization when turbulence is present; a closer analysis is needed to separate and understand the two effects.
The agreement between simulation and observation in the ion fluxes (Figure 3), ion temperature, and density (Figure 4), and consistency with previous analysis (energization near z = 0) (Greco et al. 2002;Artemyev et al. 2009) makes the case that the simulation is properly representing the observed ion acceleration, and can be used to better understand the mechanism behind this energization.To unveil the role of turbulence versus Speiser-like energization, we take a closer look at individual ion trajectories of ions with and without turbulent fields.

Ion Trajectories and the Role of Turbulence
In Figure 6, we classify the most common trajectories we observe in the no-turbulence simulation by (1) their entering and exiting faces of the ions and (2) whether they cross z = 0 more than once.This classification yields six types of orbits.
(1) Jet-forming ions: These enter through a z face, drift toward the magnetic null, and exit through an x face, often accelerated to near-Alfvén speed (Figure 6(a), blue trace).From a 2D magnetic reconnection point of view, most particles forming the ion jet partake in this trajectory; however, we find that this ion jet trajectory forms roughly 30% of the ion jet in the 3D case.(2) Fermi-reflected ions: These enter and exit through the ±x faces (Figure 6(a), violet trace), and make a The blue squares and red squares represent MMS observations extracted from the event in Figure 1; the red squares fit a power-law index of −3.9.The dashed green lines show the ion flux for a simulated x-line twice as long (4R E ).With respect to the red line, the high-energy flux increases by more than a factor of 3 when the x-line is doubled in size.similar (30%) contribution to the ion jet.(3) x-entering Speiser orbit: Ions that enter through the ±x faces and migrate along y with a serpentine orbit, gain additional energy, and leave the +y face (Figure 6(b), green trace).(4) Speiser-enhanced jetforming ion: These enter the −y face, gain energy while partaking in a serpentine orbit, and then exit the ±x face (Figure 6(b), red trace).These orbits are also common (likely due in part to the limited y extent of the simulation domain) and dominate (>60%) the highest-energy part of the ion jet.(5) zentering Speiser orbit: These enter ±z, and exit +y (Figure 6(c), purple trace).They are the most common +y exiting trajectory.(6) Speiser orbit: These enter from −y and exit through the +y face, while forming a characteristic serpentine trace (Figure 6, orange trace).They are relatively rare in our simulation, comprising only about 5% of the net energy flux exiting the +y face in the no-turbulence case.
To study how turbulence changes this picture, we do two 1000 s runs: one with turbulence and one without (δB and δE set to zero).For each run, we count the ions leaving the +y face, while recording their energy and whether they cross the z = 0 plane.The results are presented in Figure 7.In order to highlight the differences between the plotted energy fluxes, we normalize the result such that the reference 4 keV Maxwellian entering energy flux peaks at unity.A primary difference between the no-turbulence and turbulence cases is the dramatic increase in >80 keV fluxes when turbulence is present (Figure 7, right panel).In both cases, the most energetic ions cross the z = 0 plane before leaving the +y face (green traces in Figure 7), but the turbulence further separates this energetic population into a secondary population of highly energized ions, effecting a bimodality in the flux-energy distribution.We identify a quasithermal peak at an energy of 9 keV (red trace, right panel) and a highly energetic peak at 85 keV (green trace, right panel).
The peaks in quasi-thermal cores of the exiting fluxes (red traces in Figure 7) have nearly the same energy.The turbulence case has a slightly higher flux level.The increase in energy of the quasi-thermal core of the distribution (from 8-9 keV) is consistent, roughly, with moving along Ey a distance of an ion gyroradius, which is the most likely displacement in the noturbulence case (Taktakishvili et al. 2003;Russell et al. 2016).However, a displacement along E y cannot explain the energies gained at the energetic end in the turbulent case, which is considerably greater than what the E y field can yield, even when considering displacements across the entire simulation domain.
We can further examine the relationship between Speiserlike energization and energization by turbulent fields by plotting energization as a function of the displacement Δy.In Figure 8(a), the dots represent the ions that leave the +y face in the same 1000 s run that we use for Figure 7.The abscissa and ordinate values correspond to the ion's final energy gain and displacement, respectively, and the blue and black dots correspond to the simulations with and without turbulence.The relationship shown suggests that the acceleration process that generates the tail of the energy distribution still has a linear trend with displacement, the final energization of the ions being a combination of the effects of the Speiser and turbulent fields.Although it is expected in the no-turbulence case to see a monotonically increasing trend between the displacement in y and the energization of the particle, the slope of this dependence steepens when turbulence is present.
Figure 8(b) displays the net exiting particle energy flux summed over all boundaries as a function of the turbulence amplitude.Zero on the horizontal axis represents no turbulence, and unity represents turbulent amplitudes (δE and δB) as measured by MMS.The exiting energy fluxes are normalized to the total particle energy flux entering the simulation domain.With no turbulence, Speiser and ion jet acceleration combine to increase the net exiting particle energy flux (exiting minus entering) by 0.73 of the entering particle flux (Figure 8(b), black circle at 0).As the turbulence amplitude increases, the net exiting energy flux increases quadratically, as expected from second-order stochastic energization (see Somov 2013, Chapter 14).At the full turbulence level, the net particle energy flux is doubled.Most noticeably, the total exiting energy flux of ions  As discussed earlier, the degree of ion energization is sensitive to the dawn-dusk dimension.Nevertheless, we find that the ratio of turbulent to Speiser-like energization does not vary as strongly with x-line length.The energization of both types decreases with a shorter x-line, and both increase with a longer x-line.On the other hand, the ratio of turbulent to Speiser-like energization strongly depends on the ratio of the turbulence amplitude to the amplitude of E y , which can dramatically vary from event to event.
As a consequence of the interplay between Speiser-like and turbulent energization, the net energization of +y exiting ions, in particular, is greatly enhanced by the turbulent fields.This stems from the fact that ions that cross z = 0 once tend to cross it multiple times, and the energization rate is also significantly larger for particles near z = 0 (see Figure 5).Therefore, ions that get energized due to their proximity to z = 0, tend to have  trajectories that favor more subsequent energization, which leads to a runaway that generates the highly energized population seen in Figure 7.
That electromagnetic and electrostatic turbulence does not affect the average dwell time of particles complements earlier results (Veltri et al. 1998;Greco et al. 2002;Taktakishvili et al. 2003), which show that magnetostatic turbulence, contrastingly, increases the dwell time and the average Δy of particles, effecting a net acceleration.The combination of electromagnetic and electrostatic turbulence, on the other hand, decreases and increases the dwell time of individual ions by roughly the same amount on average.Examining the trajectories in Figure 9, we can see both effects on the dwell time at play.We show four trajectories of ions in runs without (top) and with (bottom) turbulence, the two trajectories shown for each case having the same initial conditions.The blueyellow trajectory is a known Speiser orbit with no turbulence: it enters the x face, drifts into the current sheet (E × B drift), and then exits the +y face.On the other hand, the green-red orbit contributes to the ion jet with no turbulence: it enters the −y face, then follows the magnetic field out of the +x face.The dwell times for these orbits are represented in the color bars at the lower-right corner of the plots.While the dwell time of the blue-yellow trajectory increases, due to turbulence (and its energy), the green-red trajectory accelerates through the region faster also due to turbulence.We find that when applying MMS-based electromagnetic and electrostatic turbulence, the turbulent acceleration exceeds the Speiser energization mechanism, and together they can energize particles by over a factor of ∼20, achieving energies greater than 80 keV.
Figure 10 shows the fluxes out of each of the faces of the simulation domain, without turbulence (top row), and with turbulence (bottom row).The red region highlights the >80 keV energies of the flux-energy distribution.Fluxes exiting the ±x faces in the no-turbulence case show the expected ion jets common of magnetotail reconnection (e.g., Angelopoulos et al. 1994).Speiser-enhanced ions, such as the one depicted in the green-red trajectory in Figure 9 (top),  .Top: two simulated ion trajectories without turbulence.Blue-yellow trace: Speiser-like trajectory that drifts toward the magnetic null plane and then is directed toward the +y face.Green-red trace: a particle that contributes to the ion jet that follows the hyperbolic magnetic field and exits the +x face.Bottom: the same ions but with turbulence.Blue-yellow trace: a nonthermally accelerated ion.When the particle crosses the magnetic null plane, it engages in a Speiser-like orbit, and the dwell time and energization increase.The energized particle is redirected toward the +y face, due to the background electric field.Green-red trace: the ion is redirected toward the +y face, due to turbulence, and the dwell time decreases.contribute the majority of the ions in the most energetic part of the jet (>30 keV).Without turbulence, we find the expected result that the majority of the energetic ions leave through the combined ±x faces; however, a substantial number of energetic ions exit the +y face as well, likely due to the relatively short extent of the turbulence in the y-direction.In the no turbulence case, very few ions are energized to over 80 keV.
Once the MMS-based electromagnetic and electrostatic turbulence is applied, the net energy flux out of the ±x faces increases, but the majority of the energetic ions exit through the +y face (Figure 10, bottom row).The net energy flux out of the +y face dramatically increases, and the energetic tail (>80 keV) forms.
Therefore, as shown by the trajectories in Figure 9, to understand the ion energization process, one needs a 3D simulation with open boundaries, given that the ion jet acceleration and the Speiser acceleration and turbulence interact and enhance the total current leaving tailward.The inclusion of this feature in our simulation, together with the MMS-measured turbulent field, provides the necessary component to produce the measured power law.Comparisons between our scheme and others will be discussed in the next section.

Discussion
In modeling an observed turbulent magnetic reconnection event, we have found that accelerated ion jets, Speiser-like energization, and turbulence are all necessary ingredients to explain the net acceleration of ions in the magnetic reconnection current sheet of the magnetotail.Veltri et al. (1998), Greco et al. (2002), andTaktakishvili et al. (2003) did not include the electric field turbulence because they theorized that the measured magnetic turbulence by Geotail was due to tearing instabilities.However, MMS measured a turbulent electrostatic field that is significantly stronger than expected when considering the inductive field of the magnetic turbulence alone.
The impact of the stronger electric fields is clear when considering that in Artemyev et al. (2009), where the only electric component comes from this inductive field, and no ions have energies of 100 keV or more.The main difference between the work of Artemyev et al. (2009) and this work is the introduction of an electrostatic component of δE to the power spectrum of the turbulent field.Secondarily, Artemyev et al. (2009) use a cluster-measured spectrum for the plane waves with a single index of - 7 4 , and their turbulence is based on a measured correlation length, while we use MMSmeasured indices (see Appendix) and correlation length and time in our reconstruction of the turbulent field.We believe that the principal difference in the ion energies in these two simulations stems from the inclusion of the electrostatic component of the turbulence.
The effects of the electrostatic turbulence in the energies also become apparent when comparing to the magnetostatic-only turbulence case studied by Veltri et al. (1998), Greco et al. (2002), and Taktakishvili et al. (2003).The highest ion energies in our simulation are an order of magnitude higher compared to the ones obtained in their simulation.Nevertheless, we confirm some of the findings in Veltri et al. (1998) and Greco et al. (2002): the density depletion centered at x = 0 and the thickening of the Speiser layer.Unlike Veltri et al. (1998) and Greco et al. (2002), however, we do not see a reverse current layer forming on top and below the quasi-current sheet.Such a structure appears to be hard to maintain in an environment with high electric field fluctuations (δE).Note that the density depletion has been measured by MMS and other spacecraft (Artemyev et al. 2010;Ergun et al. 2020a), and while the pressure remains mostly balanced by the corresponding increase in the temperature within the turbulent region, the existence of extra turbulent magnetic field pressure is still important to balance the system (this is consistent with Greco et al. 2002).
The only relatively unconstrained free parameter in our model is the length of the y-axis.All other parameters are constrained by either MMS measurement or theoretical consistency.This prompts the question of how the energization changes if we change the size of the turbulent region.The dashed green line in Figure 3 indicates that energization increases significantly with the y dimension.Figure 8(a) suggests a linear relation can represent the dependence between the displacement and the energization in the turbulent case, which invites speculation on the energization that can occur in larger regions in other astrophysical contexts.
There is no theoretical limit on the scale of turbulence.On the other hand, the extent of an x-line of magnetic reconnection or the dimension of a current sheet could be limited or discontinuous, given the increased pressure from ion acceleration.Nonetheless, widely scattered turbulence could extend over scales far greater than that in Earthʼs magnetotail.In supernova remnants, the high-energy emission region has been estimated to be of the order of 10 −2 pc, and the Neumann layer (where turbulence is expected to be strong) is estimated to extend 10 −4 pc (Zhang et al. 2018).Even if strong turbulence has a limited filling factor in these extensive regions, ion acceleration could plausibly contribute to the solution of what Fermi called the ion acceleration problem: ions need to already be accelerated to 200 MeV to partake in the Fermi acceleration mechanism in supernova shock fronts (Fermi 1949;Davis 1956).
From MMS observations, it appears that the strong turbulence is enabled by magnetic field annihilation due to magnetic reconnection.The substantial ion energy flux that exits the +y face of the domain also opens the question of how strong turbulence and the resulting ion energization can influence magnetic reconnection.It has been observed that the x-line length in the magnetotail (which extends over 20R E ) rarely exceeds ∼3R E (Nakamura et al. 2004).The test-particle simulation results may provide an explanation.The ion energy density (pressure) is significantly higher on the +y side of the domain (e.g., Figure 8(a)).The increased ion pressure may limit the extent of an x-line in turbulent magnetic reconnection.

Conclusions
In this work, we employed a test-particle simulation that recreates an observed turbulent magnetic reconnection region of Earthʼs magnetotail.The plasma parameters, including density profiles, ion temperatures, and magnetic fields, are based on measurements.The turbulent electromagnetic and electrostatic fields are reproduced by taking advantage of the four-spacecraft MMS mission, which allows for the estimation of correlation lengths, correlation time, and wave speeds, as well as detailed spectral properties and accurate parallel and perpendicular amplitudes.The dimensions of the turbulent region are estimated from observation and/or constrained by theory, the only exception being the extent along the x-line, which is treated as a free parameter.The test-particle simulation employs fully open boundaries allowing particles to enter and exit the domain as they would in nature.
The test-particle simulation is able to reproduce the measured ion distributions, density depletion, and ion temperatures, which corroborates the realistic reproduction of the turbulent magnetic reconnection region.Even though testparticle simulations are not self-consistent, the boundary conditions are designed to achieve approximate pressure balance.The ability to reproduce measured ion properties lends support that one can derive meaningful conclusions on the energization process.
The test-particle simulation in this article focuses on the case of strong turbulence in a magnetic reconnection region.We found that (1) electromagnetic and electrostatic turbulence is largely responsible for the generation of the power-law tail in the ion distribution.We also found that (2) all acceleration mechanisms (the ion jets, the Speiser drift, and the turbulence energization) contribute significantly to the overall dynamics of the ions; the largest energy input comes from the turbulence itself.Comparing several simulations that have been done before supports several of their findings, but highlights the necessity of imposing realistic electromagnetic and electrostatic fields and employing open boundaries.Moreover, we found that (3) the presence of turbulence significantly enhances the number flux and energy flux of ions out of the +y face.The most energetic part of this flux has crossed z = 0 multiple times, and hence, is dwelling in the region where we find the highest energization rates in our simulation.
We speculate the consequence of scaling this picture to a larger reconnection region, specifically, those generated by the shocks of supernovae shells.When considering such scales, we find that the reconnection mechanism described here can accelerate ions to the order of hundreds of MeV, which can contribute to generating an ion population energetic enough to be further accelerated by Fermi acceleration and generate the energetic part of the cosmic-ray spectrum.

Appendix Description of 3D Test-particle Simulation and the Field Reconstruction
We designed the test-particle simulation to recreate an event in the magnetotail of Earth by imposing observed conditions and tacking ions as they pass through the domain.Importantly, ions can enter and exit the do+main as they would in the magnetotail (fully open boundaries).In treating acceleration, escape is critical (Blandford & Eichler 1987).The 3D testparticle code used in this study was based on a 3D quasi-selfconsistent simulation of the Parker Solar Probe interaction with the solar wind done by Ergun et al. (2010), whose major results were ultimately verified with fully self-consistent simulations (Marchand et al. 2014).The quasi-self-consistent code also was used to characterize a Maven instrument (Ergun et al. 2021) The simulation domain is shown in Figure 2. The conditions at the boundaries, T i = 4 keV, n i = 0.1 cm −3 , and |B Lobe | = 20 nT are based on observations (d i ∼ 1000 km; ρ i ∼ 320 km).The domain includes a magnetic reconnection plane that extends 16R E ( ∼ 120d i ; ∼ 400ρ i ) in the x-direction (the direction of the reconnecting magnetic field) and 4R E ( ∼ 25d i ; ∼ 80ρ i ) in the z-direction (current sheet normal).The y-direction along the nominal magnetic reconnection x-line is adjusted to 2R E ; the simulation domain extends 4R E .
The simulation contains a background B x that varies as ( ) l z tanh z and B z that varies as ( ) l x tanh x (Equation ( 4)).The asymptotic value of |B x | is 20 nT, and of |B z | is 2.5 nT, matching the observations.B y = 0, and there is no change in the background B as a function of y.This background B mimics that of magnetic reconnection (shown in Figure 2).The background electric field, E y , is set to 2.7 mV m −1 in the turbulent region to match observations (Ergun et al. 2022).It is reduced by a factor of two outside of the turbulent region.
The simulation has open boundaries.Ions routinely enter from all boundaries with fluxes of a 4 keV Maxwellian distribution with a density of 0.1 cm −3 at z = 0 (n sheet ).The density varies as a function of z (shown in Figure 2), lowering to 0.025 cm −3 (n lobe ) at z = ± 2R E as described in Equation (6).The values of n sheet and n lobe are based on observations.

A.2. Turbulent Electromagnetic Fields
The turbulent region is confined in x at ±7R E , in y at ±2R E , and in z at ±0.75R E .The boundaries of the turbulence region are not abrupt.The turbulence is ramped to full amplitude following a cos 2 shape starting with no turbulence at 0.25R E outside of the turbulent region, it reaches reaching full strength at 0.25R E inside of the turbulent region.This subdomain in the

Figure 1 .
Figure 1.MMS 2 observations of a magnetic reconnection event in 2017 July 28.Panel (a) displays high-energy ions measured by FEEPS; these include ions with energies >80 keV.Below are the turbulent magnetic fields (b), electric fields (c), and density depletion in the region (d).Panel (e) shows the ion fluxes extracted from the measurements of the FEEPS instrument (red squares), and the FPI instrument (blue squares).The flux measurements at energies <3 keV (gray squares denote contamination with high-energy electrons and protons).

Figure 2 .
Figure 2. Sketch of the simulation box.The turbulent region is shaded in tan.The boundaries of the turbulent region are gradual, extending 0.25 R E into and outside of the stated boundary.The blue lines depict the background magnetic fields.The diagram at the right shows the boundary density as a function of z.

Figure 3 .
Figure 3. MMS data compared to the baseline simulation.The black dots represent the omnidirectional ion flux as a function of energy taken in the turbulent region of the simulation, and the yellow line represents a power-law fit with an index of −4.6, for which we only use ion energies up to 200 keV.The blue squares and red squares represent MMS observations extracted from the event in Figure1; the red squares fit a power-law index of −3.9.The dashed green lines show the ion flux for a simulated x-line twice as long (4R E ).With respect to the red line, the high-energy flux increases by more than a factor of 3 when the x-line is doubled in size.

Figure 4 .
Figure 4. Top: evolution of the temperature inside the turbulent region.The simulation is initially loaded with particles at 4 keV and particles are injected at the boundaries also at 4 keV.The temperature quickly reaches a steady state at 14 keV.Bottom: density profile in the x-direction for the baseline simulation (normalized to n sheet ).

Figure 5 .
Figure 5. 2D histogram of the change in energy (gain or loss) over 0.1 s of each ion in the baseline simulation, as a function of their position in z.The vast majority of ions see a change in energy of less than 0.1 keV in 0.1 s, which appears as a horizontal red bar.However, a subset of ions see strong changes in energy, both positive and negative, which indicates a strong stochastic contribution from turbulence.When averaged, the net energization is positive (purple trace).

Figure 6 .
Figure 6.Types of orbits extracted from the simulation for the no-turbulence case.The background magnetic field is drawn as dashed black arrows, while the background electric field (E y ) appears as blue dashed arrows.(a) blue trace: a common trajectory of an ion jet, which enters a z face and leaves through the x face.(a) Violet trace: Fermi-reflected ion which enters and exits the +x face.(b) Green trace: an x-entering Speiser orbit, which leaves out of the +y face after crossing the z = 0 plane and engaging in a serpentine motion.(b) Red trace: a jet-forming ion that enters the −y-axis and leaves through the +x face.(c) Magenta trace: a trajectory that drifts toward z = 0, due to the E × B drift, engages in serpentine motion, and exits the +y face.(c) Orange trace: serpentine trajectory (yellow) moving parallel to the field E y and crossing the z = 0 plane multiple times.

Figure 7 .
Figure7.Left: normalized log-linear plot of the exiting energies at the +y face over a 1000 s run.Particles that cross z = 0 tend to partake in Speiser-like orbits and exit with higher energy than those that do not cross the magnetic null.Right: the bimodal distribution, results from turbulent δE and δB acting on the ions that cross z = 0.The green histogram only counts the particles that cross the z = 0 plane at least once.The magenta histogram counts the particles that cross z = 0 only once.The observed energization of low-energy peaks (red traces) is from 8-9 keV.Adding the turbulent results in a a similar low-energy profile profile, but the extra population of heated particles emerges.

Figure 8 .
Figure 8.(a) Energization of all particles leaving the face as a function of their displacement in the y-direction in a 1000 s run.The black and blue dots represent the ions in the no-turbulence and turbulence runs, respectively.The linear trend due to the electric field E y is still present in the energy gained by the ions affected by turbulence (after Figure 4 in Dolgonosov et al. 2013).(b) Exiting energy out of the simulation box as a function of turbulence amplitude.The turbulence amplitude is parametrized as the rms of the simulated δB and δE, over the observed rms values B rms = 8.1 nT and E rms = 18 mVm −1 , respectively.An ordinate value of 1 signifies that we have used the observed value of the rms in the simulation.We see that the energy outflow increases quadratically with turbulence amplitude.The orange line represents the inflowing Ponyting flux into the region, computed at the boundaries of the simulation box where the turbulent field vanishes.The blue line represents the fraction of the energy leaving in low-energy ions (>80 KeV), while the blue line represents the fraction of high-energy ions leaving the simulation.

Figure 9
Figure9.Top: two simulated ion trajectories without turbulence.Blue-yellow trace: Speiser-like trajectory that drifts toward the magnetic null plane and then is directed toward the +y face.Green-red trace: a particle that contributes to the ion jet that follows the hyperbolic magnetic field and exits the +x face.Bottom: the same ions but with turbulence.Blue-yellow trace: a nonthermally accelerated ion.When the particle crosses the magnetic null plane, it engages in a Speiser-like orbit, and the dwell time and energization increase.The energized particle is redirected toward the +y face, due to the background electric field.Green-red trace: the ion is redirected toward the +y face, due to turbulence, and the dwell time decreases.

Figure 10 .
Figure 10.Exiting fluxes for each face of the simulation, with and without turbulent fields.The red-shaded region indicates the >80 keV energies of the flux-energy distribution.A tail forms in the +y face when we apply the turbulent fields.Moreover, the reconnection jets get energized as well.
. As done in this work, Ergun et al. (2020b) and Ergun et al. (2022) employed realistic δB and δE signals (in 1D and 3D) to study electron acceleration in Earthʼs magnetosphere.The heart of the code has a relativistic Boris advancement algorithm (Boris 1970), which has excellent energy conservation and can accept input of realistic B and E signals.The code in this article was tested for long simulation periods and displayed less than a few percent change in energy.A.1.Simulation Domain, Boundary Conditions, and Background Magnetic Field

Figure 11 .
Figure 11.(a) The measured (black line) and the reproduced (blue line) δB spectra.(b) The measured (black line) and the reproduced (red and blue lines) δE spectra.(c) The measured speeds of δB vs. frequency derived from time delays from the four MMS spacecraft.The blue dots represent measurements.The black squares denote the averages; the red line indicates the fit of these averages.The formula denoted by the red line is used in the simulation.(d) The measured correlation times of δB signals vs. frequency derived from the four MMS spacecraft (same legend as in the δE case).