Proton and heavy ion acceleration by stochastic fluctuations in the Earth ’ s magnetotail

Spacecraft observations show that energetic ions are found in the Earth’s magnetotail, with energies ranging from tens of keV to a few hundreds of keV. In this paper we carry out test particle simulations in which protons and other ion species are injected in the Vlasov magnetic field configurations obtained by Catapano et al. (2015). These configurations represent solutions of a generalized Harris model, which well describes the observed profiles in the magnetotail. In addition, three-dimensional time-dependent stochastic electromagnetic perturbations are included in the simulation box, so that the ion acceleration process is studied while varying the equilibrium magnetic field profile and the ion species. We find that proton energies of the order of 100 keV are reached with simulation parameters typical of the Earth’s magnetotail. By changing the ion mass and charge, we can study the acceleration of heavy ions such as He and O, and it is found that energies of the order of 100–200 keV are reached in a few seconds for He, and about 100 keV for O.


Introduction
One of the unsolved problems of magnetospheric plasma physics concerns the generation mechanisms of energetic electrons and ions.Energies from tens of keV to a few hundreds of keV are found in the Earth's magnetotail (Christon et al., 1989;Keiling et al., 2004;Imada et al., 2007;Haaland et al., 2010;Artemyev et al., 2014).Heavy ions like singly charged oxygen O + are also observed in the magnetotail, with energies often reaching several hundred keV (Keika et al., 2013;Kronberg et al., 2014Kronberg et al., , 2015)).Further, relevant levels of energetic electrons and ions are also observed during periods of low geomagnetic activity (Christon et al., 1989), leaving open the search for the acceleration mechanism.In addition, the observed proton temperatures, of the order of 5-10 keV, are also much larger than the energy corresponding to their original source, but the mechanism that can generate such heated particles is not fully understood (Runov et al., 2006;Artemyev et al., 2011).Grigorenko et al. (2009) pointed out that 100 keV ions are accelerated near reconnecting regions, possibly by a process related to time-dependent, patchy reconnection.Spacecraft observations suggest that both electrons and ions are accelerated not only in the vicinity of the reconnection X line but also in a larger area around the reconnection region (e.g., Imada et al., 2007;Grigorenko et al., 2009;Ono et al., 2009).Indeed, during disturbed periods, strong electric (Cattel and Mozer, 1982), velocity and magnetic fluctuations (e.g., Hoshino et al., 1994;Borovsky et al., 1997;Zimbardo et al., 2010), as well as energetic ions, are observed in the magnetotail.
Also, for energetic particles in solar flares, only a limited number of particles are believed to be accelerated directly by the reconnection electric field (e.g., Ambrosiano et al., 1988).Therefore, a second process, in which either firstorder Fermi acceleration (e.g., Tsuneta and Naito, 1998) or stochastic Fermi acceleration (e.g., LaRosa et al. , 1996) is involved, is also required.The possible ion heating and acceleration due to reconnection jets in the solar corona was Published by Copernicus Publications on behalf of the European Geosciences Union.
considered by Artemyev et al. (2014), where heavy ions were studied too.
With regard to the geospace environment once more, ion acceleration at dipolarization fronts in the geomagnetic tail has been studied, among others, by Ashour-Abdalla et al. (2011), Ukhorskiy et al. (2013), Birn and Hesse (2014), Greco et al. (2014) and Greco et al. (2015).Perri et al. (2009) and Greco et al. (2010) considered the combined effect of a steady-state dawn-dusk electric field and of stochastic Fermi acceleration due to the presence of transient magnetic structures: by performing a two-dimensional (2-D) test particle simulation, they showed that proton energies of up to 80-100 keV can be reached, in agreement with the above observations.Further, Perri et al. (2011) developed a threedimensional (3-D) model which takes into account the equilibrium magnetic structure of the current sheet (CS).A Harris-like profile of a magnetic field reversal with the presence of a normal magnetic field component was used.It was found that the stochastic Fermi acceleration is efficient to explain the acceleration of protons, although somewhat smaller energies than in the 2-D case are obtained.
While the Harris magnetic field profile is very popular in space plasmas, many observations in the Earth's magnetotail show that the current sheet has a non-Harris-like profile, with the current sheet being often embedded in a thicker plasma sheet (e.g., Runov et al., 2006).Recently, Catapano et al. (2015) generalized the well-known solution of the Harris current sheet to the case where several current carrying populations, i.e., multiple electron and ion populations, are present.Those solutions allow for adjustment of the temperature profiles and the density profiles of the plasma sheet, in a wide range of configurations, with the magnetic field profile being obtained self-consistently.A comparison with data from the Cluster spacecraft shows that different current sheet crossings can be well reproduced by the new solutions (Catapano et al., 2015).Thus, it is interesting to understand influence of the new magnetic field profiles on the particle acceleration process for different ion species, also in view of the fact that the current sheet thickness and the current profiles are observed to change significantly during magnetospheric activity (Hoshino et al., 1996;Sergeev et al., 2003;Vasko et al., 2015).
In this work we carry out test particle simulations in which protons and heavier ions are injected in the magnetic field configurations obtained by Catapano et al. (2015).The 3-D stochastic electromagnetic perturbations of Perri et al. (2011) are included in the model, so that the ion acceleration process is studied while varying the equilibrium magnetic field profile, the fluctuation level, and the ion species.By changing the ion mass and charge, we can study the acceleration of protons and heavy ions such as He ++ and O + .
In Sect.2, we present the self-consistent Vlasov solutions for the generalized Harris current sheet.In Sect. 3 we present the numerical model adopted to describe the magnetotail electromagnetic perturbations and discuss the numerical re-sults on acceleration of particles, both for protons and for heavier ions.In Sect. 4 we discuss the energy gain for the different ion species, and finally, in Sect.5, we give our conclusions.
2 Self-consistent current sheet model Catapano et al. (2015) developed a current sheet model which allows for regulating the level of plasma temperature and density inhomogeneities across the sheet.These models generalize the classical Harris model (Harris, 1962) via including two-temperature current-carrying plasma populations and one background plasma population not contributing to the current density.We use the notation f j v j , T j to describe a Maxwellian distribution function with the corresponding drift velocity v j and temperature T j , where j indicates the species (j = i for ions and j = e for electrons).Any linear combination of f j is solution of the stationary Vlasov equation.In the velocity plane the Maxwellian distribution function can be presented as a peak of phase space density with contour lines represented by concentric circles, centered at the origin.The shifted Maxwellian has the same shape in the velocity plane, but the center is shifted from the coordinate origin to the distance equal to the drift velocity.The spacecraft observations show that the phase space density of the particles velocity distribution in thin current sheet is often represented by a ring distribution (see Artemyev et al., 2009).The simplest way to construct such a distribution function consists in using a combination of two Maxwellian distributions with different temperatures.Therefore, to develop our model, we use the following distribution function: where the parameter f ring represents the ring distribution, so that F j represents a linear combination of ring distribution and background plasma.The ring distribution is obtained by the combination of two shifted Maxwellians with different temperatures, where the parameter γ describes the relative density of the second current-carrier population.f j (0, T 0j ) represents the background plasma with zero bulk velocity.The parameter δ represents the relative density of the background population and its value is chosen to have a positive phase space density everywhere in the phase space.For γ = 0 we obtain the distribution function presented in Yoon and Lui (2004).
We consider the coordinate system with the x axis directed along Sun-Earth direction, the z axis oriented across the CS (i.e., along the direction of the CS inhomogeneity), and the y axis directed along the electric current flow.Thus, in this model, the magnetic field has one component B x (z) corresponding to the vector potential A y (z), while the electric field has only the component E z (z) corresponding to the scalar potential φ(z).The Maxwell equations reduce to equations for the scalar potential φ and vector potential A y as in the classic Harris equilibrium.Using the quasi-neutrality condition (e.g., ∂ 2 φ/∂z 2 ∼ = 0) we obtain the following equation for the dimensionless vector potential (a = A y q i v 1i /cT 1i , where q i is the ion charge and c is the speed of light and where the Boltzmann constant is set to 1), and the dimensionless scalar potential ( = φq i /T 1i ), where the function G, obtained from the quasi-neutrality condition, is (Catapano et al., 2015) Above, t e , t i , and t b are the normalized temperatures with respect to T 1i , of electrons and ions of first, second and background populations.The most important parameter is , where we assume that We notice that T i ∼ 5T e most frequently in the magnetotail (e.g., Baumjohann et al., 1989;Runov et al., 2006), and this condition has an influence on the current carrying populations.Here, the condition T 0e = T 0i can be used because the background populations do not contribute to the current density.Usually, electrons are colder than ions and we use this condition for the first and second populations that are carrying the current.
For α = 1 we obtain the class of solution proposed by Harris with the null scalar potential.We have investigated different classes of solutions varying the parameters δ, γ , α and t b to describe the temperature inhomogeneity in the CS.At variance with Harris solutions, an analytical solution for the vector potential is not possible in general, so that we need to carry out a simple numerical integration along z of Eqs. ( 2)-( 4) to obtain the various quantities (Catapano et al., 2015).Several populations of ions and distribution of scalar potential φ(z) result in realistic distribution of plasma parameters across the magnetotail.To illustrate the usefulness of the our CS model, we show the comparison of modeled and observed distributions of ion temperature.From the Cluster database we picked up four crossings with distinguished variations in the ion temperature across the sheet (see Fig. 5 in Artemyev et al., 2011).In Fig. 1 (adapted from Catapano et al., 2015) we report the comparison of T i, eff as a function of B x spacecraft observations showing the model capability of describing different configurations.The ratio of pressure and density gives the effective ion temperature T i, eff .Indeed, one of the basic features of the Harris solution is that the temperature is constant across the CS thickness.By changing the parameters, we obtain different profiles of the magnetic field B x (z) and electric field E z (z).In this work we use the two self-consistent solutions that better describe observations, as shown in Fig. 1.In particular, we focus on profiles (c) and (d) of Fig. 1, and the corresponding magnetic and electric field profiles are shown in Fig. 2.These two solutions are those which have the lowest and the highest slope, respectively, in the magnetic field profile among solutions in Fig. 1.In Fig. 2 the black line represents the Harris solution with the null electric field, the red line represents the magnetic field profile that has the lowest slope, or slower variation in the magnetic field across the CS, and the blue line represents the profile with the highest slope.The parameters used to obtain the profiles in Fig. 2 are shown in Table 1.
3 Numerical simulations In the test particle numerical simulation, we overimpose time-dependent electromagnetic fluctuations (Perri et al., 2011), which represent fluctuations frequently observed in the magnetotail even in quiet periods (Borovsky et al., 1997), on the self-consistent solutions shown in Fig. 2. We use a three-dimensional simulation box with L x = L y = 10 5 km and a size along the z direction that ranges from −L z to L z , where L z = 2.5 × 10 4 km.This has to be considered a local simulation box; indeed, its size of about 15 R E is only a fraction of the actual magnetotail extension, so that large-scale variations in the magnetotail are neglected.The characteristic thickness of the CS is set to λ = L z /5 = 5 × 10 3 km.We use the coordinate system as in Catapano et al. (2015).N = 25 time-dependent electromagnetic fluctuations are located in the (x, y) plane; such perturbations create an oscillating motion in the plane.Thus, the total magnetic field is given by where the first term is the magnetic field from self-consistent numerical solutions of Eq. ( 2) (shown in Fig. 2), the second term is the out-of-plane magnetic field, B n , that simulates a remaining part of the Earth dipole magnetic field.The third term represents the time-dependent fluctuations, δB(r, t) = × A(r, t).The vector potential perturbations are modeled as (Perri et al., 2011) where r is the position of a particle, r i (t) is the position of the fluctuation center at time t, and the sum is made on the number of fluctuations i.The parameter represents the decreasing scale of the vector potential, which we can imagine as the size of a magnetic cloud that interacts with the particle (see Perri et al., 2007).The positions r i are fixed randomly in the (x, y) plane and are oscillating with velocity V = 400 km s −1 (i.e., the typical value for the Alfvén speed in the magnetotail; Hoshino et al., 1994;Nakamura et al., 2004;Vörös et al., 2007) and with random phases.The electric field is obtained as E = − φ − ∂A(r, t)/∂t.Thus, we have with E 0y the constant dawn-dusk electric field and where E z (z) is the electric field coming from the self-consistent solutions of Eq. (3) (see Fig. 2), and δE(r, t) = ∂A(r, t)/∂t represents the time-dependent fluctuating term of the electric field coming from Eq. ( 6).All the equations are normalized using the following normalization quantities: L = 10 5 km, B 0 = 2 nT and E 0 = 40 mV m −1 (see Perri et al. (2011) for more details).Further, the fluctuating potential decrease scale is set to = 8000 km (e.g., Nakamura et al., 2004;Greco et al., 2010;Perri et al., 2011); the amplitude of electromagnetic fluctuations, A 0 / , is of the order of 10 nT, which is consistent with observations (Borovsky et al., 1997); the dawn dusk electric field is E 0y = 0.2 mV m −1 ; the normal component of magnetic field is B n = 3 nT; and the asymptotic value of the magnetic field in the lobes is B max = 20 nT = 10 B 0 (Sergeev et al., 2003).The magnetic field B x (z) reaches its maximum (minimum) value for z = L z (z = −L z ), but most of the variation happens on scale λ.We simultaneously inject N p = 10 4 particles in the simulation box at t = 0, at z = 0, and randomly distributed in the (x, y) plane; the starting coordinates x and y vary from 0 to L, where L is the size of the simulation box.Particles exiting the simulation box are substituted with freshly injected particles.The initial velocities for protons are extracted from a Maxwellian distribution with v th = 120 km s −1 ; heavier ions are injected with different thermal velocities in order to have the same initial temperature of protons.We numerically solve the equations of motion, for each particle via a fourth-order Runge-Kutta scheme.The integration step is fixed to t = 0.001 −1 p , while t 0 = −1 p = 5 s.

Results for protons
To study the particle dynamics we have integrated the trajectories of a few particles using the profiles 1 and 2, reported in Table 1. Figure 3 shows the projections of trajectories in two different planes of the box for profile 1 and the energy gained  by particles as a function of time (upper panels) and the trajectories and the energies as a function of time for profile 2 (lower panels).We can notice that the parabolic structure of the large-scale magnetic field forces the particles to leave the box along the positive x direction (see, for example, the middle panels in Fig. 3) and the ∇B drift along y determines the quasi-cycloid orbits.At z ∼ 0, when the magnetic field B x (z) reaches the zero value, particles are less magnetized and the probability of being accelerated by the time-dependent fluctuations increases.Indeed, in this region particles experience frequent interactions with the electromagnetic fluctuations in the CS via a stochastic Fermi process (see, for example, the blue trajectory in Fig. 3).Particles can also be accelerated by the constant E 0y electric field.It is interesting to note that, as a result of the interaction with magnetic fluctuations, particles show a chaotic behavior and exhibit a meandering motion (see red lines in the middle panels in Fig. 3).Particles undergoing Fermi acceleration reach energies of up to 120 keV.Between the simulations with profile 1 and 2 we observed some different behavior in the individual particle trajectories (see, for example, particle 2 in the lower panel).The starting points, denoted by triangles in the figure, are the same as the case with profile 1.Although there are differences among particle orbits, they have no influence on the statistical behavior of energization.Results from a statistical analysis are shown in Fig. 4, which displays comparison among the probability density functions (PDFs) of the proton energies at t = 500 t 0 for profile 1, profile 2, and injection in the case with A 0 / = 10 nT (large A 0 ), A 0 / = 3 nT (medium A 0 ) and A 0 / = 1 nT (small A 0 ).In all of these three cases the protons have the same initial distribution (red line in Fig. 4).As expected, the particle energization increases with the perturbation strength and protons gain energies of up to 100 keV in the case with large A 0 .This means that the electromagnetic perturbations create an efficient accelerator for protons in the magnetotail (also shown in Perri et al., 2009;Greco et al., 2010).On the other hand, the influence of the magnetic field profile B x (z) is very small, as shown by the PDFs being nearly overlapped.It is interesting to note that, for A 0 / = 1 nT, the PDF(E) reaches energies of 6.5 keV.
Considering that the proton temperature in the magnetotail plasma sheet is observed to be of the order of 5 keV (e.g., Borovsky et al., 1997;Runov et al., 2006;Artemyev et al., 2011), a very small fluctuation level with δB ∼ 1 nT, almost always present in the magnetotail (Hoshino et al., 1994), is  sufficient to heat the cold plasma coming from the solar wind to a temperature of 3-5 keV (Runov et al., 2006).

Results for heavy ions, He ++ and O +
A completely new study is now performed by integrating heavier ions' trajectories within the above simulation model.The equations and the numerical integration are the same as the ones we used for protons, except that we change the ratio q/m in Eq. ( 9).We perform numerical investigations using He ++ and O + ions and compare the cases with the profile and 2 of B x (z) (see Fig. 2). Figure 5 reports trajectories of He ++ and O + ions, and energies vs. time for profile 1.The black lines represent the He ++ ions and red lines represent the O + ions.We use the solid lines for particle 1, dashed lines for particle 2, and dashed-dotted lines for particle 3. Again, the triangles denote the starting point of the ions.At injection, the He ++ ions have a Larmor radius r He++ L = r L p ∼ 600 km, which is 0.078 times the size of fluctuations .Because of the equal Larmor radius, He ++ ions interact with the fluctuations as protons do, but they can reach larger energy because of their double charge.This effect is discussed extensively in Sect. 4. The Larmor radius of O + is r L O+ = 4r L p ∼ 2400 km, which is ∼ 0.3 times .Because of the larger Larmor radius, the electromagnetic fluc-tuations have a different effect on O + ions.For example, we can compare the trajectories represented by the solid black line (for He ++ ) and red dashed line (for O + ) in Fig. 5.In the left panel it is shown that these two particles of different species have similar trajectories.However, if we observe the energies vs. time (middle and right panels), we notice that the He ++ ions have been energized up to 300 keV, while O + ions reach an energy of around 75 keV in a few interactions corresponding to stochastic Fermi acceleration.Results for profile 2 are not reported since they are similar to the ones for profile 1. Figure 6 shows the PDFs of particle energy of He ++ and O + ions for different values of A 0 / .As for protons, the energization grows with the perturbation strength A 0 .We can see that in the case with A 0 / = 10 nT, He ++ ions reach energies as large as 200 keV, while O + ions reach energies of about 100 keV.

Energy gain
Figure 7 shows the PDFs of protons along with heavy ions in the presence of profile 1 and A 0 / = 10 nT: we can see the different effect of the electromagnetic fluctuations on the different species, with He ++ being the most energized particles.It is interesting to estimate the rate of energy gain, W/ t, for each species.A rough a priori estimate of W/ t can be obtained by considering the work done on a particle as L = qδE • s.Assuming that δE = V δB, we can obtain W/ t = qδBV v, with V the perturbation velocity and v the particle speed.In order to compare the results, consider the energy W (t i ) at each step of integration t i of N = 500 particles and we calculate the corresponding energy gain: The stochastic Fermi-like process implies multiple interactions with fluctuations that can produce an increase or a decrease in energy (see the right panel in Fig. 3 and the middle and right panels in Fig. 5).For each interaction, if ( W/ t) i < 0, particles lose energy; conversely, if ( W/ t) i > 0, particles gain energy.To estimate the rate of energy gain, we take into account only the energy gains larger than 1 keV s −1 .We count the energy gain that exceeds this threshold for each particle and then we average over the ensemble of particles.The average value of energy gain for protons is found to be around 3.52 keV s −1 , while for He ++ ions it is 4.62 keV s −1 , that is, it is similar for the two species.This is consistent with the fact that these two species have the same Larmor radius and perform similar interactions with the time-dependent electromagnetic fluctuations.Helium has a double charge but a smaller average speed v; thus, these effects compensate for each other in the estimate of W/ t ∼ qδBV v.For O + ions the average energy gain is of 1.98 keV s −1 , that is smaller than protons and He ++ .Indeed, they have a larger Larmor radius but a single charge and a smaller speed, so that the interaction with the fluctuations is less efficient.It is interesting to note that the proton acceleration rates are larger than other acceleration mechanisms, like those considered at shock fronts in the heliosphere (Perri andZimbardo, 2012, 2015;Zimbardo and Perri, 2013).
Previous studies have pointed out that, in the near-Earth magnetotail, heavy ions as O + are accelerated more efficiently than protons, especially during periods of strong geomagnetic activity.Generally, this is observed during substorms or strong reconnection, as has been reported in the literature (e.g., Nosé et al., 2000a, b;Ono et al., 2009;Keika et al., 2013;Luo et al., 2014).Since the energization process has been observed to happen during substorms and localized dipolarization of the magnetic field, the acceleration process can be considered "local".O + ions can be accelerated up to 100 keV and more.In our numerical model, instead, only the He ++ ions exceed the 100 keV.This is probably due to the different acceleration mechanism: particles continuously interact with stochastic fluctuations and diffuse in the simulation box, and the perturbed magnetic field is at most 10 nT; conversely, in dipolarization fronts the peak magnetic field δB can be as large as 20-40 nT (Ono et al., 2009;Sergeev at al., 2009), and this is clearly influencing the acceleration rate.

Conclusions
In this paper we have investigated the dynamics and the acceleration of protons and heavier ions in a CS model that includes transient electromagnetic perturbations.The equilibrium magnetic field profile is obtained from a new class of self-consistent solutions of the Vlasov-Maxwell equations which extends the Harris equilibrium to the case where several current carrying populations, with different temperatures and bulk velocities, are present (Catapano et al., 2015).In particular, we have chosen three different solutions for the magnetic field profile across the CS: one corresponding to the Harris sheet, one shallower than the Harris sheet (profile 1), and another one steeper (profile 2).We choose to use these profiles obtained by a self-consistent solution of the Vlasov-Maxwell equations system, because these can well reproduce the observed profile of temperature in the CS (see Fig. 1).The asymptotic lobe magnetic field is kept constant at B max = 20 nT, and the normal magnetic field is B n = 3 nT.The three-dimensional electromagnetic perturbations are obtained from a model developed by Perri et al. (2011), in which the size and the strength of the perturbations can be changed.By injecting ions in a local simulation box in Cartesian coordinates, we have studied the acceleration of protons, He ++ ions, and O + ions.
Regarding the numerical results for protons, we find that the different equilibrium magnetic fields affect the particle dynamics in a similar way.The degree of energization is determined by the amplitude of the electromagnetic fluctuations.This means that when the perturbations have a large amplitude, their influence dominates and the acceleration process is not very sensitive to the equilibrium magnetic field profile.For both profiles, proton energies of up to 100 keV are obtained, in agreement with observations.When a small magnetic fluctuation level is present, A 0 / 1 nT, the proton F. Catapano et al.: Proton and heavy ion acceleration by stochastic fluctuations in the Earth's magnetotail energy can easily reach 5-6 keV and He ++ and O + can reach 10 keV, in agreement with typical ion temperatures observed in the magnetotail (e.g, Baumjohann et al., 1989;Runov et al., 2006;Artemyev et al., 2011).For He ++ ions, notably, energies of up to 200 keV are obtained (see Fig. 6).These energies are reached in 5-10 s, indicating a high efficiency of the acceleration process.We consider that, for He ++ ions, the dynamics are very similar to that of protons, but the double charge allows for stronger and faster acceleration by the fluctuating electric field, as also shown by the estimate of the average energy gain W/ t.
For O + ions we get similar results as far as the influence of profile 1 and profile 2 is concerned, although the Larmor radius is now substantially larger and the trajectories look less magnetized.However, the energy gained by O + ions is lower than for helium, and it is below 100 keV.We can argue that the larger oxygen Larmor radius modifies the particle interaction, decreasing the energy gain.This effect might also be due to the fact that, just because of the larger Larmor radius, oxygen ions spend less time in the quasi-neutral sheet (z ∼ 0), where the electromagnetic fluctuations are the strongest.In order to check this possibility -i.e., in order to investigate whether the ion energization depends on the finite size of the simulation box, namely on the residence time in the simulation box, we have performed additional simulations (not shown).In these simulations we set the normal component of the magnetic field B n to 6 nT instead of 3 nT.The results show that all the species can reach higher energy in this configuration but that He ++ remains the most energized species.In order to understand the cause of this behavior, we estimate the residence time in the simulation box of each ion species, for different values of B n and of the amplitude of fluctuations.In all the cases we found that O + ions spend a longer time in the box (also due to the lower speed) but that He ++ is the most energetic species.We can conclude that the numerical results are not due to a finite size effect but to the stochastic interaction.
It is interesting to compare our results with spacecraft observations in the magnetotail.In the presence of substorms and local dipolarization of magnetic field, the most energized species is found to be oxygen, reaching energies larger than 200 keV (Nosé et al., 2000a, b;Ono et al., 2009;Luo et al., 2014).During those periods the magnetic fluctuations are generated by reconnection jets and produce a spatially localized acceleration of ions.Those observations are in agreement with the results obtained by Greco et al. (2015), where multiple ion species are interacting with a single, high-speed dipolarization front corresponding to a reconnection jet, and where it is found that the energy gain grows with the ion mass, with oxygen reaching more than 200 keV.Also, Grigorenko et al. (2015) carried out a numerical study of ion dynamics in a plasmoid-like configuration with the presence of electromagnetic turbulence, and it is found that ion energization also depends on the resonant interaction of ions and wave harmonics: with strong enough waves, O + ions can be rapidly accelerated up to several hundreds of keV.
The acceleration process considered in the present work turns out to be more efficient for He ++ ions than for protons and oxygen ions, and therefore it is not able to explain the most energetic O + observations.The main difference between the mentioned studies and the presented model consists in the particle interaction with multiple magnetic perturbations moving "randomly".In our model the acceleration is not localized as for dipolarization fronts, and ions are diffusing while interacting with the transient fluctuations.The proposed model is more similar to a stochastic Fermi acceleration process and is probably more relevant for quiet geomagnetic periods in the magnetotail (Borovsky et al., 1997).We would also like to stress that a similar stochastic process can be at work jointly with other mechanisms; a multi-step acceleration process (high-speed dipolarization front and Fermilike acceleration) will be investigated in a future work.

Data availability
Data supporting all the figures reported in this work come from test particle simulations performed on computers at the Physics Department, University of Calabria.Data can be accessed by writing to the following address: filomena.catapano@unical.it.

Figure 1 .
Figure 1.Comparison of four profiles of ion temperature measured in the magnetotail CSs (shown by blue circles) with model profiles (shown by black curves).See Catapano et al. (2015) for details.

Figure 2 .
Figure 2. Model profiles of B x (z) and E z (z), the black line represents the Harris solution, the red line the solution with lower slope in B x , and the blue line the profile with higher slope.

Table 1 .
Parameters used to obtain the profiles in Fig.2

Figure 3 .
Figure 3. Trajectories and kinetic energies of protons in the presence of B x (z) as in profile 1 (upper panels) and as in profile 2 (lower panels).The fluctuation level is set to A 0 / = 10 nT.The colored lines represent the particle orbits and their starting points are denoted by triangles.Axes not to scale.Particles are injected at t = 0 randomly within the (x, y) plane and at z = 0. Their initial velocities are extracted from the same Maxwellian distribution (see text for further details).

Figure 4 .
Figure 4. Comparison among the PDFs of energies of protons in the presence of profile 1 and profile 2, by varying A 0 / = 10 nT (large A 0 ), A 0 / = 3 nT (medium A 0 ) and A 0 / = 1 nT (small A 0 ).

FFigure 5 .
Figure 5. Trajectories and energies of He ++ and O + ions in the presence of B x (z) as in profile 1. Black lines represent the He ++ ions and red lines the O + ions (solid line for particle 1, dashed line for particle 2, dashed-dotted line for particle 3).The triangles denote the starting point.

Figure 6 .
Figure 6.Comparison of the PDFs of energies of O + and He ++ ions in the presence of profile 1 and profile 2.

Figure 7 .
Figure 7. PDF of protons and heavy ions in the presence of profile 1 and A 0 / = 10 nT.