Proceedings of ITC/ISHW2007

Monte-Carlo Simulation of Neoclassical Transport in Magnetic
Islands and Ergodic Region
Ryutaro KANNO1,2) , Masanori NUNAMI3) , Shinsuke SATAKE1) , Hisanori TAKAMARU4) ,Yukihiro TOMITA1,2) , Keisuke NAKAJIMA2) , Masao OKAMOTO4) and Nobuyoshi OHYABU1,2)
1)2)

National Institute for Fusion Science, Toki 509-5292, Japan
Department of Fusion Science, Graduate University for Advanced Studies, Toki 509-5292, Japan
3)
Institute of Laser Engineering, Osaka University, Suita 565-0871, Japan
4)
Department of Computer Science, Chubu University, Kasugai 487-8501, Japan

It is shown in Large Helical Device experiments that the transport modeling based only on the fluid
description is not sufficient for expressing edge transport phenomena in a magnetic island. Existence of abootstrap current around the island is strongly suggested. On the other hand, in recent tokamak experiments
it is found that so-called stochastic diffusion theory based on the “field line diffusion” over estimates theradial energy transport in the collisionless edge plasma affected by resonant magnetic perturbations, though theperturbations induce a chaotic behavior in the field lines. These results imply that the conventional modelingof the edge transport should be reconsidered for a lower-collisionality case. In order to take a new look at
the modeling of the edge transport, we investigate neoclassical effect on the transport in magnetic islands
and ergodic region. By using a drift kinetic equation solver without an assumption of existence of nested
flux surfaces (the KEATS code), it is possible for us to execute the investigation. The simulation results
show that the radial energy flux of ions for a lower-collisionality case is quite small compared to the prediction of the stochastic diffusion theory, while the flux for a high-collisionality case is consistent with the prediction.
Keywords: neoclassical transport, edge plasma, magnetic island, ergodic region, Monte-Carlo simulation

1

Introduction

It is shown in Large Helical Device (LHD) experiments that the transport modeling based only on the fluid
description is not sufficient for expressing edge transportphenomena in a magnetic island [1, 2, 3]. This result isgiven in the experiments of observing the healing of the
m/n = 1/1 magnetic island in the edge, where m and n are
the poloidal and toroidal mode numbers, respectively. A
current depending on the pressure gradient (i.e. the bootstrap or Pfirsch-Schlüter current) is expected to explain thehealing in the experiments. In results of a simulation study
based on the fluid description [4], the healing phenomenon
is not explained by the Pfirsch-Schlüter current only. The
important role of a bootstrap current in the edge region isstrongly suggested, and thus the kinetic modeling of theedge plasma is needed for understanding of the edge transport phenomena, where in the LHD experiments the temperature is &500 eV and the plasma density ∼ 1019 m−3 inthe island.
On the other hand, in recent tokamak experiments it
is found that so-called stochastic diffusion theory basedon the “field line diffusion” [5, 6] over estimates the radial energy transport in the edge region added resonant
magnetic perturbations (RMPs) [7, 8]. This fact is discovered in the experiments of ELMs (edge localized modes)
author’s e-mail: kanno@nifs.ac.jp

suppression by means of RMPs in collisionless tokamak
plasmas [7, 8]. (Historically, the idea of suppressing the
ELMs and controlling the edge transport by using RMPs
has been proposed about 20 years ago [9].) When the
RMPs induce a chaotic behavior in the field lines, thestochastic diffusion theory predicts that a thermal diffusivity is given as χaql = χak |δBr /Bt |2 or vath πRax q|δBr /Bt |2 forthe collisional or collisionless limit, where a is a particle
species, χak = 3.91T a τa /ma the parallel diffusivity, T a thetemperature, τa the collision time, ma the particle mass,
vath the thermal velocity, Rax the major radius of the magnetic axis, q the safety factor, δBr the strength of RMPs,
and Bt the toroidal component of the magnetic field. Thisprediction has been demonstrated in experiments on highcollisional tokamak plasmas [10]. However, in collisionless plasmas, the experimental thermal diffusivity χex is inconsistent with the prediction of the stochastic diffusion
theory, e.g. χeql /χeex  10 for the electron thermal diffusivity [8]. Small RMPs cause the complete suppression of
the ELM events, and have a negligible effect on the energy
confinement.
The above experimental results in torus plasmas imply
that the conventional modeling of the edge transport should
be reconsidered for a lower-collisionality case. There isno established theory describing radial transport in magnetic islands and ergodic region. In order to take a newlook at the modeling of the edge transport, we investigate

Proceedings of ITC/ISHW2007
neoclassical effect on the transport in magnetic islands andergodic region. Here, even in the three dimensional field
line structure disturbed by RMPs, the Coulomb collisioncauses the transition between a passing particle orbit and atrapped particle orbit in toroidal and helical ripples (localized and/or blocked particle orbits) [11]; in the present paper we call it the neoclassical effect on the edge transportphenomena. Recently, we develop a new transport simulation code without an assumption of existence of nested
flux surfaces; the code is named “KEATS” [12, 13]. The
code is programmed by expanding a well-known MonteCarlo particle simulation scheme based on the δ f method[14, 15]. By using the KEATS code, it is possible for us to
execute the investigation. In this paper we show the simulation results, applying the code to a torus plasma having
the ergodic region in the edge.

2

In general, effects of neutrals and an electric field are important in the edge transport, but in the present paper theseeffects are neglected for simplicity. (The modeling of afluctuating field in the KEATS code is described in Refs.
[13, 17].)
To solve Eq. (1) by Monte-Carlo techniques, we adopt
the two-weight scheme of the δ f formulation [14, 15]. In
evolution of the δ f part, the background fM is assumed tobe fixed because the background is in a quasi steady-statefrom the viewpoint of the δ f part. The Monte-Carlo simulation code, KEATS, is programmed in an Eulerian coordinate system, i.e. so-called helical coordinates [4], thus
the code does not need magnetic flux coordinates. Simulation results (e.g. estimation of particle and energy fluxes)of the KEATS code for a case of a simple tokamak field areagreed with ones of the “FORTEC-3D” code [15] whichuses magnetic flux coordinates.

Simulation Model

We consider that a guiding center distribution function of plasma f is separated into an equilibrium-likebackground f0 and a kinetic part δ f of the distribution,f = f0 + δ f , where the kinetic part δ f is considered asa small perturbation from f0 . The zeroth-order distribution function f0 is given as a local Maxwellian distribution
f0 = fM (x, ξ, v) = n{m/(2πT )}3/2 exp{−mv2 /(2T )}, whereξ = vk /v is the cosine of the pitch angle, vk = v · b,b = B/B the unit vector along a field line, B a magnetic
field, B = |B|, v = |v|, n = n(x) the density, m the particle
mass, and T = T (x) the temperature. Applying the decomposition f = fM + δ f to the drift kinetic equation, we havethe following equation of the kinetic part δ f :D
δ f = − {vd · ∇ fM − CF fM } ,Dt(1)

where the operator D/Dt is defined as D/Dt := ∂/∂t +(vk + vd ) · ∇ − C T , vk = vk b the parallel velocity, and vd thedrift velocity of guiding center motion. The test particle
collision operator C T is given, for simplicity, as
CT =

 ∂
νdef ∂ 
1 − ξ2
,
2 ∂ξ
∂ξ

(2)

m
v · u0 ,T(3)and it can be implemented numerically by random kicks invelocity space [16], which represents the Coulomb scattering process, where νdef is the deflection frequency. Here,
we should note statistical accuracy of the operator C T expressed by the Monte-Carlo method [16], in particular,
around |ξ| ≈ 1. The operator C F is the field particle collision term, which represents local momentum conservation(CF is needed to treat accurately the parallel transport):
C F = νdef

and u0 is given as
,Z
Z
mv2
3
u0 =
d v νdef v δ f
fM .
d3 v νdef
3T

(4)

3

Simulation Results

For the investigation of neoclassical effect on the transport in the ergodic region, we use a magnetic configuration
which is formed by adding RMPs into a simple tokamak
field having concentric circular flux surfaces, where themajor radius of the magnetic axis Rax = 3.6 m, the minor radius of the plasma a = 1 m, and the magnetic fieldstrength on the axis Bax = 4 T. The Poincaré plots of themagnetic field lines on a poloidal cross section are shownin Fig.1. Onepcan see the ergodic region in r/a = 0.7 ∼ 1,
where r = (R − Rax )2 + Z 2 . In the KEATS code, thenumber of test particles is NTP = 16, 000, 000.
To investigate effect of the existence of the ergodic
region on the transport phenomena, we evaluate the ion
energy flux Qi in two cases, i.e., in the configurations (a)having lower edge temperature T edge ∼ 200 eV at a centerof the ergodic region and (b) having higher edge temperature T edge ∼ 1 keV. The temperature profile is given as
T i = T ax {0.02 + 0.98 exp[−4(r/a)α ]} with T ax = 2 keV andα = 2.5 (case (a)) or 7.86 (case (b)), which neglects theexistence of the ergodic region. The density profile is set
homogeneous, ni = const. = 1×1019 m−3 . The backgroundfM is fixed in the calculations. The radial profiles of the energy flux estimated from the KEATS computations for twocases (a) and (b) are shown in Figs.2a and b, respectively.
For simplicity, the radial energy fluxes are given neglecting the existence of the ergodic region, because we haveno magnetic coordinate system including several magnetic
field structures as the core and ergodic regions. The energy
flux Qi is averaged over concentric circular shell region in
the whole toroidal angles as if there were nested flux surfaces. Here, in the KEATS computations the energy flux Q
is given as
Q(t, x) =

Z

d3 v

mv2
(vk + vd )δ f ,2(5)

Proceedings of ITC/ISHW2007

Fig. 1 Poincaré plots of magnetic field lines on a poloidal crosssection.

where · · · means the time-average, and the averaging time
is longer than the typical time-scale of δ f . The heat flux
predicted by the stochastic diffusion theory is given as
qql = nχql ∇T . From the results of the KEATS code, wefind that the ion energy flux is affected by the RMPs, but
the flux in the lower-collisionality region around r/a ≈ 0.8
is quite small compared with the prediction of the stochastic diffusion theory, shown in both cases (a) and (b). On
the other hand, in the high-collisionality region aroundr/a ≈ 1, the flux is consistent with the prediction. These
results can be explained theoretically, see Appendix.

4

Summary

We have been developing the neoclassical transport
code, KEATS, to study the transport phenomena in theislands and ergodic region. We apply the code to theedge disturbed by resonant magnetic perturbations, andfind that the ion energy flux estimated by the KEATS
code for a lower-collisionality case is quite small compared to the prediction of the stochastic diffusion theory
based on the “field line diffusion,” while the flux for a highcollisionality case is consistent with the prediction.

Acknowledgements
This work was performed with the support andunder the auspices of the NIFS Collaborative Research Program NIFS06KLDD007, NIFS06KTAT018 and
NIFS07KTAD001, and was also supported by MEXT.
HAITEKU(2004) and the 2006-2007 JIFT exchange program.

Fig. 2 Radial profiles of ion energy/heat flux for (a) lowerand (b) higher edge temperatures, where r =p(R − Rax )2 + Z 2 . The collisionless (or collisional) quasilinear (QL) model means the prediction of the stochastic
diffusion theory for the collisionless (or collisional) case.

Appendix. Stochastic Analysis of Radial Transport in a Perturbed Field
A fluid equation in a steady-state corresponds toa stochastic differential equation described as dXti =jγU i (X t )dt + cij (Xt )dWt and i, j = 1, 2, 3 [18, 19], whereγ is a constant (e.g. γ = 5n/2 for the heat balance equation if n = const.), U = (U 1 , U 2 , U 3 ) a flow in a steadystate, Di j = cik gk` c`j a diffusion coefficient, gk` a metric coefficient, X t = (Xt1 , Xt2 , Xt3 ) a diffusion process, andW t = (Wt1 , Wt2 , Wt3 ) a Brownian process. It is assumed thata fluid is exposed to noise caused by resonant magnetic
perturbations (RMPs), and that a fluid particle motion is
described by an Itô process dYti = γŨ i (t, ω)dt + cij (Y t )dWtj
instead of the process X t , where the flow is represented
as Ũ(t, ω) = U(Y t ) + “noise”, ω is a label of a fluid
particle, “noise” is a random function having zero mean
and finite strength, and the definition of an Itô process isgiven in Ref. [20]. It is known that an Itô process Y t coincides in law with a diffusion process X t if and only if
E x0 [Ũ(t, ω)|PYt ] = U(Y t ) [20], where X 0 = Y 0 = x0 is a
starting point of a fluid particle at t = 0, PYt is the σ-algebra
generated by the set {Y s ; 0 ≤ s ≤ t}, and E x0 [· · · |PYt ] denotes the conditional expectation with respect to PYt . This

Proceedings of ITC/ISHW2007
theorem means that the “noise” cannot cause the diffusionof fluid particles. We should reconsider the reason why thenoise created by RMPs can affect the transport.
Let us take the following collision operator:
)
(
∂
2 ∂f
,
(A.1)
· u f + vth
C( f ) = νcol∂u∂u

where νcol = νcol (x) is the collision frequency, vth the thermal velocity, and v = U+u the velocity of a guiding center,and U = U(x) the mean velocity [21]. The operator (A.1)
is simpler, but is used only to get a rough idea of collisional effects [22]. We consider the motion of a guiding
center along a field line for estimation of radial spreading
the guiding centers by their parallel motions in a perturbed
magnetic field. The guiding center motion exposed to thecollisions (A.1) is given as an Ornstein-Uhlenbeck process:dx = vdt = (U + u)dt,
du = −νcol udt + σdW kt ,

(A.2)
(A.3)
√

where U = Uk b, u = uk b, σ = vth νcol , W kt a Brownian
process for the parallel direction, i.e. dW kt = bdWt , andb = B/B the unit vector along a field line. Here, the effectsof ripples are neglected for simplicity. The solutions ofEqs. (A.2) and (A.3) are described respectively asZ t
x = x0 +(U + u)ds,
(A.4)
0
Z t
u = e−νcol t u0 +e−νcol (t−s) σdW ks ,(A.5)
0where x0 and u0 are the initial values at t = 0.
One may consider that effect of a perturbation fieldon the motion is interpreted as noise on the motion. If theeffect is expressed by a linear operator ṽ = Ñv, thendx = (v + ṽ)dt = (v + Ñv)dt.

(A.6)

The solution of Eq. (A.6) is given asZ t
x = x0 +(1 + Ñ) U + e−νcol s u0 ds
0
Z t Z s+ds
e−νcol (s−h) σ(1 + Ñ)dW kh . (A.7)
00For the collisional limit t  1/νcol , the diffusion (causedby the perturbation field) in configuration space is derived
from Eq. (A.7):
vth
(A.8)
dx ≈ (1 + Ñ)Udt + √ (1 + Ñ)dW kt ,
νcoli.e., for the collisional limit the diffusion in velocity space
directly becomes the diffusion in configuration space. We
should note that the diffusion in configuration space originates from the collisions in velocity space. When the
RMPs are added to the original magnetic field having
nested flux surfaces, the parallel motion of a guiding center
may cause radial fluctuation in configuration space [5, 6].
If the noise ṽ = (ṽ1 , ṽ2 , ṽ3 ) is given asṽi = (Ñv)i =δBr εi jk√ θ̂ j vk φ̃(t, i),
Bt
g

(A.9)

then the radial diffusivity Dr = Dk |δBr /Bt |2 is obtained in
the fluid equations given from the drift kinetic equation
having the collisions (A.1) for the collisional limit, where θ̂
is the unit vector for the poloidal direction, δBr the strengthof the RMPs satisfying |δB r/Bt |  1, Bt the toroidal component of B, g = det gi j the square of Jacobian, εi jk
the Levi-Civita symbol, φ̃(t, i) the ith component ofh aizeromean random vector having the mean square of E φ̃2 = 1
and being independent of dW kt , and Dk = v2th /νcol the parallel diffusivity. We should note that the noise term Ñ(U+u0 )cannot cause diffusion, as shown in the first paragraph inthis section; see also Refs. [13, 17, 20].
The above discussion shows that for a high-collisional
plasma (the characteristic time t  1/νcol ), the guiding center motions become close to the prediction of thestochastic diffusion theory based on the “field line diffusion” [5, 6]. On the other hand, for a lower-collisionality
plasma (t . 1/νcol ), the motions are not interpreted as the
diffusion process predicted by the stochastic diffusion theory. These consequences are consistent with the results
shown in Fig. 2 and also ones obtained in the test particle
simulations [23].
[1] N. Ohyabu et al., Phys. Rev. Lett. 88, 055005 (2002).
[2] N. Ohyabu et al., Plasma Phys. Control. Fusion 47, 1431
(2005).
[3] Y. Nagayama et al., Nucl. Fusion 45, 888 (2005).
[4] R. Kanno et al., Nucl. Fusion 45, 588 (2005).
[5] A.B. Rechester and M.N. Rosenbluth, Phys. Rev. Lett. 40,
38 (1978).
[6] B.B. Kadomtsev and O.P. Pogutse, in Proc. 7th IAEA Int.
Conf. on Plasma Physics and Controlled Nuclear Fusion
Research (IAEA, Vienna, 1979) Vol.1, pp.649.
[7] K.H. Burrell et al., Plasma Phys. Control. Fusion 47, B37(2005).
[8] T.E. Evans et al., Nature Phys. 2, 419 (2006).
[9] N. Ohyabu et al., Nucl. Fusion 27, 2171 (1987).
[10] Ph. Ghendrih et al., Plasma Phys. Control. Fusion 38, 1653
(1996).
[11] S. Jimbo et al., Nucl. Fusion 45, 1534 (2005).
[12] M. Nunami et al., Research Report NIFS Series No. NIFS871 (2007).
[13] R. Kanno et al., accepted for publication in Contributions
to Plasma Physics.
[14] W.X. Wang et al., Plasma Phys. Control. Fusion 41, 1091
(1999).
[15] S. Satake et al., Plasma Fusion Res. 1, 002 (2006).
[16] A.H. Boozer and G. Kuo-Petravic, Phys. Fluids B 24, 851(1981).
[17] R. Kanno et al., Plasma Fusion Res. 1, 012 (2006).
[18] Y. Feng et al., J. Nucl. Mater. 266-269, 812 (1999).
[19] A. Runov et al., Nucl. Fusion 44, S74 (2004).
[20] B. Øksendal, Stochastic Differential Equations (SpringerVerlag, Berlin Heidelberg, 2003).
[21] P.C. Clemmow and J.P. Dougherty, Electrodynamics ofParticles and Plasmas (Addison-Wesley, Reading, Mass.,1969).
[22] D.R. Nicholson, Introduction to Plasma Theory (John Wiley & Sons, New York, 1983).
[23] A. Maluckov et al., Physica A 322, 13 (2003).