A PRIORI ERROR ANALYSIS OF NEW SEMIDISCRETE, HAMILTONIAN HDG METHODS FOR THE TIME-DEPENDENT MAXWELL’S EQUATIONS

. We present the first a priori error analysis of a class of space-discretizations by Hybridizable Discontinuous Galerkin (HDG) methods for the time-dependent Maxwell’s equations introduced in S´anchez et al . [ Comput. Methods Appl. Mech. Eng. 396 (2022) 114969]. The distinctive feature of these discretizations is that they display a discrete version of the Hamiltonian structure of the original Maxwell’s equations. This is why they are called “Hamiltonian” HDG methods. Because of this, when combined with symplectic time-marching methods, the resulting methods display an energy that does not drift in time. We provide a single analysis for several of these methods by exploiting the fact that they only differ by the choice of the approximation spaces and the stabilization functions. We also introduce a new way of discretizing the static Maxwell’s equations in order to define the initial condition in a manner consistent with our technique of analysis. Finally, we present numerical tests to validate our theoretical orders of convergence and to explore the convergence properties of the method in situations not covered by our analysis.


Introduction
The goal of this paper is to present the first a priori error analysis of a class Hybridizable Discontinuous Galerkin (HDG) methods for the time-dependent Maxwell's equations. This class of HDG methods is based on the formulation of the magnetic field and its vector potential. It is well-suited for symplectic time-integration and has the advantage of providing time-invariant, non-drifting approximations to the total electric and magnetic charges, and to the total energy, see [31].
To describe our results, let us begin by considering the classical formulation of the Maxwell's equations based on the electric and magnetic field: ∇¨p q " in Ωˆp0, s, ∇¨p q " 0 in Ωˆp0, s, complete with the following initial and boundary conditions: " on Γˆp0, s, Γ :" BΩ, " 0 , " 0 on Ωˆt " 0u.
Here we assume Ω Ă R is a polyhedral domain with Lipschtz boundary. In addition, we assume that and , both lying in 8 pΩq, are bounded from below by 0 ą 0 and 0 ą 0, respectively, and , P 8 pΩq are non-negative functions.
We are interested in finite element spatial discretizations of (1). A natural way is to use Nédélec's pcurlqconforming edge elements [27,28]. Compared to 1 -conforming nodal elements, they provide more robust convergence. For instance, it is well known that nodal elements can give spurious solution when the domain has re-entrant corners [12,13], or spurious modes when the domain has closed cavities [3]. Edge elements are free from these issues. For error analyses of edge elements for the time-dependent Maxwell's equations, we refer to [23,25,26].
A key reason for the above-mentioned lack of convergence by 1 -conforming elements is their heavy enforcement of the inter-element conformity along the normal direction, since this renders the discrete space not dense in the solution space. Therefore, it is not surprising to see that Discontinuous Galerkin (DG) methods are free from spurious solutions [5]. In addition, compared to edge elements, DG methods are more straightforward to implement for ℎ -adaptivty and meshes with hanging-nodes. Numerical studies also suggest DG methods reach high scaliability due to their small parallel communication footprint [24]. These advantages make DG an appealing class of methods to solve the Maxwell's equations. However, since no conformity is enforced in the discrete space, DG methods display more globally-coupled degrees of freedom for the same mesh, as compared with edge methods. To overcome this problem, the Hybridizable Discontinuous Galerkin (HDG) methods were introduced. In the framework of electromagnetic problems, they first appeared in [30] for solving the timeharmonic Maxwell's equations. Error analyses are given in [6][7][8][9]15,22] for either the time-harmonic or the static Maxwell's equations.
Relatively fewer papers exist on HDG methods for the time-dependent Maxwell's equations; see [10,19,21,29]. Most existing methods are based on a direct application of the HDG discretization of (1), namely, by using approximations of and on each element, as well as approximations of the tangential traces of and on the skeleton of the mesh. In [31], we show that the methods obtained in this way, that is, by using HDG space discretiations, are not a Hamiltonian dynamical system. As expected, we numerically observe a drifting of the total energy, the electric and magnetic charges, and the linear and angular momenta. This makes -based methods unsuitable for long-time simulation. To overcome this problem, we proposed a new class of methods based on a rewriting of the original problem (1) in terms of the electric field and the magnetic potential . The corresponding -formulation [31] reads as follows: : " ∇ˆ`9´in Ωˆp0, s, " ∇ˆin Ωˆp0, s, ∇¨p q " in Ωˆp0, s, and is complete with the following initial and boundary conditions: " on Γˆp0, s, Γ :" BΩ, " 0 , 9 "´0 on Ωˆt " 0u.
Here is the magnetic vector potential; see (2b). Note that the condition (2c) (or (1c)) only needs to be enforced in the initial condition (2e), and similarly for (1d); see [25]. Correspondingly, in the HDG methods we shall introduce later based on (2), the equation (2c) is not used in the discretization of the evolutionary equation (3), but is used in the discretization of the initial condition (5e). We proved that the HDG methods obtained based on (2) are Hamiltonian dynamical systems. Consequently, they provide non-drifting energy, charges and momenta; see Corollary 4.2 of [31] for the proof for these non-drifting properties.
In this paper, we present an a priori error analysis for the HDG methods based on the formulation (2). To the best of our knowledge, this is the first error analysis for HDG methods for the time-dependent Maxwell's equations. By adopting the techniques used in [15], we are able to analyze at once most existing spatial HDG discretizations, see Table 1. Moreover, since (2) is essentially a wave-form formulation based on the magnetic potential, we are able to reuse many existing analysis techniques from [11,14] for acoustic and elastic waves. However, for the Maxwell's equations special attention needs to be paid to the discretization of the initial condition in order to ensure that the initial error does not pollute the convergence order for the rest of the time. To achieve this, a new class of HDG methods for discretizing the static Maxwell's equations is introduced.
We present our main result in Section 2. We describe the class of HDG semi-discretizations as well as the choice of the approximation to the initial condition. We then state and and discuss our error estimates. Then in Section 3, we present its detailed proof. Finally, in Section 5, we show numerical tests which validate our theoretical results.

Notation and preliminaries
We begin by introducing some notation. Let ℎ be a collection of conforming triangulations of Ω by shaperegular, star-shaped polyhedral elements. Let ℱ ℎ denote the set of all the faces of ℎ , and ℱ the set of faces of an element P ℎ . For any element P ℎ , we denote its diameter by ℎ , and set ℎ :" max P ℎ ℎ and ℎ m :" min P ℎ ℎ .
Let be either an element P ℎ or a face P ℱ ℎ . We denote by p q and p q the polynomial spaces of degree , supported on , taking values in R and R , respectively. We let Π and Π denote the 2 -projection into p q and p q, respectively. For any vector normal to a face , , and any linear space of R -valued functions defined on , p q, we set p q :" t P p q :¨" 0u, and p q :" t P p q :ˆ" 0u. For an R -valued space p q defined on an element , we shall denote pB q :" ś Pℱ p q, where ℱ is the collection of the faces of , and p q is the restriction of the space p q on the face .
For any subdomain of R , we set p , q :" ş for scalar-valued functions and , and p , q :" ş¨v for R -valued functions and . We then write p , q ℎ :" For any natural number , we set where, as usual, | | 2 , :" Finally, for functions depending on time and space, and for any natural numbers ℓ, , we set ' ∇˜`2 [7] where p q represents the -th order time derivative of the function . We write 8 instead of 0,8 , and 2 instead of 0 . Similar definitions hold when replacing ℎ by B ℎ . We now introduce the notation for the approximation spaces of our numerical methods. Given any element P ℎ , we let p q and p q represent the approximation spaces for the magnetic vector potential and the magnetic field on , respectively. Since 9 "´, p q is also the space for the electric field . For any face P ℱ ℎ , we denote by p q the approximation space for the tangential component of the magnetic potential on . We also set pB q :" ś Pℱ p q where p q can be p q, p q, or p q. With the local spaces defined, we introduce the global approximation spaces: where 0 ℎ and B ℎ represent the subspaces of ℎ taking zero values on the boundary faces and the interior faces, respectively. Finally, we denote by the maximal polynomial degree such that p q Ă p q or p q Ă p q hold for all P ℎ or P ℱ ℎ , where can be , , or .

Spatial semi-discretization
We assume the approximation spaces , and satisfy the following assumption: Assumption 2.1 (On the local spaces). For all P ℎ , we have that t P p q : ∇ˆ" 0u Ă pB q.
Assumptions 2.1-(1), -(2) and -(3) were first used in [15] for a unified analysis of HDG methods for the static Maxwell's equations. Assumption 2.1-(4) is added in order to guarantee the existence of a class of projections needed for our error analysis; see Proposition 3.1.
Note that we do not require pB q Ă pB q. Thus, pB q can be larger than pB q. To handle this case, we use the 2 -projection : 2 pB ; R 3 q Ñ pB q. The appearance of projections of this type is a hallmark of the Lehrenfeld-Schöberl projection for HDG methods, which was first introduced in [20] and later applied to the static Maxwell's equations in [7]. A remarkable feature of these methods is that they achieve super-convergence (of their numerical traces) on general polyhedral meshes.
Our analysis holds for any HDG method satisfying the above Assumption 2.1 on the local spaces. Examples are the three HDG methods listed in Table 1. ,0 For instance, for the methods in Table 1, we choose ,0 ℎ :" ℎ X 1 0 pΩq where ℎ :" " for HDG-S and "`1 for HDG-B /B +. The equations (5) define a discretization of the following static Maxwell's equations: in Ω, where the curl operator is discretized by an HDG method while the gradient operator is discretized by Lagrange elements. Here is a pseudo-pressure. Note that ∇ˆp0q is divergence-free, which leads to " 0. This property is inherited by the discretization (5). In fact, we have the following Proposition 2.1. Its proof is included in Section 3.
In our previous work [31], to define the initial condition, we used discontinuous approximations for the pressure and for its trace over the skeleton, and weakly imposed the continuity of the numerical normal trace of ℎ . In contrast, here, we propose an alternative approach to discretize the initial condition for two reasons. The first is that, the method has less degrees of freedom since it uses continuous approximations to the pressure . The second is that, we can guarantee that ℎ " 0, see Proposition 2.1, which makes this initial discretization consistent with the time-marching scheme since the time-marching (3) does not involve the pressure.

Error estimates
Now we present the error estimates to the HDG methods defined by (3)-(5). We first present two assumptions which are only used in the duality arguments for estimating the error of ℎ . Assumption 2.3 (On the data and elliptic regularity). We assume that: (1) there exist constants reg, in p´1 2 , 1s and reg, P p 1 2 , 1s such that the following elliptic regularity inequality holds: for all P pcurl; Ωq such thatˆ" 0 and ∇¨" 0; (2) the functions , , and are piecewise 1 , namely, there exists a constant , , ą 0 such that sup P ℎ p} } 1 p q`} } 1 p q`} } 1 p q q ď , , ; (3) the domain Ω is simply connected with a connected boundary.
When is a constant function, the elliptic regularity inequality of Assumption 2.3-(1) holds by Proposition 3.7 of [1]. If in addition, Ω is a convex polyhedron, then it holds with reg, " 1 and reg, " 1. For the cases when has a lower regularity, we refer to [2,16].
Assumption 2.4 (Minimal polynomial degree). For all P ℎ , we have (1) 0 p q Ă p q, (2) 1 p q Ă p q, Table 2. Orders of convergence for any HDG method satisfying the assumptions of Theorem 2.1. We assume that the solution is very smooth and the we have the full elliptic regularity, that is, reg, " reg, " 1.
We are now ready to state our main result. It provides estimates of the quality of the approximate magnetic field, ℎ , magnetic potential, ℎ , and the electric field ℎ :"´9 ℎ .
Theorem 2.1. Consider any HDG method defined by the weak formulation (3), and the initial condition (4) and (5). Assume that its approximation spaces satisfy the Assumptions 2.1 and 2.2. Then we have .
Here the constant 1 depends only on the polynomial degree , , , the shape regularity of ℎ , and the parameters , , . The constant 2 depends additionally on the domain Ω through the elliptic regularity constant reg .
Observe that we can have an estimate for the approximation error of the magnetic vector potential ℎ without assuming the data regularity by simply using the estimate for the electric field. This would give a sub-optimal order of convergence, though. For the main choices of the stabilization parameter , this result gives the orders of convergence displayed in Table 2. The orders of convergence of the HDG methods in Table 1 can now be easily obtained, see Table 3. Let us briefly discuss them.
First, note that for " ℎ´1, the three methods considered here provide approximations to converging optimally. Next, note that we chose to consider the method HDG-S`1, and not the method HDG-S , because Table 3. Orders of convergence for the HDG methods listed in Table 1. The conditions of Theorem 2.1 are assumed to hold. In boldface are indicated the optimal orders of convergence. We assume that the solution is very smooth and that we have the full elliptic regularity, that is, reg, " reg, " 1.
this method has the same spaces than HDG-B for both the magnetic potential ℎ and its trace p ℎ . Only the local space for the magnetic field, ℎ , is bigger for HDG-S`1. We see that this difference pays off for " ℎ and " 1 as all the unknowns of the method HDG-S`1 converge with one additional order of convergence. The difference disappears for the larger stabilization " ℎ´1 though. Finally, note that all the variables of the methods HDG-B and HDG-B + converge with the same orders, even though the local space for the numerical traces, p q, is smaller for the latter. In fact, the numerical trace p ℎ provided by the HDG-B + method achieves super-convergence since it lies in a proper subspace of`1.
The estimates of Theorem 2.1 hold for general polyhedral meshes. Since using polyhedral meshes allows an automatic handling of hanging nodes, our analysis naturally extends to non-conforming meshes.

Proofs
In this section, we prove the results of Section 2. We proceed in several steps.
Step 1: The initial condition Let us begin by proving the unique solvability of the discretization of the initial condition (5) of Proposition 2.1.
It remains to prove that the approximate pressure ℎ is identically zero. Again, by Assumption 2.2-(1), we can set :" ∇ ℎ in equation (5a) to get . Finally, by Assumption 2.2-(1), ∇ ℎ P ℎ and so, p∇ ℎ q P ℎ by Assumption 2.1-(4). Since ℎ " 0 on the boundary of Ω, by Assumption 2.2-(3), we have that p∇ ℎ q P 0 ℎ . This means that we can take 0 :" p∇ ℎ q in equation (5c), to get that As in the previous argument, we can now conclude that ℎ " 0. This completes the proof of Proposition 2.1.
Step 2: The compact form of the weak formulation Here, we introduce a compact form of the weak formulation (3), which allows us to investigate its properties more concisely. So, we define Note how the lines of the definition of the bilinear form ℎ correspond to the equations (3a), (3c), and (3b), respectively. Now, the original weak formulation (3) can be rewritten into the following more compact form: find p ℎ , p ℎ , ℎ q P ℎˆℎˆℎ such that for all p , , q P ℎˆ0 ℎˆℎ , with the boundary condition enforced by for all P B ℎ . The bilinear form ℎ has a sub-skew-symmetry structure, as we see in the next result. This structure will be frequently exploited in our error analysis of the method. Lemma 3.1 (Sub-skew-symmetry). The following identities hold: ℎ p , p, ; , p, q " ℎ p , p,´; , p,´q, ℎ p , p, 0; 0, 0, q`ℎp0, 0, ; , p, 0q " 0, ℎ p , p, ; , p, q " ℎ p , p, 0; , p, 0q`ℎp0, 0, ; 0, 0, q, " x p´pq,´py B ℎ`p , q ℎ , provided p " p. Proof. We obtain the first identity by examining the structure of (6). The next two identities follow from the first one. The last identity is a consequence of the definition of ℎ . This completes the proof.
Step 3: Projections Here, we introduce some projections and study their properties. Our analysis consists in controlling approximation errors of these projections.
hold for all P ℎ . In addition, where depends only on the shape regularity of , and the polynomial degree and .
If we define Π , on each element P ℎ , as the element of p q which solves the equations pΠ´, ∇ˆq " xˆ´pˆq, y B @ P p q, (10a) where K means taking the orthogonal complement in p q, we will have that it satisfies the equation (8a), by construction. So, it remains to prove that it Π is well defined. To do that, we only have to prove that the above set of equations defines a square system. The fact that Π is well defined will then follow from the inequality (9a). Let us show that the system of equations is square. The number of unknowns is clearly dim p q. Let us count the number of linearly independent equations. By Assumption 2.1-(4), we know that, if ∇ˆ" 0, then P pB q and xˆ´pˆq, y B " 0. As a result, the enforcement by (10a) only has a rank of dimp∇ˆp qq. Considering that (10b) has a rank of dimp p qq´dimp∇ˆp qq, we see that the equations (10) define a square system. It remains to prove the inequality (9a). Let be the 2 -projection to ℎ . Now set :" Π´P ℎ . On each element P ℎ , we can write " 1`2 where 1 P ∇ˆp q and 2 P p∇ˆp qq K . By equation (10b), 2 " 0 and so, " 1 " ∇ˆfor some P p q. Then, by equation (10a), we have for any P p q such that ∇ˆ" 0. Applying the Cauchy-Schwarz inequality, we get since the operator is a bijection from the orthogonal complement in p q of its kernel to ∇ˆp q. Since ∇ˆ" 1 " , we get This completes the proof of (9a).
Step 4: The equations of the projections of the errors for ą 0 Now that the projection Π , is defined, we can obtain the equations satisfied by the projection of the errors. To do this, we first introduce an associated boundary remainder operator: The boundary remainder can be regarded as an interpolation error on the element boundary. It also serves as a quantification for the "non-conformity" of the projection operator Π , . In the following result, we place its most important property, namely, a generalized version of the weak-commutativity property.
for all P p q.
Proof. By the first equation defining the projection Π , (8a), we have The proof is complete.
Next, we obtain an identity for the bilinear form ℎ involving the projections of our unknowns. Lemma 3.3 (An identity for the bilinear form ℎ ). Suppose that Assumption 2.1 holds. Set :" Π , p , q.
for all p , , q P ℎˆ0 ℎˆℎ .
As we see in the next result, the projections of the errors satisfy a set of equations which share a similar structure to the HDG discretization defined by (7).
for all p , , q P ℎˆ0 ℎˆℎ .
Proof. The equation can be directly obtained by taking the difference between the equation of the projections of Lemma 3.3 and the compact form of the HDG semi-discretization (7). This completes the proof.
Step 5: The energy argument With the error equations just obtained, we can obtain the following energy identities. , : ℎ˘ℎ``: , : ℎ˘ℎ .
To prove these estimates, we are going to rely on the following inequality. Proof. We provide a proof for the sake of completeness. Consider the interval r0, s and let p˚q maximize in the interval. Then we have This completes the proof.
We are now ready to prove Proposition 3.3.
Then by the first energy identity of Proposition 3.2, we know , , , and satisfy the condition of Lemma 3.5. This proves the first estimate. Similarly, by letting : , : :" 2 p0q`p0q p0q.
Then we know, from the second energy identity of Proposition 3.2, that , , , and satisfy condition of Lemma 3.5. This completes the proof of Proposition 3.3.
Step 6: The equations of the projection of the errors at the initial time Lemma 3.6 (Error equations at the initial time). Suppose that Assumptions 2.1 and 2.2 hold. Then ℎ`ℎ p0q, p ℎ p0q, ℎ p0q; , ,˘" p p0q, q ℎ`x p0q,´y B ℎ , for all p , , q P ℎˆ0 ℎˆℎ .
Proof. This result can be proved in the same way as Lemma 3.4 by exploiting the similarity between the equations defining the weak formulation of the method for ą 0 (3), and the equations defining its initial condition (5). We also must use the fact that the approximate pressure ℎ is zero, by Proposition 2.1. This completes the proof.
Step 7: Estimates of the errors at the initial time The following result gives an estimate of the initial errors Er 0 p0q and Er 1 p0q.
The above identity suggests that ℎ can degrade the order of convergence for Er 1 p0q, and then consequently also } 9 ℎ p q}, } 9 ℎ p q}, and } 1{2 p 9 ℎ p q´9 p ℎ p qq} B ℎ . Since these term are used for the duality estimate of ℎ p q, which we shall present in the next subsection, the error pollution from ℎ may also have an effect on the convergence of } ℎ p q}. Our choice of the initial condition by (5) is free from these problems since ℎ " 0 by Proposition 2.1.
Step 8: Proof of the first estimate of Theorem 2.1 Combining the estimates obtained in Propositions 3.3 and 3.4, we can bound the projection of the errors , as we get that where ℎ " Π´ℎ, ℎ " Π´ℎ, and 1{2`ℎ´p ℎ˘" 1{2´p Π´ℎq´´´p ℎ¯¯.
and, by the triangle inequality, we get that The wanted estimate now follows from the approximation estimates (9): This completes the proof.
Step 9: The dual problem In what follows, we give an estimate of ℎ p q by using a duality argument, also called the "Aubin-Nitsche trick" in second-order elliptic problems.
Since we are going to use the following characterization of the 2 -norm of ℎ p q; We need to introduce the corresponding dual problem, namely, with˚p q " 0 and 9˚p q " Θ P 2 pΩq. We next define the projections of the dual solution, namely Π˚and Π˚, and their corresponding boundary remainder:˚: For any time-dependent function p q with P r0, s, we introduce the following notation for the time-integration from to : Note that, since Π , is linear,˚" Π , p˚,˚q.
Proposition 3.5. Suppose Assumption 2.3 holds. We have Proof. To obtain the first inequality, we use a simple energy argument. Indeed, if we set p q :" p 9˚p q, 9˚p qq Ω`p´1 ∇ˆ˚p q, ∇ˆ˚p qq Ω , we get, by (14), that 9 p q " 2p 9˚p q, 9˚p qq Ω ě 0. So, p q ď p q " } 1{2 Θ} 2 Ω for all P r0, s. It remains to estimate the term }˚} 8 p , 2 pΩqq . To do that, we use the Poincaré-type inequality: The above inequality holds for a bounded and simply connected domain Ω with a Lipschitz connected boundary, see Proposition 16.9 of [32]. This condition for Ω is satisfied here since we have assumed Assumption 2.3-(3). Now, since ∇¨˚" 0 andˆ˚" 0, we have }˚} Ω À }∇ˆ˚} Ω . This completes the proof of the first estimate.
To obtain the second inequality, we do use an elliptic regularity result. By integrating (14a) and (14b) from to , we obtain ∇ˆ˚p q "´9˚p q`9˚p q´˚p q, p q´∇ˆ˚p q " 0, hold for all P r0, s. Then by the elliptic regularity inequality of Assumption 2.3-(1), we obtain Now we use the first inequality. This completes the proof.
All we have to do now is to estimate each of the terms 0 , ş 0 , " 1, 2, 3, 4 and conclude. Before that, we introduce a lemma on the projection Π and Π , which will be used frequently in the estimates of the terms.
To do that,we are going to use several times the following auxiliary result.
Proof. Note that by Assumption 2.4, we know 1 p q Ă p q. In addition, recall that we define Π as the 2 projection to p q. Thus, we have pΠ´, 1 p qq " 0. To prove pΠ´, q " 0, for all P 0 p q, first recall the definition of Π , namely (10). By Assumption 2.4, we knowˆ1p qˆĂ pB q and 1 p q Ă p q. Now consider (10a), we take P 1 p q, then we have pΠ´, ∇ˆ1p qq " 0.
This with the fact that 0 Ă ∇ˆ1 completes the proof.
Step 11: Estimate of 0 Lemma 3.8. Suppose Assumptions 2.1-2.4 hold. Then we have Proof. By the definition of 0 , we have by the first inequality of Proposition 3.5 (which requires Assumption 2.3 to hold). It remains to estimate the quantity } ℎ p0q} ℎ . Since we use a duality technique to get the estimate. Thus, we begin by introducing the dual probleḿ ∇ˆp´q " 0,´∇ˆ" in Ω, " 0 on BΩ, and writing that`, by Lemma 3.3. Next we use the first identity of Lemma 3.1 to get By the error equation at the initial time, Lemma 3.6 (which requires Assumption 2.2 to hold), we obtaiǹ Setting p0q :" Π , p p0q, p0qq and 1 :" p , q, we arrive at : 01`02`03`04 .
Similarly, we have To estimate 04 , first note that So we have Finally, we combine the estimates for 0 with " 1, 2, 3, 4, and then use the first estimate of Proposition 3.4 to control the terms } ℎ p0q} ℎ and } 1{2 p ℎ p0q´p ℎ p0qq} B ℎ . This completes the proof.

So we havěˇˇˇˇż
Then we use the definition of˚, namely, equations (15) and (11), and the projection convergence properties (9) to estimate }´1 {2˚} B . Note that Assumption 2.1 is needed for˚to be well-defined. This gives ušˇˇˇˇż Next, we use Proposition 3.5 (which needs Assumption 2.3 to hold) to estimate p}˚} reg, , ℎ`}˚} reg,`1 , ℎ q.
Similarly, we carry out the estimate for | ş 0 12 | as follows:ˇˇˇˇż This completes the proof.
Step 13: Estimate of ş Proof. Recall that 2 :" p pΠ˚´˚q, ℎ´q ℎ`p pΠ´q,˚q ℎ ": 21`22 . We then use Assumption 2.3-(2) which says that is piecewise 1 . So we havěˇˇˇˇż where Proposition 3.5 is used in the last step. This completes the proof.
Finally we use Proposition 3.5 to control }˚} 8 p , 1`r eg, p ℎ qq . This completes the proof.
We finally use Proposition 3.5 to control the˚related terms. This completes the proof.
Step 16: Conclusion It only remains to prove the second inequality of Theorem 2.1. But, to do that, we simply have to gather the estimates of the terms 0 , ş 0 , " 1, 2, 3, 4, obtained in the previous steps. Then we use the second estimate of Proposition 3.3 to control sup Pr0, s Er 1 p q and Proposition 3.4 to control Er 0 p0q. This completes the proof.

Fully discrete scheme
In this section we present the fully discrete method. First, we cast the semidiscrete scheme (3) as a first order system of differential equations and then discretize using arbitrary-order, diagonally implicit Runge-Kutta methods. Alternative schemes such as explicit Runge-Kutta, partitioned Runge-Kutta, and fully implicit Runge Kutta methods can also be considered. See [18] for a comparison of explicit and implicit HDG methods for wave equation. Symplectic DIRK schemes are chosen here to showcase the energy conservation of the numerical scheme and corroborate the estimates in the error analysis. (1) , " ℎ`∆ ř´1 "1 , for " , .
(2) Solution p , ℎ , , ℎ , , ℎ , p , ℎ q P ℎˆℎˆℎ ,ˆℎ of the system: for all p , , 0 , B q P ℎˆℎˆ0 ℎˆB ℎ . : Finally, update the approximate solution`1 ℎ " ℎ`∆ ř "1 , for " , , , p . The fully discrete system is written for a general -stages DIRK scheme. Here we can consider for instance singly DIRK methods, symplectic DIRK methods, symmetric DIRK methods, explicit first-stage, just to name a few. See [17] for a complete review of DIRK methods.

Numerical examples
In this section, we present several numerical examples demonstrating the properties of the semidiscrete scheme (3), implemented with the spaces HDG-S and HDG-B (see Tab. 1), and discretized in time using DIRK methods. In the first two examples, we study numerically the error convergence; in the third example, we provide numerical support for our choice of initial condition (5); and in the last example we show the energyconserving property of the scheme. All experiments in this section were implemented using the open source finite element software NETGEN/NGSolve [33,34].

Example 1: History of convergence test
Here we validate the predictions given in Theorem 2.1 for the methods HDG-S of index and HDG-B of index combined with with a DIRK scheme of order`2. For each of the approximations ℎ , ℎ , and ℎ , we compute the maximum over the time steps of the 2 -errors of the corresponding error, and then estimate their orders of convergence (e.o.c.). For instance, for the magnetic vector potential approximation we compute errorpℎq " max } p q´ℎ} 2 pΩq 3 , e.o.cpℎq " where ℎ 1 correspond to the previous mesh size parameter used in the computations. The experiment is carried on the unit cubic domain Ω " p0, 1q 3 and using the solution of problem (2) given by p , , , q "¨c osp q sinp q sinp q sinp ? with material parameters: permittivity " 1, permeability " 1, and conductivity " 0.5. The remaining data of the problem are computed accordingly to the exact solution.
We present in Tables 4 and 5 the errors and orders of convergence for the DIRK-(HDG-S ) and DIRK-(HDG-B ) methods for " 0, 1, 2, 3, for a sequence of uniform triangulations (tetrahedra) with meshsize parameter ℎ and number of elements .
Optimal convergence of the 2 errors is observed in Table 5 for method DIRK-(HDG-B ) as predicted in Theorem 2.1 and Table 3. We observe convergence of order`2 for the magnetic vector potential ℎ , and of order`1 for the magnetic field ℎ . The convergence of the electric field ℎ lies in between the convergence orders of ℎ and ℎ , which is consistent with our theoretical prediction.
In Table 4, we observe a slightly faster convergence than the one predicted in Theorem 2.1 and Table 3, namely, that of`1{2 for the electric and magnetic fields and of`1 for the magnetic vector potential.

Example 2: Low-regularity test
We now consider the analytic solution of problem (2) with low regularity defined on the L-shaped domain Ω " pp´1, 1q 2 zpp´1, 0q 2 qˆp0, 1q by p , , , q " ∇ p , q cosp q, p , q " 4{3 sinˆ4 3˙, where is given in cylindrical coordinates p , q. Observe that and have a singularity located at the z-axis, and , P r 4{3´p Ωqs 3 , ą 0. Similar test are presented in [7,15] for the static case. We compute the numerical approximations on a sequence of uniform triangulations with meshsize ℎ and number of elements by the methods DIRK-(HDG)-S and DIRK-(HDG-B ) for " 0, 1. Errors and rates of convergence are displayed in Tables 6 and 7.  Table 6 shows optimal order approximation of the errors using the HDG-S method and with stabilization parameter " 1. Indeed, it shows order 1 for and in the case " 0, and (approximately) order 4{3 in the Table 7. History of convergence of numerical approximations to Maxwell's equations (2) by a DIRK method of order`2 and the HDG-B method of degree , with " 1{ℎ, up to " 2, CFL " 0.05.

Example 3: Initial condition test
The estimates obtained in our main result, Theorem 2.1, rely on the particular definition of the initial condition given by (5). In this numerical experiment, we provide evidence indicating that the use of this initial condition is indeed essential for the optimal convergence of the scheme. We compare the performance of the DIRK-(HDG-B ) method initialized with: a) the standard 2 projection of the initial data 0 and 0 onto the HDG spaces, and with b) the HDG-approximation proposed in (5).
In Figure 2, we consider the DIRK-(HDG-B ) method, and plot the 2 -errors, of the electric field and the magnetic vector potential, vs. the corresponding meshsizes ℎ. In red, are the results obtained by using the standard 2 -projection as initial condition. In blue are the results obtained with our new initial condition (5). We observe sub-optimal rates of convergence corresponding to the scheme initialized with the 2 -projections. This supports our preference for our choice for the initial condition.
We present, in Figure 3, a plot of the magnitude of the magnetic vector potential at time " 20 and the evolution of the energy associated to the system, specifically we plot the quantities As it is shown in [31], the method DIRK-(HDG-B ) method conserves the energy of the semidiscretization in space, ℋ ℎ , up to 10 digits of precision, and does not allow the "exact" energy, r ℋ ℎ , to drift in time.

Concluding remarks
In this paper, we provide the first a priori error analysis of a wide, new class of Symplectic-Hamiltonian finite element methods for the time-dependent Maxwell's equations, see [31]. We provide numerical results which validate our theoretical predictions and confirm that the schemes conserve exactly the energy of the discrete system and that the exact energy does not drift in time. Moreover, we propose a new way of defining the initial data which uses a continuous approximation for the "pressure". It behaves better than the simple L 2 -projection, as our numerical experiments suggest.