Ion Diffusion and Acceleration in Plasma Turbulence
Francesco Pecora Sergio Servidio Antonella Greco William H. Matthaeus David Burgess Christopher T. Haynes Vincenzo Carbone Pierluigi Veltri
LaTeX2e


























[subfigure]labelformat=empty


    
  
 
      
      
 




[-5pt]0pt20pt






































Fig. <ref>
Eq. <ref>
Eq.s <ref>


















[name=Luca, color=red]LSV
[name=Denise, color=blue]DP
[name=Oreste, color=Green]OP
[name=Sergio, color=Maroon]SS
[name=Franco, color=Gray]FV


sistema
{[ .






F. Pecora et al.]F. Pecora^1Email address for correspondence: 
francesco.pecora11@unical.it, S. Servidio^1, A. Greco^1, W. H. Matthaeus^2, D. Burgess^3, C. T. Haynes^3, V. Carbone^1  and P. Veltri^1

^1Dipartimento di Fisica, Università della Calabria, I-87036 Cosenza, Italy

^2Bartol Research Institute and Department of Physics and Astronomy, University of Delaware, Newark, DE 19716, USA

^3School of Physics and Astronomy, Queen Mary, University of London, 327 Mile End Road, London, E1 4NS 

 






Ion Diffusion and Acceleration in Plasma Turbulence
    [
    ?; revised ?; accepted ?. - To be entered by editorial office
=================================================================


Particle transport, acceleration and energisation are phenomena of major importance for both space and laboratory plasmas. Despite years of study, an accurate theoretical description of these effects is still lacking. Validating models with self-consistent, kinetic simulations represents today a new challenge for the description of weakly-collisional, turbulent plasmas. We perform two-dimensional (2D) hybrid-PIC simulations of steady-state turbulence to study the processes of diffusion and acceleration. The chosen plasma parameters allow to span different systems, going from the solar corona to the solar wind, from the Earth's magnetosheath to confinement devices. To describe the ion diffusion, we adapted the Nonlinear Guiding Center (NLGC) theory to the 2D case. Finally, we investigated the local influence of coherent structures on particle energisation and acceleration: current sheets play an important role if the ions Larmor radii are on the order of the current sheets size. This resonance-like process leads to the violation of the magnetic moment conservation, eventually enhancing the velocity-space diffusion. 



Authors should not enter PACS codes directly on the manuscript, as these must be chosen during the online submission process and will then be added during the typesetting process (see http://www.aip.org/pacs/ for the full list of PACS codes)






§ INTRODUCTION


Processes such as turbulence, diffusion and particle acceleration are ubiquitous both in astrophysical and laboratory plasmas. Understanding particle diffusion is of fundamental importance in order to characterize the distribution of the charged gas in the universe and, more specifically, in the heliosphere. The understanding of the energetic particle motion, originating, for example, from solar flares or coronal mass ejections, can help to prevent injuries for space travelers, as well as hardware damages for satellites. Moreover the distribution of heavy ions in the Earth's magnetosphere can have effects on climate changes <cit.>. Particle transport theory is fundamental for the dynamics of the solar corona <cit.> and of the interplanetary medium <cit.>. In laboratory plasma experiments, the magnetic confinement could be improved by understanding what affects particle turbulent transport <cit.>.


Charged particle dynamics depends on the stochastic motion of the magnetic field lines. The random walk of the magnetic field lines affects the diffusion both across and along the mean magnetic field <cit.>. Charged particles gyrate along the magnetic field but, if the field is turbulent, they spread in the perpendicular direction <cit.>, jumping from a field line to another. The turbulent nature of space and laboratory plasmas suggests that the best way to describe their dynamics is given by the statistical approach <cit.>. In this statistical (Lagrangian) description, particles move in a turbulent electromagnetic field, similarly to the motion of neutral tracers in atmospheric turbulence <cit.>. The turbulent nature of the fields that scatter the particles limits the analytical treatment of the subject. For this, we still lack an exact and universal theory for describing the perpendicular diffusion <cit.>. 

A common theoretical approach, applied when studying diffusion of charged particles, relies on calculating separately the diffusion coefficient in directions parallel and perpendicular to the main guiding field <cit.>. Currently the NLGC theory <cit.> gives a rather  accurate prediction of the diffusion coefficient for systems with a three-dimensional (3D) geometry. However, this theory has been tested only with test-particle simulations - simulations that do not take into account the back-reaction of the particles motion on the electromagnetic fields. Nowadays, more realistic simulations of turbulence are available and is therefore interesting to test the validity of these predictions, using, for example, kinetic models.

A full 3D description of plasma turbulence requires huge computational efforts that can be lightened by reducing the dimensionality of the problem. Magnetohydrodynamics (MHD) turbulence simulations <cit.> have shown that, if a strong guiding magnetic field is present, turbulence loses the property of isotropy. The magnetic structures are mainly present in the plane perpendicular to the main field <cit.>. These structures are essentially composed by a sea of flux tubes (magnetic islands) and current sheets, where reconnection might eventually occur <cit.>. In this complex pattern, particles diffuse and accelerate. 

MHD simulations suggest that turbulence can be considered as a network of reconnecting magnetic islands, that merge, change topology and convert energy <cit.>. Magnetic reconnection is considered one of the most effective mechanisms for particle acceleration and energisation <cit.>, being crucial for explosive events in the solar atmosphere, like solar flares <cit.> and coronal mass ejections <cit.>. Analogously, in the Earth's magnetosphere <cit.> and far away to the heliopause <cit.>, magnetic reconnection is thought to be a very active mechanism. Observations and measurements strongly relate the magnetic island merging <cit.> and the magnetic field discontinuities <cit.> with the increase of suprathermal energetic particles. This increasing number of energetic particles in zones of merging magnetic islands and of strong discontinuities can be the key interpretation to relate the magnetic reconnection with acceleration and energisation <cit.>.

Both 2D and 3D simulations of MHD turbulence with a strong guiding magnetic field <cit.> reveal power law energy spectra for test-particles. The power-law spectra can be indeed observed for anomalous cosmic rays <cit.>, energetic ions in the solar wind <cit.> and energetic electrons in solar flares <cit.>.  It is worth noting, however, that since test-particle models are not self-consistent, particles in these simulations can reach extremely high energies. These very high energies might not be commonly observed in the solar wind, in usual turbulence conditions, especially for ions, suggesting that a self-consistent kinetic model need to be used in order to correctly describe the plasma dynamics.

The present work is organized as follows. In Section <ref> we will present a global overview on the model that describes the plasmas in the low-collisionality limit. We will introduce the numerical algorithm that simulate the plasma-particle dynamics. In Section <ref>, the particle trajectory statistics will be presented, interpreting the numerical results via the plasma turbulent diffusion theories. The investigation of the acceleration process will be the central topic of Section <ref>, where we will perform local analysis of the acceleration phenomenon.  Finally, in the last Section, we will present the discussions, conclusions and future perspectives.









§ SIMULATIONS OF 2D PLASMA TURBULENCE


We perform three hybrid-PIC simulations, varying the plasma β (the ratio between kinetic and magnetic pressure), with β=5, 0.5 and 0.1, in order to cover a wide range of relevant plasma scenarios. These chosen values are such that the high value is typical of plasmas found in the Earth's magnetosheath; the value of β=0.5 is close to the typical solar wind conditions, whereas the lowest is appropriate for the the solar corona and laboratory devices.

Particles and fields are described by the Vlasov-Maxwell system eq:hybridPIC that we solve via a hybrid-PIC approach, using kinetic ions and fluid electrons. The system of equations is given by


    ẋ =  v
       v̇ =  E+ v× B
       B t = - E
       E = - ( u×B ) + 1n j×B - 1/n∇ P_e +ηj .


In the above equations, x is the particles position, v their velocity, E is the electric field, B is the magnetic field, u is the proton bulk velocity (the first moment of the velocity distribution function), n is the proton number density (the zero-th moment of the ion velocity distribution function), j is the current density. The pressure term is adiabatic P_e = β n^γ and η=0.006 is the resistivity that introduces a small scale dissipation, for numerical stability. The electric field is given by the generalised Ohm's law. In the simulations distances are normalized to c/ω_p_i, where c is the speed of light and ω_p_i is ion plasma frequency. The time is normalized to Ω_c_i^-1, that is the ion cyclotron frequency. Finally, velocities are normalized to the Alfvén speed v__A= c Ω_c_i/ω_p_i.

The three simulations have the same initial conditions: uniform density and a Maxwellian distribution of particles velocities with uniform temperature. We impose large scale fluctuations in order to mimic the motion of the large energy-containing vortices. We have chosen random fluctuations, at large scale, for both magnetic field and the ion bulk velocity field. The eq:hybridPIC are solved on a square grid of size L_0=128 d_p, where d_p is the proton skin depth defined as d_p = c /Ω_c_i, discretized with 512^2 points, with periodic boundary conditions. This initial state consists of a 2D spectrum of fluctuations, perpendicular to the main field B_0 (the latter chosen along z). The fluctuations amplitude is δ b/B_0 ∼ 0.3. To suppress the statistical noise of the PIC method we use 1500 particles per cell (about 4× 10^8 total particles).


To describe the magnetic topology, in the 2.5D approximation, it is useful to define the in-plane magnetic field as B_⊥ = ∇ a_z ×ẑ, where a_z is the magnetic potential and ẑ is the out-of-plane (axial) unit vector. The current in the axial direction j_z = (B_⊥ ) ·ẑ = -∇^2 a_z. To achieve a stationary state of fully developed turbulence, we initially let the system decay freely, and then we introduce a forcing at the time t^* at which nonlinearity reaches its peak (roughly the peak of ⟨ j_z^2⟩), namely t^* ∼ 25 Ω_cp^-1. The forcing consists of “freezing” the amplitude of the large-scale modes of the in-plane magnetic field, with 1≤ m≤ 4, with constant phases. This corresponds to a large-scale input of energy, as described in <cit.>. 


In order to have a significant statistics, we perform our analysis when a steady state has been achieved, namely for 50<tΩ_cp<250. fig:jz shows the shaded contour of the current density j_z along with the contour line of the vector potential a_z, as the system evolves toward turbulence. Panel a shows the initial state of the system, where islands and current sheets are not defined yet. As time goes on, smaller vortices and sharp current sheets develop. In particular, these magnetic structures represent magnetic islands (flux tubes in 3D). The global appearance of the system remains unchanged, in a statistical sense, when the peak of nonlinearity has been reached. The regions of big magnetic gradient appear in between reconnecting magnetic islands, with associated intense current sheets <cit.>. 




To better describe the state of fully developed turbulence, we computed the Fourier spectra, as a function of the wavenumber |k|, of both the electric and magnetic fields. The power spectra for the trace of the correlation tensor, |B̃(k, t)|^2, where B̃(k, t) are the respective Fourier coefficients, are reported in fig:EB_dev_spectra. These spectra exhibit the classical scenario of MHD turbulence, in which energy flows from large to small scales. This cascade of energy occurs over the so-called inertial subrange, where energy that scales as k^-5/3 <cit.>. In particular, fig:EB_dev_spectra shows the turbulent development. The state of fully developed turbulence is achieved after 50Ω_cp^-1. After this time the spectra are almost stationary (do not experience large fluctuations). It is important to note, that as observed in the solar wind <cit.>, the power in the electric fields is higher, at k's that correspond to characteristic ion lengths. This is in agreement with previous studies and simulations of plasma turbulence <cit.>.



In order to understand the dynamics of particles in a turbulent scenario, where structures are present at different scales, it is useful to measure different characteristic lengths. The correlation (or integral) length of a turbulent field generally corresponds approximately to the size of the largest energy-containing turbulent eddies, (in our case the typical sizes of the islands.) This length can be obtained from the two-point correlation function C(r) as

    λ__C = ∫_Ω C(r) dr =   1⟨ b^2 ⟩∫_Ω dr⟨b(x+r)·b(x) ⟩.

Here b represents the magnetic fluctuations 

and ⟨·⟩ is the average over the total volume Ω. Note that fluctuations are isotropic in the plane perpendicular to the main field, hence x and r are in-plane vectors and the above length is essentially the same, namely C(r_x)∼ C(r_y)≡ C(r). Moreover, λ_C is very similar for the bulk velocity fluctuations, because of the choice of the initial conditions and the driving. In the above definitions we suppressed the time dependence because  of stationary, and we computed these lengths averaging over time, in the steady state regime. For our system λ_C ∼ 10 d_p and does not vary for runs wit different plasma β.

Apart of the above energy-containing scale, which is a large scale characteristic length of turbulence, it is important to characterize also the smallest scales properties. In hydrodynamics the Taylor length is the scale at which the viscous dissipation term is no longer negligible <cit.>. This scale is related to the largest width of the structures,  where dissipation starts to be relevant. In our case, the magnetic Taylor scale  can be defined as

    λ__T = √(b_⊥^2⟨ j_z^2 ⟩).

In the above expression, b_⊥ is the root mean square of the magnetic field in the plane perpendicular to the axial direction (perpendicular to the main guiding field), and ⟨ j_z^2 ⟩ is the axial current averaged over the 2D space, and averaged also over the times in which the system is in the stationary state of fully developed turbulence. This scale can be considered as the scale at which the inertial range of turbulence terminates <cit.>. In our simulations λ__T∼ 1.7 d_p. Since the Taylor length λ__T gives the biggest size of the dissipative structures, we introduce now an average measure of these structures, namely the typical width of the layer cores (O-points). The current sheet width, δ_c,  is the half maximum width of the current sheet intensity profile. The average core width of the current sheets we measured as δ_c ∼ 0.3 d_p.










§ PARTICLE DIFFUSION IN PLASMA TURBULENCE






We investigate the particle motion in steady state turbulence. Among all PIC macro-particles, we have followed the path of 10^5 samples, verifying the convergence of the statistical results. From their  positions as a function of time, it is possible to compute the diffusion coefficient, measuring their mean squared displacement. After a general overview on the erratic trajectory of particles in 2D turbulence, we will see whether the plasma β affects the statistics of ions diffusion and acceleration.


In fig:2dtraj we show the 2D trajectories of some randomly selected particles, during the whole simulation, for two distinct values of β (low and high β). As it can be seen, particles that move in high plasma are less magnetized and their Larmor radii is large. Particles in low-plasma, instead, are highly magnetized and perform tight gyrations, with more abrupt change of direction. Small-β particles spread less than in the high β case. Generally, in all cases, one can notice that some particles have closed orbits as they are trapped in a magnetic island while others experience turnovers that suddenly bend the trajectory. As we shall see, particles that perform closed trajectories are mostly trapped in turbulent vortices, whereas sharply segmented trajectories can belong to particles which encounter local intense turbulent structures, as magnetic discontinuities or current sheets, that make them deviate suddenly <cit.>.




The overall motion, as represented in fig:2dtraj, seems quite random: particles scatter and diffuse in time. We computed the square displacement Δ s^2=Δ x^2+Δ y^2, finding that, for sufficiently long time intervals, the motion is diffusive, namely

    ⟨Δ s ^2 ⟩ = 2D τ.

In the above expression, D is the diffusion coefficient, τ is the time interval over which the particle moves by Δ s, and the brackets now represent an ensemble average computed over particles.
Note that τ=t-t_0, and because of stationarity we varied t_0 in the range of the steady state regime. To understand the dependence of the diffusion coefficient on particle energy, we divide the particles in energy bins and evaluate the above diffusion coefficient, collecting particles depending their average energy. The averaged mean squared displacement as a function of the time interval for several energy values is shown in fig:brown_2dt. It is possible to observe the normal diffusive behavior is achieved for for τ≳ 60 Ω_cp^-1. By fitting the curves after this interval with the eq:brown_disp, we can obtain the numerical value of D.
For shorter time intervals particles seems to not follow the linear law. Indeed, for short time intervals the particle motion cannot be stochastic since they have “memory” of their initial condition, and temporal correlations exist.  It is clear that the higher is the particles energy, the larger is the diffusion coefficient. Before reporting the numerical results on the diffusion of particles in self-consistent kinetic simulations, we will introduce the theoretical background.



 §.§ A Non Linear Guiding Center (NLGC) model in 2D

We want to develop a theory that can describe particle diffusion in 2D turbulence, and test it with the self-consistent plasma simulations. Currently, a useful description of perpendicular diffusion is given by the Non Linear Guiding Center (NLGC) theory <cit.>. This theory works well for 3D turbulence, with a uniform mean magnetic field, and it has been tested with test-particle simulations. To develop the 2D version of the theory we start from the 3D NLGC derivation, considering the particle gyromotion. Using the Taylor-Green-Kubo (TGK) formulation <cit.>, one can write the diffusion coefficient D_xx, for instance along the x direction, as

    D_xx = 1B_0^2∫_0^∞ dτ⟨ v_z(0)B_x(x(0),0) v_z(τ)B_x(x(τ),τ) ⟩,

where it is assumed that the particle motion projected on the 2D plane follows the 2D  projection of the magnetic field lines.
We assume that the magnetic field fluctuations in the perpendicular plane are completely uncorrelated with the velocity in the z direction and this allows us to write the diffusion coefficient as

    D_xx = 1B_0^2∫_0^∞ dτ⟨ v_z(0) v_z(τ) ⟩⟨ B_x(x(0),0) B_x(x(τ),τ) ⟩.

Now we encounter a deviation from the 3D NLGC, by saying that the velocity correlation function in the z direction is nothing but the square value of the z velocity i.e.

    ⟨ v_z(0) v_z(τ) ⟩∼  v_z^2.

This assumption has been directly proved by measuring the quantity ⟨ v_z(0) v_z(τ) ⟩ in our simulations. fig:cvz shows that the motion along z is given by free streaming, with ⟨Δ z^2 ⟩∼ t^2 and therefore the velocity v_z, in the direction along the main field, is almost constant.



Using this approximation, we can estimate the diffusion coefficient neglecting the parallel scattering, and get: 

    D_xx = v_z^2B_0^2∫_0^∞ dτ ⟨ B_x(x(0),0) B_x(x(τ),τ) ⟩.

Using the Corrsin's independence hypothesis to transform the integrand function as

    ⟨ B_x(x(0),0) B_x(x(τ),τ) ⟩ = ∫ dr  R_xx(r,τ) P(r,τ)

one arrives at

    D_xx = v_z^2B_0^2∫_0^∞ dτ ∫ dr  R_xx(r,τ) P(r,τ),

where R_xx(r,τ) is the Eulerian two-point two-time correlation tensor and P(r,τ) is the probability function of the particle having displacement r after a time interval τ <cit.>. At this point is convenient to express R_xx(r,τ) using its Fourier transform

    R_xx(r,τ) = ∫ dk  S_xx(k,τ) e^ik·r.

We can model S_xx(k,τ) = S_xx(k)Γ(k,τ), where Γ(k,τ) represents the time-propagator of the spectrum, and S_xx(k) is spatial spectrum. Usually, this functional form is described via the so-called sweeping decorrelation mechanism, namely Γ(k,τ) ∼ e^-τ/τ_c(k), where τ_c(k) is a characteristic decorrelation time. The latter is usually thought to be the  the “sweeping decorrelation time” <cit.>, in analogy with fluid turbulence (see below). 

At this point, we can proceed with the theoretical estimate of the diffusion coefficient, given by

    D_xx = v_z^2B_0^2∫_0^∫ dk  S_xx(k) dτΓ(k,τ)  ∫ dr  P(r, τ) e^ik·r.

Since P is Gaussian, then

    ∫ dr  P(r, τ) e^ik·r = e^-(k_x^2D_xx + k_y^2D_yy)τ.

In In axisymmetric turbulence D_xx = D_yy≡ D, S_xx = S_yy≡ S and, since we have only in-plane structures (i.e. k_z=0) k_x^2 + k_y^2 ≡ k^2, the diffusion coefficient is

    D = v_z^2B_0^2∫ dk  S(k) ∫_0^ dτ  e^-k^2Dτ e^-τ/τ_c( k).



Finally, integrating over τ, one gets

    D = v_z^2B_0^2∫ dk S(k)[τ_C(k)]^-1 + k^2D.

Note that this prediction is only a small modification to the NLGC theory, having suppressed the z-dependence by using the particle free streaming along z (eq:vz2appr). 

In order to obtain a first estimate of the diffusion coefficient D, we can make an approximation.  Assuming that the sweeping decorrelation time is rather long, (as in the case of large scale slow driving, for example) we can drop the term τ_C(k)^-1 in eq:2DNLGCD. This latter term, indeed, makes the equation more complex to solve. We therefore obtained:

    D^* ∼√(v_z^2B_0^2∫ dkS(k)k^2).

The reader may notice that this result is essentially of the same form as the so-called Field Line Random Walk (FLRW)  limit of perpendicular particle scattering <cit.>. For such cases, D= v_eff D_fl where D_fl is the Fokker Planck coefficient for field line transport, and v_eff is the effective velocity of the particle along the magnetic field <cit.>.

To evaluate the diffusion coefficient via eq:apprD and eq:2DNLGCD, we computed the average power spectrum S(k), from the simulations. fig:D_v22 shows the numerical values of D obtained by fitting the particles mean squared displacements in fig:brown_2dt. The theoretical values of the 2D NLGC, evaluated via both the exact and the approximated formulas, are reported as a function of the particles energies. The simulations results follow fairly well the theoretical prediction at low . The theory slightly deviates at very high energy (β's), although the functional monotonic behaviour is very similar. We will note that we do not observe large energy excursions like in test-particle models, due to the self-consistency of the ions treatment.








Here we briefly investigate the decorrelation mechanism in self-consistent, plasma turbulence. We performed a Fourier transform in time of the magnetic fluctuations, computing the propagator Γ(k, τ), as described in Ref.s <cit.>. As in fluid, MHD, and Hall MHD models of turbulence, this time-dependent correlation of turbulence strongly depends on the amplitude of k, as reported in fig:rab. As it can be seen, the decorrelation mechanism depends on k and drops quickly in time (only a few inertial range modes are reported). From this functional form, we computed the decorrelation time τ_C(k), represented in the panel b of the same figure. The decorrelation time scales as ∼ 1/k, indicating the clear dominance of the sweeping effect. To be more quantitative, in order to compute the diffusion coefficient for our experiment, we found that τ_C(k) ∼ 3/(δ b_⊥ k d_p), where δ b_⊥ is the rms of the in-plane magnetic fluctuations.









§ PARTICLE ACCELERATION


In agreement with previous literature <cit.>, we have found that diffusivity depends on particles energy. The behavior of diffusion is in agreement with a modified model of the NLGC theory, here adapted to collisionless 2D plasmas. Moreover, we observed that many of these particles experience variations of their momentum – namely an acceleration process. In this section, we investigate the mechanisms responsible for acceleration, using conditional statistics. We will establish the relationship between the effectiveness of the acceleration mechanism and the particle energies, and the existence of possible resonance conditions that energize the ions.

We have computed, for each particle, the Lagrangian acceleration a= ∂v/∂ t,  where v is the particle velocity. We computed the the Probability Distribution Function (PDF) of the Lagrangian acceleration, for all the particles, at different times and for different simulations. The acceleration has been computed by using a 6^th order finite difference (centered) method. The distributions are reported in fig:PDFa2, for the two extreme β, at the beginning and at the end of the simulation. As it can be observed, the PDFs of particle acceleration in high plasma are well described by the χ^2 distribution. Let { x_1, x_2, …, x_k } be a set of k independent normally distributed variables. The sum of the square of these k variables distributes according to the χ^2_k PDF that is defined as

    PDF(χ^2_k) = 12^k Γ(k/2) x^k/2-1 e^-x/2

where k are commonly referred to as degrees of freedom, and Γ(x) is the gamma function. In our case the three independent variables are the acceleration components { a_x, a_y, a_z } and the square modulus of the acceleration |a|^2 = a_x^2+a_y^2+a_z^2 distributes according to the χ^2_3 with 3 degrees of freedom.

For β=5,  in particular, this distribution does not change in time, indicating the lack of very extreme events during the evolution of the system. Whereas at low β (β = 0.1), the acceleration, which is initially randomly distributed, develops a tail at later times (blue solid line). This means that the number of particles with anomalous acceleration is higher for low plasma, suggesting that a physical process that depends on the plasma conditions is at work.




A straightforward acceleration mechanism can be due to an electric field parallel to the local magnetic field: 

    E_∥ = E·B|B|.

We identified the anomalous particles, i.e. the particles with acceleration values exceeding the variance of the followed distribution and represented their positions, at a given time, over the parallel electric field map with a contour plot of the magnetic potential in the z direction, in order to see whether they show a correlation with the turbulent structures. This map is reported in fig:aQmap. The figure shows that particles with anomalous acceleration are affected by the above field: accelerating ions are non uniformly distributed, they cluster where the parallel electric field is more intense, on the flanks of the magnetic islands.





In order to establish a more quantitative link between the parallel electric field and possible local acceleration effects, we compared the PDFs of the single terms of the electric field that can be parallel to the local magnetic field (namely ∇ P_e and η j) and of the total parallel electric field itself (E_||) 

over the positions of non-anomalous particles, and the PDFs of the same quantities measured only at the positions of high-acceleration particles. The particles were distinguished by thresholds with respect to the global acceleration PDF: the anomalous particles are those which have an acceleration value exceeding the 3σ of the global distribution, whereas normal particles have an acceleration value within 1σ. Note that the above conditional statistics would give the same distribution only if acceleration and E_|| are uncorrelated. These PDFs are shown in fig:QPDF, where the quantities with superscript ^↑ are measured at the positions of anomalously accelerated ions, whereas the fields with superscript ^↓ are related to normal particles. As it can be seen, accelerated particles have higher parallel electric field –  the population of the distribution is higher at big E_|| regions. Moreover we used the Partial Variance of Increments (PVI) method to find candidate regions likely to be identified as coherent structures. These structures contribute to non-Gaussian statistics and therefore to intermittency. The PVI time series is substantially defined as the normalized sequence of magnetic increments. By applying the PVI technique, we found that the regions of bigger E_|| occur in correspondence of magnetic discontinuities and not in smooth regions, as the PDF of E_|| conditioned on PVI values <cit.> clearly evidences in the inset of panel a of fig:QPDF.


This statistics further confirm the relation between the parallel electric field and the stochastic acceleration mechanisms in 2D turbulence. 










This acceleration mechanism in the out-of-plane direction has a global effect of elongating the ion velocity distribution function (VDF). In order to see if this typical alignment effect is present in our numerical experiments, we computed both the PDF of the angle between particles velocity and the local magnetic field and the PDF of the angle the particles velocity has with the main magnetic field (along z)

    cos(θ) =  v· B/| v|| B|,     cos(ψ) = v_z/|v|

In fig:cosine we report the distributions of both the angles, at the initial and final times of the simulation, for different values of β. At the initial time, when turbulence is very young, the distributions are quite flat, meaning particles are moving isotropically. As the simulation goes on, particles tend to align with the main magnetic field in the z direction. This effects is much more evident for the low plasma, the more magnetized one. Another feature one can notice in the low plasma is that the particles orient themselves more on the local magnetic field rather than its z component, although the z component is its main one and the difference is not statistically relevant.






It is crucial to see now whether this acceleration mechanism, locally related to intense parallel electric fields, can actually increase the particles energy, and up to which values. This correspondence is not trivial since the particles with anomalous acceleration are only a small fraction of the plasma. We have constructed the energy PDF to see whether the energy has a different behavior at different β. The energy PDFs for the low and high values are show in fig:PDFE. In the low scenario, the PDF develops a power law tail, already seen in observation and previous numerical simulations. Whereas in the high plasma the energy remains similar to the initial distribution, namely close to Maxwellian distribution. This suggests that the acceleration mechanism can energize particles, and that the process depends on the plasma β.









 §.§ Magnetic Trapping





We have seen that the acceleration mechanism involves the parallel electric field. However, there must be something else that makes this field energizing particles in low plasmas. In figure fig:aQmap we have seen that the anomalous particles lie within (or in between) magnetic islands, in agreement with previous discussions. The acceleration mechanism, evidently occurs in association with magnetic reconnection. Particles temporarily trapped inside small flux tubes or in the dynamically active region of larger flux tubes experience Fermi-like processes. In the  region near current sheets, particles can experience nearly continuous first order energisation <cit.>. In particular we followed the trajectory of one of the most energetic particles and monitored its energy, its acceleration, the energy derivative with respect to time and the parallel current density it samples during its journey (the parallel pressure term is not shown because it is too noisy for non-statistical treatment) (fig:wonder). The figure also shows that the particle remains trapped in a magnetic island and its energy grows until it reaches a high enough energy to escape.





To estimate the escaping times we calculated the Lagrangian auto-correlation time using eq:corr_time. This quantity represents the time at which each particle experience a correlated field, and can be estimated as

    T_esc = 1⟨ v_x(t_0)^2 ⟩_a∫_0^∞⟨ v_x(t_0) v_x(t_0+τ) ⟩_a dτ.

In the above definition, the operation ⟨·⟩_a is the average only over the particles that have high acceleration. For the usual isotropy hypothesis. we computed the above relation along all directions, finding similar results. We obtained the following escaping times: T_esc^anom(β=5)∼ 0.8 Ω_cp^-1, T_esc^anom(β=0.5)∼ 3.0 Ω_cp^-1 , and T_esc^anom(β=0.1)∼ 5.7 Ω_cp^-1. In addition we computed the same quantity for all the particles, obtaining these values: T_esc^all(β=5) = 0.95Ω_cp^-1, T_esc^all(β=0.5) = 4.13Ω_cp^-1, T_esc^all(β=0.1) = 9.64Ω_cp^-1.
This scenario is consistent with the fact that particles in low plasmas develop high energy tails, since they are confined within magnetic islands for longer periods, experiencing the same parallel electric field. On the other hand, particles in high plasmas easily escape from magnetic islands and are not efficiently energized. The comparison of the trapping times values at the same β for the different kind of particle populations also supports this view. The escaping times calculated over all the particles are bigger than those of the anomalous particles suggesting that more energetic particles are more likely to escape from magnetic islands while lower energy particles remain trapped for longer period and can be coherently energised. This Fermi-like process, invoked in small scale reconnection <cit.>, is now quantitatively observed in large scale plasma turbulence. It is interesting now to further characterize the energisation process by looking at the characteristic parameters of these anomalous particles.






 §.§ The Magnetic Moment





We now further inspect the acceleration process by looking at the magnetic moment μ of ions, defined as

    μ = m v_⊥^22B.


In the above expression, m is the particle mass and v_⊥ its velocity (perpendicular to the magnetic field B, measured at the particle's position). μ is an adiabatic invariant of the system, if B is slowly varying. Indeed, the orbits are like closed circles and the flux of magnetic field passing through them is almost constant. This suggests that the magnetic moment might not be a constant of the motion in a turbulent system where several spatial scales are present, and where  magnetic field variations are neither negligible nor adiabatic <cit.>. To quantify the behaviour of the magnetic moment for each particle, we computed the normalized moment

    μ̃_p = μ_p(t)-μ_p(0)μ_p(0),

where the label p now indicates a single particle. This measure gives us information about the variation of the particle magnetic moment with respect to its initial value. fig:11murt shows μ̃_p as a function of time, for some particles, randomly selected, for two plasma values. Particles that moves in low plasma have the highest magnetic moment excursions, while in the high plasma, where particles are  not so energized in time, their magnetic moment is much more conserved. It is important to notice that the magnetic moment distribution is very similar to the distribution of energies and acceleration, further confirming the relevance of this quantity for the process of plasma acceleration in turbulence.




To see whether the magnetic moment is a constant of the motion, in a more quantitative way (fig:11murt refers only to a small portion of particles), we calculated the standard deviation σ_μ(p). This can be interpreted as the dispersion of the magnetic moment, for each particle, and is defined as 

    σ_p = √(1/T∫_t_0^t_0+T[ μ_p(t')-⟨μ_p⟩_T ]^2 dt'
    ),

where μ_p(t) is the magnetic moment of the p-th particle at the time t, and the average in the integral is calculated over the whole time. In case of a perfectly conserved magnetic moment, this quantity is null.
We have then built the PDF of

    ε = σ_p⟨μ_p ⟩_t,

that indicates how much the magnetic moment deviates from its mean value – how much the magnetic moment is “broken”. The PDF(ε) is shown in panel a of fig:epsilonmuRLPDF. As expected, the violation of the magnetic moment is much more pronounced (broader distribution) in the case of β=0.1. Low-particles have small gyro-radii and they can probably interact with the local strong inhomogeneities.




 
We can relate the magnetic moment violation directly to the turbulence characteristic scales. Our analysis suggests that particles can interact with current sheets, if their Larmor radius is much shorter than the Taylor length λ_T – the largest current sheet size (eq:lambdat ). If the radius is bigger, the particle can gyrate without even “noticing” the current sheet. The other important scale is the current width δ_c, that is the mean size of the current sheet cores. Panel b of fig:epsilonmuRLPDF shows a more quantitative view of the above speculation, suggesting a kind of spatial resonance. In this figure we report the PDF of the Larmor radius along with the turbulence characteristic lengths λ__T and δ_c.  The Larmor radii of particles in high plasma are broadly distributed and are much bigger than the current sheets thickness (while they are on the order of the Taylor scale). The case with β=0.1 shows, instead, that the average R_L are smaller than (or on the order of) δ_c. In this scenario, particles with the resonant Larmor radius feel the presence of the sharp discontinuities, and undergoes a magnetic moment break. This finally leads to the acceleration mechanism, and to the energisation of ions in turbulent plasmas.





 §.§ Approaching to Velocity Space Diffusion

We have studied the process of particles acceleration and energisation. Acceleration can be interpreted as a “motion” in the velocity space, i.e. v-diffusion <cit.>. In this last section we will  briefly mention this aspect of diffusion. As for spatial diffusion, it is instructive  to observe particles trajectories in the velocity space. In particular, we represent these trajectories in a 2D space made by v_z (parallel to the global magnetic field) and the in-plane velocity v_⊥ = √(v_x^2+v_y^2). The trajectories in the velocity space are shown in fig:vtraj.



The velocity space trajectories show, from another point of view, the same behavior we have seen in the previous sections. In the high plasma, particles are accelerated but not effectively energized by the parallel electric field. This means that they can diffuse only in pitch-angle and gyrophase, namely they move on the surface of iso-energetic shells, varying only the angle the velocity forms with the magnetic field. In the low scenario, particles, instead, are accelerated and energized by the current sheets. The combined acceleration and energisation phenomena make particles trajectories in the velocity space more complex, since particles can diffuse both in pitch-angle and momentum-amplitude. In this case, ions change their velocity modulus and can acquire a continuous range of values making the trajectories almost ergodic. The ergodicity domain is bounded by the upper limit of energy a particle can gain before escaping a vortex and run into a decelerating (de-energizing) electric field. This interesting aspect of diffusion will be further inspected in future works.




§ DISCUSSION AND CONCLUSIONS

In this work we described the diffusion and the acceleration of ions in plasma turbulence, using self-consistent kinetic simulations (kinetic ions and fluid electrons). These topics are of primary importance both in astrophysical and laboratory plasmas. We simplified the problem by using a 2D geometry, which can be a valid approximation to understand the nature of strongly anisotropic (magnetized) fluctuations <cit.>. Particles have been evolved making use of a PIC algorithm. The PIC algorithm is very useful and of practical fundamental importance when treating non-equilibrium plasmas for which the wave-particle interactions are crucial. The 2D hybrid-PIC simulations have been performed for three different β's, in order to reproduce a wide range of physical systems. The different values of β used vary for more than one order of magnitude allowing us to describe plasma scenarios spanning from the solar corona to fusion devices. The 2D maps of the current density along with the magnetic field showed the development of fundamental turbulent structures, such as vortices and current sheets.

We have then studied the motion of ions, moving self-consistently in the electromagnetic field. As suggested by previous (and numerous) test-particles studies, the motion is very erratic. Particles can be trapped in magnetic vortices or  scattered away by current sheets, wandering like a pollen in the atmosphere <cit.>, or field lines in the solar corona <cit.>. This kind of trajectories are achieved after rather long time intervals and can be statistically described within the theory of diffusion. The particle motion becomes uncorrelated when the particle is no longer trapped by the same vortex. Low energy particles have been found to have longer correlation times as they cannot easily escape from vortices.

All the existing theories on the diffusion coefficient can describe 3D systems and all these theories have been verified only with test-particle models. The NLGC theory is possibly the most precise theory that gives an estimate of the diffusion coefficient. We have then reduced this theory to the 2D case. This 2D NLGC theory has been found to be valuable in describing the diffusion coefficient of particles moving in self-consistent turbulent fields.

From the acceleration process, the PDFs show that the acceleration nature depends on the plasma β. Acceleration is a stochastic variable for high plasmas, whereas it is distributed with power-law tails for low case. We have found that the electric field component parallel to the magnetic field is correlated with particle acceleration. The “anomalously” accelerated particles, i.e. the particles with acceleration values that exceed the variance of the distribution, are connected to regions with high parallel electric field. This has been seen, qualitatively, by spotting these particle on the parallel electric field map and, quantitatively, by computing the conditional statistics.

By looking at the energy PDFs, we noticed that particles in low systems, such as in the solar wind, are more effectively energized. This process can be linked to the presence of narrow current layers. The main phenomenon acting on the particle is the breaking of the particle magnetic moment. Local spatial resonances break this constant of motion, leading to acceleration and finally to energisation of particles. In this process, particle undergo a spatial resonance with the background turbulent structures: ions that have their Larmor radius on the order of the current sheet thickness experience large excursions of their magnetic moment. These particles experience a local acceleration process, while in the case with much larger β, as can be found in the magnetospheric environment, the plasma elements do not “see” the reconnecting current sheets and the embedded parallel electric field. This kind of interaction with vortices and current sheets has been proved also by looking at characteristic scales resonances. Particles characteristic scale is the Larmor radius, whereas for turbulence we computed the Taylor length and the current width. A consistent percentage of high energy particles has a Larmor radius bigger than the Taylor length, that is, qualitatively, the in-plane length of the  biggest current sheet. This means high energy particles can barely “notice” the current sheets. 

Finally, we have introduced the concept of velocity space diffusion, because above processes might be related to the stochastic motion of particles in the velocity space. We still lack a fundamental theory to determine the diffusion coefficient in velocity space <cit.>, and we leave this work for future studies. Future works can be focused on an analytical treatment of the velocity-space diffusion. Moreover, a full 3D study will be performed in the future, taking into account also the role of kinetic electrons.




§ ACKNOWLEDGMENTS

This work is partly supported by the International Space Science Institute (ISSI) in the framework of International Team 405 entitled “Current Sheets, Turbulence, Structures and Particle Acceleration in the Heliosphere", by the US NSF AGS-1156094 (SHINE), and by NASA
grant NNX14AI63G (Heliophysics Grandchallenge Theory), the MMS mission through grant NNX14AC39G, and the Solar Probe Plus science team (ISOIS/Princeton subcontract SUB0000165).


jpp ]