Progress in NUCLEAR SCIENCE and TECHNOLOGY, Vol. 2, pp.72-77 (2011)

ARTICLEA New Simulation Method of Geodesic Acoustic Mode in Toroidal Plasmas byUsing Band-Limited White Noise in a δf Neoclassical Transport Code
Shinsuke SATAKE1,2,∗ , Hideo SUGAMA1,2 , Ryutaro KANNO1,2 , Takeshi IDO1,3 ,Seikichi MATSUOKA2 and Masayuki YOKOYAMA1,2
12National Institute for Fusion Science, Oroshi 322-6, Toki, Gifu 509-5292, JapanThe Graduate University for Advanced Studies (Sokendai), Oroshi 322-6, Toki, Gifu 509-5292, Japan3
Department of Energy Engineering and Science, Nagoya University, Nagoya 464-8603, JapanFORTEC-3D code, which solves the drift-kinetic equation for torus plasmas and radial electric field using the δfMonte Carlo method, has developed to study the variety of issues relating to neoclassical transport phenomena inmagnetic confinement plasmas. Here the numerical techniques used in FORTEC-3D are briefly reviewed, and recent
progress in the simulation method to simulate GAM oscillation is also explained. A band-limited white noise term
is introduced in the equation of time evolution of radial electric field to excite GAM oscillation, which enables us toanalyze GAM frequency with fine resolution even in the case the collisionless GAM damping is fast.
KEYWORDS: Monte Carlo, neoclassical transport, pseudo-random number generator, white noise, FokkerPlanck-Landau collision operator, GAM oscillation

I. Introduction

been utilized to calculate transport coefficients in LHD plasmas with including radial electric field, and it is found that the
non-local effect alters the equilibrium electric field profile. Its
wide applicability is also being demonstrated in researches ofrecently concerned topics such as the geodesic acoustic mode(GAM) of electric field oscillation and Landau damping of
it,5, 7) and simulation study of neoclassical toroidal viscosity
in tokamaks with weak perturbation magnetic field.9)
In this paper, first we briefly review the numerical methodof the δf Monte Carlo method, which solves the DKE for deviation of distribution function from local Maxwellian, δf =f − fM . The conserved-form Monte Carlo collision operator
and implementation of a parallelized pseudo-random numbers
generator (PRNG) used in it are also explained. Then recent
progress in the simulation method is shown about the excitation of GAM by introducing a band-limited white noise in
the time evolution of the radial electric field. It enables us toanalyze GAM frequency spectrum in LHD plasma, of which
GAM damping is strong, by carrying out a long-duration
GAM oscillation simulation with the excitation scheme.

Neoclassical transport theory1, 2) treats transport phenomena of plasmas in toroidal magnetic field configurations,which are governed by guiding-center motion and Coulomb
collisions of charged particles. In torus plasmas, particles ofwhich velocity pitch |v∥ /v| is small are trapped in a region
where magnetic field strength is weak, and those trapped particles mainly contribute the radial particle and heat transport.
The plasma distribution function tends to deviate from local
Maxwellian fM as the magnetic field is non-uniform and as
the plasma collisionality is low. Then precise information ofplasma distribution in the phase space is required to evaluate
the transport. In neoclassical transport theory, the time evolution of plasma distribution function is described by the driftkinetic equation (DKE), which describes time evolution of
gyro-averaged distribution function in five-dimensional phasespace (x, v∥ , v⊥ ).
From the early researches, DKE has been solved analytically or by computers adopting several approximations
and simplifications, some of which are not always proper
to simulate realistic cases in high-performance burning plasmas in forthcoming ITER experiments,3) for example. Recent
progress in computation environment enables us to solve DKE
directly and to evaluate transport coefficients more correctly.
We have developed FORTEC-3D code4–7) which is applicableto general 3-dimensional magnetic field like helical configurations of Large Helical Device (LHD) in NIFS.8) It treats
self-consistent evolution of radial electric field and the finiteorbit-width (FOW) effect of guiding-center motion, which areusually neglected in conventional methods. FORTEC-3D has

II. The δf Monte Carlo Method
FORTEC-3D solves the drift-kinetic equation for perturbed
distribution function δf = f − fM ,D
δf
Dt

where fM is a time-independent local Maxwellian which hasconstant density n and temperature T on a flux surface, vd ≡∗ Corresponding author, E-mail:satake@nifs.ac.jp

c

(
)
∂
∂δf
δf + v∥ + vd · ∇δf + K̇− CT (δf )
∂t(
∂K
)
∂
= − vd · ∇ + K̇
fM + PfM .
(1)
∂K
≡

2011 Atomic Energy Society of Japan, All Rights Reserved.
72

A New Simulation Method of Geodesic Acoustic Mode in Toroidal Plasmas by Using Band-Limited White Noiseẋ − v∥ b is the drift velocity, b is the unit vector along themagnetic field line, and the independent velocity variables are2chosen here as K = mv 2 /2 and µ = mv⊥/(2B). FORTEC3D adopts the Boozer magnetic coordinates10) (ψ, θ, ζ), whichis the toroidal magnetic flux within a flux surface, poloidal angle, and toroidal√ angle, respectively. The effective minor radius r = a ψ/ψedge is also used as a flux-surface variable,√
where a = 2ψedge /B0 , ψedge is the value of ψ at the outermost flux surface, and B0 is the magnetic field strength at themagntic axis (r = 0). Magnetic field in Boozer coordinates isgiven byB = ∇ψ × ∇θ + ι∇ζ × ∇ψ = G∇ζ + I∇θ + β∗ ∇ψ, (2)where G(ψ) and I(ψ) are the poloidal and toroidal current,and ι(ψ) is the rotational transform of the field line, respec√tively. The Jacobian of Boozer coordinates is gB = (G +ιI)/B 2 . We assume that only electrostatic radial electric field
E = −(dΦ/dψ)∇ψ = Eψ ∇ψ exists and the magnetic field
is steady in time. Note that Φ is allowed to be time-dependent.
To solve Eq. (1), FORTEC-3D adopts the two-weight
scheme of δf Monte Carlo method.11, 12) Each simulation maker, of which distribution function is written asg(x, v∥ , v⊥ ), has two weights which are defined so as to satisfy the relations wg = δf and pg = fM , respectively. Since
each marker moves according to the l.h.s. of Eq. (1), or
Dg/Dt = 0, time evolution of marker weights are given by[]p∂ẇ = −
vd · ∇ + K̇
− P fM ,
(3a)
fM
∂K
[
]
p
∂
ṗ =
vd · ∇ + K̇
fM .
(3b)
fM
∂K
The linearized Fokker-Planck-Landau collision operator(for like-particle collisions) is decomposed into the testparticle operator CT and the field-particle operator PfM . The
former is implemented in a Monte Carlo way by a random
walk of simulation marker velocity in the (v∥ , v⊥ ) space.13)
The latter is defined so that the collision operator satisfies thefollowing relations∫ 3
d vM (CT (f ) + PfM ) = 0
for M = {1, v∥ , v 2 },(4a)
CT (f ) + PfM = 0 for δf = (c0 + c1 v∥ + c2 v 2 )fM ,ci ∈ R,(4b)which describes the conservation property and the null-spaceof the linearized operator, respectively. The accuracy of theconservation property has improved recently by introducing
new numerical scheme for PfM operator.6)
The time evolution of the radial electric field on each flux
surface is given from the time derivative of Poisson’s equationin Boozer coordinates as follows:
(
⟨
)⟩
∑
c2
∂
2
ea ΓD
ϵ0 |∇ψ| 1 + 2
Eψ (ψ, t) = −
(5)
a ,
vA
ψ ∂t
a
where ⟨· · ·⟩ denotes a flux-surface average defined asH√dθdζ gB A(ψ, θ, ζ)H⟨A⟩ψ ≡,√dθdζ gB
VOL. 2, OCTOBER 201173

vA = ϵ0 cB/(mi ni )1/2 is the Alfvén velocity, and the subscript a specifies the particle species. The radial particle andheat fluxes are evaluated from the guiding-center radial drift
velocity as follows:⟨∫⟩D3
Γa (ψ, t) =d vδfa vd · ∇ψ
.
(6a)
ψ

⟨∫
QD
a (ψ, t) =

d3 v

ma 2
v δfa vd · ∇ψ
2

⟩
.

(6b)

ψ

In neoclassical transport simulations for LHD plasmas, determination∑of the ambipolar electric field profile, which evolves
so that a ea ΓD
a = 0 is satisfied on each flux surface, is oneDof the main tasks since ΓDa and Qa in non-axisymmetric plasmas are strongly affected by E × B rotation speed. In order to reduce the computation time, FORTEC-3D only solves
ion neoclassical transport, and another code GSRAKE14) isused to evaluate ΓDe , which solves a bounce-averaged approximated DKE for helical configurations. In solving DKE
for ions, FORTEC-3D neglects the ion-electron collisions because me /mi ≪ 1.
We have applied this method so far, because the finiteorbit-width effect and the strict conservation property of likeparticle collision term, which are the merit of FORTEC-3D,
are considered to be more important for ion transport analysis
than that for electron, because of the large mass difference.
In the recent LHD experiments, however, core electron temperature becomes very high in Core Electron-Root Confinement (CERC) plasmas,15) and positive ambipolar radial electric field is considered to play an important role to achieve
good electron confinement. Since the Te gradient in CERC isvery steep and trapped electrons in such high-Te plasmas isinsensitive to Coulomb collisions and can drift widely in radial direction, it is concerned that the FOW effect on electron
neoclassical transport in CERC plasmas would not be negligible. To simulate neoclassical transport and ambipolar electricfield in CERC plasmas, FORTEC-3D is being extended forelectron transport simulations. For that purpose, the unlikeparticle (electron-ion) collisions operator as well as the likeparticle ones is included. From benchmarks of the code,16)
it is found that the radial electron particle and heat fluxes inhigh-Te LHD plasma are largely different from those evaluated from conventional neoclassical calculations especially if
the E × B rotation is weak, which suggests the importanceof precise calculation method of neoclassical transport in analyzing high-Te core plasmas.
FORTEC-3D is written in Fortran90 and parallelized byusing MPI. In the parallelization of a massive Monte Carlocode, pseudo-random number generator (PRNG) may be abottle neck. In a typical simulation run, FORTEC-3D uses
about 1 ∼ 4 × 107 markers and calculates about 104 ∼ 105
times collisions, which requires a fast, long period, and statistically good quality PRNG. We have adopted a parallel
PRNG according to Mersenne twister17) with the dynamic creation method,18) which generates independent random number sequences on each MPI process. FORTEC-3D is usually executed using 64 to 256 MPI processes with 8 parallel threads on shared memory, and the time cost of PRNG isless than 1% of the total simulation time. The parallel PRNG

Shinsuke SATAKE et al.

74

has passed a statistical independency test.19) The Mersenne
twister in FORTEC-3D can manage up to 1,024 MPI parallel
computation now, though increase of communication time of
MPI_REDUCE-type operations which are required in solving
PfM operator and Eq. (6) becomes bottle neck as the numberof MPI processes increases.

III. Excitation of GAM Using White-Noise
Zonal flows and the geodesic acoustic mode (GAM) oscillation in toroidal plasmas have been investigated recently, as
the zonal flow suppresses micro instabilities and anomaloustransport.20) Cross-device investigation21) has revealed bothsimilarities and differences in the nature of GAMs among several devices. One possibility to explain it is the magnetic fieldstructure. Theoretical studies on GAMs in non-axisymmetricmagnetic field based on the gyrokinetic model by Sugama and
Watanabe22) have shown that the GAM frequency and damping rate depend on the magnetic field spectrum, and the collisionless damping of GAM oscillation in helical plasma isfaster than that in a comparable scale tokamaks. To verify
the analytic expectation of GAM in LHD magnetic configurations, we have carried out a simulation study of the dependence of GAM frequency and damping rate using FORTEC3D.5) It has been shown there that GAM frequency ωG islower and the collisionless damping rate of GAM γG is higher
when the magnetic axis position of LHD is shifted inward. We
have also shown there that the DKE simulation of FORTEC3D corresponds to the Te /Ti → 0 limit of the gyrokinetic theory. Since executing global simulations of GAM in realistic
configurations in a gyrokinetic model is time-consuming, researching GAM using a DKE model is considered as a usefultool. The reliability of FORTEC-3D for GAM simulation has
been proved by a benchmark7) with full-f gyrokinetic Vlasovcode GT5D.23)
In the previous study, however, we have used MHD equilibrium magnetic field of LHD with the rotational transformprofile ι is fixed to a low value in surveying the GAM dependence on the magnetic axis position, ι = 0.35 to 0.5 at
r = 0 ∼ 0.8, since the GAM damping rate is roughly proportional to ι, and strong Landau damping occurs when ι ≃ 1.
In LHD experiments the ι profile is determined by the internal
plasma current profile and magnetic axis position, usually increasing with minor radius. However, even if ι ≪ 1, the helical magnetic perturbations in LHD also causes rapid dampingof GAM. Because of the strong GAM damping, it is difficult
to analyze the frequency spectrum precisely from simulation
result, since in evaluating GAM frequency by power-spectrum
analyze, long-duration data of ΓD (ψ, t) and Eψ (ψ, t) are required to obtain a fine resolution in frequency. In the DKE
model, however, there is basically no driving force of GAM,
and the initial GAM in the simulation is excited just because
the initial condition δf = Eψ = 0 is not an equilibrium stateof toroidal plasma.
To resolve this problem, we devised a new technique of exciting GAM in DKE simulation. A white-noise term is introduced in the equation of time evolution of Er on a flux surface,which is described here as a stochastic differential equation as

follows:
dEr (r, t) = −

)
e (
D
Zi ΓD
i (r, t) − Γe (r, t) dt
ϵ0 ϵ⊥+σE (r)dW (r, t),(7)where the radial coordinate is transformed from ψ to r for simplicity, and the "diffusion coefficient" of radial electric field
σE is considered as a very small factor. Equation (7) can be
regarded as a random-walk of radial electric field around an
ambipolar state which mimics the effect of any turbulence onelectrostatic potential in plasma. Concerning the noise term, aband-limited white-noise, of which frequency range contains
the GAM frequency ωl < ωG < ωu , is adopted :
√
N
2dt ∑
dW (r, t) =cos[ωk t + ϕk (r)],
(8)
N
k=1

where N is sufficiently large number (≫ 100), ωk = ωl +(k−1/2)∆ω, k = {1, 2, · · · , N }, ∆ω = ωu − ωl , and the phaseϕk (r) are random numbers in [0, 2π) given independently oneach radial simulation mesh. Since the ambipolar electric fieldis O(Ti /(ea)) and ωG ∼ vti /R0 where vti is ion thermalvelocity and R0 is the major radius, the diffusion coefficientσE is given by()2
Ti (r)
vti (r)2σE(r) = δ 2
(9)
,
ea
R0
with a control parameter δ ≪ 1. Band-limited white noise
is used so that GAM oscillation can be excited efficiently bysmall noise amplitude without bringing very high and low frequency components which are irrelevant to GAM excitation
into the simulation. The idea of application of white noise isborrowed from the numerical method used in the researches ofseismic resistant and structural analyses of buildings.24) Since
white-noise consists of continuum spectra, one can simulate
the response of a structure to general earthquakes, without
knowing details of the spectrum of "real" earthquake wave
nor the eigenfrequncies of the structure. Similarly, though thephysical driving mechanism of zonal flows and GAMs is regarded as nonlinear couplings of micro instabilities, which isbeyond the drift-kinetic approach, we can still study the behavior of excited GAM with DKE by using white-noise instead of unknown driving term.
Numerically, Eq. (7) is solved by the Euler method, while
marker motion and weight evolution are solved by the 4-th order Runge-Kutta method. To suppress the discretization error
in time evolution, the time step size dt is chosen small enough,ωG dt, ωu dt ≪ 10−2 . We have checked the simulation results(frequency and amplitude of excited GAM) are not affected by
the choice of dt and N in Eq. (8) by varying these parameters.
To demonstrate the numerical scheme, we carried out several simulations in a LHD magnetic field configuration (R0 =3.75 m , B0 = 1.5 T, without plasma current) obtained froma MHD equilibrium calculation code VMEC.25) FORTEC-3D
adopts a Fourier series expression of magnetic field spectra as
shown in Ref. 5). We used 22 components of Bm,n spectrain the present study. Core ion density and temperature are
ni0 = 2.0 × 1018 m−3 and Ti0 = 1.0 keV, respectively.
PROGRESS IN NUCLEAR SCIENCE AND TECHNOLOGYA New Simulation Method of Geodesic Acoustic Mode in Toroidal Plasmas by Using Band-Limited White Noiser/a=0.2104d=0.00d=0.01
d=0.034e+18

Γi[1/m s]

0

2e+18

2

Er[kV/m]r/a=0.210
6e+18

d=0.00d=0.01
d=0.03noise(d=0.01)
noise(d=0.03)

2

75

0

-2
-2e+18-4-4e+18
0

0.02

0.04

0.06

0.08
0.1
time [t/τi]

0.12

0.14

0.16

0

0.02 0.04 0.06 0.08 0.1
time [t/τi]

0.12 0.14 0.16

Fig. 1 Time evolution of radial electric field Er in GAM ex-

Fig. 2 Time evolution of radial particle flux Γi in GAM ex-

citation simulations observed at r/a = 0.21 surface. The
time is normalized by the ion collision time τi . The amplitude of the white-noise δ is varied 0.00, 0.01, and 0.03
in these simulations. The noise parts are plotted at the offcenter position in the graph with magnified by 200.

citation simulations observed at r/a = 0.21 surface. The
amplitude of the white-noise δ is varied 0.00, 0.01, and 0.03
in these simulations.

The ι profile is monotonically increasing, ι(r = 0) = 0.34,
ι(0.56a) = 0.50, and ι(a) = 1.2. The collisionality withthese parameters was very low, and then the effect of collisions on GAM damping is negligible. The plasma configuration is made similar to a LHD experiment with observation of
GAM.26, 27)
The time evolution of the radial electric field and particle
flux in the simulations are shown in Figs. 1 and 2, where threecases (δ = 0.00 (without noise), 0.01, and 0.03) are plotted.
In Fig. 1, the given noise portions in δ =0.01 and 0.03 cases
are also plotted. Note that the noise amplitude is very small
compared to the total Er , then they are magnified by 200 in the
figure. The GAM frequency ωG is approximately estimated as
ωG0 =√7/4(vti /R0 )(1 + 46/49q 2 )1/2

(10)

in the Te /Ti → 0 limit of simple axisymmetric tokamak
configuration. In the present cases, it ranges about 26 kHz(r = 0) > ωG0 /(2π) > 17 kHz(r = a). Therefore,
the band-limited white noise was given in the range fromωl /(2π) = 10 kHz to ωu /(2π) = 36 kHz. In Figs. 1 and
2, one can see that GAM was successfully excited continuously in time with very small amplitude of white noise. Comparing the amplitude of the noise part Ẽr and the total Er ,the ratio is about |Ẽr /Er | ∼ 1/100, and excited GAM amplitude is roughly proportional to σE . On the other hand, if
the white noise was turned off, initial GAM oscillation disappeared soon because of the strong Landau damping of GAM
in helical configurations, and only a small fluctuation fromMonte Carlo noise remained.
From the data of continuous wave form of GAM in particle
flux ΓDi (r, t) and radial electric field Er (r, t), we can analyze
the GAM frequency by taking the power spectrum of themas it is shown in Figs. 3(a), (b) and (c). For δ = 0.01 case,the peak of the spectrum appears slightly below the estimation
ωG0 , which is an approximation by neglecting the helicity ofVOL. 2, OCTOBER 2011

LHD magnetic field. The band-limited white noise has uniform power spectrum in the range from ωl to ωu as shwon
in Fig. 3(c). On the other hand, we cannot see clear peakof power spectrum if the white noise is cut as in Fig. 3(a).
Though we do not show the power spectrum of ΓDi , it has almost the same profile as that of Er . The radial position whereclear peak of GAM oscillation is r/a = 0.05 ∼ 0.20, which
seems consistent with the observation of GAM using HeavyIon-Beam-Probe (HIBP) in a LHD experiment,27) of which
profile is similar to the one in the simulation. We have also
confirmed that the GAM frequency and the peak radial position of GAM in simulation are almost unaffected by varying
the white-noise amplitude from δ=0.01 to 0.03.
Next, the peak frequencies and peak amplitudes of thepower spectrum of ΓDi and Er are compared with the gyrokinetic theory,22) which takes into account of the helical magnetic ripples of LHD configurations. Analytic solutions ofωG and γG for the gyrokinetic analysis in the Te /Ti → 0
drift-kinetic limit is shown in Eqs. (7) and (8) in Ref. 5). Figure 4(a) shows the radial profile of GAM frequency. Notethat ωG0 is an approximate estimation of GAM frequency in atokamak geometry. Both gyrokinetic analysis and FORTEC3D simulation show lower ωG than ωG0 , which is the effect ofhelical configuration of LHD magnetic field. The radial profile of the peak of power spectrum amplitudes of ΓDi and Erare plotted in Fig. 4(b). The peak of amplitude is located at
r/a ≃ 0.125, and the amplitudes decay at outboard region.
It was diffucult to find a clear peak in power spectrum andidentify ωG from simulation results at r/a > 0.35. It is because both ι and the helical ripple components of magnetic
field spectrum are increasing function of r in the LHD configuration used here, which results in the larger damping rate of
GAM, |γG |. One can see that the peak amplitudes observedin FORTEC-3D simulation tend to be proportional to |1/γG |,which means that the GAM excited by the white noise actually damps according to the kinetic mechanism that gyrokinetic theory explains. Recent analysis of LHD experiment has
found that the spacial profile of GAM power shows the sim-

Shinsuke SATAKE et al.

76ilar tendency as Fig. 4(b).27) This is simply explained usinga model of forced vibration of an oscillator of which natural
frequency is ω0 ,(a) Power Spectrum Er : d=0.00
0.25
0.2

40

0.15

30

0.1

20

f[kHz]

mẍ + γω0 ẋ + ω02 x = f (t),

0.05
10
0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

r/a

(b) Power Spectrum Er : d=0.01
7
6

40
30

4
3

20

2
1

f[kHz]

5

0
0.1

0.2

0.3

0.4

0.5

0.6

0.7

r/a

(c) White noise : d=0.01
3e-06
40

20

1e-06

f[kHz]

2e-06
30

10
0
0

0.1

0.2

0.3

0.4

0.5

0.6

where γ is the damping rate. If the driving force f (t) is whitenoise of which power spectral density is S(ω) = S0 , the variation of the displacement σx2 = E[x2 (t)] is known to be proportional to S0 /(γω03 ),24) where E[g] means the expected value
2of g. Note here that σE∝ S0 in Eq. (7). Though Eq. (11) isa simple model of GAM damping by neglecting radial correlations, it explains that E[|Er − Er |2 ] ∝ |γG |−1 . We could
have seen GAM at r/a > 0.35 if the driving term had beenmuch stronger, however, the coincidence of the radial profileof the GAM power between simulations and observation implies that the realistic choice of the GAM driving term δ is theorder of 10−2 . More investigation by comparing the GAM
oscillation amplitude between simulation and observation isplanned near future.
Finally, at the closest position to the magnetic axis r/a =0.05, difference between simulation and theory is large. It
is possibly because of the FOW effect of passing particle orbit around the magnetic axis, of which radial drift width is
wide there (so called potato orbits), or because of the effect ofboundary condition Er (r = 0) = 0 in the simulation. The gyrokinetic analysis do not assume the existence of potato orbits
nor the magnetic axis.

IV. Summary

10
0

(11)

0.7

r/a

Fig. 3 Contour plots of the power spectrum of Er for (a):

δ = 0.00 and (b): δ = 0.01 case, respectively, and (c): thepower spectrum of the white-noise part for δ = 0.01 case.
Green line represents an approximate estimation of GAM
frequency for tokamak case, Eq. (10). Note that the scaleof the color contour is different among these figures.

We have developed a new simulation scheme to simulate
GAM in LHD plasmas, in which the GAM damping rate is
essentially high compared with that in tokamaks, by applying a band-limited white noise as a driving term of GAM oscillation in FORTEC-3D DKE Monte Carlo code. From thetest simulations in a LHD magnetic configuration it has been
demonstrated that the GAM oscillation can be excited by verysmall amplitude of white noise, and the GAM frequency andthe peak amplitude of the power spectrum of GAM oscillation are consistent with a gyrokinetic analysis for helical plasmas. With the GAM excitation scheme FORTEC-3D simulation is able to simulate GAM oscillation and damping in realistic LHD plasma profile and analyze GAM power spectrumwith fine resolution. LHD is equipped with HIBP which observes electrostatic potential profile in core plasma. GAM frequency profile can be obtained from the diagnostics, and thenwe can carry out precise comparisons among observation, gyrokinetic theory, and FORTEC-3D simulation, which will be
profitable for further understanding the physics of GAMs andzonal flows in helical plasmas. Further improvement of thesimulation model is planned in near future to include the impurity heavy ion effect, which is one of the candidates thatchange the observed GAM frequency from the analytic estimation.

Acknowledgments
This work is supported in part by the Japanese Ministry of Education, Culture, Sports, Science, and Technology, Grant No. 21560861 and in part by the NIFS collabPROGRESS IN NUCLEAR SCIENCE AND TECHNOLOGYA New Simulation Method of Geodesic Acoustic Mode in Toroidal Plasmas by Using Band-Limited White Noise(a) Peak frequency
1.11
ωG/ωG0

0.9
0.8
0.7
0.6

Peak(Er)
Peak(Ion Flux)
G-K theory

0.5
0.4
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.3

0.35

r/a

Power amp., 1/|γ| (Arb. unit)(b) Peak amplitude
Peak(Er)
Peak(Ion Flux)
1/|γ|(G-K theory)

1

0.1

0.01
0

0.05

0.1

0.15

0.2

0.25

r/a

Fig. 4 (a): Comparisons of GAM frequency estimated frompeak positions of power spectrum of Er and ΓDi (symbols),
and an estimation from gyrokinetic analysis5, 22) (line). The
vertical axis is normalized by an approximate estimation fortokamak case, Eq. (10). (b): Comparisons of the peak amplitudes of power spectrum of Er and ΓDi (symbols), andthe inverse of GAM damping rate 1/|γG | from gyrokinetic
analysis (line). FORTEC-3D results for δ = 0.03 case areused here.
orative Research Programs, NIFS09KTAD006, 09KTAL022,
10KDAT020, 10KNXN186, and 10KTAT040.
References
1) R. Balescu, Transport Processes in Plasmas Vol. 1 & 2, NorthHolland, Amsterdam (1998).
2) P. Helander, D. J. Sigmar, Collisional Transport in Magnetized
Plasmas, Cambridge Univ. Press, Cambridge (2002).
3) ITER Physics Basis Editors et al., Nucl. Fusion, 39, 2137
(1999).
4) S. Satake, M. Okamoto, N. Nakajima, H. Sugama, M.
Yokoyama, C. D. Beidler, "Non-Local Simulation of the Formation of Neoclassical Ambipolar Electric Field in NonAxisymmetric Configurations," Plasma Fusion Res., 1, 002(2006).
5) S. Satake, H. Sugama, T.-H. Watanabe, "Simulation studies on
the GAM oscillation and damping in helical configurations,"
Nucl. Fusion, 47, 1258 (2007).
6) S. Satake, R. Kanno, H. Sugama "Development of a NonLocal Neoclassical Transport Code for Helical Configurations,"
Plasma Fusion Res., 3, S1062 (2008).
7) S. Satake, Y. Idomura, H. Sugama, T.-H. Watanabe, "Benchmark test of drift-kinetic and gyrokinetic codes through neoclassical transport simulations," Comput. Phys. Comm., 181, 1069
(2010).
VOL. 2, OCTOBER 2011

778) A. Iiyoshi, A. Komori, A. Ejiri et al., "Overview of the Large
Helical Device project," Nucl. Fusion, 39, 1245 (1999).
9) S. Satake, H. Sugama, R. Kanno, J.-K. Park, "Calculationof Neoclassical Toroidal Viscosity in Tokamaks with Broken
Toroidal Symmetry," Plasma Phys. Contr. Fusion, 53, 054018(2011).
10) A. H. Boozer, "Plasma equilibrium with rational magnetic surfaces," Phys. Fluids, 24, 1999 (1981).
11) W. X. Wang, N. Nakajima, M. Okamoto, S. Murakami, "A newδf method for neoclassical transport studies," Plasma Phys.
Contr. Fusion, 41, 1091 (1999).
12) S. Brunner, E. Valeo, A. Krommes "Collisional delta-f schemewith evolving background for transport time scale simulations,"
Phys. Plasmas, 6, 4504 (1999).
13) X. Q. Xu, M. N. Rosenbluth, "Numerical simulation of iontemperature-gradient-driven modes," Phys. Fluids, B3, 627(1991).
14) C. D. Beidler, H. Maaßberg "An improved formulation of theripple-averaged kinetic theory of neoclassical transport in stellarators," Plasma Phys. Contr. Fusion, 43, 1131 (2001).
15) M. Yokoyama, H. Maaßberg, C. D. Beidler et al., "Common
Features of Core Electron-Root Confinement in Helical Devices," Fusion Sci. Technol., 50, 327 (2006).
16) S. Matsuoka, S. Satake, M. Yokoyama, A. Wakasa, S. Murakami, "Neoclassical electron transport calculation by using δfMonte Carlo method," Phys. Plasmas, 18, 032511 (2011).
17) M. Matsumoto, T. Nishimura "Mersenne twister: a 623dimensionally equidistributed uniform pseudo-random number
generator," ACM Trans. Model. Comput. Simulat., 8, 3 (1998).
18) M. Matsumoto, T. Nishimura, "Dynamic Creation of Pseudorandom Number Generators," Monte Carlo and Quasi-Monte
Carlo Methods 1998, Springer, 56 (2000).
19) S. Satake, M. Okamoto, N. Nakajima, H. Takamaru, "Development of three-dimensional neoclassical transport simulation
code with high performance Fortran on a vector-parallel computer," Lecture notes in Computer Science, 4759, SpringerVerlag, Berlin (2008).
20) K. Itoh, S.-I. Itoh, P. H. Diamond et al., "Physics of zonal
flows," Phys. Plasmas, 13, 055502 (2006).
21) A. Fujisawa, T. Ido, A. Shimizu et al., "Experimental progresson zonal flow physics in toroidal plasmas," Nucl. Fusion, 47,
S718 (2007).
22) H. Sugama, T.-H. Watanabe, "Collisionless damping of zonal
flows in helical systems," Phys. Plasmas, 13, 012501 (2006).
23) Y. Idomura, M. Ida, T. Kano, N. Aiba and S. Tokuda, "Conservative global gyrokinetic toroidal full-f five-dimensional
Vlasov simulation," Comput. Phys. Comm., 179, 391 (2008).
24) M. Hoshiya, Kakuritsuron shuho ni yoru shindo kaiseki, Kajima
Kenkyujo Shuppankai, Chapters 3 and 4 (1974), [in Japanese].
25) S. P. Hirshman, O. Betancourt, "Preconditioned descent algorithm for rapid calculations of magnetohydrodynamic equilibria," J. Comput. Phys., 96, 99 (1991).
26) T. Ido, A. Shimizu, M. Nishiura et al., "Experimental study ofradial electric field and electrostatic potential fluctuation in the
Large Helical Device," Plasma Phys. Contr. Fusion, 52, 124025(2010).
27) T. Ido, A. Shimizu, M. Nishiura et al., "Potential fluctuation associated with the energetic-particle-induced geodesic acoustic
mode in the Large Helical Device," Nucl. Fusion, 51, 073046(2011).