Dynamical low-rank approximation of the Vlasov–Poisson equation with piecewise linear spatial boundary

We consider dynamical low-rank approximation (DLRA) for the numerical simulation of Vlasov–Poisson equations based on separation of space and velocity variables, as proposed in several recent works. The standard approach for the time integration in the DLRA model uses a splitting of the tangent space projector for the low-rank manifold according to the separated variables. It can also be modified to allow for rank-adaptivity. A less studied aspect is the incorporation of boundary conditions in the DLRA model. We propose a variational formulation of the projector splitting which allows to handle inflow boundary conditions on spatial domains with piecewise linear boundary. Numerical experiments demonstrate the principle feasibility of this approach.


Introduction
In this work we consider the dynamical low-rank approximation of the Vlasov-Poisson equation in d ≤ 3 spatial dimensions: with a bounded domain Ω x ⊂ R d , and Ω v = R d .This equation models the time evolution of an electron density f = f (t, x, v) of a collisionless plasma as a function of space and velocity in the presence of an electrical field E. We assume the initial condition and inflow boundary conditions of the form Here n x denote the outward normal vectors of the spatial domain Ω x .The electrical field E can either be fixed or dependent on the density f via a Poisson equation: supplemented with appropriate boundary conditions.Here ρ b is a background charge.Simulations of such systems are computationally demanding since the time evolution of an 2d-dimensional, i.e. up to six-dimensional, function has to be calculated.Applying standard discretization schemes thus leads to an evolution equation in O(n 2d ) degrees of freedom, where n is the number of grid points in one dimension and a corresponding high computational effort.To tackle the problem, methods such as particle methods [24], adaptive multiscale methods [5], and sparse grids [17] have been used.
In the seminal paper [7] dynamical low-rank approximation (DLRA) has been proposed for solving (1).DLRA is a general concept of approximating the evolution of time-dependent multivariate functions using a low-rank model, typically based on separation of variables.It originates in classic areas of mathematical physics, such as molecular dynamics [21], and has been proposed in the works [18,19] as a general numerical tool for the time-integration of ODEs on fixed-rank matrix or tensor manifolds.An overview from the perspective of numerical analysis and further references can be found in [23].The approach can also be made rigorous for low-rank functions in L 2 -spaces thanks to their tensor product Hilbert space structure; see, e.g., [2].
For the simulation of (1) using DLRA, a rather natural separation of space and velocity variables is applied.After discretizing the problem in a corresponding tensor product discretization space, one then seeks an approximate solution curve for equation (1) of the form that is, f r (t, •, •) is a rank-r function for every time point t.For the numerical time-integration of the low-rank factors in this representation, [7] adopted the so-called projector-splitting approach from [22], which is one of the work-horse algorithms for DLRA, to the case of the transport equation (1).
In the present work, we wish to study in more detail how to incorporate the inflow boundary condition (2) into such a scheme.In [7] this question was somewhat circumvented by considering periodic boundary conditions, which, however, is not always applicable in real world problems.Enforcing boundary conditions directly on low-rank representations like (4) appears to be rather impractical, or at least poses some difficulties.For the nonlinear Boltzmann equation such an approach has been considered in [16].Here, we will instead start out from a weak formulation of the transport problem (1) including a weak formulation of the boundary condition according to [11,Sec. 76.3].We then use the projector splitting approach for the time-stepping in this weak formulation, where the problem is iteratively reduced to subspaces belonging to variations of space or velocity variables only.We note that in [15] boundary conditions for a two-dimensional product domain (corresponding to d = 1 in our notation) were included in this way through discretization by a discontinuous Galerkin method.
It is important to note already here that in order to still benefit from the seperation of variables, our approach assumes the boundary function g to be a finite sum of tensor products, or at least be approximable by such, see (7) below.Moreover, we require the spatial domain to have a piecewise linear boundary so that the outer normal vectors are piecewise constant, see (6).
As will be derived in section 2, the resulting effective equations for the low-rank factors in our weak formulation of the projector splitting scheme take the form of Friedrichs' systems, that is, systems of hyperbolic equations in weak formulation, that respect the boundary conditions without violating the tensor product structure.These systems can hence be solved by established PDE solvers.We remark that our derivation of these equations remains more or less formal, and the existence of weak solutions (in the continuous setting) will not be studied.Let us also mention the work [20] which is related to our work insofar as it shows that for hyperbolic problems a continuous formulation of the projector splitting integrator prior to discretization has favorable numerical stability properties.
We then proceed by deriving corresponding discrete equations that allow to solve the system numerically.Specifically, in our implementation we will use stabilized finite elements [11] for discretization.We then apply the method to solve the weak Landau damping on a periodic domain (in order to verify the discretization) as well as a linear equation involving a piecewise linear boundary.Here, in addition to the standard projector-splitting approach we will also consider a weak formulation of the unconventional low-rank integrator proposed in [4].It consists in modifications regarding the update strategy for the low-rank factors and allows to make the whole scheme rank-adaptive [3], which is important in practice, based on subspace augmentation.Similar strategies are considered in [14] and [15].
Several recent works have been focussing on the conservation properties in dynamical low-rank integrators [8,9,6,10,13,12], such as mass, momentum and energy, based on modified Galerkin conditions.This important aspect is not yet addressed in the present paper and will remain for future work.
The paper is organized as follows.In section 2 we detail our idea of a weak formulation of the projector-splitting integrator for the Vlasov-Poisson equation which is capable of handling the inflow boundary conditions.From this, in section 3 we obtain discrete equations by restricting to finite-dimensional spaces according to the Galerkin principle, resulting in Algorithm 1.The rank-adaptive unconventional integrator is also considered (Algorithm 2).In section 4 results of numerical experiments are presented to illustrate the principle feasibility of our algorithms.

Weak formulation of the projector splitting integrator
Our goal is to develop a projector splitting scheme for DLRA of the transport problem (1) that is capable of handling the inflow boundary conditions (2).To achieve this, we start out from a weak formulation of (1) including a weak formulation of the boundary condition according to [11,Sec. 76.3].Assume for now that the electrical field is fixed.Let W be an appropriate closed subspace of L 2 (Ω) (e.g.W ⊆ H 1 (Ω) is sufficient) and denote by P W the L 2 -orthogonal projector onto W . Then the goal in the weak formulation is to find f ∈ C 1 ([0, T ], W ) with f (0) = P W f 0 and such that at every t ∈ (0, T ) it holds where Here we have taken into account that Ω v = R d is unbounded and hence the normal vectors for Ω read n = (n x , 0).While our derivations are based on the above formulation with fixed electric field E, the case where E depends on f will be discussed in section 2.5.The variational formulation (5) has the advantage that the boundary conditions are integrated in the bilinear form and right hand side.It therefore can be easily combined with the so-called time-dependent variational principle, also called Dirac-Frenkel principle, underlying DLRA.In order to benefit from the separation of variables in DLRA, it will be necessary to impose some restrictions on the boundary and the function g in the inflow condition (2).Specifically, we assume that the spatial domain has a piecewise linear boundary, so that one has a decomposition where each part Γ (ν) x has a constant outer normal vector n (ν) x .In addition, we require that g itself admits a separation into space and velocity variables by a finite sum of tensor products, or can at least be approximated in this form.
In the following sections we derive the projector splitting approach of DLRA in the weak formulation and derive effective equations for the low-rank factors in (4) in the form of Friedrichs' systems (systems of hyperbolic equations in weak formulation) that are amenable to established PDE solvers.In this way the boundary conditions will be incorporated without sacrificing the tensor product structure.

DLRA and projector splitting
We first present the basic idea of the projector splitting approach for the dynamical low-rank solution of the weak formulation (5).As required for DLRA, we consider a possibly infinitedimensional tensor product subspace where we assume that W x ⊆ L 2 (Ω x ) and W v ⊆ L 2 (Ω v ) are suitable subspaces to ensure that ( 5) is well-defined on W .For example, W x and W v could be subspaces of H 1 , which would correspond to a mixed regularity with respect to space and velocity.
The corresponding manifold M r of low-rank functions in W is with the additional conditions that the X 1 , . . ., X r and V 1 , . . ., V r are orthonormal systems in W x and W v , respectively, and the matrix S = [S ij ] has rank r.Following the Dirac-Frenkel variational principle, a dynamical low-rank approximation for the weak formulation (5) asks for a function for all t ∈ (0, T ), where T fr(t) M r is the tangent space of M r at f r (t) specified below.For the initial value one may take f r (0) = f0 , where f0 is a rank-r approximation of f 0 .The DLRA formulation (9) can be interpreted as a projection of the time derivative of the solution onto the tangent space, which implicitly restricts it to the manifold.The projector splitting integrator from [22] is now based on the rather peculiar fact that the tangent spaces of M r can be decomposed into two smaller subspaces corresponding to variations in the X and V component only.Specifically, let be in M r at time t, then it holds that where Note that these two subspaces are not complementary as they intersect in the space Correspondingly, the L 2 -orthogonal projection onto the full tangent space T fr(t) M r can be decomposed as where P V (t) , P X(t),V (t) , and P X(t) denote the L 2 -orthogonal projections onto the subspaces T X(t) , T X(t),V (t) , and T V (t) , respectively.
The projector-splitting integrator performs the time integration of the system according to the decomposition (10) of the tangent space projector.Sticking to the weak formulation (9), one time step from some point t 0 with to a point t 1 = t 0 + ∆t then consists of the following three steps: 1. Solve the system (9) restricted to the subspace Find an orthonormal system X 1 1 , . . ., X 1 r for K 1 , . . ., K r such that 2. Solve the system (9) restricted to the subspace T X 1 ,V 0 on the time interval [t 0 , t 1 ] with initial condition f , and taking into account the minus sign in the projector.At time This step is often interpreted as a backward in time integration step, which is possible if there is no explicit time dependence in the coefficients or boundary conditions.However, in our case the inflow g is in general time dependent which then does not admit for such an interpretation in general.
3. Solve the system (9) restricted to the subspace T X 1 on the time interval [t 0 , t 1 ] with initial condition f .At time t 1 one obtains As the final solution at time point t 1 one then takes A modification of this scheme is the so called unconventional integrator proposed in [4].There, the K-step and L-step are performed independently using the same initial data from f r (t 0 ) (i.e. the K-step is identical with the one above, but the L-step is performed on T X 0 with initial value f r (t 0 )).As a result, one obtains two new orthonormal sets of component functions X 1  1 , . . ., X 1 r and V 1 1 , . . ., V 1 r for the spatial and velocity domains.The S-step is then performed afterwards in a "forward" way (i.e.without flipping signs) using these new bases.Compared to the standard scheme, this decoupling the S-step from the K-and L-steps also offers a somewhat more natural way of making the scheme rank-adaptive, which is important in practice since a good guess of r might not be known.Specifically, in [3] it is proposed to first augment the new bases to X = {X 0 1 , . . ., X 0 r , K 1 1 , . . ., K 1 r } and V = {V 0 1 , . . ., V 0 r , L 1 1 , . . ., L 1 r } (i.e.including the old bases), orthonormalize these augmented bases, and perform the S-step in the space Note that this space contains functions of rank 2r (in general) and still contains f r (t 0 ) which one takes as initial condition for the S-step.The solution at time point t 1 can then be truncated back to rank r or any other rank (less or equal to 2r) depending on a chosen truncation threshold using singular value decomposition.
Regardless of whether the projector splitting or the modified unconventional integrator is chosen, the efficient realization of the above steps is based on their formulation in terms of the component functions K, S and L only, which in turn requires a certain separability of space and velocity variables in the system (9).We therefore next investigate the single steps in detail, focusing only on the projector splitting integrator and omitting the required modifications for the (rank-adaptive) unconventional integrator (it will be discussed in the discrete setting in section 3.3).A particular focus is how to ensure the required separability for the boundary value terms, which will lead to the assumptions ( 7) and ( 6) that g is separable and the boundary is piecewise linear.As we will see, this gives effective equations in the form of Friedrichs' systems for the factors K and L, and a matrix ODE for S. Their discrete counter-parts as well as the algorithmic schemes for both the projector splitting integrator and the unconventional integrator will be presented in section 3.

Weak formulation of the K-step
The time dependent solution of the first step, i.e. ( 9) restricted to T V 0 , is a function on the time interval [t 0 , t 1 ] with initial condition f (t 0 , x, v) = f r (t 0 , x, v) from (11).Here the V 0 j are fixed and the functions K 1 , . . ., K r need to be determined.Their dynamics are governed by a system of first order partial differential equations which we derive in the following.The test functions w ∈ T V 0 have the form for all j = 1, . . ., r. Testing with w in (9) specifically means We will write this equation in vector form and therefore define which can be regarded as vector-valued functions in (W x ) r .
We investigate the terms in (12) separately.The term regarding the inner product in L 2 (Ω, R) has three parts.The first involves the time derivative: where we have used the pairwise orthonormality of V 1 , . . ., V r .The second part reads with the symmetric r × r matrix .
Finally, the part involving the electrical field can be written as where For the terms in (12) involving the inflow boundary Γ − we make use of our assumption that the spatial boundary ∂Ω x can be decomposed as in ( 6) with a constant normal vector n x of the boundary.We decompose Γ − accordingly into Then the boundary term from the bilinear form can be written as where with characteristic function χ.
For the boundary term on the right hand side we recall the decomposition (7) of the function g into a sum of tensor products.Proceeding as above we obtain .
In summary ( 12) is equivalent to the following system of first order partial differential equations for K(t) in weak formulation with boundary penalty, which needs to be solved for K(t 1 ) with the initial values from the previous approximation (11).

Formulation of the S-step
The solution of the second step of the splitting integrator is a time-dependent function in the subspace space T X 1 ,V 0 , where the evolution for S on the interval [t 0 , t 1 ] is governed by (9) restricted to test functions w ∈ T X 1 ,V 0 , but taking into account the minus sign of P X 1 ,V 0 in the corresponding projector splitting (10).The initial condition reads f (t . The test function in w ∈ T X 1 ,V 0 will be written as For such test functions, we again consider the different contributions as in (12) (with f instead of f ).For the time derivative we get where (•, •) F is the Frobenius inner product on R r×r .The first two contributions of the bilinear form read The terms stemming from the boundary condition take the form Putting everything together and testing with all Σ = [Σ mn ] ∈ R r×r yields the ODE (in strong form) for obtaining S(t 1 ) = S from the initial condition S(t 0 ) = Ŝ.Note again that compared to the original problem (9) all signs (except the one of Ṡ) have been flipped due to the negative sign of the corresponding projector in the splitting (10).

Weak formulation of the L-step
The equations for the L-step are derived in a similar way as the K-step in section 2.2.Here we seek f (t, in the subspace T X 1 with fixed X 1 1 , . . ., X 1 r from the initial value f (t 0 , x, v) = f (t 1 , x, v).The test functions in T X 1 are of the form ..,r .be vector valued functions in (W v ) r .
When restricting the weak formulation as in (12) (with f instead of f ) to test functions in T X 1 , the first three resulting contributions take a similar form as in the K-step, but this time with (note that in the K-step A k x featured the velocity and K x the electric field).The boundary term in the bilinear form gives Hence, the boundary condition on the spatial domain Ω x leads to the additional multiplicative term on the whole domain Ω v .Boundary terms for the velocity variable are not present, since since we assume Ω v = R d to be unbounded.For the right hand side, we have .
The resulting equation for the L-step reads: Find L(t) ∈ (W v ) r for t ∈ [t 0 , t 1 ] such that The initial condition is , where the V 0 j are from f r (t 0 ) in the previous time step (11), and S has been determined in the S-step.

Electrical field
So far we have assumed that the electrical field does not depend on the density f .However, in the setting of a Lie splitting the other case can be included quite easily as well.At the beginning of the time step from t 0 to t 1 , as discussed in sections 2.2-2.4,the electrical field can be computed via the Poisson equation ( 3) and then simply held constant during that step.This will only introduce an error of first order in time, as does the projector splitting itself.

Discrete equations
In sections 2.2-2.4 the governing equations for the projector splitting integrator have been formulated in a continuous setting.They consist of the two Friedrichs' systems ( 14) and ( 17) for K and L, and the finite-dimensional system of ODEs (15) for S. In this section we will introduce a finite element discretization for the approximate numerical solution of the hyperbolic systems and derive the governing discrete equations that need to be solved in the practical computation.
A natural choice for the discretization is to use the discontinuous Galerkin method.However, the resulting equations for K, L, and S (Sec.2) include derivatives of the bases functions X i and V j , see for example Eqns.( 16) and (13).Hence, for a first approach we will use continuous finite element method, which of course has to be stabilized.It remains to future work to investigate the use of DG methods.

Discretization
Instead of an unbounded domain we will now use a finite domain Ω v = [−v max , v max ] d for the velocity variable with v max big enough, and impose periodic boundary conditions.
Let meshes T x/v on the domains Ω x/v be given and define H 1 -conforming discretization spaces

We are now looking for an approximate solution of the Vlasov-Poisson equation (1) in the manifold
of rank-r finite element functions, where as before the systems X h 1 , . . ., X h r and V h 1 , . . ., V h r are assumed to be orthonormal, and S = [S ij ] has rank r.To represent elements in M h r , we introduce coefficient matrices The goal is then to find the solution in the form at prescribed time points t.Note that in order to have the X h 1 , . . ., X h k orthonormal in L 2 as was always assumed above, the matrix X needs to be orthogonal with respect to the inner product of the mass matrix of the basis functions φ (x) α (the matrix M x below).A similar remark applies to V. Following the Galerkin principle, the steps of the projector splitting integrator can be formulated in terms of the matrices X, S, and V by solving the equations ( 14), (15), and ( 17) with the finite element spaces W h x/v instead of W x/v .This means that we use the finite dimensional spaces (W h x/v ) r both as ansatz and test spaces.Specifically, for the Friedrichs' systems ( 14) and ( 17) our aim is to find factors K h j and L h i of the form for i, j = 1, . . ., r.We gather the coefficients in the matrices The governing discrete equations are obtained by inserting the expressions (18) for X i and V j as well as (19) for K i and L j into ( 14), (15), and (17).In order to formulate them conveniently we will use the following additional discretization matrices: In addition, we will take into account that finite element solutions of hyperbolic partial differential equations have to be stabilized.For that purpose we use the continuous interior penalty (CIP) stabilization [11].For a mesh T , its bilinear form is given by where F h 0 is the set of all interior interfaces of T and [[•]] is the jump term with respect to the interface F .The corresponding discretization matrices are denoted by Eventually, the discrete versions of equations ( 14), (15), and ( 17) read as follows: Here we use the abbreviations for discretization matrices M. The parameter δ ≥ 0 controls the stabilization.

Discrete projector splitting integrator
Based on the above equations, the complete realization of the projector splitting scheme, including the appropriate initial conditions and the orthogonalization steps, is outlined in Algorithm 1.It also includes in lines 1 and 2 the computation of the electric field as described in section 2.5.For that purpose the density ρ has to be calculated given by f h r (t 0 , •, •) ∈ M h r by integrating out the v variable.The result is a finite element function ρ h ∈ W h x which is used on the right hand side of the Poisson equation (3).We used ansatz functions of order two to solve the resulting linear equation directly for the potential Φ h .The electrical field E h can be computed as the negative gradient of Φ h and is a discontinuous finite element function on the same mesh.Now, updating the discretization matrices M x,E k involving the electrical field, the projector splitting steps from t 0 and t 1 can be performed.
Algorithm 1 First order dynamical low rank integrator with projector splitting Data: X 0 , S 0 , V 0 ▷ update electrical field 1: Solve Poisson equation (3) for density resulting from X 0 , S 0 , V 0 2: Compute E and update

Rank adaptive algorithm
The rank-adaptive algorithm is based on the unconventional low-rank integrator proposed in [4] and its rank-adaptive extension in [3].It is displayed in Algorithm 2. Starting with a rank-r 0 function, it first performs independent K-steps and L-steps from the initial data X 0 and V 0 .The computed factors K 1 and L 1 are used to enlarge the previous bases for X 0 and V 0 to dimension 2r 0 each.This is achieved by computing orthonormal bases X and V of the augmented matrices which differs from (20b) in the signs of the right-hand side.For the initial condition S(t 0 ) = Ŝ, X 0 S 0 V 0 has to be expressed with respect to the larger bases X and V, i.e. as X Ŝ V.The new solution then has rank 2r 0 (in general) and can be truncated to a lower rank r 1 according to a tolerance, which yield the actual new factors Algorithm 2 First order rank adaptive unconventional integrator (RAUC) ▷ Compute augmented V basis 6: L 0 = V 0 (S 0 ) T 7: Solve (20c) on interval [t 0 , t 1 ] with initial condition L 0 to obtain L 1 8: Compute orthonormal basis V with respect to and monotonically decreasing σ i 12: For tolerance ϵ compute r 1 ≤ 2r 0 as the minimal number such that

Numerical experiments
We present numerical results of two experiments for testing the methods described in section 3. Since to our knowledge there does not exist an analytic solution of the Vlasov-Poisson equation on bounded domains to which our numerical results could be compared, we chose to consider the following two setups.The first experiment will be the classical Landau damping where analytic results on the decay of the electric energy are known.However, the domain is assumed to be periodic so that no boundary conditions are needed.This scenario has been treated in several previous works.Our second experiment will then include boundary conditions on a polygonal spatial domain to actually test our proposed approach for handling these.
In our tests linear finite elements are used and the matrix ODEs (lines 4, 6, and 8 in Alg. 1, and line 11 in Alg. 2) are solved using an explicit Runge-Kutta method of third order.The implementation is based on the finite element library MFEM [1] and uses its Python wrapper PyMFEM. 1 The computations are carried out using standard numerical routines from numpy and scipy.All experiments have been performed on a desktop computer (i9 7900, 128 GB RAM), the code is available online. 2

Landau damping
As a first numerical test we consider the classical Landau damping in two spatial and velocity dimensions.We use periodic domains Ω x = [0, 4π] 2 and Ω v = [−6, 6] 2 and a background density of ρ b = 1.As initial condition we choose The same setup was investigated in [7].For this case linear analysis shows that the electric field decays with a rate of γ ≈ 0.153.For the discretization we use regular grids with n x = 64 2 and n v = 256 2 degrees of freedom (level 0).Using Algorithm 1 the simulation was carried out fixed ranks r = 5, 10, 15 and a time step of ∆t = 0.005.The results are shown in Figure 1.As can be seen in the top left plot, for rank r = 5 the computed electric energy 1  2 Ωx |E(t, x)| 2 dx exhibits the analytical decay rate approximately up to t = 35.For rank r = 10 the electrical energy starts to deviate at the end of the time interval, while the solution with rank r = 15 shows the correct rate within the full simulation time.
Furthermore, it is known that the particle number, the total energy and the entropy, that is, the quantities are invariants of the exact solution.In Figure 1 we see that the mass and the total energy are almost conserved, whereas the entropy is only preserved up to a small error.
For the rank adaptive case (Algorithm 2) the system was simulated for the same discretization (level 0) and different tolerances ϵ.The corresponding results in Figure 2 show that the electric energy deteriorates at around t = 25 for both tolerances ϵ = 10 −5 , 10 −6 .Although the rank increases up to 40 (the maximal rank allowed) for the case of ϵ = 10 −6 the accuracy in the electric energy does not improve.To improve accuracy the simulation is carried out on a uniformly refined spatial and velocity mesh (level 1).The results for a time step of ∆t = 0.00125 and a tolerance of ϵ = 10 −6 are shown in Fig. 2. The electric energy shows the correct decay in electric energy up to approximately t = 40.
Investigating the computed bases X and V in more detail shows that in the rank adaptive case spurious oscillations are present.In Figure 3 two exemplary basis functions X i (x) of level 0 at time t = 50 are displayed.In the fixed rank case (r = 15) the basis function is much smoother than in the rank adaptive simulation (ϵ = 10 −5 ).
In summary spurious modes may enter in the course of the simulation for the rank adaptive simulation.However, the accuracy can be improved by refining the discretization.In contrast, the algorithm using a fixed rank seems to have a regularizing effect.It remains to future work to investigate this effect more closely.

Inflow boundary condition with constant electrical field
In the second example we focus on the boundary condition and compare the numerical to an analytical solution.In this setting, we are solving the transport equation ( 1) in 2 + 2 dimensions with a constant electrical field E = [0, 4] T on a triangular spatial domain on the boundary of Ω x will be determined by a function f , which is a solution of the same equation as for f , but on the whole domain R 2 and with initial condition f (0, x, v) = f0 (x, v).
Here, f0 has compact spatial support just outside the triangular domain Ω x and a compactly supported velocity distribution centred around v = [2, 0] T .More precisely, we set where σ x = 0.2, σ v = 0.5, and 1] and centered around 0. In consequence f0 is a product function supported in a four-dimensional cube with side lengths controlled by σ x and σ v .By the method of characteristics we obtain Restricting f on Ω x × R 2 solves the original problem and can be used to assess the quality of its numerical solution f h r .In order to use f as the inflow function in our dynamical low-rank approach, we have to work with an approximation in separable form instead.For that purpose, we use the fact that, by (25) and our choice (24), f can be written as Figure 4 shows the evolution of the computed spatial density for level ℓ = 1 at different time steps together with its error computed using the analytical solution f .At the beginning of the simulation, the density flows into the domain, is transported and finally leaves the domain.The error remains reasonably small during the whole process.
A more detailed investigation of the L 2 error is depicted in the upper part of obtained by comparing the numerical solutions f h r to the analytical solution f on a once uniformly refined grid.As one can see, in the beginning the error increases for all levels but decays again for t ≥ 0.35.At about that time a significant fraction of the density has left the domain already (see Figure 4).The maximal error decays as the level increases.
The lower part of the figure shows the ranks which were used by the rank adaptive algorithm.As more particles enter the domain and the distribution spreads out, the ranks increase.For higher levels of discretization a smaller truncation parameter is used in order to balance the low-rank approximation error with the discretization error, leading to higher ranks.

Conclusion and outlook
In this paper we have studied how to incorporate inflow boundary conditions in the dynamical lowrank approximation (DLRA) for the Vlasov-Poisson equation based on its weak formulation.The single steps in the projector splitting integrator, or the rank-adaptive unconventional integrator, can be interpreted as restrictions of the weak formulation to certain subspaces of the tangent space.The efficient solution of these sub-steps requires the separability of the boundary integrals, which is ensured for piecewise linear boundaries together with a separable inflow function.The resulting equations can be solved using FEM solvers based on Galerkin discretization.We confirmed the feasibility of our approach in numerical experiments.
As a next step, the conservation of physical invariants such as mass and momentum in the numerical schemes could be addressed, perhaps by extending methods from [8,6,10,13,12].Enforcing nonnegativity of the density function in the DLRA approach is another open issue.
A potential advantage of the weak formulation for the sub-problems in the projector splitting integrator is that in principle it should allow for a great flexibility regarding the discretization spaces.In particular, they do not need to be fixed in advance and mesh-adaptive methods could be used for solving the sub-steps, provided suitable interpolation and prolongation operations are available, as well as adaptive solvers for Friedrichs' systems.We leave this as possible future work.Other elaborate adaptive integrators have been proposed in [15], which also could be combined with our approach.

Figure 1 :
Figure1: Simulation results of the 2+2-dimensional Landau damping using fixed ranks (Alg.1); electric energy including the analytical decay rate (upper left) and relative error of the invariants(22).For the total energy and entropy the lines for all three ranks almost overlap.

Figure 2 :
Figure 2: Simulation results of the 2+2-dimensional Landau damping using the rank adaptive algorithm (Algorithm 2) for different tolerances ϵ and discretizations; electric energy including the analytical decay rate (left) and ranks (right).

Figure 5 Figure 5 :
Figure 5: Numerical solution of (1) for different levels of discretization.The upper graph shows the L 2 error computed with respect to the analytical solution f , see (25), on a uniformly refined grid.The lower graph shows the ranks used by the rank adaptive integrator.