Fusion Res.
Res. SERIES,SERIES,Vol.
Vol.6 (2004)
6 (2004)597–600
000–000
J. Plasma Fusion

Statistical Characteristics from Gyro-fluid Transport Simulation
MATSUMOTO Taro, KISHIMOTO Yasuaki and LI Jiquan
Japan Atomic Energy Research Institute, Ibaraki 311-0193, Japan(Received: 10 December 2003 / Accepted: 20 July 2004)

Abstract
Statistical characteristics are analyzed from gyrofluid simulation for sheared slab electron temperature gradientturbulence. It is found that the autocorrelation time is enhanced by zonal flow formation. For turbulent plasma, relatively high correlation dimension (~ 8) is obtained. On the other hand, the formation of the zonal flows is found tocause the reduction of the correlation dimension less than 3. Then, it is concluded that the variation of turbulent structure is well reflected in the correlation dimension. For turbulent plasma the deviation of probability density function
from Gaussian profile at high amplitude region is supposed to be the origin of intermittent large heat flux.
Keywords:
plasma turbulence, ETG (electron temperature gradient), zonal flow, statistical analysis, fractal,correlation dimension, probability density function, gyro-fluid simulation

zonal flow intensity, which characterize plasma performance,
have been systematically investigated by using gyrofluid
simulation of sheared slab ETG modes [5].

1. Introduction
For magnetic confinement fusion research, it is important to clarify and control the turbulent transport mechanismthat governs the performance of confinement plasma. In
recent years, it has been clarified that magnetic confinement
plasma autonomously generates a variety of structure such as
internal transport barrier through dynamics of the fluctuations with various temporal and spatial scales [1-2]. Then, it
is necessary to find out the rules to identify what kind of turbulent structure is generated there.
For understanding the transport mechanism in toroidal
plasmas, the statistical characteristics of the plasma turbulence have been investigated by experimental observations.
The fractal dimension of density and temperature fluctuations
was found to characterize the attractor of the turbulence [34]. It was shown that the fractal dimension is affected by thedifference of the plasma heating method [3]. The relation
between the tail component of probability density function(PDF) and the intermittency of turbulent structure have been
discussed in recent years. function (PDF) is affected by theturbulent structure. However, from a viewpoint of difficultyof preparing detail experimental condition, the knowledge ofsuch statistical characteristics enough to identify turbulent
structure has not been understood well.
Then, it is instructive to simulate microturbulence forfinding out the rules to identify turbulent structure by inductive method, in particular, from a viewpoint of feasibility ofpreparing detail condition. In this research, the relations
between statistical characteristics of turbulent structure andparameters such as temperature gradient, magnetic shear,

2. Simulation model and analysis method
2.1 Model equations
Gyrofluid models for ITG (ion temperature gradient)modes have been developed in recent years [6]. For ETG turbulence, gyro-fluid models including the Landau damping
effect are described from the symmetry to those of ITG turbulence, although the treatment for electron adiabatic
response to zonal flows is different from ion one [6]. Thus,
the gyrofluid equations for sheared slab configuration in thecollisionless limit are described asˆˆˆ 2 ) ∂ φ = [1 + (1 + η )∇
ˆ 2 ] ∂ φ + [φˆ ,∇ˆ 2 φˆ ](1 − ∇e⊥
⊥⊥ŷ
∂
∂ t̂
ˆ v̂ − µ ∇ˆ 4 φˆ,+∇|| ||⊥⊥∂ v̂ ||2ˆ
ˆ
ˆ
ˆ
ˆ
φ
φ
= ∇|| − ∇|| p̂ − [ ,v̂ ||] + η⊥∇ ⊥v̂ ||,
∂ t̂
∂ p̂
∂ φˆ
ˆ v̂
= − (1 + ηe)−[ φˆ, p̂ ] − Γ ∇|| ||∂ t̂∂ ŷˆ 2 p̂,− (Γ − 1) 8/π | k | (p̂ + φˆ ) + χ ∇||⊥⊥

(1)
(2)

(3)

where Γ = 5/3, and coordinates are normalized to characteristic scales as x̂ = x/ρe, ŷ = y/ρe, ẑ = Ln and t̂ = vtet/Ln. Here, theperturbed quantities, the electrostatic potential, the electronparallel velocity and the electron pressure, are normalized as∼φ̂ = (Ln /ρe)(eφ /Τe), v̂ || = (Ln/ρe)(ve||vte) and p̂ = (Ln/ρe)( p∼/p0),
with Ln = (d lnn/dx)-1, where ρe is the electron Larmor radius,©2004 by The Japan Society of Plasma
Science and Nuclear Fusion Research

Corresponding author’s e-mail: tmatsumo@naka.jaeri.go.jp
1597

Matsumoto T. et al., Statistical Characteristics from Gyro-fluid Transport Simulation

vte is the electron thermal velocity, Τe is the electron temperature, ve|| is the electron parallel velocity and n is the electrondensity. In this model, the adiabatic ion response is assumed,so that ηe = dlnΤe /dlnn. The magnetic shear parameter ŝ isdefined as ŝ = Ln/Ls, where Ls = B0 /(dBy /dx) and B0 is the longitudinal magnetic field.
The dynamics of fluid are computed in a three dimensional rectangular box with the Cartesian coordinates (x, y, z).
A periodic boundary condition is adopted in y and z directions, Numerical parameters are selected as Lx = 100ρe, Ly =10πρe, Lz = 2πρeLn. For numerically stable and efficient calculation, the viscosity terms are introduced in Eqs. (1) ~ (3)as a damping effect for the instabilities with wave length
smaller than electron Larmor radius, which does not change
the linear growth rate of ETG turbulence. Then, the perpendicular viscosity µ⊥, η⊥ and the perpendicular diffusivity χ⊥ areassumed to be 0.5. The detail of computation is described inRef. 5.

0 and H(r) = 0 for r < 0. For each m, C(r) would have a scale
like C(r) ~ rν. The correlation dimension D2 is obtained as
the exponent v saturated for enough embedding dimension m

3. Simulation results
We have simulated ETG turbulence in three cases. For
the case (A) where magnetic shear is relatively weak (ŝ =0.1), zonal flow can be obviously excited for large ηe. Then,electron heat conductivity is restricted. Kelvin-Helmholtz(KH) instability is also generated as a result of zonal flow
formation [5]. On the other hand, for small ηe, several modes
compete and the formation of zonal flow is restricted. For thecase (B), the magnetic shear is the same to (A). However,zonal flow component (0/0 mode) is artificially off to distinguish the pure effect of zonal flow component. For the case(C), magnetic shear is relatively strong (ŝ = 0.2). In eachcases, electron temperature gradient is changed within ηe = 3~ 6.
A typical behavior of time series data of poloidal electric field Ey for the strong ETG turbulence is shown in Fig.1.
There is no periodic behavior and large spiky bursts are intermittently observed. On the other hand, for the coherence withzonal flow, the time series data has large periodic oscillation
with small spiky fluctuation.

2.2 Autocorrelation time
Autocorrelation function of time series data, x (ti)[i =1,2,3,…, N ], is calculated as〈(x (t) − 〈x〉)(x (t + τ ) − 〈x〉)〉,(4)

where 〈…〉 means time average and τ is delay time. The autocorrelation time τac is defined as the delay time when theautocorrelation function has zero at the first time, so that 〈(x(t) − 〈x〉)(x (t + τac) − 〈x〉)〉 = 0.

3.1 Autocorrelation time
Figure 2 shows the autocorrelation time τac obtained
from Ey after the saturation of linear growing modes. For
weak shear case (A), the steep jump by one order is observed
from ηe = 3 to 4. This is because the effect of zonal flow formation, that is, the 0/0 mode becomes dominant for ηe ≤ 4.
On the other hand, τac is kept small regardless of ηe for (B)and (C). Then, it is found that the autocorrelation time isenhanced by zonal flow formation because long time scalecharacterizing KH mode appears through the dispersion relation.
We have also analyzed autocorrelation time for the radial electric field Ex, the electrostatic potential φ̂ and the heat

2.3 Correlation dimension
To evaluate the fractal dimension D of x(ti), it is necessary to reconstruct a trajectory in multi-dimensional space byembedding theorem. For one-variable time series data, it isappropriate to use a time-delay coordinate,{xi (ti), x(ti + τ ),…,x (ti + (m − 1)τ )},(5)in a viewpoint of noise reduction [7], where m is embedding
dimension and τ is embedding log (time delay). In order toobtain precise D, m should be larger than the twice of resultant fractal dimension D, i.e. m ≥ 2D + 1 [8]. The linear correlation between embedding coordinates would be too strong
for too small τ. On the other hand, the deterministic relation
between coordinates would disappear for too large τ. Then,the autocorrelation time τac would be appropriate for theembedding log τ.
To evaluate the fractal dimension, the correlation dimension D2 is useful from the point of view of computational feasibility. The correlation dimension can be calculated with
Grassberger-Procaccia algorithm [9] by taking correlation
integral of time series data asN∑C(r) = lim 12H(r − |xi − xj |)N→∞N i,j = 1(6)where H is the Heaviside function defined by H(r) = 1 for r ≥Fig. 1 Time series data of Ey for the case (A) ŝ and ηe = 3.

2598

Matsumoto T. et al., Statistical Characteristics from Gyro-fluid Transport Simulation

flux Q (= p̂Ey). The steep jump is also observed for all ofthese physical quantities. Then, hereafter, we mainly discuss
the statistical characteristics for Ey .

3.2 Correlation dimension
Figure 3 shows the correlation dimension D2 obtained
from Ey of the case (A). As embedding dimension, D2 linearly increases and saturates to the original correlation dimension of the dynamical system. For small temperature gradient(ηe = 3), relatively high dimension (D2 = 8 ~ 9) is obtained
because the strong turbulence with various temporal scales
has large degree of complexity.
However, as increasing ηe , D2 significantly decreases. In
particular, D2 becomes less than 3 for ηe = 6. At the sametime, the time series data have dominant periodic oscillation
with small fluctuations because the zonal flow excited by
ETG turbulence drives KH instability. Then, D2 is mainly
determined by the coherent oscillation rather than the turbulent fluctuations.
Figure 4 shows D2 as a function of ηe for the cases (A)Fig. 4 Correlation Dimension of Ey. The open circles andcrosses denote the case(A) ŝ = 0.1 with 0/0 mode andthe case(B) ŝ = 0.1 without 0/0 mode, respectively.

and (B). By comparing D2 for ηe = 3 and 6 of (B), it is foundthat the pure temperature gradient effect without zonal flow
is small (∆ D2 ~ 1). The increase of electron temperature gradient without 0/0 mode is not accompanied with the variationof turbulent structure. Then, the degree of complexity is not
so changed, although the electron heat conductivity isenhanced a lot.
On the other hand, by comparing D2 of (A) and (B) forηe = 6, the pure effect of zonal flow formation is relatively
large (∆ D2 ~ 7 for ηe = 6). Therefore, it is clarified that thevariation of turbulent structure is well reflected in the correlation dimension.

3.3 PDF
Figure 5 shows PDF of Ey for ηe = 3 and 6 of the case(A). In general, for strong turbulence, PDF is close toGaussian profile. Fig. 5(a) has good agreement with Gaussian
at small amplitude region. However, the deviation fromGaussian is observed at large amplitude region. This tailcomponent corresponds to the large spiky bursts observed intime series data of Fig.1.
In the PDF of heat flux for ηe = 3 of the case (A), thedeviation from Gaussian is found to be localized at positive
amplitude region. This asymmetry generates net heat flux.
Then, it is supposed that Ey and the pressure perturbation p̂
could synchronize at the intermittent burst,
On the other hand, Fig. 5(b) is also different fromGaussian because the periodic oscillation with small fluctuations is dominant in the presence of KH instability. Then, it isconcluded that the variation of turbulent structure is also wellreflected in PDF.

Fig. 2 Autocorrelation time of Ey . The open circles, crosses
and triangles denote the cases (A) ŝ = 0.1 with 0/0
mode, (B) ŝ = 0.1 without 0/0 mode and (C) ŝ = 0.2,
respectively.

4. Summary and discussion
Fig. 3 Correlation dimension of Ey. The open circles, crosses, triangles and squares denote ηe = 3, 4, 5 and 6,respectively. The broken line and dotted line mean
white noise case and the enough condition forembedding, respectively.

For sheard slab ETG turbulence, the relation between
the statistical characteristics and parameters characterizing
plasma performance have been systematically investigated byusing gyrofluid simulation. It is found that the autocorrela3599

Matsumoto T. et al., Statistical Characteristics from Gyro-fluid Transport Simulation

formation of the zonal flows causes the reduction of the correlation dimension less than 3. Therefore, it is clarified thatthe variation of turbulent structure is well reflected in the correlation dimension.
The deviation of probability density function fromGaussian profile is observed for both turbulent plasma andcoherent plasma with zonal flow. For turbulent plasma thedeviation at high amplitude region is supposed to be the origin of intermittent large heat flux. For coherent plasma, PDF
is very different from Gaussian due to the periodic oscillationexcited by KH instability.

Acknowledgements
This work was supported in part by a Grant-in-Aid for
Scientific Research from Japan Society for the Promotion ofScience.

References[1] A. Hasegawa and M. Wakatani, Phys. Rev. Lett. 59,1581 (1987).
[2] Y. Kishimoto et al., Plasma Phys. Control. Fusion, 40,A663 (1998).
[3] A. Komori et al., Phys. Rev. Lett. 73, 660 (1994).
[4] V. Budaev et al., Plasma Phys. Control. Fusion, 35, 429(1993).
[5] J. Li and Y. Kishimoto, Phys. Plasmas 9, 1241 (2002).
[6] G.W. Hammet et al., Phys. Rev. Lett. 64, 3019 (1990).
[7] N.H. Packard et al., Phys. Rev. Lett. 45, 712 (1980).
[8] F. Takens, in Dynamical Systems and Turbulence, edited by D.A. Rand and B.S. Young (Springer-Verlag,Berlin, 1981), p.366.
[9] P. Grassberger and I. Procaccia, Phys. Rev. Lett. 50, 346(1983).

Fig. 5 Probability distribution function of Ey for (a) ŝ = 0.1, ηe= 3 and (b) ŝ = 0.1, ηe = 6. The dashed lines indicate
Gaussian profile for reference.

tion time is enhanced by zonal flow formation because longtime scale characterizing KH mode appears through the dispersion relation.
For turbulent plasma, relatively high correlation dimension (~ 8) is obtained. On the other hand, it is found that the4600