Multidimensional Vlasov–Poisson Simulations with High-order Monotonicity- and Positivity-preserving Schemes

The Vlasov Equation (2) can be broken down into six one-dimensional advection equations along each dimension of the phase space: three for advection in the physical space

$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+{v}_{i}\displaystyle \frac{\partial f}{\partial {x}_{i}}=0\,\,(i=1,2,3),\end{eqnarray} \tag{ 4 }$

and another three in the velocity space

$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}-\displaystyle \frac{\partial \phi }{\partial {x}_{i}}\displaystyle \frac{\partial f}{\partial {v}_{i}}=0\,\,(i=1,2,3).\end{eqnarray} \tag{ 5 }$

We solve these six equations sequentially. The gravitational potential ϕ in Equation (5) is updated after advection Equation (4) in the physical space is solved and the mass density distribution in the physical space is updated. Clearly, an accurate numerical scheme for a one-dimensional advection scheme is a crucial ingredient of our Vlasov solver.

Let us consider a numerical scheme of the advection equation

$\begin{eqnarray}&&\displaystyle \frac{\partial f(x,t)}{\partial t}+c\displaystyle \frac{\partial f(x,t)}{\partial x}=0\end{eqnarray} \tag{ 6 }$

on a one-dimensional uniform mesh grid. In a finite-volume manner, we define an averaged value of $f(x,t)$ over the ith grid centered at $x={x}_{i}$ with the interval of ${\rm{\Delta }}x$ at a time $t={t}^{n}$ as

$\begin{eqnarray}&&{f}_{i}^{n}=\displaystyle \frac{1}{{\rm{\Delta }}x}{\int }_{{x}_{i-1/2}}^{{x}_{i+1/2}}f(x,{t}^{n})\,{dx},\end{eqnarray} \tag{ 7 }$

where ${x}_{i\pm 1/2}={x}_{i}\pm {\rm{\Delta }}x/2$ are the coordinates of the boundaries of the ith mesh grid. Without loss of generality, we can restrict ourselves to the case with a positive advection velocity $c\gt 0$.

Physical and mathematical considerations impose the following requirements on numerical solutions of the advection equation: (i) Monotonicity—since the solution of the advection Equation (6) preserves monotonicity, so should its numerical solution. (ii) Positivity—the distribution function must be positive by definition, and thus the numerical solution should be non-negative. Note that monotonicity is not sufficient to ensure positivity of the numerical solution.

2.1.1. MP Schemes with the Positivity-preserving Limiter

We first present the spatially fifth-order scheme for the advection Equation (6) based on the MP5 scheme of Suresh & Huynh (1997), in which the value of $f(x,t)$ at the interface between the $(i-1)$th and ith mesh grids is reconstructed with five stencils ${f}_{i-3}^{n}$, ${f}_{i-2}^{n}$, ${f}_{i-1}^{n}$, fⁿ_i and ${f}_{i+1}^{n}$ as

$\begin{eqnarray}&&{f}_{i-1/2}^{\mathrm{int},5}=(2{f}_{i-3}-13{f}_{i-2}+47{f}_{i-1}+27{f}_{i}-3{f}_{i+1})/60.\end{eqnarray} \tag{ 8 }$

In order to ensure monotonicity of the numerical solution, we impose the constraint

$\begin{eqnarray}&&{\bar{f}}_{i-1/2}=\mathrm{MP}({f}_{i-1/2}^{\mathrm{int},5},{f}_{i-3},{f}_{i-2},{f}_{i-1}.{f}_{i},{f}_{i+1}).\end{eqnarray} \tag{ 9 }$

The detailed prescriptions of the constraints are given in Appendix A. The MP5 scheme adopts the three-stage TVD Runge–Kutta (TVD–RK) time integration as

$\begin{eqnarray}{f}_{i}^{(1)} & = & {f}_{i}^{n}-\nu {L}_{i}({f}^{n})\\ {f}_{i}^{(2)} & = & \displaystyle \frac{3}{4}{f}_{i}^{n}+\displaystyle \frac{1}{4}({f}_{i}^{(1)}-\nu {L}_{i}({f}^{(1)}))\\ {f}_{i}^{n+1} & = & \displaystyle \frac{1}{3}{f}_{i}^{n}+\displaystyle \frac{2}{3}({f}_{i}^{(2)}-\nu {L}_{i}({f}^{(2)})),\end{eqnarray} \tag{ 10 }$

where $\nu =c{\rm{\Delta }}t/{\rm{\Delta }}x$ is the CFL number and the operator L_i(f) is defined as

$\begin{eqnarray}&&{L}_{i}(f)={\bar{f}}_{i+1/2}-{\bar{f}}_{i-1/2}.\end{eqnarray} \tag{ 11 }$

Hereafter, this scheme is referred to as RK–MP5.

As already noted above, the MP scheme does not ensure positivity of the numerical solution. We thus introduce a positivity-preserving (PP) limiter based on the method proposed by Hu et al. (2013). By denoting ${U}_{i}^{+}={f}_{i}^{n}\,-2\nu {\bar{f}}_{i+1/2}^{n}$ and ${U}_{i}^{-}={f}_{i}^{n}+2\nu {\bar{f}}_{i-1/2}^{n}$, the single-stage Eulerian time integration can be written as

$\begin{eqnarray}&&{f}_{i}^{(1)}=\,{f}_{i}^{n}-\nu ({\bar{f}}_{i+1/2}^{n}-{\bar{f}}_{i-1/2}^{n})\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{1}{2}({f}_{i}^{n}-2\nu {\bar{f}}_{i+1/2}^{n})+\displaystyle \frac{1}{2}({f}_{i}^{n}+2\nu {\bar{f}}_{i-1/2}^{n})\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{1}{2}{U}_{i}^{+}+\displaystyle \frac{1}{2}{U}_{i}^{-}.\end{eqnarray} \tag{ 14 }$

Here, we consider a linear combination of the interface value obtained with the MP scheme and that with the spatially first-order upwind scheme as

$\begin{eqnarray}&&{\hat{f}}_{i+1/2}^{n}={\theta }_{i+1/2}{\bar{f}}_{i+1/2}^{n}+(1-{\theta }_{i+1/2}){f}_{i+1/2}^{\mathrm{UP},n},\end{eqnarray} \tag{ 15 }$

where ${f}_{i+1/2}^{\mathrm{UP},n}$ is expressed as

$\begin{eqnarray}&&{f}_{i+1/2}^{\mathrm{UP},n}=\displaystyle \frac{1}{2c}[{f}_{i}^{n}(c+| c| )+{f}_{i+1}^{n}(c-| c| )]={f}_{i}^{n}\end{eqnarray} \tag{ 16 }$

and ${\theta }_{i+1/2}$ is a parameter with a condition of $0\leqslant {\theta }_{i+1/2}\leqslant 1$. Note that the spatially first-order upwind scheme preserves positivity of the numerical solution if the CFL condition $\nu \leqslant 1$ is satisfied. Then, the time integration of the combined interface values (15) can be written as

$\begin{eqnarray}&&{f}_{i}^{(1)}={f}_{i}^{n}-\nu ({\hat{f}}_{i+1/2}^{n}-{\hat{f}}_{i-1/2}^{n})\end{eqnarray} \tag{ 17 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{1}{2}({f}_{i}^{n}-2\nu {\hat{f}}_{i+1/2}^{n})+\displaystyle \frac{1}{2}({f}_{i}^{n}+2\nu {\hat{f}}_{i-1/2}^{n})\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray} & = & \displaystyle \frac{1}{2}\{{\theta }_{i+1/2}{U}_{i}^{+}+(1-{\theta }_{i+1/2}){U}_{i}^{\mathrm{UP},+}\}\\ & & +\displaystyle \frac{1}{2}\{{\theta }_{i-1/2}{U}_{i}^{-}+(1-{\theta }_{i-1/2}){U}_{i}^{\mathrm{UP},-}\},\end{eqnarray} \tag{ 19 }$

where ${U}_{i}^{\mathrm{UP},+}={f}_{i}^{n}-2\nu {f}_{i+1/2}^{\mathrm{UP},n}$ and ${U}_{i}^{\mathrm{UP},-}={f}_{i}^{n}+2\nu {f}_{i-1/2}^{\mathrm{UP},n}$. Therefore, if both of the terms on the right-hand side of Equation (19) are positive, ${f}_{i}^{(1)}$ is also positive. A simple calculation shows that it can be realized by setting the parameter ${\theta }_{i+1/2}$ as

$\begin{eqnarray}&&{\theta }_{i+1/2}=\min ({\theta }_{i+1/2}^{+},{\theta }_{i+1/2}^{-}),\end{eqnarray} \tag{ 20 }$

where

$\begin{eqnarray}{\theta }_{i+1/2}^{+}=\left\{\begin{array}{ll}\displaystyle \frac{{U}_{i}^{\mathrm{UP},+}}{{U}_{i}^{\mathrm{UP},+}-{U}_{i}^{+}} & \,\,\mathrm{if}\,{U}_{i}^{+}\lt 0\\ 1 & \,\,\mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 21 }$

and

$\begin{eqnarray}{\theta }_{i+1/2}^{-}=\left\{\begin{array}{ll}\displaystyle \frac{{U}_{i+1}^{\mathrm{UP},-}}{{U}_{i+1}^{\mathrm{UP},-}-{U}_{i+1}^{-}} & \,\,\mathrm{if}\,{U}_{i+1}^{-}\lt 0\\ 1 & \,\,\mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 22 }$

In addition, the condition of ${\theta }_{i+1/2}\geqslant 0$ can be met by setting ${U}_{i}^{\mathrm{UP},+}\geqslant 0$ and ${U}_{i+1}^{\mathrm{UP},-}\geqslant 0$ as a sufficient condition, which leads to a constraint on the CFL parameter as $\nu \leqslant 1/2$. Let us denote the above procedures to evaluate the positive interface value (15) as

$\begin{eqnarray}&&{\hat{f}}_{i+1/2}^{n}=\mathrm{PP}({\bar{f}}_{i+1/2}^{n},{f}_{i}^{n},{f}_{i+1}^{n}).\end{eqnarray} \tag{ 23 }$

It can be shown that positivity of the numerical solution is preserved even with the three-stage TVD–RK time integration by setting the CFL parameter to $\nu \lt 1/2$. In the following, we call this scheme the RK–MPP5 scheme (MPP for "monotonicity and positivity preserving").

It is tedious but straightforward to construct the spatially seventh-order schemes, RK–MP7 and RK–MPP7, by simply replacing ${f}_{i-1/2}^{\mathrm{int},5}$ in Equation (9) with the seventh-order interpolation

$\begin{eqnarray}{f}_{i-1/2}^{\mathrm{int},7} & = & (-3{f}_{i-4}+25{f}_{i-3}-101{f}_{i-2}+319{f}_{i-1}\\ & & +214{f}_{i}-38{f}_{i+1}+4{f}_{i+2})/420.\end{eqnarray} \tag{ 24 }$

2.1.2. Semi-Lagrangian Time Integration

The MP schemes described above adopt the three-stage TVD–RK time integration (10), which needs explicit Euler integration three times per single time step. It can be a source of significant computational cost. To overcome this drawback, we devise an alternative method based on the semi-Lagrangian (SL) time integration.

The averaged value of $f(x,t)$ over the ith mesh grid ${f}_{i}^{n+1}$ at a time of ${t}^{n+1}={t}^{n}+{\rm{\Delta }}t$ is written as

$\begin{eqnarray}{f}_{i}^{n+1} & = & \displaystyle \frac{1}{{\rm{\Delta }}x}{\displaystyle \int }_{{x}_{i-1/2}}^{{x}_{i+1/2}}f(x,{t}^{n+1}){dx}\\ & = & \displaystyle \frac{1}{{\rm{\Delta }}x}{\displaystyle \int }_{X({t}^{n},{t}^{n+1},{x}_{i-1/2})}^{X({t}^{n},{t}^{n+1},{x}_{i+1/2})}f(x,{t}^{n})\,{dx},\end{eqnarray} \tag{ 25 }$

where $X({t}_{1},{t}_{2},x)$ is the position of the characteristic curve at a time of t₁ originating from the phase space coordinate $(x,{t}_{2})$, and for the advection Equation (6), $X({t}_{1},{t}_{2},x)\,=x-c({t}_{2}-{t}_{1})$; see Figure 1. Denoting

$\begin{eqnarray}&&{{\rm{\Phi }}}_{i+1/2}^{n}=\displaystyle \frac{1}{c{\rm{\Delta }}t}{\int }_{X({t}^{n},{t}^{n+1},{x}_{i+1/2})}^{{x}_{i+1/2}}f(x,{t}^{n})\,{dx}\end{eqnarray} \tag{ 26 }$

and

$\begin{eqnarray}&&{{\rm{\Phi }}}_{i-1/2}^{n}=\displaystyle \frac{1}{c{\rm{\Delta }}t}{\int }_{X({t}^{n},{t}^{n+1},{x}_{i-1/2})}^{{x}_{i-1/2}}f(x,{t}^{n})\,{dx},\end{eqnarray} \tag{ 27 }$

we can rewrite Equation (25) in terms of numerical flux as

$\begin{eqnarray}&&{f}_{i}^{n+1}={f}_{i}^{n}-\nu ({{\rm{\Phi }}}_{i+1/2}^{n}-{{\rm{\Phi }}}_{i-1/2}^{n}).\end{eqnarray} \tag{ 28 }$

By interpolating the discrete values of fⁿ_j, we can compute the numerical fluxes ${{\rm{\Phi }}}_{i+1/2}$ and ${{\rm{\Phi }}}_{i-1/2}$. Actually, the PFC scheme also adopts essentially the same approach. It adopts the third-order reconstruction to interpolate fⁿ_j with a slope limiter and suppresses numerical oscillation while preserving positivity of the numerical solutions. For more detailed description of the PFC scheme, see Filbet et al. (2001) and the Appendix of YYU13.

Zoom In Zoom Out Reset image size

Our new schemes are based on the conservative semi-Lagrangian (CSL) scheme (Qiu & Christlieb 2010; Qiu & Shu 2011), which is originally coupled with a finite-difference scheme but is also applicable to the finite-volume schemes. In the fifth-order CSL scheme, ${{\rm{\Phi }}}_{i-1/2}^{n}$ can be expressed with five stencil values as

$\begin{eqnarray}&&{{\rm{\Phi }}}_{i-1/2}^{n}=\displaystyle \sum _{j=0}^{4}{C}_{j}{\zeta }^{j},\end{eqnarray} \tag{ 29 }$

where

$\begin{eqnarray}&&{C}_{0}=\displaystyle \frac{{f}_{i-3}^{n}}{30}-\displaystyle \frac{13}{60}{f}_{i-2}^{n}+\displaystyle \frac{47}{60}{f}_{i-1}^{n}+\displaystyle \frac{9}{20}{f}_{i}^{n}-\displaystyle \frac{{f}_{i+1}^{n}}{20}\end{eqnarray} \tag{ 30 }$

$\begin{eqnarray}&&{C}_{1}=-\displaystyle \frac{{f}_{i-2}^{n}}{24}+\displaystyle \frac{5}{8}{f}_{i-1}^{n}-\displaystyle \frac{5}{8}{f}_{i}^{n}+\displaystyle \frac{{f}_{i+1}^{n}}{24}\end{eqnarray} \tag{ 31 }$

$\begin{eqnarray}&&{C}_{2}=-\displaystyle \frac{{f}_{i-3}^{n}}{24}+\displaystyle \frac{{f}_{i-2}^{n}}{4}-\displaystyle \frac{{f}_{i-1}^{n}}{3}+\displaystyle \frac{{f}_{i}^{n}}{12}+\displaystyle \frac{{f}_{i+1}^{n}}{24}\end{eqnarray} \tag{ 32 }$

$\begin{eqnarray}&&{C}_{3}=\displaystyle \frac{{f}_{i-2}^{n}}{24}-\displaystyle \frac{{f}_{i-1}^{n}}{8}+\displaystyle \frac{{f}_{i}^{n}}{8}-\displaystyle \frac{{f}_{i+1}^{n}}{24}\end{eqnarray} \tag{ 33 }$

$\begin{eqnarray}&&{C}_{4}=\displaystyle \frac{{f}_{i-3}^{n}}{120}-\displaystyle \frac{{f}_{i-2}^{n}}{30}+\displaystyle \frac{{f}_{i-1}^{n}}{20}-\displaystyle \frac{{f}_{i}^{n}}{30}+\displaystyle \frac{{f}_{i+1}^{n}}{120}\end{eqnarray} \tag{ 34 }$

and

$\begin{eqnarray}&&\zeta =\displaystyle \frac{c{\rm{\Delta }}t}{{\rm{\Delta }}x}.\end{eqnarray} \tag{ 35 }$

The seventh-order CSL scheme is given by

$\begin{eqnarray}&&{{\rm{\Phi }}}_{i-1/2}^{n}=\displaystyle \sum _{j=0}^{6}{D}_{j}{\zeta }^{j},\end{eqnarray} \tag{ 36 }$

where

$\begin{eqnarray}{D}_{0} & = & -\displaystyle \frac{{f}_{i-4}^{n}}{140}+\displaystyle \frac{5}{84}{f}_{i-3}^{n}-\displaystyle \frac{101}{420}{f}_{i-2}^{n}+\displaystyle \frac{319}{420}{f}_{i-1}^{n}\\ & & +\displaystyle \frac{107}{210}{f}_{i}^{n}-\displaystyle \frac{19}{210}{f}_{i+1}^{n}+\displaystyle \frac{{f}_{i+2}^{n}}{105}\end{eqnarray} \tag{ 37 }$

$\begin{eqnarray}{D}_{1} & = & \displaystyle \frac{{f}_{i-3}^{n}}{180}-\displaystyle \frac{5}{72}{f}_{i-2}^{n}+\displaystyle \frac{49}{72}{f}_{i-1}^{n}\\ & & -\displaystyle \frac{49}{72}{f}_{i}^{n}+\displaystyle \frac{5}{72}{f}_{i+1}^{n}-\displaystyle \frac{{f}_{i+2}^{n}}{180}\end{eqnarray} \tag{ 38 }$

$\begin{eqnarray}{D}_{2} & = & \displaystyle \frac{7}{720}{f}_{i-4}^{n}-\displaystyle \frac{19}{240}{f}_{i-3}^{n}+\displaystyle \frac{7}{24}{f}_{i-2}^{n}\\ & & -\displaystyle \frac{23}{72}{f}_{i-1}^{n}+\displaystyle \frac{{f}_{i}^{n}}{48}+\displaystyle \frac{7}{80}{f}_{i+1}^{n}-\displaystyle \frac{{f}_{i+2}^{n}}{90}\end{eqnarray} \tag{ 39 }$

$\begin{eqnarray}{D}_{3} & = & -\displaystyle \frac{{f}_{i-3}^{n}}{144}+\displaystyle \frac{11}{144}{f}_{i-2}^{n}-\displaystyle \frac{7}{36}{f}_{i-1}^{n}\\ & & +\displaystyle \frac{7}{36}{f}_{i}^{n}-\displaystyle \frac{11}{144}{f}_{i+1}^{n}+\displaystyle \frac{{f}_{i+2}^{n}}{144}\end{eqnarray} \tag{ 40 }$

$\begin{eqnarray}{D}_{4} & = & -\displaystyle \frac{{f}_{i-4}^{n}}{360}+\displaystyle \frac{{f}_{i-3}^{n}}{48}-\displaystyle \frac{13}{240}{f}_{i-2}^{n}\\ & & +\displaystyle \frac{23}{360}{f}_{i-1}^{n}-\displaystyle \frac{{f}_{i}^{n}}{30}+\displaystyle \frac{{f}_{i+1}^{n}}{240}+\displaystyle \frac{{f}_{i+2}^{n}}{720}\end{eqnarray} \tag{ 41 }$

$\begin{eqnarray}&&{D}_{5}=\displaystyle \frac{{f}_{i-3}^{n}}{720}-\displaystyle \frac{{f}_{i-2}^{n}}{144}+\displaystyle \frac{{f}_{i-1}^{n}}{72}-\displaystyle \frac{{f}_{i}^{n}}{72}+\displaystyle \frac{{f}_{i+1}^{n}}{144}-\displaystyle \frac{{f}_{i+2}^{n}}{720}\end{eqnarray} \tag{ 42 }$

and

$\begin{eqnarray}{D}_{6} & = & \displaystyle \frac{{f}_{i-4}^{n}}{5040}-\displaystyle \frac{{f}_{i-3}^{n}}{840}+\displaystyle \frac{{f}_{i-2}^{n}}{336}-\displaystyle \frac{{f}_{i-1}^{n}}{252}\\ & & +\displaystyle \frac{{f}_{i}^{n}}{336}-\displaystyle \frac{{f}_{i+1}^{n}}{840}+\displaystyle \frac{{f}_{i+2}^{n}}{5040}.\end{eqnarray} \tag{ 43 }$

However, the numerical solution obtained with the CSL scheme does not preserve the monotonicity. Thus, we apply the MP constraint to ${{\rm{\Phi }}}_{i-1/2}^{n}$ and ${{\rm{\Phi }}}_{i+1/2}^{n}$ as

$\begin{eqnarray}&&{\bar{{\rm{\Phi }}}}_{i-1/2}^{n}=\mathrm{MP}({{\rm{\Phi }}}_{i-1/2}^{n},{f}_{i-3}^{n},{f}_{i-2}^{n},{f}_{i-1}^{n},{f}_{i}^{n},{f}_{i+1}^{n})\end{eqnarray} \tag{ 44 }$

and

$\begin{eqnarray}&&{\bar{{\rm{\Phi }}}}_{i+1/2}^{n}=\mathrm{MP}({{\rm{\Phi }}}_{i+1/2}^{n},{f}_{i-2}^{n},{f}_{i-1}^{n},{f}_{i}^{n},{f}_{i+1}^{n},{f}_{i+2}^{n}).\end{eqnarray} \tag{ 45 }$

The numerical solution with the monotonicity preservation at $t={t}^{n+1}$ is then computed as

$\begin{eqnarray}&&{f}_{i}^{n+1}={f}_{i}^{n}-\nu ({\bar{{\rm{\Phi }}}}_{i+1/2}^{n}-{\bar{{\rm{\Phi }}}}_{i-1/2}^{n}),\end{eqnarray} \tag{ 46 }$

which is referred to as the SL–MP5 and SL–MP7 schemes (SL for "semi-Lagrangian") for the fifth- and seventh-order accuracies, respectively.

The PP limiter described in the previous subsection can be applied as

$\begin{eqnarray}&&{\hat{{\rm{\Phi }}}}_{i-1/2}^{n}=\mathrm{PP}({\bar{{\rm{\Phi }}}}_{i-1/2}^{n},{f}_{i-1}^{n},{f}_{i}^{n})\end{eqnarray} \tag{ 47 }$

and

$\begin{eqnarray}&&{\hat{{\rm{\Phi }}}}_{i+1/2}^{n}=\mathrm{PP}({\bar{{\rm{\Phi }}}}_{i+1/2}^{n},{f}_{i}^{n},{f}_{i+1}^{n}).\end{eqnarray} \tag{ 48 }$

Then, the numerical solution at ${t}^{n+1}$ is given by

$\begin{eqnarray}&&{f}_{i}^{n+1}={f}_{i}^{n}-\nu ({\hat{{\rm{\Phi }}}}_{i+1/2}^{n}-{\hat{{\rm{\Phi }}}}_{i-1/2}^{n}).\end{eqnarray} \tag{ 49 }$

We refer to these schemes as the SL–MPP5 and SL–MPP7 schemes hereafter.

The Poisson Equation (3) is numerically solved with the convolution method (Hockney & Eastwood 1981) using the fast Fourier transform (FFT). Under periodic boundary conditions, the discrete Fourier transform (DFT) of the gravitational potential $\phi ({\boldsymbol{k}})$ can be written as

$\begin{eqnarray}&&\hat{\phi }({\boldsymbol{k}})=G({\boldsymbol{k}})\hat{\rho }({\boldsymbol{k}}),\end{eqnarray} \tag{ 50 }$

where $\hat{\rho }({\boldsymbol{k}})$ is the DFT of the density field, ${\boldsymbol{k}}=({k}_{x},{k}_{y},{k}_{z})$ is the discrete wavenumber vector, and $G({\boldsymbol{k}})$ is the Fourier-transformed Green function of the discretized Poisson equation. The second-order discretization of the Green function, $G({\boldsymbol{k}})$, is given by

$\begin{eqnarray}&&G({\boldsymbol{k}})=-\displaystyle \frac{\pi G{{\rm{\Delta }}}^{2}}{{\displaystyle \sum }_{i=x,y,z}{\sin }^{2}\left(\tfrac{{k}_{i}{\rm{\Delta }}}{2}\right)},\end{eqnarray} \tag{ 51 }$

with ${\rm{\Delta }}x={\rm{\Delta }}y={\rm{\Delta }}z={\rm{\Delta }}$. Since we adopt advection solvers with spatially fifth- and seventh-order accuracy, it is reasonable to adopt a Poisson solver with the same or similar order of accuracy, which can be simply constructed by high-order discretization of the Poisson equation. As shown in Appendix C, the DFTs of the Green functions with spatially fourth- and sixth-order accuracy are given by

$\begin{eqnarray}&&G({\boldsymbol{k}})=\displaystyle \frac{12\pi G{{\rm{\Delta }}}^{2}}{\displaystyle \sum _{i=x,y,z}\left[{\sin }^{2}({k}_{i}{\rm{\Delta }})-16{\sin }^{2}\left(\tfrac{{k}_{i}{\rm{\Delta }}}{2}\right)\right]}\end{eqnarray} \tag{ 52 }$

and

$\begin{eqnarray}&&G({\boldsymbol{k}})\,=\\ &&-\displaystyle \frac{180\pi G{{\rm{\Delta }}}^{2}}{\displaystyle \sum _{i=x,y,z}\left[2{\sin }^{2}\left(\tfrac{3{k}_{i}{\rm{\Delta }}}{2}\right)-27{\sin }^{2}({k}_{i}{\rm{\Delta }})+270{\sin }^{2}\left(\tfrac{{k}_{i}{\rm{\Delta }}}{2}\right)\right]},\end{eqnarray} \tag{ 53 }$

respectively. The gravitational potential $\phi ({\boldsymbol{x}})$ is computed with the inverse DFT of $\hat{\phi }({\boldsymbol{k}})$.

For isolated boundary conditions, the doubling-up method (Hockney & Eastwood 1981) is adopted, where the number of mesh grids is doubled for each spatial dimension and the mass density in the extended grid points is zeroed out. The Green function is set up in the real space as

$\begin{eqnarray}&&G(x,y,z)=\displaystyle \frac{G}{{({x}^{2}+{y}^{2}+{z}^{2})}^{1/2}}\end{eqnarray} \tag{ 54 }$

and is Fourier-transformed to obtain $\hat{G}({\boldsymbol{k}})$. Then, the gravitational potential can be calculated in the same manner as in the periodic boundary condition.

After obtaining the gravitational potential, we compute gravitational forces at each of the spatial grids with the finite-difference approximation (FDA). Again we adopt FDA schemes with the same order of accuracy as that of the gravitational potential calculation. The two-, four- and six-point FDAs are described in Appendix C. In what follows, we adopt the sixth-order calculation of the gravitational forces and potentials irrespective of adopted advection schemes because we are primarily interested in the accuracy of our advection equation solver.

We parallelize our Vlasov–Poisson solvers by decomposing the six-dimensional phase space into distributed memory spaces. We decompose the computational domain in the phase space in the same manner as YYU13. The physical space is divided into rectangular subdomains along x-, y-, and z-directions, while the velocity space is not decomposed; each spatial point has an entire set of velocity mesh grids. This decomposition is convenient when computing the velocity moments of the distribution function such as the mass density, mean velocities, and velocity dispersions.

In solving the advection equations in the physical space (4) across the decomposed subdomains, the values of the distribution function adjacent to the boundaries of the subdomains are transferred using the Message Passing Interface (MPI). The number of adjacent mesh grids to be transferred is two for the spatially third-order PFC scheme and three and four for spatially fifth- and seventh-order MP schemes, respectively.

We advance the distribution function from $t={t}^{n}$ to ${t}^{n+1}$ using the directional splitting method on the order of

$\begin{eqnarray}f({\boldsymbol{x}},{\boldsymbol{v}},{t}^{n+1}) & = & {T}_{{v}_{z}}({\rm{\Delta }}t/2){T}_{{v}_{y}}({\rm{\Delta }}t/2){T}_{{v}_{x}}({\rm{\Delta }}t/2)\\ & & \times {T}_{x}({\rm{\Delta }}t){T}_{y}({\rm{\Delta }}t){T}_{z}({\rm{\Delta }}t)\\ & & \times {T}_{{v}_{z}}({\rm{\Delta }}t/2){T}_{{v}_{y}}({\rm{\Delta }}t/2){T}_{{v}_{x}}({\rm{\Delta }}t/2)f({\boldsymbol{x}},{\boldsymbol{v}},{t}^{n}),\end{eqnarray} \tag{ 55 }$

where ${T}_{{\ell }}({\rm{\Delta }}t)$ denotes an operator for advection along the ℓ-direction with time step ${\rm{\Delta }}t$. It should be noted that the Poisson equation is solved after the advection operations in the physical space (${T}_{x}{T}_{y}{T}_{z}$) to update the gravitational potential and forces for the updated density field. We adopt the directional splitting method rather than the unsplit method because the former is the combination of one-dimensional advection schemes and hence rigorously assures the monotonicity and the positivity even in the multidimensional problems and also because the latter requires a larger amount of memory space to temporally store the numerical fluxes or the updated distribution function. Furthermore, the unsplit method requires a relatively smaller time step than the directional splitting method to achieve the same level of numerical accuracy. See Section 3.1.4 for the comparison of the dimensional splitting and unsplit methods in the two-dimensional advection problems. The temporal accuracy of the directional splitting method is in the second order at best. However, we find that the errors in the numerical solutions presented in the following numerical test are mostly dominated by ones that originate from the spatial discretization, and the spatial resolution is more important than the temporal accuracy in terms of the quality of the numerical solution.

As for time step width, we have several constraints on the CFL parameter. The MP schemes are subject to the constraint with $\nu \leqslant 1/5$ if the parameter α in the MP scheme is set to $\alpha =4$ (see Appendix A) and the PP limiter requires $\nu \leqslant 1/2$. As shown in Appendix B, all of our new schemes work well for $3\lt \alpha$ and $\nu \lt 0.4$. In practice, the condition $\nu =1/(1+\alpha )$ is desirable, and thus we set $\alpha =4$ and $\nu =0.2$ throughout the present paper unless otherwise stated.⁵

We set the time step for integrating the Vlasov equation as

$\begin{eqnarray}&&{\rm{\Delta }}t=\min ({\rm{\Delta }}{t}_{{\rm{p}}},{\rm{\Delta }}{t}_{{\rm{v}}}),\end{eqnarray} \tag{ 56 }$

where ${\rm{\Delta }}{t}_{{\rm{p}}}$ and ${\rm{\Delta }}{t}_{{\rm{v}}}$ are the time steps for advection equations in physical and velocity space, respectively. ${\rm{\Delta }}{t}_{{\rm{p}}}$ is given by

$\begin{eqnarray}&&{\rm{\Delta }}{t}_{{\rm{p}}}=\nu \min \left(\displaystyle \frac{{\rm{\Delta }}x}{{V}_{x}^{\max }},\displaystyle \frac{{\rm{\Delta }}y}{{V}_{y}^{\max }},\displaystyle \frac{{\rm{\Delta }}z}{{V}_{z}^{\max }}\right),\end{eqnarray} \tag{ 57 }$

where ${\rm{\Delta }}{\ell }$ denotes the spacing of mesh grids along the ℓ-direction and ${V}_{{\ell }}^{\max }$ is the maximum extent of the velocity space along the v_ℓ-direction. Similarly, we compute ${\rm{\Delta }}{t}_{{\rm{v}}}$ as

$\begin{eqnarray}&&{\rm{\Delta }}{t}_{{\rm{v}}}=\nu \,\mathop{\min }\limits_{i}\left(\displaystyle \frac{{\rm{\Delta }}{v}_{x}}{| {a}_{x,i}| },\displaystyle \frac{{\rm{\Delta }}{v}_{y}}{| {a}_{y,i}| },\displaystyle \frac{{\rm{\Delta }}{v}_{z}}{| {a}_{z,i}| }\right),\end{eqnarray} \tag{ 58 }$

where ${a}_{{\ell },i}$ is the gravitational acceleration along the ℓ-direction at the ith mesh grids in the physical space and the minimization is taken over all the spatial mesh grids.

As the first of the test suite of our new schemes, we perform numerical simulations of the advection equation. Here, we solve the one- and two-dimensional advection equations. In the one-dimensional tests, we examine the fundamental properties of our new schemes such as monotonicity, positivity, and accuracy. In the two-dimensional tests, we present the applicability of our one-dimensional schemes to multidimensional problems.

The first test run we perform is the one-dimensional advection problem. We solve

$\begin{eqnarray}&&\displaystyle \frac{\partial f(x,t)}{\partial t}+\displaystyle \frac{\partial f(x,t)}{\partial x}=0\end{eqnarray} \tag{ 59 }$

with our new schemes. Since the accuracy of the numerical solutions depends on the shape of $f(x,t)$, we consider three initial conditions and examine critically the numerical solutions.

3.1.1. Test-1a: Rectangular-shaped Function

We configure the rectangular-shaped initial condition given by

$\begin{eqnarray}f(x,t=0)=\left\{\begin{array}{ll}1 & 0.25\leqslant x\leqslant 0.75\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 60 }$

in the simulation domain over $0\leqslant x\leqslant 1$ and integrate Equation (59) using PFC, RK–MPP5, RK–MPP7, SL–MPP5, and SL–MPP7 under the periodic boundary conditions. Figure 2 shows the numerical solutions at t = 8 with the number of mesh grids ${N}_{{\rm{m}}}=64$. The numerical solutions with RK–MPP5, SL–MPP5, and SL–MPP7 are less diffusive compared with that of PFC. Note that RK–MPP7 yields undesirable asymmetric features around $x\simeq 0.2$ and 0.7. Actually these features are persistent irrespective of ${N}_{{\rm{m}}}$ and ν and are likely caused by the lack of accuracy in the time integration. In order to avoid these features, we would need the fourth- or higher-order TVD–RK integration schemes, which are computationally quite expensive (see Section 4). Thus, we exclude RK–MP7 and RK–MPP7 in the rest of our test simulations. Unlike RK–MPP7, SL–MPP7 does not exhibit such features, indicating that the SL time integration method is able to be combined with spatially seventh- or higher-order schemes. We confirm that the spatially ninth-order CSL scheme combined with the SL method works well at least in one-dimensional advection problems. Figure 3 shows the relative error of the numerical solutions to the analytic solution defined as

$\begin{eqnarray}&&{\epsilon }_{1{\rm{D}}}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{{N}_{{\rm{m}}}}| {f}_{i}-f({x}_{i},t)| }{{\displaystyle \sum }_{i=1}^{{N}_{{\rm{m}}}}| f({x}_{i},t)| }.\end{eqnarray} \tag{ 61 }$

Since there are discontinuities in the initial condition, the numerical fluxes computed with PFC and our new schemes are effectively of the first order, avoiding numerical oscillations. This is why the relative errors scale with the first-order power law.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.2. Test-1b: Sinusoidal Wave

As a good example of advection of a smooth function, we consider the initial condition given by

$\begin{eqnarray}&&f(x,t=0)=1+\sin (2\pi x),\end{eqnarray} \tag{ 62 }$

in the periodic domain of $0\leqslant x\leqslant 1$. Figure 4 shows the relative errors of numerical solutions obtained with RK–MPP5, SL–MPP5, and SL–MPP7 in the left panel and with RK–MP5, SL–MP5, and SL–MP7 without the PP limiter in the right panel as a function of ${N}_{{\rm{m}}}$. For comparison, the result obtained with PFC is also presented in both panels. The numerical solutions obtained with RK–MPP5 have nearly third-order accuracy for ${N}_{{\rm{m}}}\gtrsim {10}^{2}$, whereas those with RK–MP5 have fifth-order accuracy. This is because the PP limiter operates at the minima of numerical solutions at each stage of the TVD–RK time integration and degrades overall accuracy. Contrastingly, the accuracies of the SL schemes are not affected by the PP limiter, i.e., the results of SL–MPP5 and SL–MPP7 show fifth- and seventh-order scaling.

Zoom In Zoom Out Reset image size

3.1.3. Test-1c: Quartic Sine Function

The final example is the initial condition given by

$\begin{eqnarray}&&f(x,t=0)={\sin }^{4}(4\pi x).\end{eqnarray} \tag{ 63 }$

This case serves as a good example to examine how well the PP limiter works. Figure 5 shows the numerical solutions with ${N}_{{\rm{m}}}=64$ at t = 8 computed with PFC, RK–MP5, and RK–MPP5 (left panel) and those obtained with SL–MP5, SL–MPP5, and SL–MPP7 (right panel). The PFC solution is significantly smeared near the maxima and minima, although it preserves positivity. The numerical solutions obtained with RK–MP5 and SL–MP5 without the PP limiter exhibit noticeable negative values around the minima. On the other hand, RK–MPP5, SL–MPP5, and SL–MPP7 successfully preserve positivity and reproduce the analytic solution very well. Note that SL–MPP7 is better at reproducing the profiles around the maxima than SL–MPP5 and RK–MPP5, showing the importance of high-order accuracy. As can be seen in Figure 6, the order of accuracy for RK–MPP5, SL–MPP5, and SL–MPP7 is actually close to the fifth, and the result with SL–MPP7 is more accurate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.4. Test-1d: Two-dimensional Advection

Here we present numerical simulations of two-dimensional advection problems with our new schemes to test their numerical accuracy in multidimensional problems. We solve a two-dimensional advection equation

$\begin{eqnarray}&&\displaystyle \frac{\partial f({\boldsymbol{x}},t)}{\partial t}+{c}_{x}\displaystyle \frac{\partial f({\boldsymbol{x}},t)}{\partial x}+{c}_{y}\displaystyle \frac{\partial f({\boldsymbol{x}},t)}{\partial y}=0,\end{eqnarray} \tag{ 64 }$

where c_x and c_y are advection speed along the x- and y-direction, respectively. We set up a two-dimensional domain with $0\leqslant x,y\leqslant 1$ with a periodic boundary condition and consider three initial conditions similar to those presented in the one-dimensional tests: a rectangular-shaped wave given by

$\begin{eqnarray}f(x,y,t=0)=\left\{\begin{array}{cc}1 & 0.25\leqslant x,y\leqslant 0.75\\ 0 & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 65 }$

a sinusoidal wave

$\begin{eqnarray}&&f(x,y,t=0)=1+\displaystyle \frac{1}{2}\sin (2\pi x)+\displaystyle \frac{1}{2}\sin (2\pi y),\end{eqnarray} \tag{ 66 }$

and a quartic sine wave

$\begin{eqnarray}&&f(x,y,t=0)={\sin }^{4}(4\pi x)+{\sin }^{4}(4\pi y).\end{eqnarray} \tag{ 67 }$

In all three runs, we set the advection velocities to $({c}_{x},{c}_{y})=(1,0.5)$. The simulation domain is discretized into the Cartesian mesh grids with ${N}_{{\rm{x}}}$ and ${N}_{{\rm{y}}}$ mesh grids along the x- and y-directions, respectively.

We compute the relative error of numerical solutions to the analytic solution defined by

$\begin{eqnarray}&&{\epsilon }_{2{\rm{D}}}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{{N}_{{\rm{x}}}}{\displaystyle \sum }_{j=1}^{{N}_{{\rm{y}}}}| {f}_{i,j}-f({x}_{i},{y}_{j},t)| }{{\displaystyle \sum }_{i=1}^{{N}_{{\rm{x}}}}{\displaystyle \sum }_{j=1}^{{N}_{{\rm{y}}}}| f({x}_{i},{y}_{j},t)| },\end{eqnarray} \tag{ 68 }$

where ${f}_{i,j}$ is a numerical solution at $(x,y)=({x}_{i},{y}_{j})$ and x_i and y_j are the x and y coordinates of the ith and jth mesh grids in the x- and y-direction, respectively. Figure 7 shows the relative errors of numerical solutions obtained with PFC, SL–MPP5, and SL–MPP7 combined with the directional splitting or the unsplit methods, in which the numbers of mesh grids along each dimension are set to ${N}_{{\rm{x}}}={N}_{{\rm{y}}}=32$, 64, 128, 256, and 512.

Zoom In Zoom Out Reset image size

For the advection of a rectangular-shaped function, the scaling of the relative errors is nearly first order in space irrespective of adopted numerical schemes, while the results obtained with SL–MPP7 and SL–MPP5 are less diffusive than those obtained with PFC. On the other hand, numerical solutions of the advection of a two-dimensional sinusoidal wave obtained with PFC, SL–MPP5, and SL–MPP7 give the third-, fifth-, and seventh-order scalings of the relative error, respectively, in which the sinusoidal shape is so smooth that the MP constraint and PP limiter do not affect the overall accuracy of the schemes. In the case of a quartic sine function, the PP limiter compromises the numerical accuracy around the minima of the wave, and the numerical solutions with SL–MPP5 and SL–MPP7 exhibit nearly fifth-order accuracy, while the scaling of the relative errors obtained with the PFC scheme is not better than the third order. All these results are consistent with what we have in the one-dimensional Test-1a, Test-1b, and Test-1c, indicating that our new schemes also exhibit spatially high-order accuracies in the multidimensional problems.

As described in Section 2.4, we do not adopt the unsplit method for the time integration, due to its several disadvantages against the directional splitting method. It is, however, useful to make a comparison of numerical accuracy obtained with these two methods. The relative errors obtained with the unsplit method combined with PFC, SL–MPP5, and SL–MPP7 are also shown in each panel of Figure 7. We heuristically find that the CFL parameter must be less than 0.1 in the runs with the unsplit method to keep the monotonicity and positivity of the numerical solutions, and we set $\nu =0.02$ in computing the results with the unsplit method. It is clearly seen that the unsplit method yields almost the same level of numerical errors as the directional splitting method. We also find that the computational cost in computing with the unsplit method is larger than that of the directional splitting method by $\simeq 30$% in terms of wall-clock time per step.

Next, we perform Vlasov–Poisson simulations of spatially one-dimensional self-gravitating systems in two-dimensional phase space volume. We consider the periodic simulation domain with $0\leqslant x\leqslant L$ and set the initial condition to be

$\begin{eqnarray}&&f(x,v,t=0)=\displaystyle \frac{\bar{\rho }[1+A\cos ({kx})]}{{(2\pi {\sigma }^{2})}^{1/2}}\exp \left(-\displaystyle \frac{{v}^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 69 }$

where σ is the velocity dispersion, $\bar{\rho }$ is the mean density of the system, and A is the amplitude of density inhomogeneity. This distribution function yields the initial mass density as

$\begin{eqnarray}&&\rho (x)=\bar{\rho }[1+A\cos ({kx})],\end{eqnarray} \tag{ 70 }$

and the density inhomogeneity develops through the gravitational instability when the wavenumber k satisfies the condition $k/{k}_{{\rm{J}}}\lt 1$, where ${k}_{{\rm{J}}}$ is the Jeans wavenumber defined by

$\begin{eqnarray}&&{k}_{{\rm{J}}}={\left(\displaystyle \frac{4\pi G\bar{\rho }}{{\sigma }^{2}}\right)}^{1/2}.\end{eqnarray} \tag{ 71 }$

In this test, we set the wavenumber k to $k=4\pi /L$, the velocity dispersion σ so as to be $k/{k}_{{\rm{J}}}=0.5$, and the amplitude to A = 0.01. The velocity space is defined over $-V\leqslant v\leqslant V$, where V = L/T and T is the dynamical timescale defined by $T={(G\bar{\rho })}^{-1/2}$. We adopt outflow boundary conditions in the velocity space. Figure 8 shows the maps of the phase space density computed with PFC, RK–MPP5, SL–MPP5, and SL–MPP7 at $t=2T$, $3T$, and $4T$. The numbers of mesh grids in physical and velocity space are set to ${N}_{{\rm{x}}}={N}_{{\rm{v}}}=64$. The phase space density computed with PFC with third-order spatial accuracy is significantly smeared especially at later epochs compared to the results with other spatially higher-order schemes. The results with RK–MPP5 and SL–MPP5 schemes, both of which have the spatially fifth-order accuracy, appear almost identical to each other. SL–MPP7 is able to resolve finer spiral-arm structure of matter distribution in the phase space located at the turnaround position ($x\simeq 0.35$ and $x\simeq 0.65$).

Zoom In Zoom Out Reset image size

Figure 9 compares the results computed with RK–MPP5, SL–MPP5, and SL–MPP7 with the PP limiter (top panels) and with RK–MP5, SL–MP5, and SL–MP7 without it (bottom panels). The regions with negative phase space density depicted by white dots in the bottom panels are not present in the results computed with the PP limiter.

Zoom In Zoom Out Reset image size

To measure the accuracy of the numerical solutions obtained with various schemes, we compute the following quantity:

$\begin{eqnarray}&&{\epsilon }_{\mathrm{VP}}(t)=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{{N}_{{\rm{x}}}}{\displaystyle \sum }_{j=1}^{{N}_{{\rm{v}}}}| f({x}_{i},{v}_{j},t)-{f}_{\mathrm{ref}}({x}_{i},{v}_{j},t)| }{{\displaystyle \sum }_{i=1}^{{N}_{{\rm{x}}}}{\displaystyle \sum }_{j=1}^{{N}_{{\rm{v}}}}| {f}_{\mathrm{ref}}({x}_{i},{v}_{j},t)| },\end{eqnarray} \tag{ 72 }$

where ${f}_{\mathrm{ref}}(x,v,t)$ is a reference solution to be compared with. Since we do not have any analytical solution for this test, we construct the reference solution ${f}_{\mathrm{ref}}(x,v,t)$ in the following manner: First, we conduct runs with ${N}_{{\rm{x}}}={N}_{{\rm{v}}}=4096$ using SL–MPP7 and apply the fifth-order spline interpolation to the results so that the values of the numerical distribution function can be directly compared between different numbers of mesh grids. The left panel of Figure 10 shows the relative errors of the numerical solutions at $t=4T$ computed using PFC, RK–MPP5, SL–MPP5, and SL–MPP7 with ${N}_{{\rm{x}}}={N}_{{\rm{v}}}=64$–2048 mesh grids. The scalings of the relative errors of our new schemes and the PFC scheme are nearly in second and first order, respectively, and somewhat compromised compared with their formal orders of accuracies, because the functional shape of the distribution function is not so smooth that the MP constraint and the PP limiter effectively work. The accuracies of RK–MPP5, SL–MPP5, and SL–MPP7 for a given number of mesh grids are significantly better than PFC with the same number of mesh grids and are similar to or even better than that with four times as many mesh grids.

Zoom In Zoom Out Reset image size

We also measure the accuracies of our new schemes in the runs with $k/{k}_{{\rm{J}}}=1.1$, in which the density fluctuation damps through the collisionless damping or the free streaming. In the case with $k/{k}_{{\rm{J}}}=1.1$, we set the amplitude of density fluctuation to A = 0.1 and extend the velocity space as $-2V\leqslant v\leqslant 2V$ since the velocity dispersion is larger than the case with $k/{k}_{{\rm{J}}}=0.5$. We compute the relative errors of the distribution function at $t=4T$ in the same manner as done in the case with $k/{k}_{{\rm{J}}}=0.5$. The right panel of Figure 10 shows the relative errors of numerical solutions computed with our new schemes as well as PFC for ${N}_{{\rm{x}}}={N}_{{\rm{v}}}=64$–2048 mesh grids. In this test, the shape of the distribution function is relatively smoother than the runs with $k/{k}_{{\rm{J}}}=0.5$, and the numerical solutions are not as strongly affected by the MP constraint and the PP limiter as the runs with $k/{k}_{{\rm{J}}}=0.5$. Accordingly, the scalings of the relative errors are better than those in the case with $k/{k}_{{\rm{J}}}=0.5$ and are nearly third and second order for our new schemes and PFC, respectively. Note that the accuracies of RK–MPP5 and SL–MPP5 are almost the same as each other in both cases with $k/{k}_{{\rm{J}}}=0.5$ and 1.1, indicating that the RK schemes and SL schemes yield almost the same accuracy in these tests.

We perform simulations of Landau damping in an electrostatic plasma, in which the time evolution of the distribution function of electrons is described by the following Vlasov equation:

$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+v\cdot \displaystyle \frac{\partial f}{\partial x}-\displaystyle \frac{e}{m}\displaystyle \frac{\partial {\phi }_{{\rm{s}}}}{\partial x}\cdot \displaystyle \frac{\partial f}{\partial v}=0,\end{eqnarray} \tag{ 73 }$

where ${\phi }_{{\rm{s}}}$ is the electrostatic potential and e and m are the charge and mass of an electron, respectively. The electrostatic potential satisfies the Poisson equation

$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\phi }_{{\rm{s}}}=4\pi {ne}\left(1-{\int }_{-\infty }^{\infty }{fdv}\right),\end{eqnarray} \tag{ 74 }$

where n is the number density of ions. In this formulation, the distribution function is normalized such that its integration over the velocity space gives the fraction of electron number density relative to ions.

We consider a two-dimensional phase space with $-2\pi L\leqslant x\leqslant 2\pi L$ and $-6V\leqslant v\leqslant 6V$, where $V\equiv {\omega }_{\mathrm{pe}}L$ and ${\omega }_{\mathrm{pe}}={(4\pi {e}^{2}n/m)}^{1/2}$ is the electron plasma frequency and the periodic boundary condition is imposed on the physical space. The initial condition is given by

$\begin{eqnarray}&&f(x,v,t=0)=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^{2}}}\exp \left(-\displaystyle \frac{{v}^{2}}{2{\sigma }^{2}}\right)(1+A\cos {kx}),\end{eqnarray} \tag{ 75 }$

where the wavenumber is set to $k=1/(2L)$ and $\sigma =V$. With this initial condition, it is expected that the electric field damps exponentially as $\exp (-\gamma {\omega }_{\mathrm{pe}}t)$ with $\gamma =0.153$.

First, we show results of linear Landau damping in which we set the amplitude of the initial perturbation to A = 0.01. We set ${N}_{{\rm{x}}}={N}_{{\rm{v}}}=64$. Figure 11 shows time evolution of the amplitude of the electric field with a mode of $k=1/(2L)$ obtained with PFC and SL–MPP5. The results obtained with RK–MPP5 and SL–MPP7 are almost identical to that of SL–MPP5 and are not presented in this figure. Although the numerical solution with PFC ceases to damp at $t\gtrsim 40{\omega }_{\mathrm{pe}}^{-1}$, the one obtained with the SL–MPP5 scheme continues to damp to $t\lesssim 55{\omega }_{\mathrm{pe}}^{-1}$, just before the recurrence time ${t}_{{\rm{R}}}=2\pi /(k{\rm{\Delta }}v)=67{\omega }_{\mathrm{pe}}^{-1}$, indicating that SL–MPP5 is less dissipative than PFC.

Zoom In Zoom Out Reset image size

Next, we perform simulations of nonlinear Landau damping with A = 0.5. This problem has been investigated in a number of previous studies (e.g., Klimas 1987; Manfredi 1997; Nakamura & Yabe 1999). In this problem, the initial perturbations damp (in space) quickly, but a highly stratified distribution in the phase space develops with time through rapid phase mixing (see Figure 12). The amplitude of the electric field initially decreases, but particle trapping by the electric field prevents complete damping. Figure 12 shows the phase space density of electrons at $t=10\,{\omega }_{\mathrm{pe}}^{-1}$, $20\,{\omega }_{\mathrm{pe}}^{-1}$, $40\,{\omega }_{\mathrm{pe}}^{-1}$, and $50\,{\omega }_{\mathrm{pe}}^{-1}$ computed with PFC, RK–MPP5, SL–MPP5, and SL–MPP7. All the results are obtained with ${N}_{{\rm{x}}}=64$ and ${N}_{{\rm{v}}}=256$. Rapid phase space mixing generates a number of peaks along the velocity axis. The PFC solution exhibits a smeared distribution at $t\gtrsim 40\,{\omega }_{\mathrm{pe}}^{-1}$, while those computed with other higher-order schemes preserve features of phase mixing even at $t=50\,{\omega }_{\mathrm{pe}}^{-1}$. This is also seen in Figure 13, in which we show the velocity distribution function

$\begin{eqnarray}&&F(v,t)=\int f(x,v,t){dx}\end{eqnarray} \tag{ 76 }$

at $t=50\,{\omega }_{\mathrm{pe}}^{-1}$. There are a number of peaks in the velocity distribution function, which represent electrons trapped by electrostatic waves around their phase velocities. The peaks in the distribution function are more prominent in the results with fifth- and seventh-order schemes than in PFC, again showing the less dissipative nature of our schemes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we study gravitational collapse of a uniform-density sphere. The mass and the initial radius of the sphere are M and $2R$, respectively, and the initial distribution function is given by

$\begin{eqnarray}&&f({\boldsymbol{x}},{\boldsymbol{v}},t=0)=\displaystyle \frac{{\rho }_{0}}{{(2\pi {\sigma }^{2})}^{3/2}}\exp \left(-\displaystyle \frac{| {\boldsymbol{v}}{| }^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 77 }$

where ${\rho }_{0}=M/(4\pi {(2R)}^{3}/3)$ is the initial density of the sphere and ${\sigma }^{2}={GM}/(20R)$ is the initial velocity dispersion of the matter inside the sphere. We set the initial virial ratio to be 0.5. The extent of the simulation volume is $-3R\leqslant x,y,z\leqslant 3R$ in the physical space and $-3V\leqslant {v}_{x},{v}_{y},{v}_{z}\leqslant 3V$, where V is given by V = R/T and $T={({R}^{3}/{GM})}^{1/2}$ is the dynamical timescale of the system. In this test, we adopt the out-going boundary condition in both the physical and velocity spaces. We perform the runs with ${N}_{{\rm{x}}}={64}^{3}$ for the physical space and ${N}_{{\rm{v}}}={32}^{3}$ or ${N}_{{\rm{v}}}={64}^{3}$ mesh grids for the velocity space.

We plot the phase space density in the $(x,{v}_{x})$-plane with $y=z=0$ and ${v}_{y}={v}_{z}=0$ at $t=3T$, $4T$, $6T$, and $8T$ in Figure 14. We set ${N}_{{\rm{x}}}={N}_{{\rm{v}}}={64}^{3}$. Initially, the sphere starts to collapse toward the center and reaches the maximum contraction around $t=4T$. Then the system expands again, and a part of the matter passes through the center of the system to reach the outer boundary of the physical space defined at $x=\pm 3R$, $y=\pm 3R$, and $z=\pm 3R$. The results computed with PFC show a smeared distribution compared with those obtained with RK–MPP5, SL–MPP5, and SL–MPP7. We note that RK–MPP5 and SL–MPP5 yield almost identical results, although the computational costs for the two are significantly different (see Section 4).

Zoom In Zoom Out Reset image size

Figure 15 is the same plot as Figure 14 but shows the results with ${N}_{{\rm{v}}}={32}^{3}$. It is interesting that the results obtained by SL–MPP7 with ${N}_{{\rm{x}}}={64}^{3}$ and ${N}_{{\rm{v}}}={32}^{3}$ (bottom panels in Figure 15) are similar to or even better than those obtained with PFC with ${N}_{{\rm{x}}}={N}_{{\rm{v}}}={64}^{3}$ (top panels in Figure 14), despite the eight times smaller number of mesh grids.

Zoom In Zoom Out Reset image size

For this test, we compare our new schemes with another independent scheme. To this end, we utilize a high-resolution N-body simulation in the following manner. First, we perform an N-body simulation starting from the initial condition equivalent to Equation (77). In this N-body simulation, the gravitational potential field is dumped and stored at all the time steps. Then, we compute the phase space density at a given phase space coordinate $({\boldsymbol{x}}(t),{\boldsymbol{v}}(t))$ by tracing back the coordinate to its initial coordinate $({\boldsymbol{x}}(0),{\boldsymbol{v}}(0))$ using the previously stored gravitational potentials. Since the phase space density is constant along each trajectory, we can, in principle, compute the phase space density from the initial phase space coordinate $({\boldsymbol{x}}(0),{\boldsymbol{v}}(0))$ using Equation (77). However, trajectories in conventional N-body simulations are subject to rather grainy gravitational potential compared to what it should be in the actual physical system owing to their superparticle approximation. Thus, we adopt the self-consistent field (SCF) method (Hernquist & Ostriker 1992; Hozumi 1997), in which the gravitational potential is computed by expanding the density field into a set of basis functions, and particle trajectories are calculated for the smooth gravitational potential. The initial condition is constructed with 10⁷ equal-mass particles, and the particle orbits are computed using the basis set constructed by Clutton-Brock (1972), where the numbers of the expansion terms are 64 and 144 in the radial and angular directions, respectively. Figure 16 shows the phase space density obtained with this procedure. Comparing with Figures 14 and 15 manifests that our spatially higher-order schemes (RK–MPP5, SL–MPP5, and SL–MPP7) yield results much closer to those obtained with the SCF approach. It should be noted that, in practice, the SCF method is applicable only to self-gravitating systems with a certain spatial symmetry, such as in the present test, to which the lowest-order functions of the basis set suffice to describe the overall structure of the system. Also, storing the gravitational potential (or its expansion) at all the time steps is very costly. Therefore, the SCF approach cannot be applied to a wide range of self-gravitating systems. For further comparison with the results obtained by the SCF methods, we perform a Vlasov–Poisson run with the largest numbers of mesh grids, ${N}_{{\rm{x}}}={128}^{3}$ and ${N}_{{\rm{v}}}={128}^{3}$ using SL–MPP5, and show the phase space density at $t=3T$, $4T$, $6T$, and $8T$ in the $(x,{v}_{x})$-plane in Figure 17. With such high resolution, the phase space density is consistent with that obtained with the SCF method in detail.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We further examine the mass and energy conservation for this spherical collapse test. We show the time evolution of the total mass, kinetic energy, and gravitational potential energy of the results simulated with PFC and SL–MPP5 in Figure 18, where the results obtained with ${N}_{{\rm{v}}}={64}^{3}$ and ${N}_{{\rm{v}}}={32}^{3}$ are shown in the left and right panels, respectively. Note that the results simulated with RK–MPP5 and SL–MPP7 are almost the same as those of SL–MPP5 and hence are not shown in this figure. For comparison, we also plot the results obtained with an N-body simulation, in which we employ 128³ particles, and adopt the Particle-Mesh (PM) method with the number of mesh grids in solving the Poisson equation equal to that of the Vlasov–Poisson simulation, i.e., 64 mesh grids in each dimension, in order to realize the same spatial resolution of the gravitational field and force. The total mass in the Vlasov simulation starts to decrease around $t\simeq 8T$ because part of the matter reaches the boundary of the simulation volume at that time, as can be seen in Figures 16 and 17. In the runs with PFC, such a decrease of the total mass begins earlier and the amount of the escaped mass is larger than that in the runs with SL–MPP5 and other higher-order schemes owing to larger numerical diffusion, irrespective of ${N}_{{\rm{v}}}$. The most prominent difference can be seen around the maximum contraction epoch $t\simeq 4T$, when the kinetic and gravitational potential energies have their maximum and minimum values, respectively. As can be seen in Figures 14 and 15, the large numerical diffusion of PFC, in both the physical and velocity spaces, smears the "S"-shaped phase space distribution and enlarges the physical extent of the system at $t=4T$, yielding small kinetic and shallow gravitational potential energies seen in Figure 18. The overall accuracy can be improved for PFC by simply increasing the number of mesh grids from ${N}_{{\rm{v}}}={32}^{3}$ to 64³. On the other hand, the result of SL–MPP5 is close to the N-body simulation, irrespective of ${N}_{{\rm{v}}}$, and in the runs with ${N}_{{\rm{v}}}={64}^{3}$ both of the results with PFC and SL–MPP5 are in good agreement with the N-body simulation. This suggests that SL–MPP5 (as well as RK–MPP5 and SL–MPP7) with ${N}_{{\rm{v}}}={32}^{3}$ provides roughly the same level of accuracy as PFC with ${N}_{{\rm{v}}}={64}^{3}$ in terms of mechanical energy.

Zoom In Zoom Out Reset image size