Collisionless diffusion of particles across a magnetic field at a plasma edge in the presence of turbulence

Communications in Nonlinear Science and Numerical Simulation 13 (2008) 183–188 www.elsevier.com/locate/cnsns

Collisionless diffusion of particles across a magnetic field at a plasma edge in the presence of turbulence E. Pohn a

a,*

, M. Shoucri

b

Association EURATOM-Austrian Academy of Sciences, Kegelgasse 27, A-1030, Vienna, Austria b Institut de recherche Hydro-Que´bec (IREQ), Varennes, Qc, Canada J3X1S1 Available online 24 May 2007

Abstract We study the collisionless diffusion of electrons and ions across a magnetic field at a plasma edge, in the turbulent field produced by a Kelvin–Helmholtz instability. Ions and electrons are described by gyro-kinetic equations which include the finite Larmor radius correction and the polarization drift. Characteristic methods associated with cubic spline interpolation in two dimensions (2D) are used to integrate the gyro-kinetic equations in a slab geometry.  2007 Published by Elsevier B.V. PACS: 52.65.Ff Keywords: Vlasov equation; Eulerian codes

1. Introduction Various mechanisms can lead to charge separation and an electric field at a plasma edge. The resulting ~ E ~ B drift, which is charge and mass independent, gives rise to poloidal rotation (see [1] and the more recent work in [2]). The addition of the effect of the finite Larmor radius to this drift motion allows for a charge separation between electrons and ions to exist. Also the polarization drift has different sign for ions and electrons and accentuates a charge separation in a time varying electric field and in the presence of a strong gradient. The Kelvin–Helmholtz instability which results from this charge separation at the plasma edge can produce a low frequency turbulent spectrum. We also note that density and temperature gradients are included in the model, so other drift instabilities like the ion temperature gradient (ITG) driven instability can couple with the velocity shear instability. We study the phase–space dynamics of the particles in this turbulent field.

*

Corresponding author.

1007-5704/$ - see front matter  2007 Published by Elsevier B.V. doi:10.1016/j.cnsns.2007.05.018

184

E. Pohn, M. Shoucri / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 183–188

2. The pertinent equations We consider a two dimensional slab geometry, x being the periodic (poloidal) direction and y the non-periodic (radial) direction. The magnetic field ~ B is situated in the (x, z) plane and makes an angle h with the x axis. Both species are assumed to be homogeneous in the z direction. The gyro-kinetic equations are e  ot fi;e þ r  ð~ v? fi;e Þ þ vk cos Hox fi;e  E cos Hovk fi;e ¼ 0; ð1Þ mi;e x ~ ~ m B E? þ ð~ E? . The subscripts i and e stand for ions where ~ v? ¼ ~ vD þ~ vpi;e , ~ vD ¼ E and ~ vpi;e ¼  eBi;e2 ½ot~ vD þ~ vk Þ  r ? ~ B2 (positive sign) and electrons (negative sign). The star over a quantity is an abbreviation for an integral operator which takes into account the finite Larmor radius effect and is defined, using a Gaussian kernel G(r), by the expression Z  a ðrÞ ¼ Gðr  r0 Þaðr0 Þ dr0 : ð2Þ ~

In Fourier space each coefficient of the eik~r mode is multiplied by a factor Gk ¼ expðk 2? q2i;e =2Þ, where qi;e ¼ vti;e =xci;e is the gyroradius and k? is the component of the wave vector perpendicular to the magnetic field (k 2? ¼ k 2x sin2 h þ k 2y ). Note that the Gaussian kernel reflects the fact that the distribution function in the direction perpendicular to the R magnetic field is Maxwellian. The Poisson equation is Du ¼ 4peðni  ne Þ, ~ E ¼ ru, and ni;e ¼ fi;e dvk . The difference ni ðyÞ 6¼ ni ðyÞ is due to the finite gyroradius effect and is bigger the more important is the edge gradient. The polarisation drift ~ vpi;e accentuates this charge separation. Eq. (1) couples the drift perpendicular to the magnetic field to the kinetic effects in the direction of the magnetic field when the magnetic field is tilted with respect to the normal to the poloidal (x, y) plane. 3. Results Time is normalized to x1 pe , space to kDe and velocity to vte . We take for both species the initial density profile N(y) = 0.5(1 + tanh(1.6y)), and temperature profile Ti,e(y) = Ti,eo[0.2 + 0.4tanh(1.6y)], with Tio = Teo = 1. The initial conditions for the ions and electrons distribution function are given by m =m v2

i e k ni  fi ðx; y; vk Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2T i ðyÞ ; 2pT i ðyÞ

v2

k ne  f e ðx; y; vk Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2T e ðyÞ ; 2pT e ðyÞ

ð3Þ

pffiffiffiffiffiffiffiffi T i =T e v =x 1 where ni,e(x, y, t = 0) = N(y)(1 +  sinkox +  sin2kox +  sin3kox), kqDei ¼ Tki De ci ¼ xci =xpi ¼ 0:9 , mi/me = 1840, xci/xpi = 0.9, and  = 0.002. The computational domain is 0 6 x 6 28, 6 6 y 6 6, 5 6 vk 6 5. The number of points is Nx · Ny · Nv = 64 · 128 · 128 and the time step xpeDt = 1. The numerical scheme used is similar to what has been presented in [2,3], with the exception that the two consecutive shifts fi ðx  dx; y; vk ; tÞ and fi ðx; y  dy; vk ; tÞ are combined in one step fi ðx  dx; y  dy; vk ; tÞ, calculated by 2D interpolation using a tensor 0 0

-1 -2

-2

-3 -4

-4

-5 -6

-6

1.mode 2.mode 3.mode 4.mode

-7 -8

-8

-9 0

5000

1.mode 2.mode 3.mode 4.mode

10000 15000 20000 25000 30000

0

5000

Fig. 1. Fourier modes.

10000 15000 20000 25000 30000

E. Pohn, M. Shoucri / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 183–188

185

Fig. 2. Potential and electric field (full curve t = 0, broken curve t = 30,000).

product of cubic B-splines, preceded by an iterative step to calculate dx and dy as outlined in [4]. Fig. 1a shows, on a logarithmic scale, the time evolution of the amplitude of the first four Fourier modes of the spectrum in the poloidal x direction, obtained for h = 88.5. Only the first three modes are excited initially with the perturbation , the fourth mode is left to grow from round-off errors. A turbulent spectrum is developing, the first and second mode showing the same growth rate initially, then the second mode dominates until the end, when the first mode reaches the same level. In Fig. 1b, for h = 89, the first and second mode show the same growth rate initially, but the first mode dominates until the end. This is in agreement with the fact that for h closer to the normal, the fundamental mode should dominate (inverse cascade [5]). The simulation are carried up to xpet = 3 · 104. Fig. 2a shows (averaged over x) for the case h = 88.5 the initial profile of the potential (sinusoidal full curve) and the electric field, and the corresponding values at xpet = 3 · 104 (broken curves). The spectrum is saturating at low level, the profiles vary little. A similar plot is shown in Fig. 2b for the case h = 89, where the deformation of the initial profiles are more important. The low noise level of the Vlasov code allows to use test particles diagnostic technique to study the particles motion in the phase–space. They follow the characteristics of the Vlasov equation, dx sin h ¼ Ey þ vpx þ vk cos h dt B dy sin h ¼ Ex þ vpy dt B dvk me ¼  Ex cos h dt mi dw ¼ wr ~ v dt

ð4Þ

The last equation in Eq. (4) comes from the fact that Eq. (1) contains a source term of the form f r ~ vp , since the polarization drift is not divergence free. The quantity w(t) is a weight attached to each particle, reflecting the effect of this source term. Initially w = 1. If the polarization drift is zero, the source term f r ~ vp is zero, and w remains constant and equal to 1. We note that an invariant exists if vpy = 0, which is obtained by dividing the second and third equations in Eq. (4). We get vk 

eB y ¼ const: me;i tan h

ð5Þ

This invariant represents a straight line in the ðy; vk Þ-space, and reflects the constant generalized momentum in the homogeneous direction z. It shows that an acceleration along the magnetic field is coupled to a diffusion across the magnetic field. We show in Fig. 3a the ðy; vk Þ-space for the test electrons for h = 89 at xpet = 30,000. These electrons were initially uniformly loaded at 1 6 vk 6 1 and 2 < y < 1. They show a motion close to straight lines according to Eq. (5). The corresponding results for ions are given in Fig. 3b,

186

E. Pohn, M. Shoucri / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 183–188

Fig. 3. (v, y)-phase space for H = 89 at t = 30,000.

where the ions were initially loaded at the same position in space, and at 0:015 6 vk 6 0:015. We have calculated the mean square displacements Dy2 and Dv2k , Dy 2 ¼ Dv2k ¼

1

N part X

N part

i¼1

1

N part X

N part

i¼1

wi ðtÞ½y i ðtÞ  y i ð0Þ

2

2

wi ðtÞ½vk;i ðtÞ  vk;i ð0Þ :

Fig. 4. Diffusion of electrons (full curve H = 89, broken curve H = 88.5).

Fig. 5. Surface plots for H = 89.

E. Pohn, M. Shoucri / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 183–188

187

Fig. 6. Density (H = 89) potential (H = 88.5) charge (H = 88.5).

The time evolution of Dy2 is shown for the electrons in Fig. 4a for h = 88.5 (broken curve) and for h = 89 (full curve, showing a higher spatial diffusion). Similar results are presented in Fig. 4b for Dv2k . We show in Fig. 5a the ion density at the edge for h = 89, and a 3D view of the potential in Fig. 5b and of the charge in Fig. 5c at xpet = 30,000, dominated by the fundamental mode, while the corresponding ones in Fig. 6b and c for h = 88.5 show the signature of the first and second mode dominating at the end. Fig. 6a gives the profiles of the electrons (full curve), the ions (dotted curve), and ni ðyÞ (broken curve). 4. Conclusions A code with gyro-kinetic equations for both ions and electrons has been used to study a turbulent spectrum produced by an instability at a plasma edge related to Kelvin–Helmholtz and ITG instabilities, and the resulting collisionless diffusion of particles across the magnetic field. The integration scheme involves 2D interpolation using a tensor product of cubic B-splines to integrate the equation in (x, y) space [4]. The two cases where h = 88.5 and h = 89 show a transition from a solution dominated by a turbulent spectrum to a solution dominated by a fundamental mode. This suppression of turbulence when h gets closer to 90 is the effect of energy flowing to the longest wavelength available in the system (inverse cascade [5]), and the diffusion in space becomes more important than the diffusion in velocity space. We finally note how sensitive is the turbulent spectrum with respect to the numerical scheme. We repeat the case h = 88.5, using linear interpolation in the iterative steps to calculate dx and dy [4], keeping everything else in the code unchanged. The results in Fig. 7a show the modes having the same growth rate initially as in Fig. 1a, and saturating at the same level,

0 -1 -2 -3 -4 -5 -6 1.mode 2.mode 3.mode 4.mode

-7 -8 -9 0

5000

10000 15000 20000 25000 30000

Fig. 7. Linear interpolation in iterative steps to calculate dx and dy (H = 88.5).

188

E. Pohn, M. Shoucri / Communications in Nonlinear Science and Numerical Simulation 13 (2008) 183–188

however in the nonlinear phase the fundamental mode is continuously growing, and the higher modes declining, characteristic of numerical diffusion. In Fig. 7b we note the appearance of a big structure due to the dominant fundamental mode, which contrast with the result in Fig. 6b showing the modulation of the two dominating modes of Fig. 1a. References [1] [2] [3] [4] [5]

Burrell KH. Phys Plasmas 1977;4:1499. Shoucri M, Manfredi G, Bertrand P, et al. J Plasma Phys 1999;61:191. Shoucri M, Gerhauser H, Finken K-H. Comp Phys Commun 2003;154:65. Shoucri M. Czech J Phys 2001;51:1139. Knorr G, Pe´cseli HL. J Plasma Phys 1989;41:157.