Vlasov simulations of collisionless selfgravitating systems and astrophysical plasmas
Kohji Yoshikawa
Division of Astrophysics
Center for Computational Sciences, University of Tsukuba
CCS-EPCC Workshop
December 7 - 8, 2017, CCS, Univ. of Tsukuba

collaborators:

Satoshi Tanaka (CCS, Univ. of Tsukuba)
Takashi Minoshima (JAMSTEC)
Naoki Yoshida (Univ. of Tokyo)

N-body Simulation ofSelf-Gravitating Systems
a standard method for simulations of selfgravitating systems (galaxies, clusters of galaxies,
the large-scale structure in the universe) for morethan 40 years.
the mass distribution is sampled by particles in
the 6D phase-space volume (x, p) in a MonteCarlo manner.

sophisticated algorithms to treat large numberof particles such as Tree and TreePM methods
developedIshiyama 2013

Particle-In-Cell (PIC) Simulation of Plasmas
Particle-based approach to solve collisionless(astrophysical) plasma

E- and B-fields are computed in the finite
difference manner

3D PIC simulation of collisionless shock

Matsumoto et al. (2017)
PIC simulation of magnetic reconnection

particle acceleration in collisionless shockmagnetic reconnection(pCANS web page)

Drawbacks of Particle Simulations
intrinsic contamination of shot noise in physical quantities
not good at simulating kinetic physical processes
in which the tail of the distribution function plays
important rolesv

matter in the tail is not fairly sampled in particle
simulations
collisionless damping, two-stream instability
magneto-rotational instability starting fromhigh-beta plasmas
The grid spacing of PIC simulations should be less than the Debye length.
difficult to simulate phenomena in the macroscopic scalexVlasov Simulations
Simulate self-gravitating systems and plasmas by directly solving
the Vlasov equation in a finite volume manner: distribution function in six-dimensional phase space
Vlasov-Poisson simulation

Vlasov-Maxwell simulation
magnetized plasma

self-gravitating systemelectro-static plasma

Vlasov Simulations
Free from shot noise contamination
Good at capturing kinematic physical processes
collisionless damping, dynamical friction etc in self-gravitating systems
Landau damping, weibel instability etc in astrophysical plasma
Free from mesh size constraint to the Debye length
requires huge amount of memory and computational costs
had not been practical for a long timesolves hyperbolic PDE instead of ODE in particle simulations
inevitable numerical diffusionneed for sophisticated less diffusive numerical schemes

Vlasov-Poisson Simuation in 6D phase space
Yoshikawa, Yoshida, Umemura (2013)
directional splitting method: advection along ℓ-direction
numerical scheme for a one-dimensional advection equation
Positive and Flux Conservative (PFC) scheme

Filbet, Sonnendrucker, Bertrand (2001)

Merging of Two Spheres
offset merging of two king self-gravitating
spherescomparison with a equivalent N-bodysimulations, in which each King sphereis represented with a million particles
Vlasov simulation

N-body simulation

Grav. Instablity and Collisionless Damping
Initial condition: gravitational instability: collisionless damping
The initial density fluctuation δ (x) is given so that it has a white noise power
spectrum.

lines : theoretical prediction

New High-Order Schemes
for Vlasov Simulation

Higher-Order Advection Schemes
curse of dimensionality in Vlasov simulations
huge memory consumption due to high dimensionality of phase space
size of numerical simulation is limited by the amount of available memoryhow to overcome
adaptive mesh refinement
higher-order schemes for the advection equation
mathematical and physical requirements
monotonicity
monotonicity-preserving (MP) constraints (Suresh & Huynh 1997)
positivity
a new positivity-preserving (PP) limiter in the form of a flux limiter
spatially fifth- and seventh-order schemes with monotonicity- andpositivity-preserving features.
(c.f. The PFC scheme has a spatially third-order accuracy.)

Naming Convention of the Schemes
XX-MPP5
RK : 3-stage temporary 3rd order Runge-Kutta time integration
SL : single-stage semi-Lagrangian time integration
MPP : monotonicity and positivity preserving
MP : monotonicity preserving
final digit : the spatial order of accuracyLinear Advection
linear advection of

The result of PFC scheme is severelysmeared.
Schemes w/o PP limiter yield significant
negative value around the minima.
Only 7th-order scheme (SL-MPP7) canreproduce the maxima.

1D Self-Gravitating System
one-dimensional space with periodic boundary condition
initial condition

critical Jeans wave number

gravitational instability
collisionless damping

1D Self-Gravitating System
Nx = Nv =64
As time proceeds, numerical diffusion
takes place and smear small structuresin the lower-order schemes.

with PP limiter

w/o PP limiter

Negative regions in the lower panels
disappear in the results with positivity
preserving limiter.

Application To Cosmological Neutrinos

Cosmological Relic Neutrinos
massive neutrinos in the universe
lots of neutrinos in our universe
decoupled at early stage of theuniverse when they are relativistic.

currently non-relativistic andgravitationally interacting withcold dark matter (CDM)

Discovery of neutrino oscillation

absolute mass of neutrinos
and its hierarchy are still
unknown

dynamical effect of massive neutrinos
free streaming (collisionless damping)large velocity dispersion of neutrinos
growth of density fluctuation suppressed beyond the damping scaleCollisionless Damping on LSS

density fluctuation suppressed at scales smaller than the damping scale.
the amount of suppression depends on the mass of neutrinos
the mass (and its hierarchy) of neutrinos can be estimated with suchdamping feature.
non-linear features should be investigated with numerical simulations.

Initial Conditioncosmological parameters
PLANCK 2015 results：
curvature fluctuation：(pivot scale : k = 0.002 Mpc-1)

total neutrino mass

initial condition created at redshift of zi = 10
computaional domain: Lbox = 20000h-1 Mpc, 2000h-1 Mpc, 200h-1 Mpc
number of particles / number of mesh grids
N-body simulation :Vlasov simulation：

CDM and neutrino distribution
CDMneutrino

200 h-1 Mpc

200 h-1 Mpc

Power Spectrum:

ratio of power spectra with massive and massless neutrinos

damping scale :density fluctuation with k > 3x10-2 h/Mpc damps owing to collisionless dampingconsistent with the perturbation theory (e.g. Saito et al. 2009) in the early stages
up turn at k>1 h-1Mpc, probably due to dynamical feedback by CDM

CDM and neutrino distribution
contour : CDM

color : neutrino

diffuse distribution of neutrinos
owing to its large velocity
dispersion
regions with similar CDM density
can have significantly different
neutrino density
remarkable offset in density peaksof CDM and neutrinos
neutrino wake / dynamical friction ?
Zhu et al. (2016)200 h-1 Mpc

40 h-1 Mpc

20 h-1 Mpc

Vlasov-Maxwell Simulation

Vlasov-Maxwell Simulation
solve the distribution functions
for both of ions and electrons
difficulty to solve gyro-motionin the velocity space corrrectly
in a long term
rigid-body rotation problem
rigid-body rotation of a 2D gaussian profile

temporal variation of dispersion
CSL2 : CIP-CSL2
Takizawa et al. (2002)
BKSB :
backsubstitution method
Schmitz & Grauer (2006)

Minoshima, Matsumoto, Amano (2011)

Rigid-Body Rotation with Our Scheme
initial condition

SL-MPP7

PFC

100 rotations

PFC (3rd order)

rigid-body rotation of a gaussian profile in a 2D
plane
SL-MPP5 and SL-MPP7 schemes yields only a few
per cent increase in velocity dispersion.

SL-MPP5 (5th order)

The lower-order scheme suffers from significant
numerical heating.

SL-MPP7 (7th order)Magnetic Reconnection with VlasovMaxwell Simulations
Vlasov-Maxwell simulation in 5D phase spaceNx= 432, Ny = 216

color : Jz

Nvx = Nvy = Nvz = 32
Hall effect to trigger the fastreconnection

dissipation of magnetic field

momentum transport to upstream region

Magnetic Reconnection with VlasovMaxwell Simulations
Vlasov-Maxwell simulation in 5D phase spaceNx= 432, Ny = 216

color : Jz

Nvx = Nvy = Nvz = 32
Hall effect to trigger the fastreconnection
Deviation from MHD approximation

dissipation of magnetic field

momentum transport to upstream region

Summary
Vlasov simulations in 6-dimensional phase space are now practical.
New high-order advection schemes with monotonicity- andpositivity-preservation and with single-stage time integration
Vlasov-Poisson simulation of cosmic neutrinos in the large-scale
structure in the universe
Our new schemes are also suitable for Vlasov-Maxwell simulation
owing to good accuracy in rigid-body rotation.
Vlasov-Maxwell simulations of magnetic reconnection in 5D phasespace show kinetic features of magnetic recconection in amacroscopic MHD scale.