D Mandal

Y Elskens
email: yves.elskens@univ-amu.fr

X Leoncini
email: xavier.leoncini@univ-amu.fr

Nicolas Lemoine
email: nicolas.lemoine@univ-lorraine.fr

F Doveil

Sticky islands in stochastic webs and anomalous chaotic cross-field particle transport by ExB electron drift instability

Keywords: Stochastic web, ExB drift instability, Hall thruster, super-diffusive transport PACS : 05.45.-a Nonlinear dynamics and chaos 52.20.Dq Particle orbits 52.25.Fi Transport properties 52.75.Di Ion and plasma propulsion 05.45.Pq Numerical simulations of chaotic systems

I. INTRODUCTION

The formation of stochastic web structures and the chaotic transport of charged particles in presence of electrostatic waves and magnetic field have been investigated for several decades [1][2][3][4][5][START_REF] Balescu | Aspects of anomalous transport in plasmas[END_REF][START_REF] Bouchara | [END_REF]. In purely chaotic situations where a central limit theorem is valid, the transport process is like a discrete time random walk, and the variance grows linearly with time [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF][START_REF] Balescu | Statistical dynamics -Matter out of equilibrium[END_REF][START_REF]Stochastic differential equations[END_REF]. But in the case of mixed phase space where both chaotic and regular trajectories coexist, the transport processes are not so clear [START_REF] Chirikov | [END_REF][12][13]. Usually, trajectories spend more time near the border of the regular region. This type of dynamics can sometimes be modeled using continuous time random walks (CTRW), where the number of jumps within a time interval [0, t), and the displacement in each jump are taken from two mutually independent probability densities, and these two probability densities fully specify the probability distribution describing the random walk [14,[START_REF] Méndez | Stochastic Foundations in Movement Ecology -Anomalous Diffusion[END_REF]. Transport in such systems can be linked with Lévy flight type processes [START_REF]Lévy flights and related topics in physics[END_REF][START_REF] Del-Castillo-Negrete | [END_REF][18]. In presence of a magnetic field, due to the interaction with electrostatic waves, the dynamics of the charged particles become chaotic and, for certain parameter values, they form stochastic webs where chaotic sticky islands, inside which trajectories show regular features, coexist with a chaotic "sea" between islands. Large scale transport is possible through this chaotic domain [19][START_REF] Zaslavsky | The physics of chaos in hamiltonian systems[END_REF][START_REF] Contopoulos | [END_REF].

These web structures exhibit different shapes which depend on the wave vectors k and amplitudes of the electrostatic waves, and on ω k /ω ce the frequency ratios of electrostatic waves frequencies ω k to the electron cyclotron frequency ω ce [22]. The study of particle transport in these web structures helps to understand the anomalous collisionless transport mechanism in magnetized plasmas. In most previous studies, the formation of stochastic webs and the associated transport were investigated for high wave frequency (ω k ω ce ). Moreover, the perpendicular diffusion coefficient in that limit is calculated by invoking the linear time dependence of the variance [23]. However, deviations from the linear time dependence of variance are frequently observed in case of mixed phase space.

Here, we consider the collisionless transport mechanism of electrons due to the E × B electron drift instability. In magnetized Hall plasmas [24], the E × B electron drift, plasma density, temperature, magnetic field gradients and ion flow are the sources of the E × B drift instabilities or electron cyclotron drift instability [START_REF] Mikhailovskii | Electromagnetic instabilities in an inhomogeneous plasma[END_REF]. This instability is observed in many magnetized plasma devices like magnetrons for material processing [START_REF] Abolmasov | [END_REF], magnetic filters [27], Penning gauges [28], linear magnetized plasma devices dedicated to study cross-field plasma instabilities [29], Hall thrusters for space propulsion and many fusion devices. In Hall thrusters and other devices, this E × B drift instability plays a dominant role in anomalous particle transport. In most of these devices, the electrostatic modes generated by E×B drift instability have very small frequencies compared to the electron cyclotron frequency (ω k ω ce ). Therefore, the resonance condition with the cyclotron harmonics, ω kk v = ω ce , is not satisfied.

In previous works, the chaotic dynamics of test particle (electron) near the anode region in Hall thrusters due to inhomogeneities in magnetic field, and its effect on ionization efficiency and anomalous electron transport, were reported even in the absence of waves [30,31]. In our recent work [32], we present the anomalous transport of electrons due to wave-particle interaction in Hall thruster using a three-dimensional test particle model, even in a uniform magnetic field. But in three space dimensions, the dynamics in the presence of waves is very complicated.

In this paper, we focus on the consequence of the E×B drift instability. We discuss the formation of chaotic web structures and characterize the associated transport properties using a reduced two-degrees-of-freedom Hamiltonian which helps to simplify the original dynamics complexity. In real thrusters, due to the presence of wall boundaries, particles can reflect from the boundary, a process which may destroy the web formation. We find that, in presence of a single background electrostatic wave along with the uniform, static electric and magnetic fields, the trajectories generate web structures, and that, due to the presence of sticky islands, particle transport is super-diffusive. Then we analyse the sticking time statistics by observing a power-law decay of particle presence in each of the sticky sets. Therefore, this work complements our previous findings with a new understanding of the particle transport, which can be applicable in other systems having mixed phase space.

Section II presents the model and its two descriptions (respectively time-dependent and time-independent). Sec. III indicates the numerical method used to integrate the evolution equations. Sec. IV discusses the chaotic web structures generated by the dynamics, Sec. V analyses the transport in these structures, and Sec. VI discusses the effect of sticking to invariant islands on transport. We conclude in Sec. VII.

II. REDUCED HAMILTONIAN DYNAMICS AND THE ELEMENTARY MODEL

A. Fields acting on an electron

We consider a Cartesian coordinate system for the numerical modelling, with x-direction along the magnetic field B 0 , y-direction as E 0 × B 0 drift direction and zdirection along the constant electric field E 0 (see fig. 1 of [32]).

In Hall thruster geometry, unstable low frequency (ω k ω ce ) electrostatic waves are generated due to E × B drift instability. A 3D dispersion relation of this instability for Hall thruster has been derived by Cavalier et al. [33]. The most unstable mode [34] is given by k max ∼ (λ De √ 2) -1 and ω max ∼ ω pi / √ 3, where λ De and ω pi are the electron Debye length and ion plasma frequency respectively. Its propagation angle deviates by tan -1 (k z /k y ) ∼ 10 -15 • from the y-direction near the thruster exit plane. Hence, the wave vector along the z-direction is small (k z ∼ 0.2 k y ), and the electric field along the z-direction is dominated by the stronger constant axial electric field E 0 . Therefore, for simplicity, we remove here the z-variation of the electric field.

For this first investigation, we consider only the fastest growing mode. The total electric field acting on the particle is

E(x, y, z, t) = φ 1 k sin α(x, y, t) + E 0 e z , (1)
with the local phase α(x, y, t) := k x x + k y yω 1 t, where the wave vector k = k x e x + k y e y and the wave angular frequency ω 1 follow the dispersion relation of the E × B electron drift instability [33] and k z = 0. The origin of time is such that α = 0 for x = y = 0, t = 0. The position r = (x, y, z), velocity v, time t and the potential φ 1 are normalized with Debye length λ De , thermal velocity v the , inverse electron plasma frequency ω -1 pe and m e v 2 the /|q e |, respectively. We choose the amplitude φ 1 equal to the saturation potential [34] at the exit plane of the thruster |δφ y,rms | = T e /(6 √ 2) = 0.056 v 2 the . We consider a single mode with (k x , k y , ω 1 ) = (0.001, 0.754, 1.23 • 10 -3 ).

Earlier studies of Hall thruster [33] reported that the comb-like unstable modes found in axial-azimuthal plane get smoothed out as one increases the wave vector k x parallel to the magnetic field ; for large k x values (k x > 0.08), the dispersion relation tends to follow an asymptotic curve corresponding to a modified ion acoustic instability. At small k x values, the real part of the frequency given by the dispersion function oscillates about this asymptotic curve and crosses this asymptotic curve at each resonance (and also one time between each resonance), as can be seen in figure 2 of ref. [33], while the growth rate depends strongly on k x at that point. As our model is not self-consistent, so that the amplitude of the wave is fixed according to experimental measurements in Hall thruster [35], we only use the real part of the frequency of the unstable mode. The mode considered is the one having the largest growth rate and corresponds to one of the resonances. So the real part of the frequency at that point is correctly given by the approximated dispersion relation [34], even though the value of k x chosen here is too small to ensure a smoothed dispersion relation.

From here on, we write ω c for ω ce . In normalized units, the cyclotron frequency ω c = |q e B 0 |/m e = 0.1 ω pe , |q e E 0 |/m e = 0.04 ω pe v the , and the drift velocity v d = E 0 /B 0 = 0.4 v the . Therefore, the y-component of the mode phase velocity is small (ω 1 /k y v d ). As a result, in the Lorentz equation of motion of a particle with mass m and charge q

r = q m (E(r, t) + ṙ × B) , (2)
the electric field E(r, t) has a constant part E 0 along zdirection and a slowly time varying part in (x, y) plane. Eq. ( 2) can be written componentwise, using Eq. ( 1), as

ẍ = qE 1x m sin(k x x + k y y -ω 1 t), (3) ÿ
= qE 1y m sin(k x x + k y y -ω 1 t) + ω c ż, (4) z
= qE 0 m -ω c ẏ, (5)
where E 1x = k x φ 1 and E 1y = k y φ 1 are the amplitude of the x-and y-components of electric field, respectively. Eq. ( 5) can be integrated:

ż + ω c y = qE 0 m t + a, (6)
where a = v z0 + ω c y 0 is a constant of integration, v z0 and y 0 are the particle's initial z-component velocity and position along y-direction, respectively. Substituting ż in Eq. ( 4), and recalling the drift velocity v d = E 0 /B 0 , we reduce the equation of motion of the particle to a system of two equations,

ÿ + ω 2 c y = qE 1y m sin (k x x + k y y -ω 1 t) + v d ω 2 c t + ω c a, ẍ = qE 1x m sin (k x x + k y y -ω 1 t) . (7)
B. Time-dependent Hamiltonian System (7) derives from the Hamiltonian H(p x , p y , x, y, t)

H = p 2 x + p 2 y 2m + m 2 ω 2 c y 2 -(t + A) mv d ω 2 c y +qφ 1 cos (k x x + k y y -ω 1 t) , (8)
where A = (v z0 + ω c y 0 ) / (ω c v d ) is a constant. By means of the generating function

F (P x , P y , x, y, t) = P x x + (P y + v d ) (y -(t + A)v d ) , (9)
we change to new variables (P x , P y , X, Y ) in a frame moving with a constant velocity v d along the y-direction (which we call a "drifted frame" for figures),

X = ∂F ∂P x = x, Y = ∂F ∂P y = y -(t + A) v d , p x = ∂F ∂x = P x , p y = ∂F ∂y = P y + v d , ∂F ∂t = -(P y + v d ) v d . (10)
Using these new coordinates [START_REF]Stochastic differential equations[END_REF], the new Hamiltonian (after removing terms irrelevant to the motion) and the equations of motion read

K (P x , P y , X, Y, t) = P 2 x + P 2 y 2m + m 2 ω 2 c Y 2 + qφ 1 cos α, Ẍ = qφ 1 m k x sin α, Ÿ + ω 2 c Y = qφ 1 m k y sin α, (11)
with

α = k x X + k y Y + (v d k y -ω 1 )t + ζ and where ζ = k x v d A is constant.
The dimensionless equations of motion are obtained using the dimensionless variables

X = k x X + ζ, Y = k y Y , t = ω c t. Introducing the new notation β = k x /k y , ε = qφ 1 k 2 y /(mω 2 c ) and ν 1 = (v d k y -ω 1 )/ω c , we obtain the dimensionless equations d 2 X dt 2 = εβ 2 sin (X + Y + ν 1 t ) , d 2 Y dt 2 + Y = ε sin (X + Y + ν 1 t ) . (12)
In this paper, we solve Eqs (12) numerically using a second order symplectic scheme. The dynamics involves two degrees of freedom with a time-periodic dynamics (with period 2π/ν 1 ), so that the effective phase space is 5dimensional. The coordinate X admits periodic boundary condition (with period 2π), whereas Y runs over the real line.

The dynamics depends on three parameters, ε, β and ν 1 . For ε = 0, X is ballistic and Y is a harmonic oscillator, in agreement with the well-known solutions for particle motion in stationary, uniform fields E 0 and B 0 . Note that, in Hall thrusters, B 0 is radial and E 0 is axial, so that the drift is azimuthal. The coordinates y and Y are thus defined on circles, while x and X are actually bounded by the inner and outer cylindrical chamber walls. The origin for Y and X are determined by the initial conditions (y 0 , v z0 ) and the phase convention for the electrostatic mode, respectively.

C. Time-independent Hamiltonian

A time-independent Hamiltonian can be derived by means of a Galileo transformation along X with velocity ν 1 ω c /k x . With the generating function and change of variables

F = P x - ν 1 mω c k x X + ν 1 ω c t + ζ k x + P y Y, X = X + ν 1 ω c k x t + ζ k x , P x = P x - ν 1 mω c k x , Y = Y, P y = P y , ∂F ∂t = ν 1 ω c k x P x - ν 1 mω c k x , (13)
the Hamiltonian (up to terms irrelevant to the motion) and the equations of motion can be written as

K (P x , P y , X , Y) = P 2 x + P 2 y 2m + m 2 ω 2 c Y 2 +qφ 1 cos (k x X + k y Y) , Ẍ = qφ 1 m k x sin (k x X + k y Y) , Ÿ + ω 2 c Y = qφ 1 m k y sin (k x X + k y Y) . (14)
Setting X = k x X , Y = k y Y and t = ω c t, the equations of motion ( 14) reduce to

d 2 X dt 2 = εβ 2 sin (X + Y ) , d 2 Y dt 2 + Y = ε sin (X + Y ) . (15)
Eq. ( 15) is solved numerically for various parameters and initial conditions. This dynamics depends on two parameters only, ε and β. It involves two degrees of freedom, with the coordinate X obeying periodic boundary condition (with period 2π), whereas Y runs over the real line. As the dynamics is autonomous, it preserves the "energy" K. Therefore, trajectories stay on smooth 3-dimensional surfaces, and they may be visualised by means of 2-dimensional Poincaré sections.

While coordinates (x, y) and (X, Y ) are related with the Hall thruster chamber, coordinates (X , Y) simplify further the dynamics, provided one does not worry about boundary conditions. Therefore, we use both the timeindependent and the time-dependent representations in the following discussions.

For ε = 0, viz. in absence of electrostatic wave, the dynamics is integrable. Its dimensionless actions are the linear momentum dX /dt with angle the position X periodic with period 2π in agreement with the boundary condition, and the gyration energy R 2 /2 = (Y 2 + (dY /dt ) 2 )/2 (divided by the cyclotron frequency, which is 1) with angle the gyrophase in the (Y , dY /dt ) plane. In presence of the electrostatic wave, for small ε, these actions generate two adiabatic invariants. For β also small, the actions evolve on different time scales (in terms of the dimensionless t ), namely ε -1 for the oscillations of Y and ε -1 β -2 for the nearly-ballistic motion of X .

III. NUMERICAL METHOD

Because the right hand sides of Eqs ( 12) and ( 15) depend on space, the infinitesimal generators for both velocity and position equations do not commute, and one uses a time-splitting numerical integration scheme. Since the dynamics is Hamiltonian and we are interested in long-time evolution, we choose a symplectic scheme, which guarantees preservation of the hamiltonian structure exactly and ensures over very long time the conservation of a hamiltonian close to the original one [36,37].

The positions are advanced with the map r(t + ∆t) = T v,∆t (r(t)) = r(t) + v∆t, and the velocities are advanced in the form

v x (t + ∆t) = T Ex,∆t (v x (t)) = v x (t) + εβ 2 sin(X + Y) ∆t and v y (t + ∆t) = T Ey,∆t (v y (t)) = v y (t) + (ε sin(X + Y) -Y) ∆t.
As a result, we use a second-order symmetric leapfrog scheme, which evolves [START_REF] Méndez | Stochastic Foundations in Movement Ecology -Anomalous Diffusion[END_REF] as the map

r(t + ∆t) v(t + ∆t) = A r(t) v(t) , (16)
A = T v,∆t/2 • T E,∆t • T v,∆t/2 . (17)
We evolve Eqs ( 12) similarly, with T E,∆t evaluated at midstep t + ∆t/2.

IV. STOCHASTIC WEB STRUCTURE

The values of (ε, β 2 , ν 1 ) in a Hall thruster device for the fastest growing mode (k x = 0.001, k y = 0.754, ω = 1.23 • 10 -3 ) are ε = 3.21, β 2 = 1.75 • 10 -6 and ν 1 = 3. As ε is not small, the gyration action is definitely not conserved, as will be seen in the phase space plots. However, εβ 2 ∼ 5 • 10 -6 so that the ballistic action is almost conserved over times t ∼ 10 5 ω -1 c = 10 6 ω -1 pe . Moreover, k x = 0.01, so that εβ 2 = 5.6 • 10 -4 which also ensures conservation of the action for long time. We checked that setting k x = 0.01 with all other parameter kept constant gives similar phase space plots as k x = 0.001. For higher values like k x ≥ 0.1, the action is not conserved over long times, which changes the phase space structure.

Here we first focus on transport for web structures with three-fold rotational symmetry (ν 1 = 3) as in the Hall thruster geometry, and its harmonic the six-fold rotational symmetry (ν 1 = 6) with respect to gyration variables (Y , v y ). By "rotational symmetry" we mean that, if one rotates the (Y , v y )-space web structure by angle 2π/n around the origin, then it generates an identical phase space structure. To assess the importance of having an integer value for the forcing frequency ν 1 , we also consider the non-resonant value ν 1 = 1.39. For the timeindependent description, the initial velocity Ẋ0 plays a similar role, and we also contrast the integer values 3 and 6 with the rational value 3.5.

For small values εβ 2 1, the dynamics generate a spiral stochastic web [3] in the three-dimension space (X , Y , Ẏ ). Moreover, for ε 1, the dynamics is nearly integrable and the chaoticity vanishes. In the section cut by a plane X = const, the web is a system of concentric circles and of straight lines passing through the co-ordinate origin in the (Y , Ẏ ) plane. The cells of the web form concentric "belts" around the coordinate origin. These trajectories are called "web trajectories", due to their special structural shape. Birkhoff's theorem ensures that, for small, nonzero ε and β, the Poincaré map has O type (elliptic) fixed points and X type (hyperbolic) ones ; regular trajectories are organized in islands around the O points, while chaotic ones wander through the heteroclinic and homoclinic tangles associated with the stable and unstable manifolds of the X points [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF].

Fig. 1 presents the web network for different initial particle positions, for ε = 0.5, β 2 = 1.75 • 10 -6 and ν 1 = 3. Due to small ε value, the Hamiltonian is nearly integrable and KAM theory applies. When X varies, as a result of detuning from longitudinal cyclotron resonance (k x Ẋ = ω ce ), the separatrix network and, consequently, the stochastic web rotate about the coordinate origin with angular velocity Ẋ ∼ β 2 . But the shape and area of the shells remain constant, and a particle that is executing rapid rotation inside a cell sufficiently far from the separatrix of the averaged motion will never intersect the stochastic web and will move regularly.

For high values of β 2 or Ẋ 0 , X (t) varies rapidly and the chaotic tangle network also rotates rapidly compared to the case with slowly varying X . This precludes the presence of islands in the (Y , Ẏ ) phase space for long time evolution. Fig. 2 presents the rotation of a regular web trajectory due to the high β 2 value.

In our simulations, we evolve the dynamics of 1024 particles having initial Gaussian velocity distribution with unit standard deviation along the y-direction and covariance Ẋ 0 Ẏ 0 = 0. Along the x-direction, we consider a very small standard deviation σ x = 0.001, and choose the square of wave vector ratio β 2 ∼ 10 -6 . We first consider the time-dependent dynamics in Sec. IV A, then the time-independent cases in Sec. IV B.

A. Time-dependent Hamiltonian

For the time-dependent Hamiltonian, the dynamics depends on (ε, β 2 , ν 1 ). Parameter ε expresses the ratio of bounce frequency to cyclotron frequency, β expresses the ratio of the parallel and perpendicular components of the wave electric field, and ν 1 is the normalized frequency of the electrostatic wave in the drifted frame.

A large value of β 2 or Ẋ 0 causes the dynamics detuning from the longitudinal resonance condition, k x Ẋ = ω c , where is an integer. Therefore, the orbits in the web structure drift more rapidly away from their initial action values, covering the entire phase space inside the web and destroying islands. In our present simulation, β 2 ∼ 10 derived by Karney [22] must hold. For all three sets of parameters, this stochasticity criterion is indeed satisfied.

We first plot the stroboscopic Poincaré section in the (Y , P y ) plane of a single particle trajectory at times t = 2nπ/ν 1 , where n = 0, 1, 2.... The parameter ε determines the radius of the stochastic web. The value of ν 1 determines the shape of the web structure. For integer ν 1 , we observe a web structure with ν 1 -fold rotational symmetry (Figs 3 and4). For parameters (ε, β 2 , ν 1 ) = (0.69, 1.83 • 10 -5 , 1.39) with non-integer ν 1 , the dynamics generates a Halloween mask-like, deformed three-fold web structure (Fig. 5). Along Y , the dynamics Eq. ( 12) has two time scales, one associated with the electrostatic wave (with period 2π/ν 1 in the drifting frame) and the other associated with a simple harmonic oscillator with period 2π. Therefore, an integer value of ν 1 causes resonance between these two

-3 -2 -1 0 1 2 3 y[k -1 y ](alongE 0 × B 0 ) -4 -3 -2 -1 0 1 2 3 v y [ω c k -1 y ]
(yv y ) in drifted frame For fixed values of ε, β 2 and ν 1 , any initial condition (X 0 , Y 0 , Ẋ 0 , Ẏ 0 ), within the chaotic domain of the stochastic web, generates a similar web structure, and the particles with initial conditions outside the web structure and well inside the sticky islands (regions with no points in the web structures) generate regular trajectories.

B. Time-independent Hamiltonian

In the time-independent Hamiltonian, Ẋ = Ẋν 1 and the dynamics depends on (ε, β 2 ) only. In this case, for any initial condition (X 0 , Y 0 , Ẋ 0 , Ẏ 0 ) within the chaotic domain of the stochastic web, the shape of the web structure depends on the initial Ẋ 0 values. For different Ẋ 0 values, the trajectory lies on different energy surfaces K = constant. Thus, particles with different initial conditions generate different web structures. Since β 2 ∼ 10 -6 in this study, the motion along the X direction is almost ballistic, Ẋ ∼ = constant. Therefore, we can generate Poincaré section plots by taking sections at X = n2π, where n is an integer.

Figs. 6 and 7 display the stroboscopic plot of the timeindependent dynamics Eq. ( 15) with (ε, β 2 ) = (3.21, 1.75• 10 -6 ) and two different initial velocities along the X direction, Ẋ 0 = 3 and 3.5 respectively. For integer values of Ẋ 0 as in Fig. 6, the dynamics generates web structures similar to those generated in the Poincaré section plot for the cases of time-dependent dynamics (12) with same integer values of ν 1 . Fig. 7 presents the stroboscopic plot of a particle with v 0x = 3.5, which corresponds to a higher-order resonance (7/2) : for fractional values of Ẋ 0 , the stroboscopic plot generates different structures, because each of the different initial conditions lies on a different energy (K = constant) surface. Therefore, the web structures in the time-independent dynamics highly depend on the initial conditions of the particle.

V. TRANSPORT PROPERTIES

To characterise the transport properties, we consider a simple observable. Previous studies [START_REF] Bouchara | [END_REF][START_REF] Leoncini | [END_REF] dependent one-degree-of-freedom Hamiltonian systems focused on the norm of velocity ( ṗ, q) in phase space, where p, q are canonical co-ordinates. Here, we consider the arc length s of the trajectory in position space only, or, in dimensionless variables of Eqs (12),

S (t) = t 0 dX 2 + dY 2 . (18)
Numerically, we consider the global average speed along the trajectory of a typical particle i

vi (n) = 1 n∆t n-1 k=0 [∆X i (t k )] 2 + [∆Y i (t k )] 2 , ( 19
)
where k is the timestep index, with coordinate increments

∆X i (t k ) = X i (t k+1 ) -X i (t k ) , ( 20
) ∆Y i (t k ) = Y i (t k+1 ) -Y i (t k ) . (21)
In eq. ( 19), the discretized form of the integral ( 18) is used.

When conditions are met so that we can apply the central limit theorem, the distribution of quantity

Z i := 1 √ n n k=1 (v i (k) -v ),
gives a normal distribution function with mean 0 and variance σ 2 , where the distribution increases with rate √ n as the variance shrinks. We define ρ n (v) as the sampling density of the distribution of vi (n)'s. Therefore, good ergodic properties of the dynamics (12) would include the convergence of ρ n towards a Dirac distribution for n → ∞, in which case the support v of the limit is the time average of the vi (n)'s,

v i (k) = [∆X i (t k )] 2 + [∆Y i (t k )]
lim n→∞ ρ n (v) = δ(v -v ), (22)
and, almost surely with respect to the initial condition (viz. index i),

v := lim n→∞ vi (n). ( 23
)
In the present case, what is understood by ergodic properties is that, when considering initial conditions in the chaotic sea, they will almost surely give rise to the same natural ergodic measure, meaning that they equally sample all the accessible chaotic domain of phase space when the dynamics evolve for sufficiently long time.

One method to assess the convergence of ρ n is to look at how fast its maximum value ρ max (n) diverges towards +∞ with n. In order to characterize this speed, one can typically expect that a scaling applies to increments in (19), so that one may expect

ρ max (n) ∼ n α , (24)
where the exponent α characterises the nature of the transport. If increments in (19) are quite independent and a central limit theorem applies, transport is diffusive and α = 1/2. For α > 1/2 it is sub-diffusive, and for α < 1/2 it is super-diffusive. Indeed, instead of considering the global average quantity (speed) vi , one may consider the arc length S i (t) = n k=1

[∆X i (t k )] 2 + [∆Y i (t k )] 2
= nv i in order to characterise the transport. Here, when we can apply the central limit theorem, one expects the variance of the distribution of S i 's to grow like √ n. One can also calculate different moments M q of the distribution of S i 's and extract the characteristic exponent µ from

M q = |S i -S i | q ∼ t µ(q) . (25)
The second order exponent q = 2 may be related [START_REF] Bouchara | [END_REF] to the variance of the displacement in space, σ 2 ∼ t µ (2) . For a chaotic system of diffusive type, µ(2) = 1, while µ(2) = 1 for an anomalous transport. Moreover, the area under the distribution function is constant, which implies that the variance grows faster in presence of fat tails than in the diffusive case. Therefore, µ(2) > 1 implies super-diffusive transport and µ(2) < 1 implies sub-diffusive transport. The exponent α associated with ρ max of the speed distribution is different from this µ(2) which is associated with the moment of the position-displacement distribution, but in the same spirit and due to the fact that normalized distributions have constant area, the existence of fat tails will imply that the maximum ρ max will not grow so fast, and thus α < 1/2 for ρ max indicates super-diffusive transport, whereas α > 1/2 indicates sub-diffusive transport. Both the scaling parameter α and µ(2) characterise the diffusion property, and one may derive [START_REF] Leoncini | [END_REF] a relation α = 1µ(2)/2 between them.

The non-diffusive aspect of transport implies that the usual central limit theorem does not apply. However, the motion of trajectories in the chaotic web does lose correlations, and each trajectory typically generates a similar picture of the web. Thus, one can still view a trajectory as a sequence of pieces which are essentially mutually independent, but with duration in time and extent in space which do not meet the finite variance assumption necessary for the gaussian central limit. In other words, the successive pieces generate a process obeying a scaling law, viz. it is infinitely divisible, and their (properly rescaled) sum has a Lévy distribution characterized by the scaling exponent [START_REF] Balescu | Statistical dynamics -Matter out of equilibrium[END_REF][START_REF] Zaslavsky | The physics of chaos in hamiltonian systems[END_REF].

A. Time-dependent Hamiltonian

In Figs. 8 and9, we plot the distribution of

u i,n = vi (n) -v (26)
for two different web structures, with ν 1 = 3 and ν 1 = 6, respectively. One can calculate the arc length for a timeindependent dynamics also, but in the time-dependent dynamics, the parameter ν 1 = (v d k yω 1 )/ω c is important in Hall thrusters. It expresses the frequency of the electrostatic modes generated by the E × B instability, in a frame which moves with the drift velocity v d along the E 0 × B 0 direction. As our study is motivated by the anomalous chaotic transport in Hall thruster devices, we choose the time-dependent dynamics for characterising the transport for ν 1 = 3 and ν 1 = 6.

To characterize the transport, we consider two different stochastic webs corresponding to values (3.21, 1.75 • 10 -6 , 3) and (3.21, 1.75 • 10 -6 , 6) for parameters (ε, β 2 , ν 1 ). We evolve the equations of motion (12) for 1024 particles with all initial conditions inside the chaotic domain, we calculate the arc length for each particle trajectory for a long time evolution (10 8 ω -1 pe ) using a timestep value ∆t = 3.33•10 -3 ω -1 c = 3.33•10 -2 ω -1 pe in simulation, and we generate the distribution of the global average speed.

We plot the distribution at different times in Figs. 8 and9 for both cases. To avoid nonphysical peaks in the distribution of ρ n , the length n of the time sequence should be sufficiently long for the dynamics to reach a saturation state, i.e. for the Poincaré section of each particle's trajectory to sample the whole phase space reach of the web. Here, we construct the distribution functions ρ n at times t n ≥ 800 ω -1 pe . In the plot, the strong sharp peaks are associated with the stickiness phenomenon, by which a trajectory may remain for a long time close to the regular islands. The number of sharp peaks depends on the number of resonance generating sticky islands within the web structures, which we further discuss in the next section. We obtain more peaks for ν 1 = 6 (Fig. 9) than for ν 1 = 3 (Fig. 8). In the figure, log 10 t = log 10 (nτ P ), where τ P is the average time between two successive Poincaré sections. The relative magnitude of these sharp peaks decreases as n → ∞, because the contribution from the chaotic domain becomes large compared to the contribution from the sticky regular trajectories as we consider a longer time evolution.

In both Figs. 8 and9, the distribution after time t n ∼ 10 7 ω -1 pe (yellow line) has almost zero relative strength of the sharp peaks, compared to the height of the smoother distribution. Stickiness generates a memory effect and Lévy flights [40]. In absence of these sticky trajectories, the transport is purely diffusive and the exponent α takes the value 1/2. In the presence of these sticky trajectories, the transport will be anomalous.

To measure α, we find the value of ρ max from the local maximum of the central smooth flat peak location which corresponds to the bulk and particles evolving in the chaotic domain, and not the sticky domains. In Fig. 10, we plot the time evolution of ρ max for both cases ν 1 = 3 and ν 1 = 6. From the curve fitting, we obtain two different values α 3 = 0.17 and 0.39 for the case ν 1 = 3 (magenta dashed line and data dots), and α 6 = 0.15 and 0.33 for ν 1 = 6 (black). Both values of α in both cases are below 0.5. Thus, the diffusion is anomalous and superdiffusive.

After a longer time evolution, most of the particles spend more time exploring the chaotic region of the stochastic web and sampling the ergodic measure. It appears that the contribution from the sticky islands "decreases" in comparison with the contribution from the chaotic domain, and that the "diffusion" rate increases at longer time (t n > 10 5 ω -1 pe ). Note however that, even for large times, the exponent being smaller than 1/2 implies that the average speed fluctuations viv do not approach a Gaussian distribution, hence they do not obey the central limit theorem over this time scale, and transport is superdiffusive.

B. Time-independent Hamiltonian

Similarly, we analyse transport for a stochastic web structure generated from the time-independent Hamiltonian with the corresponding arc length

S (t) = t 0 dX 2 + dY 2 . ( 27
)
Transport in web with simple rational parameter: The trajectories of 1024 particles are computed numerically up to time 10 8 ω -1 pe with time step ∆t = 5.5 • 10 -3 . All particles are initially randomly distributed in (X , Y ) plane within -π ≤ X 0 ≤ π and -π ≤ Y 0 ≤ π. Their initial velocities along the Y -direction are drawn from a Gaussian distribution with unit standard deviation. Along the X -direction, we consider three different values Ẋ 0 = 3, 3.5 and 6 (simple rational) in order to analyse the transport in three different web structures generated from the time-independent Hamiltonian [START_REF] Méndez | Stochastic Foundations in Movement Ecology -Anomalous Diffusion[END_REF]. For all three cases, we consider ε = 3.21 and β 2 = 1.75 • 10 -6 . Figs. 11, 12 and 13 present the distribution of the average speed at different times for all three cases. Similarly to the time-dependent cases, sharp peaks in the distribution of average speed appear due to the presence of the sticky islands, and the number of sharp peaks is larger for the web structures with six-fold rotational symmetry, Ẋ 0 = 6, than for the three-fold rotational symmetry, Ẋ 0 = 3. In the case Ẋ 0 = 3.5, as seen in Fig. 7, the number and area of sticky islands are smaller than for Ẋ 0 = 3. Therefore, the height in the smooth part of the distribution, due to the chaotic domain, is larger for Ẋ 0 = 3.5.

To estimate the exponent values from Eq. ( 24), we plot in Fig. 14 the time evolution of ρ max for all three cases. For Ẋ 0 = 3 (magenta) and 6 (black), the plots are similar to the time-dependent cases. From the curve fitting, we obtain two different values of α in two different regimes of the plots, α 3 = (0.15, 0.25) and α 6 = (0.18, 0.25). In both cases, the transport is anomalous of super-diffusive type. The values for the shorter time span (t < 5 • 10 5 ω -1 pe ) are very close to the values that are recovered for the time-dependent dynamics.

For fractional values of Ẋ 0 , the number and area of the sticky islands are smaller than for the other two cases. Most of the region within the stochastic webs is part of the chaotic domain, therefore the height of the distribution increases at a higher rate than in the other two cases, as we increase the value of n. For Ẋ 0 = 3.5 = 7/2, we find higher exponent values, α 3.5 = (0.23, 0.46). At short time, the relative contribution from sticky trajectories is significantly large, which reduces the exponent value to α 3.5 = 0.23 ; in contrast, for longer time t > 5 • 10 4 ω -1 pe , the contribution from the chaotic region dominates over the contribution from sticky islands, and sharp peaks almost disappear, which increases the exponent value to α 3.5 = 0.46, so that the transport becomes closer to a diffusive type.

Transport in Halloween mask like web: For Ẋ 0 = 1.39 (further away from a simple rational) and (ε, β 2 ) = (0.69, 1.83 • 10 -5 ), we also draw 1024 initial conditions in the chaotic part of the domain defined by -π ≤ X 0 ≤ π, Ẏ 0 a Gaussian random number with expectation 0 and standard deviation 1, and Y 0 outside the islands (typically, 0.1 π ≤ |Y 0 | ≤ 0.6 π). With these parameters and initial conditions, the same analysis applies, as seen from the peaks in the distribution in Fig. 15 and from the slopes α 1.39 = (0.32, 0.42) in Fig. 16. In the next section, we discuss this change of exponent values more quantitatively.

VI. EFFECT OF STICKY ISLAND ON ANOMALOUS TRANSPORT: CHANGE OF α

In this section, we will attempt to quantify the impact of a sticky set on the transport, recalling that α must be 0.5 when the central limit theorem can be applied, whereas it is anomalous if α = 0.5. The superdiffusive or subdiffusive nature of transport is measured by the weight of the tails in the distributions : (i) if there are fat tails (compared to the diffusive bulk of the distribution), the maximum ρ max of the distribution will grow slower than power 1/2 and we expect superdiffusion, while (ii) if there are thinner tails, the maximum ρ max of the distribution will grow faster than power 1/2 and we expect subdiffusion. The presence of the sharp peaks in the distribution of ρ n (v), Figs 8,[START_REF] Balescu | Statistical dynamics -Matter out of equilibrium[END_REF][START_REF] Chirikov | [END_REF]12,13 and 15, increases the effective weight of the tail of the distribution, which makes exponent α < 1/2.

Each of the sharp peaks is related to the presence of sticky islands in the phase space. Thus, one can select the portions of trajectories contributing to each peak and locate them in the Poincaré map. Specifically, given a trajectory z i (t) (with initial condition labeled 1 ≤ i ≤ N = 1024) over a long time span T max , and a sticking time T s , we fix a delay τ , and consider the arc lengths ∆S i,m (T s ) = S i (T s +mτ )-S i (mτ ) of the portion of z i over [mτ, T s + mτ ], for 0 ≤ m ≤ M -1 for some large M = (T max -T s )/τ , where T max is the maximum value of the time evolution, T max = 10 8 . For a moderate value of τ (say, ∼ 200), this method generates M N time sequences of length T s which we analyse. Indeed we generate the distribution functions using the same technique. Therefore, T s is any time when the distribution function is generated. The minimum value of T s and τ should be such that the dynamics become sufficiently ergodic within that time, i.e. it covers the entire chaotic domain of the phase space. Since we plot T s in log 10 scale, we choose the values T s = 200, 400, 800, 1600, 3200, 6400.... The maximum value of T s is determined by the relation M ≥ 50, such that the total number of data points for each distribution M N > 5 • 10 3 , which constructs a smooth distribution function at any time T s .

In particular, we extract the points (X , Y , Ẏ ) associated with each vi that generates a sharp peak in the distribution and we plot their Poincaré section in the (Y , dY /dt ) plane, which generates the particular sticky set. This is done in Fig. 17 for the web structures with three-fold symmetry for the trajectories with sticking duration T s such that log 10 (ω pe T s ) = 2.92. There are three peaks, Pk 1 at u n = 1.0, Pk 2 at u n = 3.0 and Pk 3 at u n = 5.0, in Fig. 11, each with a finite width. For each peak, we identify the trajectories contributing to the peak, and plot their Poincaré section, whereby the sticky regions emerge. In Fig. 17, the sticky regions denoted by blue, red and black dots are associated with the peaks Pk 1 , Pk 2 and Pk 3 , respectively.

In Fig. 13, the distribution of average speed for the web with six-fold symmetry has seven peaks, namely Pk 1 at u n = -3.6, Pk 2 at u n = -1.5, Pk 3 at u n = 0.5, Pk 4 at u n = 2.3, Pk 5 at u n = 4.2, Pk 6 at u n = 6.2 and Pk 7 at u n = 8.1. In a similar way, we locate the sticky region in the phase space for each of these peaks for the same duration T s , such that log 10 (ω pe T s ) = 2.92. In Fig. 18, the blue, red, green, magenta, cyan, yellow and black dots identify the sticky regions associated with the peaks Pk 1 to Pk 7 , respectively. Thus, all peaks in the distribution plots are associated with different sticky sets. In Fig. 19, we similarly identify the sticky sets associated with the peaks for the stochastic web with Halloween mask like structure of Eq. ( 15) for ε = 0.69, β 2 = 1.83 • 10 -5 and Ẋ 0 = 1.39. Sticky sets are special in the dynamics, because they influence transport in the chaotic domain. Therefore, they are not islands in which the trajectory remains forever : their trajectories leave them to wander further through the chaotic domain. This implies that a group of trajectories initially in the sticky set will leak into the main chaotic sea -but slowly enough for their stickiness to show up. Geometrically, sticky regions correspond to tangles bordering islands. If the initial condition of a particle is well inside an island, the dynamics of the par-FIG. 18: Localization of seven different sticky regions in the stochastic web with six-fold symmetry of Eq. ( 15) for ε = 3.21, β 2 = 1.75 • 10 -6 and Ẋ 0 = 6.0. These different sticky regions are associated with seven different peaks in the distribution plot of average speed (Fig. 13) with sticking duration Ts such that log 10 (ωpeTs) = 2.92, which are Pk1 at un = -3.6 (blue dots), Pk2 at un = -1.5 (red dots), Pk3 at un = 0.5 (green dots), Pk4 at un = 2.3 (magenta dots), Pk5 at un = 4.2 (cyan dots), Pk6 at un = 6.2 (yellow dots) and Pk7 = 8.1 (black dots), respectively.

ticle will be essentially regular : it never exits from the islands chain and generates a sharp peak in the distribution which never becomes empty with time. Those regular trajectories do no contribute to the anomalous transport. We exclude all the regular trajectories to remove those peaks by considering all the initial condition inside the chaotic domain.

Each sticky set leaks into the main chaotic domain at its own rate, which we discuss below. However, because our dynamics is hamiltonian, we also know from the Poincaré recurrence theorem that any trajectory will, after a sufficiently long but finite time, return to a state arbitrarily close to its initial state. Therefore, the leakage from a sticky set is a transient process and, in general, the trajectories having left such a set behave afterwards like any other trajectory of the chaotic domain. It will almost surely return, at a later time, close to its initial state, but the statistics of this Poincaré recurrence time and its sensitivity to the initial data are another issue, which we do not discuss.

To understand the influence of stickiness on anomalous transport, we now consider the web structure with three-fold symmetry, and investigate the change of each peak for increasing evolution time t = n∆t. This analysis makes it possible to see the time evolution of the particles trapped in the corresponding islands. From the distribution ρ n (v) of average speed, one can count the number of data points that contribute to each specific peak at different times t. Then one estimates how long the particles are sticking to each specific island, by monitoring the change of area localized under each of those peaks as a function of n. Therefore, this area yields the weight of sticking to that particular island, until at least t 1 = n∆t, which can be written as

w Pk (t 1 , T max ) = (T max -t 1 ) -1 Tmax t1
ρ Pk (t)dt, (28) where Pk is the index for each peak, ρ PK (t) is the area under the peak Pk at time t, and the statistics are gathered from a "very long" run [0, T max ]. Under an ergodic assumption [40,[START_REF] Meziani | Chaos, Complexity and Transport[END_REF], this weight would enable one to estimate the probability that a trajectory would stick to island Pk for at least the duration t 1 . For large sticking time, a self-similar behaviour in the small scales in phase space near the island will be associated with a power law decay with an exponent γ, w Pk (T s , T max ) ∼ T 1-γPk s ,

where we consider ρ Pk (t) ∝ t -γ and replace t 1 with T s . This integration is valid only for γ PK > 1.

In order to analyze the sticking-times statistics, we count the number of data points sticking to each island and plot them in logarithmic scale versus the duration.

One way to do is by integrating the distribution along the u axis, at the particular sharp peak location. But when we locate, in the (Y , Ẏ ) Poincaré section, the dots associated with the sharp peaks, we found some dots actually situated within the main chaotic domain, and moreover the sticky sets overlap with each other on this one-variable (u) axis. Therefore, to find the actual ρ Pk value associated with the sticky sets only, one should count the dots situated only near the chaotic tangle bordering the island. Each sticky set leaks with time and mixes into the more regular chaotic domain. So, we define the border of the web network by an invariant surface with dimensionless radius R = Y 2 + (dY /dt ) 2 , which is associated with the gyration action. The number density within each sticky set decreases rapidly away from the surface R = constant. Therefore, to count the dots belonging only to the sticky sets, we denote the maximum and minimum value of R by R out and R in respectively for each peak. We consider three different annular domains in phase space, one for each peak, using the radius R . In the (Y , dY /dt ) section plane, Pk 1 is associated with the annulus with inner radius R in = 9.4 and outer radius R out = 13.5 ; similarly, Pk 2 with R in = 12.23 and R out = 16.76, and Pk 3 with R in = 15.72 and R out = 20.0. From the data set associated with each peak, we identify those points which satisfy R in ≤ k y Y 2 + Ẏ2 /ω 2 c ≤ R out . We perform this counting for all the n values to obtain W Pk (n) ∼ = M N w Pk , we plot log 10 (W Pk ) vs. log 10 (ω pe T s ), and read the exponent γ Pk from the slope according to Eq. (29). During the counting process, one should keep in mind that the length of analysed data set M N should be constant for all the time T s . For the weakest peak Pk 3 , the exponent changes to γ Pk3 = 2.15 when log 10 (ω pe T s ) = 5.7. Since the exponent value for the strongest peak changes when log 10 (ω pe T s ) = 5.3, the strength of this peak starts to decrease (leak) faster, which helps to increase the maximum ρ max of the average distribution ρ n (v) at a faster rate. Therefore, the value of the transport exponent α, from Eq. ( 24), increases after log 10 (ω pe T s ) = 5.3, which is also observed in Fig. 14.

FIG. 1 :

 1 FIG. 1: Web network generated by Poincaré section of different trajectories with different initial Ẏ 0 values and X 0 = 0 for ε = 0.5, β 2 = 1.75 • 10 -6 and ν1 = 3.

  -6 and | Ẋ 0 | 1 which induces a slow drift of the trajectory. We evolve Eqs (12) numerically for three different sets of parameters (ε, β 2 , ν 1 ) = (3.21, 1.75 • 10 -6 , 6.0), (3.21, 1.75 • 10 -6 , 3.0) and (0.69, 1.83 • 10 -5 , 1.39), respectively. The frequency of the electrostatic mode ω 1 is very small with respect to ω c . In the drifted frame, one has ν 1 > ω c , and the chaoticity criterion ε > 0.25 ν 2/3 1

  FIG.3: Poincaré section of a single trajectory of (12) for ε = 3.21, β 2 = 1.75 • 10 -6 and ν1 = 6.

FIG. 5 :FIG. 6 :

 56 FIG.5: Poincaré section of a single trajectory of (12) for ε = 0.69, β 2 = 1.83 • 10 -5 and ν1 = 1.39.

FIG. 7 :

 7 FIG. 7: Stochastic web of a single trajectory of (15) for ε = 3.21, β 2 = 1.75 • 10 -6 and Ẋ 0 = 7/2 = 3.5.

6 FIG. 10 :

 610 FIG. 10: Evolution of ρmax versus n for the time-dependent Hamiltonian with ν1 = 3 (magenta) and ν1 = 6 (black).

44 FIG. 11 : 44 FIG. 12 :

 44114412 FIG. 11: Distribution of average speed for the stochastic web with Ẋ 0 = 3 and (ε, β 2 ) = (3.21, 1.75 • 10 -6 ) at ωpe t = 8 • 10 2 , 13 • 10 3 , 21 • 10 4 , 85 • 10 4 , 34 • 10 5 , 13 • 10 6 and 5.4 • 10 7 .

44 FIG. 13 :

 4413 FIG. 13: Distribution of average speed for the stochastic web with Ẋ 0 = 6 and (ε, β 2 ) = (3.21, 1.75 • 10 -6 ) at ωpe t = 8 • 10 2 , 6.7 • 10 3 , 1.0 • 10 5 , 4.2 • 10 5 , 1.7 • 10 6 , 6.9 • 10 6 and 2.7 • 10 7 .

6 FIG. 14 :

 614 FIG.14: Evolution of ρmax versus n for the time-independent Hamiltonian for (ε, β 2 ) = (3.21, 1.75 • 10 -6 ) with Ẋ 0 = 3 (magenta), Ẋ 0 = 6 (black) and Ẋ 0 = 3.5 (green).

log10(t) = 7. 47 FIG. 15 : 42 FIG. 16 :

 47154216 FIG. 15: Distribution of average speed for the stochastic web with Ẋ 0 = 1.39 and (ε, β 2 ) = (0.69, 1.83 • 10 -5 ) at ωpe t = 7 • 10 3 , 1.1 • 10 5 , 4.6 • 10 5 , 1.8 • 10 6 , 3.7 • 10 6 , 1.4 • 10 7 and 2.9 • 10 7 .

FIG. 17 :

 17 FIG. 17: Localization of three different sticky regions in the stochastic web with three-fold symmetry of Eq. (15) for ε = 3.21, β 2 = 1.75 • 10 -6 and Ẋ 0 = 3.0 as in Fig.6. These different sticky regions are associated with three different peaks in the distribution plot of average speed (Fig.11) when log 10 (ωpet) = 2.92, which are Pk1 at un = 1.0 (blue dots), Pk2 at un = 3.0 (red dots) and Pk3 at un = 5.0 (black dots), respectively.

FIG. 19 :

 19 FIG. 19: Localization of six different sticky regions in the stochastic web with Halloween mask like structure of Eq. (15) for ε = 0.69, β 2 = 1.83 • 10 -5 and Ẋ 0 = 1.39. These different sticky regions are associated with six different peaks in the distribution plot of average speed (Fig. 15) with sticking duration Ts such that log 10 (ωpeTs) = 3.86, which are Pk1 at un = -1.75 (red dots), Pk2 at un = -7.5 (blue dots), Pk3 at un = 0.4 (magenta dots), Pk4 at un = 1.05 (green dots), Pk5 at un = 1.22, 1.33, 1.38 (yellow dots), Pk6 at un = 1.6, 1.78 (black dots), respectively.

FIG. 20 :

 20 FIG.20: Time evolution of number of points WPk in each of the three peaks Pk1 at un = 1.0 (blue dots), Pk2 at un = 3.0 (red dots) and Pk3 at un = 5.0 (black dots), respectively, in the distribution of average speed (Fig.11). Here t denotes the sticking time Ts.

FIG. 21 :

 21 FIG. 21: Three sticky regions in the stochastic web with three-fold symmetry for ε = 3.21, β 2 = 1.75 • 10 -6 and Ẋ 0 = 3.0 as in Fig. 17 with sticking duration Ts such that log 10 (ωpeTs) = 5.63.

Fig. 20

 20 Fig.20presents these results. The time evolutions of Pk 1 , Pk 2 and Pk 3 are presented by blue, red and black dot respectively. Among all three peaks, Pk 1 is the strongest peak in the distribution. Initially, the exponents in Eq. (29) for all three peaks have very small values |1γ Pk | 1, namely γ Pk1 = 1.23, γ Pk2 = 1.17 and γ Pk3 = 1.06. For log 10 (ω pe T s ) = 5.3, the exponent value for the strongest peak Pk 1 increases to γ Pk1 = 2.11. The cross-over of the second strongest peak Pk 2 occurs for log 10 (ω pe T s ) = 4.2, when the exponent value changes to γ Pk2 = 1.64, which still implies |1γ Pk2 | < 1. For the weakest peak Pk 3 , the exponent changes to γ Pk3 = 2.15 when log 10 (ω pe T s ) = 5.7. Since the exponent value for the strongest peak changes when log 10 (ω pe T s ) = 5.3, the strength of this peak starts to decrease (leak) faster, which helps to increase the maximum ρ max of the average distribution ρ n (v) at a faster rate. Therefore, the value of the transport exponent α, from Eq. (24), increases after log 10 (ω pe T s ) = 5.3, which is also observed in Fig.14.Fig.21presents the three sticky regions for the sticking time T s such that log 10 (ω pe T s ) = 5.63. Comparing with Fig.17(with log 10 (ω pe T s ) = 2.92) extracted from the same M N time sequences, we see that the strength (number of dots) of the sticky set associated with the peak Pk 1 decreases by a very large amount and starts to become empty.

Fig. 21

 21 Fig.20presents these results. The time evolutions of Pk 1 , Pk 2 and Pk 3 are presented by blue, red and black dot respectively. Among all three peaks, Pk 1 is the strongest peak in the distribution. Initially, the exponents in Eq. (29) for all three peaks have very small values |1γ Pk | 1, namely γ Pk1 = 1.23, γ Pk2 = 1.17 and γ Pk3 = 1.06. For log 10 (ω pe T s ) = 5.3, the exponent value for the strongest peak Pk 1 increases to γ Pk1 = 2.11. The cross-over of the second strongest peak Pk 2 occurs for log 10 (ω pe T s ) = 4.2, when the exponent value changes to γ Pk2 = 1.64, which still implies |1γ Pk2 | < 1. For the weakest peak Pk 3 , the exponent changes to γ Pk3 = 2.15 when log 10 (ω pe T s ) = 5.7. Since the exponent value for the strongest peak changes when log 10 (ω pe T s ) = 5.3, the strength of this peak starts to decrease (leak) faster, which helps to increase the maximum ρ max of the average distribution ρ n (v) at a faster rate. Therefore, the value of the transport exponent α, from Eq. (24), increases after log 10 (ω pe T s ) = 5.3, which is also observed in Fig.14.Fig.21presents the three sticky regions for the sticking time T s such that log 10 (ω pe T s ) = 5.63. Comparing with Fig.17(with log 10 (ω pe T s ) = 2.92) extracted from the same M N time sequences, we see that the strength (number of dots) of the sticky set associated with the peak Pk 1 decreases by a very large amount and starts to become empty.

  Fig.20presents these results. The time evolutions of Pk 1 , Pk 2 and Pk 3 are presented by blue, red and black dot respectively. Among all three peaks, Pk 1 is the strongest peak in the distribution. Initially, the exponents in Eq. (29) for all three peaks have very small values |1γ Pk | 1, namely γ Pk1 = 1.23, γ Pk2 = 1.17 and γ Pk3 = 1.06. For log 10 (ω pe T s ) = 5.3, the exponent value for the strongest peak Pk 1 increases to γ Pk1 = 2.11. The cross-over of the second strongest peak Pk 2 occurs for log 10 (ω pe T s ) = 4.2, when the exponent value changes to γ Pk2 = 1.64, which still implies |1γ Pk2 | < 1. For the weakest peak Pk 3 , the exponent changes to γ Pk3 = 2.15 when log 10 (ω pe T s ) = 5.7. Since the exponent value for the strongest peak changes when log 10 (ω pe T s ) = 5.3, the strength of this peak starts to decrease (leak) faster, which helps to increase the maximum ρ max of the average distribution ρ n (v) at a faster rate. Therefore, the value of the transport exponent α, from Eq. (24), increases after log 10 (ω pe T s ) = 5.3, which is also observed in Fig.14.Fig.21presents the three sticky regions for the sticking time T s such that log 10 (ω pe T s ) = 5.63. Comparing with Fig.17(with log 10 (ω pe T s ) = 2.92) extracted from the same M N time sequences, we see that the strength (number of dots) of the sticky set associated with the peak Pk 1 decreases by a very large amount and starts to become empty.

  for time-Distribution of average speed for the stochastic web with ν1 = 3 at ωpe t = 8 • 10 2 , 13 • 10 3 , 21 • 10 4 , 85 • 10 4 , 34 • 10 5 , 13 • 10 6 and 5.4 • 10 7 .

		8
			log10(t) = 2.92
		7	log10(t) = 4.13 log10(t) = 5.33
			log10(t) = 5.93
		6	log10(t) = 6.53
			log10(t) = 7.14
		5	log10(t) = 7.74
	ρ ′ n (u)	4
		3
		2
		1
		0
		-4	-2	0	2	4
				u n
	FIG. 8:

  2 and v is defined below. The difference between the quantity vi and Z i is the rescaling by 1/ Distribution of average speed for the stochastic web with ν1 = 6 at ωpe t = 8•10 2 , 6.7•10 3 , 1.0•10 5 , 4.2•10 5 , 1.7•10 6 , 6.9 • 10 6 and 2.7 • 10 7 .

		0.7
			log10(t) = 2.92
			log10(t) = 3.83
		0.6	log10(t) = 5.03
			log10(t) = 5.63
		0.5	log10(t) = 6.23
			log10(t) = 6.84
			log10(t) = 7.44
	ρ ′ n (u)	0.4
		0.3
		0.2
		0.1
		0.0
		-20	-15	-10	-5	0	5	10	15
					u n
	FIG. 9:
									√ n in vi . The
									variance of the distribution of vi shrinks with rate σ 2 ∼ 1/ √ n as n increases. Since the area under the distribu-
									tion function is constant (equal to unity), the height of

Acknowledgements

This work is part of IFCPRA project 5204-3. We acknowledge the financial support from CE-FIPRA/IFCPRA. The Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. We are grateful to Professors Dominique Escande and Abhijit Sen for many fruitful discussions, and to anonymous reviewers for their comments.

VII. CONCLUSIONS

In this paper, we discuss the transport due to electrostatic waves generated by the E × B electron drift instability. The original time-dependent 3-degrees-offreedom problem is reduced to a 2-degrees-of-freedom time-dependent model and a 2-degrees-of-freedom autonomous model. Due to the wave-particle interaction, the dynamics become chaotic, and trajectories form stochastic web structures with different shape for different parameters, which we investigated for both the timedependent and time-independent descriptions.

Along with each web structure, there occur sticky islands where the trajectory spends more time compared to the purely chaotic domain, which affects the diffusion rate [32]. We use a scaling exponent for characterising the particle transport, and find that the transport is anomalous, of super-diffusive type. The presence of sticky islands generates sharp peaks in the distribution of average speed (a phase space observable) which increases the effective weight of the tail in the distribution. Considering the Poincaré recurrence theorem and Kac' lemma for the sticky sets, we estimate a power law decay for the probability that a trajectory would stick to an island. With increasing duration of the time evolution, sticky sets start to become empty and they decay with a higher exponent value. This change in the exponent γ values also affects the transport-coefficient exponent α values.

In real Hall thrusters, the E × B instability generates many unstable modes, with different frequencies and wavevectors. In this case, even for small amplitude waves, the dynamics cannot be reduced to a timeindependent 2-degrees-of-freedom model. However, each wave will typically bear its own dimensionless parameters (ε i , β i , k zi /k yi ), with small values for β i and k zi /k yi . Therefore, the several-wave dynamics will exhibit resonance overlap between the structures generated by these individual waves, resulting in smaller islands (if any survive [START_REF] Elskens | [END_REF][43][44]) and more regular transport [32]. In case of presence of harmonics modes, different types of web structures (mixed phase space) are formed in the phase space [3], which therefore will generate super-diffusive transport.

Beside the effect of several waves, three issues must also be considered. First, the thruster chamber is a cylinder, where the intensity of the radial magnetic field decreases for larger radius (here x) and the azimuthal coordinate (here y) is periodic. Second, electrons do not stay forever in the chamber, which implies that tools of transient chaos [45] will be relevant. Third, the electrons charge and current generate electromagnetic fields, so that the system needs a self-consistent many-body description.