PRL 103, 015003 (2009)PHYSICAL REVIEW LETTERS

week ending
3 JULY 2009

Nonlinear Phase Mixing and Phase-Space Cascade of Entropy in Gyrokinetic Plasma Turbulence
T. Tatsuno,1 W. Dorland,1 A. A. Schekochihin,2 G. G. Plunk,1 M. Barnes,1,2,3 S. C. Cowley,3 and G. G. Howes4
1

Department of Physics, IREAP and CSCAMM, University of Maryland, College Park, Maryland 20742, USA
2
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3NP, United Kingdom
3
Euratom/UKAEA Fusion Association, Culham Science Centre, Abingdon OX14 3DB, United Kingdom
4Department of Physics and Astronomy, University of Iowa, Iowa City, Iowa 52242, USA(Received 16 November 2008; published 30 June 2009)
Electrostatic turbulence in weakly collisional, magnetized plasma can be interpreted as a cascade ofentropy in phase space, which is proposed as a universal mechanism for dissipation of energy inmagnetized plasma turbulence. When the nonlinear decorrelation time at the scale of the thermal
Larmor radius is shorter than the collision time, a broad spectrum of fluctuations at sub-Larmor scales
is numerically found in velocity and position space, with theoretically predicted scalings. The results areimportant because they identify what is probably a universal Kolmogorov-like regime for kinetic
turbulence; and because any physical process that produces fluctuations of the gyrophase-independent
part of the distribution function may, via the entropy cascade, result in turbulent heating at a rate that
increases with the fluctuation amplitude, but is independent of the collision frequency.
DOI: 10.1103/PhysRevLett.103.015003

PACS numbers: 52.30.Gz, 52.35.Ra, 52.65.Tt

Introduction.—Turbulence is inherently nonlinear anddynamically complicated. In the general case, a broad
spectrum of fluctuations is excited, in both wave number
and frequency. For turbulent, magnetized plasma, the equations of magnetohydrodynamics provide a pedagogicallyrich description of the dynamics. However, for those turbulent eddies whose parallel wavelengths (relative to themagnetic field) are comparable to or smaller than the collisional mean free path and whose perpendicular wavelengths are comparable to or smaller than the Larmor
radius of one of the constituent species of the plasma,magnetohydrodynamic theory breaks down. In such cases,
the gyrokinetic (GK) theory [1,2] represents a rigorous
limit of plasma kinetics for anisotropic (kk  k? ), lowfrequency (!  , the ion cyclotron frequency) fluctuations. In this Letter, we present a GK description of turbulence in a simplified situation, chosen to isolate a novel
phenomenon which is a generic component of all GK
turbulence: the simultaneous cascade of entropy to smallerscales in both real space and velocity space. This phasespace cascade is the mechanism by which turbulent energy
associated with fluctuating fields is brought to small scalesin velocity space, where even very infrequent collisions are
sufficient to provide irreversibility and thus heating.
Below, we present the theory and first-principles simulations of the phase-space cascade in a homogeneous, electrostatic, magnetized plasma.
It is well known that Landau and Barnes damping ofelectromagnetic plasma fluctuations lead to the generationof small-scale structures in fðvk Þ, where f is the oneparticle distribution function, and vk is the velocity coordinate along the background magnetic field [3,4]. This isassociated with the free-streaming of particles along thefield. As t increases, a single Fourier harmonic of the
distribution function fkk  eikk vk t gets progressively more
0031-9007=09=103(1)=015003(4)oscillatory in vk -space. Eventually, even infrequent collisions are sufficient to smooth these oscillatory features,
since the collision operator is roughly a diffusion operator
in velocity space. As long as collisions are sufficiently
infrequent, the damping rate depends not on the collisionrate, but on the nature of the wave and its phase velocity
relative to the thermal speeds of the plasma species.
Physically, Landau damping is the smearing of spatial
perturbations that occurs when there is a spread in the
distribution of parallel velocities. We recall for futurereference that this generation of velocity-space structure
is independent of the fluctuation amplitudes.
Besides this linear parallel phase mixing, there exists a
nonlinear phase mixing process [5] that, in a strongly
turbulent plasma and at spatial scales smaller than the
Larmor radius, drives the formation of structure in fðv? Þmuch more rapidly than parallel phase mixing drives fðvk Þ.
Physically, this nonlinear phase mixing is the smearing of
spatial perturbations due to the spread in the distribution of
gyroaveraged E  B velocities (see Fig. 1). Unlike for theparallel phase mixing, the rate of generation of v-spacestructure by this process is proportional to the fluctuation
amplitude. In this Letter, we present a study of this nonlinear process, which we interpret as a turbulent cascade ofentropy in phase space [6]. As such, it represents a conceptually novel nonlinear phenomenon, where generationof small scales in the position and velocity space occurs in
an intertwined way. This process, which is likely to be afundamental and ubiquitous feature of magnetized plasma
turbulence, has never been numerically diagnosed andanalyzed before, although Krommes [3] did point out thegeneral possibility of the coupling between position andvelocity space.
Gyrokinetics in 2D.—Let the distribution function be
f ¼ F0 þ f, where F0 is a Maxwellian with density n0015003-1

Ó 2009 The American Physical Society

FIG. 1 (color online). Schematic view of the nonlinear phasemixing superimposed on the potential from the Run (iii) (seeTable I) at t=init ¼ 10 and the largest wavelength mode takenout. When the fluctuation scale ‘ & , the gyroaverage of theelectric field induces a decorrelation of the distribution function
at the velocity-space scale corresponding to the difference inLarmor radii ‘v ¼ v=  ‘ [see (7)].

and temperature T0 , and f ¼ h  q’F0 =T0 , where q is
the particle charge and ’ is the electrostatic potential. To
keep the focus on the nonlinear process, we consider
electrostatic GK turbulence in slab geometry with kk ¼
0. Then the non-Boltzmann part h of the perturbed ion
distribution function satisfies [1]invariant W is proportional toRminus the perturbed partof the entropy of the system,  f lnfdrdv [3,4]. Here wewill refer to W as ‘‘entropy’’ to emphasize this connection.
The second invariant E is conserved in the 2D electrostaticcase only.
Scalings.—A scaling theory of the entropy cascade in
the sub-Larmor scale range can be developed in a way
reminiscent of the Kolmogorov-style turbulence theories[6]. Assume that at (perpendicular) scales ‘  , thetransfer of entropy is local in scale. On dimensional
grounds, the entropy flux isv2th hv3th 2¼ const(5)
‘ n 0
until it reaches the collisional dissipation scale, where vth
is the thermal speed and ‘ is the nonlinear decorrelation
time at scale ‘. The neglect of the ’2 term in W [see (3)] is
justified post hoc due to its smallness in the ‘   regime[see (9) and Fig. 3(a)]. There is a self-consistent electrostatic potential at the scale ‘: from (2),q’  ‘ 1=2 v? 1=2 hv3th ‘
:
(6)


T0

vth
n0 
Here we have assumed that the nonlinear phase mixing
produces velocity-space structures correlated with the spatial scale via (see Fig. 1 and Refs. [6,8])

qF @h’iR@h cz^  rh’iRþ; (1) rh ¼ hC½hiR þ 0
B0
T0 @t
@t
where B0 is the background magnetic field aligned with thez axis and hiR is the gyroaverage holding the guiding
center position R constant. The collision operator C½hused in our simulations contains pitch-angle scatteringand energy diffusion with proper conservation properties[7]. The quasineutrality condition yields
Z
Z k v X
Q’ ¼ q hhir dv ¼ q eikr J0 ? ? hk dv; (2)kwhere hir denotes the gyroaverage at fixed
P particle position r, J0 is the Bessel function, Q ¼ s q2s n0s =T0s for
Boltzmann-response (3D) electrons or Q ¼ q2i n0i =T0i forno-response (2D) electrons, and s and i are the species
indices. Our results are not affected by the choice of theelectron response. For concreteness, we henceforth use noresponse electrons since electrons cannot contribute to thepotential if kk ¼ 0 exactly. In the absence of collisions, thesystem has two positive definite conserved integrals [6,8]:

ZZ T0 f2
Z Z T0 hh2 ir
Q
W¼
drdv ¼
dv ’2 dr; (3)
2
2F0
2F0
E ¼

week ending
3 JULY 2009PHYSICAL REVIEW LETTERS

PRL 103, 015003 (2009)QX
ð1  0 Þj’k j2 ;
2 k
2

(4)

v?
‘
 :
vth


(7)

This has allowed us to estimate the velocity integral in (2)as a random-walk-like accumulation of the integrand represented by the product of hk , which is a random functionof v? whose ‘‘step size’’ is given by (7) with ‘  k1? , andof the Bessel function, which introduces a reduction factorof ð‘=Þ1=2 .
The decorrelation time ‘ may be estimated by balancing the @t term with the nonlinear term in (1), leading to 1=2
‘2

‘2

:
(8)
‘ 
‘
ch’iR =B0c’=B0
Substituting (6) into (8) and (8) into (5) yields h  ‘1=6 and’  ‘7=6 . Therefore, the spectra of h and ’ are;E’ ðk? Þ  k10=3;(9)Eh ðk? Þ  k4=3?
?
RP
where Eh ðk? Þ ¼ jk? j¼k? T0 jhk j2 =2F0 dv and E’ ðk? Þ ¼
P2
2
0 . Note that the total entropy (3) can
jk? j¼k? q n0 j’k j =2TR
be expressed as W ¼ ½Eh ðk? Þ  E’ ðk? Þdk? .
Dissipation cutoff.—From the balance between the nonlinear decorrelation time (8) and the collision time 1 ,one obtains an estimate of the dissipation cutoff scales inboth the velocity space and the real space [see (7)]. Using
C½h  v2th h=v2? , we find the cutoffs

2

where 0 ¼ I0 ðk2? 2 =2Þek?  =2 , I0 is the modified Bessel
function and  is the ion thermal Larmor radius. The
015003-2

v?c1
 D3=5 ;

vth
k?c 

D¼

1
;


(10)

PRL 103, 015003 (2009)

week ending
3 JULY 2009PHYSICAL REVIEW LETTERS

where  is the nonlinear decorrelation time measured at‘ ¼ . We have introduced a new dimensionless number D
to characterize the scale separation in gyrokinetic turbulence: analogous to the Reynolds number in fluid turbulence, large D corresponds to a broader scaling range overwhich the entropy cascade extends, and to dissipation at
smaller scales. Here, however, the smallest spatial scaleobserved is determined by the v-space scale for which
diffusion in velocities becomes important, through thecorrelation between real and velocity space given by (7).
The fact that D increases with the amplitude of the fluctuations at the Larmor scale clearly distinguishes this processfrom linear Landau damping. We note that for 3D gyrokinetic turbulence, the nonlinear phase mixing is a much
faster process than the linear one if the fluctuation amplitude is sufficiently large.
Numerical simulations.—We now report the first-of-akind numerical investigation of the entropy cascade inphase space, carried out with the GK code AstroGK. The
code uses a Fourier pseudospectral scheme for the realspace dimensions perpendicular to the background magnetic field and a Legendre collocation scheme for thevelocity-space integrations. The velocity space is discretized in energy " ¼ v2 and  ¼ v2? =". In the absence ofcollisions, AstroGK conserves the invariants (3) and (4)with a high precision.
The results reported below were obtained in three runs atdecreasing collision frequency  and correspondingly increasing spatial and velocity resolution. They are indexedin Table I, where Nx  Ny is number of collocation points
in the real space and N"  2N is the number of grid points
in velocity space—the factor of 2 corresponds to the sign ofpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
vk ¼  "ð1  Þ. Our highest-resolved run required 36
wallclock hours on 8192 processors.
The code evolves g ¼ h  qF0 h’iR =T0 and ’ via
Eqs. (1) and (2). We take the box size Lx ¼ Ly ¼ 2
and start from the initial condition ginit ¼ g0 ½cosð2x=Þ þcosð2y=Þ þ ðx; yÞF0 , where g0 is a constant and ðx; yÞ
is a small-amplitude white noise superimposed on all
Fourier modes. From (2), we can calculate ’init .
Time evolution.—The initial jkx j, jky j ¼ 2 configuration is unstable: the amplitudes of ’ corresponding tojkx j, jky j ¼ 1 grow and then saturate around t=init ’
9, where init ¼ 2B0 =ðck2? k h’init iR kÞ is the turnovertime associated with the initial condition and k h’iR k¼
RR½ð1=n0 Þ jh’iR j2 F0 dvdR1=2 . The nonlinear interactions
between modes produce smaller scales down to a cutoff
determined by D [see (10)]. The turbulent spectra fill up byt=init ’ 10, then decay with time.
The time evolution of the collisionless conserved quantities W and E [see (3) and (4)] is shown in Fig. 2. During
the initial growth of the instability of the jkx j, jky j ¼ 1,
W decays very slowly at a rate , consistent with acollisional decay rate associated with the large-scale
phase-space variation of ginit . Once turbulence develops,

TABLE I. Index of the runs.
Run

Nx  NyN"  2Ninit

D

k?c 

(i)
(ii)
(iii)

642
1282
2562

322
642
1282

5:6  103
1:9  103
7:4  104

48
118
440

20
35
77

W decays more rapidly as the entropy cascade transfers it
nonlinearly to smaller scales in phase space, until thefluctuations of the distribution function are thermalized(dissipated) at the collisional cutoff.
The decrease of W from its initial value corresponds tothe amount of entropy (heat) production due to the irreversible collisional smearing of the distribution function.
The turbulence that follows the initial instability enhances
the heating, suggesting that small-scale velocity-spacestructure is generated (this is confirmed below). As expected, the rate of dissipation is not strongly affected by thecollision frequency, i.e., there is a finite amount of dissipation even as the collision frequency tends to zero. The
dissipation rate is determined instead by the nonlinear
cascade rate.
While W decays, E stays almost constant. If we increase
the size of the simulation box, the jkx j, jky j ¼ 1 modes
themselves become unstable to even longer-wavelength
modes. We attribute both this instability and the failureof E to decay to the intrinsic tendency of E to have an
inverse cascade [8,9], which we do not discuss here.
Spectra and scalings.—The wave-number spectra of the
decaying developed turbulence are given in Fig. 3(a). They
are angle integrated over wave-number shells jk? j ¼ k? ,normalized by WðtÞ at each time and then averaged overtime for 10  t=init  15. As resolution is increased, thespectra appear to converge to the theoretically predictedscalings (9), which supports the validity of our dimensional
and physical considerations of the entropy cascade.
To characterize the entropy cascade in the velocity
space, Plunk et al. [8] introduced velocity-space spectra

FIG. 2 (color online). Time evolution of W and E [Eqs. (3)and (4)] normalized to initial W. The runs (i)–(iii) are indexed inTable I. Evolution of E does not differ among runs significantly,
and is given for run (iii).

015003-3

PRL 103, 015003 (2009)PHYSICAL REVIEW LETTERS

week ending
3 JULY 2009

FIG. 3 (color online). Time-averaged
normalized (a) wave-number (Fourier)spectra Eh ðk? Þ=W and E’ ðk? Þ=W[cf. (9)], and (b) velocity-space (Hankel)
spectrum E^ g ðpÞ=W for the runs indexedin Table I. Theoretically predicted slopes
are given for comparison.

R
PE^ g ðpÞ ¼ k pjg^ k ðpÞj2 , where g^ k ðpÞ ¼ J0 ðpv? Þgk ðvÞdv
is a Hankel transform. The theoretical expectation is thatE^ g ðpÞ  p4=3 because the real- and velocity-space scales
should be related according to (7), which, in terms of thedual variable p, becomes k?   pvth . The time-averaged
Hankel spectrum E^ g ðpÞ obtained in our simulations isshown in Fig. 3(b). This again shows approximate consistency with the theoretical prediction and confirms that
small-scale structure is formed in the velocity space.
Dissipation cutoff.—In Table I, we show for each of ourruns the dimensionless number D ¼ ð Þ1 , where  ¼
2B0 =ðck2? k h’0 iR kÞ measured at t=init ¼ 10 and ’0 is’ with the jkx j, jky j ¼ 1 modes taken out (see also
Fig. 1). Also shown is the theoretical estimate (10) for
the wave-number cutoff k?c  ¼ D3=5 , where  ¼ 2 is anempirical value that corresponds to our particular set up.
Comparing with the wave-number and velocity-space
spectra in Fig. 3, we see that (10) describes the resolution
requirements quite well. With fewer velocity grid points,
we find shallower wave-number spectra than the resolved
ones, while with more, we resolve below the velocity cutoff
without any change in the wave-number spectra. Thus D isa good indicator of necessary and sufficient resolution in
full 4D phase space.
Conclusions.—We have presented electrostatic, decaying turbulence simulations for weakly collisional,
magnetized plasmas using the gyrokinetic model in 4D
phase space (two real-space and two velocity-spacedimensions). Landau damping was removed from the system by ignoring variation along the background magnetic
field. Nonlinear interactions introduce an amplitudedependent perpendicular phase mixing of the gyrophaseindependent part of the perturbed distribution function andcreate structure in v? which is finer for higher k? . We havefound that the wave-number (Fourier) and velocity-space(Hankel) spectra of the perturbed distribution function andthe resulting electrostatic fluctuations at sub-Larmor scales
agree well with theoretical predictions based on the interpretation of the nonlinear phase mixing as a cascade ofentropy in phase space [6,8]. We have introduced a dimensionless number D (analogous to Reynolds number) thatcharacterizes the scale separation between the thermal
Larmor scale and the collisional cutoff in phase space [see(10)], and showed that this number correctly predicts theresolution requirements for our simulations.
We note that there are, in general, entropy cascades foreach 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 [6,8]. Similar considerations are also
possible for ion-scale electromagnetic turbulence [6] andfor minority species.
The small-scale phase-space structure that we have discovered is likely to be a universal feature of strong, magnetized plasma turbulence. Understanding it theoreticallyand diagnosing it numerically is akin to the inertial-range
studies for Kolmogorov turbulence, extended to the kinetic
phase space. One should expect rich and interesting physics to emerge and it is likely that predicting large-scale
dynamics will require effective models for the small-scalecascade. An immediate key physical implication of theexistence of the entropy cascade is a turbulent heating
rate independent of collisionality in weakly collisional
plasmas.
Numerical computations were performed at NERSC,NCSA and TACC. This work was supported by the U.S.
DOE Center for Multiscale Plasma Dynamics, STFC, andthe Leverhulme Trust Network for Magnetized Plasma
Turbulence.

[1] E. A. Frieman and L. Chen, Phys. Fluids 25, 502 (1982).
[2] G. G. Howes et al., Astrophys. J. 651, 590 (2006).
[3] J. A. Krommes, Phys. Plasmas 6, 1477 (1999); J. A.
Krommes and G. Hu, Phys. Plasmas 1, 3211 (1994).
[4] T.-H. Watanabe and H. Sugama, Phys. Plasmas 11, 1476
(2004).
[5] W. Dorland and G. W. Hammett, Phys. Fluids B 5, 812(1993).
[6] A. A. Schekochihin et al., Astrophys. J. Suppl. Ser. 182,
310 (2009); A. A. Schekochihin et al., Plasma Phys.
Controlled Fusion 50, 124024 (2008).
[7] I. G. Abel et al., Phys. Plasmas 15, 122509 (2008); M.
Barnes et al., Phys. Plasmas (to be published).
[8] G. G. Plunk et al., J. Fluid Mech. (to be published);
arXiv:0904.0243.
[9] Y. Idomura, Phys. Plasmas 13, 080701 (2006); J. Candyet al., Plasma Phys. Controlled Fusion 49, 1209 (2007).

015003-4