Exponential Runge-Kutta methods for Vlasov-Poisson equation

Vlasov equation plays an important role in the modeling of plasma physics, it is a non-linear transport in $(x,v)$ of a density distribution $f = f(t,x,v)$; coupled with Poisson equation for estimate electric field $E=E(t,x)$:

High order methods are needed in phase space $(x,v)$ to study filamentation.

We study the stability problem on a simplify model: $$ \partial_tu + a\partial_xu = 0. $$ We discretize in space with WENO method [1], and we test various Runge-Kutta methods for the time integrator. We apply von Neumann analysis on space discretization to study amplification factor as a curve in complex plane $\mathbb{C}$.

We define stability domain of an explicit Runge-Kutta method RK($s$,$n$) as: $$ \mathcal{D}_{(s,n)} = \left\{ z\in\mathbb{C} : |p_{(s,n)}(z)|\leq 1 \right\} $$ with $s$ the number of stages, $n$ the order and $p_{(s,n)}$ the stability function of RK($s$,$n$).

It is possible to interpret CFL number between time integrator and space method as the biggest homothety ratio that wedges all the amplification factor curve into the stability domain of considered Runge-Kutta methods. We estimate numerically the CFL $\sigma$ (respectively $\sigma^\ell$) for various couples RK($s$,$n$) with WENO (respectively linearized WENO). All CFL have been approved with a long time simulation of ∗.

We apply a Fourier transform in $x$ direction of VP: $$\partial_t\hat{f} + i\kappa v\hat{f} + \widehat{(E\partial_vf)} = 0$$ of the form $$\partial_t u + Lu + N(u) = 0.$$ Next we use an exponential integrator (Duhamel formula): $$\partial_t(e^{i\kappa vt}\hat{f}) + e^{i\kappa vt}\widehat{(E\partial_vf)} = 0$$ of the form $$\partial_t(e^{Lt}u) + e^{Lt}N(u) = 0.$$ With this exponential form we can use a classical Runge-Kutta method, in this context we call them Lawson methods or Integrating factor Runge-Kutta methods [4].

Stability function of Lawson schemes can be expressed in terms of the underlying Runge-Kutta method: $$ p_{\text{Lawson}(s,n)}(z) = p_{(s,n)}(z)e^{L\Delta t} $$ with $L=i\kappa v\in i\mathbb{R}$ so stability domain is the same. Our numerical CFL study, on simplified model, still works on Duhamel formulation of Vlasov-Poisson.

We are interested in the numerical cost $\frac{\Delta t}{s}$ of RK($s$,$n$). To compare each time integrator, we compute total energy in Vlasov-Poisson system: $$ H(t) = \int_{\Omega}\int_{\mathbb{R}}v^2f\,\mathrm{d}v\mathrm{d}x + \int_{\Omega}E^2\,\mathrm{d}x $$ which is preserved in time. We propose to select the best method by considering:

$$ h_{(s,n)}:\frac{\Delta t}{s} \mapsto \left\| \frac{H(t)-H(0)}{H(0)} \right\|_{\infty}\!\! $$ Tested strategy:

• Spectral scheme in $x$,
• WENO method in $v$,
• Lawson method of the underlying RK($s$,$n$) in time.

• Promising approach for Vlasov-Poisson simulation, less stages than splitting strategy.
• Selected Runge-Kutta method is Dormand-Prince RK(6,5), this is an embedded Runge-Kutta method (adaptative stepsize, relaxation of the CFL).
• Extension to multi-dimensional Vlasov-Poisson is easier than splitting strategy.

Simulation of the bump on tail test case : $$f(t=0,x,v) = \left(\frac{0.9}{\sqrt{2\cdot\pi}}e^{-\frac{|v|^2}{2}} + \frac{0.2}{\sqrt{2\cdot\pi}}e^{-2|v-4.5|^2}\right)\left(1+0.04\cdot\cos(0.3\cdot x)\right)$$