A posteriori error estimates for the one and one-half Dimensional Relativistic Vlasov-Maxwell system

This paper concerns a posteriori error analysis for the streamline diffusion (SD) finite element method for the one and one-half dimensional relativistic Vlasov-Maxwell system. The SD scheme yields a weak formulation, that corresponds to an add of extra diffusion to, e.g. the system of equations having hyperbolic nature, and convection-dominated convection diffusion problems. A procedure that improves the convergence of finite elements for this type of problems. The a posteriori error estimates relay on a dual problem formulation and yields an error control based on the, computable, residual of the approximate solution. The lack of dissipativity enforces us considering negative norm estimates. To derive these estimates, the error term is split into an iteration and an approximation error where the iteration procedure is assumed to converge. The computational aspects and implementations, which justify the theoretical results of this part, are the subject of these studies and addressed in [5].


Introduction
This paper concerns a posteriori error analysis for approximate solution of the Vlasov-Maxwell (VM) system by the streamline diffusion (SD) finite element methods. Our main objective is to prove a posteriori error estimates for the SD scheme in the H −1 (H −1 ) and L ∞ (H −1 ) norms for the Maxwell equations and L ∞ (H −1 ) norm for the Vlasov part. The VM system lacking dissipativity exhibits severe stability draw-backs and the usual L 2 (L 2 ) and L ∞ (L 2 ) errors are only bounded by the residual norms. Thus, in order not to rely on the smallness of the residual errors, we employ the negative norm estimates to pick up convergence rates also involving powers of the mesh parameter h and having optimality properties due to the maximal available regularity of the exact solution. Both Vlasov and Maxwell equations are of hyperbolic type and for the exact solution in the Sobolev space H r+1 , the classical finite element method for hyperbolic partial differential equations will have, an optimal, convergence rate of order O(h r ), where h is the mesh size. On the other hand, with the same regularity (H r+1 ) the optimal convergence rate for the elliptic and parabolic problems is of order O(h r+1 ). This phenomenon, and the lack of diffusivity in the hyperbolic equations which cause oscillatory behavior in the finite element schemes, sought for constructing modified finite element schemes that could enhance stability and improve the convergence behavior for hyperbolic problems. In this regard, compared to the classical finite element, the SD schemes, corresponding to the add of diffusion term to the hyperbolic equation, are more stable and have an improved convergence rate viz, O(h r+1/2 ). Roughly, the SD method is based on a weak formulation where a multiple of convection term is added to the test function. With this choice of the test functions the variational formulation resembles to that of an equation which, compared to the original hyperbolic equation, has an additional diffusion term of the order of the multiplier.
A difficulty arises deriving gradient estimates for the dual problems, which are crucial for the error analysis for the discrete models in both equation types in the VM system. This is due to the lack of dissipative terms in the equations. An elaborate discussion on this issue can be found in the classical results, e.g., [17], [9] and [23] as well as in relatively recent studies in [10] and [22].
We use the advantage of low spatial dimension that, assuming sufficient regularity, yields existence and uniqueness through d'Alembert formula. This study can be extended to higher dimensional geometries, where a different analytical approach for the well-posedness is available in the studies by Glassey and Schaeffer in, e.g., [12] and [13]. Numerical implementations for this model will appear in the second part: [5]. We also mention related studies [18] and [19] for the Maxwell's equations where stabilized interior penalty method is used.
Problems of this type have been considered by several authors in various settings. In this regard, theoretical studies for the Vlasov-Maxwell system relevant to our work can be found in, e.g. [8] for treating the global weak solutions, [14] for global weak solutions with boundary conditions and more adequately [11]- [13] for relativistic models in different geometries. SD methods for the hyperbolic partial differential equations have been suggested by Hughes and Brooks in [15]. Mathematical developments can be found in [16]. For SD studies relevant to our approach see, e.g., [2] and the references therein some containing also further studies involving discontinuous Galerkin schemes and their developments.
An outline of this paper is as follows: In the present Section 1, following the introduction, we comment on particular manner of various quantities in the Maxwell equations and introduce the relativistic one and one-half dimensional model with its well-posedness property. In Section 2 we introduce some notations and preliminaries. Sections 3 is devoted to stability bounds and a posteriori error estimates for the Maxwell equations in both H −1 (H −1 ) and L ∞ (H −1 ) norms. Sections 4 is the counterpart of Section 3 for the Vlasov equation which is now performed only in L ∞ (H −1 ) norm.
Finally, in our concluding Section 5, we summarize the results of the paper and discuss some future plans.
Throughout this note C will denote a generic constant, not necessarily the same at each occurrence, and independent of the parameters in the equations, unless otherwise explicitly specified.
The Vlasov-Maxwell (VM) system which describes time evolution of collisionless plasma is formulated as Here f is density, in phase space, of particles with mass m, charge q and velocitŷ Further, the charge and current densities are given by ρ(t, x) = 4π qf dv, and j(t, x) = 4π qfv dv, (1.2) respectively. For a proof of the existence and uniqueness of the solution to VM system one may rely on mimicking the Cauchy problem for the Vlasov equation through using Schauder fixed point theorem: Insert an assumed and given g for f in (1.2). Compute ρ g , j g and insert the results in Maxwell equations to get E g , B g . Then insert, such obtained, E g and B g in the Vlasov equation to get f g via an operator Λ: f g = Λg. A fixed point of Λ is the solution of the Vlasov equation.
For the discretized version employ, instead, the Brouwer fixed point theorem. Both these proofs are rather technical and non-trivial. The fixed points argument, rely on viewing the equations in the Maxwell's system as being valid independent of each others, but the quantities f , B, E, j and ρ are physically related to each others by the Vlasov-Maxwell system of equations and it is not the case that some of them are given to determine the others. However, in one and one-half geometry, relying on d'Alembert formula Schauder/Brouwer fixed point approach, is unnecessary. The fixed point approach, which was first introduced by Ukai and Okabe in [24] for the Vlasov-Poisson system, is performed for the Vlasov-Maxwell system in [21] in full details and therefore is omitted in here.
1.1. Relativistic model in one and one-half dimensional geometry. Our objective is to construct and analyze SD discretization schemes for the relativistic Vlasov-Maxwell model in one and one-half dimensional geometry (x ∈ R, v ∈ R 2 ), which then can be generalized to higher dimensions: The system (1.3) is assigned with the Cauchy data and with This is the only initial data that leads to a finite-energy solution (see [11]). In (1.3) we have for simplicity set all constants equal to one. The background density ρ b (x) is assumed to be smooth, has compact support and is neutralizing. This yields To carry out the discrete analysis, we shall need the following global existence of classical solution due to Glassey and Schaeffer [11]. Theorem 1.1 (Glassey, Schaeffer). Assume that ρ b , the background density, is neutralizing and we have Then, there exists a global C 1 solution for the Relativistic Vlasov-Maxwell system.
Note that for the well-posedness of the discrete solution the existence and uniqueness is due to [21], whereas the stability of the approximation scheme is justified throughout Sections 3 and 4.

Assumptions and notations
Let Ω x ⊂ R and Ω v ⊂ R 2 denote the space and velocity domains, respectively. We shall assume that f (t, x, v) , E 2 (t, x), B(t, x) and ρ b (x) have compact supports in Ω x and that f (t, x, v) has compact support in Ω v . Since we have assumed neutralizing background density, i.e. ρ(0, x)dx = 0, it follows that E 1 also has compact support in Ω x (see [11]). Now we will introduce a finite element structure on Ω Before we define our finite dimensional spaces we need to introduce some function spaces, viz where In the discretization part, for k = 0, 1, 2, . . ., we define the finite element spaces where P k (·) is the set of polynomial with degree at most k on the given set. We shall also use some notation, viz To proceed, we shall need to perform an iterative procedure: starting with f h,0 we compute the fields E h, 1 1 , E h,1 2 and B h,1 and insert them in the Vlasov equation to get the numerical approximation f h,1 . This will then be inserted in the Maxwell equations to get the fields E h, 2 1 , E h,2 2 and B h,2 and so on. The iteration step i yields a Vlasov equation for f h,i with the fields E h,i 1 , E h,i 2 and B h,i . We are going to assume that this iterative procedure converges to the analytic solution of the Vlasov-Maxwell system. More specifically, we have assumed that the iteration procedure generates Cauchy sequences.
Finally, due to the lack of dissipativity, we shall consider negative norm estimates. Below we introduce the general form of the function spaces that will be useful in stability studies and supply us the adequate environment to derive error estimates with higher convergence rates. In this regard: Let Ω be a bounded domain in R N , N ≥ 2. For m ≥ 0 an integer, 1 ≤ p ≤ ∞ and G ⊆ Ω, W m p (G) denotes the usual Sobolev space of functions with distributional derivatives of order ≤ m which are in L p (G). Define the seminorms and the norms We shall only use the L 2 -version of the above norm.

A posteriori error estimates for the Maxwell equations
Our main goal in this section is to find an a posteriori error estimate for the Maxwell equations. Let us first reformulate the relativistic Maxwell system, viz Then, the Maxwell equations can be written in compact (matrix equations) form as The streamline diffusion method on the ith step for the Maxwell equations can now be formulated as: find W h,i ∈Ṽ h such that for m = 1, 2, . . . , M , and δ is a multiple of h (or a multiple of h α for some suitable α), see [10] for motivation of choosing δ. Now we are ready to start the a posteriori error analysis. Let us decompose the error into two parts where W i is the exact solution to the approximated Maxwell equations at the ith iteration step: a posteriori error analysis for the Maxwell equations. We will start by estimating the numerical errorẽ i . To this end, we formulate the dual problem: The idea is to use the dual problem to get an estimate on the H −1 -norm of the errorẽ i . Multiplying (3.4) withẽ i and integrating overQ T we obtain Likewise, due to the fact that all involved functions have compact support in Ω x , we can write (3.7) Inserting (3.6) and (3.7) into the error norm (3.5), we get Let nowφ be an interpolant of ϕ and use (3.3) with g =φ to get Now, to proceed we introduce the residuals Further, we shall use two projections, P and π, for our interpolantsφ. These projections will be constructed from the local projections ∀u ∈Ṽ h m and π m u| Sm = 1 h Im u(t, ·) dt. Now we define P and π, slab-wise, by the formulas respectively. See Brezzi et al. [6] for the details on commuting differential and projection operators in a general setting. Now we may choose the interpolants as ϕ = P πϕ = πP ϕ, and write an error representation formula as (3.10) To estimate J 1 and J 2 we shall use the following identity We estimate each term in the error representation formula separately: where in the last estimate the, piecewise time-constant, residual is moved inside the time integration. As for the J 2 -term we have that (3.12) Thus we can derive the estimate To estimate the second term in (3.9) we proceed in the following way (3.14) Combining (3.9)-(3.14) yields To get an estimate for the H −1 -norm we need to divide both sides by χ H 1 (QT ) and take the supremum over χ ∈ [H 1 (Q T )] 3 . We also need the following stability estimate.
Lemma 3.1. There exists a constant C such that Proof. To estimate the H 1 -norm of ϕ we first write out the equations for the dual problem explicitly:

(3.15)
We start by estimating the L 2 -norm of ϕ. Multiply the first equation by ϕ 1 and integrate over Ω x to get

Standard manipulations yields
The second integral vanishes because ϕ 1 is zero on the boundary of Ω x . We therefore have the following inequality Applying Grönwall's inequality and then integrating over (0, T ) we end up with the stability estimate Similarly we estimate the second and third component of ϕ as follows: We multiply the second and the third equations of (3.15) with ϕ 2 and ϕ 3 , respectively. Adding the resulting equations and integrating over Ω x , yields the equation We may rewrite (3.16) as (3.17) Note that the third term on the left hand side of (3.17) is identically equal to zero because both ϕ 2 and ϕ 3 vanish at the boundary of Ω x . We therefore have the following inequality . Integrating over (t, T ) we get that Applying Grönwall's inequality and then integrating over (0, T ) we end up with the stability estimate Next we need to prove that the L 2 -norms of the derivatives of ϕ are bounded by χ H 1 (QT ) . To do this we first note that ϕ has analytical solutions, see [21], (3.18) Let us start by estimating the x-derivative of ϕ 1 . By the above formula for ϕ 1 we have that

Cauchy-Schwartz inequality and a suitable change of variables yields
Integrating both sides of the inequality over (0, T ), gives the estimate Now we can use this inequality together with the first equation in (3.15) to get an estimate for the time derivative of ϕ 1 : Similar estimates can be derived for the derivatives of ϕ 2 and ϕ 3 . We omit the details and refer to the estimations of the derivatives for ϕ 1 .
Summing up we have proved following estimate for the numerical errorẽ i .

Theorem 3.2 (A posteriori error).
There exists a constant C such that As for the iterative errorẼ i we assume that W i converges to the analytic solution, so that, for sufficiently large i,ẽ i is the dominating part of the error W − W h,i , see [21] for motivation of the iteration assumption. Therefore, for large enough i, we have that This together with Theorem 3.2 yields the following result: There exists a constant C such that 3.2. L ∞ (H −1 ) a posteriori error analysis for the Maxwell equations. In this part we perform a L ∞ (H −1 ) error estimate. The interest in this norm is partially due to the fact that the Vlasov part is studied in the same environment. To proceed we formulate a new dual problem as where χ ∈ [H 1 (Ω x )] 3 . We multiplyẽ i (T, x) by χ and integrate over Ω x to get Using (3.6) and (3.7) the above identity can be written as With similar manipulations as in (3.8), this equation simplifies to Following the proof of Theorem 3.2 we end up with the following result: Theorem 3.4. There exists a constant C such that In the proof of this theorem we use the stability estimate: Lemma 3.5. There exists a constant C such that The proof is similar to that of Lemma 3.1 and therefore is omitted. With the same assumption on the iteration errorẼ i as in (3.19), the numerical errorẽ i will be dominant and we have the following final result: Corollary 3.6. There exists a constant C such that

A posteriori error estimates for the Vlasov equation
The study of the Vlasov part rely on a gradient estimate for the dual solution. Here, the L 2 -norm estimates, would only yield error bounds depending on the size of residuals, with no h α -rates. Despite the smallness of the residual norms this, however, does not imply concrete convergence rate and smaller residual norms require unrealistically finer degree of resolution. The remedy is to employ negative norm estimates, in order to gain convergence rates of the order h α , for some α > 0. In this setting a H −1 (H −1 )-norm is inappropriate. Hence, this section is devoted to L ∞ (H −1 )-norm error estimates for the Vlasov equation in the Vlasov-Maxwell system. 4.1. L ∞ (H −1 ) a posteriori error estimates for the Vlasov equation. The streamline diffusion method on the ith step for the Vlasov equation can be formulated as: find f h,i ∈ V h such that for m = 1, 2, . . . , M , where the drift factor is computed using the solutions of the Maxwell equations. As in the Maxwell part we decompose the error into two parts where f i is the exact solution of the approximated Vlasov equation at the ith iteration step: To estimate the numerical error we formulate a corresponding dual problem as Since G(f h,i−1 ) is divergence free (i.e., we have a gradient field), we may manipulate the sum above as in (3.6) and (3.7), ending up with Adding and subtracting appropriate auxiliary terms, see (3.8), this simplifies to where we have used (4.2). Let nowΨ i be an interpolant of Ψ i and use (4.1) with g =Ψ i to get

(4.4)
In the sequel we shall use the residuals Finally, we introduce the projections P and π defined in a similar way as in the Maxwell part, where the local projections The main result of this section is as follows: There exists a constant C such that Proof. We choose the interpolants so thatΨ The termsJ 1 andJ 2 are estimated in a similar way as J 1 and J 2 , ending up with the estimate . Hence, it remains to bound the second and third terms in (4.4). To proceed, recalling the residuals R i 1 and R i 2 , we have that , where we used that δ = Ch. Summing up we have the estimate . Together with the following stability estimate this completes the proof.
Proof. We start estimating the L 2 -norm: multiply the dual equation (4.3) by Ψ i and integrate over Ω to get The second integral is zero, since Ψ i vanishes at the boundary of Ω, so that we have Integrating over (t, T ) yields Ψ i (t, ·, ·) 2 L2(Ω) ≤ χ 2 L2(Ω) . Once again integrating in time, we end up with It remains to estimate ∇Ψ i L2(QT ) . To this approach we rely on the characteristic representation of the solution for (4.3), see e.g., [20]: with X(s, t, x, v) and V (s, t, x, v) being the solutions to the characteristic system where Hence, we have Thus, it suffices to estimate the gradients of X and V . Below we shall estimate the derivatives of X, V 1 and V 2 with respect to x. The estimates with respect to v i , i = 1, 2 are done in a similar way. Differentiating (4.5) with respect to x we get Integrating these equations over [s, t] and then taking the absolute values give Summing up we have that

Now an application of the Grönwall's lemma yields
By similar estimates for derivatives with respect to velocity components we have The estimates (4.6) and (4.7) would result to the key inequalities ∂ x Ψ ≤ C T ∇χ and ∂ vj Ψ ≤ C T ∇χ , j = 1, 2, (4.8) which together with the equation for Ψ gives Summing up we have shown that which proves the desired result.
Using the assumption that the iteration error E i converges and is dominated by the numerical error e i , together with Theorem 4.1, we get the following result:

Conclusions and future works
We have presented an a posteriori error analysis of the streamline diffusion (SD) scheme for the relativistic one and one-half dimensional Vlasov-Maxwell system. The motivation behind our choice of the method is that the standard finite element method for hyperbolic problems is sub-optimal. The streamline diffusion is performed slab-wise and allows jump discontinuities across the time grid-points. The SD approach have stabilizing effect due to the fact that, adding a multiple of the streaming term to the test function, it corresponds to an automatic add of diffusion to the equation.
Numerical study of the VM system has some draw-backs in both stability and convergence. The VM system lacks dissipativity which, in general, affects the stability. Further, L 2 (L 2 ) a posteriori error bounds would only be of the order of the norms of residuals. In our study, in order to derive error estimates with convergence rates of order h α , for some α > 0, the H −1 (H −1 ) and L ∞ (H −1 ) environments are employed. However, because of the lack of dissipativity, the H −1 (H −1 )-norm is not extended to the Vlasov part, where appropriate stability estimates are not available. Therefore the numerical study of the Vlasov part is restricted to the L ∞ (H −1 ) environment.
The computational aspects and implementations, which justify the theoretical results of this part, are the subject of a forthcoming study which is addressed in [5].
Future studies, in addition to considering higher dimensions and implementations, may contain investigations concerning the assumption on the convergence of the iteration procedure, see end of Section 2.
We also plan to extend this study to Vlasov-Schrödinger-Poisson system, where we rely on the theory developed by Ben Abdallah et al. in [3] and [4] and consider a novel discretization procedure based on the mixed virtual element method, as in Brezzi et al. [7].