On hp-Streamline Diffusion and Nitsche schemes for the Relativistic Vlasov-Maxwell System

We study stability and convergence of $hp$-streamline diffusion (SD) finite element, and Nitsche's schemes for the three dimensional, relativistic (3 spatial dimension and 3 velocities), time dependent Vlasov-Maxwell system and Maxwell's equations, respectively. For the $hp$ scheme for the Vlasov-Maxwell system, assuming that the exact solution is in the Sobolev space $H^{s+1}(\Omega)$, we derive global {\sl a priori} error bound of order ${\mathcal O}(h/p)^{s+1/2}$, where $h (= \max_K h_K)$ is the mesh parameter and $p (= \max_K p_K)$ is the spectral order. This estimate is based on the local version with $h_K=\mbox{ diam } K$ being the diameter of the {\sl phase-space-time} element $K$ and $p_K$ is the spectral order (the degree of approximating finite element polynomial) for $K$. As for the Nitsche's scheme, by a simple calculus of the field equations, first we convert the Maxwell's system to an {\sl elliptic type} equation. Then, combining the Nitsche's method for the spatial discretization with a second order time scheme, we obtain optimal convergence of ${\mathcal O}(h^2+k^2)$, where $h$ is the spatial mesh size and $k$ is the time step. Here, as in the classical literature, the second order time scheme requires higher order regularity assumptions. Numerical justification of the results, in lower dimensions, is presented and is also the subject of a forthcoming computational work [20].


Introduction
We study stability and convergence for some specific finite element schemes for a model problem for the three dimensional, relativistic, Vlasov-Maxwell (VM) system with 3-dimensional spatial domain (x ∈ Ω x ⊂ R 3 ) and 3-dimensional velocities domain (v ∈ Ω v ⊂ R 3 ). The objective is two-fold: i) Numerical investigations of the hp -version of the streamline diffusion (SD) finite element method for VM where both Maxwell's and Vlasov equations are discretized using a space-velocity-time scheme both in h (mesh size) and in p (spatial order) versions. In this part we derive optimal a priori error bounds for a SD scheme in a L 2 -based norm.
ii) The study of the combined effect of Nitsche's symmetrization (cf [21] and [6]) in the spatial scheme for a Galerkin method and a time discretization, for a second order pde obtained through the combined Maxwell's fields.
The SD method was suggested by Hughes and Brooks in [18] for the fluid problems. The method was further developed (by T. Houghs and co-workers) to include several engineering problems. A through mathematical analysis was first given by Johnson et al in [19] in a study of the Navier-Stokes equations and was extended to most pdes with hyperbolic nature, where, e.g., [2]- [5] and [24] are relevant in the present study. In the SD method the test function is modified by adding a multiple of the streaming part in the equation, in terms of the test function, to it. Then, in the weak formulation we obtain a multiple of streaming terms in test and trial functions. This can be viewed as an extra diffusion term in the streaming direction in the original equation. Hence, the name of the method (the streamline diffusion). Such an extra diffusion would improve both the stability and convergence properties of the underlying Galerkin scheme. It is well known that the standard Galerkin method has a weaker convergence property for the hyperbolic problems: O(h s−1 ) versus O(h s ) for the elliptic and parabolic problems with exact solution in the Sobolev space H s (Ω). The SD method improves this weak convergence to O(h s−1/2 ) (see [1] for Sobolev spaces of non-integer order) and also, having an upwinding character, enhances the stability.
These two properties are achieved by discontinuous Galerkin as well (see, e.g. [7]). The hp-approach is to capture local behavior in the sense that: in the vicinity of singularities refined mesh h is combined with the lower order (small p) polynomial approximations, whereas in more smooth regions higher order polynomials (large p) and non-refined (large h) meshes are used. In a sense the hp-approach may be interpreted as a kind of automatic adaptivity.
The Vlasov-Maxwell (VM) system which describes the time evolution of collisionless plasma is formulated as Here f is the density, in phase space, time of particles with charge q, mass m and velocitŷ v = (m 2 + c −2 |v| 2 ) −1/2 v (v is momentum).
Further, c is the speed of light and the charge and current densities ρ and j are given by ρ(t, x) = 4π qf dv and j(t, x) = 4π qfv dv.
The phase-space variables may have different dimension: The Vlasov-Maxwell equations arises in several branches of continuum physics, e.g. in astrophysics or rarefied gas dynamics. The main assumption underlying the model is that collisions are rare and therefore negligible. In this setting the above system describes the motion of a collisionless plasma, e.g., a high-temperature, low-density, ionized gas.
For a thorough mathematical study of VM models we refer to DiPerna and Lions [12] and a most recent work by Glassey and co-workers [14]- [15] and the references therein. The results in [14] are for a lower dimensional model where the interest lies in classical solutions, and are based on compactness and regularity assumptions on the initial density and fields.
The main mathematical concern in dealing with the Vlasov-Maxwell system is related to the nonlinear term (E+v×B)·∇ v f which can be written in the divergence [12] the nonlinear form Numerical approaches for the VM system have been considered by several authors in different setting. The most relevant studies to this work are given by Gamba and co-workers [9] devoted to a discontinuous Galerkin approach, and Standar in [22] where the stability and a priori error estimates for the h version of SD method for VM are derived. As some related studies we mention the analysis of a one dimensional model problem for the relativistic VM system in an interval given by Filbet and co-workers in [13]. Also in a very recent work [11] Degond and co-workers study a particle-in-cell method for the Vlasov-Maxwell system. An outline of this paper is as follows. We gather notation and assumptions in Section 2. In Section 3 we formulate the SD schemes for both Maxwell's equations and the Vlasov-Maxwell system. Section 4 is on stability and convergence of the hp SD finite element method of the Maxwell's equations, based on a space-time iterative scheme. We insert such approximated field function in the drift term in Vlasov-Maxwell equation and prove stability and derive optimal convergence rates in the SD phase-space-time discretization scheme. Section 5 is devoted to a the study of a Nitsche scheme combined with a time discretization for a second order pde obtained from the Maxwell's equations. Here, we rely on a modified form of the Ritz projection and derive optimal error estimates for the Nitsche scheme in spatial discretization. Assuming additional regularities in time we also prove optimal convergence of a second order time scheme for the fields. Finally, in our concluding Section 6 we present some numerical tests of the studied schemes in lower dimensional geometry.

Notation and Assumptions
The divergence equations in (1.1) can be derived from the rest of the equations, assuming that the initial data E 0 and B 0 satisfy corresponding divergence equations. Hence we will consider the following relativistic Vlasov-Maxwell system in R d (in the paper we focus on the dimension d = 3, but one can easily obtain the analogous results for d = 2) withv = (1 + |v| 2 ) −1/2 v and j(t, x) = fv dv, where for simplicity we set the charge q and all constants equal to one.
Our objective is to use an iterative scheme to approximate the solution of the Vlasov-Maxwell (henceforth referred as VM) equations. First we take a guess for the density f and then calculate the corresponding j. Next, we plug these quantities into the Maxwell's equations and solve these equations. Finally, we solve the Vlasov equation with the such approximated E and B as coefficients.
We start from the Maxwell's part.
Then the Maxwell's equations in (2.1) can be written in the following form: where ∂ i denotes the derivative with respect to x i . Hence defining the matrices the Maxwell's equations can be written as the system Now we return to the Vlasov equation given by For simplicity, we introduce the notation and define the total gradient so that, we can rewrite the Vlasov equation in compact form as Note that G is divergence free Throughout this paper 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.

hp-Streamline Diffusion Method
Let Ω x ⊂ R 3 and Ω v ⊂ R 3 denote the space and velocity domains, respectively. We assume that f (t, x, v), E i (t, x), B i (t, x) for i = 1, 2, 3 have compact supports in Ω x and that f (t, x, v) has compact support in Ω v . Now we will introduce a finite element structure on Ω = Ω  × Ω x × Ω v , associated with initial and corresponding boundary data. Here one may assume that f has compact support in the velocity space R d v , and hence assume homogeneous Dirichlet boundary condition for Ω v .
Thus to define an adequate finite element space we let where Here Ω stands for either Ω x or Ω x × Ω v . For k = 0, 1, 2, . . ., we define the finite element spaces for the Maxwell's equations (resp. Vlasov equation) as the space of piecewise polynomials which are continuous in x (resp. in x and v) and with possible discontinuities at the interior time levels t m , m = 1, . . . , M : with the extension for the Vlasov part and with τ = τ x × τ v , viz: where P p K (·) is the set of polynomial of degree at most p K on the given set. In this setting we allow the degree of polynomial to vary from cell to cell, hence we define the piecewise constant function p(t, x, v) := p K . We shall also use some notation, viz. ( where g = (g 1 , . . . , g 6 ) T , g ± (t, x) = lim s→0 ± g(t + s, x). The problem (3.1) is equivalent to: find W h,i ∈Ṽ h such that where the bilinear form is defined as [W ], g + m + W + , g + 0 and the linear form bỹ Now let (·, ·) K denote the L 2 -inner product over K and define a non-negative piecewise constant function δ by i.e., δ K is a non-negative constant on element K. Counting for the local character of the parameters h K , p K and δ K , to formulate a finite element method based on the local space-time elements, the problem (3.2) would have an alternative formulation where we replace in the definitions forB andL the sum of the inner products (·, ·) m involving δ K by the corresponding sum K∈C h (·, ·) K and all δ by δ K . Thus we have the problem (3.2) where in the bilinear, and linear, forms the first sum is replaced by K∈C h as, e.g., We also have that the solution W of equation (2.2) satisfies Subtracting (3.2) from this equation, we end up with the following relatioñ which is of vital importance in the error analysis. Now assuming jump discontinuities at the time levels t = t m , m = 1, . . . , M − 1, the suitable norm for stability and convergence would read as follows: 3.2. Vlasov-Maxwell Equations. The hp-streamline diffusion method on the ith step for the Vlasov part (2.3) can be formulated as follows: find f h,i ∈ V h such that for m = 0, 1, . . . , M − 1, where the trilinear form B is defined as and the linear form L is given by Analogously as for the Maxwell's equations we reformulate (3.5) considering phase-space-time finite element discretization. This yields replacing the first sum in (3.6) by a sum over the prismatic elements K ∈ C h of the form K∈C h and thus have the terms with Therefore, the adequate norm to derive stability and convergence estimates for the Vlasov equation will be the following triple norm:

Maxwell Equations.
Lemma 4.1 (M-coercivity). The bilinear formB(·, ·) satisfies the coercivity rela-tionB Proof. By definition ofB we have that Integrating by parts we get that Then, the proof follows immediately through adding all above terms.
For any positive constant C we have that for g ∈H 0 , Proof. For t m < t < t m+1 , we may write where in the second equality we used (4.1). Finally, by Grönwall's lemma we have that Now, integrating over [t m , t m+1 ] we obtain the desired result.
We proceed to the error analysis. First, letW be an interpolant of W in the finite dimensional discrete function spaceṼ h and denote by W h,i a solution to (3.2). Then we represent the error as the following split To estimate the convergence rate for both the Maxwell's and the Vlasov-Maxwell equations we use the Theorem 3.2 from [5], which are based on classical interpolation estimates by [10]. As a consequence of this Theorem we have the following bounds for the interpolation error where Ω stands for Ω x in the Maxwell's equations and Ω x × Ω v for the Vlasov case) and its gradient: where the sums are taken over all space-time elements of the triangulation of the domain, [0, T ] × Ω, 0 ≤ s K ≤ min(p K , k), with p K being the local spectral order. Closed formulas for Φ 1 and Φ 2 are given in Theorem 3.2 of [5]. A less involved formula for Φ 1 can be found in [16]. Now we state the following convergence theorem for the Maxwell's equations.
Partial integration gives the identities sinceη andξ have compact support in Ω x . Inserting these equations into the expression forB(η,ξ) we end up with the following equality Further, using some standard inequalities it follows that Now let us estimate the second term The above two inequalities and Lemma 4.2, with properly chosen C and the bound on p K h K , imply the estimate Hiding theξ-term on the right hand side in theξ-term on the left hand side, gives us the following inequality Thus, we have estimated |||ξ||| 2 M . This implies that (4.4) Now we have to estimate the interpolation error terms: For the term J 1 we use (4.2) and (4.3) to get To estimate the term J 2 we use the trace estimate combined with the inverse inequality and get to obtain where for all terms in J 2 we combined (4.6) with (4.2) and (4.3). Now, moving the triple norm ofẽ on the right hand side of (4.4) to the left hand side and using estimates (4.5) and (4.7) it follows that The next step is to apply a discrete Grönwall lemma of the form: If h ≤ 1/2C 2 , then we have that The discrete Grönwall's lemma yields . . , M . Plugging these inequalities into (4.8) will give the stated error estimate and the proof is complete.
Lemma 4.5 (V-coercivity). We have that Proof. Taking into account that ∇G(f h,i−1 ) = 0, g is zero on ∂Ω and following the proof of Lemma 4.1 we get the desired result.
Lemma 4.6 (Poincaré-type V-estimate). For any constant C we have for g ∈ H 0 , The proof is similar to that of Lemma 4.2 and therefore is omitted. Now we proceed with the error analysis. First we letf be an interpolant of f . Then we set We state the following convergence theorem.  10) and the parameter δ K on each K satisfies δ K = C 1 h K p K for some positive constant C 1 with p K h K ≤ C 2 < 1 for some constant C 2 > 0. Then there exists a constant C > 0 independent of p K , s K and h K such that where 0 ≤ s K ≤ min(p K , k) and the subscript V in the triple norm above, as well as a subscript for Φ is to emphasize that these quantities are in the Vlasov part.
Proof. By (3.5) and Lemma 4.5 we get that where and We start with the term T 1 . Integrating by parts and using the facts that η and ξ are zero on ∂Ω and ∇G(f h,i−1 ) = 0, we get Now using Cauchy-Schwarz inequality we obtain the estimate where for the last term we have the bound Next we estimate T 2 : To proceed we need to estimate . Hence using (4.9) we obtain where we denote . Now combining the estimates for T 1 and T 2 together with (4.12), (4.13) and (4.10) we have Moving the triple norm to the left hand side and estimating ξ 2 Q T ≤ C( e 2 Q T + η 2 Q T ), will give us the following inequality We now estimate |||e||| V as follows Using Lemma 4.6 for the term e 2 Q T with an appropriately chosen constant C and the bound on p K h K we get Recalling (4.2) and (4.3) for all η-terms we obtain (in a similar way as in the proof of Theorem 4.3 for the terms J 1 and J 2 ) Now we use the discrete Grönwall lemma stated in the proof of Theorem 4.3 and we get the following inequality which gives the desired result and completes the proof.
Proof. Using Lemma 4.6 with a properly chosen constant C and the bound on where S denotes the terms with sums from the right hand side of (4.11). Using this inequality repeatedly will give us f − f h,i 2 s K +1,K ≤ CS + C(ph) i , which ends the proof.
Remark 4.9 (The DG approach). The whole theory developed in the previous sections work for the streamline diffusion based discontinuous Galerkin (SDDG) method as well. In this method our approximations are also allowed to have jump discontinuities across the inter-element boundaries. Then the norms, including the sum over such jump terms, are much more involved. However, as we mentioned above, the analysis although lengthy follow the same path.

Nitsche's method for Maxwell equations
Alternative, yet more desirable numerical scheme for the Maxwell's equations can be obtained using a symmetrizing penalty approach known as Nitsche's method. This, however cannot be extended to the Vlasov part due to the hyperbolic nature of the Vlasov equation. Nevertheless, the advantages of the symmetrizing are overwhelming. Therefore below we include analysis of the Nitsche's approach for the Maxwell part.
Recall that we have the Maxwell's equations given by We differentiate the first equation with respect to time to get Now recall the Green's formula where n = n(x) is the unit outward normal to the boundary at the point x ∈ ∂Ω x . Apply Green's formula to (5.3) This relation being non-symmetric (see the contribution from the boundary terms) causes severe restrictions in, e.g., deriving stability estimates. To circumvent such draw-backs Nitsche introduced a symmetrized scheme for elliptic and parabolic problems (see [21]), which is also known as the penalty method. This can be seen in, e.g., [8] and [23]. In our case Nitsche's method is performed by the add of extra boundary terms making the bilinear form symmetric and coercive, viz.
Here γ is a constant that will be specified later. Now, we define the symmetric bilinear form and the element space of piecewise linear polynomials Thus, we can formulate the semi-discrete problem as: for each fixed t, It is straightforward to observe the consistency of the method.
Lemma 5.1 (Consistency). The exact solution E of (5.2) satisfies Now, defining the mesh dependent discrete norm where in the last inequality we used the following trace estimate: and the trivial inequality Now, if we choose the constants γ and α, such that γ > 4α and α >C, then we have proved the following coercivity result.
Lemma 5.2 (Coercivity). If γ is large enough, then there exists a constant C > 0, such that a(g, g) ≥ C g 2 h ∀g ∈ V x h . To have continuity of the form a(·, ·) we need to define a mesh dependent triple norm as |||g||| 2 h := g 2 h + h 1/2 ∇ x × g 2 Γx . Then we have the following lemma. Lemma 5.3 (Continuity). The bilinear form a(·, ·) is continuous with respect to the triple norm ||| · ||| h and we have the following estimate Proof. Using inequality (5.5) and simple algebra we get the result: For the triple norm ||| · ||| h we have the following inverse estimate which holds for all g ∈ V x h . We note that the trace estimate implies the coercivity of the form a(·, ·) also in the triple norm.

A modified Ritz projection. Let us define a projection
We have the following error estimates for the projection Q h : Proof. Consider the stationary problem The Nitsche formulation for this problem is given by: Further, we have the Galerkin orthogonality Hence, the projection Q h can be seen as the solution operator of (5.7). We therefore need an a priori error estimate of (5.7). To this end we split the error into two terms ϕ−ϕ h = (ϕ−I h ϕ)+(I h ϕ−ϕ h ) = η+ξ, where I h is the standard nodal interpolation operator. By coercivity, continuity of a(·, ·) and the Galerkin orthogonality (5.8) we have that |||ξ||| 2 h ≤ Ca(ξ, ξ) = −Ca(ξ, η) ≤ C|||ξ||| h |||η||| h . It follows that |||ξ||| h ≤ C|||η||| h , so it remains to estimate |||η||| h . Below we estimate each term in |||η||| h separately. For the interpolation error we have . As for the boundary integrals, by trace inequality we have the estimates and similarly where, in both estimates, in the last inequalities we have used the interpolation estimates. Summing up we end up with It remains to estimate the error in the L 2 -norm. To this end we consider the auxiliary problem Multiplying the first equation with ϕ − ϕ h and integrating over Ω x yields By the stability of the elliptic problem Ωx the estimate for the L 2 -norm follows.
Remark 5.6. For the magnetic field B we get a slightly different system of equations, but the same error estimates will hold.

Time discretization.
Let {t m } M m=0 be a uniform partition of [0, T ] of step size k = T /M . Before formulating the fully discrete problem, we introduce the following notations of difference quotients where u m = u(t m ). Then a fully discrete problem reads as follows: for m = 2, 3, . . . , M , find E m such that The choices of the first two approximations E 0 and E 1 will be discussed later. We split the error as We use Lemma 5.4 to estimate ρ m , hence it remains to estimate θ m . To do so we note that Then we have (5.14) We define the discrete energy Now, if in the left hand side of (5.14) we use the second form of χ in the first term and the third form in the second term, and in the right hand side we use the second form of χ in both terms, then we get We estimate the right hand side using the continuity of a(·, ·), with C a = (9 + γ), and the inverse inequality for the triple norm to get

After cancellation it follows that
Iterating the above inequality leads to Now we estimate the terms on the right hand side. Let us begin with ω j and split it as We write ω j 1 in the following way Summing over j and using Lemma 5.4 gives As for ω j 2 we use Taylor expansion of E j and E j−2 in polynomials of degree 2 about t j−1 . Then, due to cancellations, we end up with Once again summing over j leads to The third term in (5.15) will be estimated as follows It remains to estimate E 1 , which depends on how E 0 and E 1 are chosen. Let E 0 = Q h E 0 and assume that E 1 is chosen such that θ 1 Ωx ≤ C(h 2 + k 2 ). Then we have We combine the above estimates to get ∂ t θ m Ωx ≤ C(h + k) for m = 2, 3, . . . , M , assuming that k is proportional to h. Finally, if we use the following estimate then we also have that θ m Ωx ≤ C(h + k). We sum up the results in the following theorem.

Numerical Results
Here we present some numerical results justifying the accuracy of our method. We performed the calculations for the simplified case of one space variable and two velocities variables, i.e. the one and one-half dimensional Vlasov-Maxwell system (cf [22]), which takes the following form: We assume here the non-relativistic case of the Vlasov-Maxwell system, since there is a wider literature available to compare the numerical test with. We note that our theoretical results are also valid in this case. The initial conditions are given by which corresponds to the streaming Weibel instability (cf [9]) with β = 0.01 and b = 0.001. We perform the calculations for two sets of values of parameters:   Table 1 lists the errors for the fixed mesh set H 1 and increasing degree of finite elements polynomial approximation, whereas in Table 2 we list the errors for the fixed degree p = 1 of polynomial approximation and different mesh sizes.  We can see from the tables the convergence of our method for all functions. We present the results for one value of the stability parameter δ = 0.05, since its choice does not influence importantly (in some reasonable interval of values) the accuracy, but only the stability of the method.
6.2. Streaming Weibel instability tests. In this section we present the preliminary results for the streaming Weibel instability tests. More results will be included in the forthcoming paper [20].
We present the time evolution of the magnetic, electric and kinetic energies for both cases of parameters values. The calculations were carried out for the mesh sizes h t = 1/20, h x = 1/30, h v = √ 2 · 11/300 with p = 1. We plot the components of electric energy E i = 1 2L L 0 E 2 i dx 1 , i = 1, 2, and the magnetic energy B = 1 2L L 0 B 2 dx 1 in Figure 1. The kinetic energy is showed in Figure 2 as the separate components defined by K i = 1 2L L 0 Ωv v 2 i f dvdx 1 , i = 1, 2. The qualitative behaviour of the time evolution of both the electromagnetic and kinetic energies is in agreement with the theory and the results presented in [9] for different numerical methods.

Conclusion
This paper concerns two approaches in the numerical investigation for the Vlasov-Maxwell system. The first study is devoted to the hp streamline diffusion method for the relativistic Vlasov-Maxwell system in fully 3-dimensions. Our objective is to present unified phase-space (for Maxwell's equations) and phase-space-time (for the Vlasov-Maxwell system) discretization schemes that have optimal order convergence for the hyperbolic problems ( O(h s−1/2 ) for solutions in the Sobolev space H s (Ω) ) with strong stability properties and adaptivity features. The adaptivity in a priori regime is based on refining in the vicinity of singularities combined with lower order approximating polynomials and non-refined mesh with higher spectral order in smooth regions. In this way we have constructed a finite element mesh with several improving properties, e.g. stability, convergence, and adaptivity, gathered in it. To our knowledge, except in some work in convection-diffusion problems, see e.g. [17] and our study in [5], such approach is not considered for this type of equations elsewhere. The second study concerns a penalty method for the Maxwell's equations, which is based on a certain Nitsche type symmetrization scheme. In this part we have combined the field equations to a second order pde. For this equation we derive a second order spatial approximation for the Nitsche's scheme. We also prove a second order temporal discretization, assuming a somewhat more regular-in-time field functions. Even this approach is not considered in any other works for the VM system.
The results are justified in lower dimensional cases through the accuracy and the streaming Weibel instability tests presented in this paper and through implementing some numerical examples in the forthcoming paper, see [20]. The full-dimensions are to expensive to experiment. However, the theoretical analysis and numerical justifications in low dimensions are indicating the robustness of the considered schemes.