A mass, momentum, and energy conservative dynamical low-rank scheme for the Vlasov equation

The primary challenge in solving kinetic equations, such as the Vlasov equation, is the high-dimensional phase space. In this context, dynamical low-rank approximations have emerged as a promising way to reduce the high computational cost imposed by such problems. However, a major disadvantage of this approach is that the physical structure of the underlying problem is not preserved. In this paper, we propose a dynamical low-rank algorithm that conserves mass, momentum, and energy as well as the corresponding continuity equations. We also show how this approach can be combined with a conservative time and space discretization.


Introduction
Solving kinetic equations efficiently is important in applications ranging from plasma physics to radiative transfer. The main challenge in this context is the up to six-dimensional phase space and the associated unfavorable scaling of computational cost and memory requirements, usually referred to as the curse of dimensionality. Because of this, particle methods, such as the particle in cell (PIC) scheme, have been and are still widely used. However, it also well known that particle methods miss or do not resolve certain physical phenomena (see, e.g., [2]) or require an immense number of particles thereby negating their advantage. In addition, complexity reduction techniques such as sparse grids have been investigated. While they can provide some advantage compared to a full grid simulation the gain is usually modest [20].
More recently, using dynamical low-rank approximations to solve kinetic problems has received considerable interest. Such methods have been developed for both the Vlasov equation [19,11,12,13] and radiation transport problems [8,28,10]. Dynamical low-rank integrators approximate a six-dimensional Vlasov equation by a set of only three dimensional advection problems. Moreover, they have a range of properties that makes them well suited for performing kinetic simulations. In particular, such methods are able to resolve filamentation and are almost exact if the dynamics can be well represented by a linearized equation [13]. In addition, dynamical low-rank schemes can capture the limiting fluid or diffusive regime [9,8,10,29]. Therefore, dynamical low-rank integrators can drastically decrease the numerical effort that is required to solve a number of kinetic problems. This enables the simulation of such problems on desktop computers or small clusters that are otherwise either unfeasible or require large supercomputers.
A major disadvantage of dynamical low-rank integrators, however, is that they do not respect the physical structure of the underlying equations. In particular, mass, momentum, and energy are not preserved by the low-rank approximation. This is in stark contrast to Eulerian and semi-Lagrangian Vlasov solvers, where at least mass and momentum are conserved [27] and methods with good long-time behavior with respect to energy have been obtained [7], and particle methods, where usually mass and either momentum or energy can be conserved [32].
Some approaches to improve this deficiency have been proposed. In [28] the solution is simply rescaled such that the mass is preserved. Such an approach, however, can not be extended to simultaneously conserve momentum or energy. It also does not respect the underlying continuity equation for the mass density, which could actually be more important for the long time behavior of the integrator [12]. In [12] a method based on Lagrange multipliers is proposed. This approach succeeds in improving the conservative properties of the low-rank approach, but it does not allow us to simultaneously conserve both the continuity equations for density and momentum as well as the corresponding invariants. It also does not change the dynamical lowrank approximation. Instead, this method adds a correction to each step of the low-rank integrator. In [29] the fluid moments are integrated explicitly. These moments are then coupled to a low-rank approximation that resolves the kinetic dynamics. To enforce conservation, a correction is added to the low-rank part of the algorithm that requires, similar to [12], the solution of a linear system of equations.
Neither of the conservative or quasi-conservative methods developed in the literature solve the fundamental problem. Namely, that the classic dynamical low-rank approximation does not take the structure of the equations that we are trying to solve into account. In this paper we introduce such an approach that conserves both mass, momentum, and energy and ensures that the corresponding continuity equations are satisfied as well. The proposed method is based on the observation that if certain functions of velocity belong to the approximation space, then the desired conservation follows. To accomplish this two steps are necessary. First, we have to formulate the dynamical low-rank scheme in a different function space than the L 2 space that is usually used (we will use an appropriately weighted L 2 space instead). Second, we have to constrain the low-rank factors such that the desired functions belong to the approximation space at all times. This is done by fixing certain basis functions and using a modified Petrov-Galerkin condition in such a way that this is compatible with the remainder of the dynamical low-rank approximation. This results in equations of motions that are somewhat different from the ones used in the classic dynamical low-rank algorithm introduced in [18] (and that subsequently was used heavily in the literature; see, e.g, [26,25,11,12,28,13,10]). We also introduce a time and space discretization of the resulting equations of motions that conserves mass and momentum as well as a discrete version of the corresponding continuity equations up to machine precision.
The remainder of the paper is structured as follows. In section 2 we introduce the Vlasov-Poisson equation and discuss the invariants and corresponding physical structure of this model. This is followed by a recollection of the classic dynamical low-rank integrator in section 3. The proposed conservative dynamical low-rank integrator is then introduced in section 4 and its properties are discussed in section 5. The fully discretized integrator, i.e. time and space discretization, is discussed in section 6. Finally, a number of numerical simulations are presented in section 7.

Vlasov-Poisson equations and conservation
In this work we consider the Vlasov-Poisson equations The sought after quantities are the particle-density function f and the electric field E. This equation models the dynamics of electrons that are subject to a homogeneous ion background distribution. We consider the single species Vlasov-Poisson equation here, but it should be emphasized that the described algorithm can be easily extended to the multi-particle case as well as to more complicated models, such as the Vlasov-Maxwell equations.
The dynamics of the Vlasov-Poisson equations conserves a number of physically relevant invariants. In particular, mass M , momentum J and energy E are invariants of the solution to the continuous equations. We note that in the single species setting mass conservation is equivalent to charge conservation (since charge and mass are proportional) and momentum conservation is equivalent to current conservation. For each of these invariants, a continuity or moment equation, that is posed in physical space only, is satisfied by an associated density. In the case of mass, the mass density Conservation of mass can be easily derived from the continuity equation by integrating in x. It is in fact these continuity equations that our dynamical low-rank integrator satisfies and from which conservation of the global invariants follows. We note that this is a much stronger result than simply conserving the invariants, as the present approach also preserves the underlying physical structure of the equation. The associated density for momentum conservation is the momentum density j(t, x) which satisfies the following continuity equation We can derive conservation of momentum by recognizing that E(1 − ρ) = ∇ · (E ⊗ E − 1 2 E 2 ) and integrating in physical space. Note that due to the normalization of the particle-density function we have E dx = 0. The energy density satisfies the following continuity equation We derive global conservation, i.e. conservation of energy, by The last equality follows since ∂ t E(t, x) = j(t, x) is just the electrostatic version of Ampere's law. A more detailed discussion can be found in [30]. Looking at the conserved quantities in light of their corresponding continuity equations also relates the kinetic equation considered here to their corresponding fluid model. For more details we refer the reader to [16,10,9].

The classical dynamical low-rank scheme
For a dynamical low-rank scheme we seek an (approximate) solution of the form where r is the rank of the approximation. The dynamics is represented by the low-rank factors X i and V j , that only depend on physical space x and velocity space v, respectively, and the low-rank factor S ij , which carries no spatial or velocity dependence. From a physical point of view, the decomposition into physical space and velocity space is natural. If the dynamics is integrable then there are conserved action variables j(x, v) and conjugate angle variables θ(x, v).
In this case, the natural basis functions are Θ k (θ, t) = e ikθ and J k (j, t) = δ(j −j k ). Note, however, that there is a nonlinear coordinate transformation that must be determined in order for this representation to hold. The separation into physical space and velocity used here does hold for linearized equations, such as the Case-Van Kampen eigenfunctions [31] for the linear Vlasov-Poisson equation. For a mathematical perspective in the context of the dynamical low-rank approximation, see [13]. In the collisional case, the nonlinear collision operator generates scattering in v and is considered local in x. This simplification would not hold in generalized phase space coordinates {θ, j}. Hence, when collisions are dominant, the decomposition between physical space and velocity space is still the natural one. Due to the fact that sharp structures, such as shocks and boundary layers, can form in the x direction, good basis functions are localized in x. For the v direction, good basis functions include functions orthogonal over a Maxwellian distribution and eigenfunctions of the collision operator. The low-rank approximation can accommodate this without any difficulty, as has been demonstrated in [9,8,10]. From now on, we will, for reasons of simplicity, write f for f (t, x, v), X i for X i (t, x), etc. The equations of motions for the low-rank factors X i , S ij , and V j are then determined by imposing a Galerkin condition, see [18]. Following this argument for the PDE case we obtain [11] where For the Vlasov-Poisson equation, (3). We note that the low-rank factors X i and V j are orthonormal; that is, they satisfy ( The main utility of the low-rank approximation is that instead of a problem in dimension 2d, as in equation (1), we only have to solve 2r equations of dimension d and r ordinary differential equations. This is the reason why the dynamical low-rank approximation drastically reduces the memory and computational effort required to solve such problems, assuming that r is not too large, which is often true in practice. For a more detailed discussion we refer the reader to [11,13] and the subsequent discussion in section 4. Equations (4)-(6), however, do not preserve the invariants of the original system. That is, even if the system is solved exactly (i.e. no time or space error is introduced) conservation of mass, momentum, and energy is lost. This is not very surprising as the dynamical low-rank approximation described above does not take the physical structure of the equations into account. Thus, there is no mechanism that prevents, e.g., mass to be removed from the system due to the truncation performed by the low-rank approximation.

The conservative dynamical low-rank scheme
In this section, we propose a novel dynamical low-rank integrator with equations of motions, that will take the place of (4)- (6), that ensure mass, momentum, and energy conservation.
For the classical dynamical low-rank integrator we have and thus by integrating over v where we have used equations (4)- (6). Now, if we could ensure that 1 ∈ V , where V = span{V j }, then we could find coefficients c j such that 1 = j c j V j . Plugging this into equation (7) we get which would imply mass conservation. However, 1 ∈ V is clearly not true as constant functions do not lie in L 2 (R d ). Thus, the continuity equation is not satisfied and the classic dynamical low-rank integrator does not conserve mass. A very similar argument can be made for the momentum, with v ∈ V , and energy, with v 2 ∈ V . One might object at this point that in a numerical simulation we necessarily use a truncated domain and thus the functions 1, v, and v 2 lie in L 2 (Ω v ). However, on a finite domain the approximation spaces have to be equipped with appropriate boundary conditions (usually either periodic or homogeneous Dirichlet conditions are imposed). The functions 1, v, and v 2 do not satisfy boundary conditions that are compatible with the problem. In the periodic case it is clear that v and v 2 are not periodic and assuming so would incur a large numerical error. In the case of homogenous Dirichlet boundary conditions 1, v, v 2 are not zero at the boundary and thus do not lie in the desired approximation space.
To obtain a conservative dynamical low-rank integrator we proceed as follows. First, we will use a weighted function space that includes the possibility of representing constant functions in V . However, this on its own is not yet sufficient. In fact, the dynamic low-rank algorithm automatically chooses appropriate basis functions in order to minimize the overall error, according to some Galerkin condition, in the particle-density function f . There is no guarantee that such a choice satisfies 1 ∈ span{V }. Thus, we have to constrain the approximation space in such a way that constant functions are always part of the basis. This, in particular, requires us to modify the Galerkin condition that is used in the classic integrator. Those two ideas combined allow us to formulate a conservative dynamical low-rank algorithm. The details of this procedure will be the topic of the remainder of this section.
For the conservative dynamical low-rank scheme the solution is approximated by a function of the following form where r is the rank of the approximation and f 0 ( ) is a, yet to be chosen, weight function. The low-rank factors X i (t, ·) and V j (t, ·) are assumed to lie in the L 2 spaces weighted by f 0x and f 0v , respectively. That is, . The corresponding weighted inner products are denoted by respectively. In the following we will choose f 0x (t, The choice of f 0v must guarantee that 1, v , and v 2 lie in L 2 (Ω v , f 0v ) in order to obtain conservation of mass, momentum, and energy, respectively. It is also important to choose the temperature large enough such that the relevant part of the phase space is captured. For many problems that model plasma instabilities is a good choice as this represents the equilibrium distribution and guarantees that all powers of v lie in the approximation space. However, other than the constraints outlined above the choice of f 0v is arbitrary.
The second crucial ingredient of the algorithm is that some of the functions V j are held fixed as the system evolves in time. We write where the U a (v) are fixed, i.e. they are not changed by the dynamical low-rank integrator, and the W p (t, v) are allowed to vary in time according to the low-rank algorithm. In the following, indices i, j, k, l span 1, . . . , r, indices p, q span m + 1, . . . , r, and the indices a, b span 1, . . . , m. For example, to obtain mass, momentum, and energy conservation in 1+1 dimension and for f 0v It is clear that then the orthogonality constraint U a , U b v = δ ab is satisfied. As usual, we also impose the orthogonality conditions W p , V j = δ pj . Note that this implies that the W p are orthogonal both with respect to each other and with respect to the U a . The equations of motion for the low-rank factors are derived by considering the low-rank manifold of all functions of a given rank r. In our case this manifold M can be written as As usual we will impose the gauge conditions ∂ t X i , X j x = 0 and ∂ t W p , W q v = 0 in order to make sure that the dynamics of the low-rank factors is uniquely determined. The tangent space of the manifold is then where we have used the fact that ∂ t U a = 0. While this construction is similar to the classic dynamical low-rank integrator, it is crucial that we are cognizant for which low-rank factors the gauge conditions are imposed and that we make sure that the entire set of the V j are orthogonal to each other (and not only the W p ).
We now derive the equations of motion for the low-rank factors X i , S ij , and V ij that satisfy the following Petrov-Galerkin condition where (f, g) xv = (f g) xv and ∂ t f ∈ T f M. Note that (·, ·) x , (·, ·) v , and (·, ·) xv denote the usual L 2 inner products, while ·, · x , ·, · v , and ·, · xv denote the L 2 inner products weighted with f 0x , f 0v , and f 0 , respectively. Since we have assumed f 0x (t, x) = 1, it holds that (·, ·) x = ·, · x . The Petrov-Galerkin condition is different from the classic dynamical low-rank scheme; an additional factor of 1/f 0v has been introduced. The reason for this is that in the following our goal is to choose an appropriate v ∈ T f M such that the equation of motion for a specific low-rank factor is isolated. This change in the Petrov-Galerkin condition is what makes this possible in the weighted approximation spaces that we consider here.
Equations for X i : In this case there is little difference to the classic algorithm. We consider a family of test The Petrov-Galerkin condition (9) then becomes We can rewrite this as and thus by using the orthogonality and gauge conditions as well as the fact that χ is arbitrary we obtain which is precisely the first equation of motion that is obtained for the classic algorithm (see, e.g., [11]). Now, we can plug the right hand-side of the Vlasov-Poisson equation into the right-hand side of equation (10). This yields We note that due to the weight function f 0v the coefficients c 1 kl and c 2 kl are changed compared to the classic algorithm in [11]. Equations for W p : In this case we consider the family of test functions ν q = f 0v ζ(v) i X i S iq , where ζ is an arbitrary function of v. Those ν q lie in the tangent space as we can write ν q = f 0v ip X i S ipẆp witḣ W p = δ pq ζ(v). The Petrov-Galerkin condition (9) then becomes and thus (since ζ is arbitrary and using the gauge conditions) On the left-hand side we have the matrix T qp = i S iq S ip . Since S has full rank the same is true for T and thus we can invert it in order to obtain the equations of motion for the W p .
For the Vlasov-Poisson equation we have

Equation for S ij :
In this case we consider ν kl = f 0v X k V l , which lies in the tangent space as ν kl = f 0v ij X iṠij V j withṠ ij = δ ki δ jl . The Petrov-Galerkin condition (9) becomes For the Vlasov-Poisson equation Together equations (10)- (12) are the equations of motions for the proposed conservative dynamical low-rank integrator. They take the place of equations (4)-(6) that have been used in the classic algorithm. The primary difference lies in the equation for the W p , where it has to be ensured that the W p are updated in such a way that they remain orthogonal to not only the other W p but also to the fixed U a (this and the desire to keep the U a fixed is the reason why we can not simply choose ν q = χ(v)X q in deriving the equations for the W p above). In addition, due to the use of the weighted L 2 spaces all coefficients are changed as well and f 0v and its derivatives make an appearance in the equations of motion.

Mass, momentum, and energy conservation for the proposed low-rank approximation
In this section we will show that the numerical scheme derived in section 4 is indeed mass, momentum, and energy conservative.
To ensure mass conservation we choose, as is discussed in the previous section, U 1 ∝ 1. Then, using K j = i X i S ij , the density is given by This can be easily seen as Using the orthogonality condition U 1 , V j = δ 1j we obtain the desired relation.
Deriving the corresponding continuity equation is now straightforward. We have From equation (10) we This is precisely the relation that is obtained for the continuous evolution, i.e. the solution of the Vlasov-Poisson equation without any low-rank approximation present. Since From the continuity equation conservation of mass follows at once (by integrating in space).
For conservation of momentum we proceed in precisely the same manner. For simplicity we only consider the 1+1 dimensional case here. By choosing U 2 such that v = v U 2 , i.e. U 2 ∝ v, momentum is given by We thus have which is the desired continuity equation that implies conservation of momentum. The extension to the multi-dimensional case is straightforward. We have to use one fixed basis function for each of the directions in which we want to conserve momentum (e.g.
For energy conservation we can not choose U 3 ∝ v 2 since v 2 is not orthogonal to 1. Thus, we use U 3 ∝ v 2 − 1.
In this setting we have U 3 , U 2 v = U 3 , U 1 v = 0, as desired. We then can represent v 2 as follows Using this relation we have Thus, for the energy density we have Our goal is once again to derive the corresponding continuity equation. We have Evaluating the integral and using the relations for the electric field derived in section 2 we get the desired continuity equation which implies conservation of energy.
We have here outlined an approach with m = 3, U 1 ∝ 1, U 2 ∝ v, and U 3 ∝ v 2 − 1 which conserves mass, momentum, and energy. However, it is also possible to use this method to only conserve one or two out of these three invariants. If only energy conservation is desired, e.g., we can simply choose U 1 ∝ v 2 and the argument above is even simpler as we do not have to take the orthogonality between U 1 and U 3 into account.
We can clearly see that for our algorithm the argument with respect to conservation made for the Vlasov-Poisson equations in section 2 carries over to the low-rank approximation. This is by design and highlights the ability of our approach to preserve the physical structure of the Vlasov-Poisson equations.

Discretization
In the previous section we have established that the proposed dynamical low-rank integrator conserves mass, momentum, and energy. This is done completely within a continuous formulation. Thus, the only approximation made is due to the fact that the proposed scheme uses a low-rank representation of the solution. However, to actually implement the conservative dynamical low-rank integrator on a computer, we have to also introduce a time and space discretization. Devising a numerical method that discretizes the equations of motions for the proposed dynamical low-rank integrator, while maintaining conservation, is the main purpose of this section.
In section 6.1 we will introduce an explicit integrator for equations (10)-(12) that preserves mass and momentum up to machine precision. For dynamical low-rank approximations integrators that are robust to the presence of small singular values have been developed recently [24]. Unfortunately, it turns out that the projector splitting integrator can not be easily adapted to the present situation. In section 6.2 we propose a robust dynamical low-rank integrator that fits within the framework considered in this paper.

Conservative Euler scheme
Simply applying a classic Runge-Kutta scheme to the evolution equations (10)-(12) does destroy the conservative properties of the method. As an example, let us consider the classic explicit Euler scheme where τ is the time step size and the upper indices denote the value of the corresponding quantity at the discrete time t n . The reason why this scheme fails to be conservative is that we use the inverse of S n , i.e. S at time t n , to compute X n+1 . This inconsistency means that there is no well defined K n i and K n+1 i and thus the argument in section 5 can not be applied. However, we can write equation (10) in the following conservative form Applying the explicit Euler scheme to that equation yields No change is made to the other two equations of motions. Putting this together we obtain In the following we will call this method the conservative Euler scheme. This is still a fully explicit method as we can first compute S n+1 using equation (14), which is then used in the computation of X n+1 . The scheme satisfies a discrete version of the continuity equation as can be easily seen from equation (15). Integrating this equation in x we get M n+1 = M n and thus conservation of mass. The conservative Euler scheme also yields the following discrete continuity equation which implies conservation of momentum. For energy we have The reason why we have rewritten the electric field in this particular way should become clear shortly.
Integrating the above relation in space yields We again write the electric field using its potential, i.e. E n = −∇φ n , to get Thus, the first in term on the right-hand side of (17) vanishes. That leaves the second term which can clearly not be zero (except for trivial solutions that satisfy E n = const for all n). Thus, the conservative Euler scheme commits a first order error in the energy. The reason for this is that in the discrete setting the relation ∂ t (E 2 ) = 2E∂ t E does not hold true. This introduces the second term in the calculations above and thus destroys the conservation of energy in the discrete setting. Fundamentally, the issue is that the Euler scheme is not symmetric under time reversal. If one allows for implicit methods this deficiency can be remedied. For example, using E n+1/2 = (E n+1 +E n )/2 in the kinetic update results in Integrating in physical space, as above, then shows conservation of energy. However, since the electric field depends on the particle-density function this approach results in a fully implicit scheme that has to be solved up to machine precision, if energy conservation to the same level of accuracy is desired. We consider the construction of efficient energy conservative time integrators a subject of future research. Let us now turn our attention to the discretization of space. Fortunately, this is relatively straightforward. As long as the discrete approximation of the derivatives and the quadrature rule chosen to define the invariants allows us to perform integration by parts without introducing any error, the calculations made in this section carry over to the fully discretized case. This is true for a number of space discretization strategies. For example, using fast Fourier techniques (FFT) or the standard second-order centered finite difference scheme to compute the derivatives in combination with the trapezoidal rule to evaluate the integrals satisfies this property. The former will be used in the numerical results in section 7. Moreover, for a number of more advanced and higher-order finite difference, finite volume, and discontinuous Galerkin schemes this property is satisfied as well (see, e.g., [1,14,33]). Such methods can then also be used in situations where, e.g. , homogeneous Dirichlet boundary conditions are required. For non-homogeneous boundary conditions additional modifications have to be made to the low-rank integrator (see, e.g., [21]). There is one additional concern in an actual implementation that, by necessity, operates in finite precision arithmetic (i.e. with doubles or floats on a computer). To show conservation we employed the orthogonality condition W p , U a = 0. This is true for the dynamical low-rank integrator by construction (see section 4). The conservative Euler discretization also preserves this property. This can be seen explicitly from equation (16) by taking the inner product with U a Since W 0 p , U a v = 0 is true for the initial value we can assume that W n p , U a v = 0 and thus However, in the course of the time integration round-off errors can accumulate and it can thus still happen that the W p fail to be exactly orthogonal to the U a . This can cause a small linear drift in the error of mass and momentum. To avoid this we perform an orthonormalization of the W p at the end of each step. This does not negatively impact the accuracy of the numerical method and the incurred computational cost is negligible.

Unconventional integrator
The integrator described in the previous section is not robust if the matrix S has small singular values. The reason for this is that inverting the matrix S is then numerically ill-conditioned. Small singular values commonly occur if the rank of the solution is lower than the rank chosen to conduct the numerical simulation. We refer to [17] for a more detailed discussion, but note that robustness is generally considered a desirable property especially if the rank is adaptively changed during the simulation. For the classic dynamical low-rank approximation a projector splitting integrator has been proposed by Lubich & Oseledets [24] that remedies this deficiency. The method has also been extended to a variety of tensor formats [22,23,11,5,4]. The main utility of the projector splitting integrator is that by using a QR decomposition it avoids the inversion of S. Unfortunately, this projector splitting integrator can not be used in the present situation, because if we solve equation (11) to obtain ip S iq S ip W p , while holding the X i constant, it is not possible to use a QR decomposition to extract the low-rank factors S and W . Instead of the projector splitting integrator, we consider an approach based on the the recently developed unconventional integrator by Ceruti & Lubich [3]. The main idea is that the dynamical low-rank approximation has two parts. On the one hand, the low-rank factors X i and V j determine the subspaces X = span{X i } and V = span{V j } in which an approximation is sought. However, it does not matter how exactly the X i are chosen as long as they span the same approximation space X. On the other hand, the low-rank factor S contains the coefficients that combine the basis functions in an appropriate way in order to obtain a good approximation. Inverting S is only required in equations (10) and (11). That is, it is only required to obtain the low-rank factors X i and V j . We thus proceed as follows. First, we discretize equation (13) to compute and discretize equation (11) to compute The approximation spaces V n+1 and X n+1 are uniquely defined by L n+1 q and K n+1 k , respectively. Thus, we perform a QR decomposition The R parts of the QR decomposition, i.e. R 2 pq and R 1 ik , are simply discarded. We then determine S n+1 as follows This procedure is a discretization of equation (12), where we have used the fact that That is, we transform the coefficient matrix S n ij to the new basis spanned by the previously determined V n+1 j and X n+1 i and then solve equation (12) with the basis held fixed over one time step. This gives us a robust dynamical low-rank integrator as the matrix S needs not be inverted. The downside of this integrator, however, is that mass and momentum conservation up to machine precision is lost. The reason for this is that project the subspaces X n and V n onto the subspaces X n+1 and V n+1 , respectively. Since this projection is not exact, conservation can be lost in the process. We will also see this in the numerical experiments conducted in section 7.

Numerical experiments
In this section we illustrate the conservative dynamical low-rank scheme using a number of numerical examples. The numerical algorithm outlined in sections 4 and 6 will be employed. In order to compute the derivatives in space we employ techniques based on the fast Fourier transform (FFT).
The results of a numerical simulation with the conservative Euler dynamical low-rank integrator and rank r = 6 is shown in Figure 1. For this problem the linearized decay rate can be determined analytically and matches the observed numerical results well. To study the conservation properties we run simulations using m = 0 (no conservation), m = 1 (U 1 ∝ 1 and thus mass conservation), m = 2 (U 1 ∝ 1, U 2 ∝ v and thus mass and momentum conservation), and m = 3 (U 1 ∝ 1, U 2 ∝ v, U 3 ∝ v 2 − 1 and thus mass, and momentum conservation as well as energy conservation in absence of any time integration error). The numerical results in Figure 1 show clearly that mass and momentum are conserved up to machine precision. Thus, the results match perfectly with the theory in sections 5 and 6. Energy is interesting as mandating only mass or momentum conservation can increase the error in energy. This is not entirely unexpected as similar behavior can be observed even for the case of Hamiltonian ordinary differential equations, see e.g. [Example 4.3 in Chap. IV. 4 15]. However, one should note that for m = 3 the error in energy is significantly reduced compared to all the other configurations. Thus, using the energy conservative dynamical low-rank integrator is clearly beneficial, even if the time integration error is taken into account. Reducing the time step size further also improves energy conservation, indicating that the energy error in the m = 3 configuration is limited by the time integration error, as expected. We also note that for the m = 3 configuration the fidelity of the simulation is somewhat worse. The reason for this that the dynamical low-rank integrator is only able to choose 3 basis functions (all others are already determined by imposing conservation) and has thus less freedom to decrease the error in the particle density function (see also the discussion in section 7.3).

Maxwellian with position-dependent velocity
As our second example, we consider a Maxwellian particle density with a position-dependent velocity u(x) = α cos(x) on the domain (x, v) ∈ [0, 4π] × [−6, 6]. The density is given by ρ = 1 + cos( 1 2 x). In the simulation the parameters are chosen as α = 0.2 and = 10 −2 . Note that the initial value given in (19) is not low-rank. However, for small to moderate u, as is the case here, we can use the following expansion (see, e.g., [6,9]) Thus, for the initial value of the dynamical low-rank integrator we use the following rank 6 approximation The results of a numerical simulation with the conservative Euler dynamical low-rank integrator and rank r = 6 is shown in Figure 2. We again observe excellent agreement between theory and the numerical results.
In particular, mass and momentum are conserved up to machine precision (and are between 9 and 14 orders of magnitude smaller than for the m = 0 configuration) and the error in energy is reduced and shown to be dominated by the time integration error.  (19) is shown. To determine the error in the electric energy a reference solution is computed using a full grid semi-Lagrangian scheme that uses 1024 degrees of freedom in both the x and v direction. For mass and energy the relative error and for momentum the absolute error is reported. The simulation is conducted using the conservative Euler dynamical low-rank integrator with step size τ = 10 −4 (except if otherwise indicated), rank r = 6, and 128 grid points in both the spatial and velocity direction are used. The configurations considered are m = 0, m = 1 ( . The time evolution of the electric energy is shown with a dashed black line (axis on the right).

Landau damping using the unconventional dynamical low-rank integrator
The conservative Euler dynamical low-rank integrator can not be used if the rank with which the simulation is run is larger than the rank of the solution (see the discussion in section 6.2). This is often the case for a number of commonly considered plasma instabilities, where the initial value has rank 1. However, for the unconventional dynamical low-rank integrator (see section 6.2) this is not an issue. To demonstrate this we will consider the following linear Landau damping problem which is rank 1. For the simulation we use α = 10 −2 on the domain (x, v) ∈ [0, 4π] × [− 6,6]. The results of the numerical simulation with rank r = 10 is shown in Figure 3. We observe that the analytic decay rate is reproduced accurately by all configurations. In particular, no reduction in fidelity is observed for the m = 3 configuration (as is the case in section 7.1). Thus, having a certain number of basis functions that the algorithm can choose freely is clearly beneficial. As explained in section 6.2, the unconventional integrator is not conservative up to machine precision. The error in mass, momentum, and energy is dominated by the time integration error and thus reducing the time step size reduces the error in the conserved quantities, as is illustrated for mass in Figure 3.   (21) is shown. For mass and energy the relative error and for momentum the absolute error is reported. The simulation is conducted using the unconventional dynamical low-rank integrator with step size τ = 10 −3 , rank r = 10, and 128 grid points in both the spatial and velocity direction are used. The configurations considered are m = 0, m = 1 (U 1 ∝ 1), m = 2 (U 1 ∝ 1, U 2 ∝ v), and m = 3 (U 1 ∝ 1, U 2 ∝ v, U 3 ∝ v 2 − 1). The results for all configurations nearly overlap in the plot.  7,7]. This initial value has rank 1. However, the rank of the solution increases significantly as the system is evolved in time and nonlinear effects become stronger. Nevertheless, it is known that the dynamical low-rank integrator resolves the corresponding dynamics well, at least, up to saturation [11,12]. From the numerical results in Figure 4 we see that this is also the case for the conservative dynamical low-rank integrator. As in the previous section the error in mass, momentum, and energy is dominated by the time integration error in all cases.  (22) is shown. For mass and energy the relative error and for momentum the absolute error is reported. The simulation is conducted using the unconventional dynamical low-rank integrator with step size τ = 10 −3 , rank r = 10, and 128 grid points in both the spatial and velocity direction are used. The configurations considered are m = 0, m = 1 (U 1 ∝ 1), m = 2 (U 1 ∝ 1, U 2 ∝ v), and m = 3 (U 1 ∝ 1, U 2 ∝ v, U 3 ∝ v 2 − 1). The results for all configurations nearly overlap in the plot. The full grid solution is computed using a semi-Lagrangian scheme that uses 1024 degrees of freedom in both the x and v direction. Note that initially the error in mass is dominated by round-off and thus smaller time step sizes lead to a larger error.