¹

1 Introduction↩︎

Graphene is a monolayer of carbon atoms arranged in a hexagonal lattice and has demonstrated exceptional thermal, electrical and optical properties, as well as structural robustness. This makes it an attractive candidate for reinforcement in materials, as well as for applications in organic electronics and optoelectronics. There is significant interest in both the scientific community and industrial sectors in graphene. In recent years, numerous other 2D materials have emerged, with the term referring to crystalline solids that have a reduced thickness, usually consisting of a single or only a few atomic layers. This unique characteristic gives them exceptional properties. Another important class of composite 2D materials are the semiconducting Transition Metal Dichalcogenides (TMDCs) such as or which also come in single layers, albeit with a more complicated unit cell than graphene. These materials have a wide range of applications, including catalysis, spintronics, and optoelectronics; see the reviews Cho.L.S.Et.2010?, Ted.L.O.2016?. Numerical simulations are vital for studying such materials.

The optical properties of such materials can be studied by depositing a sheet on a thin dielectric layer on top of a metal plate and exciting it with light pulses. The interaction between the pulses and the material is described by Maxwell equations coupled to quantum mechanical models, see, e.g., Cha.E.Z.Et.2008?, Kan.M.2005?, Sta.2014?. A simpler way of modeling the interaction of the 2D material in the Maxwell equations is to use conductivity surfaces or current sheets. Here it is assumed that the material sheet has zero thickness, thus is truly two-dimensional, and the constitutive equation \[\vvJ^{}_{\operatorname{surf}}(\omega) = \sigma_{\operatorname{surf}}(\omega) \vvE^{}(\omega)\] holds along the plane of the material in the frequency domain, where \(\sigma_{\operatorname{surf}}\) describes the surface conductivity of the 2D material. We refer to Chapter 1 in Dep.2016.GOES? for details about the modelling.

As a first step toward the full model, we consider linear and isotropic time-dependent Maxwell equations on a polyhedral domain \(Q\subset \bbR^3\) with plane faces. Furthermore, we assume that \(Q\) is composed of two plane-faced sub-polyhedra \(Q_-\) and \(Q_+\) which share a common face \({F_{\operatorname{int}}}= \overline{Q_-} \cap \overline{Q_+}\). We choose the coordinate frame such that \({F_{\operatorname{int}}}\) is a subset of the \(x_2x_3\)-plane and the normal vector \(\boldsymbol{n}_{\operatorname{int}}\) on \({F_{\operatorname{int}}}\) is aligned with the \(x_1\)-axis. See 1 for different examples we have in mind. We assume that the surface current \(\vvJ^{}_{\operatorname{surf}}\) supported on \({F_{\operatorname{int}}}\) is a given function. The governing equations read \[\tag{1} \begin{align} \partial_t\vvH^{}_\pm&= -\mu^{-1}_\pm\mathop{\mathrm{curl}}\vvE^{}_\pm, & \relax\Bigl(\mu^{}_\pm\vvH^{}_\pm\Bigr) & = 0, \tag{2}\\ \partial_t\vvE^{}_\pm &= \varepsilon^{-1}_\pm\mathop{\mathrm{curl}}\vvH^{}_\pm - \varepsilon^{-1}_\pm \vvJ^{}_\pm, & \relax\Bigl(\varepsilon^{}_\pm\vvE^{}_\pm\Bigr) &= \rho^{}_\pm, \tag{3} \end{align}\] on \(Q_\pm\), for \(t\geq 0\). We denote by \(f_\pm = \restrict[\big]{f}{Q_\pm}\) the restriction of a function \(f\in\lebfuns{2}{Q}\). Here, for \(x\in Q_\pm\), \(\vvE^{}(t,x),\vvH^{}(t,x) \in \bbR^3\) denote the electric and magnetic field, \(\vvJ^{}(t,x)\in\bbR^3\) the volume current density and \(\rho^{}(t,x)\in\bbR\) the charge density, respectively. We assume that the material parameters \(\mu^{}_\pm,\varepsilon^{}_\pm\geq \delta > 0\) are constant on \(Q_\pm\). The equations are equipped with perfectly conducting boundary conditions \[\label{eq:PECBoundary} \mu\vvH^{}\cdot \boldsymbol{\nu}= 0, \vvE^{}\times \boldsymbol{\nu}= 0,\tag{4}\] on \(\partial Q\), for \(t\geq 0\) with outer unit normal vector \(\boldsymbol{\nu}\), see, e.g., Dau.L.1990.MANM?. At the interface \({F_{\operatorname{int}}}\), the conditions \[{\tag{5}} \begin{align} [\![\mu\vvH^{}\cdot\boldsymbol{n}_{\operatorname{int}}]\!]_{F_{\operatorname{int}}}&= 0, & [\![\varepsilon\vvE^{}\cdot \boldsymbol{n}_{\operatorname{int}}]\!]_{F_{\operatorname{int}}}&= \rho^{}_{\operatorname{surf}}, \tag{6}\\ [\![\vvH^{}\times\boldsymbol{n}_{\operatorname{int}}]\!]_{F_{\operatorname{int}}}&= \vvJ^{}_{\operatorname{surf}}, & [\![\vvE^{}\times\boldsymbol{n}_{\operatorname{int}}]\!]_{F_{\operatorname{int}}}&= 0,\tag{7} \end{align}\] hold for \(t\geq 0\), where \(\boldsymbol{n}_{\operatorname{int}}\) denotes the inner unit normal vector on \({F_{\operatorname{int}}}\) pointing from \(Q_-\) to \(Q_+\) and \([\![f]\!]_{F_{\operatorname{int}}}= \restrict[\big]{f_+}{{F_{\operatorname{int}}}} - \restrict[\big]{f_-}{{F_{\operatorname{int}}}}\) denotes the jump on \({F_{\operatorname{int}}}\) whenever the functions \(f_\pm\) admit well-defined traces on the interface.

By 7 , we obtain the condition \[\vvJ^{}_{\operatorname{surf}} \cdot \boldsymbol{n}_{\operatorname{int}}= [\![\vvH\times\boldsymbol{n}_{\operatorname{int}}]\!] \cdot \boldsymbol{n}_{\operatorname{int}}= 0,\] since the tangential trace of the magnetic field belongs to the tangential plane of \(\boldsymbol{n}_{\operatorname{int}}\). With the specific coordinate frame chosen above, this reads as \[\vvJ^{}_{\operatorname{surf}}(t,x) = \bigl(0,J^{}_{\operatorname{surf},2}(t,x),J^{}_{\operatorname{surf},3}(t,x)\bigr).\]

Finally, by \(\rho^{}_{\operatorname{surf}}(t,x) \in \bbR\) we denote the surface charge density, see, e.g., for details Dau.L.1990.MANM?.

Discretization↩︎

The discontinuous Galerkin (dG) time-domain method is a well-establish method for Maxwell equations, see, e.g., Hes.W.2002?. We briefly recall the construction of the dG space discretization and refer to ¿sec:subsec:DiscreteSetting? ¿sec:subsec:SpatialDiscretization? for details.

Assume that \(\calT_h\) is the union of suitable meshes for \(Q_\pm\) with elements \(K\) and matching faces at \({F_{\operatorname{int}}}\). The set of all element faces \(F\) is denoted by \(\calF_h\). The broken polynomial space of degree at most \(k\geq 1\) is defined as \[\tag{8} \begin{equation}\tag{9} \PolyAtMost{3}{k}[\calT_h] = \setc[\big]{v_h\in\lebfuns{2}{Q}}{\restrict{v_h}{K} \in \PolyAtMost{3}{k}[K] \text{ for all } K\in\calT_h}, \end{equation} which are polynomials on every element \(K\) and, in general, discontinuous{} across element faces \(F\). The vector valued ansatz space for the magnetic and electric field is given by the broken finite element space of degree \(k\) defined as \begin{equation}\tag{10} V_h= \PolyAtMost{3}{k}[\calT_h]^3. \end{equation}\]

The dG method is a non-conforming method in the sense that a function \(\vvU_h \in V_h\) does not admit a \(\mathop{\mathrm{curl}}\) on the whole domain \(Q\), i.e., \(V_h\not\subset \curlfuns[]{Q}\). Therefore, we need to introduce a discretized \(\mathop{\mathrm{curl}}\) operator acting on \(V_h\). One way to do so is given by means of the central flux discretization, see, e.g., Hes.W.2002?. The notation is based on Hoc.S.2016?. We define the discrete operator \(\mathop{\mathrm{curl}}_h {:}V_h\to V_h\) such that for all \(\boldsymbol{\phi}_h\in V_h\) it holds \[\label{eq:DiscreteCurlIdea} \begin{align} \int_{Q}\mathop{\mathrm{curl}}_h\vvU_h\cdot\boldsymbol{\phi}_h\ \mathrm{d}x =& \sum_{K\in\calT_h} \int_{K}\mathop{\mathrm{curl}}\vvU_h \cdot\boldsymbol{\phi}_h\ \mathrm{d}x\\ &- \sum_{F\in\calF_h} \int_{F}[\![\vvU_h\times \vvn_{F}]\!]_F\cdot \{\{\boldsymbol{\phi}_h\}\}_F\ \mathrm{d}s. \end{align}\tag{11}\] Here, \([\![\vvU_h\times\vvn_F]\!]_F\) denotes the tangential jump of \(\vvU_h\) and \(\{\{\boldsymbol{\phi}_h\}\}_F\) a weighted average of \(\boldsymbol{\phi}_h\) on \(F\) defined below 28 . The first sum on the right hand side of 11 acts locally on single elements, thus decoupling the action of the curl while the second sum couples neighboring elements through tangential jumps. The coupling terms are referred to as numerical fluxes and they penalize non-zero tangential jumps across faces. We recall that functions \(\vvH^{}\in\curlfuns{Q}\) have zero tangential jumps across faces, i.e., \[[\![\vvH^{}\times \vvn_{F}]\!]_F= 0.\] By \(\mathop{\mathrm{curl}}_{h,0} {:}V_h\to V_h\) we denote a discrete operator that additionally enforces homogeneous tangential boundary conditions of the electric field 4 .

In order to incorporate the surface current \(\vvJ^{}_{\operatorname{surf}}\), we follow the idea of 11 , but instead of penalizing zero tangential jumps, we apply the inhomogeneous interface condition 7 for all faces \(F\in\calF_h^\circ\) with \(F\subset{F_{\operatorname{int}}}\), i.e., \[[\![\vvH^{}\times \vvn_F]\!]_F= \restrict[\big]{\vvJ^{}_{\operatorname{surf}}}{F}.\] Equation 11 motivates to define an extension \(\vvJ^{}_{\operatorname{surf},h} \in V_h\) via \[\label{eq:DiscreteLiftIdea} \int_Q\vvJ^{}_{\operatorname{surf},h}\cdot \boldsymbol{\phi}_h\ \mathrm{d}x = \sum_{F\in\calF_h, F\subset{F_{\operatorname{int}}}} \int_{F} \vvJ^{}_{\operatorname{surf}}\cdot \{\{\boldsymbol{\phi}_h\}\}\ \mathrm{d}s\tag{12}\] for all \(\boldsymbol{\phi}_h \in V_h\). Note that \(\vvJ^{}_{\operatorname{surf}}\) is defined only on the interface \({F_{\operatorname{int}}}\) whereas \(\vvJ^{}_{\operatorname{surf},h}\) is defined on the whole domain \(Q\), with support only on elements adjacent to \({F_{\operatorname{int}}}\). We end up with the following spatially discrete system of differential equations \[\label{eq:FirstHandwavingDiscretization} \begin{align} \partial_t\vvH^{}_h(t) &= -\mu^{-1}\mathop{\mathrm{curl}}_{h,0}\vvE^{}_h(t),\\ \partial_t\vvE^{}_h(t) &= \varepsilon^{-1}\mathop{\mathrm{curl}}_h\vvH^{}_h(t) - \vvJ^{}_{h}(t) - \vvJ^{}_{\operatorname{surf},h}(t), \end{align}\tag{13}\] for \(t\geq 0\). Here, the inhomogeneous interface conditions 7 are incorporated by \(\vvJ^{}_{\operatorname{surf},h}\) that acts like an artificial current on the evolution of the electric field. We refer to ¿sec:subsec:DiscreteSetting? ¿sec:subsec:SpatialDiscretization? for the precise definitions.

We integrate the spatially discrete system in time by the second order leapfrog method. Let \(\tau > 0\) be the time step size and \(t_n = n\tau\) for \(n\in\bbN\). The fully discrete scheme then reads \[\begin{align} \vvH^{n+1/2}_h - \vvH^{n}_h &= -\frac{\tau}{2}\mu^{-1}\mathop{\mathrm{curl}}_{h,0} \vvE^{n}_h,\\ \vvE^{n+1}_h - \vvE^{n}_h &= \tau\varepsilon^{-1}\mathop{\mathrm{curl}}_h \vvH^{n+1/2}_h - \frac{\tau}{2}(\vvJ^{n}_{h} + \vvJ^{n+1}_{h}) - \frac{\tau}{2}(\vvJ^{n}_{\operatorname{surf},h} + \vvJ^{n+1}_{\operatorname{surf},h}),\\ \vvH^{n+1}_h - \vvH^{n+1/2}_h &= -\frac{\tau}{2} \mu^{-1}\mathop{\mathrm{curl}}_{h,0} \vvE^{n+1}_h, \end{align}\] for \(n\geq 0\), starting from appropriate initial values \(\bigl(\vvH^{0}_h, \vvE^{0}_h\bigr)\in V_h^2\). Other time integration schemes can be applied to 13 as well.

It is well known that the explicit leapfrog scheme exhibits a step size restriction, which is also known as the Courant-Friedrichs-Lewy (CFL) condition. The scheme is only stable for time step sizes \(\tau < \tau_{\operatorname{CFL}}\), with \(\tau_{\operatorname{CFL}}\sim h_\text{min}\), where \(h_\text{min}\) denotes the diameter of the smallest element \(K\in\calT_h\).

Contributions of the paper↩︎

The challenges associated with interface problems have been thoroughly investigated from both analytical and numerical perspectives, albeit within a geometric framework different from the one specified earlier. In that context, it is assumed that a positive distance exists between the interface and the domain boundary. For example, well-posedness and regularity of quasilinear Maxwell equations for such a geometric setting is found in Sch.S.2022?. From a numerical point of view, finite element methods have been explored in Che.Z.1998? concerning elliptic and parabolic problems, and in Dek.2017?, Dek.S.2012? for hyperbolic equations. However, these results are not applicable to the problem described in 1 , 4 and 5 .

In a recent study by Dörich and Zerulla Dor.Z.2023?, a different technique is employed to establish well-posedness and regularity for the model problem 1 , 4 and 5 . on the cuboid domain depicted in [fig:NormalDomain]. From a numerical perspective, the discontinuous Galerkin time-domain method was successfully applied to an interface problem concerned with Graphene sheets in the above mentioned setting, see, e.g., Wer.W.M.Et.2015?, Wer.K.B.Et.2016?. There, the focus is on the physical modelling of such sheets. Their excellent numerical results motivate a thorough mathematical error analysis.

In this paper, we provide a mathematical framework that is suitable for both, analysis and numerics of the problem at hand. We prove well-posedness and stability for the governing equations building up on the techniques in Dor.Z.2023?. Transferring the ideas from analysis, a rigorous spatial and full discretization error analysis is provided for the numerical scheme. Under suitable regularity conditions on the exact solution, we prove that the error of the scheme is of second order in time and of \(k\)th order in space with respect to the \(L^2\)-norm, i.e., \[\norm{\bigl(\vvH^{},\vvE^{}\bigr)(t_n) - (\vvH^{n}_h, \vvE^{n}_h)}_{\lebfuns{2}{Q}^3 \times\lebfuns{2}{Q}^3} \leq C(\tau^2 + h^{k}), \quad 0 \leq t_n \leq T.\] Furthermore, the results are underpinned by several numerical examples showing the sharpness of the estimates with respect to spatial regularity.

Note that the results are consistent with the case where the surface current vanishes, i.e., \(\vvJ^{}_{\operatorname{surf}} = 0\). However, new techniques are required for \(\vvJ^{}_{\operatorname{surf}} \neq 0\). We address the problems in brief.

One of the challenges is that by the interface condition 7 , the state-space \(\curlfuns{Q}\times \hcurlfuns{Q}\), typically used for the evolution of linear Maxwell equations, is no longer suitable for the problem described by 1 , 4 and 5 . Circumventing this, we enlarge the state-space with functions \(\vvV\in\lebfuns{2}{Q}^3\) that only possess a weak \(\mathop{\mathrm{curl}}\) on each sub-polyhedra, i.e., \(\vvV_\pm \in \curlfuns{Q_\pm}\). This causes several problems both from an analytical and numerical perspective. Analytically, \(C^0\)-semigroup techniques are no longer applicable and numerically, we must handle a non-consistent discretization. Motivated by the treatment of inhomogeneous Dirichlet boundary conditions, we modify the problem that we can treat it in a standard way. Nonetheless, since the interface \({F_{\operatorname{int}}}\) intersects with the boundary, special care is necessary to treat the perfectly conducting boundary conditions 4 correctly in the modified problem. For this, we extend the ideas from Dor.Z.2023? to the more general geometric setting described above. They turn out to be essential for both analysis and numerics.

Structure of the paper↩︎

In 2, we first introduce a suitable analytical framework for the problem described by 1 , 4 and 5 . We proceed by presenting the main result of this section concerning an existence and stability result for the analytical problem. With the strategy of proof outlined, we introduce an important extension result that is later used frequently. The section is closed with the proof of the main result.

3 is concerned with space discretization. We first provide a standard description of the discontinuous Galerkin method and point out in detail how inhomogeneous interface problems are treated. After that, the main result of this section provides an error bound on the semi-discretization. The section carries on with the discussion of an important extension of the semi-discrete scheme utilizing a nodal interpolation on the interface. The section again closes with the remaining proofs.

3 is concerned with space discretization. We provide a detailed description of the discontinuous Galerkin method and point out in detail how inhomogeneous interface problems are treated. We present the main result of this section consisting of error bounds for the semi-discretization as well as an important extension of the semi-discrete scheme utilizing a nodal interpolation on the interface. The section again closes with the remaining proofs.

The main result of 4 is concerned with an error bound on the full discretization. We first prove stability of the scheme and provide afterwards the proof of the main result.

In 5, we provide several numerical experiments in a two-dimensional setting that confirm our theoretical findings.

2 Well-posedness↩︎

General setting and notation↩︎

The volume charge density \(\rho^{}\) is determined by the volume current \(\vvJ^{}\) through \[\label{eq:VolumeChargeCurrentRelation} \rho^{}_\pm(t) = \rho^{}_\pm(0) + \int_0^t \relax\vvJ^{}_\pm(s)\;\mathrm{d}s\tag{14}\] on \(Q_\pm\), for \(t\geq 0\). Equation 14 is called the continuity relation for electricity. It is well known that the divergence conditions 2 3 and the magnetic boundary condition in 4 hold, if they are valid for \(t=0\) and 14 holds. We refer to Dau.L.1990.MANM? for details.

A similar relation exists for the surface charge density \(\rho^{}_{\operatorname{surf}}\). It is determined by the volume current \(\vvJ^{}\) and the surface current \(\vvJ^{}_{\operatorname{surf}}\) through \[\label{eq:SurfaceChargeCurrentRelation} \rho^{}_{\operatorname{surf}}(t) = \rho^{}_{\operatorname{surf}}(0) + \int_0^t \relax_{F_{\operatorname{int}}}\vvJ^{}_{\operatorname{surf}}(s) - [\![\vvJ^{}(s)\cdot\boldsymbol{n}_{\operatorname{int}}]\!]_{F_{\operatorname{int}}}\;\mathrm{d}s\tag{15}\] on \({F_{\operatorname{int}}}\), for \(t\geq 0\), where \(\relax_{F_{\operatorname{int}}}\) denotes the two-dimensional divergence on \({F_{\operatorname{int}}}\). It is shown in Sch.S.2022? that equations 6 are valid for \(t\geq 0\) if they are valid for \(t=0\) and 15 holds. Thus, it remains to solve the \(\mathop{\mathrm{curl}}\)-equations in 1 subject to the boundary conditions 4 and the tangential interface conditions 7 .

The speed of light is denoted by \(c_\pm=(\mu^{}_\pm\varepsilon^{}_\pm)^{-1/2}\) and we use the notation \[\eta_\infty = \max\monosetc{\eta_-,\eta_+}, \quad \eta \in \monosetc{\varepsilon^{},\mu^{},c}\] for piecewise defined constants. We employ the weighted inner \(L^2\)-products \[\label{eq:weightedL2} \sprodmu{\cdot}{\cdot} = \sprod{\mu\cdot}{\cdot}_{\lebfuns{2}{Q}},\sprodeps{\cdot}{\cdot} = \sprod{\varepsilon\cdot}{\cdot}_{\lebfuns{2}{Q}},\sprodmueps{\cdot}{\cdot} = \sprodmu{\cdot}{\cdot} + \sprodeps{\cdot}{\cdot}\tag{16}\] and their induced norms \(\normmu{\cdot}\), \(\normeps{\cdot}\) and \(\normmueps{\cdot}\). Note that they are equivalent to the standard \(\lebfuns{2}{Q}\)-norm.

The spaces \(\sobfuns{s}{Q}\) for \(s>0\) denote fractional Sobolev spaces. We write \(\norm{\cdot}_{\sobfuns{s}{Q}}\) and \(\singlenorm{\cdot}_{\sobfuns{s}{Q}}\) for their associated norms and semi-norms and refer to Edm.E.2023.FSSI? for details. On the interface \({F_{\operatorname{int}}}\), we denote by \(H^{1/2}({F_{\operatorname{int}}})\) the fractional Sobolev space, and by \(H^{-1/2}({F_{\operatorname{int}}})\) its dual with pivot \(L^2({F_{\operatorname{int}}})\).

We introduce the space of functions that exhibit a weak variation \(\mathop{\mathrm{curl}}\) \[\curlfuns[]{Q} = \setc[\big]{\vvV\in\lebfuns{2}{Q}^3}{\mathop{\mathrm{curl}}\vvV \in\lebfuns{2}{Q}^3},\] as well as the subspace that contains functions with a vanishing tangential trace \[\hcurlfuns{Q} = \setc[\big]{\vvV \in \curlfuns{Q}}{\restrict{\vvV\times \boldsymbol{\nu}}{\Gamma} = 0}.\] Since we deal with solutions that lag regularity across the interface \({F_{\operatorname{int}}}\), we employ the notion of piecewise spaces. For \(s > 0\) we denote the piecewise Sobolev spaces \[PH^{s}(Q) = \setc[\big]{v\in\lebfuns{2}{Q}}{v_\pm\in\sobfuns{s}{Q_\pm}},\] and we define analogously \[PH(\mathop{\mathrm{curl}}, Q)= \setc[\big]{\vvV\in\lebfuns{2}{Q}^3}{\mathop{\mathrm{curl}}\vvV_\pm \in \lebfuns{2}{Q_\pm}^3}.\] Note that we use the same symbol for both \(\mathop{\mathrm{curl}}\)-operators. The following connection between \(\curlfuns{Q}\) and \(PH(\mathop{\mathrm{curl}}, Q)\) is a simple corollary of Green’s formula; see, for example, Gir.R.1986.FEMNa?.

\[\tag{17} We define the Maxwell operators acting on the magnetic and electric field respectively as \begin{alignat}{3}\tag{18} \operatorname{\widehat{\boldsymbol{\calC}}_H}&{:}D(\operatorname{\widehat{\boldsymbol{\calC}}_H})= PH(\mathop{\mathrm{curl}}, Q)\to \lebfuns{2}{Q}^3,\; &&\vvH &&\mapsto \varepsilon^{-1}\mathop{\mathrm{curl}}\vvH,\\ \operatorname{\boldsymbol{\calC}_E}&{:}D(\operatorname{\boldsymbol{\calC}_E})= \hcurlfuns{Q} \to \lebfuns{2}{Q}^3,\; &&\vvE &&\mapsto \mu^{-1}\mathop{\mathrm{curl}}\vvE. \end{alignat} The operator acting on the combined field \(\vvu = \bigl(\vvH,\vvE\bigr)\) is defined as \begin{equation}\tag{19} \operatorname{\widehat{\boldsymbol{\calC}}}{:}D(\operatorname{\widehat{\boldsymbol{\calC}}})= D(\operatorname{\widehat{\boldsymbol{\calC}}_H})\times D(\operatorname{\boldsymbol{\calC}_E})\to \lebfuns{2}{Q}^6,\; \operatorname{\widehat{\boldsymbol{\calC}}}= \begin{pmatrix} 0 && -\operatorname{\boldsymbol{\calC}_E}\\ \operatorname{\widehat{\boldsymbol{\calC}}_H}&& 0 \end{pmatrix}. \end{equation} We emphasize that \(\curlfuns{Q} \subset PH(\mathop{\mathrm{curl}}, Q)\) and define the restricted operators \begin{align}\tag{20} \operatorname{\boldsymbol{\calC}_H}{:}D(\operatorname{\boldsymbol{\calC}_H})&= \curlfuns{Q} \to \lebfuns{2}{Q}^3, & \operatorname{\boldsymbol{\calC}_H}&= \restrict{\operatorname{\widehat{\boldsymbol{\calC}}_H}}{D(\operatorname{\boldsymbol{\calC}_H})},\\ \operatorname{\boldsymbol{\calC}}{:}D(\operatorname{\boldsymbol{\calC}})&= \curlfuns{Q} \times \hcurlfuns{Q} \to \lebfuns{2}{Q}^6, & \operatorname{\boldsymbol{\calC}}&= \restrict{\operatorname{\widehat{\boldsymbol{\calC}}}}{D(\operatorname{\boldsymbol{\calC}})}. \end{align} Note that operators with a hat are always associated with piecewise{} domains. We stick to this notation throughout the paper.\]

The Maxwell equations now read: seek \(\bigl(\vvH(t),\vvE(t)\bigr)\in D(\operatorname{\widehat{\boldsymbol{\calC}}_H})\times D(\operatorname{\boldsymbol{\calC}_E})\) such that \[\tag{21} \begin{align} \partial_t\vvH &= -\operatorname{\boldsymbol{\calC}_E}\vvE &&\text{in }[0,T]\times Q,\\ \partial_t\vvE &= \operatorname{\widehat{\boldsymbol{\calC}}_H}\vvH - \varepsilon^{-1} \vvJ^{}&&\text{in }[0,T]\times Q,\\ \vvH(0) &= \vvH^0,\quad \vvE(0) = \vvE^0 &&\text{in }Q,\\ [\![\vvH\times \boldsymbol{n}_{\operatorname{int}}]\!]_{F_{\operatorname{int}}}&= \vvJ^{}_{\operatorname{surf}}&&\text{on }[0,T]\times {F_{\operatorname{int}}},\tag{22} \end{align}\]

Surface current↩︎

As mentioned in the introduction, the normal component of the surface current vanishes, i.e., \(\vvJ^{}_{\operatorname{surf}} \cdot \boldsymbol{n}_{\operatorname{int}}= 0\). Therefore, we can always identify \(\vvJ^{}_{\operatorname{surf}}\) with a two-dimensional vector-valued function in the tangential plane described by \(\boldsymbol{n}_{\operatorname{int}}\). We will therefore use the same notation \(\vvJ^{}_{\operatorname{surf}}\) interchangeably as a three-dimensional vector-valued function perpendicular to the normal vector \(\boldsymbol{n}_{\operatorname{int}}\) and as a two-dimensional function on \({F_{\operatorname{int}}}\) identified as a subset of \(\bbR^2\).

The surface current needs to satisfy certain compatibility conditions for the problem to be well-posed. This boils down to whether \(\vvJ^{}_{\operatorname{surf}}\) is extensible by zero at the boundary of \({F_{\operatorname{int}}}\) in the sense of \(H^{1/2}\). This leads to a delicate notion of spaces, typically denoted with \(H^{1/2}_{00}\) in the literature, caused by the difficulty of assigning a meaning to traces of functions with Sobolev-regularity \(1/2\).

We choose here a definition that is intrinsic to \({F_{\operatorname{int}}}\) interpreted as a bounded polygon in \(\bbR^2\) and, therefore, independent of both polyhedra \(Q_\pm\). We follow the construction and notation in Buf.C.2001?, Buf.C.2001b?, and we stick to the notation of this paper.

Doing so, we denote the edges of \({F_{\operatorname{int}}}\) as \(\Gamma_i \subset \bbR^2\) for \(i\in\monosetc{1,\ldots,N_{{F_{\operatorname{int}}}}}\). The outer unit normal of \(\Gamma_i\) is denoted with \(\nu_i\in\bbR^2\) and the tangential vector of \(\Gamma_i\) as \(\tau_i\in\bbR^2\). See 2 for an example.

We further define the distance functions \[\delta_{\Gamma_i}(x) = \operatorname{dist}(x,\Gamma_i), \text{ for } x\in\bbR^2.\]

The space of functions on \({F_{\operatorname{int}}}\) with Sobolev-regularity \(1/2\) which are extensible by zero to the whole of \(\bbR^2\) is defined as

\[\begin{align} \hsobfunszeroo[]{1/2}{{F_{\operatorname{int}}}} &= \setc[\big]{f\in\sobfuns[]{1/2}{{F_{\operatorname{int}}}}}{\delta_{\Gamma_i}^{-1/2}f \in \lebfuns[]{2}{{F_{\operatorname{int}}}}, i=1,\ldots,N_{{F_{\operatorname{int}}}}},\\ \norm{f}_{\hsobfunszeroo[]{1/2}{{F_{\operatorname{int}}}}}^2 &= \norm{f}_{\sobfuns[]{1/2}{{F_{\operatorname{int}}}}}^2 + \sum_{i=1}^{N_{{F_{\operatorname{int}}}}} \norm{\delta_{\Gamma_i}^{-1/2}f}_{\lebfuns[]{2}{{F_{\operatorname{int}}}}}^2 \end{align}\] The additional constraints thus penalize functions that do not decay rapidly enough to zero in close neighborhood to the boundary of \({F_{\operatorname{int}}}\). The dual space of \(\hsobfunszeroo[]{1/2}{{F_{\operatorname{int}}}}\) with respect to the pivot space \(\lebfuns[]{2}{{F_{\operatorname{int}}}}\) is denoted with \(\hsobfunszeroo[]{-1/2}{{F_{\operatorname{int}}}}\) and equipped with the introduced dual-norm.

Since we deal with vector-valued functions, we further introduce the space of vector fields with Sobolev-regularity \(1/2\) that are extensible by zero parallel to the edges as \[\begin{align} \boldsymbol{H}^{1/2}_{\parallel,00}({F_{\operatorname{int}}})&= \setc[\big]{\vvF\in\sobfuns[]{1/2}{{F_{\operatorname{int}}}}^2}{\delta_{\Gamma_i}^{-1/2}\vvF\cdot\tau_i \in \lebfuns[]{2}{{F_{\operatorname{int}}}}^2, i=1,\ldots,N_{{F_{\operatorname{int}}}}},\\ \norm{\vvF}_{\boldsymbol{H}^{1/2}_{\parallel,00}({F_{\operatorname{int}}})}^2 &= \norm{\vvF}_{\sobfuns[]{1/2}{{F_{\operatorname{int}}}}^2}^2 + \sum_{i=1}^{N_{{F_{\operatorname{int}}}}} \norm{\delta_{\Gamma_i}^{-1/2}\vvF\cdot \tau_i}_{\lebfuns[]{2}{{F_{\operatorname{int}}}}^2}^2 \end{align}\] The dual with respect to \(\lebfuns[]{2}{{F_{\operatorname{int}}}}^2\) is denoted with \(\boldsymbol{H}^{-1/2}_{\parallel,00}({F_{\operatorname{int}}})\) and also equipped with the introduced dual-norm.

We introduce yet another space of functions with Sobolev-regularity \(3/2\), defined by \[\hsobfunszeroo[]{3/2}{{F_{\operatorname{int}}}} = \setc{f\in \hsobfuns[]{1}{{F_{\operatorname{int}}}}}{ \nabla_{{F_{\operatorname{int}}}} F \in \boldsymbol{H}^{1/2}_{\parallel,00}({F_{\operatorname{int}}}) },\] endowed with the associated graph norm. Note that this notation suggests that there are no compatibility conditions in the direction of the normal derivatives. And indeed, this definition is equivalent with norms to the space \(\hsobfuns[]{1}{{F_{\operatorname{int}}}}\cap H^{3/2}({F_{\operatorname{int}}})\). We refer to Buf.C.2001? for details. The dual of \(\hsobfunszeroo[]{3/2}{{F_{\operatorname{int}}}}\) with respect to \(\lebfuns[]{2}{{F_{\operatorname{int}}}}\) is denoted with \(\sobfuns[]{-3/2}{{F_{\operatorname{int}}}}\). We define \(\relax_{{F_{\operatorname{int}}}} {:}\boldsymbol{H}^{-1/2}_{\parallel,00}({F_{\operatorname{int}}})\to \sobfuns[]{-3/2}{{F_{\operatorname{int}}}}\) via duality by \[\sch{\relax_{{F_{\operatorname{int}}}} \vvF}{f}_{3/2} = -\sch{\vvF}{\nabla_{{F_{\operatorname{int}}}}f}_{\parallel,1/2}, \quad f\in \hsobfunszeroo[]{3/2}{{F_{\operatorname{int}}}},\] with the dual pairing in \(\hsobfunszeroo[]{3/2}{{F_{\operatorname{int}}}}\) and in \(\boldsymbol{H}^{1/2}_{\parallel,00}({F_{\operatorname{int}}})\) on the left- and right-hand side, respectively.

We finally define the space \[\label{def:TraceSpaceParallel} \boldsymbol{H}^{-1/2}_{\parallel,00}(\relax_{F_{\operatorname{int}}},{F_{\operatorname{int}}}) = \setc{\vvF\in \boldsymbol{H}^{-1/2}_{\parallel,00}({F_{\operatorname{int}}})}{\relax_{F_{\operatorname{int}}}\vvF \in \hsobfunszeroo[]{-1/2}{{F_{\operatorname{int}}}}},\tag{23}\] endowed with the graph norm. This space will suffice as a domain for the spatial component of the surface current, since the tangential jump is surjective onto this space. The following result follows from Buf.C.2001b?.

Main result↩︎

Our first main result states existence, uniqueness and stability under appropriate regularity assumptions in weak variational spaces.

Strategy of proof↩︎

It is well known that the operator \(\operatorname{\boldsymbol{\calC}}{:}D(\operatorname{\boldsymbol{\calC}})\to\lebfuns{2}{Q}^6\) is the generator of a unitary \(C^0\)-semigroup, whereas the operator \(\operatorname{\widehat{\boldsymbol{\calC}}}{:}D(\operatorname{\widehat{\boldsymbol{\calC}}})\to \lebfuns{2}{Q}^6\) does not inherit any good properties, see, e.g., Dor.Z.2023?. Thus, we define a lifting \(\vvJ^{}_{\vvH} = \mathcal{R}_{{F_{\operatorname{int}}}}\vvJ^{}_{\operatorname{surf}}\) and introduce a shifted magnetic field \(\vtH = \vvH - \vvJ^{}_{\vvH}\) with a vanishing tangential jump \([\![\vtH\times\boldsymbol{n}_{\operatorname{int}}]\!]_{F_{\operatorname{int}}}= 0\). The shifted system then reads: seek \(\bigl(\vtH(t),\vvE(t)\bigr) \in D(\operatorname{\boldsymbol{\calC}_H})\times D(\operatorname{\boldsymbol{\calC}_H})\) such that \[\label{eq:ShiftedMaxwellSystem} \begin{align} \partial_t\vtH &= -\operatorname{\boldsymbol{\calC}_E}\vvE - \vtJ_{1} &&\text{in }[0,T]\times Q,\\ \partial_t\vvE &= \operatorname{\boldsymbol{\calC}_H}\vtH -\vtJ_{2} &&\text{in }[0,T]\times Q,\\ \vtH(0) &= \vtH^0,\quad \vvE(0) = \vvE^0 &&\text{in }Q, \end{align}\tag{25}\] with \(\vtJ_1 = -\partial_t\vvJ^{}_{\vvH}\), \(\vtJ_2 = -\varepsilon^{-1}\vvJ^{}+ \varepsilon^{-1}\mathop{\mathrm{curl}}\vvJ^{}_{\vvH}\) and \(\vtH^0 = \vvH^0 - \vvJ^{}_{\vvH}(0)\). We then use \(C^0\)-semigroup theory to obtain unique existence and stability for the shifted problem, and, therefore also for 21 .

Proof of main result↩︎

We write 25 in the compact form: seek \(\vtu(t) \in D(\operatorname{\boldsymbol{\calC}})\) such that \[\label{eq:ShiftedMaxwellSystemCombined} \begin{align} \partial_t\vtu(t) &= \operatorname{\boldsymbol{\calC}}\vtu(t) + \vtj(t), \text{ for } t\in [0,T],\\ \vtu(0) &= \vtu^0, \end{align}\tag{26}\] with \(\vtu(t) = \bigl(\vtH(t),\vvE(t)\bigr)\), \(\vtj(t) = \bigl(-\partial_t\vvJ^{}_{\vvH}(t), -\varepsilon^{-1}\vvJ^{}(t) + \varepsilon^{-1}\mathop{\mathrm{curl}}\vvJ^{}_{\vvH}(t)\bigr)\), \(\vvj^{}_{\vvH}(t) = \bigl(\vvJ^{}_{\vvH}(t), 0\bigr)\) and \(\vtu^0 = \vvu^0 - \vvj^{}_{\vvH}(0)\). This problem fits into the framework of Cauchy problems and standard theorems for existence and stability can be applied.

We are now able to prove the first main result.

3 Spatial discretization↩︎

In this chapter, we introduce a concrete space discretization and rigorously derive the discrete curl in 11 and the discrete extension in 12 . We first present our two main results involving rigorous error bounds for the spatially discrete scheme and an extended scheme. The chapter then proceeds with a spatially discrete analogue of the stability bound 2 and is concluded with the proofs of the main results.

Discrete setting↩︎

We denote with \(\calT_h\) matching simplicial meshes of the domain \(Q\), generated by a reference element \(\widehat{K}\). The subscript \(h\) indicates the mesh size, defined as \(h = \max_{K\in\calT_h} h_K\), where \(h_K\) denotes the diameter of a mesh element \(K\). Furthermore, we assume that the mesh sequence is shape regular in the sense of Ern.G.2021.FEAI?. Thus, there exists \(\sigma > 0\) independent of \(h\) such that \(h_K\leq \sigma\rho_K\), where \(\rho_K\) denotes the diameter of the largest inscribing ball of \(K\).

We collect the faces \(F\) of all mesh elements in the set \(\calF_h= \calF_h^\circ\cup \calF_h^{\partial}\), where \(\calF_h^\circ\) denotes the set of all faces in the interior of \(Q\) and \(\calF_h^{\partial}\) the set of all faces on the boundary \(\partial Q\). Refer to Ern.G.2021.FEAI? for a precise definition of mesh faces.

The outer unit normal vector of \(K\) is denoted by \(\vvn_{K}\). Every interior face \(F\in\calF_h^\circ\) intersects two elements \(K_{F,l}\) and \(K_{F,r}\). The order of the elements is arbitrary, but fixed. We choose the unit normal \(\vvn_{F}\) to \(F\) pointing from \(K_{F,l}\) to \(K_{F,r}\). For boundary faces \(F\in \calF_h^{\partial}\), we choose the unit normal \(\vvn_{F}\) to \(F\) as the outer unit normal vector \(\vvn_{K}\) of the associated element \(K\).

Let \(F\) be an interior face and \(v{:}Q\to \bbR\) be a function that admits a well-defined trace on \(F\). The weighted average of \(v\) on the face \(F\) is defined as \[\begin{equation}\label{eq:WeightedAverage} \{\{v\}\}^\omega_F= \frac{ \omega_{K_{F,l}}\restrict[\big]{(\restrict{v}{K_{F,l}})}{F} + \omega_{K_{F,r}} \restrict[\big]{(\restrict{v}{K_{F,r}})}{F}}{\omega_{K_{F,l}} + \omega_{K_{F,r}}}, \end{equation} where \(\omega {:}Q\to (0,\infty)\) denotes a positive weight function that is piecewise{} constant, i.e.{}, \(\restrict{\omega}{K} \equiv \omega_K\) for all \(K\in\calT_h\). Analogously, we define the jump of \(v\) on \(F\) as \begin{equation} [\![v]\!]_F= \restrict[\big]{(\restrict{v}{K_{F,r}})}{F} - \restrict[\big]{(\restrict{v}{K_{F,l}})}{F}. \end{equation}\tag{28}\] For vector fields, both definitions hold component-wise.

The following assumption is necessary to resolve the interface conditions 5 .

Similar to the definition of the broken polynomial space 9 , we introduce for \(s\geq 0\) the broken Sobolev space on \(\calT_h\) defined by \[\begin{equation}\label{eq:BrokenSobolev} \sobfuns{s}{\calT_h} = \setc{v\in\lebfuns{2}{Q}}{\restrict{u}{K}\in\sobfuns{s}{K} \text{ for all } K\in \calT_h} . \end{equation} The piecewise{} semi-norm{} on \(\sobfuns{s}{\calT_h}\) is denoted by \(\singlenorm{\cdot}_{\sobfuns{s}{\calT_h}}\) and we define \begin{equation} \norm{\cdot}_{\sobfuns[]{s}{\calT_h}}^2 = \norm{\cdot}_{\lebfuns[]{2}{Q}}^2 + \singlenorm{\cdot}_{\sobfuns{s}{\calT_h}}^2. \end{equation}\tag{29}\]

In the following, more regularity of the solution is assumed such that it admits classical traces on element faces. Therefore, we define the spaces \[\label{eq:MoreRegularSolutionSpaces} \begin{equation} \widehat{V}_\ast^{\vvH}= D(\operatorname{\widehat{\boldsymbol{\calC}}_H})\cap \sobfuns{1}{\calT_h}^3, \quad V^{\vvE}_\ast= D(\operatorname{\boldsymbol{\calC}_E})\cap \sobfuns{1}{\calT_h}^3, \quad\widehat{V}_\ast= \widehat{V}_\ast^{\vvH}\times V^{\vvE}_\ast \end{equation} and the restricted spaces \begin{equation} V^{\vvH}_\ast= D(\operatorname{\boldsymbol{\calC}_H})\cap \sobfuns{1}{\calT_h}^3, \quad V_\ast= V^{\vvH}_\ast\times V^{\vvE}_\ast. \end{equation}\tag{30}\] Since functions of the approximation space \(V_h\), defined in 9 , do not admit a well-defined \(\mathop{\mathrm{curl}}\), we introduce the following spaces containing both the analytical solution and the approximation \[\label{eq:SpacesContainingBothAnaAndApprox} \begin{equation} \widehat{V}_{\ast,h}^{\vvH}= \widehat{V}_\ast^{\vvH}+ V_h, \quad V_{\ast,h}^{\vvE}= V^{\vvE}_\ast+ V_h, \widehat{V}_{\ast,h}= \widehat{V}_{\ast,h}^{\vvH}\times V_{\ast,h}^{\vvE}, \end{equation} and similarly \begin{equation} V_{\ast,h}^{\vvH}= V^{\vvH}_\ast+ V_h, \quad V_{\ast,h}= V_{\ast,h}^{\vvH}\times V_{\ast,h}^{\vvE}. \end{equation}\tag{31}\]

Spatial discretization↩︎

\[\tag{32} As motivated in the introduction with \ref{eq:DiscreteLiftIdea}, we define the discrete lift operator \begin{equation}\tag{33} {\boldsymbol{\frakL}}_{\operatorname{int}}{:}\lebfuns{2}{{F_{\operatorname{int}}}}^3 \to V_h,\quad \sprodeps{{\boldsymbol{\frakL}}_{\operatorname{int}}\vvV}{\boldsymbol{\phi}^{}_h} = -\sum_{F\in\calF_h^{\operatorname{int}}} \sprod{\restrict{\vvV}{F}}{\{\{\boldsymbol{\phi}^{}_h\}\}^{\overline{\mu c}}_F} \end{equation} for \(\boldsymbol{\phi}^{}_h\in V_h\), and the discrete magnetic Maxwell operator \(\operatorname{\widehat{\boldsymbol{\frakC}}_{H}}{:}\widehat{V}_{\ast,h}^{\vvH}\to V_h\) \begin{equation}\tag{34} \sprodeps{\operatorname{\widehat{\boldsymbol{\frakC}}_{H}}\vvH}{\boldsymbol{\phi}^{}_h} = \begin{multlined}[t] \sum_{K\in\calT_h} \sprod{\vvH}{\mathop{\mathrm{curl}}\boldsymbol{\phi}^{}_h}_K\\ -\sum_{F\in\calF_h^{\partial}} \sprod{\vvH\times\vvn_{F}}{\boldsymbol{\phi}^{}_h}_F -\sum_{F\in\calF_h^\circ} \sprod{\{\{\vvH\}\}^{\mu c}_F}{[\![\boldsymbol{\phi}^{}_h]\!]_F\times \vvn_{F}}_F. \end{multlined} \end{equation} Analogously, we define the electric Maxwell operator \(\operatorname{\boldsymbol{\frakC}_E}{:}V_{\ast,h}^{\vvE}\to V_h\) for \(\vvE\in V_{\ast,h}^{\vvE}\) and \(\boldsymbol{\psi}^{}_h\in V_h\) by \begin{equation}\tag{35} \sprodmu{\operatorname{\boldsymbol{\frakC}_E}\vvE}{\boldsymbol{\psi}^{}_h} = \sum_{K\in\calT_h} \sprod{\vvE}{\mathop{\mathrm{curl}}\boldsymbol{\psi}^{}_h}_K - \sum_{F\in\calF_h^\circ} \sprod{\{\{\vvE\}\}^{\varepsilon c}_F}{[\![\boldsymbol{\psi}^{}_h]\!]_F\times \vvn_{F}}_F. \end{equation} This definition incorporates the perfectly conducting boundary condition for the electric field. The discrete operator acting on the combined field is defined as \begin{equation}\tag{36} \operatorname{\boldsymbol{\frakC}}{:}\widehat{V}_{\ast,h}\to V_h^2,\quad \operatorname{\boldsymbol{\frakC}}= \begin{pmatrix} 0 & - \operatorname{\boldsymbol{\frakC}_E}\\ \operatorname{\widehat{\boldsymbol{\frakC}}_{H}}& 0 \end{pmatrix}. \end{equation} Analogously to \ref{eq:NonPiecewiseOperators}, we define the restrictions \begin{align}\tag{37} \operatorname{\boldsymbol{\frakC}_H}{:}V_{\ast,h}^{\vvH}&\to V_h, & \operatorname{\boldsymbol{\frakC}_H}&= \restrict{\operatorname{\widehat{\boldsymbol{\frakC}}_{H}}}{V_{\ast,h}^{\vvH}},\\ \operatorname{\boldsymbol{\frakC}}{:}V_{\ast,h}&\to V_h^2, & \operatorname{\boldsymbol{\frakC}}&= \restrict{\operatorname{\widehat{\boldsymbol{\frakC}}}}{V_{\ast,h}}. \end{align}\]

The semi-discrete problem now reads: seek \(\bigl(\vvH_h(t),\vvE_h(t)\bigr)\in V_h^2\) such that \[\label{eq:SemiDiscreteMaxwell} \begin{align} \partial_t\vvH_h(t) &= -\operatorname{\boldsymbol{\frakC}_E}\vvE_h(t) &&\text{for } t\in [0,T],\\ \partial_t\vvE_h(t) &= \operatorname{\widehat{\boldsymbol{\frakC}}_{H}}\vvH_h(t) -\vvJ^{}_{h}(t) - \vvJ^{}_{\operatorname{surf},h}(t) &&\text{for } t\in [0,T],\\ \vvH_h(0) &= \vvH_h^0,\quad \vvE_h(0) = \vvE_h^0, \end{align}\tag{38}\] where \(\vvJ^{}_{\operatorname{surf},h}= {\boldsymbol{\frakL}}_{\operatorname{int}}\vvJ^{}_{\operatorname{surf}}\), \(\vvH_h^0 = \Pi_h\vvH^0\), \(\vvE_h^0 = \Pi_h\vvE^0\) and \(\vvJ^{}_{h}= \Pi_h\varepsilon^{-1}\vvJ^{}\). We denote with \(\Pi_h{:}\lebfuns{2}{Q}\to\PolyAtMost{3}{k}[\calT_h]\) the broken \(L^2\)-orthogonal projection defined by \[\label{eq:BrokenOrthProjection} \sprod{v - \Pi_hv}{\phi_h}_{\lebfuns{2}{Q}} = 0\qquad \text{for all } \phi_h \in \PolyAtMost{k}{3}[\calT_h].\tag{39}\] For typical properties of this projection, compare Ern.G.2021.FEAI? or Di.E.2012.MADG?.

The second main result gives an error bound on the spatially discrete solution \(\vvu_h=\bigl(\vvH_h,\vvE_h\bigr)\) of 38 . For a sufficiently regular problem, we obtain convergence in the mesh parameter \(h\). The proof is given below.

Note that this agrees with the results obtained for the special case \(\vvJ^{}_{\operatorname{surf}}= 0\), see, e.g., Hes.W.2002?.

Nodal interpolation↩︎

The calculation of the lift operator 33 involves the evaluation of integrals over mesh faces. In practice, those integrals are approximated by quadrature formulas. This can be quite expensive, since the evaluation may be required at every time step. Moreover, in the future, we want to consider currents that depend on the electric field, i.e., \(\vvJ^{}_{\operatorname{surf}}\approx \vvJ^{}_{\operatorname{surf}}(\vvE^{}_h)\). In such a situation, the evaluation of the lift would cause evaluations of the finite element function at every quadrature point which is quite expensive and should be avoided. In such cases, nodal discontinuous Galerkin methods are attractive since they allow for a fast evaluation of integrals and functions, see, e.g., Di.E.2012.MADG? for a detailed discussion. In the following, we construct a scheme that makes use of nodal interpolation of \(\vvJ^{}_{\operatorname{surf}}\) and provide error bounds.

We specify the construction from ¿sec:subsec:DiscreteSetting? and choose \(\calN_k = \dim{\PolyAtMost{3}{k}}\) nodes \(\Sigma_{\widehat{K}} = \monosetc{\sigma_{\widehat{K},1},\ldots,\sigma_{\widehat{K},\calN_k}}\) in the closure of the reference element \(\widehat{K}\). Then, the Lagrange polynomials, defined by \(\theta_{K,i}(\sigma_{K,j}) = \delta_{ij}\) for \(i,j\in\monosetc{1,\ldots,\calN_k}\), form a basis of \(\PolyAtMost{k}{3}[K]\). Thus, we can define for \(\ell>3/2\) the local interpolation operator \[\calI^h_K{:}\sobfuns{\ell}{K} \to \PolyAtMost{k}{3}[K], \quad \calI^h_Kv = \sum_{j=1}^{\calN_k} \restrict{v}{K}(\sigma_{K,i}) \theta_{K,i}\] and hence, the global interpolation operator by restriction, i.e., \[\calI^h{:}\sobfuns{\ell}{\calT_h} \to \PolyAtMost{k}{3}[\calT_h], \quad \restrict{\calI^hv}{K} = \calI^h_Kv, \quad \text{for } K\in\calT_h.\] Note that the interpolation operator acts component-wise for vector fields.

The surface current \(\vvJ^{}_{\operatorname{surf}}\) is only supported on the interface \({F_{\operatorname{int}}}\) and hence we construct an interpolation operator on the sub-mesh \(\calF_h^{\operatorname{int}}\). Therefore, we need the following two assumptions.

We also need to make sure how the nodes of neighboring elements come in contact with each other.

Given [ass:FaceUnisolvence,ass:FaceMatching], we are able to define for \(\kappa > 1\) the local interpolation operator\[\frakI^h_F{:}\sobfuns{\kappa}{F} \to \PolyAtMost{2}{k}[F], \quad\frakI^h_Fv = \sum_{j = 1}^{\frakN_k} \restrict{v}{F}(\sigma_{F,j}) \theta_{F,j}\] and the global interpolation operator \[\frakI^h{:}\sobfuns{\kappa}{\calF_h^{\operatorname{int}}} \to \PolyAtMost{2}{k}[\calF_h^{\operatorname{int}}],\quad \restrict{\frakI^hv}{F} = \frakI^h_Fv, \quad \text{for } F\in\calF_h^{\operatorname{int}}.\]

The problem now reads: seek \(\bigl(\vcH(t), \vcE(t)\bigr) \in V_h^2\) such that \[\label{eq:InterpolateSemiDiscreteMaxwell} \begin{align} \partial_t\vcH_h(t) &= -\operatorname{\boldsymbol{\frakC}_E}\vcE_h(t) &&\text{for } t\in [0,T],\\ \partial_t\vcE_h(t) &= \operatorname{\widehat{\boldsymbol{\frakC}}_{H}}\vcH_h(t) -\vvJ^{}_{h}(t) - \vcJ^{}_{\operatorname{surf},h}(t) &&\text{for } t\in [0,T],\\ \vcH_h(0) &= \vvH_h^0,\quad \vcE_h(0) = \vvE_h^0, \end{align}\tag{41}\] with \(\vcJ^{}_{\operatorname{surf},h} = {\boldsymbol{\frakL}}_{\operatorname{int}}\frakI^h\vvJ^{}_{\operatorname{surf}}\). Note that the semi-discrete solutions of 38 and 41 only differ in the fact that we use nodal interpolation under the lift operator. Our third main result concerns the error introduced by this additional approximation. The proof is given below.

The following corollary follows immediately from 5 7.

Stability↩︎

We proceed by proving a discrete analogue to the stability bound 24 .

The broken \(L^2\)-projection \(\Pi_h\), defined in 39 , has the following piecewise approximation properties, see, e.g., Ern.G.2021.FEAI?.

The following lemma shows an important relation between the Maxwell operators 17 and their discrete counterparts 32 .

The result is stated in Hoc.S.2016?. Thus, we only prove involving the new domain special to the inhomogeneous interface problem.

A proof of this statement is included in Hoc.S.2016?. We emphasize that all estimates there hold since they are local to every element \(K\in\calT_h\) and, thus, do not depend on the domain \(D(\operatorname{\widehat{\boldsymbol{\calC}}})\).

The following Lemma is essential for the well-posedness of the semi-discrete problem. A proof is provided in Hoc.S.2016?.

We infer from the skew-adjointness that \(\operatorname{\widehat{\boldsymbol{\frakC}}}\) is a generator of a unitary \(C^0\)-semigroup on \(V_h^2\). Therefore, the semi-discrete problem 38 has a unique solution \(\vvu_h(t) = \bigl(\vvH_h(t),\vvE_h(t)\bigr) \in V_h^2\) given by the variation-of-constants formula \[\vvu_h(t) = e^{t\operatorname{\widehat{\boldsymbol{\frakC}}}}\vvu_h^0 + \int_0^t e^{(t-s)\operatorname{\widehat{\boldsymbol{\frakC}}}}\bigl(\vvj^{}_{h}(s) + \vvj^{}_{\operatorname{surf},h}(s)\bigr)\mathrm{d}s,\] with \(\vvu_h^0 = \bigl(\vvH_h^0,\vvE_h^0\bigr)\), \(\vvj^{}_{h}= (0,-\vvJ^{}_{h})\) and \(\vvj^{}_{\operatorname{surf},h}= (0,-\vvJ^{}_{\operatorname{surf},h})\).

The following stability bound holds true for the semi-discrete problem and is a discrete analogue to 24 .

Error analysis↩︎

We proceed by proving the main error bounds of this section.

Interpolation error↩︎

The following local estimates for the nodal interpolation hold, compare for example Ern.G.2021.FEAI?.

In order to obtain approximation properties for the local interpolation operator \(\frakI^h\) on the sub-mesh, we need to ensure that \(\calF_h^{\operatorname{int}}\) does not degenerate, i.e. that the sub-mesh is again shape regular. Recall the following notation. For \(F\in\calF_h^{\operatorname{int}}\), we denote with \(h_F\) the largest diameter of \(F\) and with \(\rho_F\) the diameter of the largest inscribing ball of \(F\). It is clear from the definition that \(h_F\leq h_K\). Furthermore, Gan.L.1996? shows that \(\rho_K\leq \rho_F\), i.e., the diameter of the largest inscribing ball of \(K\) is always less or equal to the diameter of the largest inscribing ball of \(F\). Therefore, \[h_K\leq \sigma \rho_K \implies h_F\leq \sigma \rho_F,\] i.e., the sub-mesh \(\calF_h^{\operatorname{int}}\) inherits the shape regularity from \(\calT_h\). We infer again from Ern.G.2021.FEAI? the following approximation properties.

Similar to 4, we obtain an approximation result under the discrete lift operator.

With 8, we have all ingredients to prove the second main result of this section.

4 Full discretization↩︎

In time, we discretize 38 with the explicit leapfrog scheme with step size \(\tau > 0\) and set \(t_n = n\tau\) for \(n\in\bbN\). The fully discrete scheme reads \[\label{eq:LeapfrogScheme} \begin{align} \vvH^{n+1/2}_h - \vvH^{n}_h &= -\frac{\tau}{2}\operatorname{\boldsymbol{\frakC}_E}\vvE^{n}_h,\\ \vvE^{n+1}_h - \vvE^{n}_h &= \tau\operatorname{\widehat{\boldsymbol{\frakC}}_{H}}\vvH^{n+1/2}_h - \frac{\tau}{2}(\vvJ^{n}_{h} + \vvJ^{n+1}_{h}) - \frac{\tau}{2}(\vvJ^{n}_{\operatorname{surf},h} + \vvJ^{n+1}_{\operatorname{surf},h}),\\ \vvH^{n+1}_h - \vvH^{n+1/2}_h &= -\frac{\tau}{2} \operatorname{\boldsymbol{\frakC}_E}\vvE^{n+1}_h, \end{align}\tag{56}\] for \(n\geq 0\) and \(\vvH^{0}_h = \Pi_h\vvH^{0}\), \(\vvE^{0}_h = \Pi_h\vvE^{0}\). It is well-known that the leapfrog scheme is stable if for some \(\theta\in (0,1)\), the CFL condition \[\label{def:CFLCondition} \tau < \tau_{\operatorname{CFL}}= \frac{2\theta}{\normeps{\operatorname{\widehat{\boldsymbol{\frakC}}_{H}}\operatorname{\boldsymbol{\frakC}_E}}} \lesssim \theta\min_{K\in\calT_h} h_K\tag{57}\] is satisfied. Here, \(\normeps{\cdot}\) denotes the induced operator norm, cf.@eq:eq:weightedL2 . For more details on the constant within the CFL condition, we refer to Hoc.S.2016?.

The main result of this section is the following bound on the full discretization error.

Our analysis is inspired by Hoc.S.2016?, where the locally implicit method for linear Maxwell equations is considered. With the discrete Maxwell operators from 32 and \(\vvu^{n}_h = \bigl(\vvH^{n}_h,\vvE^{n}_h\bigr)\), we write 56 in the following one-step formulation \[\tag{59} \begin{equation}\tag{60} \widehat{\boldsymbol{\frakR}}_{-}\vvu^{n+1}_h = \widehat{\boldsymbol{\frakR}}_{+}\vvu^{n}_h + \frac{\tau}{2}\bigl(\vvj^{n+1}_{h} + \vvj^{n}_{h}\bigr) + \frac{\tau}{2}\bigl(\vvj^{n+1}_{\operatorname{surf},h} + \vvj^{n}_{\operatorname{surf},h}\bigr) \end{equation} for \(n\geq 0\) with operators \begin{equation}\tag{61} \widehat{\boldsymbol{\frakR}}_{\pm}{:}V_h^2 \to V_h^2, \quad \widehat{\boldsymbol{\frakR}}_{\pm}= I \pm \frac{\tau}{2}\operatorname{\widehat{\boldsymbol{\frakC}}}- \frac{\tau^2}{4}\operatorname{\widehat{\boldsymbol{\frakD}}} \end{equation} and perturbation operator \begin{equation}\tag{62} \operatorname{\widehat{\boldsymbol{\frakD}}}{:}V_h^2 \to V_h^2, \quad \operatorname{\widehat{\boldsymbol{\frakD}}}= \begin{pmatrix} 0 & 0\\ 0 & \operatorname{\widehat{\boldsymbol{\frakC}}_{H}}\operatorname{\boldsymbol{\frakC}_E} \end{pmatrix}. \end{equation}\] For \(\operatorname{\widehat{\boldsymbol{\frakD}}}=0\), the scheme 59 is equivalent to the Crank–Nicolson method and thus, one can interpret the leapfrog scheme as perturbation of it. We further use the operator \(\widehat{\boldsymbol{\frakR}}{:}V_h^2\to V_h^2\), defined as \(\widehat{\boldsymbol{\frakR}}= \widehat{\boldsymbol{\frakR}}_{-}^{-1}\widehat{\boldsymbol{\frakR}}_{+}\).

Note that for the special choice of \(\boldsymbol{\mathcal{C}_H}^i=0\) and \(\boldsymbol{\mathcal{C}_E}^i=0\) in Hoc.S.2016?, one obtains the leapfrog scheme on the whole spatial domain. Thus, we can use bounds on the operators in 59 from that work.

Stability↩︎

The following theorem provides stability for the fully discrete scheme and is a discrete analogue of the bound provided in 2 for the exact solution.

Error analysis↩︎

As usual, we split the full discretization error into \[\label{eq:FullDisErrorSplit} \vvu^{n} - \vvu^{n}_h = \vvu^{n} - \Pi_h\vvu^{n} + \Pi_h\vvu^{n}- \vvu^{n}_h = \vve_\Pi^n + \vve_h^n.\tag{69}\] Here, \(\vve_\Pi^n = \vvu^{n} - \Pi_h\vvu^{n}\) is the best approximation error and \(\vve_h^n = \Pi_h\vvu^{n}- \vvu^{n}_h\) is the dG-leapfrog-error at time \(t_n\). Since the best approximation error is covered by projection results, cf. 2, we determine the defect \(\vvd^n\) by inserting the projected exact solution into the scheme 60 , i.e., \[\label{eq:ProjectExactSolutionInMethod} \widehat{\boldsymbol{\frakR}}_{-}\Pi_h\vvu^{n+1} = \widehat{\boldsymbol{\frakR}}_{+}\Pi_h\vvu^{n} + \frac{\tau}{2}\bigl(\vvj^{n+1}_{h} + \vvj^{n}_{h}\bigr) + \frac{\tau}{2}\bigl(\vvj^{n+1}_{\operatorname{surf},h} + \vvj^{n}_{\operatorname{surf},h}\bigr) - \vvd^n.\tag{70}\] This allows us to infer the error recursion for the dG-leapfrog scheme.

With this representation of the defect, we are able to prove the main result of this section.

5 Numerical experiments↩︎

In this section, we present several numerical experiments that underline our theoretical findings. The software was build with the Maxwell toolbox TiMaxdG² which is build upon the finite element library deal.II³ dealII95?. The full software with executables for reproduction purposes can be found under

All experiments are conducted in transverse electric (TE) polarization, i.e., \(H_{1} = H_2 = E_3 = 0\), to reduce the computational effort, with material parameters \(\mu^{}_\pm = \varepsilon^{}_\pm = 1\). Thus, we solve the system \[\begin{align} \partial_t H_{3,\pm} &= -\partial_1 E_{2,\pm} + \partial_2 E_{1,\pm}, &&\text{in}\;Q_\pm,\\ \partial_t E_{1,\pm} &= \partial_2 H_{3,\pm} - J_{1,\pm}, &&\text{in}\;Q_\pm,\\ \partial_t E_{2,\pm} &= -\partial_1 H_{3,\pm} - J_{2,\pm}, &&\text{in}\;Q_\pm,\\ H_{3}(0) &= H_{3}^0,\;E_{1}(0) = E_{1}^0,\; E_{2}(0) = E_{2}^0, &&\text{in}\;Q,\\ [\![H_3]\!]_{F_{\operatorname{int}}}&= J_{1,\operatorname{surf}}, &&\text{on}\;{F_{\operatorname{int}}}. \end{align}\]

Cavity solution↩︎

In this experiment, we use the well-known cavity solution, cf. Hoc.K.2022?, to construct a regular reference solution of the interface problem 21 . On each cuboid \(Q_-\), \(Q_+\) we make the ansatz \[\label{eq:CavityAnsatz} \begin{align} H_{3,\pm}(x_1, x_2, t) &= \frac{1}{\omega^{\pm}} (k^{\pm}_1 A_2 - k_2 A^{\pm}_1)\cos{(k_2 x_2)}\cos{(k^{\pm}_1 (x_1+1))}\sin{(\omega^{\pm} t)},\\ E_{1,\pm}(x_1, x_2, t) &= -A^{\pm}_1 \sin{(k_2 x_2)} \sin{(k^{\pm}_1 (x_1+1))} \cos{(\omega^{\pm} t)},\\ E_{2,\pm}(x_1, x_2, t) &= -A_2 \cos{(k_2 x_2)} \sin{(k^{\pm}_1 (x_1+1))} \cos{(\omega^{\pm} t)}, \end{align} with spatial wave numbers \begin{equation} k^{\pm}_1 = \frac{\pi k^\pm}{2}, \quad k_2 = \pi m, \quad k^\pm, m\in\bbN, \end{equation} temporal wave numbers \begin{equation} \omega^{\pm} = \sqrt{\bigl(k^{\pm}_1\bigr)^2 + k_2^2}, \end{equation} and amplitudes \begin{equation} A^{\pm}_1 = - A_2 \frac{ k_2}{k^{\pm}_1},\quad A_2 \in \bbR. \end{equation}\tag{77}\] We choose the data such that \[\begin{align}\tag{78} J_{\operatorname{surf}}(x_2, t) &= \lim_{x_1\to 0^+} H_{3,+}(x_1,x_2,t) - \lim_{x_1\to 0^-} H_{3,-}(x_1,x_2,t),\\\tag{79} J_{1,\pm} &= J_{2,\pm} = 0. \end{align}\] The remaining constants are chosen such that \(J_{\operatorname{surf}} \neq 0\). In our simulation we used the specific values \[k^- = 2, \qquad k^+ = 4, \qquad m = 1, \qquad A_2 = 1.\]

We chose a mesh sequence of 20 meshes with mesh sizes in the range of \(1\cdot 10^{-1}\) and \(1\cdot 10^{-2}\), a fixed time-step width \(\tau = 1\cdot 10^{-4}\), and polynomial degrees between one and four. For each mesh size, we calculated two different numerical solutions differing in the treatment of the surface current 78 . One series of simulations is done with the lifting defined in 33 and one series is done with the interpolation of the surface current described in 41 . At several time steps, we calculated the \(L^2\)-error against the reference solution 77 . 3 depicts the different mesh sizes on the \(x\)-axis and the maximal \(L^2\)-error obtained on the \(y\)-axis. We observe for \(k\)-th order ansatz polynomials \(k\)-th order spatial convergence until a plateau is reached where the error of the time discretization dominates. This agrees with both, 5 3 Additionally, 3 shows that the interpolation of the surface current leads to the same spatial error. This is expected since the surface current 78 is smooth.

Low regularity surface current↩︎

The aim of this experiment is to show the effect of spatial regularity of surfaces currents on the spatial convergence order. We follow the ideas in Hoc.L.O.2020? and construct for \(\alpha \geq 0\) trigonometric polynomials \[\label{eq:TrigonometricPolynomial} f_\alpha(x) = \sum_{j=-M/2+1}^{M/2} \nu_{\alpha,j} e^{ijx}, \qquad x\in [-\pi,\pi], \qquad M = 2^m, m\in\bbN,\tag{80}\] with coefficients \[\begin{align} \nu_{\alpha,0} &= \nu_{\alpha,M/2} = 0,\\ \nu_{\alpha,j} &= \frac{i \cdot r_j}{(1+j^2)^{\frac{1}{2}(\frac{1}{2}+\alpha)}} &&\text{for } j=1,\ldots,\frac{M}{2}-1,\\ \nu_{\alpha,j} &= -\nu_{\alpha,j+M/2} &&\text{for } j=-\frac{M}{2}+1,\ldots,-1. \end{align}\] The factors \((r_j)_{j=1}^{M/2-1}\) are uniformly sampled numbers from the interval \([-1,1]\). In the limit \(M\to\infty\), the sequence of trigonometric polynomials converges to a function in the Sobolev space \(H^{\alpha}_{\operatorname{per}}(-\pi,\pi)\) with norm \[\norm{g}_{\alpha}^2 = 2\pi \sum_{j\in\bbZ} (1+j^2)^\alpha |\widehat{g}_j|^2, \quad g(x) = \sum_{j\in\bbZ} \widehat{g}_j e^{ijx}.\] The coefficients are chosen in such a way that the trace of \(f_\alpha\) vanishes on the boundary \(\{-\pi,\pi\}\) and, thus, satisfy homogenous Dirichlet boundary conditions. By construction, the norm \(\norm{f_\alpha}_{\eta}\) is bounded uniformly in \(M\) for \(\eta\leq \alpha\) and grows otherwise. 4 demonstrates this behavior.

The surface current is defined as \[J_{\operatorname{surf}}(t,x_2) = \frac{1}{\norm{f_\alpha}_0} f_\alpha(2\pi x_2 - \pi) \sin^2(\pi t)\] for \(x_2\in [0,1]\) and \(t\in[0,T]\).

The fully discretize scheme 56 is used without an additional interpolation on the interface. In our experiment, we chose \(M = 2^m\) for \(m = 22\) and a series of regularity coefficients \(\alpha \in [0,4]\). For every \(\alpha\), we calculated a reference solution on fine mesh with polynomial degree \(k_{\operatorname*{ref}} = 3\), mesh size \(h_{\operatorname*{ref}} = 1\cdot 10^{-3}\), and step size \(\tau_{\operatorname*{ref}} = 5\cdot 10^{-5}\). We then compared the \(L^2\)-error of a sequence of solutions on 8 different meshes with mesh sizes in the range between \(1\cdot 10^{-1}\) and \(5\cdot 10^{-3}\) at the end time \(T=1\) against the reference solution and estimated the order of convergence (EOC). The experiment was performed for first and second order polynomials with the fixed time step size \(\tau = 2.5\cdot 10^{-4}\).

5 shows the dependence of the convergence order on the spatial regularity of the surface current. For \(\alpha < 1/2\) little to no convergence is observed. The order then grows linearly for \(\alpha\in [0.5,1.5]\) until it stagnates.

This agrees with 5 provided that one can improve on the results in Dor.Z.2023? to solutions with piecewise regularity \(PH^{s}(Q)\) for \(s=\min\monosetc{\alpha+1/2,2}\), \(\alpha > 1/2\). In particular, looking at the proofs one would expect that the regularity requirements on \(\vvJ^{}_{\operatorname{surf}}\) can be reduced by one order of Sobolev regularity.

Polynomial solution↩︎

In this example, we investigate the temporal errors by constructing a polynomial solution which does not create spatial errors. We construct a solution that is polynomial in space in order to isolate the error introduced by time discretization. \[\label{eq:pol95sol} The ansatz polynomials in space are given by \begin{equation} q_-(x_1) = 2 + x_1, \quad q_+(x_1) = -1 + x_1, \quad r(x_2) = x_2(1-x_2) \end{equation} and the ansatz function in time by \begin{equation} p(t) = \sin(2\pi t). \end{equation} We define on \(Q_\pm\) the fields \begin{align} H_{3,\pm}(x_1, x_2, t) &= p(t)q_\pm(x_1)r'(x_2),\\ E_{1,\pm}(x_1, x_2, t) &= p'(t)q_\pm(x_1)r(x_2),\\ E_{2,\pm}(x_1, x_2, t) &= 0 \end{align} and define the surface current again as \begin{equation} J_{\operatorname{surf}}(x_2, t) = \lim_{x_1\to 0^+} H_{3,+}(x_1,x_2,t) - \lim_{x_1\to 0^-} H_{3,-}(x_1,x_2,t). \end{equation} Additionally, we define the volume current on \(Q_\pm\) as \begin{align} J_{1,\pm}(x_1,x_2,t) &= p(t)q_\pm(x_1)r''(x_2),\\ J_{2,\pm}(x_1,x_2,t) &= -p(t)q_\pm'(x_1)r'(x_2). \end{align}\tag{81}\]

The fully discretized scheme 56 is used without an additional interpolation on the interface. We chose a mesh sequence with 5 different mesh sizes between \(5\cdot 10^{-1}\) and \(5\cdot 10^{-2}\). For every mesh in the sequence, we compared the \(L^2\)-error between the reference solution and the numerical scheme at several time steps for a total of 40 different time step sizes in the range between \(1\cdot 10^{-1}\) and \(1\cdot 10^{-4}\). Throughout all calculations, we used a polynomial degree of 3 in order to discretize the reference solution exactly in space. 6 shows on the x-axis the time stepsize \(\tau\) and the maximal \(L^2\)-error on y-axis. The method converges with second order in time if the CFL condition 57 is satisfied.

Acknowledgement↩︎

We thank Kurt Busch for the fruitful discussions on plasmonic nanogaps and the insights into the exciting physical phenomena, and Constantin Carle for many helpful discussions on the implementation.

---

1. Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 258734477 – SFB 1173↩︎

2. https://gitlab.kit.edu/kit/ianm/ag-numerik/projects/dg-maxwell/timaxdg↩︎

3. https://www.dealii.org↩︎