A low-rank projector-splitting integrator for the Vlasov--Maxwell equations with divergence correction

TL;DR

This paper introduces a low-rank integrator for the Vlasov--Maxwell equations that reduces computational cost by approximating the solution within a low-rank manifold, effectively capturing plasma dynamics and enforcing divergence constraints.

Contribution

It develops a dynamic low-rank integrator for the Vlasov--Maxwell system with divergence correction, enabling efficient and accurate plasma simulations.

Findings

The scheme accurately captures key plasma features in numerical tests.

It maintains divergence constraints with machine precision.

The low-rank approach reduces computational complexity.

Abstract

The Vlasov--Maxwell equations are used for the kinetic description of magnetized plasmas. As they are posed in an up to 3+3 dimensional phase space, solving this problem is extremely expensive from a computational point of view. In this paper, we exploit the low-rank structure in the solution of the Vlasov equation. More specifically, we consider the Vlasov--Maxwell system and propose a dynamic low-rank integrator. The key idea is to approximate the dynamics of the system by constraining it to a low-rank manifold. This is accomplished by a projection onto the tangent space. There, the dynamics is represented by the low-rank factors, which are determined by solving lower-dimensional partial differential equations. The proposed scheme performs well in numerical experiments and succeeds in capturing the main features of the plasma dynamics. We demonstrate this good behavior for a range of test problems. The coupling of the Vlasov equation with the Maxwell system, however, introduces additional challenges. In particular, the divergence of the electric field resulting from Maxwell's equations is not consistent with the charge density computed from the Vlasov equation. We propose a correction based on Lagrange multipliers which enforces Gauss' law up to machine precision.

Figures

Equations

$${\partial }_{t}f(t,x,v)+v\cdot {\nabla }_{x}f(t,x,v)-F(f)\cdot {\nabla }_{v}f(t,x,v)=0.$$

$${\partial }_{t}f(t,x,v)+v\cdot {\nabla }_{x}f(t,x,v)-\left(E(t,x)+v\times B(t,x)\right)\cdot {\nabla }_{v}f(t,x,v)=0$$

$${\partial }_{t}E(t,x)=c^2\nabla \times B(t,x)-j(t,x)$$

$$\nabla \cdot E(t,x)=n_0-\rho (t,x)$$

$${\partial }_{t}B(t,x)=-\nabla \times E(t,x)$$

$$\nabla \cdot B(t,x)=0.$$

$$\rho (t,x)={\int }_{{\Omega }_v}f(t,x,v)\mathrm{d}v\text{and}j(t,x)={\int }_{{\Omega }_v}vf(t,x,v)\mathrm{d}v,$$

$${\partial }_t\rho +\nabla \cdot j=0.$$

$$f(t,x,v)=\sum _{i,j=1}^r{X}_i(t,x)S_{ij}(t)V_j(t,v).$$

$${\partial }_t{f}^{(1)}(t,x,v)+v\cdot {\nabla }_x{f}^{(1)}(t,x,v)-(E^{(1)}(t,x)+v\times B^{(1)}(t,x))\cdot {\nabla }_v{f}^{(0)}(v)=0$$

$${\partial }_t{E}^{(1)}(t,x)=\nabla \times B^{(1)}(t,x)-j^{(0)}-{\int }_{{\Omega }_v}v{f}^{(1)}(t,x,v)\mathrm{d}v$$

$${\partial }_t{B}^{(1)}(t,x)=-\nabla \times E^{(1)}(t,x),$$

$${\partial }_t{\hat{f}}_k^{(1)}(t,v)+i(v\cdot k){\hat{f}}_k^{(1)}(t,v)-\left({\hat{E}}_k^{(1)}(t)+v\times {\hat{B}}_k^{(1)}(t)\right)\cdot {\nabla }_v{\hat{f}}^{(0)}(v)=0$$

$${\partial }_t{\hat{E}}_k^{(1)}(t)=ik\times {\hat{B}}_k^{(1)}(t)-j^{(0)}{\delta }_{k,0}-{\int }_{{\Omega }_v}v{\hat{f}}_k^{(1)}(t,v)\mathrm{d}v$$

$${\partial }_t{\hat{B}}_k^{(1)}(t)=-ik\times {\hat{E}}_k^{(1)}(t).$$

$$f^{(1)}(0,x,v)=\sum _{i=1}^m{\hat{f}}_{k_i}^{(1)}(0,v)e^{i{k}_{i}x},B^{(1)}(0,x)=\sum _{i=1}^m{\hat{B}}_{k_i}^{(1)}(0)e^{i{k}_{i}x}.$$

$${\hat{E}}_k^{(1)}(0)=\frac{ik}{k^2}{\int }_{{\Omega }_v}{\hat{f}}_k^{(1)}(0,v)\mathrm{d}v=0.$$

$${\hat{f}}_k^{(1)}(t,v)=0,{\hat{E}}_k^{(1)}(t)=0,{\hat{B}}_k^{(1)}(t)=0.$$

$$f(t,x,v)=f^{(0)}(v)+\sum _{i=0}^m{\hat{f}}_{k_i}^{(1)}(t,v)e^{i{k}_{i}x}$$

$$f(t,x,v)=\sum _{i,j=1}^r{X}_i(t,x)S_{ij}(t)V_j(t,v),$$

$$\mathcal{M}=\{$$

$$X_i(x)\in L^2({\Omega }_x),V_j(v)\in L^2({\Omega }_v),⟨X_i,X_k{⟩}_x={\delta }_{ik},⟨V_j,V_l{⟩}_v={\delta }_{jl}\},$$

$${\mathcal{T}}_f\mathcal{M}=\{\dot{f}$$

$$\text{with}\dot{S}=({\dot{S}}_{ij})\in {\mathbb{R}}^{r\times r},{\dot{X}}_i(x)\in L^2({\Omega }_x),{\dot{V}}_j(v)\in L^2({\Omega }_v),⟨X_i,{\dot{X}}_k{⟩}_x=0,⟨V_i,{\dot{V}}_j{⟩}_v=0\}.$$

$$P(f)g=\sum _j⟨V_j,g{⟩}_v{V}_j-\sum _{i,j}{X}_i⟨X_i{V}_j,g{⟩}_{x,v}{V}_j+\sum _i{X}_i⟨X_i,g{⟩}_x.$$

$$P(f)g=P_{\bar{V}}g-P_{\bar{V}}{P}_{\bar{X}}g+P_{\bar{X}}g.$$

$${\partial }_{t}f=-P(f)\left(v\cdot {\nabla }_{x}f-(E+v\times B)\cdot {\nabla }_{v}f\right).$$

$${\partial }_{t}f=-P_{\bar{V}}\left(v\cdot {\nabla }_{x}f-(E+v\times B)\cdot {\nabla }_{v}f\right)$$

$${\partial }_{t}f=P_{\bar{V}}{P}_{\bar{X}}\left(v\cdot {\nabla }_{x}f-(E+v\times B)\cdot {\nabla }_{v}f\right)$$

$${\partial }_{t}f=-P_{\bar{X}}\left(v\cdot {\nabla }_{x}f-(E+v\times B)\cdot {\nabla }_{v}f\right).$$

$$f(t,x,v)=\sum _j{K}_j(t,x)V_j(t,v)\text{with}{K}_j(t,x)=\sum _i{X}_i(t,x)S_{ij}(t).$$

$$\sum _j{\partial }_t{K}_j(t,x)V_j(t,v)+\sum _j{K}_j(t,x){\partial }_t{V}_j(t,v)=$$

$$-\sum _j⟨V_j(t,\cdot ),\left[v\mapsto v\cdot {\nabla }_{x}f(t,x,v)-\left(E(x)+v\times B(x)\right)\cdot {\nabla }_{v}f(t,x,v)\right]{⟩}_v{V}_j(t,v).$$

$${\partial }_t{K}_j(t,x)=-\sum _l{c}_{jl}^1\cdot {\nabla }_x{K}_l(t,x)+\sum _l{c}_{jl}^2\cdot E(t,x)K_l(t,x)+\sum _l{c}_{jl}^3\cdot B(t,x)K_l(t,x).$$

$$c_{jl}^1={\int }_{{\Omega }_v}v{V}_j^0(v)V_l^0(v)\mathrm{d}v,c_{jl}^2={\int }_{{\Omega }_v}{V}_j^0(v){\nabla }_v{V}_l^0(v)\mathrm{d}v.$$

$$c_{jl}^3={\int }_{{\Omega }_v}{V}_j^0(v)\left({\nabla }_v{V}_l^0(v)\times v\right)\mathrm{d}v.$$

$$K_j^1(x)=\sum _i{X}_i^1(x)S_{ij}^1,$$

Full text

A low-rank projector-splitting integrator for the Vlasov–Maxwell equations with divergence correction

Lukas Einkemmer

Alexander Ostermann

Chiara Piazzola

Abstract

The Vlasov–Maxwell equations are used for the kinetic description of magnetized plasmas. As they are posed in an up to 3+3 dimensional phase space, solving this problem is extremely expensive from a computational point of view. In this paper, we exploit the low-rank structure in the solution of the Vlasov equation. More specifically, we consider the Vlasov–Maxwell system and propose a dynamic low-rank integrator. The key idea is to approximate the dynamics of the system by constraining it to a low-rank manifold. This is accomplished by a projection onto the tangent space. There, the dynamics is represented by the low-rank factors, which are determined by solving lower-dimensional partial differential equations. The proposed scheme performs well in numerical experiments and succeeds in capturing the main features of the plasma dynamics. We demonstrate this good behavior for a range of test problems. The coupling of the Vlasov equation with the Maxwell system, however, introduces additional challenges. In particular, the divergence of the electric field resulting from Maxwell’s equations is not consistent with the charge density computed from the Vlasov equation. We propose a correction based on Lagrange multipliers which enforces Gauss’ law up to machine precision.

1 Introduction

A kinetic description of plasmas is required in many applications; for example, in magnetic confined fusion. In this paper, we are interested in the Vlasov equation

[TABLE]

The Vlasov equation describes a collisionless plasma, where $f$ is the density of particles in $\mathbb{R}^{d_{x}+d_{v}}$, with $d_{x},d_{v}=1,2,3$. The physically interesting case is $d_{v}\geq d_{x}$. The particles are subjected to a force $F$, which gives rise to a nonlinear term. When considering a plasma interacting with an electromagnetic field the Lorentz force has to be considered. It is given by $F=E+v\times B$, where $E$ and $B$ denote the electric and magnetic fields, respectively. Moreover, the Vlasov equation has to be coupled with Maxwell’s equations, which describe the time evolution of the electric and magnetic fields subject to two divergence constraints; Gauss’ laws for electricity and magnetism, respectively.

The numerical solution of the Vlasov–Maxwell equations faces several challenges. Those difficulties are primarily caused by the high dimensionality of the system. A full discretization of phase space with $n_{x}$ grid points per direction in $x$ and $n_{v}$ points per direction in $v$ leads to a storage cost of $\mathcal{O}(n_{x}^{d_{x}}n_{v}^{d_{v}})$, which is extremely expensive. This is often referred to as the curse of dimensionality. To alleviate this problem, particle in cell (PIC) methods have been widely used. In this approach, only the physical space is discretized, and a conglomerate of particles is advanced by following the characteristic curves of the Vlasov equation; see, e.g., the review paper [36]. Numerical integrators based on fully discretizing phase space, i.e. the Eulerian approach, have become affordable only in recent years, both due to improved algorithms and the greatly increased performance of computer systems (see, e.g., [5, 33, 13, 12, 17, 39, 35, 37]). Algorithms based on splitting methods have been considered. The idea in this setting is to split the problem into a sequence of simple advection equations and then to solve those by semi-Lagrangian or spectral methods. We refer, e.g., to [2, 19, 3, 4, 34, 29].

Solvers for the Vlasov–Poisson equations based on low-rank approximations were proposed in the last few years in order to reduce computation and storage cost [10, 22]. In both papers, the authors discretize first, perform a time step, and then use tensor truncation to project back to the low-rank approximation space. Recently, a different approach was proposed in the context of the Vlasov–Poisson equations; see [15] and the follow-up papers [11, 14]. The key idea is to approximate the dynamics of the system by constraining it to the low-rank manifold. The solution is then represented in low-rank form already at the continuous level, before any space or time discretization is performed. Differential equations for the low-rank factors of the solutions are derived, such that the dynamics of the systems is shifted to the low-rank factors representing the solution. These partial differential equations can then be discretized by any suitable numerical method. A convergence proof for a simplified model, the radiative transfer equation, has been conducted in [9].

An integrator based on this approach was first considered in [25] in the context of matrix differential equations. It is known as the projector-splitting integrator for the dynamic low-rank approximation. More theoretical insights were given in [20, 18] and extensions for various algorithms and tensor formats can be found in [21, 26, 27, 28]. Moreover, the dynamic low-rank approximation was employed for different types of problems, originally in quantum mechanics; see [31, 30, 23, 24].

In the context of kinetic equations, the low-rank integrator has the drawback that it does not preserve the physical structure of the system. That is, physical invariants are in general not preserved leading to severe limitations on the long time integration; see [15]. In [14] the authors studied this issue and proposed a quasi-conservative scheme for the preservation of mass and momentum for the Vlasov–Poisson equations.

In this paper, we propose a dynamic low-rank integrator for the Vlasov–Maxwell equations. In contrast to the Vlasov–Poisson equations, we need to take into account the influence of the magnetic field. The obtained differential equations, which describe the dynamics on the low-rank manifold, are advection equations in a lower-dimensional subspace and therefore well suited for applying semi-Lagrangian schemes or Fourier techniques. The proposed method never requires the computation of quantities that lie outside of the low-rank approximation space and therefore has a computational advantage compared to the method in [22].

The coupling of the Vlasov equation with Maxwell’s equations poses additional challenges. In the Vlasov–Poisson equations the values of the density function and the electric field are consistent for the entirety of the time integration process. The value of the charge density obtained from the Vlasov equation is directly used to compute the electric field. In the Vlasov–Maxwell equations, this approach is not possible, since the electric and magnetic fields are themselves coupled. Thus, the charge density obtained from the Vlasov equation is, in general, not consistent with the divergence of the electric field obtained from Maxwell’s equations. This leads to a violation of Gauss’ law. We propose to include a correction strategy in the low-rank algorithm, such that the divergence constraint for the electric field is enforced at each integration step. We follow the ideas in [1] and associate a Lagrange multiplier with the constraint on the electric field.

The paper is structured as follows. In section 2 we introduce the Vlasov–Maxwell equations. We then motivate the use of a low-rank approach in section 3. The proposed dynamic low-rank integrator is explained in detail in section 4. In section 5 we explain a correction technique that enforces the divergence constraint at each time step. The description of the setting for the numerical experiments and the implementation details can be found in section 6. In sections 7, 8, 9, and 10 we illustrate the method by means of some numerical experiments. We consider a Landau-type problem, a two-stream instability, a bump-on-tail instability, and a Weibel instability, respectively. While the basic low-rank integrator succeeds in capturing the main dynamics of the system, we observe that a number of invariants are not preserved and, in particular, Gauss’ law is strongly violated. The numerical comparison between the basic low-rank integrator and the divergence preserving low-rank integrator is also carried out here.

2 The Vlasov–Maxwell equations

In this section, we introduce the Vlasov–Maxwell equations in detail. We restrict our discussion to one particle species, in our case electrons. However, the extension to multiple species is straightforward.

We denote the phase space by $\Omega_{x}\times\Omega_{v}$ and write $(x,v)\in\Omega_{x}\times\Omega_{v}$. We consider here the dimensionless form of the Vlasov equation coupled with Maxwell’s equations

[TABLE]

The unknown $f(t,x,v)$ denotes the electron density at time $t$ at the location $(x,v)$, $E$ and $B$ are the electric and magnetic fields, respectively. In the following we set the speed of light $c$ to the value 1. The charge and current densities are defined as

[TABLE]

respectively. Since we consider only one particle species here, the value $n_{0}=\frac{1}{\lvert\Omega_{x}\rvert}\int_{\Omega_{x}}\rho(t,x)\,\mathrm{d}x$ normalizes the electron charge density with respect to the total charge of the surrounding ions. This corresponds to assuming a uniform background distribution of ions.

The evolution of the electromagnetic field is restricted by two divergence constraints, the so called Gauss laws for electricity (3) and magnetism (5). These equations are automatically satisfied in the continuous case provided that the initial values $E(t=0,x)$ and $B(t=0,x)$ satisfy the same properties. In fact, Gauss’ law (3) follows from the continuity equation

[TABLE]

In other words, Gauss’ law is a direct consequence of the conservation of charge. Further, recalling that $\nabla\cdot(\nabla\times\nobreak\cdot)=0$, the zero divergence condition (5) for $B$ holds. Hence, most Maxwell solvers take into account only the time evolution of the two fields, i.e., (2) and (4), and neglect the divergence constraints. However, one usually has to be careful that these physical constraints are still satisfied in the numerical method. Otherwise, the solution can become inaccurate, and in some situations, even numerical instabilities are known to develop. We refer the reader to [7] and references contained therein. The zero divergence constraint for the magnetic field can be fulfilled by selecting an appropriate space discretization. The discrete differential operator $\widetilde{\nabla}$, a discretization of $\nabla$, has to be chosen such that $\widetilde{\nabla}\cdot(\widetilde{\nabla}\times\cdot)=0$. One possibility is to use spectral approximations. A staggered grid can also be employed, where the components of the electric and magnetic fields are considered at different grid points, the so called Yee discretization; see [40]. However, this is not as easily possible for Gauss’ law for electricity as the charge density couples to the Vlasov equation.

Solving the Vlasov equation is very expensive due to its high dimensionality, amongst other factors. Recently, a low-rank approach was proposed for the Vlasov–Poisson equations; see [15, 14]. The key assumption is that the solution has low-rank structure, i.e., it can be approximated as a combination of a relatively small number of basis functions that only depend either on $x$ or on $v$. Then, a function $f$ of rank $r$ can be written as

[TABLE]

Such an approximation will be employed here.

3 Low-rank representation of the Vlasov–Maxwell equations

Our goal in this section is to analytically investigate a situation in which the solution of the Vlasov–Maxwell equations admits a low-rank representation. This will be helpful for the interpretation of the numerical results that are obtained in sections 7, 8, 9, and 10.

We start by linearizing the Vlasov–Maxwell equations around an equilibrium, as done, for example, in [8]. To that end we introduce quantities $f(t,x,v)=f^{(0)}(v)+f^{(1)}(t,x,v)$, where $f^{(0)}$ is the unperturbed equilibrium, $E(t,x)=0+E^{(1)}(t,x)$, and $B(t,x)=0+B^{(1)}(t,x)$. This assumes that our initial value is a perturbation of the equilibrium density given by $f^{(0)}(v)$, say a Gaussian, and that no electric or magnetic fields are required to maintain the equilibrium. The linearized Vlasov–Maxwell equations are then given by

[TABLE]

where $j^{(0)}=\int_{\Omega_{v}}vf^{(0)}(v)\,\mathrm{d}v$ is the (constant) current of the equilibrium. Here, we have neglected higher order terms. Then we perform a Fourier transform in $x$ and denote the corresponding frequencies by $k$. This gives

[TABLE]

Thus, all of the Fourier modes decouple.

Now, let us assume that

[TABLE]

That is, the initial perturbation of the equilibrium $f^{(0)}(v)$ and of the magnetic field can be written as a linear combination of $m$ Fourier modes. This is not a very restrictive assumption. In fact, most of the initial values for the well known test problems can be written in this form (often even for $m=1$).

Now, if $k\neq 0$ and $k\neq k_{i}$ for all $i\in\{1,\dots,m\}$ we have $\hat{f}_{k}^{(1)}(0,v)=0$. We initialize the electric field using Gauss’ law in Fourier space, i.e.

[TABLE]

Then, the unique solution of (6)–(8) is given by

[TABLE]

Thus, no new Fourier modes are excited (except possibly the constant mode $k=0$) and it follows that the solution of the linearized Vlasov–Maxwell equations can be written as

[TABLE]

with $k_{0}=0$. This is a low-rank approximation of rank at most $m+1$, which is the desired result.

The full nonlinear dynamics of the Vlasov–Maxwell equations can, of course, excite new Fourier modes. If the new modes reach a significant amplitude, the linear analysis conducted above is no longer accurate. However, the analysis is still useful since it explains a behavior that we observe in a number of numerical experiments (see also sections 7, 8, 9, and 10). Namely, that a very small rank is sufficient to almost exactly represent the dynamics of the Vlasov equation in the linear regime. This is the case, for example, for the two-stream instability before the electric energy saturates or for linear Landau damping.

It is also interesting to note that the analysis conducted does not exploit the full potential of the low-rank approximation considered in this paper. More specifically, the low-rank structure of (9) assumes that $X_{i}(x)=e^{ik_{i}x}$. In contrast, the low-rank approximation considered in this paper makes no such stipulation. Consequently, in our algorithm the $X_{i}$ can be arbitrary functions of $x$. More details will be given in the next section. Hence, the class of functions that can be represented by a low-rank approximation is much wider than what the Fourier analysis above would indicate. Thus, being in the linear regime is a sufficient but not necessary condition to obtain a low-rank approximation with small rank. As an example, we note that the the saturation of the electric energy (in the two-stream or Weibel instability) is resolved very well even for a small rank (see sections 8 and 10). Saturation is a nonlinear phenomenon which is not at all captured by the linear theory (i.e. by the dispersion relation).

It is also interesting to note at this point that the solution can have a very high degree of filamentation, but still be low-rank. This fact was explained and then used in [15] to efficiently solve the plasma echo problem in the electrostatic regime. In the above analysis, filamentation is implicitly contained in functions such as $f_{k_{i}}^{(1)}(t,v)$, which can be oscillatory (and thus have large derivatives) in $v$. This is a significant advantage of the low-rank approach, compared to other dimension reduction techniques.

In the next section, we summarize the framework for the construction of a low-rank integrator, we define the low-rank manifold in which the low-rank solution lives and describe its tangent space. We then use these two structures to derive a low-rank integrator.

4 The low-rank projector-splitting integrator

In this section, we explain the low-rank integrator for the Vlasov–Maxwell equations for the 3+3 dimensional case in full detail, following the ideas in the seminal paper [25] and the recent works [15, 11, 14].

First, we summarize the framework for the construction of the integrator. We assume that the particle density function $f$ can be represented in low-rank format in the following way

[TABLE]

where $r$ indicates the rank. The functions $X_{i}(x),\ i=1,\dots,r$ express the dependence of $f$ on the space variable $x\in\Omega_{x}$ and the functions $V_{j}(v),\ j=1,\dots,r$ the dependence of $f$ on the velocity variable $v\in\Omega_{v}$. In the following we drop the range of the summation indices and always assume that they run from 1 to $r$.

We seek an approximation of $f$ in the low-rank manifold, defined as follows

[TABLE]

where $\langle\cdot,\cdot\rangle_{x}$ and $\langle\cdot,\cdot\rangle_{v}$ denote the $L^{2}$ inner products in $\Omega_{x}$ and $\Omega_{v}$, respectively. The corresponding tangent space is

[TABLE]

We aim to approximate the dynamics of the Vlasov equation (1) by constraining it to the manifold $\mathcal{M}$. This is done by means of an orthogonal projection onto the tangent space $\mathcal{T}_{f}\mathcal{M}$. The projector $P(f)$ is given by the following expression

[TABLE]

By introducing the following two vector spaces $\bar{X}=\text{span}\{X_{1},\dots,X_{r}\}$ and $\bar{V}=\text{span}\{V_{1},\dots,V_{r}\}$, the projection can be written as

[TABLE]

In [25] the authors proposed a splitting integrator for solving the projected equation

[TABLE]

We split the right hand side into three terms and solve the arising subproblems:

[TABLE]

In the following we explain in detail the first-order scheme, the so called Lie splitting. The second-order integrator is described in algorithm 1. At first we consider the projector $P_{\bar{V}}$ and look at subproblem (12). We represent the approximate particle density $f$ as

[TABLE]

The initial value is represented by $X_{i}^{0}(x)=X_{i}(0,x)$, $S_{ij}^{0}=S_{ij}(0)$ and $V_{j}^{0}(v)=V_{j}(0,v)$. Then, equation (12) can be rewritten as

[TABLE]

By making use of the orthogonality conditions $\langle V_{j},V_{l}\rangle_{v}=\delta_{jl}$ and $\langle V_{i},\dot{V}_{j}\rangle_{v}=0$ we obtain the following differential equation for the coefficients in the basis expansion

[TABLE]

Moreover, from the previous calculation it also follows that $V_{j}(t,v)=V_{j}^{0}(v)$. The coefficients $c_{jl}^{1,2}\in\mathbb{R}^{3}$ are given by the following expressions

[TABLE]

The last term in (15) was rewritten by means of the relation between scalar and vector products $(a\times b)\cdot c=(c\times a)\cdot b$. The coefficient $c_{jl}^{3}\in\mathbb{R}^{3}$ is

[TABLE]

Thus, we integrate equation (15) with initial value $K_{j}(0,x)=\sum_{i}X_{i}^{0}(x)S_{ij}^{0}$ to obtain $K_{j}(\tau,x)=K_{j}^{1}(x)$, where $\tau$ denotes the integration step. After performing a QR-decomposition

[TABLE]

we obtain orthonormal $X_{j}^{1}(x)$, and $S_{ij}^{1}$.

The subproblem (13) is treated in a similar way. In this step the particle function is represented as $f(t,x,v)=\sum_{i,j}X_{i}(t,x)S_{ij}(t)V_{j}(t,v)$ and only $S_{ij}$ is updated. The factors $X_{i}(t,x)$ and $V_{j}(t,v)$ are unchanged, i.e., $X_{i}(t,x)=X_{i}^{1}(x)$ and $V_{j}(t,v)=V_{j}^{0}(v)$. The corresponding differential equation is

[TABLE]

which has to be integrated with initial value $S_{ij}(0)=S_{ij}^{1}$. After one integration step we obtain $S_{ij}(\tau)=S_{ij}^{2}$. The coefficient $d_{ik}^{2}\in\mathbb{R}^{3}$ is given by

[TABLE]

whereas the coefficients $d_{ik}^{1,3}\in\mathbb{R}^{3}$ are defined as follows

[TABLE]

For the solution of (14) we represent the particle function as

[TABLE]

The corresponding differential equation is

[TABLE]

We integrate it with initial value $L_{i}(0,x)=\sum_{j}S^{2}_{ij}V_{j}^{0}(v)$ and obtain $L_{i}(\tau,v)=L_{i}^{1}(v)$. After a QR-decomposition we get orthonormal $V_{j}^{1}(v)$, and $S_{ij}^{3}$ such that

[TABLE]

The value of the particle function after one time step $f(\tau,x,v)$ is approximated by

[TABLE]

Note that the three differential equations (15), (16) and (17) are solved by a favorable numerical scheme, where the electric and magnetic fields are evaluated at the beginning of the time step.

Further, the solution of the evolution equations for the electric and magnetic fields, (2) and (4) respectively, has to be computed. We remark that the space discretization has to be chosen such that the zero divergence constraint (5) is satisfied. We suggest either a spectral discretization or the Yee space discretization, see [40], and denote the discretization of the operator $\nabla$ by $\widetilde{\nabla}$. A first-order scheme is sufficient for the time integration, for example

[TABLE]

where $E(0,\cdot)=E^{0}$, $B(0,\cdot)=B^{0}$, and the current density is computed as follows

[TABLE]

In algorithm 1 we illustrate a second-order scheme. The solution of the Vlasov equation is approximated by means of the Strang splitting scheme. The solution of Maxwell’s equations is computed on a staggered grid in time; see e.g. [40].

5 A divergence correction technique

The low-rank integrator illustrated in the previous section does not explicitly take into account Gauss’ law  (3). Such an approach is motivated by the fact that this divergence constraint is automatically satisfied in the continuous case. This follows immediately from the continuity equation, as already explained in section 2. However, at the discrete level the charge is not conserved and the divergence constraint is not satisfied. Thus, we propose to explicitly include Gauss’ law in the numerical scheme.

We employ an approach based on Lagrange multipliers, following the ideas in [1]. The evolution of the electric field is then described by

[TABLE]

The function $\nabla\phi$ is a correction potential, which is interpreted as the Lagrange multiplier. This approach turns out to be a particular case of a more general idea proposed in, e.g., [32, 7]. There, the evolution equation for the electric field (20) is coupled to the corresponding constraint via the function $\phi$. The Gauss law is modified to

[TABLE]

where $\mathcal{D}(\phi)$ is a differential operator. Depending on the choice of such an operator, different techniques arise. In particular, if $\mathcal{D}(\phi)\equiv 0$ the above system is obtained.

The numerical solution of (20) coupled with (21) is computed by an operator splitting technique. We first solve the subflow

[TABLE]

This is equivalent to the original equation and thus any suitable numerical method can be applied. We explain here only a second-order scheme where we make use of a staggered grid in time; see [40]. As before we denote the discretization of $\nabla$ by $\widetilde{\nabla}$, such that $\widetilde{\nabla}\cdot(\widetilde{\nabla}\times\cdot)=0$. We compute first $\bar{E}^{n+1}$ as the solution of (22) after one time step of size $\tau$ and get

[TABLE]

The intermediate value of the electric field $\bar{E}^{n+1}$ is then updated by solving the second subflow

[TABLE]

coupled to Gauss’ law. We obtain the equation

[TABLE]

and

[TABLE]

Next, we can determine $\phi$ by solving the Poisson equation

[TABLE]

This equation is obtained by applying the discrete divergence to (25) and substituting the discrete divergence condition (26). Then, the electric field is updated as in (25). This correction step for the electric field has to be incorporated in the low-rank scheme described in algorithm 1. The resulting second-order integrator is summarized in algorithm 2.

Note that the correction potential $\phi$ is constructed by taking into account the specific numerical integrator we use for the solution of (22). The error in Gauss’ law made in (23) is then corrected. It is crucial to consider the numerical value of the mass instead of the analytic one in the discrete divergence constraint (26), since the employed low-rank integrator is not mass preserving. Therefore, the divergence constraint has to be determined at the discrete level in terms of the current numerical values of mass and charge density.

Note that the divergence correction procedure and the quasi-conservative scheme proposed in [14] are complementary approaches. The conservation of mass is enforced on the distribution function, whereas the correction proposed here only acts on the electric field; the distribution function only enters as an input via the charge density and the current density. Thus, both procedures can be combined.

6 Implementation details

In the numerical experiments reported in the following sections, we consider the reduced model obtained by taking into account one dimension in the physical space and two in the velocity space. We describe this model here. To do that, we write $v=(v_{1},v_{2})$ and denote by $v_{1}$ and $v_{2}$ the two components of the velocity. The phase space is then $\Omega_{x_{1}}\times\Omega_{v_{1}}\times\Omega_{v_{2}}$. The electric field $E=E(t,x_{1})=(E_{1},E_{2})(t,x_{1})$ is two dimensional, where $E_{1}(t,x_{1})$ and $E_{2}(t,x_{1})$ denote its components. The magnetic field is orthogonal to the electric field and is one dimensional. We thus have

[TABLE]

The Vlasov equation reduces to

[TABLE]

where the matrix $\mathcal{J}$ is given by

[TABLE]

Further, the Maxwell system is given by

[TABLE]

Note that the zero divergence condition for the magnetic field is automatically fulfilled.

The Vlasov–Maxwell systems conserves the total energy. In the next sections, we study the behaviour of the proposed integrator with respect to this matter. In this respect, we recall here that the electric and magnetic energies are the energies stored in the electric and magnetic fields, respectively. They are defined as

[TABLE]

The kinetic energy is given by

[TABLE]

The projector-splitting integrator illustrated in section 4 is slightly simpler in the present case. The coefficients $c_{jl}^{1,3}$ are real numbers, whereas $c_{jl}^{2}$ is a vector in $\mathbb{R}^{2}$. Similarly, $d_{ik}^{1,3}$ are real numbers and $d_{ik}^{2}$ is in $\mathbb{R}^{2}$. The equations for the substeps can then be easily obtained from (15), (16), and (17).

The numerical discretization of the spatial and velocity operators is carried out in spectral space. We use the FFTW library [16] for computing the discrete Fourier transform in one dimension. Given $n$ real values $X_{k}$, $k=0,\dots,n-1$ the discrete Fourier transform is defined as follows

[TABLE]

where $\widehat{X}_{k}=\widehat{X}_{n-k}^{*}$ with $*$ denoting the complex conjugate. Due to this symmetry, only half of the output data has to be stored in computer memory. Moreover, it follows that $\hat{X}_{0}$ and $\hat{X}_{n/2}$, for $n$ even, are real numbers. The imaginary part of these numbers is not stored. This introduces an additional numerical error in the solution of the Poisson equation as we explain in the following.

For the time integration of (15), (16), and (17) we employ the boost::odeint library. The DOPRI5 (Dormand–Prince) method of order $5$ is used. In order to keep the error small and to guarantee stability, we use $5$ substeps for each step of the splitting scheme.

For computing the correction to the electric field the Poisson equation (27) has to be solved first. In spectral space the solution is given by

[TABLE]

where $g$ denotes the right-hand side of (27). Then the space derivative of $\phi$ has to be computed. In spectral space it is done by multiplying $\hat{\phi}_{k}$ by $ik$, i.e.

[TABLE]

Therefore, for $n$ even, the entry $\widehat{\phi}_{n/2}$ is a purely imaginary number. However, due to the symmetry this information can not be stored in memory and thus a numerical error is introduced by the inverse transformation. Normally these errors can be neglected. Here, however, the errors grow as the solution evolves in time. An easy remedy is to use an odd number of grid points in the real space $\Omega_{x}$. We will do this in all the simulations conducted in the following sections.

Let us discuss the computational complexity and storage requirements of the proposed integrator. If $r$ denotes the chosen rank, then for each time step we have to solve

• •

a system of $r$ advection equations in $x$ (this requires the computation of at most $3d_{v}r^{2}$ integrals over $\Omega_{v}$),

• •

a system of $r^{2}$ ordinary differential equations (this requires the computation of at most $3d_{v}r^{2}$ integrals over $\Omega_{x}$),

• •

a system of $r$ advection equations in $v$.

Additionally, one has to consider the cost for solving Maxwell’s equations. These are equations posed in $x$ only. In case of the divergence correction the cost of solving the Poisson equation has to be counted also. However, since we do this in Fourier space this constitutes only a negligible amount of the total run time of the algorithm. If $n_{x}$ and $n_{v}$ degrees of freedom are used in each coordinate direction in $x$ and $v$, respectively, then the storage cost is reduced from $\mathcal{O}(n_{x}^{d_{x}}n_{v}^{d_{v}})$ to $\mathcal{O}(r(n_{x}^{d_{x}}+n_{v}^{d_{v}}))$. As we are using spectral methods, the cost of solving the evolution equations is proportional up to a logarithm to the number of degrees of freedom. In this case our low-rank algorithm requires $\mathcal{O}(r(n_{x}^{d_{x}}+n_{v}^{d_{v}}))$ arithmetic operations for the evolution equations and $\mathcal{O}(r^{2}(n_{x}^{d_{x}}+n_{v}^{d_{v}}))$ operations for computing the integrals.

7 Numerical results for the Landau-type problem

We consider first a Landau-type problem (in the sense that this configuration shows Landau damping when the speed of light goes to infinity). This is the same setup as in [4]. We consider the following initial value

[TABLE]

where we choose $k=0.4$ and $\alpha=0.01$. The electric field is initialized by solving Gauss’ law (31). This yields

[TABLE]

and the magnetic field at $t=0$ is chosen as

[TABLE]

The space domain here is $\Omega_{x_{1}}=[0,2\pi/k]$ and the velocity domain is chosen as $\Omega_{v_{1}}\times\Omega_{v_{2}}=[-5,5]^{2}$.

In figure 1, on the top, we show the results obtained with algorithm 1 (without correction). On the left, we plot the values of the electric and magnetic energies for rank = 15. Our results are very similar to those obtained with an Eulerian Vlasov solver in [4]. On the right we plot the time evolution of the smallest singular value of the distribution function $f$. The smallest singular value of the matrix $S=(S_{ij})$ gives a measure of the error introduced by truncating the rank. We observe that this error stays bounded in time for all the approximation ranks we tested. Moreover, the choice of a higher approximation rank gives better results, as expected.

In the third plot in figure 1 we display the time evolution of the electric and magnetic energies for the divergence preserving scheme given in algorithm 2. We observe that the results are very close to the uncorrected results.

In figure 2 we compare the behaviour of the integrator given in algorithm 1 (dashed-dotted lines) and the divergence preserving scheme of algorithm 2 (full lines) in terms of mass and energy conservation, and preservation of Gauss’ law. The first plot shows the relative error in mass

[TABLE]

where the mass $m$ is defined as

[TABLE]

The relative error in the total energy $\mathcal{E}_{tot}=\mathcal{E}_{e}(t)+\mathcal{E}_{m}(t)+\mathcal{E}_{k}(t)$, i.e.

[TABLE]

is shown in the second plot. Conservation of mass and energy is quite good for both algorithms. For the case of algorithm 1 we observe that mass conservation and the error in Gauss’ law improve as we increase the rank indicating that the error is dominated by the low-rank approximation. We also observe that among the three invariants considered here, the largest error, approximately $10^{-3}$ for rank $15$, is committed in Gauss’ law. When applying the divergence preserving scheme of algorithm 2 Gauss’ law is satisfied up to machine precision. We further note that the divergence correction does not cause a significant difference in the behavior of the other invariants.

8 Numerical results for the two-stream instability

In this section we consider the two-stream instability. We employ the problem found in [3]. The initial particle density is given by

[TABLE]

where we choose $\beta=2\cdot 10^{-3}$, $\Omega_{x_{1}}=[0,2\pi]$ and $\Omega_{v_{1}}\times\Omega_{v_{2}}=[-0.4,0.4]^{2}$. The instability is driven by a magnetic perturbation only. The initial magnetic field is chosen as

[TABLE]

with $\alpha=10^{-3}$. The electric field is initialized to zero.

In figure 3, on the top, we plot the time evolution of the energies obtained with algorithm 1 (without correction) for different values of the approximation rank. Initially, we observe an oscillatory behavior of the electric and magnetic energies. Later, the electric energy increases exponentially and then saturates. The results for rank = 5 are on the left, on the right the results for rank = 15 are displayed. We clearly notice the benefit of choosing a higher approximation rank. In particular, the electric and kinetic energies show a better match after saturation. This might indicate that a higher approximation rank is necessary to resolve the nonlinear regime. On the other hand, the plot of the smallest singular values brings a different conclusion. The error in the distribution function in the nonlinear regime is relatively large (independent of the choice of the rank). Despite this large error in the distribution function, we observe that averaged quantities of the system, such as the electric energy, match the physical dynamics quite well.

In figure 3 we also show that the splitting integrator described in algorithm 1 converges in time with the expected order. We display the time error of the scheme for different values of the time step $\tau$ with respect to a reference solution.

In figure 4 we show the behavior of the integrator with respect to the conservation of mass, energy, and the preservation of Gauss’ law. The numerical error of the proposed integrator is an interplay between different factors. We consider the influence of the choice of the approximation rank, the spatial discretization, and the time integration error due to the projector splitting integrator. On the left, the error in mass is plotted. First, let us focus on the curves corresponding to different approximation ranks obtained with time step size $\tau=0.1$. We observe that increasing the rank gives more accurate results. In particular, the method preserves the mass up to machine precision for the full rank (i.e. rank = 33 for the configurations with $n_{x_{1}}=33$) until saturation. Note, however, that in the nonlinear regime the error increases in time. This is due to the fact that even if the dynamics can be exactly represented by the chosen rank, there is en error comitted by the projector splitting integrator. In this regime, conservation of mass can be enforced by choosing a smaller time step (we consider $\tau=0.01$ in the figure).

Similar conclusion can be drawn by looking at the error in energy shown on the right. In this case we can further observe that a finer space discretization is beneficial in the nonlinear regime. In fact, the red and black curves overlap until saturation. Afterwards, the error for the fine grid is smaller compared to the other configurations. The plot of the error in Gauss’ law, also, confirms the previous observations. If the approximation rank is high enough, the time integration error can be clearly observed. This error decreases by reducing the time step size from $\tau=0.1$ and $\tau=0.01$ by two orders of magnitude, as expected for a second-order integrator.

In figure 5, top left, we observe that the time evolution of the energies for the divergence preserving scheme are close to what is obtained with algorithm 1. As for the Landau-type problem, the divergence correction does not cause a significant difference in the behavior of the other invariants.

9 Numerical results for the bump-on-tail instability

The aim of this section is the study of the bump-on-tail instability. We choose an initial value similar to the one in [6]:

[TABLE]

The following numerical results are obtained with $\alpha=9/10$, $\beta=2/10$, $\gamma=0.03$ and $k=0.3$. We consider the phase space $\Omega_{x_{1}}\times\Omega_{v}$ to be $[0,20\pi]\times[-9,9]^{2}$. The electric field is initialized according to Gauss’ law and no initial magnetic field is considered. The perturbation of the equilibrium results in an instability which consists of three vortices traveling periodically in phase space. The numerical simulation of this phenomenon is challenging. The three vortices must be kept separate and the filamentation must be accurately resolved. We report our results for algorithm 1 (without correction) in figure 6 where we compare the distribution function at different times. In the upper row, the results obtained for rank = 5 are shown. The qualitative behavior of the solution improves significantly when choosing rank = 15. In both configurations, the formation of the vortices is resolved very well.

10 Numerical results for the Weibel instability

In this section we consider an instance of the Weibel instability; see [3]. The initial particle density is

[TABLE]

where we choose $v_{th}=0.02$, $T_{r}=12$, $k=1.25$, and $\alpha=10^{-4}$. The magnetic field at $t=0$ is chosen as

[TABLE]

with $\beta=10^{-4}$. The electric field $E_{1}$ is initialized by (31). The space domain here is $\Omega_{x_{1}}=[0,2\pi/k]$, and the velocity domain is chosen as $\Omega_{v_{1}}\times\Omega_{v_{2}}=[-0.3,0.3]^{2}$.

As explained in section 6, we approximate the spatial derivatives in Fourier space. This results in a numerical scheme that does not have any dissipation. This is the correct behavior from a physical point of view. However, in some configurations this can cause numerical instabilities. In order to avoid this problem, we add some artificial dissipation into Maxwell’s equations as follows

[TABLE]

By tuning the value of the parameter $\varepsilon$, we can adjust the strength of the numerical dissipation. We use $\varepsilon=10^{-2}h^{2}$ in all simulations in this section, where $h$ is the spatial grid size. This is sufficient to remove the numerical instability. Such an approach is commonly used and dates back to [38]. The need of introducing artificial dissipation is due to the fact that here, contrary to the previous examples, the magnetic field is strong and thus non-physical oscillations are amplified by the numerical scheme. This is not an issue with the low-rank approximation. Eulerian schemes based on spectral methods, e.g. the scheme proposed in [3], also require artificial dissipation. In the implementation, we treat the diffusion implicitly.

In figure 7, on the top, the first Fourier mode of the electric and magnetic fields, respectively, are plotted for algorithm 1 (without correction). We compare the results obtained for rank = 5 (left) and rank = 15 (right). We observe that the integrator works very well in both the linear regime and to simulate saturation. More specifically, it reproduces the exponential growth of the magnetic field until the saturation point. We refer to [3, sect. 5.2] for some reference results and conclude that both choices of the approximation rank give good qualitative results also in the nonlinear regime. On the other hand, the behavior of the smallest singular value (bottom right figure) indicates that the error in the distribution function is quite large in the nonlinear regime. Thus, we once again have a situation in which the averaged quantities are resolved much better than the distribution function itself.

Similar conclusions can be drawn from figure 8. Mass and energy are well conserved for both the schemes in the case of short time integration, but a significant error is incurred in the nonlinear regime. However, the divergence correction influences the overall behavior of the scheme in the nonlinear regime. In the top left figure we observe worse qualitative behavior of the Fourier modes, in particular for the electric field. Enforcing Gauss’ law results in a deterioration of the quality of the other invariants we considered. Let us explain this behavior. In general, there is no guarantee that a small correction is sufficient to obtain a point on the low-rank manifold which satisfies Gauss’ law. Thus, there is a competition between the approximation error and the divergence correction. This is precisely what happens in the nonlinear regime. Therefore, introducing the divergence correction deteriorates the qualitative behavior of the numerical solution. This is primarily a problem in situations where the low-rank error in the distribution is large. We also remark that the divergence correction procedure explained in section 5 has to be adapted to the damped equations for the electromagnetic field (35)–(37).

Acknowledgements

We would like to thank Bruno Després (Laboratoire Jacques-Louis Lions) for the helpful discussions. We also thank the referees for their helpful comments, which improved the presentation of this paper.