Stable and convergent fully discrete interior–exterior coupling of Maxwell’s equations

Stable and convergent fully discrete interior–exterior coupling of Maxwell’s equations

Balázs Kovács B. Kovács Mathematisches Institut, University of Tübingen,
Auf der Morgenstelle 10, 72076 Tübingen, Germany
1Ch. Lubich Mathematisches Institut, University of Tübingen,
Auf der Morgenstelle 10, 72076 Tübingen, Germany
   Christian Lubich B. Kovács Mathematisches Institut, University of Tübingen,
Auf der Morgenstelle 10, 72076 Tübingen, Germany
1Ch. Lubich Mathematisches Institut, University of Tübingen,
Auf der Morgenstelle 10, 72076 Tübingen, Germany
b e
b e

Maxwell’s equations are considered with transparent boundary conditions, for initial conditions and inhomogeneity having support in a bounded, not necessarily convex three-dimensional 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 time-dependent boundary integral operator that is shown to satisfy a coercivity property. The stability of the numerical method relies on this coercivity and on an anti-symmetric structure of the discretized equations that is inherited from a weak first-order 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.

transparent boundary conditions Calderon operator discontinuous Galerkin boundary elements leapfrog scheme convolution quadrature
35Q61 65M60 65M38 65M12 65R20

1 Introduction

Maxwell’s equations on the whole three-dimensional 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 time-dependent boundary integral operators abboud2011coupling ; leapfrog . All the above approaches, except the last one, are inadequate for non-convex domains. The local methods fail because waves may leave and re-enter a non-convex domain. Inclusion of a non-convex domain in a larger convex domain is computationally undesirable in situations such as a cavity or an antenna-like structure or a far-spread non-connected 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 first-order 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 ; Lubich-multistep . 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 Hiptmair-acta .

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 time-harmonic 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 anti-symmetric 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 continuous-time and discrete-time coercivity is obtained from the Laplace-domain coercivity using an operator-valued 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 semi-discretization of the symmetrized weak first-order 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 continuous-time Maxwell’s equations later, in Section 4.2, via an operator-valued 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 continuous-time 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 semi-discrete 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 second-order 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 anti-symmetric 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 Hiptmair-acta , 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 space-dependent 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 time-harmonic Maxwell’s equation and its boundary integral operators

2.1 Preliminaries and notation

Let us consider the time-harmonic Maxwell’s equation, obtained as the Laplace transform of the second order Maxwell’s equation (with constant permeability and permittivity ):


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):


We introduce an important notation, the

which appears on the right-hand side of (2). We note that the relation holds, cf. BuffaHiptmair ; KirschHettlich .

Let us now set , which provides

Moreover, if satisfies (1) and , then


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 non-homogeneous boundary conditions or, even more important, to lay the foundations for boundary integral equations we need to identify a proper trace space”111Quoted 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


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 anti-symmetric 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 time-harmonic 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




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


These estimates can be shown by adapting the arguments of (BallaniBanjaiSauterVeit, , Section 4.2). In particular, by using the anti-symmetric 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:


3 Coercivity of a Calderon operator for the time-harmonic Maxwell’s equation

An important role will be played by the following operator on , to which we refer as a Calderon operator:


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

There exists such that the Calderon operator (8) satisfies

for and for all , with


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 functional-analytic setting. The structural similarity becomes possible thanks to the anti-symmetric 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 anti-symmetric 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 time-harmonic version (1).

4.1 Recap: Temporal convolutions and Herglotz theorem

We recall an operator-valued continuous-time 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 time-dependent Maxwell’s equation. We describe the result in an abstract Hilbert space setting.

Let be a complex Hilbert space, with dual and anti-duality . 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 operator-valued 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:

  1. , for any , ;

  2. , 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:

Using the temporal convolution operators of Section 4.1, the solution of the exterior problem is given as

with boundary densities

which satisfy the equation

Here is the temporal convolution operator with the distribution whose Laplace transform is the Calderon operator defined in (8).

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 ):


with the coupling through the Calderon operator as

where and . In addition, by we obtain


where we also used (6). Hence, and

4.4 Coercivity of the time-dependent 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 time-harmonic Maxwell’s equation as given by Lemma 4 together with the operator-valued continuous-time Herglotz theorem as stated in Lemma 2.3 yields coercivity of the time-dependent 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 ,


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


hold for arbitrary , and .

While this weak formulation is apparently non-standard 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 .

Testing with , and , in (12), by using (10) we obtain

and summing up the three equations yield

For , the coercivity of the continuous-time Calderon operator, as stated in Lemmas 6 and 7, yields that the electromagnetic energy

satisfies the energy estimate (11) (with ) for arbitrary .

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 quasi-uniform and contact-regular 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),


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 ,