2 Equations and numerical methods

In astro- and plasma physics the behaviour of a collisionless plasma is modelled by the Vlasov equation

(2.1) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}f(t,x,v)+v\boldsymbol{\cdot }\unicode[STIX]{x1D735}f(t,x,v)+F\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}f(t,x,v)=0,\end{eqnarray}$$

which is posed in an up to $3+3$ -dimensional phase space (although lower-dimensional models are often employed as well). The variable $x$ denotes the position and variable $v$ denotes the velocity. The density function $f$ is the sought-after particle distribution and the (force) term $F$ describes the interaction of the plasma with the electromagnetic field. Depending on the application this force term can include the full Lorentz force or only the force due to the electric field.

In this paper we will restrict ourselves to the Vlasov–Poisson equation, i.e. the force term is given by $F=E$ and the electric field is determined by solving

(2.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }E=\unicode[STIX]{x1D70C}(x)-1,\quad \unicode[STIX]{x1D735}\times E=0\end{eqnarray}$$

with charge density

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(x)=\int f(t,x,v)\,\text{d}v.\end{eqnarray}$$

In addition, in all our simulations we impose periodic boundary conditions in both the $x$ - and $v$ -directions.

To advance the numerical solution in time we use the splitting approach introduced in Cheng & Knorr (Reference Cheng and Knorr1976). That is, we consider (the first part of the splitting algorithm)

(2.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}f(t,x,v)+v\boldsymbol{\cdot }\unicode[STIX]{x1D735}f(t,x,v)=0,\quad f(0,x,v)=g(x,v)\end{eqnarray}$$

and denote the corresponding solution at time $\unicode[STIX]{x1D70F}$ by

(2.5) $$\begin{eqnarray}\text{e}^{\unicode[STIX]{x1D70F}A}g(x,v)=f(\unicode[STIX]{x1D70F},x,v).\end{eqnarray}$$

For the second part of the splitting algorithm we then consider

(2.6a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}f(t,x,v)+E^{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v}f(t,x,v)=0,\quad f(0,x,v)=g(x,v),\end{eqnarray}$$

where the electric field $E^{n}$ is determined from the density

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D70C}^{n}(x)=\int g(x,v)\,\text{d}v\end{eqnarray}$$

and is thus taken constant during each time step. We denote the corresponding solution by

(2.8) $$\begin{eqnarray}\text{e}^{\unicode[STIX]{x1D70F}B}g(x,v)=f(\unicode[STIX]{x1D70F},x,v).\end{eqnarray}$$

Note that, contrary to what the notation suggests, $B$ is a nonlinear operator as the electric field depends on the density function $f$ . Using the introduced notation we can easily formulate a time step of the splitting algorithm. The Lie splitting

(2.9) $$\begin{eqnarray}f^{n+1}(x,v)=\text{e}^{\unicode[STIX]{x1D70F}B}\text{e}^{\unicode[STIX]{x1D70F}A}f^{n}(x,v),\end{eqnarray}$$

where $f^{n}$ is an approximation of $f(t_{n},x,v)$ and $\unicode[STIX]{x1D70F}=t_{n+1}-t_{n}$ is the time step size, is a first-order method. In all our simulations we will use the Strang splitting scheme

(2.10) $$\begin{eqnarray}f^{n+1}(x,v)=\text{e}^{(\unicode[STIX]{x1D70F}/2)A}\text{e}^{\unicode[STIX]{x1D70F}B}\text{e}^{(\unicode[STIX]{x1D70F}/2)A}f^{n}(x,v),\end{eqnarray}$$

which is a second-order method and has approximately the same computational cost as the Lie splitting scheme. Let us also note that high-order splitting methods have been constructed (Crouseilles, Faou & Mehrenberger Reference Crouseilles, Faou and Mehrenberger2011a ).

The main computational advantage of the splitting scheme is that the resulting substeps, given by (2.4) and (2.6), are in the form of an advection equation, where the advection speed is independent of the variable in the direction of the advection. Thus, we can immediately solve (2.4) to obtain

(2.11) $$\begin{eqnarray}\text{e}^{\unicode[STIX]{x1D70F}A}g(x,v)=g(x-\unicode[STIX]{x1D70F}v,v)\end{eqnarray}$$

and (2.6) to obtain

(2.12) $$\begin{eqnarray}\text{e}^{\unicode[STIX]{x1D70F}B}g(x,v)=g(x,v-\unicode[STIX]{x1D70F}E(x)),\end{eqnarray}$$

where $E$ is the electric field determined from the charge density corresponding to $g(x,v)$ . Thus, by using the splitting approach outlined above we have reduced the task of computing a numerical approximation to the Vlasov–Poisson equation to computing two translations in phase space.

Now, up to this point phase space is still continuous. However, in order to implement the numerical scheme on a computer, we have to perform a space discretization that represents the numerical solution, at a fixed time step, using a finite number of degrees of freedom. The most straightforward approach is to choose a uniform grid. However, for any grid $(x_{i},v_{j})$ the translations $x_{i}-\unicode[STIX]{x1D70F}v_{j}$ and $v_{j}-\unicode[STIX]{x1D70F}E(x_{i})$ , in general, do not coincide with a grid point. Thus, we have to use an interpolation scheme in order to evaluate $g(x_{i}-\unicode[STIX]{x1D70F}v_{j},v_{j})$ and $g(x_{i},v_{j}-\unicode[STIX]{x1D70F}E(x_{i}))$ . In the Vlasov community interpolation based on cubic splines is very popular. This method is illustrated in figure 1 and proceeds in two steps:

1. (i) construct a cubic spline using the known value of the function on the grid points (this involves the solution of a tridiagonal system of equations);

2. (ii) translate each grid point according to (2.11) or (2.12) and use the cubic spline to obtain the new value at each grid point.

Since we follow the characteristics backward in time (i.e. track particles along their trajectories) but use an Eulerian grid for the space discretization, these methods are usually referred to as semi-Lagrangian and are free of a Courant–Friedrichs–Lewy condition. In addition, they do not exhibit the numerical noise that is prevalent in particle methods (such as the popular particle in cell scheme).

Figure 1.

An illustration of the semi-Lagrangian method based on a quadratic spline interpolation is shown. In this case a continuous piecewise polynomial interpolant is constructed that is used for all evaluations. We illustrate three evaluations of the interpolant (denoted by the variable $x$ ).

For the Vlasov–Poisson equation it is vital that charge (or mass; the two quantities are equivalent as all particles in the physical model carry the same amount of electric charge) is conserved. It is, a priori, not clear that the semi-Lagrangian scheme based on spline interpolation satisfies this constraint. In order to show this, we reformulate the numerical method as a forward semi-Lagrangian scheme (as is done in Crouseilles, Respaud & Sonnendrücker (Reference Crouseilles, Respaud and Sonnendrücker2009b )). Since our interpolant is a cubic spline, we can expand $f^{n}(x,v)$ in the B-spline basis (for simplicity we restrict ourselves to the $1+1$ -dimensional case here)

(2.13) $$\begin{eqnarray}f^{n}(x,v)=\mathop{\sum }_{kl}\unicode[STIX]{x1D714}_{kl}^{n}S(x-x_{k})S(v-v_{l}),\end{eqnarray}$$

where $S$ is given by

(2.14) $$\begin{eqnarray}6S(x)=\left\{\begin{array}{@{}ll@{}}4-6(x/h)^{2}+3|x/h|^{3}\quad & 0\leqslant |x|\leqslant h\\ (2-|x/h|)^{3}\quad & h\leqslant |x|\leqslant 2h\\ 0\quad & \text{otherwise}\end{array}\right.\end{eqnarray}$$

and $h$ is the grid spacing. The coefficients $\unicode[STIX]{x1D714}_{kl}^{n}$ are determined uniquely by the function values at the grid points but in order to obtain them we have to solve a linear system of equations. The forward semi-Lagrangian method based on cubic spline interpolation for (2.11) is then given by

(2.15) $$\begin{eqnarray}f^{n+1}(x_{i},v_{j})=\mathop{\sum }_{kl}\unicode[STIX]{x1D714}_{kl}^{n}S(x_{i}-\unicode[STIX]{x1D70F}v_{j}-x_{k})S(v_{j}-v_{l}).\end{eqnarray}$$

For the constant advection case considered here this is merely a reformulation of the previously introduced (backward) semi-Lagrangian method. However, the interpretation is different (translating the basis function versus following the characteristics backward in time). That the forward semi-Lagrangian scheme is mass conservative can be easily deduced from the present formulation. First, we note that the B-splines form a partition of unity. That is, for all $x$ we have

(2.16) $$\begin{eqnarray}\mathop{\sum }_{k}S(x-x_{k})=1.\end{eqnarray}$$

Therefore, we can write the integral over a spline $u$ as follows

(2.17) $$\begin{eqnarray}\int u(x)\,\text{d}x=\mathop{\sum }_{k}\unicode[STIX]{x1D714}_{k}\int S(x-x_{k})\,\text{d}x=h\mathop{\sum }_{k}\unicode[STIX]{x1D714}_{k}=h\mathop{\sum }_{i}\mathop{\sum }_{k}\unicode[STIX]{x1D714}_{k}S(x_{i}-x_{k})=h\mathop{\sum }_{i}u(x_{i}),\end{eqnarray}$$

which shows that the discrete charge is in fact equal to the continuous charge given by the cubic spline interpolation. Finally, we have

(2.18) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{i}u(x_{i}-\unicode[STIX]{x1D70F}v) & = & \displaystyle \mathop{\sum }_{i}\mathop{\sum }_{k}\unicode[STIX]{x1D714}_{k}S(x_{i}-\unicode[STIX]{x1D70F}v-x_{k})=\mathop{\sum }_{k}\unicode[STIX]{x1D714}_{k}\mathop{\sum }_{i}S(x_{i}-\unicode[STIX]{x1D70F}v-x_{k})\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{k}\unicode[STIX]{x1D714}_{k}=\mathop{\sum }_{i}u(x_{i}),\end{eqnarray}$$

which shows that the semi-Lagrangian scheme based on the cubic spline interpolation is charge conservative (up to machine precision). We further remark that, assuming a sufficient regular solution, this method is fourth-order accurate in space.

We now turn our attention to the semi-Lagrangian discontinuous Galerkin scheme. To that end we divide our domain into cells $C_{i}=I_{i}^{(1)}\times \cdots \times I_{i}^{(d)}$ , where the $I_{i}$ are one-dimensional intervals of length $h$ , and assume that a function $u$ is given such that $u|_{C_{i}}$ , i.e. the restriction of $u$ to the $i$ th cell, is a polynomial of degree $k$ . Then the function $u$ lies in the approximation space (note that we do not enforce a continuity constraint across cell interfaces). However, in general, this is not true for the translated function

(2.19) $$\begin{eqnarray}\displaystyle T_{\unicode[STIX]{x1D70F}}u(x)=u(x-\unicode[STIX]{x1D70F}v), & & \displaystyle\end{eqnarray}$$

where we have introduced the translation operator $T_{\unicode[STIX]{x1D70F}}$ . Thus, we perform an approximation by applying a projection operator $P$ to obtain

(2.20) $$\begin{eqnarray}u^{n+1}=PT_{t}u^{n}.\end{eqnarray}$$

Note that except for this projection no approximation has been made so far. The function $u^{n+1}$ constitutes the sought-after approximation of $u^{n}(x-\unicode[STIX]{x1D70F}v)$ . The operator $P$ is the $L^{2}$ projection on the (finite-dimensional) subspace of cell-wise polynomials of degree $\ell$ . This approach is illustrated in figure 2.

Figure 2.

Illustration of the semi-Lagrangian discontinuous Galerkin approximation in a single space dimension. The discontinuous piecewise polynomial approximation is translated in space and then projected back to the finite-dimensional subspace of piecewise polynomial functions.

In order to implement the numerical scheme, we still have to choose a basis for the space of cell-wise polynomials of degree $\ell$ . This then also determines the degrees of freedom stored in computer memory. Our implementation is based on the Lagrange basis at the Gauss–Legendre quadrature nodes. Thus, the degrees of freedom correspond to the value of the function at the Gauss–Legendre nodes in each cell. However, in certain situations, it is more convenient to express the piecewise polynomial function in the Legendre basis. Thus, we can write (for an $L^{2}$ function $u$ )

(2.21) $$\begin{eqnarray}u(x)\approx Pu(x)=\mathop{\sum }_{i=0}^{N-1}\mathop{\sum }_{j=0}^{\ell }u_{ij}p_{ij}(x),\end{eqnarray}$$

where $p_{ij}$ is the (appropriately scaled) $j$ th Legendre polynomial in the $i$ th cell and $N$ is the number of cells. The coefficients in the expansion are then given by

(2.22) $$\begin{eqnarray}u_{ij}=\int u(x)p_{ij}(x)\,\text{d}x.\end{eqnarray}$$

Since the Legendre polynomials form an orthogonal basis, we have

(2.23) $$\begin{eqnarray}\int u(x)\,\text{d}x=\mathop{\sum }_{i}u_{i0}=\int Pu(x)\,\text{d}x\end{eqnarray}$$

and thus the numerical scheme is charge conservative.

Before proceeding, let us make two remarks. First, the semi-Lagrangian discontinuous Galerkin method is fully explicit. That is, no linear system of equations has to be solved. Second, the order of the space discretization is $\ell +1$ and thus we can easily construct numerical schemes of arbitrary order. Note, however, that the degrees of freedom are directly proportional to the order of the method. We will revisit this issue in some detail in § 4.

Up to this point we have described the numerical schemes used to solve the Vlasov equation. However, what we have not discussed so far is how to solve Poisson’s equation in order to determine the electric field from the charge density. Since the electric field is usually a smooth function, we can use the fast Fourier transformation (FFT). This is straightforward for the cubic spline interpolation as the degrees of freedom already reside on an equidistant grid. However, for the discontinuous Galerkin approach the degrees of freedom are located at the Gauss–Legendre quadrature nodes and thus have to be transferred to an equidistant grid before an FFT algorithm can be applied. In our implementation we simply evaluate the piecewise polynomial in order to obtain an equidistant grid. For example, for a piecewise linear approximation and the cell $[x_{i-1/2},x_{i-1/2}+h]$ we evaluate the approximant at $x_{i-1/2}+(1/3)h$ and $x_{i-1/2}+(2/3)h$ which yields an equidistant grid. In addition, once the electric field has been computed it has to be transferred back from an equidistant grid to the Gauss–Legendre quadrature nodes. In our implementation we use a polynomial interpolation that is of the same order as the discontinuous Galerkin scheme. Since this part of the algorithm takes place in a lower-dimensional space, we could, in principle, use a more elaborate interpolation scheme without significantly increasing the computational cost. In §§ 3 and 4 we will further investigate the effect of this interpolation on the conservation properties of the numerical scheme.

3 Conservation

The purpose of this section is to analyse the conservation properties from a theoretical point of view. The conserved variables (for the continuous system) which we will consider here are: charge, $L^{1}$ norm of the density function, current, energy, entropy and the $L^{2}$ norm of the density function. We will first consider the error that is made due to the time discretization (the splitting approach outlined in the introduction) and then consider the error made due to the discontinuous Galerkin or the cubic spline interpolations that are used to discretize space. An overview of the properties of each scheme is given in table 1.

Table 1.

This table provides an overview of the conserved quantities for the splitting approach in time and the two space discretization strategies. The superscript $\dagger$ means the conservation is only true up to an error made in the computation of the electric field. The superscript $\ast$ emphasizes that positivity preserving semi-Lagrangian discontinuous Galerkin (dG) schemes have been developed that do preserve the $L^{1}$ norm as well as the mass. The degree of the polynomial in each cell (for the semi-Lagrangian discontinuous Galerkin scheme) is denoted by $\ell$ and the corresponding numerical method is consistent of order $\ell +1$ in space.

3.1 Time discretization

3.1.1 Charge and $L^{p}$ norms

Since we can write both parts of the splitting as a translation along a coordinate axis in phase space, we immediately conclude that the splitting conserves the mass as well as all $L^{p}$ norms.

3.1.2 Current

More interesting is the conservation of the (total) current. The current is defined by the following relation

(3.1) $$\begin{eqnarray}\displaystyle j=\int vg(x,v)\,\text{d}(x,v). & & \displaystyle\end{eqnarray}$$

For the first part of the splitting (see (2.11)) we immediately have

(3.2) $$\begin{eqnarray}\int vg(x-vt,v)\,\text{d}(x,v)=\int vg(x,v)\,\text{d}(x,v),\end{eqnarray}$$

which is the desired conservation of the current. For the second part of the splitting (see (2.12)) we obtain

(3.3) $$\begin{eqnarray}\displaystyle \int vg(x,v-E(x)t)\,\text{d}(x,v) & = & \displaystyle \int (v+E(x)t)g(x,v)\,\text{d}(x,v)\nonumber\\ \displaystyle & = & \displaystyle \int vg(x,v)\,\text{d}(x,v)+t\int E(x)\unicode[STIX]{x1D70C}(x)\,\text{d}x.\end{eqnarray}$$

Then from

(3.4) $$\begin{eqnarray}\displaystyle \int E\unicode[STIX]{x1D70C}\,\text{d}x & = & \displaystyle \int E(\unicode[STIX]{x1D735}\boldsymbol{\cdot }E)\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \int E(\unicode[STIX]{x1D735}\boldsymbol{\cdot }E)+\frac{1}{2}\unicode[STIX]{x1D735}E^{2}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle -\int (\unicode[STIX]{x1D735}\times E)\times E\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle 0\end{eqnarray}$$

the conservation of the current follows immediately. Note that we have assumed here that the electric field is computed exactly and that it is curl free. Both of these assumptions are obviously satisfied if we neglect the space discretization. However, we will have to consider these assumptions in more detail in the next section.

3.1.3 Energy

The total energy is given by

(3.5) $$\begin{eqnarray}E=\int \frac{1}{2}v^{2}g(x,v)\,\text{d}(x,v)+\int \frac{1}{2}E^{2}(x)\,\text{d}x,\end{eqnarray}$$

where the first term represents the kinetic energy and the second term represents the energy stored in the electric field. The first part of the splitting conserves the kinetic energy

(3.6) $$\begin{eqnarray}\int v^{2}g(x-tv,v)\,\text{d}(x,v)=\int v^{2}g(x,v)\,\text{d}(x,v)\end{eqnarray}$$

and the second part of the splitting conserves the electric energy (since the charge density is not changed). However, in both cases the reverse statement is not true. Thus, the second part of the splitting does not conserve the electric energy and the first part of the splitting does not conserve the kinetic energy. Most importantly, the complete splitting procedure does not conserve the total energy.

This is to be expected. The splitting scheme considered here is a Hamiltonian splitting (see, for example, Casas et al. Reference Casas, Crouseilles, Faou and Mehrenberger2017). Such numerical methods have been studied extensively in the case of ordinary differential equations (see, for example, Hairer, Lubich & Wanner Reference Hairer, Lubich and Wanner2006). While these numerical methods do not conserve energy, they show excellent long-time properties with respect to energy conservation. The results for ordinary differential equations cannot be directly transferred to the present case. The hope is, however, that similar long-time properties can also be observed for the Vlasov–Poisson equation.

3.1.4 Entropy

The entropy is defined as

(3.7) $$\begin{eqnarray}S=-\int g\log g\,\text{d}(x,v).\end{eqnarray}$$

For physically reasonable initial values (i.e.  $g\geqslant 0$ ) the entropy is well defined and conserved by the continuous system. Entropy is maximized for a density function which has equal probability for each point in phase space. Since diffusive systems naturally tend to this limit, we use entropy as a measure of the amount of diffusion that is introduced by the numerical scheme.

To analyse the splitting scheme, we define $f(t,x,v)=g(x-tv,v)$ . Then, for the first part of the splitting algorithm, we have

(3.8) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\int f\log f\,\text{d}(x,v)=-\int (\log f)\unicode[STIX]{x1D735}\boldsymbol{\cdot }(vf)\,\text{d}(x,v)-\int \unicode[STIX]{x1D735}\boldsymbol{\cdot }(vf)\,\text{d}(x,v).\end{eqnarray}$$

Integration by parts yields

(3.9) $$\begin{eqnarray}\int (\log f)\unicode[STIX]{x1D735}\boldsymbol{\cdot }(vf)\,\text{d}(x,v)=-\int fv\boldsymbol{\cdot }\unicode[STIX]{x1D735}\log f\,\text{d}(x,v)=-\int \unicode[STIX]{x1D735}\boldsymbol{\cdot }(vf)\,\text{d}(x,v)=0\end{eqnarray}$$

and thus entropy is conserved for the first part of the splitting algorithm.

Now, for the second step in the splitting algorithm we define $f(t,x,v)=g(x,v-tE(x))$ and find

(3.10) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\int f\log f\,\text{d}(x,v)=-\int (\log f)\unicode[STIX]{x1D735}_{v}\boldsymbol{\cdot }(Ef)\,\text{d}(x,v)-\int \unicode[STIX]{x1D735}_{v}\boldsymbol{\cdot }(Ef)\,\text{d}(x,v)=0,\end{eqnarray}$$

which shows that entropy is conserved for the splitting scheme.

3.2 Space discretization

3.2.1 Charge and $L^{1}$ norm

The discontinuous Galerkin method is charge conservative (see § 2). The same is true for the spline interpolation. However, both numerical methods can produce negative values. This, in principle, is not an issue for the algorithm considered here. The only difficulty that might arise is the physical interpretation of the resulting particle density function. Note, however, that as long as the magnitude of the negative values are smaller than the accuracy required, this issue is essentially void. However, an interesting aspect is that with respect to conservation of the $L^{1}$ norm we can formulate the following result.

Theorem 1. A numerical algorithm that produces negative values and is mass conservative cannot conserve the $L^{1}$ norm.

Proof. We assume that the initial value $f(0,x,v)$ is non-negative. Then

(3.11) $$\begin{eqnarray}\int |f(0,x,v)|\,\text{d}(x,v)=\int f(0,x,v)\,\text{d}(x,v)=\int f(t,x,v)\,\text{d}(x,v)\leqslant \int |f(t,x,v)|\,\text{d}(x,v).\end{eqnarray}$$

The inequality becomes strict if the particle density function $f(t,x,v)$ is negative. In this case the $L^{1}$ norm is not conserved which is the desired result.

The above result shows that conservation of $L^{1}$ norm can be used as a measure of positivity preservation. Both the discontinuous Galerkin scheme and the method based on cubic splines, at least in their most basic formulation, can yield negative values. Thus, they do not preserve the $L^{1}$ norm. For the discontinuous Galerkin scheme (mass conservative) limiters are available that restrict the density function in such a manner that negative values are avoided (Qiu & Shu Reference Qiu and Shu2011; Rossmanith & Seal Reference Rossmanith and Seal2011). The resulting schemes are then both mass and $L^{1}$ conservative. One should note, however, that while the proposed positivity limiters can be applied locally, associated with them is a certain computational cost. What is even more important is that these procedure might degrade the performance of the scheme with respect to the other conserved quantities. For example, the limiters suggested in Qiu & Shu (Reference Qiu and Shu2011) and Rossmanith & Seal (Reference Rossmanith and Seal2011) operate by shifting and scaling the output from the translation in order to avoid negative values. However, in order to conserve mass, a positive shift must result in a compression of the polynomial function which in turn introduces additional diffusion into the numerical scheme.

On the other hand, for the cubic spline interpolation implementing a positivity limiter is significantly more complicated. Nevertheless, a number of limiters have been suggested in the literature (Zerroukat, Wood & Staniforth Reference Zerroukat, Wood and Staniforth2005, Reference Zerroukat, Wood and Staniforth2006; Crouseilles, Mehrenberger & Sonnendrücker Reference Crouseilles, Mehrenberger and Sonnendrücker2010). It has been well established in this case that adding positivity limiters to the numerical scheme degrades the conservation of the other invariants.

We will not focus on positivity limiters in this paper. However, we will investigate the appearance of negative values in some detail in § 4. In fact, the numerical results obtained seem to suggest that even without such modifications we observe good long-time behaviour for the semi-Lagrangian discontinuous Galerkin scheme.

3.2.2 Current

In the following, we will use $g(x,v)=(Pf(\cdot -vt,v))(x)$ . From § 2 we know that both the discontinuous Galerkin and the cubic spline based schemes are mass conservative and thus

(3.12) $$\begin{eqnarray}\int g(x,v)\,\text{d}x=\int f(x-vt,v)\,\text{d}x.\end{eqnarray}$$

For the first part of the splitting algorithm we thus have

(3.13) $$\begin{eqnarray}\displaystyle \int vg(x,v)\,\text{d}(x,v) & = & \displaystyle \int v\left(\int g(x,v)\,\text{d}x\right)\text{d}v\nonumber\\ \displaystyle & = & \displaystyle \int v\left(\int f(x-vt,v)\,\text{d}x\right)\text{d}v=\int vf(x,v)\text{ d}(x,y),\qquad\end{eqnarray}$$

which implies that the current is conserved for that step of the algorithm.

The situation is more involved for the second part of the splitting algorithm. This is due to the fact that the translation is now in the velocity direction. We will thus consider the two numerical schemes separately, starting with the discontinuous Galerkin method. To simplify the notation let us define $g(x,v)=(Pf(x,\cdot -tE(x)))(v)$ . Expanding $v$ in the Legendre basis (in a given cell) yields

(3.14) $$\begin{eqnarray}v=v_{0}p_{0}+v_{1}p_{1}(v),\end{eqnarray}$$

where $p_{k}$ is the $k$ th Legendre polynomial in the corresponding cell. Thus, we obtain

(3.15) $$\begin{eqnarray}\displaystyle \hspace{-12.0pt}\int vg(x,v)\,\text{d}(x,v) & = & \displaystyle \int \!\left(\int (v_{0}p_{0}+v_{1}p_{1}(v))g(x,v)\,\text{d}v\right)\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \int \!\left[v_{0}p_{0}\!\int f(x,v-tE(x))\,\text{d}v+v_{1}\!\int p_{1}(v)g(x,v)\,\text{d}v\right]\text{d}x.\qquad\end{eqnarray}$$

Now, since we can expand the function $v\mapsto g(x,v)$ into

(3.16) $$\begin{eqnarray}g(x,v)=\mathop{\sum }_{k=0}^{\ell }g_{k}(x)p_{k}(v)\end{eqnarray}$$

we have (for $\ell \geqslant 1$ )

(3.17) $$\begin{eqnarray}\int p_{1}(v)g(x,v)\,\text{d}v=g_{1}(x)=\int p_{1}(v)f(x,v-tE(x))\,\text{d}x\end{eqnarray}$$

and thus

(3.18) $$\begin{eqnarray}\int vg(x,v)\,\text{d}(x,v)=\int vf(x,v-tE(x))\,\text{d}(x,v).\end{eqnarray}$$

The only remaining issue is that we do not use the exact electric field in our computation but replace it by a numerical approximation. In principle, we can use an FFT based approach to compute the electric field up to machine precision. However, the discontinuous Galerkin approach requires an interpolation procedure in order to transfer the data from and to the equidistant grid required to compute the FFT. This transformation is not exact up to machine precision and does introduce an error in the total current. We will investigate this issue in § 4.

Now, let us consider the semi-Lagrangian method based on the cubic spline interpolation. That the first part of the splitting algorithm conserves the current is an immediate consequence of the fact that that scheme is mass conservative. For the second part of the splitting algorithm we have

(3.19) $$\begin{eqnarray}\mathop{\sum }_{j}v_{j}(Pf(x,\cdot -tE(x)))(v_{j})=\mathop{\sum }_{j}v_{j}\,f(x,v_{j}-tE(x))\end{eqnarray}$$

since the projection does not change the value of the spline at the grid points. Thus, we follow that the cubic spline interpolation preserves the current under the assumption that the electric field is determined up to machine precision. Note that satisfying this assumption is easier in the case of the spline based semi-Lagrangian scheme, as the degrees of freedom already reside on an uniform grid and can thus be directly fed into an FFT based Poisson solver.

3.2.3 Energy

We have seen in the previous section that the splitting scheme does not conserve the total energy. Thus, the question posed in this section is if the space discretization introduces an additional error source.

First, let us consider the discontinuous Galerkin scheme. For the first part of the splitting we define $g(x,v)=(Pf(\cdot -vt,x))(x)$ and obtain (using conservation of mass)

(3.20) $$\begin{eqnarray}\int v^{2}g(x,v)\,\text{d}(x,v)=\int v^{2}\left(\int g(x,v)\,\text{d}x\right)\text{d}v=\int v^{2}\left(\int f(x-vt,v)\,\text{d}x\right)\text{d}v.\end{eqnarray}$$

Now, we can rewrite the integrals on the right-hand side as follows

(3.21) $$\begin{eqnarray}\int v^{2}\left(\int f(x-vt,v)\,\text{d}x\right)\text{d}v=\int v^{2}f(x-vt,v)\,\text{d}(x,v),\end{eqnarray}$$

which shows that the kinetic energy is preserved. On the other hand, the first splitting step does not conserve the electric energy. Thus, we consider here only if the space discretization modifies the electric energy that we obtain from the time integrator. This is indeed the case as in general

(3.22) $$\begin{eqnarray}\int f(x-vt,v)\,\text{d}v\neq \int g(x,v)\,\text{d}v,\end{eqnarray}$$

which shows that the space discretization modifies the charge density and thus the electric field. Therefore, the space discretization introduces an additional error source in the electric energy.

For the second part of the splitting we define $g(x,v)=(Pf(x,\cdot -tE(x)))(v)$ and obtain in each cell (by using (3.16))

(3.23) $$\begin{eqnarray}\displaystyle \int v^{2}g(x,v)\,\text{d}(x,v) & = & \displaystyle \int \left(\mathop{\sum }_{k=0}^{2}\mathop{\sum }_{j=0}^{\ell }v_{k}g_{j}(x)\int p_{k}(v)p_{j}(v)\,\text{d}v\right)\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \int \mathop{\sum }_{k=0}^{\min (2,\ell )}v_{k}g_{k}\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \int \left(\mathop{\sum }_{k=0}^{\min (2,\ell )}v_{k}\int p_{k}(v)f(x,v-tE(x))\,\text{d}v\right)\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \int v^{2}f(x,v-tE(x))\,\text{d}(x,v).\end{eqnarray}$$

The last equality only holds for $\ell \geqslant 2$ and shows that in this case, no additional error in the kinetic energy is introduced by the space discretization. Now, the second part of the splitting does not change the charge density and thus the electric energy remains invariant. Since

(3.24) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(x)=\int f(x,v-tE(x))\,\text{d}v=\int g(x,v)\,\text{d}v,\end{eqnarray}$$

the same property holds true for the discontinuous Galerkin approximation and thus the electric energy is conserved during that step (i.e. no additional error is introduced by the space discretization).

Most of the results obtained for the semi-Lagrangian discontinuous Galerkin scheme can be transferred immediately to the case of spline interpolation. Since the projection on the spline subspace does not change the value of the approximant at the grid points, we have

(3.25) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{j}v_{j}^{2}g(x,v_{j})=\mathop{\sum }_{j}v_{j}^{2}f(x,v_{j}-tE(x)) & & \displaystyle\end{eqnarray}$$

and thus the second part of the splitting does not introduce an additional source of error (independent of the order of the underlying polynomials).

3.2.4 Entropy and $L^{2}$ norm

For both the semi-Lagrangian discontinuous Galerkin scheme and the cubic spline based scheme the entropy and the $L^{2}$ norm are not conserved; even though both are conserved by the splitting algorithm. We will evaluate the relative error made in the conservation of these quantities in the next section. However, for now, we will discuss one additional issue not directly related to conservation.

In the Vlasov–Poisson equations no diffusion mechanism is included. Thus, the $L^{2}$ norm and the entropy are conserved for the continuous model. However, once diffusion is added, the $L^{2}$ norm will decrease in time. This is due to the fact that diffusion decreases the variance of the density function $f$ and we can write $\Vert f\Vert _{2}^{2}=1+\text{Var}(f)$ . In addition, diffusion increases the entropy (which is a measure of the disorder of a system). The latter property is, of course, just a statement of the second law of thermodynamics. Thus, the decrease in the $L^{2}$ norm and the increase in entropy can be considered as a measure of the amount of diffusion in a given physical system. In our case they are used to quantify the amount of diffusion that a numerical scheme (wrongly) introduces. This point is interesting because in the numerical simulations conducted in the next section we observe that the semi-Lagrangian discontinuous Galerkin scheme always decreases/increases the $L^{2}$ norm/entropy while the spline based approach shows an oscillating behaviour for some problems. Thus, the spline based semi-Lagrangian method cannot be considered as the exact solution of a perturbed physical system (as no such system would be allowed to decrease the entropy; a consequence of the second law of thermodynamics).

Analytically, we can show that the $L^{2}$ norm decreases in time for the semi-Lagrangian discontinuous Galerkin scheme.

Theorem 2. The $L^{2}$ norm is a decreasing function of time for the semi-Lagrangian discontinuous Galerkin scheme.

Proof. Due to the structure of the Vlasov–Poisson equation it is sufficient to consider a one-dimensional advection problem. Thus, we consider a translated function $g(x)$ which is piecewise polynomial and thus lies in $L^{2}$ . We can then find a Legendre series

(3.26) $$\begin{eqnarray}\mathop{\sum }_{k=0}^{\infty }g_{k}p_{k}(x),\end{eqnarray}$$

which converges in the $L^{2}$ norm to $g$ (Pollard Reference Pollard1972). Now, we have

(3.27) $$\begin{eqnarray}\Vert Pg\Vert _{2}^{2}-\Vert g\Vert _{2}^{2}=\mathop{\sum }_{k=0}^{\ell }g_{k}^{2}-\mathop{\sum }_{k=0}^{\infty }g_{k}^{2}\leqslant 0,\end{eqnarray}$$

which implies that the $L^{2}$ norm decreases in each time step.

Based on an extensive numerical investigation, we conjecture that a similar result holds true for the entropy. In fact, for piecewise constant polynomials this immediately follows from the convexity of the logarithm and for piecewise linear polynomials a similar (although more tedious) argument can be given. However, we were not able to prove the corresponding result for piecewise polynomials of arbitrary degree.

There is one additional complication in the case of the entropy. As the numerical scheme can generate negative values, the entropy as given by (3.7) is not even defined. Since we interpret negative values as zero particle density, we use

(3.28) $$\begin{eqnarray}S=-\int \max (g,0)\log g\,\text{d}(x,v)\end{eqnarray}$$

instead. Of course, using this formula can still result in a decrease of entropy. However, we will see in the next section that this issue only plays a negligible role in the numerical simulations.

4 Numerical results

Before we discuss the results of the medium and long-time simulations, which are the main focus of this paper, let us start by establishing some baseline results for relatively small final times in case of the linear Landau damping (see § 4.4 for more details on the parameters used). In this case we have computed the error for both the electric field $E$ and the particle distribution $f$ at final times $T=12.5$ and $T=100$ . The corresponding results are shown in figure 3.

Figure 3.

The errors for linear Landau damping in the electric field $E$ and the particle distribution $f$ are shown as functions of the degrees of freedom (per dimension) for the second-, fourth- and sixth-order discontinuous Galerkin method and the cubic spline interpolation. The results for final time $T=12.5$ are shown in (a,b) and the results for final time $T=100$ are shown in (c,d).

For final time $T=12.5$ we observe convergence of the discontinuous Galerkin schemes with the expected order. In addition, the sixth-order scheme is always superior to the fourth-order scheme which in turn is superior to the second-order scheme. The comparison with the cubic spline interpolation is a bit more nuanced. In particular, we see a clear advantage of the fourth- and sixth-order discontinuous Galerkin methods for the error in the electric field. On the other hand, with respect to the error in the distribution function, the cubic spline interpolation performs best (except for very stringent accuracy requirements).

For the final time $T=100$ the analytic solution is almost identical to zero. However, numerical methods have significant difficulties due to the recurrence effect (see § 4.4 for more details). We can see this in the numerical results very clearly as initially no convergence is observed (up to $128$ degrees of freedom). Since the recurrence effect depends mostly on the grid spacing, the spline interpolation performs best in this case (this is due to the fact that an equidistant grid minimizes the grid spacing). We also note that for low tolerances the second-order discontinuous Galerkin (dG) method performs better than the higher-order variants if one is only interested in resolving the error in the electric field. Some of this conclusions have been pointed out in our earlier work. For example, it was observed in Crouseilles et al. (Reference Crouseilles, Mehrenberger and Vecil2011b ), Einkemmer & Ostermann (Reference Einkemmer and Ostermann2014b ) that higher-order discontinuous Galerkin methods only have a weak influence on the recurrence effect.

The remainder of this section is divided into two parts. In the first part we mainly discuss the properties of the semi-Lagrangian discontinuous Galerkin scheme and compare its efficiency to the cubic spline interpolation. We will consider both classic test problems (linear and nonlinear Landau damping, bump-on-tail instability) and somewhat less well-known test problems (expansion of a plasma blob, plasma echo). In the second part we will consider the long-time behaviour of the semi-Lagrangian discontinuous Galerkin scheme. These simulations are conducted with the highly parallel implementation described in Einkemmer (Reference Einkemmer2016b ).

Since one of our goals is to compare the discontinuous Galerkin scheme with the cubic spline interpolation, a measure of the cost of a given numerical scheme has to be established. In a sequential implementation both computational efficiency as well as memory consumption are possible choices. However, comparing arithmetic operations on a modern (vectorized and multi-core) central processing unit (CPU) requires a very delicate analysis that depends on both algorithmic and implementation details (for example, taking into account how the linear system is solved in the cubic spline interpolation or whether such a scheme could be vectorized). In addition, on a distributed memory system (i.e. on a cluster), especially for the cubic spline interpolation, further problems occur due to the all-to-all communication required in solving the linear system. Some remedies have been proposed in Crouseilles et al. (Reference Crouseilles, Latu and Sonnendrücker2009a ). However, these usually change the numerical scheme in a non-trivial way (which might also have a detrimental effect on the conserved quantities). In fact, this problem is one of the major motivations why we consider the discontinuous Galerkin approach.

In the present study we will exclusively use the number of degrees of freedom, i.e. the memory consumed, as the measure of efficiency. This is reasonable because for higher-dimensional simulations the amount of memory consumed is often the most important criterion for Vlasov simulations. But perhaps even more important is the fact that such algorithms are bandwidth limited on all modern architectures. Consequently, how many bytes have to be transferred to and from the CPU determines the run time of the algorithm. We have implemented our discontinuous Galerkin approach on a dual socket Intel Xeon E5-2630 v3 workstation with a (theoretical) memory bandwidth of $59~\text{GB}~\text{s}^{-1}$ . For the second-order implementation (dG2)/fourth-order (dG4) implementation we achieve approximately 94 %/89 % of the theoretically possible performance on that machine. That means that for the advection step with contiguous memory access we require approximately 0.3 ns per degree of freedom. Let us also note that this means, in particular, that no spline implementation is able to outperform our algorithm by more than 10 % (in all likelihood, even an optimized spline implementation will probably be significantly slower due to the fact that extra data structures are required to store the spline coefficients). We further note that the discontinuous Galerkin method has been implemented on graphic processing units and the Intel Xeon Phi. The corresponding performance measurements, which show a significant speed-up on both of these architectures, are given in Einkemmer (Reference Einkemmer2016a ).

In the remainder of this section we will conduct two-dimensional simulations of the Vlasov–Poisson equations using 128 degrees of freedom in each dimension (the coarse problem) and 512 degrees of freedom in each dimension (the fine problem).

4.1 Nonlinear Landau damping

We begin with the classic problem of nonlinear Landau damping. The initial value

(4.1) $$\begin{eqnarray}f(0,x,v)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\text{e}^{-v^{2}/2}(1+\unicode[STIX]{x1D6FC}\cos kx)\end{eqnarray}$$

is posed on the domain $[0,4\unicode[STIX]{x03C0}]\times [-6,6]$ with $k=1/2$ and $\unicode[STIX]{x1D6FC}=1/2$ . Periodic boundary conditions are imposed in both the $x$ - and the $v$ -direction. A large body of physics literature has been devoted to this problem and a number of numerical simulations have been conducted (see, for example, Manfredi (Reference Manfredi1997) and Brunetti, Califano & Pegoraro (Reference Brunetti, Califano and Pegoraro2000)). As such, the behaviour, at least for moderate times, is fairly well understood. Initially, we observe a decay in the electric energy (the Landau damping phenomenon) that after some time stabilizes and gives way to oscillations of the electric energy. The nonlinear Landau damping is numerically challenging because even though the initial value is perfectly smooth, small-scale structures (so-called filaments) develop in the two-dimensional phase space as the system evolves in time. Thus, in order to integrate the system up to a given precision an extremely high resolution in velocity space is needed. However, what we are interested in is to what extent a relatively coarse discretization can give qualitatively correct results.

The numerical simulations for the cubic spline and the discontinuous Galerkin approach are given in figure 4 (for the coarse discretization) and in figure 5 (for the fine discretization). We observe that for the discontinuous Galerkin method the error in the current is close to machine precision (even though this quantity is not exactly conserved). The cubic spline interpolation has an advantage with respect to the energy error (which, however, is quite accurately conserved for this test problem). We also observe a distinct difference between the second-order method ( $p=1$ ) and the higher-order methods with respect to energy conservation (as predicted by the results derived in § 3).

Figure 4.

This figure shows the particle density $f$ at $T=400$ and the time evolution of the electric energy for the nonlinear Landau damping problem. In addition, the error in the current, energy, entropy, $L^{1}$ norm and $L^{2}$ norm are shown. For all numerical schemes $128$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

Figure 5.

This figure shows the particle density $f$ at $T=400$ and the time evolution of the electric energy for the nonlinear Landau damping problem. In addition, the error in the current, energy, entropy $L^{1}$ and $L^{2}$ norm are shown. For all numerical schemes $512$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

On the other hand, the discontinuous Galerkin scheme has an advantage with respect to the $L^{1}$ norm (an error of 5 % versus 1–2 % for the coarse problem). Thus, the discontinuous Galerkin scheme produces fewer negative values. If we consider the qualitative features of the numerical solution we might be tempted to conclude that the spline interpolation reproduces more clearly the features of the exact solution. However, upon closer inspection of the results for the fine discretization (where the maximal amplitude is significantly reduced for the spline interpolation) and taking into account the error in the $L^{1}$ norm, we conclude that in fact the spline interpolation produces significant overshoot into the numerical simulation.

Let us also discuss the effect of the order of the discontinuous Galerkin approximation. For the coarse problem the order only plays a minor role (although energy conservation is better once we consider a numerical scheme with piecewise polynomials of degree $2$ or above). Other than that the performance of the numerical scheme is almost independent of the order.

4.2 Bump-on-tail instability

The second problem we consider is the so-called bump-on-tail instability. In this case we prescribe the initial value

(4.2) $$\begin{eqnarray}f(0,x,v)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}(\unicode[STIX]{x1D6FC}\text{e}^{-v^{2}/2}+\unicode[STIX]{x1D6FD}\text{e}^{-4(v-2.5)^{2}}(1+\unicode[STIX]{x1D6FE}\cos x))\end{eqnarray}$$

on the domain $[0,4\unicode[STIX]{x03C0}]\times [-6,6]$ and with $\unicode[STIX]{x1D6FC}=0.8$ , $\unicode[STIX]{x1D6FD}=0.2$ and $\unicode[STIX]{x1D6FE}=0.1$ . Periodic boundary conditions in both the $x$ - and the $v$ -direction are imposed. Thus, we have a stationary thermal plasma in which 80 % of the mass is concentrated. In addition, a smaller beam with average velocity $v=2.5$ penetrates the stationary plasma. This configuration is an equilibrium of the Vlasov–Poisson equation. However, due to the perturbation that is added, a travelling vortex appears in phase space. The bump-on-tail instability shows a quiet phase that is succeeded by the relative rapid development of the vortex which then remains stable for a long time. Since the bump-on-tail instability exhibits chaotic behaviour the onset of the vortex might be different for different numerical schemes (even for relatively fine discretization).

The numerical simulations for the cubic spline and the discontinuous Galerkin approach are given in figure 6 (for the coarse discretization) and in figure 7 (for the fine discretization). The current in case of the discontinuous Galerkin scheme is only conserved up to $10^{-5}$ which clearly confirms the behaviour discussed in § 3. However, this fact does not seem to diminish the overall performance of the numerical scheme. For the higher-order discontinuous Galerkin scheme and for the cubic spline interpolation all the remaining conserved quantities show a similar behaviour. The exception being the second-order discontinuous Galerkin scheme which exhibits significantly worse performance across the board (contrary to the results obtained in the previous section for the nonlinear Landau damping).

Figure 6.

This figure shows the particle density $f$ at $T=400$ and the time evolution of the electric energy for the bump-on-tail instability. In addition, the error in the current, energy, entropy, $L^{1}$ norm and $L^{2}$ norm are shown. For all numerical schemes $128$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

Figure 7.

This figure shows the particle density $f$ at $T=400$ and the time evolution of the electric energy for the bump-on-tail instability. In addition, the error in the current, energy, entropy $L^{1}$ and $L^{2}$ norm are shown. For all numerical schemes $512$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

4.3 Expansion into a uniform ion background

In this section we consider the following initial value

(4.3) $$\begin{eqnarray}f(0,x,v)=\frac{1}{2\unicode[STIX]{x03C0}}\text{e}^{-v^{2}/2}\text{e}^{-(x-2\unicode[STIX]{x03C0})^{2}/2}\end{eqnarray}$$

on the domain $[0,4\unicode[STIX]{x03C0}]\times [-6,6]$ . As before, we impose periodic boundary conditions in both directions. This is the only test problem where the initial value is localized in space as well as in velocity. This problem can be interpreted as the expansion of an electron blob (in thermodynamic equilibrium) into a region of uniform ion density. This is a challenging problem because there is no physical mechanism that would hold the initial blob in place. Thus, each particle with a given velocity expands at its own speed which due to the periodicity of the problem yields a very irregular phase space distribution (similar to the nonlinear Landau damping). What sets this problem apart from the nonlinear Landau damping, however, is that the filamentation is driven by the free streaming part of the Vlasov equation and is not merely superimposed on top of a Maxwellian distribution in velocity space.

The numerical simulations for the cubic spline and the discontinuous Galerkin approach are given in figure 8 (for the coarse discretization) and in figure 9 (for the fine discretization). This problem seems to be especially challenging for the cubic spline interpolation. For the coarse discretization we observe an error of almost 50 % in the $L^{1}$ norm, which renders any interpretation of the particle density function impossible. In addition, we observe two unphysical peaks in the electric energy.

Figure 8.

This figure shows the particle density $f$ at $T=400$ and the time evolution of the electric energy for the expansion into a uniform ion background. In addition, the error in the current, energy, entropy, $L^{1}$ norm and $L^{2}$ norm are shown. For all numerical schemes $128$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

Figure 9.

This figure shows the particle density $f$ at $T=400$ and the time evolution of the electric energy for the expansion into a uniform ion background. In addition, the error in the current, energy, entropy, $L^{1}$ norm and $L^{2}$ norm are shown. For all numerical schemes $512$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

On the other hand, none of these difficulties appear for the semi-Lagrangian discontinuous Galerkin scheme. With respect to the $L^{1}$ norm the second-order dG method is most favourable although that method only conserves energy up to an accuracy of $10^{-3}$ (compared to $10^{-5}$ for the higher-order methods). This favourable behaviour is despite the fact that entropy conservation is somewhat better for the cubic spline approximation and the $L^{2}$ norm conservation is significantly better.

Let us continue to discuss the cubic spline interpolation. For this method the entropy and the $L^{2}$ norm show an oscillating behaviour in time. Clearly this is a numerical artefact as has been discussed in § 3. None of these artefacts are present for the discontinuous Galerkin approach (as is expected based on the theoretical considerations in § 3).

4.4 Linear Landau damping

The linear Landau damping problem is just the strong Landau damping introduced earlier but with $\unicode[STIX]{x1D6FC}$ chosen sufficiently small such that nonlinear effects are small. In this section we consider $\unicode[STIX]{x1D6FC}=10^{-2}$ . We have postponed this classic problem until now because in some sense the problem is very easy while in another sense it is very difficult. The numerical results for the coarse grid are shown in figure 10. All of the conserved quantities look very reasonable. In fact, $L^{1}$ is conserved up to machine precision and the $L^{2}$ norm and the entropy are conserved up to approximately an error of $10^{-5}$ (by far the lowest error among the problems considered so far). Nevertheless, the result of the numerical simulation is completely wrong (even considering a qualitative analysis). From analytical results it is well known that the electric energy decays exponentially as a function in time. What we observe for the numerical simulation, however, is that it oscillates. This is the so-called recurrence effect; a numerical artefact which implies that the number of grid points has to be chosen proportional to the final time in order to obtain reasonable results. This was already recognized in Cheng & Knorr (Reference Cheng and Knorr1976) and is the reason why despite the favourable conservation properties the linear Landau damping problem is difficult for any numerical scheme.

Figure 10.

This figure shows the time evolution of the electric energy for the linear Landau damping. In addition, the error in the current, energy, entropy, $L^{1}$ norm and $L^{2}$ norm are shown. For all numerical schemes $128$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

Unfortunately, a number of interesting problems in plasma physics (such as plasma echoes; see Galeotti & Califano (Reference Galeotti and Califano2005) or Hou, Ma & Yu (Reference Hou, Ma and Yu2011) for a good exposition) suffer from the recurrence effect. The problem can be overcome to some extent by filtering high frequencies in phase space (a technique referred to as filamentation filtration). In fact, a significant body of literature has been developed on this topic (Klimas & Farrell Reference Klimas and Farrell1994; Eliasson Reference Eliasson2002; Einkemmer & Ostermann Reference Einkemmer and Ostermann2014a ). However, even though such methods can alleviate the problem to some extent, in many situations a fine velocity discretization is still required in order to resolve the interesting dynamics of the system. However, for such a fine discretization the conservation behaviour is even better than for the relative coarse discretization used here.

At this point one might rightfully object that there is no guarantee that filamentation filtration will respect the invariants (even if they are conserved by the basic numerical scheme). Thus it seems prudent to investigate the effect of such a modification on the algorithm outlined in § 2. The results for the method introduced in Einkemmer & Ostermann (Reference Einkemmer and Ostermann2014a ) are shown in figure 10. There is some effect on the entropy and the $L^{1}$ norm which, however, are still conserved up to an accuracy of $10^{-2}$ . However, despite the fact that we essentially remove high frequencies neither the mass nor any of the conserved quantities, except for the $L^{1}$ norm and the entropy, are adversely affected. This behaviour is due to the fact that the magnitude of the high frequencies is very small and thus negative values appear. Since we cutoff a couple of frequencies this averages out and leaves the mass intact (while affecting the $L^{1}$ norm).

Before proceeding, let us return to the question of how coarse a space discretization is sufficient in order to obtain reasonable results. To do this we consider a model of the plasma echo phenomenon. That is, we prescribe the initial value

(4.4) $$\begin{eqnarray}f(0,x,v)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\text{e}^{-v^{2}/2}(1+\unicode[STIX]{x1D6FC}\cos k_{1}x)\end{eqnarray}$$

with $\unicode[STIX]{x1D6FC}=10^{-3}$ and $k_{1}=12\unicode[STIX]{x03C0}/100$ . This initial perturbation is damped away by Landau damping. Then at time $t=200$ we excite a second perturbation by adding

(4.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FC}}{\sqrt{2\unicode[STIX]{x03C0}}}\text{e}^{-v^{2}/2}\cos k_{2}x,\end{eqnarray}$$

where $k_{2}=25\unicode[STIX]{x03C0}/100$ , to the particle density function. This perturbation is also damped away by Landau damping. However, all the information is still stored in the particle density function and due to constructive interference an echo appears at $t=400$ . In fact, multiple echoes can be observed. We will investigate how well the first echo is resolved by the semi-Lagrangian discontinuous Galerkin method. The numerical results are shown in figure 11. We note that in this situation, increasing the order in the velocity direction while keeping the number of degrees of freedom constant significantly degrades the numerical solution. This is particular prominent for approximations of order three and higher, where the plasma echo is barely recognizable. In the spatial direction the results are comparable no matter what order is chosen. Note, however, that decreasing the number of degrees of freedom rapidly diminishes the numerical solution and the wave echo disappears. Thus, there is no point in studying the conservation properties at low resolution in the spatial direction as a certain number of cells are required even if one is only interested in resolving the first echo accurately.

4.6 Long-time behaviour (bump-on-tail instability)

In this section we investigate the long-time behaviour of the semi-Lagrangian discontinuous Galerkin approach for the bump-on-tail instability by integrating the Vlasov–Poisson equation until $T=2\times 10^{4}$ . The numerical results are shown in figure 13. Also in this case the second-order method produces qualitatively wrong results. However, the conservation properties are also markedly worse for the second-order method as compared to the fourth- and sixth-order methods. This is especially true for the entropy and $L^{2}$ norm (which measures diffusion) and we are thus not surprised that the second-order method washes out the vortex in phase space. We also observe that the error in the current is quite large (conserved only up to approximately $10^{-2}$ ). This is true independent of the order of the method. However, this behaviour is not detrimental to the other conserved quantities nor is it detrimental to the qualitative properties of the numerical solution.

Figure 11.

This figure shows the evolution of the electric energy for the plasma echo problem. The order of the dG method is indicated (dGab denotes a method with order a in the $x$ -direction and order b in the $v$ -direction) and the number of cells are given in parenthesis.

Figure 12.

This figure shows the particle density $f$ at $T=2\times 10^{4}$ and the time evolution of the electric energy for the nonlinear Landau damping. In addition, the error in the current, energy, entropy, $L^{1}$ norm and $L^{2}$ norm are shown. For all numerical schemes $128$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

4.7 Vortex merging in the bump-on-tail instability

In the previous section only one vortex appears in the bump-on-tail instability. However, for a certain regime of parameters multiple vortices can be observed in phase space. The difficulty for a numerical scheme is then to preserve these vortices as long as possible. However, most numerical schemes eventually lead to a secondary instability that results in the three vortices merging into a single one. In the following we will consider the problem from Crouseilles et al. (Reference Crouseilles, Mehrenberger and Vecil2011b ) for which three vortices are observed in phase space. The numerical results are shown in figure 14. We note that vortex merging starts at approximately $t=650$ for the fourth-order discontinuous Galerkin method and at approximately $t=1200$ for the spline interpolation. The sixth-order discontinuous Galerkin schemes performs best in this test (vortex merging starts at approximately $t=1250$ ).

Figure 13.

This figure shows the particle density $f$ at $T=2\times 10^{4}$ and the time evolution of the electric energy for the bump-on-tail instability. In addition, the error in the current, energy, entropy, $L^{1}$ norm and $L^{2}$ norm are shown. For all numerical schemes $128$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.

Figure 14.

This figure shows the particle density $f$ at different times for the three vortex bump-on-tail instability. For all numerical schemes $128$ degrees of freedom are employed per space dimension. The order of the dG method is indicated and the number of cells are given in parenthesis.