J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

Gyrokinetic simulation of entropy cascadein two-dimensional electrostatic turbulence
T. TATSUNO1) , M. BARNES1,2,3) , S. C. COWLEY3,4) , W. DORLAND1) , G. G. HOWES5) ,R. NUMATA1) , G. G. PLUNK1,6,7) , and A. A. SCHEKOCHIHIN2,4)1)

IREAP, University of Maryland, College Park, Maryland 20742, USA
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3NP, UK
3)
Euratom/UKAEA Fusion Association, Culham Science Centre, Abingdon OX14 3DB, UK
4)
Plasma Physics, Blackett Laboratory, Imperial College, London SW7 2AZ, UK
5)Department of Physics and Astronomy, University of Iowa, Iowa City, Iowa 52242, USA
6)
Department of Physics and Astronomy, University of California, Los Angeles, California 90095, USA
7)
Wolfgang Pauli Institute, University of Vienna, Vienna A-1090, Austria
2)

February 8, 2010

(Received: 20 November 2009 / Accepted: 24 February 2010)

Two-dimensional electrostatic turbulence in magnetized weakly-collisional plasmas exhibits a cascade ofentropy in phase space [Phys. Rev. Lett. 103, 015003 (2009)]. At scales smaller than the gyroradius, this cascadeis characterized by the dimensionless ratio D of the collision time to the eddy turnover time measured at the scaleof the thermal Larmor radius. When D � 1, a broad spectrum of ﬂuctuations at sub-Larmor scales is found inboth position and velocity space. The distribution function develops structure as a function of v⊥ , the velocity
coordinate perpendicular to the local magnetic ﬁeld. The cascade shows a local-scale nonlinear interaction inboth position and velocity spaces, and Kolmogorov’s scaling theory can be extended into phase space.
Keywords: electrostatic turbulence, gyrokinetic theory, nonlinear phase mixing, kinetic dissipation, decaying
turbulence simulation, Kelvin-Helmholtz instability

1. Introduction

and the system is in principle reversible in a thermodynamic sense [10]. Wave resonances or Landau damping do
not by themselves constitute dissipation since no entropy
is generated unless we include the eﬀect of collisions. Irreversibility becomes the key issue here.
In this context, the following question arises. Collisions are usually described by a diﬀusion operator in
the velocity space [11]. How do collisions in a weakly
collisional system determine the wave-number as well as
velocity-space cutoﬀ of the inertial range? In this paper, following Ref. [12], we show detailed numerical evidence of the coupling between position- and velocity-spacestructures, and the consequent achievement of the collisional dissipation described by the velocity-space diﬀusion operator. We introduce a new dimensionless number
D, which, analogously to the Reynolds number in NavierStokes equations, characterizes the scale separation between the outer scale and the dissipation scale in the kinetic theory. We emphasize that we are concerned onlywith scales smaller than the Larmor radius.
In order to focus on the nonlinear interaction, we omit
Landau damping by ignoring variation along the ﬁeld line.
The system retains two collisionless invariants as the 2D
Navier-Stokes equations do. These two invariants cannot
share the same local-interaction space in a Kolmogorovlike phenomenology, and thus lead to a dual cascade (direct
and inverse cascades) [3,13–16]. In this paper we focus on
the direct cascade only and the discussion of the inverse
cascade is left for future publication.
This paper is organized as follows. As a reduced ki-

Fluid turbulence can be described mostly by NavierStokes equations, which consist of partial diﬀerential equations in the position space — thus in 2D or 3D space. Sucha description is possible because we may regard the ﬂuid
to be suﬃciently collisional and to be locally in thermodynamic equilibrium. From the outer driving scales to thesmall dissipative scales, the energy (or the enstrophy in the
2D case) cascades in the wave-number space [1]. The inertial range dynamics are basically governed by the nonlinear term. The ﬂuxes are ﬁnally dissipated by the viscosity
described by a diﬀusion operator, in the dissipation rangeof wave-number space. There is a dimensionless number
called the Reynolds number, which characterizes the scale
separation of the turbulent dynamics. The ratio of the dissipation scale to the outer scale is in general determined bya fractional power of the Reynolds number.
Kinetic turbulence is less well understood. We havesome knowledge especially from the extensive study of fusion plasmas [2–4] and recently of space plasmas [5–9],
but not to the same extent as ﬂuid turbulence is known, particularly in the sense described above. Since most plasmas
are in a collisionless or weakly-collisional environment,
they are not in local thermodynamical equilibrium. In thiscase we have to invoke kinetic description in the phasespace, which enables us to use Vlasov or Boltzmann theory, or their magnetized low-frequency limit, the gyrokinetics. Dissipation in the kinetic system occurs becauseof collisions. Without collisions, the entropy is conserved
author’s e-mail: tatsuno@umd.edu

509©2010 by The Japan Society of Plasma
Science and Nuclear Fusion Research

T. Tatsuno et al., Gyrokinetic Simulation of Entropy Cascade in Two-Dimensional Electrostatic Turbulence

where �·� r denotes the gyroaverage at ﬁxed particle posi�tion r, Q = s q2s n0s /T 0s for Boltzmann-response (3D)electrons or Q = q2i n0i /T 0i for no-response (2D) electrons,
n0 is the density corresponding to F0 , and s and i are thespecies indices. In this paper we use no response electrons
since electrons cannot contribute to the potential if k� = 0
exactly [22]. We have also made a separate two-species
simulation with real mass ratio and obtained almost identical results in the regime k⊥ ρe � 1 with negligible density
ﬂuctuations, where ρe is the electron Larmor radius.
There are two collisionless invariants,��Q
T 0 h2
ϕ2 dr,=dR dv −
(4)
2F0
2
Q�
(1 − Γ0 )|ϕ k |2 ,
(5)
=
2 k

netic model for plasma turbulence, we ﬁrst introduce thegyrokinetic (GK) equation brieﬂy and describe its basic
nature in Sec. 2. Scaling relations of the turbulent cascade are also brieﬂy summarized. Section 3 shows theevidence of kinetic turbulent cascade by means of the numerical simulation of the decaying turbulence. In addition
to the reproduction of the theoretical prediction, we show
how much resolution is needed in both position and velocity space to achieve the proper dissipation and show thetrend that the amount of dissipation is asymptotically independent of the collision frequency. We also present thecharacteristics of the nonlinear triad interaction in wavenumber space and show that the assumption of local-scaleinteraction is supported. We ﬁnally summarize our resultsin Sec. 4.

where the subscript k is for the Fourier coeﬃcients, Γ0 =I0 (b)e−b , I0 is the modiﬁed Bessel function, b = k⊥2 ρ2 /2,and ρ is the ion thermal Larmor radius. The invariant
is proportional to the negativeof the perturbed part of the�entropy of the system, − f ln f dr dv [2, 23], where f is
the full distribution function. Here we will refer to as“entropy” to emphasize this connection, although in thermodynamic terms, it is better interpreted as a generalized
free energy [24, 25]. The second invariant implies theexistence of a dual cascade in the 2D system, which we
do not discuss in this paper. Note that these quantities areconserved only in the collisionless limit. Any decrease of
due to collisions corresponds to the creation of entropy
and heating, hence the irreversibility.

2. Gyrokinetics
2.1 Equations

We ﬁrst introduce the GK model brieﬂy [17–20].
Since we are going to look at the turbulence in the magnetized plasmas, the dynamics of interest is much slower
than particles’ gyromotion. We may average over the gyromotion and the gyroangle can be ignored due to this gyroaveraging. A GK system has 3 spatial coordinates (x, y,z), and 2 velocity coordinates (v⊥ , v� ), where � and ⊥ standfor parallel and perpendicular directions to the backgroundmagnetic ﬁeld. It is convenient to distinguish between theparticle coordinate r and the gyrocenter coordinate R [17].
These coordinates are connected by
R= r+v × ẑ
,
Ω

(1)

2.2

where ẑ is the unit vector along the background magnetic
ﬁeld and Ω is the gyrofrequency.
We further reduce the GK equation into 2D in position
space, or 4D in phase space, by ignoring variation along
the mean ﬁeld (k� = 0). This is useful because we not only
reduce the dimension but also remove the Landau damping
from the system. The resultant equation isc∂h
qF0 ∂�ϕ� R+ {�ϕ� R , h} =+ �C[h]� R ,∂t B0
T 0 ∂t

Turbulent scaling

A scaling theory of the entropy cascade in the subLarmor scale range (� � ρ) can be developed in a way
reminiscent of the Kolmogorov-style turbulence theories.
From the conservation of in the collisionless limit, we
may assume that the entropy is going to cascade scale byscale without dissipation in the inertial range. Note thatthe quantity which cascades forward (to smaller scales)
is , not . This has been demonstrated by applying
Fjørtoft’s argument on the dual cascade [13, 25].
The nonlinear term [the Poisson bracket in (2)] thatcauses the entropy transfer introduces velocity-space structure simultaneously. For small-scale electric ﬁelds, particles with diﬀerent gyroradii execute diﬀerent E × B motions because they see diﬀerent eﬀective potentials; thisleads to nonlinear phase mixing and other novel phenomena [22,26]. It can be argued that the nonlinear phase mixing produces the following relationship between positionspace and velocity-space scales:(2)

where h is the non-Boltzmann part of the perturbed ion
distribution function, ϕ is the electrostatic potential, B0 is
the background magnetic ﬁeld aligned with the z-axis, �·� R
is the gyroaverage holding the guiding center position R
constant, { f, g} = (ẑ × ∇ f ) · ∇g is the standard Poissonbracket, q is the charge, and T 0 is the temperature of the
background Maxwellian F0 . The collision operator C[h]we use in our simulations describes pitch-angle scatteringand energy diﬀusion with proper conservation properties(see [21]).
To calculate the self-consistent potential at these small
scales, we use the quasineutrality condition�Qϕ = q �h� r dv
(3)

δv⊥
�
∼ .
vth
ρ

(6)

Dimensional arguments then lead to the following
scalings of the spectra (see Refs. [12, 25, 27] for details),
Eh (k⊥ ) ∝ k⊥−4/3 ,510

Eϕ (k⊥ ) ∝ k⊥−10/3 ,(7)

T. Tatsuno et al., Gyrokinetic Simulation of Entropy Cascade in Two-Dimensional Electrostatic Turbulencelast one is for the solar wind turbulence for which the datafrom Cluster spacecraft at 19 earth radii is used [5]. In
the solar wind, the collision frequency is extremely small
because of the small density, leading to extremely large D.

Table 1 Estimated dimensionless number D in various physicalsystems. Ion scale turbulence at LAPD [28], DIII-D [29]and TFTR [30], ETG turbulence at NSTX [31], and solar
wind observation [5]. Eddy turnover time τ for the ETG
case is evaluated from simulation data [32].

ν [sec−1 ]2.1 · 1063300
120
10, 800
2.3 · 10−7

LAPD
DIII-D
TFTR
ETG @ NSTX
Solar wind

τ [sec]� 5.0 · 10−6
7.5 · 10−6
1.5 · 10−5
3.5 · 10−7190

3. Numerical simulation
3.1 Code: AstroGK

D
� 0.1
40
500
270
2.3 · 104

AstroGK is an Eulerian initial value solver for the GK
equations in 5D phase space in general, although we only
use 4 dimensions in this study. Namely, we discretize thevelocity space into grids as standard ﬂuid codes do in theposition space. It evolves the function
g=h−

where
Eh (k⊥ ) =Eϕ (k⊥ ) =� �|k⊥ |=k⊥T0|h k |2 dv,
2F0� q2 n0|ϕ k |2 .
2T0|k |=k
⊥

(8)
(9)

=
Note
that the total entropy (4) can be expressed as�[Eh (k⊥ ) − Eϕ (k⊥ )] dk⊥ .
At small position-space scales, collisions tend to dominate as velocity-space structure becomes ﬁner due to (6).
It can be shown [12, 25, 27] that the cutoﬀs are
(10)

where
D=

1
,
ντρ

(12)

according to (2) and (3). In order to make a 4D simulation using the 5D code, we used only 3 grid points along
the ﬁeld line (in z) with a very long box size and a periodic boundary condition. Starting from the homogeneous
initial condition along z, we conﬁrmed that no structuredeveloped in this direction.
The code uses a Fourier pseudo-spectral scheme forreal space dimensions perpendicular to the mean magnetic ﬁeld and Legendre collocation scheme in the twodimensional velocity space integration. The velocity-space
grid is taken with respect to particle energy ε = v2 andλ = v2⊥ /ε, which makes our grid points radially distributedin the v� -v⊥ plane [33]. For the collision operator we use aconservative ﬁnite diﬀerence scheme [34].
Time integration is made using the 3rd order AdamsBashforth scheme for the nonlinear term. The linearized
collision operator is treated by the ﬁrst-order implicit Euler
scheme with Sherman-Morrison formula for the momentsconserving corrections [34, 35].
AstroGK is written in Fortran 95 and is parallelizedusing MPI. Our biggest run (D2 in Table 2) cost about 36
wall-clock hours using 8,192 processor cores.

⊥

1
δv⊥c
∼ D−3/5 ,
∼
vth
k⊥c ρ

qF0 �ϕ� R
T0(11)

the quantity ν is the ion collision frequency and τρ is aturnover time at the scale ρ.
We introduced a dimensionless number D which represents the ratio between the outer scale, in this case thethermal Larmor radius, and the dissipation scale valid forthe kinetic theory. This is simply the ratio of the collision time and the eddy turnover time at the ion gyroscale.
It describes how many times eddies at the Larmor radiusscale turn over before the cumulative eﬀect of small-angle
collisions smears out the resulting structures in velocity
space at this scale. Analogous to the Reynolds number,
D becomes large when dissipation is weak, i.e., when thecollision frequency is small. We emphasize here that thisnumber is amplitude dependent since turnover time is determined by the E × B velocity where E is a ﬂuctuation.
Table 1 shows D values in more realistic and geometrically complicated physical systems. The ﬁrst three are
the ion-scale turbulence in the basic plasma device LAPD
[28] and in the fusion devices DIII-D [29] and TFTR [30].
The fourth one is the electron-temperature-gradient (ETG)
driven turbulence in NSTX, where ν is estimated from experimental conditions [31] while τ is from the streamerdominated results of simulation [32] at k⊥ ρe = 0.2. The

3.2

Initial condition

Here we describe the setup of the direct cascade simulations. We use the simple straight slab geometry in a
homogeneous background in the box size L x = Ly = 2πρ.
A series of decaying turbulence simulations were carriedout, with the initial condition put in |k x ρ|, |ky ρ| = 2 scale:
�
�
2y
2x
+ cos
+ χ(x, y) F0 .
(13)
ginit = g0 cos
ρρ

where F0 is a Maxwellian. Here, g0 is a constant andχ(x, y) is small-amplitude white noise superimposed on all
Fourier harmonics. From (3) and (12), we can calculate
ϕinit . The initial condition (13) is unstable to the kinetic
version of the Kelvin-Helmholtz (KH) instability, which
self-consistently evolves into turbulence. It is a completely
autonomous system and there is no entropy input during
the time evolution. We emphasize that there is no non-

511

T. Tatsuno et al., Gyrokinetic Simulation of Entropy Cascade in Two-Dimensional Electrostatic Turbulenceντinit9.3 · 10−3
5.6 · 10−3
1.9 · 10−3
1.9 · 10−3
1.9 · 10−3
7.4 · 10−4
7.4 · 10−4

D
32
48
118
118
118
440
440

k⊥c ρ
16
20
35
35
35
77
77

N x × Ny
642
642
1282
1282
1282
1282
2562

Nε × 2Nλ
322
322
162
642
1282
1282
1282

(dW/dt) / (Winit/τinit)

Run
A
B
C1
C2
C3
D1
D2

(W, E) / Winit

Table 2 Index of the runs described in Sec. 3. Dimensionless
number D is given for the initial stage of developed turbulence (t/τinit = 10), and cutoﬀ wave number is estimatedby (16) (see Sec. 3.3 and Fig. 3).

Maxwellian velocity-space structure in the initial condition. All ﬁne structure later observed in the velocity spaceis created from nonlinear interaction described in Sec. 2.
The results reported below were obtained in the seriesof runs indexed in Table 2, where N x × Ny is number ofcollocation points in the position space and Nε × 2Nλ is
the number of grid points in velocity space — the factor√of 2 corresponds to the sign of v� = ± ε(1 − λ). For eachcollision frequency, D is calculated using the turnover time
taken from the simulation (see Sec. 3.3 and Fig. 3). The
cutoﬀ scale is listed for each case here, and we determine
the grid resolution so that it resolves the smallest positionand velocity-space structures.

|ϕk| / W1/2
init

Time evolution of the collisionless invariants and theamplitude of several lower Fourier modes |ϕ k | are shownin Figs. 1–2, where τinit is the initial turnover time of the
gyroaveraged E × B motion deﬁned by
2πB0,
ck x ky ||�ϕinit � R ||with |k x ρ|, |ky ρ| = 2 and
� �
�1/2
1
2
|�ϕ�
.
||�ϕ� R || =R | F 0 dv dR
n0

0
-0.02
-0.04
-0.06
-0.08
-0.1
-0.12
-0.14
-0.16
-0.18
(b)
-0.2
0 2

D2
C2
B
A
W

6

8 10 12 14 16 18 20
t / τinit

D2
C2
B
A
4

6

8 10 12 14 16 18 20
t / τinit

Fig. 1 Time evolution of (a) and [(4) and (5)], and (b) rateof change of , normalized to initial . Labels A–D2
correspond to runs indexed in Table 2. Evolution ofdoes not diﬀer among runs signiﬁcantly, and is represented by the run D2.

3.3 Time evolution

τinit =

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
E
0.2
0.1
(a)
0
0 2 4

(14)

101
0
10
10-1
10-2
-3
10
-4
10
10-5
10-6
-7
10
10-8
-9
10
-10
10

kx=0; ky=1
kx=1; ky=0
kx=0; ky=2
kx=2; ky=0
0

2

4

6

8 10 12 14 16 18 20
t / τinit

Fig. 2 Time evolution of several lower Fourier modes ϕ k fromthe run D2 (See Table 2).

(15)

perfectly by collisional entropy production — i.e. entropy
balance [2] is satisﬁed. This is guaranteed under the controlled velocity-space numerical scheme [33]. In all runs
except for the one with the smallest D (run A in Table 2;
D = 32), the rate of change of is almost identical in the
nonlinear stage (t/τinit � 9). The slight horizontal shift ofthe lines is due to the diﬀerent timing of linear mode saturation because of the random-noise admixture in the initialconditions. This plot suggests that the asymptotic amountof dissipation is ﬁnite as collision frequency tends to zero,so the dissipation rate tends to a value asymptotically independent of the collisionality, but rather determined by
the nonlinear dynamics of the turbulence. Also, the fallingphase of d /dt corresponds to the development of the turbulence, and we may say that the turbulence has fully de-

In Fig. 1, the decrease of corresponds to the creation ofentropy since is proportional to the negative of the perturbed entropy. Initially, while instability grows (t/τinit �
8, see also Fig. 2), decays slowly with a rate proportional
to the collision frequency ν. Then as the instability saturates around t/τinit � 9, starts to decay very rapidly. This
rapid decay of suggests nonlinear creation of small scale
structures in velocity space as well as in position space.
Some time after that (t/τinit � 20), starts to decay exponentially suggesting the end of turbulence. The constancy
of implies the presence of an inverse cascade [27], which
will be discussed elsewhere.
Time evolution of the rate of change ofis shown
in Fig. 1(b). Note that the decay of is accounted for512

T. Tatsuno et al., Gyrokinetic Simulation of Entropy Cascade in Two-Dimensional Electrostatic Turbulence

1000

200

Eh(k⊥) / W

500

D

100

D2
C2
B

100
50
20
10

10

12

14
t / τinit

16

18

−4/3

-2

10

Eϕ(k⊥) / W

0

-1

10

k−4/3
⊥

-2

10

10-3

t/τinit= 9t/τinit=11t/τinit=14t/τinit=17t/τinit=20

10-4
10-5
-6

1

2

5

10
k⊥ρ

20

50

100
-1
10
10-2
10-3
-4
10
10-5
10-6
-7
10
-8
10

k⊥

B
C1
C2
D1
D2

10-3
10

20

10
Eh(k⊥) / Winit

10-1

-4

8

Fig. 3 Time evolution of the dimensionless number D [see (11)].
Labels B–D2 correspond to runs indexed in Table 2. Values of D in Table 2 is evaluated at t/τinit = 10.

10

(a)

1

2

5

10
k⊥ρ

20

50

100

(b)

B
C1
C2
D1
D2
1

2

−10/3

k⊥

5

10
k⊥ρ

20

50

100

Fig. 5 Time-averaged normalized wave-number (Fourier) spectra (a) Eh (k⊥ )/ and (b) Eϕ (k⊥ )/ [cf. (7)] for the runs
indexed in Table 2. Theoretically predicted slopes are
given for comparison.

100

Fig. 4 Time evolution of the wave-number spectra of Eh (k⊥ )from the run D2 (see Table 2).

overall steeper slope as it has not yet entered the fully developed turbulence stage. From t/τinit = 11 to 14 the spectra show a self-similar shape consistent with the inertialrange scaling (7). The spectra at t/τinit = 17 and 20 clearly
have an inertial range; the steeper falloﬀ at large k⊥ suggests a dissipation range. Note that the implied cutoﬀ scale
is consistent with the estimate from the time-evolving dimensionless number D [see Fig. 3 and (16)].
Figure 5 shows the wave-number spectra of Eh (k⊥ )and Eϕ (k⊥ ) [see (8) and (9)] of the asymptotic runs (B–D in Table 2) which are averaged over wave-number shelland normalized by the value of total at each time, andthen averaged over time 10 ≤ t/τinit ≤ 15. It is hard to
identify the dissipation range as the time average is takenover the period when the cutoﬀ is near the highest k⊥ end
so that inertial range can be taken as wide as possible. Unclear dissipation range in the time-averaged spectra is acceptable as exponential fall-oﬀ in this regime implies exponentially small dissipation apart from the cutoﬀ. Three
resolved cases B, C2 and D2 lie on top of one another andare consistent with the theoretical slopes (Eh ∝ k⊥−4/3 andEϕ ∝ k⊥−10/3 ). In contrast, the under-resolved runs in thevelocity space (C1) and in the position space (D1) haveshallower slopes than others. With the increase of D thedissipation cutoﬀ extends to smaller scales, thus higher resolution is required both in position and velocity space. The
simulation shows a consistent result with the scaling theory
including the dissipation cutoﬀ.

veloped by t/τinit = 10. In the following sections, we willconcentrate on the asymptotic cases, which have similar
rates of change of (runs B–D).
The dimensionless number D is deﬁned as follows.
As is seen in Fig. 2, the simple use of a single representative wave number k⊥ ρ = 2 yields signiﬁcant ﬂuctuationof D which is not necessarily physical. On the other hand,
the k⊥ ρ = 1 component is created by the KH instability
as a part of the inverse cascade, and is also irrelevant forcharacterizing direct cascade. Thus we introduced a definition of the turnover time at scale � = πρ, denoted byτπρ , which uses all k⊥ ρ ≥ 2 components of the ﬂuctuating
potential. Using τπρ , the time evolution of D is shown in
Fig. 3. The values are shown for three resolved cases fromt/τinit ≥ 8 since the turbulent spectra start to develop after
that. At t/τinit = 10 the dimensionless number D takes values shown for each run in Table 2. The cutoﬀ wave number
is then obtained fromk⊥c ρ = αD3/5 .

(16)

The value α = 2 corresponds to the particular initial condition and domain size used here.

3.4 Wave-number spectra

We ﬁrst show the time evolution of the wave-number
spectra of Eh (k⊥ ) averaged over wave-number shell [deﬁned by (8)] in Fig. 4. At t/τinit = 9 the spectrum shows an513

T. Tatsuno et al., Gyrokinetic Simulation of Entropy Cascade in Two-Dimensional Electrostatic Turbulence

log10 [E^g(k⊥, p) / W]

100

20

-3

10

-4

5

-5

2

-6

1

-7
1

2

5

10 20
k⊥ρ

E^g(k⊥,p) / W

-2

10

-3

10

-4

10

50 100Fig. 6 Time-averaged normalized 2D spectrum Êg (k⊥ , p) [deﬁned by (17)] taken from the run D2 (See Table 2).

In order to characterize the velocity-space structure
quantitatively, Plunk et al. [27] proposed a dual space variable p to the velocity variable v⊥ and introduced a 2D spectrum in (k⊥ , p) space:�Êg (k⊥ , p) =p|ĝ k (p)|2 ,
(17)

k−1/3
⊥

pvth=20k−2⊥pvth=30
pvth=40
pvth=50125k⊥ρ=10

p

-2

10

-3

10

-4

10
k⊥ρ

20

50

100

(b)

−1/3

k⊥ρ=20k⊥ρ=30k⊥ρ=40p−2k⊥ρ=50

10-5
10-6

3.5 Velocity-space spectra(a)pvth=1010-5
10-6E^g(k⊥,p) / W

pvth

50

10-2

-1

1

2

5

10
pvth

20

50

100

Fig. 7 Averaged asymptotic spectra from the run D2 (See Table2) (a) as a function of k⊥ for various ﬁxed p and (b) as afunction of p for various ﬁxed k⊥ .
100

|k⊥ |=k⊥E^g(p) / W�
where ĝ k (p) = J0 (pv⊥ )g k (v) dv is the Hankel transform.
Figure 6 shows the time average of the normalized
2D spectrum Êg (k⊥ , p) taken from the run D2 (see Table2). We ﬁrst note that there is a gap in the low-k⊥ , highp region. This gap appears because high-p component is
not easily created with a low-k⊥ mode, while the high-k⊥
component can be created even for low-p as Navier-Stokes
turbulence does with Maxwellian velocity proﬁle. This is
likely due to the fact that nonlinear phase mixing acts morestrongly at sub-Larmor scales. The spectral contour develops along the diagonal pvth = k⊥ ρ, which reﬂects the factthat the correlation of position- and velocity-space structure follows our conjecture (6) very well. Thus, comparable resolution in velocity and position space is required,especially when we investigate smaller scales than the Larmor radius.
In Ref. [27] the asymptotic spectra are derived fork⊥ ρ � pvth and k⊥ ρ � pvth limits, namely⎧⎪⎪⎨k⊥−2 p−1/3 (k⊥ ρ � pvth ).
(18)Êg (k⊥ , p) ∝ ⎪⎪⎩ p−2 k−1/3 (k⊥ ρ � pvth )

10

-1

10

-2

10

-3

10-4

p−4/3
B
C2
C3
D2
1

2

5

10
pvth

20

50

100

Fig. 8 Time-averaged� normalized velocity-space (Hankel) spectrum Êg (p) = Êg (k⊥ , p) dk⊥ [see also (17)] for the runsin Table 2.

is steeper than the prediction (18) for the k⊥ ρ = 10 casebecause of the remnant of the gap in Fig. 6. The agreement with the asymptotic expectation becomes good fork⊥ ρ � 30 even though we don’t have a wide asymptotic
range in this region. The overall spectra show a reasonable symmetry in k⊥ and p except for the structure comingfrom the gap, which indicates consistency with the theoretical prediction (18).
Figure 8 shows� the time-averaged spectra in the Hankel space Êg (p) = Êg (k⊥ , p) dk⊥ normalized to at eachtime. The theoretical expectation is [27]⊥Figure 7 shows various (a) horizontal and (b) vertical slices
of Fig. 6. Depending on the value of p in Fig. 7(a), we observe a slope approaching k⊥−2 in the high-k⊥ regime. On
the other hand, the k⊥−1/3 slope for k⊥ ρ � pvth seems diﬃcult to obtain, reﬂecting the gap in the low-k⊥ , high-p region seen in Fig. 6. Figure 7(b) shows a rather good slope
consistent with that predicted for the low-p limit ( ∝ p−1/3 )even for the smallest k⊥ in the ﬁgure. The high-p slopeÊg (p) ∝ p−4/3 ,(19)which logically follows from the ﬁrst spectrum in (7)and pvth ∼ k⊥ ρ [equivalent to (6)]. The numerical result again shows approximate consistency with the theo-514

T. Tatsuno et al., Gyrokinetic Simulation of Entropy Cascade in Two-Dimensional Electrostatic Turbulence100

retical prediction and conﬁrms that small-scale structure
is formed in the velocity space. The large hump in the
low-p regime (p � 3) is due to long-wave-length modes,
which have signiﬁcantly larger amplitudes than the rest.
These long-wave-length modes have velocity-space structure
whose Hankel transform yields
� close to 2Maxwellian,22 2
J0 (pv⊥ ) e−v⊥ /vth v⊥ dv⊥ /v2th = e−p vth /4 /2. The gradualsteepening at the high-p region may be related to the gap
in Fig. 6 (Notice that Fig. 8 corresponds to taking the horizontal sum of Fig. 6). The wiggles in the high-p end of theruns C2 and D2 come from the slight lack of resolution,
which goes away by increasing the velocity resolution as
the high velocity-resolution run C3 shows. This increaseof the velocity-space resolution does not aﬀect the wavenumber spectra, while failure to resolve the scaling regime
does (see line C1 in Fig. 5).

0.002
0.0015

50

0.001

Kρ

20

0.00050-0.000510
5-0.001-0.0015-0.0022
1
1

2

5 10 20Qρ

50 100Fig. 9 Time-averaged normalized entropy transfer function
T F (K, Q) [see also (22)] for the run D2 in Table 2.

3.6 Entropy transfer

Pvth

In the scaling theory (Sec. 2.2) we assumed local-scaleinteraction by following Kolmogorov’s argument. The applicability of such an assumption may be directly investigated in the numerical simulation.
In order to make this diagnostic we have introduced a
Fourier ﬁltered function [36] deﬁned by�gK (R, v⊥ ) :=g k (v⊥ )eik·R ,
(20)

50

0.004
0.003

20

0.002
0.001

10

0
-0.001-0.002-0.003-0.004

5
2
1
1

k∈K

where K = {� : Kρ − ∞/∈ ≤ |�|ρ < Kρ + ∞/∈}. The
ﬁltering is orthogonal, and gK denotes the component of
the distribution function with scale K −1 . Then the entropy
amount contained in this scale obeys the evolution equation��d|gK |2dR dv =T F (K, Q)−collisions, (21)
dt
2F0
Q

2

5 10
Svth

20

50

Fig. 10 Time-averaged normalized entropy transfer function
T H (P, S ) [see also (24)] for the run C2 in Table 2.

where P = {p : Pvth − 1/2 ≤ pvth < Pvth + 1/2}. Then,similarly to the Fourier ﬁltering, entropy transfer functionin Hankel space is described by�
T H (P, S ) := −ǧP v E×B · ∇ǧS v⊥ dv⊥ dR,(24)

where T F (K, Q) denotes the nonlinear transfer of entropy
from scale Q−1 to K −1�
gK v E×B · ∇gQ
T F (K, Q) := −v⊥
dR dv⊥ , (22)F0√where ǧ(v⊥ ) := g(v⊥ )/ F0 is introduced to make entropy
a bilinear form, and T H (P, S ) = −T H (S , P).
Figure 10 shows the normalized, time-averaged transfer function in the velocity space. For numerical purposes
the integral in (23) is approximated by a single representative Hankel mode. Nevertheless the transfer function is localized along the diagonal very well, which, together with
Fig. 9, may suggest the entropy transfer along the diagonalin (k⊥ , p)-plane (see Fig. 6).

which is by deﬁnition symmetric with respect to the exchange of two arguments.
The numerical result obtained from the run D2 (seeTable 2) is shown in Fig. 9, which is normalized by ateach time and averaged over 10 ≤ t/τinit ≤ 15. Figure
9 shows a remarkable locality of the interaction, which isclearly seen even at each time. The color denotes the direction of the entropy transfer, which indicates that the entropy is transferred from large scales to small scales.
The scaling theory does not require locality in velocity
space scales, however the velocity-space transfer function
may still give some useful information. It may be monitored by Hankel ﬁltered function deﬁned by�gP (v⊥ ) :=ĝ k (p)J0 (pv⊥ )p d p,
(23)

4. Summary

We presented electrostatic, decaying turbulence simulations for weakly collisional, magnetized plasmas using
the gyrokinetic model in 4D phase space (two positionspace and two velocity-space dimensions; the extension to
three spatial dimensions is left for future work). Landau
damping was removed from the system by ignoring variation along the background magnetic ﬁeld.
Nonlinear interactions introduce an amplitude-

p∈P

515

T. Tatsuno et al., Gyrokinetic Simulation of Entropy Cascade in Two-Dimensional Electrostatic Turbulence

dependent perpendicular phase mixing of the gyrophaseindependent part of the perturbed distribution function
and create structure in v⊥ which is ﬁner for higher k⊥ (seeFig. 6). We found that the wave-number (Fourier) spectraof the perturbed distribution function and the resulting
electrostatic ﬂuctuations at sub-Larmor scales agreed wellwith theoretical predictions based on the interpretationof the nonlinear phase mixing as a cascade of entropy
in phase space (see Figs. 5 and 7) [12, 25, 27].
The
velocity-space (Hankel) spectra show a rough consistency
with the theoretical scaling, although the agreement is not
as good as that of the wave-number spectra (see Figs. 8
and 6).
We introduced a dimensionless number (analogous toReynolds number) that characterized the scale separation
between the thermal Larmor scale and the collisional cutoﬀ in phase space [see (11)], and showed that this number
correctly predicted the resolution requirements for nonlinear gyrokinetic simulations (see Table 2). We also showed
the trend that the entropy generation rate (or irreversibility)
is independent of the collision frequency in the asymptotic
limit of the weak collisionality (see Fig. 1). Finally we
presented diagnostics of nonlinear transfer functions andshowed that local-scale cascade of entropy is supported
very well in both Fourier and Hankel spaces.
We note that there are, in general, entropy cascades
for each plasma species. Equations for the gyrokinetic turbulence at and below the electron Larmor scale are mathematically similar to the model simulated here and identicalarguments apply [25, 27]. Similar considerations are also
possible for ion-scale electromagnetic turbulence [25] andfor minority ion species (with some diﬀerences to be discussed elsewhere).
The structure of the small scales in phase space that
we have presented is likely to be a universal feature ofmagnetized plasma turbulence. Understanding it theoretically and diagnosing it numerically is akin to the inertialrange studies for Kolmogorov turbulence, extended to thekinetic phase space. One should expect rich and interesting physics to emerge and it is likely that, just like in thecase of ﬂuid turbulence, predicting large-scale dynamics
will require eﬀective models for the small-scale cascade.
An immediate physical implication of the existence of theentropy cascade is a turbulent heating rate independent ofcollisionality in weakly collisional plasmas.

[3] Y. Idomura, Phys. Plasmas 13, 080701 (2006).
[4] J. Candy, R. E. Waltz, M. R. Fahey et al., Plasma Phys.
Control. Fusion 49, 1209 (2007).
[5] S. D. Bale, P. J. Kellogg, F. S. Mozer et al., Phys. Rev. Lett.
94, 215002 (2005).
[6] G. G. Howes, W. Dorland, S. C. Cowley et al., Phys. Rev.
Lett. 100, 065004 (2008).
[7] S. P. Gary, S. Saito, and H. Li, Geophys. Res. Lett., 35,
L02104 (2008); S. Saito, S. P. Gary, H. Li et al., Phys. Plasmas 15, 102305 (2008).
[8] F. Sahraoui, M. L. Goldstein, P. Robert et al., Phys. Rev.
Lett. 102, 231102 (2009).
[9] O. Alexandrova, J. Saur, C. Lacombe et al., Phys. Rev. Lett.
103, 165003 (2009).
[10] L. Boltzmann, Sitsungsber. Akad. Wiss. 66, 275 (1872) [inKinetic Theory 2 (Pergamon, Oxford, 1966) 88].
[11] P. Helander and D. J. Sigmar, Collisional Transport in
Magnetized Plasmas (Cambridge Univ., Cambridge, 2002).
[12] T. Tatsuno, W. Dorland, A. A. Schekochihin et al., Phys.
Rev. Lett. 103, 015003 (2009).
[13] R. Fjørtoft, Tellus 2, 225 (1953).
[14] R. H. Kraichnan, Phys. Fluids 10, 1417 (1967).
[15] A. Hasegawa and K. Mima, Phys. Fluids 21, 87 (1978).
[16] W. Horton, P. Zhu, G. T. Hoang et al., Phys. Plasmas 7,1494 (2000).
[17] P. Catto, Plasma Phys. 20, 719 (1978).
[18] T. M. Antonsen and B. Lane, Phys. Fluids 23, 1205 (1980).
[19] E. A. Frieman and L. Chen, Phys. Fluids 25, 502 (1982).
[20] G. G. Howes, S. C. Cowley, W. Dorland et al., Astrophys.
J. 651, 590 (2006).
[21] I. G. Abel, M. Barnes, S. C. Cowley et al., Phys. Plasmas,
15, 122509 (2008).
[22] W. Dorland and G. W. Hammett, Phys. Fluids B 5, 812(1993).
[23] J. A. Krommes, Phys. Plasmas 6, 1477 (1999).
[24] K. Hallatschek, Phys. Rev. Lett. 93, 125001 (2004);
B. D. Scott, E-print arXiv:0710.4899 (2007).
[25] A. A. Schekochihin, S. C. Cowley, W. Dorland et al., Astrophys. J. Suppl. Ser. 182, 310 (2009); A. A. Schekochihin, S. C. Cowley, W. Dorland et al., Plasma Phys. Control.
Fusion 50, 124024 (2008).
[26] K. Gustafson, D. del-Castillo-Negrete, and W. Dorland,
Phys. Plasmas, 15, 102309 (2008).
[27] G. G. Plunk, S. C. Cowley, A. A. Schekochihin et al., J.
Fluid Mech., submitted (2009); E-print arXiv:0904.0243.
[28] T. A. Carter and J. E. Maggs, Phys. Plasmas 16, 012304(2009).
[29] M. Jakubowski, R. J. Fonck, and G. R. McKee, Phys. Rev.
Lett. 89, 265003 (2002).
[30] R. D. Durst, R. J. Fonck, J. S. Kim et al., Phys. Rev. Lett.
71, 3135 (1993).
[31] E. Mazzucato, D. R. Smith, R. E. Bell et al., Phys. Rev.
Lett. 101, 075001 (2008).
[32] F. Jenko and W. Dorland, Phys. Rev. Lett. 89, 225001(2002).
[33] M. Barnes, W. Dorland, and T. Tatsuno, to be published inPhys. Plasmas (2010); E-print arXiv:0907.4413.
[34] M. Barnes, I. G. Abel, W. Dorland et al., Phys. Plasmas 16,
072107 (2009).
[35] W. H. Press, S. A. Teukolsky, W. T. Vetterling et al., Numerical Recipes (Cambridge Univ., Cambridge, 1988).
[36] A. Alexakis, P. D. Mininni, and A. Pouquet, Phys. Rev. E
72, 046301 (2005).

Acknowledgments

TT, WD, RN and GGP thank the Leverhulme Trust
Network for Magnetised Plasma Turbulence for travel support. AAS was supported by STFC. Numerical computations were performed at NERSC, NCSA, and TACC.
[1] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR 30, 299(1941); ibid. 32, 16 (1941) [Proc. Roy. Soc. London A 434,
9 (1991); ibid. 434, 15 (1991)].
[2] T.-H. Watanabe and H. Sugama, Phys. Plasmas 11, 1476
(2004).

516