Nonlinear Processes in Geophysics Acceleration and transport of ions in turbulent current sheets : formation of non-maxwelian energy distribution

The paper is devoted to particle acceleration in turbulent current sheet (CS). Our results show that the mechanism of CS particle interaction with electromagnetic turbulence can explain the formation of power law energy distributions. We study the ratio between adiabatic acceleration of particles in electric field in the presence of stationary turbulence and acceleration due to electric field in the case of dynamic turbulence. The correlation between average energy gained by particles and average particle residence time in the vicinity of the neutral sheet is discussed. It is also demonstrated that particle velocity distributions formed by particle-turbulence interaction are similar in essence to the ones observed near the far reconnection region in the Earth’s magnetotail.


Introduction
The observation of non-maxwellian particle energy distributions in the Earth's magnetosphere has a long history (Vasyliunas, 1968;Sarris et al., 1976;Christon et al., 1989).In particular, the observed distributions can be quite well approximated by the so called kappa-distribution (ε is energy): (1) Due to wide occurrence of such a distribution in nature its structure has been investigated in a number of papers on current sheet (CS) modeling and CS stability (Fu and Hau, 2005;Yoon et al., 2006;Hu et al., 2008).Theoretical investigation of the relationship between the system's entropy and anomalous diffusion with kappa-distribution (Milovanov and Zelenyi, 2000;Leubner, 2004Leubner, , 2008) ) has deepened the understanding of transport processes occuring in astrophysics and Correspondence to: A. V. Artemyev (ante0226@gmail.com)plasma physics.However, the problem of the mechanisms responsible for particle acceleration and kappa-distribution formation is still unsolved.Some mechanism of particle acceleration in limited spatial regions (such as planetary magnetospheres and the solar corona) like the famous Fermi acceleration mechanism (Fermi, 1949) is needed to explain the formation of power tail (f ∼ε −1−κ ) of energy distribution and presence of groups of particles with large energy values (ε ε ).Different possible mechanisms were proposed such as acceleration near the CS X-line (Hoshino, 2005;Pritchett, 2006;Drake et al., 2006), ion acceleration by electric field due to the growth of tearing instability (Zelenyi et al., 1984;Taktakishvili et al., 1998), quasiadiabatic ion acceleration in the vicinity of the magnetotail neutral sheet by the dawn-dusk electrostatic electric field (Speiser, 1967;Lyons and Speiser, 1982;Ashour-Abdalla et al., 1993, Litvinenko andSomov, 1993;Vainchtein et al., 2005;Zelenyi et al., 2007), particle acceleration by MHD turbulence in the solar corona (Kobak and Ostrowski, 2000;Dmitruk et al., 2004), acceleration due to dipolarization in the Earth's magnetosphere (Delcourt and Sauvaud, 1994;Delcourt, 2002;Apatenkov et al., 2007;Ono et al., 2009).In this paper we suggest another possible mechanism of acceleration which can operate in dynamic regions of CS.This mechanism assumes the presence of turbulent electromagnetic field (TEMF) in the central region of CS (Cattel and Mozer, 1982;Hoshino et al., 1994;Bauer at al., 1995) that can effectively interact with the charged particles and energize them.

Models of turbulent electromagnetic field
In this section, a brief description of the existing TEMF models is presented and the general equations of the considered model are discussed.
A. V. Artemyev et al.: Acceleration and transport of ions in turbulent current sheets MHD models of TEMF are often used to simulate the conditions of particle acceleration in the solar corona (Galsgaard and Nordlund, 1996;Kobak and Ostrowski, 2000;Dmitruk et al., 2004).In such models, TEMF is presented as a result of independent numerical modeling and the process of acceleration is described for test particle ensembles.Sometimes the magnetic field of CS is taken as an initial condition (Dmitruk et al., 2004;Onofri et al., 2006) and models of stationary reconnection with background MHD turbulence can be developed as well (Kobak and Ostrowski, 2000).Spatial scales of MHD structures lead to considerable difference between ion and electron acceleration mechanisms.Electrons mainly gain energy moving along field lines while ions are usually accelerated as a result of perpendicular drift (Dmitruk et al., 2004).
Some recently published papers have been devoted to particle acceleration by means of interaction with spatially localized structures, namely, magnetic clouds (Perri et al., 2007(Perri et al., , 2009)).The model used in these papers could be considered as the electromagnetic analogs of the stochastic billiards' mechanical model (Karlis et al., 2007 and references therein).
Models of electron acceleration by means of conservation of adiabatic invariants during sufficiently slow magnetic field variation can also be applied either to the solar corona (Somov and Kosugi, 1997;Bogachev and Somov, 2007) or to the Earth's magnetosphere (Smets et al., 1999;Apatenkov et al., 2007).The essence of these methods is the following: provided that the electron magnetic moment µ=T ⊥ B and the second invariant T L 2 (L is the length of a flux tube) are conserved, electron perpendicular T ⊥ and parallel T temperatures change in the course of evolution of magnetic field B.
The mechanism of particle acceleration by electrostatic field in CS geometry (Speiser, 1967;Lyons and Speiser, 1982) may also be generalized in case of CS with stationary turbulence produced by a wave ensemble with ω k = 0 (Veltri et al., 1998;Greco et al., 2002;Zimbardo et al., 2004).Such a turbulence is characterized by the absence of inductive electric field.However, acceleration by an electrostatic field becomes more efficient due to modification of particle dynamics under the influence of magnetic turbulence.Spatial diffusion in the stationary turbulence model (SM) has been analytically studied in neutral sheet geometry (Chiaravalloti et al., 2006).It was found that, the power exponent value ν ( r 2 ∼ t ν ) can vary from ν≈0.6 (subdiffusion) to ν=1 (Brownian diffusion).
Another model of particle acceleration may be obtained in case of the TEMF created by propagating waves with ω k =0.Such dynamical models (DM) used to be discussed in 2-D neutral sheet geometry.For this geometry, either electrostatic (Carbone et al., 2004) or electromagnetic waves (Zelenyi et al., 2008a) have been taken into account.The model that we consider here takes into account the TEMF produced by an ensemble of propagating electromagnetic waves, which may be regarded as a generalization of earlier models.
At the nonlinear stage, the unstable modes (Lui, 2004) interact with each other and their growth saturates (Zelenyi and Milovanov, 2004).Tearing instability causes oscillations of the normal component of magnetic field in the CS (Zelenyi et al., 2008b and ref. there).Dynamics of TEMF is mainly connected with various drift modes like "sausage" (Buechner and Kuska, 1999;Silin et al., 2002), "kink" (Daughton, 1999;Zelenyi et al., 2009) and LHDI (Daughton, 2003).However, taking into consideration all the structural peculiarities of electromagnetic waves formed as eigenmodes of instability is not an easy matter owing to the facts that only the linear theory of such modes has been thoroughly worked out, and that all these modes have essentially nonlocal structure.
Thus, in this paper TEMF is simplistically modelled as an ensemble of plane waves: . (2) Here One can easily check that the condition ∇δB=0 is naturally satisfied.Initial phases φ 1 k and φ 2 k are chosen randomly.The frequency ω k is chosen for each harmonic so that all the waves possess the same phase velocity ω k k=v φ (here waves are assumed to be drift and v φ is less then thermal ion velocity).The wave magnitude δB (k) is defined by the following expression: Here l is the vector of correlation length.The power index α=7 8 was taken from earlier papers (Veltri et al., 1998;Zelenyi and Milovanov, 2004), where it had been defined on the basis of experimental observations (Hoshino et al., 1994;Petrukovich, 2005, including references).
Using the system (1) and the condition of free charges' absence one can find from Maxwell equations the components of electric field: Now we can set the systems' parameters and geometry.The constant background magnetic field is taken according to the Harris CS modified model (Harris, 1962;Lembege and Pellat, 1982): B=B 0 tanh z L e x +B z e z and the main dimensionless parameter of the system is b n =B z B 0 .Magnitude of turbulence is set by the dimensionless parameter δ= √ δB • δB B 0 .Wavenumbers are distributed in the range 2π L −1 [0.05, 4], and |l| = 0.5L.In the simulations, we launched five hundred harmonics into the system.The

Numerical simulation scheme and particle injection
In this section, we describe the modeling scheme of particle acceleration due to interaction with TEMF is described in this section.Maxwellian velocity distributions with nonzero shift v D along magnetic field lines are generated at the top and bottom edges of CS (z=±3L): (5) Such distributions are often used for the modelling of thin CS formation as well as for construction of analytic models (Zelenyi et al., 2000, including references).Particles penetrate into the central region of the CS along field lines, interact with the TEMF and subsequently leave the system.In order to understand better the properties of this system, one may first analyze the typical particle trajectories without TEMF (B=B 0 tanh z L e x +B z e z ).As shown in Fig. 2a, a particle approaches the CS neutral sheet, makes a half-turn around B z and leaves the CS.Such trajectory is called "Speiser" or transient orbit (Speiser, 1965).In addition to Speiser orbits in such "regular" configuration there are also fully trapped particles with circular orbits and quasi-trapped particles in CS (Speiser, 1965;Chen and Palmadesso, 1986;Buechner and Zelenyi, 1989;Zelenyi et al., 2000;Vainchtein et al., 2005).On the other hand, an electric field E y pointing duskward appears in the CS as a result of interaction of the Earth's magnetosphere with the magnetized flow of the solar wind.At equator in the distant magnetotail, a typical value for this electric field is of the order of 0.1 mV/m (Volland, 1978;Kan, 1990).After the first interaction with CS, a particle at Speiser orbit experiences an energy gain ∼ E y B z 2 (Lyons and Speiser, 1982;Ashour-Abdalla et al., 1993;Litvinenko and Somov, 1993).As can be seen in Fig. 2b, the presence of SM turbulence can substantially increase the particle residence time in the CS.As a result of longer particle drift in the direction of the electrostatic field E y , particles gain more energy (Veltri et al., 1998;Greco et al., 2002).This is one of the possible mechanisms of particle heating.However, here, we focus on particle energization due to the inductive electric field produced by DM turbulence as illustrated in Fig. 2c.The comparison of these two mechanisms is discussed below.

Influence of accelerated particles on CS structure
Using a test particle method, one may estimate the current density, plasma density and temperature that are formed in a self-consistent configuration.Though the variations obtained via such a method are not self-consistent, they allow us to get some insights into TEMF influence on CS structure.In addition, these profiles could be regarded as a zero order approximation to the problem.
Let's consider the plasma density profile n p (z).In the absence of TEMF, this profile exhibits a maximum in the CS centre (maxn p =n m ) and possess nonzero values (n 0 ) at the boundary (Zelenyi et al., 2000).Gradient ∂n p (z) ∂z determines the value of diamagnetic drift velocity v DM ∼ ∂n p (z) ∂z (Baumjohann and Treumann, 1996), which is a part of the dispersion equation for CS drift modes: Reω=kv DM (Daughton, 1999;Buechner and Kuska, 1999;Zelenyi et al., 2009).As the turbulence level increases, Fig. 3 displays that the ratio n m n 0 is lessening while the gradient ∂n p (z) ∂z decreases with the increase of turbulence level δ.That means, the n m n 0 value obtained in the self-consistent modeling as a result of TEMF development should be less than the value obtained in the absence of turbulence.Also waves forming TEMF (v φ = ω k ≈ v DM ) will slow down as the turbulence magnitude increases since the density profile n(z) is flattering.

Comparison of acceleration mechanisms
One can obtain particle energy distribution in the outer CS region (for |z|>L) after the interaction of the particles with TEMF.Energy distributions for various values of δ parameter and wave phase velocity v φ are shown in Fig. 4. As one can see in the figure, the particles are effectively heated in the considered system.After particle-turbulence interaction particles with energy larger than ε 0 =v 2 0 2 by a factor of hundreds occur in the system.Unfortunately, we do not obtain power law tail for intervals larger than a decade in energy (as well as in most paper dealing with particle acceleration, see Pritchett, 2006;Drake et al., 2006;Perri et al., 2009).However, as the figure shows, the obtained distribution can have the power law tail: starting from energy 500ε 0 we get f ∼ε −κ−1 and κ<5.One can also see that the turbulence level δ influences particle heating more than the value of v φ .Moreover, observations in magnetosphere plasma sheet performed by means of various spacecraft suggest κ∼3-5 (Vasyliunas, 1968;Sarris et al., 1976;Christon et al., 1989), which is in agreement with the estimate obtained above.
The possible formation of non-maxwellian distribution due to particle interaction with E y field in the SM has already been demonstrated (Taktakishvili et al., 2003).We will not discuss it further in detail but compare the two acceleration mechanisms.Let us introduce the dimensionless ratio Ēy = cE y v 0 B 0 .In case of B 0 =20nT and v 0 ≈500km s we obtain the following value of Ēy : 0.2.Thus, Fig. 5 presents the values of average energy ε gained by particles in case of δ=0.3 as a function of the parameters v φ and Ēy .
Apparently, the dependence ε ∼E 2 y (see Fig. 5) corresponds to the nonadiabatic motion regime (Speiser, 1967).At the same time, the energy dependence from v φ is approximately linear, i.e. ε ∼v φ .This result cannot be explained in terms of single particle-wave interaction (the so called surfatron acceleration mechanism where one has ε ∼v 2 φ , Neishtadt et al., 2009).On the other hand, here we numerically consider particle interaction with an ensemble of waves propagating at different angles with respect to each other.The analytical investigation of the obtained result should be carried out in a future.
Figure 6 thus suggests that TEMF increases the average particle energy much more effectively than the dawndusk electric field with realistic amplitudes ( Ēy =0.2 or E y <2µV m).It appears that the combination of these two mechanisms may lead to particle energization by a factor of 20-30.

The influence of TEMF anisotropic spectrum on particle acceleration
All the results presented above are obtained for wave ensembles with identical phase velocity v φ spreading uniformly in (x,y) plane.On the other hand, in the Earth's magnetosphere the majority of waves mainly propagate in y directions (dawn-dusk) (Runov et al., 2005;Zelenyi et al., 2009).To take this effect into account, a modified dispersion equation has been used.On the plane (k x ,k y ) we introduce the angle α which is measured from the axis k x .We assume that the phase velocity is equal to zero in the following sectors: 2π−α<ϕ<α, π−α<ϕ<π+α (ϕ=arctan k y k x ), while it has the same value v φ in the intervals:α<ϕ<π−α and π+α<ϕ<2π−α (Fig. 6b).According to available experimental data, the values of α is of the order of 20-40 • (Zelenyi et al., 2009).As the value of α increases, a growing number of waves become stationary (ω k =0) and the particle acceleration becomes accordingly less effective.Taking into account the shorter residence time of the particles in the central CS region, the increasing α value also reduces the efficiency of this ion heating mechanism.The modeling performed here indicates that the maximum value of particle energy decreases with the growth of α and the distribution converges towards the maxwellian distribution (Fig. 6a).However, even for α=45 • the value of maximum particle energy is ∼300ε 0 (thus, the considered mechanism of acceleration is still effective).

Ion residence time in the field reversal region
The ion residence time in the field reversal region is an important characteristic of particle dynamics.The thickness of this region (L c ) can be estimated by the well-known method: to assume that the averaged Larmor radii in the inner region is equal to its size v T mc eB 0 tanh L c L = L c (here v T is the particle thermal velocity).This equation demonstrates that in the region |z|<L c = Lv T mc eB 0 particles are not affected by magnetic field B x (z) in the lowest order approximation, their Larmor radius being larger than the L c (Dobrowolny, 1968).Also, the time during which a Speiser particle travels in the region |z|<L c in the absence of turbulence is about the half period of Larmor oscillation in the field B z : T =2π mc eB z =2π ω n (Speiser, 1965(Speiser, , 1967)).
For the trapped particles the residence time value is much larger (Chen and Palmadesso, 1986;Buechner and Zelenyi, 1989;Ashour-Abdalla et al., 1991) and therefore, ensemble average time is T ω n >1 (here we use normalized time T =T 2π ).
As they travel in a central CS region, particles gain energy very efficiently because of the absence of B x field (unmagnetized particles can be accelerated by turbulent electric field while magnetized ones only drift in crossed magnetic and electric fields).Thus, Fig. 7 illustrates the correlation between the turbulence level δ and particle residence time.As one can see in this figure in the absence of turbulence (δ=0), all particles remain in the central CS region during a time interval T ω n ∼2.In the presence of turbulence (δ<0.5), this time interval increases, yielding T ω n >3.Such an impact of turbulence has been first found in SM (Veltri et al., 1998).With the increase of turbulence level δ, particles gain more and more energy; hence, velocities become large.As a result, for δ>0.5, no particle remain for a long time in the central CS region (in contrast to SM where with the increase of δ average particle residence time grows).

Correlation between particle energization and residence time in the CS
As shown above, there is a dependence between the turbulence level δ and the average particle residence time inside CS T ω n .There is also a similar correlation between   turbulence level δ and the average value of energy gained by particles ε .This is more apparently indicated in Fig. 8.One can see in this figure that with the growth of the δ and v φ parameters the value of T ω n decreases while the ε value increases correspondingly.In case of stationary (frozen) turbulence (v φ ∼0) the value of ε remains constant while T ω n increases with δ becoming larger.On the contrary, in case of fast turbulence v φ ∼v 0 , energy ε increases while the time T ω n decreases as δ increases.
Except for several cases of the ε and T ω n dependence on δ and v φ , the correlation between residence time and average energy can be expressed by the following heuristic for- mula: Here C 1 and C 2 are constants and ε 0 =v 2 0 2 is the initial energy.
In accordance with this formula, for any values of δ and v φ the average particle residence time in the CS centre depends on the average energy gained by particles.In order to prove this assumption the least-squares method has been used (Fig. 9).The constants obtained are such that the minimum possible value of T ω n approximately equals to the time required for one turn in B z (C 1 ≈ω −1 n ).Typical energy is ∼14ε 0 .Still, note that the equation ε T ω n −C 1 =ε 0 C 2 cannot be applied to the case of small turbulence level (δ 1).
Equation ( 6) may be explained in the following way: with the growth of the average energy value ( ε ε 0 1), the particle residence time converges toward some constant C 1 ω n that corresponds to a situation in which all the particles become transient and spend approximately one period of rotation about B z within the CS.On the other hand, provided that particle energy values are small ( ε ε 0 ∼1) the majority of particles are trapped in the CS and stay there for more than one period of quasiperiodic motion about B z ( T ω n C 1 ).

Velocity distribution of ions leaving the CS
Finally, we examine the structure of ion velocity distributions obtained.A typical distribution at the CS boundaries is a lima bean shape distribution with nonzero flow along v (Lyons and Speiser, 1982;Ashour-Abdalla et al., 1991).These distributions directly follow from non-adiabatic ion dynamics and acceleration along E y during Speiser motion (Lyons and Speiser, 1982;Ashour-Abdalla et al., 1991).Figure 10 reveals that the results obtained with the present model are in agreement with in-situ observations of f v ,v ⊥ (Raj et al., 2002;Ball et al., 2005).This figure shows the results of calculation both for the case when E y is absent (only TEMF effect was included) and when E y field is included to simulations.The figures correspond to the case v ⊥ >0 (for v ⊥ <0, the distributions are mirror symmetric).In-situ measurements also reveal that, in the vicinity of active CS regions, the ion distributions may display larger energies as well as larger velocity dispersion (particle temperature) than the particles could get due to electrostatic effect even if all cross-tail potential drop would be converted to particle energy (Grigorenko et al., 2009).According to our model, the observed effects result from particle acceleration by inductive electric fields in the DM case (Fig. 10e).

Conclusions
In this paper we carried out a numerical study of particle acceleration and transport in turbulent current sheets in three dimensions.The average field is modeled as a neutral sheet with a normal magnetic field component B z .The E y electric field due to the solar wind-magnetosphere coupling is included.The turbulent electromagnetic field is modeled as a spectrum of traveling waves with phase speed v φ .Several model parameters are varied in order to investigate the acceleration process in turbulent current sheet.The main results are presented below.
Particles energization has been obtained by considering the model of CS with TEMF formed by electromagnetic wave ensembles.In addition, a tentative formation of distributions with power law energy tail f ∼ε −5 has been ob-served.The proposed mechanism can compete with other mechanisms of particle acceleration near X-line (Hoshino, 2005;Pritchett, 2006;Drake et al., 2006).
The particle energization rate and residence time obtained allow us to derive a correlation between average energy and average residence time for the case with δ>0.2: ε ≈14ε 0 T ω n −1 −1 .Thus, typical value of residence time is T ω n ∼1 and energization rate may reach a factor of 14.
Particle acceleration under the combined effect of electrostatic (E y ) and electromagnetic fields (DM) results in the formation of velocity distributions similar to the ones observed at the edges of CS (Lyons and Speiser, 1982;Ashour-Abdalla et al., 1991).Moreover, the presence of TEMF allows to explain additional particle accelerations (different from acceleration due to E y ) often observed in velocity distribution on open field lines (Grigorenko et al., 2009).

Fig. 1 .
Fig. 1.Profiles of background magnetic field components with the presence of TEMF (black line) and without it (dotted line) and inductive electric field.

Fig. 2 .
Fig. 2. Trajectories and energy values of particles and magnetic field components along the trajectories for CS (a); for SM with δ=0.5 (b); for DM with δ=0.5 and v φ =v 0 (c).

Fig. 5 .
Fig. 5. Average particle energy as a function of phase velocity and electric field.

Fig. 7 .
Fig. 7. Particle residence time distribution for various turbulence levels in the CS centre.

Fig. 8 .
Fig. 8. Maps of the values of average energy and residence time of particles in the CS centre.

Fig. 9 .
Fig. 9.Average residence time T ω n as a function of average energy obtained by modelling (grey points) and curve ε T ω n −C 1 =ε 0 C 2 with C 1 , C 2 determined by means of the least-squares method (black curve).

Fig. 10 .
Fig. 10.Particles velocity distribution in log scale (log 10 N N , N is a number of particles with corresponding velocities v , v ⊥ and N -general number of particles).All the figures have mirror symmetric for the change-over: v ⊥ →−v ⊥ .