Stable and convergent fully discrete interior–exterior coupling of Maxwell’s equations
Abstract
Maxwell’s equations are considered with transparent boundary conditions, for initial conditions and inhomogeneity having support in a bounded, not necessarily convex threedimensional domain or in a collection of such domains. The numerical method only involves the interior domain and its boundary. The transparent boundary conditions are imposed via a timedependent boundary integral operator that is shown to satisfy a coercivity property. The stability of the numerical method relies on this coercivity and on an antisymmetric structure of the discretized equations that is inherited from a weak firstorder formulation of the continuous equations. The method proposed here uses a discontinuous Galerkin method and the leapfrog scheme in the interior and is coupled to boundary elements and convolution quadrature on the boundary. The method is explicit in the interior and implicit on the boundary. Stability and convergence of the spatial semidiscretization are proven, and with a computationally simple stabilization term, this is also shown for the full discretization.
Keywords:
transparent boundary conditions Calderon operator discontinuous Galerkin boundary elements leapfrog scheme convolution quadratureMsc:
35Q61 65M60 65M38 65M12 65R20∎
1 Introduction
Maxwell’s equations on the whole threedimensional space are considered with initial conditions and inhomogeneity having support in a bounded domain that is not required to be convex (or in a finite collection of such domains). The study of such problems leads to transparent boundary conditions, which yield the restriction of the solution to the domain. Such boundary conditions are nonlocal in space and time, for both acoustic wave equations and Maxwell’s equations. There is a vast literature to tackle this problem in general for wave equations: fast algorithms for exact, nonlocal boundary conditions on a ball GroteKeller ; Hagstrom , local absorbing boundary conditions EngquistMajda ; HagstromMarOrGivoli , perfectly matched layers, which were originally considered for electromagnetism in Berenger , and numerical coupling with timedependent boundary integral operators abboud2011coupling ; leapfrog . All the above approaches, except the last one, are inadequate for nonconvex domains. The local methods fail because waves may leave and reenter a nonconvex domain. Inclusion of a nonconvex domain in a larger convex domain is computationally undesirable in situations such as a cavity or an antennalike structure or a farspread nonconnected collection of small domains.
The main objective of the present work is to transfer the programme of leapfrog from acoustic wave equations to Maxwell’s equations: to propose and analyze a provably stable and convergent fully discrete numerical method that couples discretizations in the interior and on the boundary, without requiring convexity of the domain. Like Abboud et al. abboud2011coupling (and later also leapfrog ) for the acoustic wave equation, we start from a symmetrized weak firstorder formulation of Maxwell’s equations. In the interior this is discretized by a discontinuous Galerkin (dG) method in space DiPietroErn ; HesthavenWarburton ; HochbruckSturm together with the explicit leapfrog scheme in time HairerLubichWanner . The boundary integral terms are discretized by standard boundary element methods in space and by convolution quadrature (CQ) in time LubichCQ ; Lubichmultistep . This yields a coupled method that is explicit in the interior and implicit on the boundary. The choice of a CQ time discretization of the boundary integral operators is essential for our analysis, and to a lesser extent also the choice of the leapfrog scheme in the interior. However, our approach is not specific to the chosen space discretizations which could, in particular, be replaced by conformal edge elements Hiptmairacta .
While the general approach of this paper is clearly based on leapfrog , it should be emphasized that the appropriate boundary integral formulation requires a careful study of the timeharmonic Maxwell’s equation. This is based on BuffaCostabelSheen_traces ; BuffaHiptmair ; BuffaHiptmairvonPetersdorffSchwab ; Costabel2004 ; BallaniBanjaiSauterVeit , with special attention to the appropriate trace space on the boundary and to the corresponding duality. Due to the analogue of Green’s formula for Maxwell’s equations, the duality naturally turns out to be an antisymmetric pairing.
The Calderon operator for Maxwell’s equation, which arises in the boundary integral equation formulation of the transparent boundary conditions, differs from the acoustic case to a large extent, and therefore the study of its coercivity property is an important and nontrivial point. Similarly to the acoustic case, the continuoustime and discretetime coercivity is obtained from the Laplacedomain coercivity using an operatorvalued version, given in leapfrog , of the classical Herglotz theorem Herglotz . Both the second and first order formulation of Maxwell’s equations are used.
The spatial semidiscretization of the symmetrized weak firstorder formulation of Maxwell’s equations has formally the same matrix–vector formulation as for the acoustic wave equation studied in leapfrog , with the same coercivity property of the Calderon operator. Because of this structural similarity, the stability results of leapfrog , which are shown using the matrix–vector setting, remain valid for the Maxwell case without any modification. On the other hand, their translation to the functional analytic setting differs to a great extent. Therefore further care is required in the consistency analysis.
In Section 2 we recapitulate the basic theory for Maxwell’s equation in the Laplace domain. Based on Buffa and Hiptmair BuffaHiptmair , and further on BuffaHiptmairvonPetersdorffSchwab ; BallaniBanjaiSauterVeit , we describe the right boundary space, which allows for a rigorous boundary integral formulation for Maxwell’s equations. Then the boundary integral operators are obtained in a usual way from the single and double layer potentials.
In Section 3 we prove the crucial technical result of the present work, a coercivity property of the Calderon operator for Maxwell’s equation in the Laplace domain. This property translates to the continuoustime Maxwell’s equations later, in Section 4.2, via an operatorvalued Herglotz theorem.
In Section 4 we study the interior–exterior coupling of Maxwell’s equations, resulting in an interior problem coupled to an equation on the boundary with the Calderon operator. We derive a first order symmetric weak formulation, which is the Maxwell analogue of the formulation of abboud2011coupling for the acoustic wave equation. Together with the continuoustime version of the coercivity property of the Calderon operator, this formulation allows us to derive an energy estimate. Later on this analysis is transfered to the semidiscrete and fully discrete settings.
Section 5 presents the details of the discretization methods: In space we use discontinuous Galerkin finite elements with centered fluxes in the domain DiPietroErn ; HesthavenWarburton , coupled to continuous linear boundary elements on the surface. Time discretization is done by the leapfrog scheme in the interior domain, while on the boundary we use convolution quadrature based on the secondorder backward differentiation formula. An extra term stabilizes the coupling, just as for the acoustic wave equation leapfrog . The matrix–vector formulation of the semidiscrete problem has the same antisymmetric structure and the same coercivity property as for the acoustic wave equation, and therefore the stability results shown in leapfrog can be reused here.
In Sections 6 and 7 we revise the parts of the results and proofs of leapfrog where they differ from the acoustic case, which is mainly in the estimate of the consistency error. Finally, we arrive at the convergence error bounds for the semi and full discretizations.
To our knowledge, the proposed numerical discretizations in this paper are the first provably stable and convergent semi and full discretizations to interior–exterior coupling of Maxwell’s equations. We believe that the presented analysis and the techniques, which we share with leapfrog , can be extended further: to other discretization techniques for the domain, such as edge element methods Hiptmairacta , higher order discontinuous Galerkin methods, and different time discretizations in the domain, together with higher order Runge–Kutta based convolution quadratures on the boundary BanjaiLubichMelenk .
For ease of presentation we consider only constant permeability and permittivity. However, it is only important that the permeability and permittivity are constant in the exterior domain and in a neighbourhood of the boundary. In the interior these coefficients may be spacedependent and discontinuous. In the latter case the equations can be discretized in space with the dG method as described in HochbruckSturm .
In this paper we focus on the appropriate boundary integral formulation and on the numerical analysis of the proposed numerical methods. Numerical experiments are intended to be presented in subsequent work.
Concerning notation, we use the convention that vectors in are denoted by italic letters (such as ), whereas the corresponding boldface letters are used for finite element nodal vectors in , where is the (large) number of discretization nodes. Hence, any boldface letters appearing in this paper refer to the matrix–vector formulation of spatially discretized equations. Functions defined in the domain are denoted by letters from the Roman alphabet, while functions defined on the boundary are denoted by Greek letters.
2 Recap: the timeharmonic Maxwell’s equation and its boundary integral operators
2.1 Preliminaries and notation
Let us consider the timeharmonic Maxwell’s equation, obtained as the Laplace transform of the second order Maxwell’s equation (with constant permeability and permittivity ):
(1) 
where is the boundary of a bounded piecewise smooth domain (or a finite collection of such domains) , not necessarily convex, with exterior normal .
We shortly recall some useful concepts and formulas regarding the above problem, based on BuffaHiptmair and KirschHettlich . For the usual trace we will use the notation . The tangential and magnetic traces are defined, respectively, as
These traces are also often called Dirichlet trace and Neumann trace, motivated by the analogue of Green’s formula for Maxwell’s equations (for sufficiently regular functions):
(2)  
We introduce an important notation, the
which appears on the righthand side of (2). We note that the relation holds, cf. BuffaHiptmair ; KirschHettlich .
2.2 Function spaces
We collect some results on function spaces, which will play an important role later on. All of the results in the present subsection can be found in Section 2 of BuffaHiptmair .
Let us start by recalling the usual definition of the Sobolev space corresponding to the operator:
with corresponding norm .
Clearly, the above integral relations hold for functions .
Now we are turning to trace spaces. However, even though is a continuous mapping, is not the right choice for boundary integral operators. As it was emphasized by Buffa and Hiptmair BuffaHiptmair : the study of the continuous mapping is “actually sufficient for the understanding of homogeneous boundary conditions for fields in the Hilbert space context. However, to impose meaningful nonhomogeneous boundary conditions or, even more important, to lay the foundations for boundary integral equations we need to identify a proper trace space”^{1}^{1}1Quoted from Buffa and Hiptmair, BuffaHiptmair , Section 2.2.. In the following, we briefly summarize the definition of such a trace space, together with some related results.
The Hilbert space collects the traces of functions, for , i.e., . This space is equipped with an inner product such that is continuous. In particular, the space has the dual space , defined with respect to the (extended) duality .
Then, the above mentioned proper trace space is given as:
with norm
(4) 
The tangential trace satisfies the following analogue of the trace theorem.
Lemma 1 (BuffaCostabelSheen_traces , Theorem 4.1)
The trace operator is continuous and surjective.
The following lemma clarifies the role of the antisymmetric pairing .
Lemma 2 (BuffaCostabelSheen_traces , Lemma 5.6 and BuffaHiptmair , Theorem 2)
The pairing can be extended to a continuous bilinear form on . With this pairing the space becomes its own dual.
The above results clearly point out that a natural choice of trace space is , which fits perfectly to the analogue of Green’s formula (2) and to the boundary integral formulation of Maxwell’s equations. This trace space is appropriate for the analysis of boundary integral operators.
2.3 Boundary integral operators
On potentials and boundary integral operators we follow Buffa and Hiptmair BuffaHiptmair , and we also refer to BuffaHiptmairvonPetersdorffSchwab ; Costabel2004 .
The usual boundary integral potentials for the timeharmonic Maxwell’s equation
are obtained, based on BuffaHiptmair and BallaniBanjaiSauterVeit : the (electric) single layer potential is given, for , as
while the (electric) double layer potential is given, for , as
where the fundamental solution is given, for , as
The solution then has the representation
(5) 
where
(6) 
Here denotes the jumps in the boundary traces. A further notation is the average of the inner and outer traces on the boundary: . On vectors both operations are acting componentwise.
For every and , formula (5) defines . Because of the jump relations
and are reconstructed from by (6).
Let us now define the boundary integral operators. As opposed to the general second order elliptic case, due to additional symmetries of the problem, they reduce to two operators and , see (BuffaHiptmair, , Section 5). They satisfy
In (BuffaHiptmair, , Section 5) the continuity of these operators was proven, without giving an explicit dependence on . Such bounds are crucial in the analysis later, therefore we now show explicit estimates for the boundary integral operators. Our result is based on BallaniBanjaiSauterVeit .
Lemma 3
For the boundary integral operators are bounded as
Proof
These estimates can be shown by adapting the arguments of (BallaniBanjaiSauterVeit, , Section 4.2). In particular, by using the antisymmetric pairing instead of the usual inner product, the results of (BallaniBanjaiSauterVeit, , Theorem 4.4) transfer from to the estimates stated here.
Furthermore, using the potential representation of the solution (5), the averages of the traces can be expressed using the operators and in the following way:
(7)  
3 Coercivity of a Calderon operator for the timeharmonic Maxwell’s equation
An important role will be played by the following operator on , to which we refer as a Calderon operator:
(8) 
The extra factor appears unmotivated here, but will turn out to be convenient later. This operator satisfies the following coercivity result, which is the key lemma of this paper.
Lemma 4
Proof
The proof has a structure similar to the proof of the corresponding result for the acoustic Helmholtz equation (leapfrog, , Lemma 3.1), although it now uses a different functionalanalytic setting. The structural similarity becomes possible thanks to the antisymmetric duality pairing that replaces the symmetric duality pairing of the acoustic case.
For given , we define by the representation formula (5). We can then express and in terms of by (6). We note that (7) yields
Now, using the properties of the antisymmetric pairing (acting componentwise on ), the analogue of Green’s formula (3) and using the definition of the traces, we obtain
We further obtain
and for we use the fact that :
where, for the first inequalities in both estimates, we used the trace inequality of Lemma 1. Extraction of factors and dividing through completes the proof.
4 Boundary integral formulation of Maxwell’s equations
Let us consider the first order formulation of Maxwell’s equations, in the following form:
with appropriate initial and boundary conditions. If the initial conditions satisfy the last two equations, then they hold for all times, see Monk ; ChenMonk , therefore these conditions are assumed to hold. The permeability and permittivity is denoted by and , respectively, and they are assumed to be positive constants, while denotes the electric current density.
Using the relation , the above equation can be written as the second order problem
with .
Setting , applying Laplace transformation, and writing instead of , we obtain the timeharmonic version (1).
4.1 Recap: Temporal convolutions and Herglotz theorem
We recall an operatorvalued continuoustime Herglotz theorem from (leapfrog, , Section 2.2), which is crucial for transferring the coercivity result of Lemma 4 from the Maxwell’s equation in the Laplace domain to the timedependent Maxwell’s equation. We describe the result in an abstract Hilbert space setting.
Let be a complex Hilbert space, with dual and antiduality . Let and be both analytic families of bounded linear operators for , satisfying the uniform bounds:
For any integer , we define the integral kernel
For a function with vanishing initial data, , we let
that is, is the distributional convolution of the inverse Laplace transform of with .
The operatorvalued version of the classical Herglotz theorem from leapfrog yields the following result.
Lemma 5 (leapfrog , Lemma 2.2)
In the above setting, the following two statements are equivalent:

, for any , ;

, for all , with , and for all .
4.2 Calderon operator for Maxwell’s equations
Consider the second order formulation of Maxwell’s equations in three dimensions:
Let be a bounded Lipschitz domain, with boundary , and further assume that the initial values and are supported within .
We rewrite this problem as an interior problem over :
and as an exterior problem over :
The two problems are coupled by the transmission conditions:
4.3 First order formulation
From now on, we use Maxwell’s equations in their first order formulation on the interior domain (and we omit the omnipresent superscript ):
(9) 
with the coupling through the Calderon operator as
where and . In addition, by we obtain
(10)  
where we also used (6). Hence, and
4.4 Coercivity of the timedependent Calderon operator
In the same way as in (leapfrog, , Lemma 4.1) for the acoustic wave equation, the coercivity of the Calderon operator for the timeharmonic Maxwell’s equation as given by Lemma 4 together with the operatorvalued continuoustime Herglotz theorem as stated in Lemma 2.3 yields coercivity of the timedependent Calderon operator .
Lemma 6
With the constant from Lemma 4 we have that
for arbitrary and for all and all with and , and with constant .
A Gronwall argument then yields the following energy estimate; see (leapfrog, , Lemma 4.2).
Lemma 7
Let the functions , , and with , , be such that for all
Then, with ,
(11)  
4.5 Weak formulation and energy estimate
Analogously to abboud2011coupling ; leapfrog , a symmetric weak form of (9) is obtained on using
and using (10) for the boundary term. Here denotes the standard inner product.
The coupled weak problem then reads: find and such that
(12)  
hold for arbitrary , and .
While this weak formulation is apparently nonstandard for Maxwell’s equations, we will see that it is extremely useful, in the same way as the analogous formulation proved to be for the acoustic case in abboud2011coupling ; leapfrog .
5 Discretization
5.1 Space discretization: dG and BEM
For the spatial discretization we use, as an example, the central flux discontinuous Galerkin (dG) discretization from HochbruckSturm (see also DiPietroErn ; HesthavenWarburton ) in the interior and continuous linear boundary elements on the surface.
We triangulate the bounded polyhedral domain by simplicial triangulations , where denotes the maximal element diameter. For our theoretical results we consider a quasiuniform and contactregular family of such triangulations with ; see e.g. DiPietroErn for these notions. We adopt the following notation from (HochbruckSturm, , Section 2.3): The faces of , decomposed into boundary and interior faces: . The normal of an interior face is denoted by . It is kept fixed and is the outward normal of one of the two neighbouring mesh elements. We denote by that neighbouring element into which is directed. The outer faces of are used as the triangulation of the boundary .
The dG space of vector valued functions, which are elementwise linear in each component, is defined as
The boundary element space is taken as
The corresponding nodal basis functions are denoted by and , respectively. Jumps and averages over faces are denoted analogously as for trace operators on , see Section 2.3:
where is the usual trace onto the face . We often omit the subscript as it will always be clear from the context.
The discrete operator with centered fluxes was presented in (HochbruckSturm, , Section 2.3):
By the arguments of the proof of Lemma 2.2 in HochbruckSturm , we obtain that the discrete curl operator satisfies the discrete version of Green’s formula (2),
(13) 
The operator is well defined on , with the broken Sobolev space
which is a Hilbert space with natural norm and seminorm and , respectively.
Using the above discrete operator, the semidiscrete problem reads as follows: Find and such that for all and ,