Simulation of collisionless plasma with the Vlasov method
Takayuki Umeda
Solar-Terrestrial Environment Laboratory, Nagoya University, Japan1Introduction

Space plasma is a collisionless, multi-scale, and highly
nonlinear medium. There are numerous types of selfconsistent computer simulations that treat space plasma
according to various approximations. The global-scale
dynamics are commonly described by magneto-hydrodynamic (MHD), Hall-MHD and multi-ﬂuid models, whileelectron-scale processes are described by the kinetic model,i.e., the Maxwell equations and either the Newton-Lorentz
equation for charged particles or the Vlasov (collisionlessBoltzmann) equation. Hybrid methods treat ions as particles and electrons as a ﬂuid for ion-scale processes.
Conventionally, MHD simulations have been used fornumerical modeling of global-scale problems such as magnetospheres of stars and planets. However, the MHD simulations need diﬀusion coeﬃcients, which are essentially due
to kinetic processes that are eliminated in the frameworkof the MHD approximation. Recent high-resolution in-situ
observations have also suggested that ﬂuid scale and kinetic scale in space plasma are strongly coupled with eachother, which is called cross-scale coupling. To understandthe cross-scale coupling in space plasma, it is important
to include full kinetics in global-scale simulations, which isthe goal of this study.
2

Overview of Numerical Schemes

The Vlasov model solves the kinetics equations of spaceplasma, i.e., the Maxwell equations (1) and the Vlasov(collisionless Boltzmann) equation (2),⎫1 ∂E ⎬
∇ × B = μ0 J + 2c ∂t(1)∂B⎭∇×E = −∂t
∂fs
∂fs
∂fsqs
[E + v × B]+v+=0(2)∂t∂rms∂v
where E, B, J, ρ, μ0 , 0 and c represent electric ﬁeld,magnetic ﬁeld, current density, charge density, magnetic
permeability, dielectric constant and light speed, respectively. The Vlasov equation (2) describes the development of the distribution functions by the electromagnetic(Lorentz) force, with the collision term in the right handside set to be zero. The distribution function fs (r, v, t)
is deﬁned in position-velocity phase space with the subscript s being the species of singly-charged particles (e.g.,s = i, e for ions and electrons, respectively). The Maxwell
equations and the Newton-Lorentz equations or the Vlasov
equation are coupled with each other via the current density J that satisﬁes the continuity equation for charge∂ρ+∇·J =0(3)∂tThese equations are regarded as the “ﬁrst principle” of thecollisionless plasma.

The Vlasov equation (2) consists of two advection equations with a constant advection velocity and a rotation
equation by a centripetal force without diﬀusion terms.
To simplify the numerical time-integration of the Vlasov
equation, we adopt a modiﬁed version of the operator splitting [1],∂fs
∂fs+v=0(4)
∂t
∂r
∂fs
qs ∂fsE+=0∂t
ms ∂v

(5)

qs
∂fs
∂fs+=0[v × B]
∂t
ms∂v

(6)

Equations (4) and (5) are scalar (linear) advection equations in which v and E are independent of r and v, respectively. We adopt a multidimensional conservative semiLagrangian scheme [1] for solving the multidimensional
advection equations. In the full electromagnetic method,
it is essential to use conservative schemes for satisfying
the continuity equation for charge. With the multidimensional conservative semi-Lagrangian scheme, the continuity equation for charge (3) is exactly satiesﬁed. In thepresent study, we compute the numerical ﬂux by using
the multi-dimensional advection scheme [1] with a positive, non-oscillatory and conservative limiter [2, 3] for stable time-integration of advection equations. Equation (6),on the other hand, is a multi-dimensional rotation equation which follows a circular motion of a proﬁle at constant
speed by a centripetal force. For stable rotation of the proﬁle on the Cartesian grid system, the “back-substitution”
technique [4] is applied. In addition, Maxwell’s equations
are solved by the implicit Finite Diﬀerence Time Domain(FDTD) method.
The velocity distribution function has both
conﬁguration-space and velocity-space dimensions,and deﬁned as a hyper-dimensional (>3D) array, There
are some additional communications overhead in parallelizing over the velocity-space dimensions since areduction operation is required to compute the charge andcurrent densities (the zeroth and ﬁrst moments) at a given
point in conﬁguration space. Thus we adopt the “domaindecomposition” only in conﬁguration space, where thedistribution functions and electromagnetic ﬁelds are decomposed over the conﬁguration-space dimensions. This
involves the exchange of ghost values for the distributionfunction and electromagnetic ﬁeld data along boundaries
of each procession element. The non-oscillatory andconservative scheme [2, 3] uses six grids for numericalinterpolation, and three ghost grids are exchanged byusing the “Mpi_Sendrecv() ” subroutine in the standard
message passing interface (MPI) library for simplicity
and portability [5]. Note that the code allows threadparallelization over the velocity-space dimensions viaOpenMP.

Fig. 1. Performance of diﬀerent supercomputer systems with scalar processors. (left) Computational speed and (right)
Scalability as a function of the number of cores. The star, asterisk, square, circe, and downward-pointing-triangle marks
represent the results on the HA8000, FX1, HX600, RX200S6, and R815. respectively.

3

Performance Evaluation

We conduct the performance measurement test of
our parallel Vlasov code with a phase-space grid of(Nvx , Nvy , Nvz , Nx , Ny ) = (30, 30, 30, 40, 20) on one core,
which corresponds to a weak-scaling test with 1GB/core.
The tests were completed on various scalar-type parallel supercomputers, Hitachi HA8000 at the University ofTokyo, Fujitsu FX1 at Nagoya University and JAXA, Fujitsu HX600 at Nagoya University, Fujitsu RX200S6 at
Kyushu University, and DELL PowerEdge R815 at SolarTerrestrial Environment Laboratory, Nagoya University.
Figure 1 shows the computational performance on thesesystems. As seen, we obtained a high scalability of over
80% with 1,000 cores. However, the scalability becomesworse when the all computational resources (3,072 and12,032 cores on FX1, and 2,304 cores on R815) are used.
It should be noted that the internode-connection deviceof the HA8000 system is Myrinet-10G, and the scalability
becomes worse with more than 1,024 cores because of thenetwork bandwidth capacity.
4Fig. 2. Structures of bow shock and wake tail formed asa result of the interaction between a weakly-magnetizedobject and the solar wind. Spatial proﬁle of ion density
together with the conﬁguration of magnetic ﬁeld lines obtained by a global Vlasov simulation.

Conclusion

For studying multi-scale processes in space plasma, computer simulations with the ﬁrst-principle (kinetic) modelare essential. We develop numerical schemes for Vlasov
simulations for practical use on currently-existing supercomputer systems. The weak-scaling benchmark testshows that our parallel Vlasov code achievs a high scalability. Currently, we use 256-1024 cores for parallel
computations and apply the present parallel Vlasov code
to “global” simulation on the interaction between solar/stellar winds and a weakly-magnetized small body(Fig.2).
Acknowledgement
The author is grateful to Y. Ito, K. Fukazawa, andY. Nariyuki for discussions. This work was supported
by MEXT/JSPS under Grant-in-Aid for Young Scientists(B) No.23740367. The computational resources are provided as a STEL computational joint research program, aNagoya University HPC program, and a JHPCN program.

References[1] T. Umeda, K. Togano, and T. Ogino, Two-dimensional
full-electromagnetic Vlasov code with conservative scheme
and its application to magnetic reconnection, Comput.
Phys. Commun., 180, 365–374, 2009.
[2] T. Umeda, A conservative and non-oscillatory scheme for
Vlasov code simulations, Earth Planets Space, 60, 773–779, 2008.
[3] T. Umeda, Y. Nariyuki, and D. Kariya, A non-oscillatory
and conservative semi-Lagrangian scheme with fourthdegree polynomial interpolation for solving the Vlasov
equation, Computer Physics Communications, 183, 1094–1100, 2012.
[4] H. Schmitz and R. Grauer, Comparison of time splitting and backsubstitution methods for integrating Vlasov’s
equation with magnetic fields, Comput. Phys. Commun.,
175 86–92, 2006.
[5] T. Umeda, K. Fukazawa, Y. Nariyuki, and T. Ogino, A
scalable full electromagnetic Vlasov solver for cross-scale
coupling in space plasma, IEEE Trans. Plasma Sci., inpress.