Numerical study of conforming space-time methods for Maxwell’s equations

Time-dependent Maxwell’s equations govern electromagnetics. Under certain conditions, we can rewrite these equations into a partial diﬀerential equation of second order, which in this case is the vectorial wave equation. For the vectorial wave, we investigate the numerical application and the challenges in the implementation. For this purpose, we consider a space-time variational setting, i.e. time is just another spatial dimension. More speciﬁcally, we apply integration by parts in time as well as in space, leading to a space-time variational formulation with diﬀerent trial and test spaces. Conforming discretizations of tensor-product type result in a Galerkin–Petrov ﬁnite element method that requires a CFL condition for stability. For this Galerkin–Petrov variational formulation, we study the CFL condition and its sharpness. To overcome the CFL condition, we use a Hilbert-type transformation that leads to a variational formulation with equal trial and test spaces. Conforming space-time discretizations result in a new Galerkin–Bubnov ﬁnite element method that is unconditionally stable. In numerical examples, we demonstrate the eﬀectiveness of this Galerkin–Bubnov ﬁnite element method. Furthermore, we investigate diﬀerent projections of the right-hand side and their inﬂuence on the convergence rates. This paper is the ﬁrst step towards a more stable computation and a better understanding of vectorial wave equations in a conforming space-time approach.


Introduction
The time-dependent Maxwell's equations and their discretization are required in many electromagnetic applications.One application is the modeling of an electric motor which is governed by the Ampère-Maxwell equation.The Ampère-Maxwell equation relates the time-dependent electric field with the current density and the magnetic field.If the corresponding space-time domain is star-shaped with respect to a ball, then we can rewrite Maxwell's system into the vectorial wave equation.Application of this technique in the static case is investigated in [34].Extensive studies of the static case can be found not only for analytic methods, but also for the approximation by finite elements, see e.g.[7], or in a more applied static nonlinear example [19].Moreover, the quasi-static case is examined, see [6].
However, the complete space-time setting of these equations has been studied in less detail.Most approaches for time-dependent Maxwell's equation are either dealing with the frequency domain, see [28], or using different types of time-stepping methods, see [21,Chapter 12.2].For the latter, proving stability is delicate.One solution is to stabilize the systems which stem from time-stepping methods.The development of such stabilized methods is an active research field, see e.g.[35] on energy-preserving meshing for FDTD schemes or [4] for the Cole-Cole model which involves polarization.Other classical time-stepping methods are the so-called leapfrog and the Newmark-beta methods, where the first is a special case of the second.Comparison and improvement of these methods can be found in [8].Another class of time-stepping approaches is locally implicit time integrators, see e.g.[18] and references there.
An alternative to time-stepping methods is a space-time approach.Most approaches apply discontinuous Galerkin methods to Maxwell's equations, see [1,10,11,24,36] and references there.Note that these methods are non-conforming in general.
In this paper, we derive conforming space-time finite element discretizations for the vectorial wave equation without introducing additional unknowns, i.e. without a reformulation of the vectorial wave equation as a first-order system.The advantage of staying with the second-order formulation and of conforming methods is the need for a smaller number of degrees of freedom compared with first-order formulations or discontinuous Galerkin methods.First, we state the variational formulation for different trial and test spaces, i.e. a Galerkin-Petrov formulation.Second, using the so-called modified Hilbert transformation H T , introduced in [31,39], we present a variational formulation with equal trial and test spaces, i.e. a Galerkin-Bubnov formulation.The main goal of this paper is to derive and describe in detail these two conforming space-time methods, the assembling of the resulting linear systems and their comparison concerning stability and convergence.In more detail, we investigate the Galerkin-Petrov formulation and the corresponding CFL condition.Additionally, we elaborate on a Galerkin-Bubnov formulation using the modified Hilbert transformation, where all necessary mathematical tools are developed.In the end, we compare the results of the conforming Galerkin-Petrov and Galerkin-Bubnov finite element methods concerning their stability and convergence behavior by considering numerical examples.Moreover, we investigate the convergence behavior of both conforming space-time approaches when the right-hand side is approximated by two different projections.
Before we introduce the space-time approach, we derive the vectorial wave equation from Maxwell's equations to get a better understanding of the properties of the equation.We state the main ideas but refer to [16,33] and references for a more detailed derivation.To derive the vectorial wave equation, we need to rewrite Maxwell's equations using differential forms in a Lipschitz domain Q ⊂ R 4 .Hence, we define the Faraday 2-form F := e ∧ dt + b, where e is the 1-form corresponding to the electric field E and b the 2-form corresponding to the magnetic flux density B. Additionally, we define the Maxwell 2-form G := h ∧ dt − d and the source 3-form J := j ∧ dt − ρ, where h is the 1-form corresponding to the magnetic field H, d is the 2-form corresponding to the electric flux density D, j is the 2-form corresponding to the given electric current density j and ρ is the 3-form corresponding to the given charge density ρ.Following the derivations in [33], Maxwell's equations allow for the representation where the four-dimensional exterior derivative d is deployed.Here, the Hodge star operator ,−µ −1 is weighted by for the components of (dx 01 , dx 02 , dx 03 ) and by (−µ −1 ) for the components of (dx 23 , dx 31 , dx 12 ) , considering the Euclidean metric in R 4 .
Next, we derive the space-time vectorial wave equation from the 4D representation (1) of Maxwell's equations.Since we assume that the domain is star-shaped with respect to a ball, the Poincaré lemma [22,Theorem 4.1] is applicable.Hence, the closed form F is exact and there is a potential A, which is a 1-form such that dA = F. From this, we also derive the following relations in the Euclidean metric Here, the scalar function A 0 is the time component and the vector-valued function A := (A 1 , A 2 , A 3 ) is the spatial component of A. If we insert this result into the second equation of (1), combined with the third of (1), we get the second-order differential equation We emphasize that there is no uniqueness of the potential A. Indeed, adding to A the exterior derivative of any 0-form φ results again in a suitable potential Ã = A + dφ which solves equation (2).This is equivalent to adding the space-time gradient ∇ (t,x) φ of a function φ ∈ H 1 0 (Q) to (A 0 , A) , with the usual Sobolev space H 1 0 (Q).In physics, this is called the gauge invariance, which 'is a manifestation of (the) nonobservability of A', see [20, p. 676].To make the potential A unique, we use gauging.A survey of different gauges can be found in [20].The choice of the gauge determines the resulting partial differential equation.In this paper, we deploy the Weyl gauge which is also named the temporal gauge, [20].The Weyl gauge is characterized by choosing A 0 = 0, i.e. the component of A in the time direction is zero.Applying it to (2) results in the vectorial wave equation.
Before we state the vectorial wave equation in terms of the Euclidean metric, we introduce some notation that is used in this paper.First, for d ∈ {2, 3}, the spatial bounded Lipschitz domain Ω ⊂ R d with boundary ∂Ω has the outward unit normal n x : ∂Ω → R d .Second, for a given terminal time T > 0, we define the space-time cylinder Q := (0, T ) × Ω ⊂ R d+1 and its lateral boundary Σ := [0, T ] × ∂Ω ⊂ R d+1 .Next, we fix some notation, which differs for d = 2 and d = 3.For this purpose, let f : Q → R d be a sufficiently smooth, vector-valued function with components f ι , ι = 1, . . ., d, and g : Q → R sufficiently smooth, scalar-valued function.We introduce the following notation: the spatial curl of f by the scalar function 1 as well as the spatial curl of g by the vector-valued function curl x g = (∂ x 2 g, −∂ x 1 g) .
• d = 3: In this case, the spatial curl of f is given by the vector-valued function where × denotes the usual cross product of vectors in R 3 .
With this notation and the Weyl gauge, we rewrite equation ( 2) into the vectorial wave equation for d = 2, 3 simply by inserting the definition of the exterior derivative in the Euclidean metric following [2, Chapter 2].Hence, we want to find a function A : where γ t is the tangential trace operator, j : Q → R d is a given current density, : Ω → R d×d is a given permittivity, and µ : Ω → R for d = 2 and µ : Ω → R 3×3 for d = 3 is a given permeability.However, equation ( 2) is a four-dimensional equation for a four-dimensional vector potential (A 0 , A) .Under the Weyl gauge, the fourth equation in ( 2) is div x ( (∂ t A)) = −ρ in Q.On the other hand, this equation holds as long as [16] for more details.So, in case of ∂ t A(0, •) = 0 in Ω we are limited to examples where ρ(0, •) = 0 in Ω.However, using the same technique used for inhomogeneous Dirichlet boundary conditions in finite element methods for the Poisson's equation, see e.g.[13], we can extend the results of this paper to inhomogeneous initial data.
With this knowledge, we give an outline of this paper.In Section 2, we recall well-known function spaces, which are needed to state the trial and test spaces of the variational formulations.In Section 3, we introduce the modified Hilbert transformation which allows for a Galerkin-Bubnov formulation.Then, we state different variational formulations and their properties in Section 4. Section 5 is devoted to the finite element spaces.These approximation spaces are used for the space-time continuous Galerkin finite element discretizations in Section 6, including investigations of a resulting CFL condition.Numerical examples for a two-dimensional spatial domain are presented in Section 7, which show the sharpness of the CFL condition and the different behavior of the numerical solutions when changing the projection of the right-hand side.Finally, we draw conclusions in Section 8.

Classical function spaces
For completeness, we recall classical function spaces, see [3,12,28,39] for further details and references.For this purpose, let a bounded Lipschitz domain D ⊂ R m for m ∈ N be fixed.We denote by L p (D), 1 ≤ p ≤ ∞, the usual Lebesgue space with its norm • L p (D) , and by H k (D), k ∈ N, the Sobolev space with Hilbertian norm • H k (D) .Further, the subspace , which is actually a norm in H 1 0 (D) due to the Poincaré inequality.For the interval D = (0, T ), we write L 2 (0, T ) := L 2 ((0, T )), H k (0, T ) := H k ((0, T )), and the subspaces , which again, is actually a norm in H 1 0, (0, T ) and H 1 ,0 (0, T ) due to Poincaré inequalities.By interpolation, we introduce H s 0, (0, T ) := [H 1 0, (0, T ), L 2 (0, T )] s and H s ,0 (0, Further, the aforementioned spaces are generalized from real-valued functions v : D → R to vector-valued functions v : D → X for a Hilbert space X with inner product (•, •) X .In particular, for X = R n with n ∈ N, the space L 2 (D; R n ) is the usual Lebesgue space for vector-valued functions v : D → R n endowed with the inner product (v, w) With this notation, for a bounded Lipschitz domain Ω ⊂ R d with d ∈ {2, 3}, we define the vector-valued Sobolev spaces , where div x is the (distributional) divergence.Similarly, with the (distributional) curl, we set . Next, we introduce the tangential trace operator γ t for d ∈ {2, 3}, see [3], [12,Section 4.3], [14,Theorem 2.11] and [28,Subsection 3.5.3]for details.For this purpose, we denote by (•) |∂Ω the usual trace of a function as a linear, surjective mapping from Here, H 1/2 (∂Ω; R m ), m ∈ N, is the fractional-order Sobolev space equipped with the Sobolev-Slobodeckij norm, see [12, Subsection 2.2.2], and H −1/2 (∂Ω; R m ) is its dual space.To define the tangential trace operator γ t , we distinguish two cases: • d = 2: For a function w ∈ H 1 (Ω; R 2 ), we define the pointwise trace γ t w = w |∂Ω • τ x , where τ x is the unit tangent vector satisfying τ x • n x = 0.Then, the continuous mapping γ t : H(curl; Ω) → H −1/2 (∂Ω; R) is the unique extension of γ t .

Modified Hilbert transformation
In this section, we introduce the modified Hilbert transformation H T as developed in [31,39], and state its main properties, see also [30,32,40,41].The modified Hilbert transformation acts in time only.Hence, in this section, we consider functions v : (0, T ) → R, where a generalization to functions in (t, x) is straightforward because of the tensor product structure of the domain Q.

Space-time variational formulations
In this section, we state space-time variational formulations of the vectorial wave equation (3).For this purpose, we recall a variational setting with different trial and test spaces as well as establish the tools to design a new space-time variational formulation of the vectorial wave equation ( 3) with equal trial and test spaces.The latter space-time variational formulation is derived by using the modified Hilbert transformation H T for the temporal part as introduced in Section 3.
To deploy the existence and uniqueness result of [16] to the stated variational setting, we make the following assumptions, which are assumed for the rest of this work.• Ω is a bounded Lipschitz domain, • and Q := (0, T ) × Ω is star-shaped with respect to a ball B, i.e. the convex hull of B and each point in Q is contained in Q.
ess inf Note that in Assumption 4.1, the given functions j, and µ are real-valued and that the functions , µ do not depend on the temporal variable t.Further, to handle the case that Q is not star-shaped with respect to a ball, we introduce the vectorial wave equation ( 3) in an extended space-time domain Q, which is star-shaped with respect to a ball and fulfills Q ⊂ Q.
Here, the functions j, and µ have to be extended appropriately by considering the material, e.g.assuming Q\Q is filled with air.

Different trial and test spaces
In this subsection, we derive a space-time variational formulation of the vectorial wave equation (3), where the trial and test spaces are different.With the aforementioned assumptions on the functions and µ the function spaces L 2 (Ω; R d ) and H 0 (curl; Ω) are endowed with the inner products respectively, where We denote the induced norms by as well as the norm Next, we define the space-time Sobolev spaces, which are used for the variational formulation.We consider which are endowed with the Hilbertian norm , where [39, Lemma 3.4.5] is applied.
With these space-time Sobolev spaces, we motivate a space-time variational formulation of the vectorial wave equation (3).For this purpose, we multiply the vectorial wave equation ( 3) by a test function v, integrate over the space-time domain Q and then use integration by parts with respect to space and time, resulting in the space-time variational formulation: for all v ∈ H curl;1 0;,0 (Q).Note that the first initial condition A(0, •) = 0 is fulfilled in the strong sense, whereas the second initial condition ∂ t A(0, •) = 0 is incorporated in a weak sense in the right side of (10).On the other hand, the boundary condition γ t A = 0 is satisfied in the strong sense with A(t, •) ∈ H 0 (curl; Ω) for almost all t ∈ (0, T ).
The unique solvability of the variational formulation ( 10) is proven in [17], [15, Theorem 3.8], and in [16] as a special case of Theorem 2, which is summarized in the next theorem.Theorem 4.2.Let Assumption 4.1 be satisfied.Then, a unique solution A ∈ H curl;1 0;0, (Q) of the variational formulation (10) exists and the stability estimate Note that the trial and test spaces of the variational formulation (10) are different.To get equal trial and test spaces, the modified Hilbert transformation H T of Section 3 is used, which is investigated in the next subsection.

Equal trial and test spaces
The goal of this subsection is to state a space-time variational formulation of the vectorial wave equation ( 3) with equal trial and test spaces.For this purpose, we extend the definition of the modified Hilbert transformation H T of Section 3 to the vector-valued functions in H curl,1 0;0, (Q) and H curl,1 0;,0 (Q), where the extension is denoted again by H T .More precisely, we define the mapping where the given Here, the set • an orthogonal basis of H 0 (curl; Ω) with respect to (•, •) H 0, ,µ (curl;Ω) , and where we refer to [16, Subsection 1.1.1]and [3, Subsection 8.3.1.2]for the construction of such a basis.The inverse operator where the given function Next, we prove properties of this definition of the modified Hilbert transformation H T .These properties are then used for the derivation of a Galerkin-Bubnov method for the spacetime variational formulation of the vectorial wave equation (3).Recall the definition of the vector-valued space-time Sobolev spaces H curl;1 0;0, (Q), H curl;1 0;,0 (Q) defined in ( 8), (9).
As H T acts only with respect to time, the properties ( 5), (7) of Section 3 remain valid.For completeness, we prove the following lemma.
for almost all x ∈ Ω.Then we derive for almost all (t, x) ∈ Q, and hence, the assertion.

Finite element spaces
As both variational formulations, the Galerkin-Petrov formulation (10) and the Galerkin-Bubnov formulation (17), are uniquely solvable, we investigate the discrete counterparts of these formulations.For that purpose, we first introduce the finite element spaces, which are used for discretizations of the variational formulations (10) and (17).For comparison reasons, we aim at a similar notation as in [16,Subsection 3.1], where, however, piecewise quadratic instead of piecewise linear functions are used for the temporal discretizations.
Let us start with the spatial discretization.In this section, we assume that the bounded Lipschitz domain Ω ⊂ R d is polygonal for d = 2, or polyhedral for d = 3.Let ν ∈ N 0 be the refinement level and (T x ν ) ν∈N 0 the corresponding mesh sequence of admissible decompositions T x ν of the spatial domain Ω.On each level ν we decompose the spatial domain Note that also other types of elements ω ι , e.g.rectangles (d = 2), are possible.We emphasize that we define the local mesh sizes by h x,ι = ωι dx 1/d , ι = 1, . . ., N x , since this particular choice occurs in the context of the CFL condition of the proposed Galerkin-Petrov method.Further, we define h x := h x,max (T x ν ) := max ι=1,...,N x h x,ι , and h x,min (T x ν ) := min ι=1,...,N x h x,ι .
For the investigations of CFL conditions, we consider for the rest of the paper a sequence (T x ν ) ν∈N 0 , satisfying a shape-regularity and a global quasi-uniformity.As these terms are not used in a uniform manner, we recall their definition.The mesh sequence (T x ν ) ν∈N 0 is Here, r ω > 0 is the inradius of the element ω, and x − y 2 is the Euclidean distance of the points x, y in R d .The constant c F influences the conditional stability, namely a CFL condition.To derive this CFL condition, we also need an inverse inequality for the spatial curl operator.The inverse inequality on the other hand is affected by the global quasi-uniformity for finite elements in H(curl; Ω).
With these conditions on the spatial mesh in mind, we introduce three spatial approximation spaces of vector-valued type, see [12,Section 19.2], [25, Section 3.5.1]or [28, Section 5.5] for more details.For this purpose, for p ∈ N 0 , n ∈ {1, 2, 3} and K ⊂ R n , the polynomials of total degree p are denoted by P p n (K).First, we define the space of vector-valued piecewise constant functions where the N 0 x basis functions ψ 0 are the componentwise characteristic functions with respect to the spatial elements.Second, we set the lowest-order Raviart-Thomas finite element space by where the N RT x basis functions ψ RT are attached to the edges for d = 2 and the faces for d = 3 of the spatial mesh T x ν and is the local polynomial space for ω ∈ T x ν , see [12,Section 14.1].Last, we introduce the lowest-order Nédélec finite element space of the first kind by and its subspace having zero tangential trace by where the N N x basis functions ψ N are attached to the edges of the spatial mesh T x ν , and the local polynomial spaces for ω ∈ T x ν are given by and Next, we investigate temporal finite element spaces.For this purpose, we consider a sequence of meshes {T t α } α∈N 0 in time.Here, α ∈ N 0 is the level of refinement and is the mesh defined by the decomposition T of the time interval (0, T ), where N t α =: N t is the number of temporal elements τ l = (t l−1 , t l ) ⊂ R, l = 1, . . ., N t , with mesh sizes h t,l = t l − t l−1 , l = 1, . . ., N t , and h t := max l=1,...,N t h t,l .First, the space of piecewise constant functions in time is given by with the elementwise characteristic functions ϕ 0 l as basis functions.Second, we write the space of piecewise linear, globally continuous functions as where the hat function ϕ 1 l is related to t l , i.e. ϕ 1 l (t k ) = 0 for k = l and ϕ 1 l (t l ) = 1.As we need also approximation spaces that incorporate initial or terminal conditions, we define the subspace S 1 0, (T t α ) := S 1 (T t α ) ∩ H 1 0, (0, T ) = span{ϕ 1 l } N t l=1 , fulfilling the homogeneous initial condition, and the subspace ,0 (T t α ) := S 1 (T t α ) ∩ H 1 ,0 (0, T ) = span{ϕ 1 l } N t −1 l=0 , satisfying the homogeneous terminal condition.
In other words, the space-time cylinder Last, to approximate the right-hand side of the vectorial wave equation (3), we consider two different L 2 (Q) projections of j ∈ L 2 (Q; R d ).In Section 7 on the numerical results we investigate the influence of the two types of projection on the convergence rates.Let us take a look at the definition of the projections for a given g ∈ L 2 (Q; R d ).First, the projection onto the space of piecewise constant functions Π 0 for all w h ∈ S 0 (T t α ) ⊗ S 0 d (T x ν ).Second, the projection onto the space of piecewise linear functions in time and Raviart-Thomas functions in space Π RT ,1 h for all w h ∈ S 1 (T t α ) ⊗ RT 0 (T x ν ).Note that both L 2 (Q) projections possibly have inhomogeneous Dirichlet, initial or terminal conditions.

FEM for the vectorial wave equation
Following the introduction to the finite element spaces in the last section, we derive conforming space-time discretizations for the variational formulations (10), (17), using the notation of Section 5.For this purpose, let the bounded Lipschitz domain Ω ⊂ R d be polygonal for d = 2, or polyhedral for d = 3.Additionally, we assume Assumption 4.1 and the assumptions on the mesh stated in Section 5.

Galerkin-Petrov FEM
In this subsection, we state the conforming discretization of the variational formulation (10), using the tensor-product spaces ( 18), (19), by a Galerkin-Petrov finite element method: for all w h ∈ S 1 ,0 (T t α ) ⊗ N 0 I,0 (T x ν ).Hence, we approximate the solution A ∈ H curl;1 0;0, (Q) of the variational formulation (10) by with coefficients A k κ ∈ R. The operator Π h in ( 22) is either the L 2 (Q) projection Π 0 h onto the space of piecewise constant functions defined in (20) or the L 2 (Q) projection Π RT ,1 h onto piecewise linear functions in time and Raviart-Thomas functions in space given in (21).The reason to replace the right-hand side j with Π h j in ( 22) is the comparison with the Galerkin-Bubnov finite element method (40), using the modified Hilbert transformation H T .An alternative is the use of formulas for numerical integration applied to (j, w h ) L 2 (Q) .Note that such formulas for numerical integration have to be sufficiently accurate to preserve the convergence rates of the finite element method, see the numerical examples in Section 7.
Next, we use the representation (23) to rewrite the discrete variational formulation (22) into the equivalent linear system Here, we order the degrees of freedom first in time then in space so that the coefficient vector in ( 23) can be written as where Further, for the linear system (24), we define the spatial matrices by for , κ = 1, . . ., N N x and the temporal matrices by for l = 0, . . ., N t − 1, k = 1, . . ., N t as well as the right-hand side by x for l = 0, . . ., N t − 1 and For the stability analysis of the linear system ( 24), we take a closer look at the system matrix.The system matrix of the linear system ( 24) is sparse and allows for a realization as a two-step method, due to the sparsity pattern of the temporal matrices A 1 ht , M 1 ht in (27).Further, this sparsity pattern forms the basis of classical stability analysis of the Galerkin-Petrov finite element method (22) in the framework of time-stepping or finite difference methods.See [31,Section 5] for the scalar wave equation, where the ideas can be extended to the vectorial wave equation, see the discussion in [15, Section 3.5.2]or [17].Here, we skip details and only state that the Galerkin-Petrov finite element method ( 22) is stable if the temporal mesh is uniform with mesh size h t and the CFL condition is satisfied with the constant c I > 0 of the spatial inverse inequality Bounds for the constant c I > 0 are given in [15,Lemma A.2] or [16,17].
As an example domain, which is also used in Section 7, we consider the unit square Ω = (0, 1) × (0, 1) ⊂ R 2 .Using uniform triangulations with isosceles right triangles as in Figure 1 yields the constant c I = 18.This result is proven by adapting the arguments of the proof of (30), given in [15,Lemma A.2] or [17], of the general situation to the case of isosceles right triangles.More precisely, in the proof of [15,Lemma A.2], the matrix J l J l is diagonal for isosceles right triangles, where J l denotes the Jacobian matrix of the geometric mapping from the reference element to the physical one.Thus, the eigenvalues of J l J l are known explicitly.So, in this situation, the CFL condition (29) reads as Numerical examples indicate the sharpness of the CFL condition (31), see Subsection 7.1.6.1.1Using the L 2 (Q) projection Π h = Π 0 h for the right-hand side J We present the calculation of the right-hand side J of the linear system (24) when Π h is the Then the projection Π 0 h j is the solution of the linear system with the matrices and vectors Thus, using this representation for the right-hand side (28) of the space-time Galerkin-Petrov method (24) yields for = 1, . . ., N N x , l = 0, . . ., N t − 1 with the matrix where for the right-hand side J Analogously to Subsection 6.1.1,we present the calculation of the right-hand side J of the linear system ( 24) when Π h is the L 2 (Q) projection onto S 1 (T t α ) ⊗ RT 0 (T x ν ) given in (21).The projection Π RT ,1 h j is the solution of the linear system with the matrices and vectors by the representation Thus, using this representation for the right-hand side (28) yields for = 1, . . ., N N x , l = 0, . . ., N t − 1 with the matrix where Note that the matrices 39) are submatrices of the matrix M 1 ht ∈ R (N t +1)×(N t +1) in (36).

Galerkin-Bubnov FEM
In this subsection, we discretize the variational formulation (17) where we apply the modified Hilbert transformation.We use a tensor-product ansatz and the conforming finite element space (18).In particular, we consider the Galerkin-Bubnov finite element method to find for all given in (21).The discrete variational formulation ( 40) is equivalent to the linear system with the spatial matrices A N hx , M N hx given in ( 26) and the temporal matrices for l, k = 1, . . ., N t , where H T , defined in Section 3, acts solely on time-dependent functions.
As in Subsection 6.1, we use again the representation (23) of A h and the vector A ∈ R N t N N x of its coefficients (25).The right-hand side of the linear system ( 41) is given by where Let us take a look at the solvability of the system matrix of the linear system (41).The temporal matrices A H T ht , M H T ht in (42) are positive definite due to property (6).In addition, the spatial matrix M N hx in ( 26) is also positive definite, whereas the spatial matrix A N hx in ( 26) is only positive semi-definite.Hence, the Kronecker product A H T ht ⊗ M N hx is positive definite.On the other hand, the product M H T ht ⊗ A N hx is positive semi-definite.Adding both results in a positive definite system matrix of the linear system (41).In other words, the linear system ( 41) is uniquely solvable.Note that, compared to the static case, we do not need any stabilization to get unique solvability of the linear system (41) and hence of the corresponding discrete variational formulation (40).However, further details on the numerical analysis of the Galerkin-Bubnov finite element method (40) are far beyond the scope of this contribution, we refer to [26,27] for the case of the scalar wave equation.
In addition, note that the temporal matrices A H T ht , M H T ht in (42) are dense, whereas the spatial matrices A N hx , M N hx in (26) are sparse.Thus, the linear system (41) does not allow for a realization as a multistep method.However, the application of fast (direct) solvers, as known for heat and scalar wave equations [23,37,38], is possible, which is the topic of future work.

Using the
h for the right-hand side J H T In this subsection, we present the calculation of the right-hand side J H T of the linear system (41) given in (20).Since this subsection is similar to Subsection 6.1.1,we skip the details.
Then the entries (43) of the right-hand side J H T admit the representation for = 1, . . ., N N x , l = 1, . . ., N t with the matrix where M N ,0 hx is defined in (33), J is given in (34) and for the right-hand side J H T Analogously to Subsection 6.2.1, we present the calculation of the right-hand side J H T of the linear system (41 given in (21).As for Subsection 6.2.1, we skip the details, which are analogous to Subsection 6.1.2.
Then the entries (43) of the right-hand side J H T admit the representation for = 1, . . ., N N x , l = 1, . . ., N t with the matrix where M N ,RT hx is defined in (37), J is given in (38) and Note that the matrix

Numerical examples for the vectorial wave equation
In this section, we give numerical examples of the conforming space-time finite element methods ( 22) and ( 40) for a spatially two-dimensional domain Ω ⊂ R 2 , i.e. d = 2.For this purpose, we consider the unit square Ω = (0, 1)×(0, 1), and set (x) = 1 0 0 1 , µ(x) = 1 for x ∈ Ω.The spatial meshes T x ν are given by uniform decompositions of the spatial domain Ω into isosceles right triangles, where a uniform refinement strategy is applied, see Figure 2. We investigate the terminal times T ∈ √ 2, 3 2 to examine the CFL condition (31) of the Galerkin-Petrov finite element method (22).The temporal meshes T t α are defined by t l = T l N t α for l = 0, . . ., N t α , where N t α = 5 • 2 α , α = 0, . . ., 4. Note that for this choice of spatial and temporal meshes, we are in the framework of the CFL condition (31).In the following numerical examples, we measure the error of the space-time finite element methods ( 22) and (40) in the space-time norms In particular, we state the numerical results for A − A h L 2 (Q) and , where A ∈ H curl;1 0;0, (Q) is the solution of the variational formulation (10), and is the solution of space-time finite element method (22) or (40).For the vectorial wave equation (3) and its variational formulation (10), we use the manufactured solution for (t, x 1 , x 2 ) ∈ Q, which is, for comparision, also investigated in [16].This solution satisfies the homogeneous Dirichlet condition γ t A(t, x 1 , x 2 ) = A |Σ (t, x 1 , x 2 )×n x (x 1 , x 2 ) = 0 for (t, x 1 , x 2 ) ∈ Σ as well as the homogeneous initial conditions A(0, x 1 , x 2 ) = ∂ t A(0, x 1 , x 2 ) = 0 for (x 1 , x 2 ) ∈ Ω.The related right-hand side j is given by For the computations, we use high-order quadrature rules to calculate the integrals for computing the projections Π 0 h j, Π RT ,1 h j in (32), (35).The temporal matrices involving the modified Hilbert transformation H T , e.g. the matrices (42), are assembled as proposed in [40, Subsection 2.2], see also [41] for further assembling strategies.All spatial matrices, e.g. the matrices (26), are calculated with the help of the finite element library Netgen/NGSolve, see www.ngsolve.organd [29].We solve the linear systems using the sparse direct solver UMFPACK 5.7.1 [9] in the standard configuration.All calculations, presented in this section, were performed on a PC with two Intel Xeon E5-2687W v4 CPUs 3.00 GHz, i.e. in sum 24 cores and 512 GB main memory.
Last, in

Galerkin-Petrov FEM
In this subsection, we investigate numerical examples for the Galerkin-Petrov finite element method (22) in the situation described at the beginning of this section.
, the unit square Ω and A in (45).
right-hand side Π 0 h j ∈ S 0 (T t α ) ⊗ S 0 2 (T x ν ) of Subsection 6.1.1.In Table 5 and Table 6, we report the results for Π RT ,1 h j ∈ S 1 (T t α ) ⊗ RT 0 (T x ν ) of Subsection 6.1.2.All tables show conditional stability, i.e. the CFL condition (31) is required for stability.Note that the ratio of the mesh sizes h t = 0.0177 and h x = 0.0221, i.e. the last column and second last row of Tables 3 to 6, is given by h t h x ≈ 0.0177 0.0221 ≈ 0.801 and thus, fulfills the CFL condition (31) resulting in a stable method.In the case of stability, we observe first-order convergence in |•| H curl;1 (Q) , where the errors in Table 1, Table 3 and Table 5 are within the same range.When considering the errors in • L 2 (Q) , Table 4 reports only first-order convergence, whereas in Table 6, second-order convergence is observed as in Table 2 on the interpolation error.Next, we elaborate on this convergence behavior by investigating the difference |A(T, x) − A h (T, x)| for x ∈ Ω.For the mesh sizes h t = 0.0707, h x = 0.0884, we plot the function Ω x → |A(T, x) − A h (T, x)| in Figure 3 for both projections Π 0 h and Π RT ,1 h .In the literature, we find the term spurious solutions, e.g. in [21], or spurious modes, e.g. in [5].Spurious modes are parts of the numerical solution, which correspond to eigensolutions of the discrete differential operator.They are oscillations that should not be part of the computed solution.Spurious modes can be understood as numerically generated noise and have no physical meaning.An example can be found in [21,Figure 5.8].In Figure 3, we obtain a similar noise behavior that occurs for the projection Π 0 h j of the right-hand side.Here, we observe some of the zero eigensolutions of the curl-curl operator added to the solution.Hence, we get a worse L 2 (Q)error, see Table 4 and Table 6, but a similar error behavior in the |•| H curl;1 (Q) -norm for both projections if the CFL condition is met, see Table 3 and Table 5.
To summarize, the approximation Π 0 h j ∈ S 0 (T t α ) ⊗ S 0 2 (T x ν ) of j is good enough, when the Second, we examine the terminal time T = 3 2 .In Table 7, Table 8, the numerical results for the Galerkin-Petrov finite element method (22) with the approximate right-hand side Π RT ,1 h j ∈ S 1 (T t α ) ⊗ RT 0 (T x ν ) of Subsection 6.1.2show that slightly violating the CFL condition (31) leads to instability.More precisely, the ratio of the mesh sizes h t = 0.0188 and h x = 0.0221, i.e. the last column and second last row of

Galerkin-Bubnov FEM
In this subsection, we report on numerical results for the Galerkin-Bubnov finite element method (40) using the modified Hilbert transformation in the situation described at the be-  We observe unconditional stability in Table 9, Table 10, Table 11, Table 12, i.e. no CFL condition is needed.This is the main difference between the results in this subsection and the results of Subsection 7.1, where for stability, a CFL condition is required.
Besides the stability issue, the errors in Table 9, Table 10, Table 11, Table 12 are comparable with the errors of the previous Subsection 7.1 regarding the projection of the right-hand side.Thus, the approximation Π 0 h j ∈ S 0 (T t α ) ⊗ S 0 2 (T x ν ) of j is good enough to get first-order convergence in |•| H curl;1 (Q) , see Table 9.In Table 12, we again see that the better approximation Π RT ,1 h j ∈ S 1 (T t α ) ⊗ RT 0 (T x ν ) of j is needed for optimal convergence rates in • L 2 (Q) .

Conclusion
In this paper, we presented two conforming space-time finite element methods for the vectorial wave equation.First, we stated a variational formulation with different trial and test spaces.Its conforming discretization with piecewise multilinear functions leads to a Galerkin-Petrov method, which is only conditionally stable, i.e. a CFL condition is required.For a particular choice of spatial meshes, we stated the CFL condition, where numerical examples showed its sharpness.Second, to tackle the problem of a CFL condition, we introduced a variational formulation for the vectorial wave equation with equal trial and test spaces using the modified Hilbert transformation H T .We gave a rigorous derivation of this variational setting.A conforming discretization with piecewise multilinear functions of this new variational approach results in a space-time Galerkin-Bubnov method.All numerical examples showed the unconditional stability of this conforming space-time method.Further, we investigated the influence of projections of the right-hand side on the convergence.In numerical examples for both the Galerkin-Petrov method and the Galerkin-Bubnov method, we observed only first-order convergence in • L 2 (Q) when the right-hand side was projected onto piecewise constants.On the other side, second-order convergence was obtained when the right-hand side was projected onto piecewise linear functions in time and lowest-order Raviart-Thomas functions in space.
For both presented conforming space-time approaches, a generalization to piecewise polynomials of higher-order is possible, see [16].Further, the fast solvers [23,37,38] for the resulting linear systems, where some allow for time parallelization, are also applicable.In addition, other time-dependent problems in electromagnetics can be handled with the presented variational frameworks, which are topics of future work, see e.g.[16].

Figure 1 :
Figure 1: Uniform triangulations of the unit square with isosceles right triangles.

Figure 2 :
Figure 2: Spatial meshes T x ν : Initial mesh T x 0 and the mesh T x 2 after two uniform refinements.

Table 1 :
1, Table 2, we report the interpolation error in the norms |•| H curl;1 (Q) and • L 2 (Q) for the unit square Ω and T = √ 2 for the function A defined in (45), where first-order convergence is observed in |•| H curl;1 (Q) and second-order convergence is obtained in • L 2 (Q) .Interpolation errors in |•| H curl;1 (Q) for T = √2, the unit square Ω and A in (45).

Table 4 :
(22)rs in • L 2 (Q) of the Galerkin-Petrov FEM(22)with the approximate right-hand In Table3and Table 4, we present the numerical results for the Galerkin-Petrov finite element method(22)with the approximate

Table 7 ,
Table 8, is given by

Table 8 :
(22)rs in • L 2 (Q) of the Galerkin-Petrov FEM(22)with the approximate right-handside Π RT ,1 h j ∈ S 1 (T t α ) ⊗ RT 0 (T x ν ) for T = 3 2, the unit square Ω and A in (45).ginning of this section.We only show numerical results for the terminal time T = √ 2, as those for T = 3 2 are similar.