J. Plasma Fusion Res. SERIES, Vol. 8 (2009)

Vlasov simulation of nite amplitude
magnetohydrodynamic waves in the solar wind :Development of Vlasov-Hall-MHD code
Takashi KUMASHIRO, Tohru HADA, Yasuhiro NARIYUKI

1)

2)
and Takayuki UMEDA

E.S.S.T., Kyushu University, 6-1, Kasuga-koen, Kasuga City, Fukukoka 816-8580, Japan
1)
Kochi National College of Technology, Kochi 783-8508, Japan

2)

Solar-Terrestrial Environment Laboratory, Nagoya University, Nagoya, Aichi 464-8601, Japan(Received: 6 September 2008 / Accepted: 16 January 2009)

Vlasov simulation is a method to solve time evolution of a plasma by directly time advancing
the distribution function in the position-velocity phase space. Vlasov simulation is free from thermal(numerical) noise and thus is advantageous in analyzing ne details of nonlinear plasma phenomena.
With this background in mind, we have developed a new Vlasov simulation code (1-d in the real space,
3-d in the velocity space), in order to study basic properties of nonlinear evolution of large amplitude
magnetohydrodynamic (MHD) waves in the solar wind. In contrast to traditional Vlasov simulations
in which electron waves are of major concern, our simulation code focuses on solving plasma behavioraround the ion scales, assuming a massless electron uid. Since we mainly deal with low frequency
MHD waves propagating quasi-parallel to the background magnetic eld, cyclotron coupling can beassumed to be weak for parameters typical to the solar wind. Thus, in our code, the Vlasov equation
is solved only along the longitudinal direction, whereas the Hall-MHD equations are solved for thetransverse directions. Propagation properties of Alfven and ion acoustic waves in the simulation agreewith analytically obtained dispersion relations. Some results on the parametric decay instability of
Alfven waves are also reported.
Keywords: Vlasov simulation, Hall-MHD, Alfven waves, Parametric instability, solar windapplications in the near future.
1. Introduction
With this background in mind, we have develVlasov simulation is a numerically stable method
oped a new Vlasov simulation code (1-d in con guto solve time evolution of a plasma by directly time
ration space, 3-d in the velocity space), in order toadvancing the distribution function in the positionstudy basic properties of nonlinear evolution of magvelocity phase space [1]. While particle simulations,netohydrodynamic (MHD) waves, in particular, quasisuch as the PIC (particle-in-cell) simulations, areparallel (to the background magnetic eld) propagatwidely used for studies of kinetic properties of aing Alfven waves in the solar wind. Quasi-parallel
plasma, these methods cannot avoid rather large stapropagating Alfven waves are thought to be ubiquitistical uctuations (numerical noise) arising from thetous in many astrophysical situations, for instance in
use of nite number of particles. It is particularly imthe solar wind, solar corona, interstellar matter, magportant to retain the noise level low for treating longnetosphere, among others, and are thought to playtime evolution of nonlinear systems (such as the paraimportant roles in the physical processes in these remetric instabilities of Alfven waves) since the presence
gions. [2]. Because they are typically robust for linearof (large) noise can qualitatively alter the eventual
ion-cyclotron damping due to their small wave frestate of system evolution. To keep numerical noise
quencies and for linear Landau damping due to theirlow, a large number of particles are needed. Unlike
small propagation angle relative to the backgroundconventional PIC simulations, the Vlasov simulation
magnetic eld, Alfven waves are thought to propagate
is free from thermal (numerical) noise, and is thus
very long distance [3]. Thus, the dissipation process ofadvantageous in analyzing ne details of nonlinear
these Alfven waves is important from a viewpoint ofplasma phenomena. On the other hand, Vlasov simthe momentum and energy transport in the solar wind.
ulations require huge computational resource, espeAmong various possible scenarios for the Alfven waves
cially for multi-dimensional models. Recent advances
to dissipate, the parametric instabilities are regarded
both in computing power and computational schemes,
as one of the most robust and important processes rehowever, make it possible to execute relatively large
sponsible for the dissipation of the solar wind Alfvenscale calculations. It is expected that the Vlasov simwaves.
ulations will be used for increasingly wider range ofWhile propagation of MHD waves can be underauthor's e-mail: kuma-taka@esst.kyushu-u.ac.jp831©2009 by The Japan Society of Plasma
Science and Nuclear Fusion Research

T. Kumashiro et al., Vlasov Simulation of Finite Amplitude Magnetohydrodynamic Waves in the Solar Wind: Development of Vlasov-Hall-MHD Codestood in a uid model, inclusion of ion kinetic e ects
is essential in correctly describing the properties ofthe solar wind Alfven waves since the ion pressure isabout the same order as the magnetic pressure [4]. To
be more precise, the ion kinetic e ects have to be taken
into account unless i is very small (for Alfven waves)or e = i is very large (for ion acoustic waves). In theabove, i and e are the ratio of the plasma pressure and the magnetic pressure for ions and electrons.
Also, electron kinetic e ects can be important in caseof nearly perpendicular propagating Alfven waves (kinetic Alfven waves), although we do not discuss themin the present paper.
An ultimate goal of the Alfven wave simulation
may be to perform an extremely large scale simulation including the entire 3+3=6 dimensions in the
con guration and the velocity phase space. However,computations in such a system require huge memory
size and computation time, and are not realistic atpresent. Thus we focus on parallel propagating solar wind Alfven waves and reduce required memory
size and computational time by assuming the following three points. (1) As mentioned above, dynamics of the quasi-parallel propagating MHD waves ismainly controlled by ions, and so we treat electrons
as a massless charge neutralizing uid, in a way similar to the well-known hybrid code [5, 6], and a Vlasovcode recently proposed by Valentini[7]. In our
simulation, the Hall-MHD equation set is used to include the dispersiveness of the MHD waves for transverse directions (see point (3) below). (2) We assume
that wave propagating directions do not change during the time scale of our interest, and limit ourselvesto the 1-dimensional problem (1-d in the real, 3-d inthe velocity space). Underlying here is the fact thatparallel Alfven waves stay dominant in the course ofparametric instabilities of Alfven waves we treat in
the present paper [8]. Many kinds of important andoutstanding problems could be considered even in the
restricted phase space assumed in the 1-d system. (3)Since we mainly deal with low frequency Alfven waves
propagating parallel to the background magnetic eld,
cyclotron coupling can be assumed to be weak for typical parameters of the solar wind [3]. In our simulation
code (Vlasov-Hall-MHD code), the Vlasov equation issolved only along the longitudinal direction whereas
the Hall-MHD equations are solved for the transverse
directions.
By employing these assumptions, ion distributionfunction can be represented only in the 1+1=2 dimensional phase space. Transverse ion velocities andelectron dynamics are described in terms of uid (averaged) variables. In order that the above assumptions
remain valid throughout the simulation runs, conditions below have to be satis ed. First, parameters
for plasma and waves are within the range in whichenergy exchange by cyclotron resonance can be negligible. Further, magnetic eld direction at any point
has to remain roughly parallel to the main axis (thex direction). For this sake, amplitude of the Alfvenwaves should not be too large.
In this sense, the concept of the present VlasovHall-MHD model is similar to Landau uid model, in
which leading order nite Larmor radius (FLR) corrections to the pressure are retained while contributions of non-gyrotropic heat ux components are neglected. It is reported that numerical integration of
the Landau uid model can accurately describe Alfvenwave instabilities [9, 10]. In a similar way, we expect that the present Vlasov-Hall-MHD model cancorrectly describe longitudinal characteristics of the
Alfven wave evolution.
The paper is organized as follows: in the next section, we describe in some detail the Vlasov-Hall-MHD
model. In Section 3, we present our simulation results on dispersion relation of the ion acoustic and the
Alfven waves. In Section 4, we discuss Alfven wave
parametric instabilities, and mention brie y an interesting ion acceleration took place during the courseof the Alfven wave decay instability. Section 5 is asummary and conclusions.
2. Vlasov-Hall-MHD code

In this section we brie y describe our VlasovHall-MHD model. Let us consider a phase space
which is 1-d in the con guration space (x direction) and 3-d in the velocity space. The backgroundmagnetic eld (=constant) is along the x direction,
and we assume that local orientations of the magnetic eld remain 'nearly' parallel to the x-axis everywhere as well,transverse magnetic eld due
to waves/ uctuations is small. Then, the ion distribution function F (x; t; v) may be separated into
the longitudinal and the transverse directions,
F (x; t; v) = f (t; x; vx )g (t; x; vy ; vz ). Note that thelongitudinal and the x directions are interchangeably
used here. Longitudinal and transverse distributionfunctions are normalized as follows,Zf (t; x; vx )dvx = (t; x)
(1)
and Z
g (t; x; vy ; vz )dvy dvz = 1(2)where (t; x) is the plasma density, and the integration is for the entire velocity space (for each variable).
Hence,Zf (t; x; vx ) = F (t; x; v)dvy dvz(3)
From now on, we normalize all the physical variables
using

et. al.

i.e.,

i.e.,

832

T. Kumashiro et al., Vlasov Simulation of Nite Amplitude Magnetohydrodynamic Waves in the Solar Wind: Development of Vlasov-Hall-MHD code

Alfven velocity, background magnetic eldstrength, and the ion cyclotron angular frequency(with ion sense of rotation being the negative frequency), all de ned at a certain reference point.
Integration of the Vlasov equation over vy and vz
yields,@f+ vx @f+ (ex + uy bz@t@x

uz b y )

@f
=0
@vx

3. Dispersion relation of ion acoustic
waves and Alfven waves

To ascertain whether the present Vlasov-HallMHD code runs correctly, we rst made thermal runs
and evaluated the dispersion relations, and comparethem with analytical dispersion relations. In caseof parallel propagation, MHD waves decouples into
longitudinal (ion acoustic) and transverse (Alfven)waves. Strictly speaking, non-propagating entropy
waves should also be included, but we do not consider
them here.
Assuming Maxwellian for ions,
1 exp  vx f (vx ) = p
(10)
ai
ai
where ai = 2Ti is ion thermal velocity squared, dispersion relation of the ion acoustic waves is given [14].
2T(11)Z (i ) = i ;Tewhere i = !=kai, ! is the angular wave frequency, kis the wave number, and Z () is the plasma dispersion
function. Assuming !=k < ai, expansion of Eq.(11)yields!r = vs k = (Te + 3Ti ) = k
(12)
3




   21 Te 2
3 (13)
Te
!i
=
exp
j!r j
8 Ti
2Ti 2
where ! = !r + i!i. Unless ion temperature is muchless than that of electrons, ion acoustic waves su erfrom Landau damping Eq.(13).
Next we consider Alfven waves in the Vlasov-HallMHD system. If we consider parallel propagating, circularly polarized waves, longitudinal parameters including the ion distribution function remain constant,and so we only need to consider the transverse uid
variables. Then the linear dispersion relation of the
Alfven waves is derived from Eqs.(6), (8).
! = k (1 + ! )(14)
Since the parallel circular polarized Alfven waves completely decouple from longitudinal ion dynamics, theunmodulated, parallel propagating Alfven waves do
not su er from the Landau damping. Since in our normalization the ion cyclotron angular frequency is 1,
the positive and negative ! corresponds to whistlermode waves (right-hand polarized) and ion-cyclotron
wave (left-hand polarized), respectively.
For the thermal run, the initial ion distributionfunction f was given as a superposition of Eq.(10)and small perturbations. The results of thermal run
are shown in Fig.1. Used parameters are Te = 0.16,
vs = 0.4, e = 0.32, ai = 0.1. Figures show (Top)time evolution of the transverse magnetic eld and(Bottom) time evolution of the density. Small amplitude Alfven waves and ion acoustic waves are seen to(4)

where e = (ex ; ey ; ez ) and b = (bx ; by ; bz ) are the normalized electric and magnetic eld, respectively. The
ion bulk velocity, u = (ux ; uy ; uz ), is given by
1 Z vF (t; x; v)dvx dvy dvz
(5)
u=

2

0



2

The perpendicular velocity evolves according to,@u= ux @u + 1 @b
(6)
?

@t

?

?

@x

0

 @x

where for brevity we have introduced complex perpendicular magnetic eld, b = by + ibz , and complex
perpendicular velocity, u = uy + iuz . Then, Maxwell
equations and the electron equation of motion lead to
1 @ ( jb j + pe ) (7)ex = uy bz + uz by @x 2
?

?

?

@b?
=
@t

@
(u b
@x x ?

u? +

1 2

2

i @b?
)
 @x

(8)

In the above, electrons are treated as a massless,isotropic uid, quasi-charge neutrality is assumed, andthe displacement current in neglected. Further assuming that the electron uid is isothermal, pe = Te ,from Eqs.(4) and (7) we have@f@f
= vx @f+ 1 @ ( jb j + Te ) @v(9)@t@x  @x 2xThen, Eqs.(1), (5), (6), (8), (9) constitute a closed set(bx is constant), which is the basic equation set of ourmodel.
Next we brie y describe computational techniqueused in the Vlasov-Hall-MHD code. The method we
employ for time advancement of the Vlasov equationis the PIC scheme (Positive Interpolation for hyperbolic Conservation laws) [12, 13], which is a recently
developed high order interpolation scheme proposed
for Vlasov simulations. The PIC scheme is shown tobe particularly accurate and ecient for highly nonlinear problems [13], and is suited for the present study
since long time evolution of Alfven wave parametricinstabilities and associated ion dynamics are strongly
nonlinear and delicate problems (small numerical errors can lead to qualitatively di erent results). For thepresent study, we used periodic boundary conditions
for x and free boundary conditions for vx .
?

2

2

2

2

0

833

T. Kumashiro et al., Vlasov Simulation of Nite Amplitude Magnetohydrodynamic Waves in the Solar Wind: Development of Vlasov-Hall-MHD code

propagate with constant velocities and in both positive and negative directions. By Fourier transforming
these variables above, we obtain the dispersion relation as shown in Fig.2. The results are in good agreement with the analytical expressions, Eqs.(12), (14).
4. Parametricwaves

instabilityofAlfven

Alfven waves play essential roles in many spaceplasma phenomena. In particular, they are believed tocontribute in heating and acceleration of plasmas via
parametric instabilities[15, 16, 17]. The parametricinstabilities are categorized into decay type (k < k)and self-modulational type (k > k) instabilities according to the wave number of excited ion acoustic
daughter waves, where k is the wave number of theparent Alfven wave and k is that of the ion acoustic
wave. In this paper we chose the set of simulation
parameters in such a way that the decay type is thedominant instability. As initial conditions, we givemonochromatic, right-hand polarized, parallel propagating Alfven waves in the system,b = b exp[i (! t k x)](15)
u = u exp[i (! t k x)](16)
where ! = k (1+! ), u = b =vp , and vp = ! =k
(Walen relation). Also we let initially  =  = 1,ux = 0, and f = f as given by Eq.(10). The parent
wave propagates toward the positive x direction.
We rst check that the present Vlasov-Hall-MHD
code correctly solves nonlinear evolution of Alfvenwaves by running it as well as a Hall-MHD code
( uid dynamics both in the longitudinal and transverse directions) and a hybrid code ( uid electronsand particle-in-cell ions) using the same set of physical parameters (typical for the decay instability), andby making comparison the decay instability growth
rates. Physical parameters used are: electron betae = 0.72, ion beta i = 0.01 (almost cold), (ion thermal velocity ai = 0.1, electron temperature
Te = 0.36, ion temperature Ti = 0.005), normalized
parent Alfven wave amplitude b = 0.5, wave number of the parent wave k = 0:417, (which leads to)the parent wave frequency ! = 0:513. Other important parameters are: number of grids is 1024 in
the real space and 1000 in the velocity space, system size = 256 (number of parent waves in the system m = L=(2=k ) = 17), velocity space grid sizev = 0:008, time step t = 0:01, and the maximum
time duration computed is tmax = 250.
From the Vlasov-Hall-MHD simulation run we
rst con rmed that the ion acoustic waves (k) and the
backward propagating Alfven waves (Stokes waves,k k < 0) were generated at respective wave number regimes expected from the theory. Higher harmonic waves of the Alfven waves were also observed.
000

Fig. 1 Results of the thermal run. (Top) Time evolutionof one transverse component of the magnetic eld,by , which represents presence and propagation of
Alfven waves. (Bottom) Time evolution of the density, , which basically represents the presence andpropagation of the ion acoustic waves.

?

0

0

0

?

0

0

0

2
0

2
0

0

0

0

0

0

0

0

0

0

i.e.,

0

0

0

0

Fig. 2 Dispersion relations obtained from the thermal run
shown in Fig.1. (Top) Dispersion relation obtained
for the Alfven waves. (Bottom) The same for theion acoustic waves. In both gures, analytical dispersion relations Eq.(12) and Eq.(14)) are superposed as squares.

0834

T. Kumashiro et al., Vlasov Simulation of Nite Amplitude Magnetohydrodynamic Waves in the Solar Wind: Development of Vlasov-Hall-MHD code

to the ion acceleration reported in hybrid simulations[11], will be discussed in a separate paper.

In addition to the Vlasov-Hall-MHD simulation, we
have also performed Hall-MHD and hybrid simulations. Table 1 compares the maximum growth rates
evaluated from the three di erent types of simulations
as well as the ( uid) theory, using the same set of physical parameters. The growth rates for all the simulation runs are approximately equal. We conclude fromthe present test simulation that our new Vlasov-HallMHD code can provide quantitatively reliable resultsfor our study on Alfven wave evolution. The computational load of the Vlasov-Hall-MHD code is almost
the same as the Hall-MHD code and much less thanthe hybrid code or full-Vlasov code. Despite this, time
evolution of the ion distribution function is clearly obtained, which is one of the greatest advantages of the
Vlasov-Hall-MHD code.
Type of code
Maximum growth rate ( )Hall-MHD
0.0690hybrid
0.0646
Vlasov-Hall-MHD
0.0632Analytical solution
0.0682

Fig. 3 (a) The ion distribution function (the upper part)and the longitudinal electric eld (the lower part)plotted at an early stage of the simulation, t = 270.
Ion acoustic waves generated by the decay instability strongly modulate the ion distribution function.
(b) The ion distribution function (the upper part)and the magnetic eld magnitude (the lower part)plotted at a later stage, t = 770. Some ions arestrongly accelerated (x  150), in association ofenhanced magnetic eld perturbation.

Table 1Comparison of the decay instability maximum
growth rates evaluated by various simulation runs
and the analytical theory of Vlasov-Hall-MHD system.

Next we brie y discuss evolution of the ion distribution function due to the decay instability of righthand polarized Alfven waves (Fig.3). Parameters used
are the same as before, except that e = 0.32, b =0.4, k = 0:196 (! = 0:216), and m = L=(2=k ) = 8.
The upper part of Fig.3 (a) shows a snapshot of theion distribution at an early stage of the decay instability (t = 270). Modulation of the ion distribution isclearly seen due to the longitudinal electric eld (thebottom part of the gure) of the ion acoustic waves.
In particular, slashing, steepening, and trapping ofions are evident in association with the large amplitude longitudinal electric eld. Upper part of Fig.3(b) shows a snapshot of ion distribution at the later
phase of the evolution (t = 770), where the rst decay instability (.e., decay of the initial parent Alfvenwave: successive decay of daughter waves may be inprogress but not seen evidently) is saturated. Aroundx  150, some ions are strongly accelerated toward
negative vx direction. This acceleration is associatedwith an enhanced magnetic eld modulation generated by the beating between parent and the backward
propagating daughter Alfven waves (the lower part ofFig.3 shows the magnetic eld). The same type of
ion acceleration can be found also for the decay ofthe left-hand polarized Alfven waves. Details of theion acceleration mechanism, as well as its relationship
0

0

0

0

5. Summary

In this paper we have introduced our newly developed Vlasov simulation code (Vlasov-Hall-MHDcode), in which the Vlasov equation is solved only
for the longitudinal direction and the Hall-MHD setof equations are solved for the transverse directions.
The electrons are treated as a massless uid. The codeis particularly suited for solving nonlinear long time
evolution of quasi-parallel MHD waves, where transverse kinetic e ects (such as cyclotron damping) canbe neglected. Below we summarize results of some of
our initial simulation runs.
(1) By performing thermal runs, we con rmedthat the ion acoustic and the Alfven waves, which satisfy the linear dispersion relation, are present in thesimulation system,(2) By giving circularly polarized, monochromatic
Alfven waves plus a very small random uctuations as
initial conditions, we simulated the evolution of Alfvenwave decay instability. All the daughter waves are
generated, as well as the higher harmonics.
(3) Comparison of growth rates was made, by835

T. Kumashiro et al., Vlasov Simulation of Nite Amplitude Magnetohydrodynamic Waves in the Solar Wind: Development of Vlasov-Hall-MHD coderunning the present Vlasov-Hall-MHD code, the HallMHD code, and the hybrid code, using the same setof physical parameters. The present code is con rmed
to yield reliable result.
(4) It is demonstrated that ne details of the ion
distribution function can be captured by the present
simulation code. In particular, an interesting ion acceleration was found to take place in association withthe decay instability, for both right- and left-hand polarized parent waves.
Details of the ion acceleration, modulational instabilities, and conprehensive arguments on the Alfvenwave parametric instabilities (with a particular emphasis on longitudinal ion kinetic e ects) will be discussed in a forthcoming paper.

ble Alfven-cyclotron waves, Phys. Rev. Lett. (2008),
100, 125003.
[12] T. Umeda, M. Ahsour-Abdalla and D. Schriver, Comparison of numerical interpolation schemes for onedimensional electrostatic Vlasov code, J. Plasma Phys.
(2006), 72, 1057-1060.
[13] T. Umeda A conservative and non-oscillatory scheme
for Vlasov code simulations Earth Planets Space,(2008), 60(7), 773-779.
[14] F. F. Chen Introduction to plasma physics andcontrolled fusion, vol.1, "Plasma Physics", (1984),
Plenum Press, New York[15] A. A. Galeev, V. N. Oraevskii, The instability of
Alfven waves, Sov. Phys. Dokl. (1963), 7, 988-1003.
[16] M. L. Goldstein, An instability of nite amplitude circularly polarized Alfven waves, Astrophys. J. (1978),219 (2), 700.
[17] Y. Nariyuki and T. Hada, Parametric instabilities ofparallel propagating incoherent Alfven waves in a nite
ion beta plasma, Phys. Plasma, (2007) 14, 122110.

Acknowledgment
This work was in part supported by Grant-inAid for Young Scientists (Start-up) No.19840024 fromJSPS.
[1] C. Z. Cheng and G. Knorr, The integration of the
Vlasov equation in con guration space, J. Comp. Phys.
(1976), 22, 330.
[2] R. Bruno and V. Carbone, The solar wind as a turbulence laboratory, Living Rev. Solar Phys. (2005), 2.
4.
[3] S. P. Gary and J. E. Borovsky, Alfven-cyclotron uctuations: Linear Vlasov Theory, J. Geophys. Res. (2004),
109, doi:10.1029.
[4] Y. Nariyuki and T. Hada, Consequences of nite ion
temperature e ects on parametric instabilities of circularly polarized Alfven waves, J. Geophys. Res. (2007),VOL. 112, A10107[5] D. Winske, Hybrid simulation codes with application
to shocks and upstream waves, Space Sci. Rev. (1985),
42 (1-2), doi. 10.107/BF00218223.
[6] T. Terasawa, M. Hoshino, J.-I. Sakai and T. Hada, Decay instability of nite-amplitude circulary polarized
Alfven waves: A numerical simulation of stimulated
Brillouin Scattering, J. Geophys. Res. (1986), 91, 4171
-4187.
[7] F. Valentini, P. Travnicek, F. Califano, P. Hellinger, A.
Mangeney, A hybrid-Vlasov model based on the current advance method for the simulation of collisionlessmagnetized plasma, J. Comp. Phys. (2007), 225, 753770.
[8] Y. Nariyuki, S. Matsukiyo and T. Hada, Parametric instabilities of large-amplitude parallel propagating Alfven waves : 2D PIC simulation New journal of
Phys. (2008), doi. 10.1088/1367-2630/10/8/083004[9] G. Bugnon, T.Passot, and P. L. Sulem, Landau- uid
simulations of Alfven-wave instability in a warm collisionless plasma, Nonl. Proc. Geophys. (2004), 11, 609618.
[10] T.Passot, and P. L. Sulem, A landau- uid modelfor dispersive magnetohydrodynamics, Phys. Plasma,(2004), 11, 5173-5189.
[11] J. A. Araneda, E. Marsch, A. F. Vinas, Proton coreheating and beam formation via parametrically unsta-

836