PRL 100, 065004 (2008)

PHYSICAL REVIEW LETTERS

week ending
15 FEBRUARY 2008

Kinetic Simulations of Magnetized Turbulence in Astrophysical Plasmas
G. G. Howes,1,* W. Dorland,2 S. C. Cowley,3,4 G. W. Hammett,5 E. Quataert,1 A. A. Schekochihin,4 and T. Tatsuno6
1

2

Department of Astronomy, University of California, Berkeley, California 94720, USA
Department of Physics, CSCAMM, and IREAP, University of Maryland, College Park, Maryland 20742-3511, USA
3
Department of Physics and Astronomy, UCLA, Los Angeles, California 90095-1547, USA
4
Plasma Physics, Blackett Laboratory, Imperial College, London SW7 2AZ, United Kingdom
5
Princeton Plasma Physics Laboratory, Princeton, New Jersey 08543, USA
6
Department of Physics, University of Maryland, College Park, Maryland 20742-3511, USA.
(Received 21 November 2007; published 14 February 2008)
This Letter presents the first ab initio, fully electromagnetic, kinetic simulations of magnetizedturbulence in a homogeneous, weakly collisional plasma at the scale of the ion Larmor radius (iongyroscale). Magnetic- and electric-field energy spectra show a break at the ion gyroscale; the spectral
slopes are consistent with scaling predictions for critically balanced turbulence of Alfvén waves above theion gyroscale (spectral index 5=3) and of kinetic Alfvén waves below the ion gyroscale (spectral indicesof 7=3 for magnetic and 1=3 for electric fluctuations). This behavior is also qualitatively consistent
with in situ measurements of turbulence in the solar wind. Our findings support the hypothesis that thefrequencies of turbulent fluctuations in the solar wind remain well below the ion cyclotron frequency bothabove and below the ion gyroscale.
DOI: 10.1103/PhysRevLett.100.065004

PACS numbers: 52.35.Ra, 52.65.Tt, 96.50.Tf

Introduction.—A wide variety of astrophysical plasmas—the solar wind and corona, the interstellar and intracluster medium, accretion disks around compact objects—
are magnetized and turbulent. The turbulence is damped atsmall scales where the plasma is weakly collisional, requiring a kinetic description. It is often a good approximation to consider small fluctuations about an equilibrium
state with a uniform (or large-scale) dynamically strong
mean magnetic field (the Kraichnan hypothesis [1]). The
resulting (subsonic) magnetohydrodynamic (MHD) turbulence is believed to be a Kolmogorov-like cascade of
spatially anisotropic Alfvénic fluctuations [2]. Such anisotropy is observed in laboratory plasmas [3], the solar wind[4], and numerical simulations [5]. Assuming a criticalbalance between the linear frequencies and nonlinear decorrelation rates [2,6], the anisotropy is scale dependent
with wave numbers parallel and perpendicular to the local
mean field related by kk / k2=3
? . This implies that in mostastrophysical plasmas, the frequencies of the Alfvénic
fluctuations remain below the ion cyclotron frequency,!  kk vA  i , even as the perpendicular wavelength
reaches the ion gyroscale, k? i  1.
Such fluctuations are described by gyrokinetics (GK), a
rigorous low-frequency anisotropic limit of kinetic theory[7–10], which systematically averages out the particle
cyclotron motion. GK orders out the MHD fast wave andcyclotron resonances but retains finite Larmor radius effects and the collisionless Landau resonance. GK enables
numerical studies of kinetic turbulence with today’s computational resources because the gyroaveraging eliminates
fast time scales and reduces the dimensionality of phasespace from six to five. GK has been used to study electrostatic turbulence in fusion plasmas for decades, but there
0031-9007=08=100(6)=065004(4)have been few GK treatments of astrophysical turbulence.
GK is not applicable to large-scale, roughly isotropic fluctuations, such as are directly excited in the interstellar
medium by supernovae. However, the fluctuations in magnetized turbulence become lower amplitude and moreanisotropic at smaller scales. GK theory and simulations
are thus appropriate, and hold considerable promise, forstudies of microscopic phenomena such as turbulent heating and magnetic reconnection, and for interpreting observations of short-wavelength turbulent fluctuations. This
Letter reports the first ab initio, fully electromagnetic,
kinetic simulations of turbulence in a magnetized weakly
collisional astrophysical plasma.
The study of turbulence in weakly collisional plasmas
benefits from in situ measurements of the near-Earth solar
wind, illuminating the properties of turbulence from thelarge (energy-containing) scales to the small, kinetic scales
at which fluctuations are damped. The one-dimensionalfrequency spectrum of magnetic fluctuations typicallyshows a power-law behavior with a 5=3 slope at low
frequencies [11], a break at a few tenths of a Hz, and asteeper power law at higher frequencies with a slope thatvaries between 2 and 4 [12]. It is generally agreed thatthe 5=3 range is an MHD inertial range, while the break
and the dissipation-range slope have been attributed toproton cyclotron damping [13], Landau damping of kinetic
Alfvén waves (KAW) [14], or dispersion of whistler waves[15]. Recent simultaneous magnetic- and electric-fieldmeasurements found an increase in the wave phase velocity
above the spectral break [16], consistent with the conversion to a KAW cascade but inconsistent with cyclotron
damping [10]. The GK simulations presented below capture all of these spectral features with magnetic- and065004-1© 2008 The American Physical Society

PRL 100, 065004 (2008)

PHYSICAL REVIEW LETTERS

electric-energy spectra similar to those reported empirically in [16]. Our simulation results suggest that the turbulent spectra observed in the solar wind are a consequenceof the transition from an Alfvén-wave to a KAW cascade.
The code.—We have used the new code, ASTROGK,
developed to study astrophysical turbulence. ASTROGK isessentially a slab version of the publicly available code
GS2, used to study plasma turbulence in fusion devices[17]. We now give a brief overview of the code.
The simulation domain is a periodic flux tube with astraight uniform mean magnetic field B0 and no equilibrium gradients. All particle species have Maxwellian equilibrium distributions. The code solves the GK equation [8],evolving the perturbed gyroaveraged distribution function
hs x; y; z; "s ;  of the guiding centers for each species s—ions (protons) and electrons with the correct mass ratio
mi =me  1836. Spatial dimensions perpendicular to themean field x; y are treated pseudospectrally; a conservative finite-difference scheme is used in the parallel direction z. A gyroaveraged pitch-angle-scattering collision
operator [9] is used, with pitch-angle derivatives computedfrom second-order finite differences. The electromagnetic
field is represented by the scalar potential ’, parallel vector
potential Ak , and parallel magnetic field perturbation Bk .
These are determined from the quasineutrality condition
and Ampère’s law [8], where the charge densities andcurrents are calculated as velocity-space moments of the
distribution function. These velocity-space integrals (overparticle energies "s  ms v2 =2 and pitch angles   vk =v)
are done with spectral accuracy, using high-order
Gaussian-Legendre integration rules. The linear terms in
the GK system, including the field equations, are advanced
implicitly in time; for the nonlinear terms, an explicit,
third-order Adams-Bashforth scheme is used.
Linear benchmarks.—GS2 has been verified to describe
correctly the linear kinetic physics in parameter regimes
relevant to astrophysical plasmas [8]. ASTROGK has beenchecked to agree with GS2 exactly and benchmarked
against linear kinetic theory, as illustrated by Fig. 1: for

FIG. 1. Normalized frequencies !=kk vA and damping rates
=kk vA vs normalized perpendicular wave number k? i for aplasma with i  1 and Ti =Te  1. ASTROGK (squares) correctly reproduces the analytic results from the linear collisionlessgyrokinetic dispersion relation (line) [8].

week ending
15 FEBRUARY 2008k? i  1, we find Alfvén waves with !  kk vA andlittle damping; for k? i  1, these become kinetic Alfvénp
waves with !  kk vA k? i = i 2=1 Te =Ti  andstronger damping,in agreementwith linear theory [8,18].
p

Here vA  B0 = 4mi ni is the Alfvén speed, ni the ion
number density, Te and Ti the ion and electron temperatures, and i  8ni Ti =B20 .
Driving.—The driving and dissipation scales in astrophysical turbulence are widely separated —in the solar
wind, the effective driving scale is L  1011 cm and theion gyroscale is i  106 cm [10]—beyond the range
accessible numerically. Our simulation domain is understood to be much smaller than the driving scale. We model
the energy influx from larger scales by adding to Ampère’s
law a parallel ‘‘antenna’’ current jak;k . For each chosen
driving wave vector ka , the antenna amplitude is calculated from a Langevin equation whose solutions are Alfvénwaves with wave vector ka , frequency !  kak vA , and adecorrelation rate comparable to !. This driving is motivated by the expectation that turbulence in the inertial
range (at scales i  k1  L) is Alfvénic and criticallybalanced [2].
Dissipation.—The driving injects power into the system;in steady state, this power must be dissipated into heat. By
Boltzmann’s H theorem, no entropy increase, and therefore no heating, is possible in a kinetic system without
collisions. If the collision rate is low, the distributionfunction develops a small-scale structure in velocity space[8,9]. This makes the velocity derivatives in the collision
integral large so the collisions can act, a situation analogous to the emergence of small spatial scales in neutral
fluids with small viscosity (Kolmogorov cascade). In GK
turbulence, the cascades in position and velocity space arelinked, creating a kinetic cascade in five-dimensional
phase space [9]. Collisionless Landau damping of theelectromagnetic fluctuations leads to particle heating inthe sense that it transfers the electromagnetic fluctuationenergy into fluctuations of the particle distribution function(the kinetic entropy cascade [9]), which are then converted
into heat by collisions.
A detailed analysis of the kinetic cascade will be presented in a separate study, but the lesson is that kinetic
turbulence simulations must include collisions and havesufficient velocity-space resolution for the correct relationship to be established between small-scale structures invelocity and position space. Accomplishing this with aphysical collision operator simultaneously for ions andelectrons is difficult. To ease the resolution requirements,
we employ a hypercollisionality (analogous to hyperviscosity in fluid simulations) with the form of a pitch-anglescattering operator with a wave-number-dependent collision rate h k? =k?max 8 , where k?max is the grid-scalewave number. This artificially enhanced collision term
terminates the cascade and produces positive-definite heating close to the grid scale, but leaves essentially collision-

065004-2

PRL 100, 065004 (2008)

PHYSICAL REVIEW LETTERS

less physics at larger scales. For the ions, the importance ofthe hypercollisionality is marginal, while for the electrons
we need a large value of h . As a result, electron heating (atthe electron gyroscale e ) is not well modeled, but this isan acceptable sacrifice because our focus is on the turbulent cascade through the ion gyroscale at i  e .
Results.—The physical parameters in GK simulationsof plasma turbulence are i and Ti =Te . Here both are set
to 1, sensible characteristic values for the solar wind at 1
AU; a full parameter scan is desirable in the future. By
varying the driving wave number ka and the (hyper)collision rate, we may focus on various scale ranges. Here
we present results obtained for the inertial range (k? i 1) and the ion gyroscale (k? i  1). The normalized
magnetic-energy
R spectrum 2 is defined EB? k?  23
Lz =L? 2k? dzhjAk;k? zj i=8ni Ti , where k? is measured in units of 1
i , Lz and L? are parallel and perpendicular box dimensions, and the angle brackets denoteangle averaging over a wave number shell centered at
jk? j  k? with width 2=L? . The normalized electricenergy spectrum EE? k?  is defined similarly in terms of’k? , with an extra factor of c=vA 2 , where c is the speedof light.
In the inertial range, k? i  1, the reduced MHD
equations are the rigorous limit of GK for Alfvénic fluctuations [9], so kinetic turbulence in this regime must be
consistent with the numerical results obtained in MHD
simulations [5]. Figure 2 shows the normalized magnetic
and electric-energy spectra calculated gyrokinetically inthis regime. As expected for critically balanced Alfvénic
turbulence [2], these spectra are coincident, showing a. This is the first demonstrascaling consistent with k5=3?
tion of an MHD turbulence spectrum in a kinetic simulation. While not surprising, this result can be viewed as anonlinear benchmark.
Our main focus is on scales near k? i  1, a regime thatcannot be simulated by any fluid model. We know, however, that low-frequency Alfvénic turbulence is rigorously
described by reduced MHD equations for k? i  1 and

FIG. 2 (color online). Magnetic (solid line) and electric(dashed line) energy spectra in the MHD regime (k? i < 1).
The box size is L? =2  10i . Electron hypercollisionality isdominant for k? i 1 (dotted line).

week ending
15 FEBRUARY 2008

by a similarly reduced version of the electron MHD equations for k? i  1 [9]. The latter system supports kinetic
Alfvén waves (see Fig. 1). If one assumes a turbulent
cascade of KAW-like fluctuations decorrelating on a timescale comparable to the linear KAW period (critical balance), scaling arguments predict that the magnetic-energy
to k7=3, while the electricspectrum steepens from k5=3?
?
energy spectrum flattens to k1=3[9,10,19].
A spectral?
break at the transition between Alfvén-wave and KAW
turbulence is expected at k? i  1. Figure 3 shows theenergy spectra in our simulation around this transition.
Both the spectral break (at k? i ’ 2) and the steepening(flattening) of the magnetic-(electric-)energy spectra areobserved. The spectra at wave numbers below and above
the transition are consistent with the predictions for critically balanced Alfvén-wave and KAW cascades[2,9,10,19].
There is a striking similarity between the simulated
spectra shown in Fig. 3 and the magnetic- and electricenergy spectra in the solar wind reported in [16]. The
increase in phase velocity in the dissipation range (k? i >1), shown by both measurement and simulation, is compelling evidence that the observed breaks in the spectra arecaused by a transition to a KAW cascade, not by the onsetof ion cyclotron damping [10].
The scaling predictions for KAW turbulence are madeassuming negligible Landau damping. In our simulations,the damping is weak, so it is reasonable that the scaling
predictions are well satisfied. However, this will not be truein all real astrophysical situations. We have argued [10]often observed inthat the spectra much steeper than k7=3?

FIG. 3 (color online). Bold lines: normalized energy spectrafor B? (solid line), Bk (dash-dotted line), and E? (dashedline). Thin lines: solution of the turbulent cascade model of [10].
Dimensions are Nx ; Ny ; Nz ; N" ; N ; Ns   64; 64; 128; 8; 64; 2,requiring ’ 0:5 109 computational mesh points, with box size
L? =2  2:5i . Electron hypercollisionality is dominant fork? i 8 (dotted line).

065004-3

PRL 100, 065004 (2008)

PHYSICAL REVIEW LETTERS

the solar wind [12] may be due to significant Landau
damping. The effect of this damping on the energy spectra
can be estimated, as proposed in [10] (see also [18]), using
a spectral model of the turbulent cascade based on threeassumptions: (i) spectrally local energy transfer, (ii) criticalbalance, (iii) the applicability of the linear damping rates.
Using this model, the energy spectrum EB? k?  can bepredicted over the entire simulation range, given one‘‘Kolmogorov’’ constant, which quantifies the linear
damping rate relative to the nonlinear cascade rate. In
Fig. 3, we show that this analytical model reproduces theentire shape of the numerical spectrum. Using linear GK
eigenfunctions for KAWs enables the determination ofenergy spectra for the electric-field fluctuations (E? ) andfor the fluctuations of the magnetic field strength (Bk ).
The model works well without fine-tuning, for a range ofvalues of the constant; this is because the damping is small
in this simulation and our model captures the transition
from Alfvénic to KAW turbulence. The agreement between the analytical model and the simulations is a nontrivial result: it suggests that the linear damping rate does
not significantly underestimate the rate at which the electromagnetic energy is dissipated in the nonlinear simulations. Future simulations will determine whether stronger
linear damping can account for the steeper spectra often
observed in the solar wind.
Conclusions.—We have presented first-of-a-kind kinetic
simulations of turbulence in a weakly collisional, magnetized plasma. The ion-gyroscale turbulent fluctuations simulated here represent the fate of a larger-scale MHD
cascade. The qualitative agreement between our simulations and solar-wind measurements [16] supports theoretical models in which the turbulent fluctuations in the solar
wind have frequencies well below the ion cyclotron frequency even when the cascade reaches the (perpendicular)
scale of the ion Larmor radius. The observed break in themagnetic-energy spectrum in the solar wind is inferred tocorrespond to a transition to kinetic-Alfvén-wave turbulence, not to the onset of ion cyclotron damping. Although
half a billion mesh points were used in the case of Fig. 3,
the resolution in velocity space is still not fully sufficient todraw detailed conclusions about the turbulent heating.
Nonetheless, the agreement between the simulations andan analytical cascade model based on linear damping rates
implies that the latter do not significantly underestimate thetrue damping in a turbulent collisionless plasma. Future
simulations will probe a range of plasma parameters, in-week ending
15 FEBRUARY 2008

cluding more heavily damped regimes, that will allow amore quantitative study of the role of collisionless damping
in turbulent plasmas. The first results reported in this Letter
demonstrate that such kinetic simulations of plasma turbulence may be undertaken with some confidence, using
existing computational resources.
This work was supported by the U.S. DOE Center for
Multiscale Plasma Dynamics (G. G. H., T. T.), David andLucile Packard Foundation (E. Q.), U.K. STFC (A. A. S.),
Leverhulme Trust Network for Magnetized Plasma
Turbulence, and Aspen Center for Physics.

*ghowes@astro.berkeley.edu[1] R. H. Kraichnan, Phys. Fluids 8, 1385 (1965).
[2] P. Goldreich and S. Sridhar, Astrophys. J. 438, 763 (1995).
[3] D. C. Robinson and M. G. Rusbridge, Phys. Fluids 14,
2499 (1971); S. J. Zweben et al., Phys. Rev. Lett. 42,1270 (1979).
[4] J. W. Bieber et al., J. Geophys. Res. 101, 2511 (1996);
A. Bigazzi et al., Astrophys. J. 638, 499 (2006);
L. Sorriso-Valvo et al., Europhys. Lett. 75, 832 (2006);
K. T. Osman and T. S. Horbury, Astrophys. J. 654, L103
(2007).
[5] J. V. Shebalin et al., J. Plasma Phys. 29, 525 (1983); J. Choand E. T. Vishniac, Astrophys. J. 539, 273 (2000);
J. Maron and P. Goldreich, Astrophys. J. 554, 1175 (2001).
[6] J. C. Higdon, Astrophys. J. 285, 109 (1984).
[7] E. A. Frieman and L. Chen, Phys. Fluids 25, 502 (1982).
[8] G. G. Howes et al., Astrophys. J. 651, 590 (2006).
[9] A. A. Schekochihin et al., arXiv:0704.0044 [Astrophys. J.
Supp. (to be published)].
[10] G. G. Howes et al., arXiv:0707.3147 [J. Geophys. Res. (tobe published)].
[11] M. L. Goldstein et al., Annu. Rev. Astron. Astrophys. 33,
283 (1995).
[12] R. J. Leamon et al., J. Geophys. Res. 103, 4775 (1998);
C. W. Smith et al., Astrophys. J. Lett. 645, L85 (2006).
[13] P. J. Coleman, Astrophys. J. 153, 371 (1968).
[14] R. J. Leamon et al., J. Geophys. Res. 104, 22331 (1999).
[15] O. Stawicki et al., J. Geophys. Res. 106, 8273 (2001).
[16] S. D. Bale et al., Phys. Rev. Lett. 94, 215002 (2005).
[17] M. Kotschenreuther et al., Comput. Phys. Commun. 88,
128 (1995); W. Dorland et al., Phys. Rev. Lett. 85, 5579
(2000); F. Jenko et al., Phys. Plasmas 7, 1904 (2000).
[18] E. Quataert and A. Gruzinov, Astrophys. J. 520, 248(1999).
[19] D. Biskmap et al., Phys. Plasmas 6, 751 (1999); J. Cho and
A. Lazarian, Astrophys. J. 615, L41 (2004).

065004-4