Proceedings of ITC/ISHW2007

ICNTS - Impact of Incompressible E × B Flow in Estimating
Mono-Energetic Transport Coefficients
Craig D. BEIDLER1) , Maxim Yu. ISAEV1) , Sergei V. KASILOV3,4) , Winfried KERNBICHLER3) ,Henning MAASSBERG1) , Sadayoshi MURAKAMI5) , Victor V. NEMOV3,4) , Donald A. SPONG6) andVictor TRIBALDOS7)
1)

Max-Planck Institut für Plasmaphysik, Assoc. EURATOM, Greifswald, Germany
2)
Russian Research Center ”Kurchatov Institute”, Moscow, Russia
3)
Institut für Theoretische Physik, Assoc. EURATOM-ÖAW, TU-Graz, Austria
4)Institute of Plasma Physics, NSC-KhIPT, Kharkov, Ukraine
5)
Department of Nuclear Engineering, Kyoto University, Kyoto, Japan
6)
Oak Ridge National Laboratory, U.S. Department of Energy, USA
7)
Laboratorio Nacional de Fusion, CIEMAT, Assoc. EURATOM-CIEMAT, Madrid, Spain
October 11, 2007
In the neoclassical ordering, an incompressible E × B drift is assumed in the 1 st order drift-kinetic equation(i.e. the approximation E × B/hB2 i is used). Then, the terms with (E × B) · ∇B disappear in both equations for ṗand v̇ allowing for a monoenergetic treatment (with a conservative formulation of the equations of motions), as
it is used e.g. in DKES. As a consequence, the magnetic moment is not an invariant in this approximation. For
large radial electric fields, however, this ordering scheme is violated since the (E × B) · ∇B term in ṗ becomescomparable with the mirror term (∝ B · ∇B). The impact of this simplification is tested by strictly local δ-fMonte Carlo techniques (i.e. with ṙ = 0 for the equation of motion in 1 st order) where the E × B drift is treatedin compressible and incompressible form (both completely conservative) for large radial electric fields.
Keywords: ICNTS, neoclassical transport, stellarators, radial electric field

1

Introduction

The International Collaboration on Neoclassical Transport in Stellarators (ICNTS) was initiated in 2000. Thereis a subtle inconsistency in the way the effect of large radial
electric fields are handled in obtaining neoclassical estimations of transport coefficients. Rather than being a limitation of the theory itself, the complication arises from some
simplifying approximations used. The electric field enters
the calculation through the E×B drift, VE = E×B/B2 , andthe small gyroradius ordering or drift ordering assumes
that the ratio of this drift to the thermal speed of particles
is of 1 st order, i.e. VE /vth ∼ ρ/L ∼ ωth /Ω = δ << 1, whereρ is the thermal gyroradius, L a characteristic plasma scale
length and ωth and Ω the transit frequency and gyrofrequency respectively. This approximation implies that theradial electric field is not so large as to distort gyrationand varies slowly in time. The drift ordering is not very
stringent, and is easily fulfilled even at the low-order electric field resonances (where the poloidal component of thedrift speed vanishes) E res = ιvBr/R since VE res /vth = ιr/R,
where r/R is the inverse aspect ratio of the considered flux
surface and ι its rotational transform.
In the drift ordering approximation particle motion isaveraged over the the small gyro-scale and the resulting kiauthor’s e-mail: V.Tribaldos@ciemat.es

netic equation, now describing the distribution of guiding
centers, is known as the drift kinetic equation (DKE) [1, 2].
Unfortunately, even with the drift ordering approximation
solving the DKE is a daunting task. The inhomogeneity ofthe magnetic field along with the non-linearity of the collisional term makes the DKE a non-linear partial differential equation in a six dimensional space f = f (r, v, p, t),r being the guiding center position in 3D space, v andp = v·B/v its speed and pitch angle in velocity-space, and t
the time. The central point of neoclassical theory is solving
this equation to obtain different flux-surface-averaged moments of its solution. To this end the distribution of guiding
centers is linearized around a local Maxwellian distribution
and several approximations are made in the drift motion of
guiding centers. In the drift ordering approximation, thetotal energy conservation, E = mv2 /2 + qΦ with E = −∇Φ,translates into an equation for the variation of the kinetic
energy which is proportional to the divergence of the E × B
drift, v̇ = v(1 + p2 )/4 ∇ · VE (notice that drift ordering precludes rapid variation of the magnetic and electric fields.
Whenever the E × B drift is incompressible, ∇ · VE = 0, thetotal energy and the kinetic energy are conserved, whichallows for quite some simplification in solving the DKE.
Not only the speed of particles can be considered as a parameter, allowing for a mono-energetic treatment, but also
the collision operator can be approximated by only its pitchProceedings of ITC/ISHW2007
angle scattering part. However, a direct calculation showsthat ∇ · VE = −2VE ∇ ln B which in general is different
from zero, thus to retain the benefits of incompressibility
usually the exact drift VE = E × B/B2 is approximated by
VE = E × B/hB2 i. Since the magnetic field gradient is determined by the device’s configuration the incompressible
approximation restricts the usual neoclassical treatment tosmall radial electric fields.
The purpose of this work is to shed some light on theeffect of the incompressible E × B flow for a wide rangeof radial electric fields and its impact on the determination of the mono-energetic particle diffusion coefficient,D∗ . In Sec. 2 the basic equations of motion are discussed
and a local δ f Monte Carlo [3, 4] (MC) technique capable of dealing with compressible and incompressible flows(both completely conservative) is presented. Results ofthe comparison between the DKES (Drift Kinetic Equation
Solver) [5, 6] and the new δ f MC codes for different configurations will be shown is Sec. 3. Finally, Sec. 4 contains
some discussion of the results.

2

Basics

The starting point of neoclassical transport theory is the
DKE, that can be formally written as [2]df∂f∂f= ṙD · ∇ f + ṗ+ v̇= C( f, f )
dt
∂p
∂v

(1)

where ṙD = vD , ṗ and v̇ are the drift speed and thetime derivatives of the pitch (related to the conservationof the magnetic moment µ = mv2⊥ /2B) and the kinetic
energy (derived from the conservation of the total energy
E = mv2 /2 + qΦ) given by:
vD = pv
ṗ = −v̇ = −

B
mv2(1 + p2 ) B × ∇B
+ VE +
B
2qB3

(2)

p
v
(1 − p2 ) B · ∇B −(1 − p2 ) VE · ∇B(3)
2
2B
2B

v
(1 + p2 ) VE · ∇B
2B

(4)and C( f, f ) is a collision operator.
Since VE · ∇B = −B/2 ∇ · VE the second term in the
r.h.s. in Eq.3 and Eq.4 depend explicitly on the compressibility of the E × B flow. The procedure used to solve the
DKE consists of linearising the distribution of guiding centers, f , in Eq. 1 with respect to the drift ordering small
parameter δ = ρ/L << 1 as f = f0 − δ f1 (∂ f0 /∂rr ). For stationary conditions, i.e. neglecting the explicit time dependence, ∂ f /∂t, the solution to the zero order, δ0 , equation
is identically satisfied by the local Maxwellian f0 = f M .
Therefore, the goal is to find the solution, f1 , to the firstorder DKE:vDs · ∇rs f1 + ṗ

∂ f1
∂ f1
+ v̇
− C( f M , f1 ) = vrD
∂p
∂v

(5)

Notice that the full first order distribution depends on theradial gradients of the zero order Maxwellian distribution,and that the inhomogeneous term is the radial drift of guiding centers vrD , thus the radial dependence enters the equation just like a parameter. Therefore, Eq. 5 describes a
diffusion process in phase space rather than in real space.
Different methods are usually applied to solve the Eq. 5,
each with its strengths and drawbacks: i) analytical calculations [1]; ii) explicit spectral procedures, expanding f1 ina base of eigenfunctions (like the DKES code [5, 6]) andiii) Monte Carlo techniques [7]. Analytical solutions areusually restricted to simplified magnetic fields and collisionality regimes, on the other hand explicit and MC numerical methods can deal with realistic configurations forbroad collisionality ranges, but are extremely time consuming and sometimes not very accurate (e.g. non diagonal transport matrix elements in MC.
The usual approximation, as for example is done in
DKES code, for computing the diffusion transport coefficients consists in neglecting ∇ · VE in equations 3 and 4,
i.e. considering the E×B flow incompressible. Since v̇ = 0,
the energy in Eq. 5 enters only as a parameter, once the collision operator is approximated by just its pitch-angle scattering part, and the coefficients are obtained by the convolution of mono-energetic solutions with the Maxwellian
distribution. The price to pay in this approximation is thatthe radial electric field cannot be very large, and that themagnetic moment µ is not conserved (see Eq. 3). To checkthe approximation made by DKES for the mono-energetic
diffusion coefficient for large radial electric fields it would
be desirable to benchmark it against a MC type calculation with and without assuming ∇ · VE = 0. In usual fullf [7] MC calculations, based on fitting the slope of the timedependence of the radial broadening of a test particle ensemble (diffusion in real space), it is easy to include thefull set of equations (conserving µ and E. However, when
the radial electric field becomes large particles can win orlose kinetic energy as they move radially because of theirdrifts, thus making the calculation non-mono-energetic.
The method proposed here to make such comparison is using the method of characteristic to solve partial differential
equations like Eq. 5, which is at the base of the δ f MC
technique.
Formally, the solution to an equation of the typea(x, y, ...)∂ f /∂x + b(x, y, ...)∂ f /∂y + ... = g(x, y, ..., f ) is asurface S such that at each point (x, y, ...) on S , the vector V = (a(x, y, ..), b(x, y, ...), ..., g(x, y, .., f )) lies in the tangent plane. Such surface can be constructed by the unionof curves C parametrized by s such that at each pointon the curve, the vector V is tangent to the curve. In
particular C = {(x(s), y(s), ..., g(s))} will satisfy the following system of ordinary differential equations:dx/ds =a(x, y, ...); dy/ds = b(x, y, ...); ...; du/ds = g(x, y, ..., f ),called characteristic curves, and the solution is f (x, y, ...) =u(x, y, ...). Rephrasing this method for the collisionlessProceedings of ITC/ISHW2007

dp
dv
du
= vDs ;= ṗ ;= v̇ ;= vrD (6)
dt
dt
dt
dt
f1 = u being its solution. Please notice that even though
the system of equations 6 are the equations of motion of
the guiding center, and the characteristic curves are the trajectories, thus clarifying the link with a MC particle simulation view, there is no equation for the radial drift rDr .
Particles cannot escape the birth flux surface and no radialbroadeningR calculation can be done. However, the solution
f1 = u ∼ vrD dt directly depends on the bounce-averageof the radial drift speed along the parallel motion. The MC
implementation is straightforward; the system of equations
is integrated (taking care to retain the conservative natureof the motion) for an ensemble of markers (since no radial motion is included they cannot be regarded as particles
anymore), all sharing the same energy and flux surface butwith random toroidal and poloidal position and pitch. The
effect of collisions in the DKE is simulated by applying apitch angle collision operator [7].

3Results

Three different magnetic configurations were studied withquite different B structures, namely: a tokamak with thesame aspect ratio and rotational transform as the W7-AS
stellarator, the LHD heliotron configuration with major radius Raxis = 3.75 m, and the W7-X standard configuration.
The impact of the incompressibility of the E × B flow wasstudied by evaluating the diffusion coefficient for a wide
range of radial electric fields. For each case two different collisionalities were chosen corresponding to the beginning of the ripple regime, ν/v = 10−3 m−1 , and one at
the beginning of the lmfp, ν/v = 10−4 m−1 . The δf MC integrates 1024 markers for three collisional times divided ineight groups of 128 markers each. The error bar is obtainedwith the standard error from the eight estimations.
As was noted long ago, e.g. [8] Eq. 5 has a singularity
when the poloidal component of the drift speed vanishes,
which corresponds to resonant radial electric field values.
In the following the radial electric field has been normalized to the first toroidal resonance Eres = ιvBr/R.
In Fig. 1 the diffusion coefficient at half radius, r/a =0.5, for an ideal tokamak configuration, with B(r, θ) =B0 (1 + r/R cos θ), R = 2m, a = 0.2m, and ι = 0.51. Thiswas selected because radial excursions from the flux surfaces are small and only one electric field resonance exists.
There is a good agreement between DKES code results andthe δf MC incompressible calculation; both displaying apeaked feature around the resonance followed by a strong
decrease because of the disappearance of banana orbits.
When compressible effects are included in the calculation
the resonance is smoothed as well as the later sharp diffusion decrease flattened. The reason being the variation in
the kinetic energy of the particles, v̇ in Eq. 4, as can be seenin Fig. 2. The broadening of the kinetic energy spectra, due
to v̇, reduces the number of particles at the resonant field,
Er /Eres = 1 but also pushes some particles to energies that
resonate at larger values of Er /Eres . At very large radial
electric fields there is a systematic difference between the
DKES and MC results which is attributed to numerical diffusion. Nevertheless, please notice the rather small valuesof the diffusion coefficient.
The results for the other two devices (see Figs. 3 and4) are similar to the tokamak result, apart from the different resonance structure due to their broader magnetic field
spectra. The helical resonance peak at Eres = (N periods −mι)vBr/R can be clearly identified. As in the tokamak casethe flattening is more pronounced at smaller collisionalities, pitch angle scattering is less efficient in moving particles out of the resonance. The incompressible approachis rather good for small Er , but unexpectedly it is also
working at large radial electric fields, provided that there is
no resonance overlapping; see the almost mono-energetic
structure of the PDF in Fig. 2 at Er /Eres = 4. Finally, there
is an inconsistency in comparing the mono-energetic diffusion coefficients of the DKES and incompressible δf codeswith the compressible result since the latter, as has beenshown, is not really mono-energetic. Moreover, considering only the pitch-angle part of the collision operator close
to the resonance is not justified [8, 7].

D*

drDs

10

1

10

0

!/v=1e-3

10

-1

10

-2

10

-3

10

-4

DKES3
Incompressible
Compressible

-5

10
1
10
10

D*

DKE means solving the system of equations:!/v=1e-4

DKES3
Incompressible
Compressible

0

10

-1

10

-2

10

-3

10

-4

10

-5

10

-1

0

10
E /E
r

10

1

res

Fig. 1 Diffusion coefficient, normalized to the tokamak plateau,
for the tokamak configuration versus the radial electricfield, normalized to the first toroidal resonance Eres , forν/v = 10−3 m−1 (top) and ν/v = 10−4 m−1 (bottom) obtained with the DKES code (circles) and the δf MC incompressible (squares) and compressible (triangles) approaches.

Proceedings of ITC/ISHW2007
2

10

0.04
1.33
1.65
4

1

10

0

-1

10

D*

PDF

10

-2

10

-3

10

-4

10

1

10

0

10

-1

10

-2

10

-3

10

-4

10

-5

10

-6

!/v=1e-3

DKES3
Incompressible
Compressible
-7

-1

0

1

(v-v )/v
0

2

3

10
1
10

4

0

10

D*

Fig. 2 Probability distribution function of velocities for ν/v =10−3 m−1 of Fig. 1 in the compressible approach.

Discussion[1] AA Galeev and RZ Sagdeev, Sov. Phys. JETP 26, 233 (1968)[2] RD Hazeltine, Phys. Fluids 15, 77 (1973)[3] M Sasinowsky and A Boozer, Phys. Plasmas 2, 610 (1995)[4] AM Dimits and WW Lee, J. Comput. Phys. 107, 309 (1993)[5] SP Hirshman et al., Phys. Fluids 29, 2951 (1986)[6] WI van Rij and SP Hirshman, Phys. Fluids B 1, 563 (1989)[7] A Boozer and G Kuo Petravic, Phys. Fluids 24, 851 (1981)[8] EF Jaeger, CL Hedrik and DA Spong, Nuc. Fusion 19 Vol.
12, 1627 (1979)

10

-1

10

-2

10

-3

10

-4

10

-5

10

-6

10

-7

10

DKES3
Incompressible
Compressible
-1

10

0

10

E /E
r

1

res

Fig. 3 Diffusion coefficient, normalized to the tokamak plateau,for LHD configuration with Raxis = 3.75 m versus theradial electric field, normalized to the first toroidal resonance , and ν/v = 10−3 m−1 (top) and 10−4 m−1 (bottom)
for the DKES code (circles) and the δf MC incompressible (squares) and compressible (triangles) approaches.

100!/v=1e-3

-1

10

10-2
D*

The aim of this work was to check the accuracy of local
mono-energetic transport coefficient calculations for large
radial electric fields. On the one hand DKES estimations
are mono-energetic and local and transport is due to apitch-angle diffusion process in phase space. But, the E×Bterm is treated as incompressible, ∇ · VE = 0, to keep the
kinetic energy constant. On the other hand, common full- f
MC based calculations are not exactly local, because radial
space broadening is used to describe the diffusion process.
And, again their applicability is limited to consider small
radial electric fields to keep particles’ kinetic energy almost constant during their radial drifts. Here, a fully conservative δf MC code was introduced where local diffusion is computed in phase-space with particles (markers)
remaining indefinitely on their birth flux surface (solving ageneral MC problem with particle losses).
The impact of the incompressible E × B flow approximation was tested comparing DKES results with thestrictly local δf MC where the E × B drift is treated in compressible and incompressible form for large radial electricfields and several magnetic field configurations.
It was found that considering ∇ · VE = 0 is indeed
a good approximation for small radial electric fields, as
expected, but also for large Er far from the resonances.
This also serves as a benchmark between DKES and the δf
MC. The discrepancy between the results with and without∇ · VE = 0 close to the resonance is due to the variation of
the kinetic energy. This results calls into question the usual
mono-energetic calculations, which, to be safe, should be
restricted to Er < 0.5 − 0.7 Eres .

!/v=1e-4

0

10-3
10-4

DKES3
Incompressible
Compressible

10-5
100
!/v=1e-4
10-1
10-2
D*

4

10

10-3
10-4
10-5
10-1

DKES3
Incompressible
Compressible
100

101
E /E
r

res

Fig. 4 Diffusion coefficient , normalized to the tokamak plateau,for W7-X standard configuration versus the radial electric field, normalized to the first toroidal resonance, andν/v = 10−3 m−1 (top) and 10−4 m−1 (bottom) for the DKES
code (circles) and the δf MC incompressible (squares)
and compressible (triangles) approaches.