Constraint-preserving hybrid finite element methods for Maxwell's equations

TL;DR

This paper introduces a family of hybrid finite element methods for Maxwell's equations that strongly preserve physical constraints like charge conservation and magnetic monopole nonexistence, improving accuracy and reliability.

Contribution

The authors develop a domain decomposition-based hybrid finite element framework that ensures strong preservation of Maxwell's physical constraints, including a hybridized Nédélec method with enhanced constraint enforcement.

Findings

Methods strongly preserve divergence constraints in simulations.

Numerical experiments confirm improved constraint preservation.

Superconvergence observed with hybrid post-processing enhances magnetic field estimates.

Abstract

Maxwell's equations describe the evolution of electromagnetic fields, together with constraints on the divergence of the magnetic and electric flux densities. These constraints correspond to fundamental physical laws: the nonexistence of magnetic monopoles and the conservation of charge, respectively. However, one or both of these constraints may be violated when one applies a finite element method to discretize in space. This is a well-known and longstanding problem in computational electromagnetics. We use domain decomposition to construct a family of primal hybrid finite element methods for Maxwell's equations, where the Lagrange multipliers are shown to correspond to a numerical trace of the magnetic field and a numerical flux of the electric flux density. Expressing the charge-conservation constraint in terms of this numerical flux, we show that both constraints are strongly preserved. As a special case, these methods include a hybridized version of N\'ed\'elec's method, implying that it preserves the constraints more strongly than previously recognized. These constraint-preserving properties are illustrated using numerical experiments in both the time domain and frequency domain. Additionally, we observe a superconvergence phenomenon, where hybrid post-processing yields an improved estimate of the magnetic field.

Tables

Table 1. Table 1. Convergence of the hybridized method for the ω 2 = 2 superscript 𝜔 2 2 \omega^{2}=2 eigenmode of Ω = ( 0 , π ) 2 Ω superscript 0 𝜋 2 \Omega=(0,\pi)^{2} , using a uniform triangle mesh on an N × N 𝑁 𝑁 N\times N grid and degree- ( r − 1 ) 𝑟 1 (r-1) BDM edge elements. The post-processed solution H ^ h subscript ^ 𝐻 ℎ \widehat{H}_{h} exhibits superconvergence relative to H h subscript 𝐻 ℎ H_{h} , while the errors and convergence rates of D ^ h subscript ^ 𝐷 ℎ \widehat{{D}}_{h} are comparable to those of D h subscript 𝐷 ℎ {D}_{h} .

$r$	mesh	$\parallel H_h-H\parallel$		$\parallel {\hat{H}}_h-H\parallel$		$\parallel D_h-D\parallel$		$\parallel {\hat{D}}_h-D\parallel$
	$N$	error	rate	error	rate	error	rate	error	rate
2	2	7.591e-01	—	3.648e-01	—	4.324e-01	—	4.644e-01	—
	4	3.778e-01	1.007	1.070e-01	1.770	1.182e-01	1.872	1.342e-01	1.791
	8	1.862e-01	1.021	2.753e-02	1.958	3.009e-02	1.974	3.512e-02	1.934
	16	9.271e-02	1.006	6.926e-03	1.991	7.558e-03	1.993	8.906e-03	1.979
	32	4.630e-02	1.002	1.734e-03	1.998	1.892e-03	1.998	2.236e-03	1.994
3	2	2.090e-01	—	3.500e-02	—	7.521e-02	—	8.055e-02	—
	4	5.517e-02	1.922	2.750e-03	3.670	9.817e-03	2.938	9.960e-03	3.016
	8	1.400e-02	1.978	1.827e-04	3.912	1.225e-03	3.002	1.220e-03	3.029
	16	3.515e-03	1.994	1.159e-05	3.978	1.526e-04	3.005	1.512e-04	3.013
	32	8.796e-04	1.999	7.270e-07	3.995	1.903e-05	3.003	1.882e-05	3.006
4	2	4.614e-02	—	4.327e-03	—	1.281e-02	—	1.316e-02	—
	4	6.121e-03	2.914	1.250e-04	5.114	7.958e-04	4.008	8.629e-04	3.931
	8	7.769e-04	2.978	3.759e-06	5.055	4.913e-05	4.018	5.500e-05	3.972
	16	9.749e-05	2.994	1.155e-07	5.024	3.048e-06	4.011	3.454e-06	3.993
	32	1.220e-05	2.999	3.582e-09	5.011	1.898e-07	4.006	2.160e-07	3.999
5	2	8.100e-03	—	4.419e-04	—	1.737e-03	—	1.761e-03	—
	4	5.354e-04	3.919	6.307e-06	6.131	5.553e-05	4.968	5.321e-05	5.048
	8	3.394e-05	3.980	9.434e-08	6.063	1.743e-06	4.993	1.642e-06	5.018
	16	2.129e-06	3.995	1.447e-09	6.027	5.449e-08	5.000	5.105e-08	5.007
	32	1.332e-07	3.999	2.404e-11	5.911	1.702e-09	5.000	1.592e-09	5.003

Equations

$$\dot{B}$$

$$\dot{D}+J$$

$$\mathrm{div}B$$

$$\mathrm{div}D$$

$$E:=-(\dot{A}+\mathrm{grad}\phi ),B:=\mathrm{curl}A.$$

$$\dot{A}=-{ϵ}^{-1}D,\dot{D}+J=\mathrm{curl}({\mu }^{-1}\mathrm{curl}A),$$

$$-{\partial }_t(ϵ\dot{A})+J=\mathrm{curl}({\mu }^{-1}\mathrm{curl}A).$$

$$\ddot{A}+\mathrm{curl}\mathrm{curl}A=0.$$

$$H^1(\Omega )$$

$$H(\mathrm{curl};\Omega )$$

$$H(\mathrm{div};\Omega )$$

$${\mathring{H}}^1(\Omega )$$

$$\mathring{H}(\mathrm{curl};\Omega )$$

$$\mathring{H}(\mathrm{div};\Omega )$$

$$L(A,\varphi,\dot{A},\dot{\varphi})\coloneqq\int_{\Omega}\biggl{(}\frac{1}{2}E\cdot{D}-\frac{1}{2}B\cdot H+A\cdot J-\varphi\rho\biggr{)}.$$

$$\displaystyle\int_{\Omega}\bigl{(}A^{\prime}\cdot(\dot{{D}}+J)-\operatorname{curl}A^{\prime}\cdot H\bigr{)}$$

$${\int }_{\Omega }(\mathrm{grad}{\phi }^{\prime }\cdot D+{\phi }^{\prime }\rho )$$

$$L(A,\dot{A})=\int_{\Omega}\biggl{(}\frac{1}{2}E\cdot{D}-\frac{1}{2}B\cdot H+A\cdot J\biggr{)},$$

$$\int_{\Omega}\bigl{(}A_{h}^{\prime}\cdot(\dot{{D}}_{h}+J)-\operatorname{curl}A_{h}^{\prime}\cdot H_{h}\bigr{)}=0,\qquad\forall A_{h}^{\prime}\in V_{h}^{1},$$

$${\dot{B}}_h=-\mathrm{curl}{E}_h,\mathrm{div}{B}_h=0,$$

$${\int }_{\Omega }(\mathrm{grad}{\xi }_h\cdot {\dot{D}}_h+{\xi }_h\dot{\rho })=0.$$

$$H^1({\mathcal{T}}_h)$$

$$H(\mathrm{curl};{\mathcal{T}}_h)$$

$$H(\mathrm{div};{\mathcal{T}}_h)$$

$$\mathring{H}^{1}(\Omega)=\bigl{\{}u\in H^{1}(\mathcal{T}_{h}):\textstyle\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}u\lambda\cdot\mathbf{n}=0,\text{ for all }\lambda\in H(\operatorname{div};\Omega)\bigr{\}}.$$

$$H(\operatorname{div};\Omega)=\bigl{\{}u\in H(\operatorname{div};\mathcal{T}_{h}):\textstyle\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}u\lambda\cdot\mathbf{n}=0,\text{ for all }\lambda\in\mathring{H}^{1}(\Omega)\bigr{\}}.$$

$$\sum _{K\in {\mathcal{T}}_h}{\int }_{\partial K}(u\times \lambda )\cdot \mathbf{n}$$

$$={\int }_{\Omega }(\mathrm{curl}u\cdot \lambda -u\cdot \mathrm{curl}\lambda )$$

$$={\int }_{\partial \Omega }(u\times \lambda )\cdot \mathbf{n}$$

$$=0,$$

$$\displaystyle\biggl{\lvert}\int_{\Omega}u\cdot\operatorname{curl}\lambda\biggr{\rvert}$$

$$\displaystyle=\Biggl{\lvert}\sum_{K\in\mathcal{T}_{h}}\int_{K}\operatorname{curl}u\cdot\lambda\Biggr{\rvert}$$

$$\displaystyle\leq\Biggl{(}\sum_{K\in\mathcal{T}_{h}}\lVert\operatorname{curl}u\rVert^{2}_{L^{2}(K,\mathbb{R}^{3})}\Biggr{)}^{1/2}\lVert\lambda\rVert_{L^{2}(\Omega,\mathbb{R}^{3})},$$

$$L(A,\varphi,\widehat{H},\widehat{{D}},\dot{A},\dot{\varphi},\dot{\widehat{H}},\dot{\widehat{{D}}})=\sum_{K\in\mathcal{T}_{h}}\Biggl{[}\int_{K}\biggl{(}\frac{1}{2}E\cdot{D}-\frac{1}{2}B\cdot H+A\cdot J-\varphi\rho\biggr{)}\\ +\int_{\partial K}(A\times\widehat{H}+\varphi\widehat{{D}})\cdot\mathbf{n}\Biggr{]}.$$

$$\displaystyle\int_{K}\bigl{(}A^{\prime}\cdot(\dot{{D}}+J)-\operatorname{curl}A^{\prime}\cdot H\bigr{)}+\int_{\partial K}(A^{\prime}\times\widehat{H})\cdot\mathbf{n}$$

$${\int }_K(\mathrm{grad}{\phi }^{\prime }\cdot D+{\phi }^{\prime }\rho )-{\int }_{\partial K}{\phi }^{\prime }\hat{D}\cdot \mathbf{n}$$

$$\sum _{K\in {\mathcal{T}}_h}{\int }_{\partial K}(A\times {\hat{H}}^{\prime })\cdot \mathbf{n}$$

$$\sum _{K\in {\mathcal{T}}_h}{\int }_{\partial K}\phi {\hat{D}}^{\prime }\cdot \mathbf{n}$$

$${\int }_{\partial K}(A^{\prime }\times \hat{H})\cdot \mathbf{n}={\int }_K(\mathrm{curl}{A}^{\prime }\cdot H-A^{\prime }\cdot \mathrm{curl}H)={\int }_{\partial K}(A^{\prime }\times H)\cdot \mathbf{n},\forall A^{\prime }\in H(\mathrm{curl};K),$$

Full text

We use domain decomposition to construct a family of primal hybrid finite element methods for Maxwell’s equations, where the Lagrange multipliers are shown to correspond to a numerical trace of the magnetic field and a numerical flux of the electric flux density. Expressing the charge-conservation constraint in terms of this numerical flux, we show that both constraints are strongly preserved. As a special case, these methods include a hybridized version of Nédélec’s method, implying that it preserves the constraints more strongly than previously recognized. These constraint-preserving properties are illustrated using numerical experiments in both the time domain and frequency domain. Additionally, we observe a superconvergence phenomenon, where hybrid post-processing yields an improved estimate of the magnetic field.

1. Introduction

Maxwell’s equations consist of two vector evolution equations, together with two scalar constraint equations, $\operatorname{div}B=0$ and $\operatorname{div}{D}=\rho$, where $B$ is magnetic flux density, $D$ is electric flux density, and $\rho$ is charge density. These constraints are automatically preserved by the evolution, so given initial conditions satisfying the constraints, one can simply evolve forward in time without needing to “enforce” the constraints in any way.

However, if one applies a finite element method in space, then the semidiscretized evolution equations no longer necessarily preserve these constraints, at least not strongly. Nédélec [29] showed that, if one uses curl-conforming edge elements for the electric field $E$ and divergence-conforming face elements for $B$, then the semidiscretized equations preserve $\operatorname{div}B=0$ strongly. On the other hand, $\operatorname{div}{D}=\rho$ holds only in the Galerkin sense (i.e., when both sides are integrated against certain continuous, piecewise-polynomial test functions). Observing this, Christiansen and Winther [15] commented that strong preservation of both divergence constraints “appears to be necessary for many applications in electromagnetics,” and Houston et al. [20] call this “one of the main difficulties in the numerical solution of Maxwell’s equations.” For this reason, alternative approaches have been developed that enforce the constraints strongly—for instance, using Lagrange multipliers [4, 13]—instead of attempting to preserve them automatically but weakly, as Nédélec’s method does. In cases where $\rho=0$, another idea is to use divergence-free elements to construct nonconforming methods [8, 10] or discontinuous Galerkin methods [16, 20, 9].

In this paper, we attack the problem of constraint preservation from a different perspective. We perform domain decomposition of the Lagrangian (i.e., primal) variational principle for Maxwell’s equations, in terms of the vector potential $A$ and scalar potential $\varphi$, using Lagrange multipliers $\widehat{H}$ and $\widehat{{D}}$ to enforce inter-element continuity and boundary conditions. These Lagrange multipliers are shown to correspond to boundary traces of the magnetic field $H$ and electric flux density $D$. After using gauge symmetry to fix $\varphi=0$, we show that the evolution of $(A,\widehat{H})$ automatically preserves the constraints $\operatorname{div}B=0$ and $\operatorname{div}\widehat{{D}}=\rho$. Finally, we semidiscretize this domain-decomposed variational principle, obtaining primal hybrid finite element methods that preserve this formulation of the constraints in a strong sense. As a special case, we give a hybridized formulation of Nédélec’s method, implying that it preserves the constraints in a stronger sense than previously recognized.

The paper is organized as follows:

• •

In Section 2, we review Maxwell’s equations, the Lagrangian variational principle, and semidiscretization using edge elements.

• •

In Section 3, we domain decompose the Lagrangian variational principle, relate solutions to the classical (non-domain-decomposed) formulation of Maxwell’s equations, and study the domain-decomposed version of the constraints and their preservation.

• •

In Section 4, we consider primal hybrid finite element methods for semidiscretizing the domain-decomposed evolution equations, showing that constraints are preserved in a strong sense.

• •

Finally, in Section 5 we conduct numerical experiments demonstrating the behavior of the hybridized Nédélec method. In addition to the constraints being preserved to machine precision, these results illustrate a superconvergence phenomenon for the post-processed magnetic field $\widehat{H}_{h}$, similar to that observed for other hybridized mixed methods (cf. Arnold and Brezzi [2], Brezzi et al. [11]).

2. Maxwell’s equations

2.1. Maxwell’s equations

We begin by reviewing the classical formulation of Maxwell’s equations, first in terms of the electric and magnetic fields and flux densities, and then in terms of the vector and scalar potentials. We postpone the discussion of regularity until the introduction of the weak formulation, in Section 2.2; for the moment, everything may be assumed to be smooth.

2.1.1. Standard formulation

In their most familiar form, Maxwell’s equations consist of the vector evolution equations,

[TABLE]

together with the scalar constraint equations,

[TABLE]

Here, $E$ and $H$ denote the electric field and magnetic field, ${D}=\epsilon E$ and $B=\mu H$ denote the electric flux density and magnetic flux density, $\epsilon$ and $\mu$ are the electric permittivity and magnetic permeability tensors, and $J$ and $\rho$ are current density and charge density, respectively. We use the “dot” notation $\dot{u}\coloneqq\partial_{t}u$ to denote partial differentiation with respect to time.

The evolution equations (1) automatically preserve the constraints (2). Indeed, taking the divergence of (1a) implies $\operatorname{div}\dot{B}=0$, so (2a) is preserved. Similarly, taking the divergence of (1b) implies $\operatorname{div}\dot{{D}}+\operatorname{div}J=0$, so (2b) is preserved if and only if $J$ and $\rho$ satisfy $\dot{\rho}+\operatorname{div}J=0$, which is the law of conservation of charge. We refer to (2b) as the charge-conservation constraint, since it is equivalent to this condition.

2.1.2. Formulation in terms of potentials

Alternatively, Maxwell’s equations may be expressed in terms of a vector field $A$, called the vector potential, and a scalar field $\varphi$, called the scalar potential. Given $A$ and $\varphi$, we define the electric field and magnetic flux density by

[TABLE]

Note that (1a) and (2a) are automatically satisfied, so we may restrict our attention entirely to the single evolution equation (1b), which we have already seen preserves (2b).

However, Maxwell’s equations do not uniquely determine the evolution of $(A,\varphi)$. Observe that if $\xi$ is any time-dependent scalar field, then the transformation $(A,\varphi)\mapsto(A+\operatorname{grad}\xi,\varphi-\dot{\xi})$ leaves $E$, $B$, $D$, $H$ unchanged. Such transformations are called gauge transformations, and the invariance of Maxwell’s equations under gauge transformations is called gauge symmetry. In particular, any solution $(A,\varphi)$ may be transformed into one of the form $(A+\operatorname{grad}\xi,0)$ by taking $\xi$ to be a solution of $\dot{\xi}=\varphi$. Therefore, we may restrict our attention to solutions with $\varphi=0$.

Remark \theremark.

This procedure of restricting to particular solutions, which are related to a general solution by some gauge transformation, is called gauge fixing. The choice $\varphi=0$, called temporal gauge, is the most convenient for our purposes, but there are other choices as well. Note that there is still some remaining gauge symmetry, even after performing temporal gauge fixing: we may transform $A\mapsto A+\operatorname{grad}\xi$ for any $\xi$ constant in time.

After temporal gauge fixing, we can write (1b) as either a first-order system in $A$, $D$,

[TABLE]

or as a second-order equation in $A$ alone,

[TABLE]

In the special case where $\epsilon$ and $\mu$ are simply positive constants with $\epsilon\mu=1$ (as in vacuum, with units chosen so that the speed of light is $1$) and $J=0$, the latter equation just becomes

[TABLE]

Taking the Fourier transform with respect to time (the so-called frequency domain or time-harmonic approach), this latter equation transforms into the eigenvalue problem for the $\operatorname{curl}\operatorname{curl}$ operator.

2.2. Weak formulation

We next discuss the weak formulation of Maxwell’s equations, first using a Lagrangian variational principle in terms of the potentials $A$ and $\varphi$, and then fixing the temporal gauge $\varphi=0$ to arrive at a weak formulation in terms of $A$ alone.

2.2.1. Function spaces and regularity

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded Lipschitz domain, and define the function spaces

[TABLE]

We also define the following subspaces, with boundary conditions imposed:

[TABLE]

Here, $\mathbf{n}\rvert_{\partial\Omega}$ denotes the outer unit normal to $\partial\Omega$, and restrictions to $\partial\Omega$ are interpreted in the trace sense.

Let $A\colon t\mapsto A(t)$ be a $C^{1}$ curve in $\mathring{H}(\operatorname{curl};\Omega)$ and $\varphi\colon t\mapsto\varphi(t)$ be a $C^{0}$ curve in $\mathring{H}^{1}(\Omega)$. It follows that $E$ is a $C^{0}$ curve in $\mathring{H}(\operatorname{curl};\Omega)$, that $B$ is a $C^{1}$ curve in $\mathring{H}(\operatorname{div};\Omega)$, and that (1a) and (2a) hold strongly in $L^{2}$. We also assume that both $\epsilon=\epsilon_{ij}(x,t)$ and $\mu=\mu_{ij}(x,t)$ are $L^{\infty}$, symmetric, and uniformly elliptic. In particular, this implies that ${D}$ and $H$ are both $C^{0}$ curves in $L^{2}(\Omega,\mathbb{R}^{3})$. Henceforth, we restrict our attention to $(A,\varphi)$ such that ${D}$ is in fact a $C^{1}$ curve in $L^{2}(\Omega,\mathbb{R}^{3})$.

Finally, let the current density $J$ be a given $C^{0}$ curve in $H(\operatorname{div};\Omega)$ and the charge density $\rho$ be a given $C^{1}$ curve in $L^{2}(\Omega)$, satisfying the charge conservation condition $\dot{\rho}+\operatorname{div}J=0$.

2.2.2. The Lagrangian and Euler–Lagrange equations

For $(A,\varphi)$ as above, define the Lagrangian

[TABLE]

The Euler–Lagrange equations are

[TABLE]

which are weak expressions of (1b) and (2b), respectively.

These Euler–Lagrange equations imply that solutions have additional regularity properties. Since $\operatorname{curl}H=\dot{{D}}+J$ is $C^{0}$ in $L^{2}$, we have that $H$ is $C^{0}$ in $H(\operatorname{curl};\Omega)$. Likewise, since $\operatorname{div}{D}=\rho$ is $C^{1}$ in $L^{2}$, we have that ${D}$ is $C^{1}$ in $H(\operatorname{div};\Omega)$. Hence, solutions to this weak problem are in fact strong solutions of Maxwell’s equations.

Remark \theremark.

When $\epsilon$ and $\mu$ are constant in time, the electric and magnetic fields have precisely the same regularity assumed by Monk [27, eqs. (7)–(8)], namely: $E$ is $C^{1}$ in $L^{2}(\Omega,\mathbb{R}^{3})$ and $C^{0}$ in $\mathring{H}(\operatorname{curl};\Omega)$, while $H$ is $C^{1}$ in $L^{2}(\Omega,\mathbb{R}^{3})$ and $C^{0}$ in $H(\operatorname{curl};\Omega)$.

As in Section 2.1, this formulation is symmetric with respect to gauge transformations $(A,\varphi)\mapsto(A+\operatorname{grad}\xi,\varphi-\dot{\xi})$, where $\xi$ is now an arbitrary $C^{1}$ curve in $\mathring{H}^{1}(\Omega)$. Fixing the temporal gauge $\varphi=0$, the Lagrangian becomes

[TABLE]

and the Euler–Lagrange equations are just (3a). This again implies that $H$ is $C^{0}$ in $H(\operatorname{curl};\Omega)$, so (1b) holds strongly. By the same argument as in Section 2.1, this automatically preserves the charge-conservation constraint.

Remark \theremark.

Preservation of the charge-conservation constraint may also be seen as a consequence of the remaining gauge symmetry $A\mapsto A+\operatorname{grad}\xi$, mentioned in Section 2.1.2, where $\xi\in\mathring{H}^{1}(\Omega)$ is constant in time. This is a particular instance of Noether’s theorem, which relates symmetries to conservation laws. See Marsden and Ratiu [24, Section 1.6] for an account of the $J=0$ case, as well as the discussion in Christiansen and Winther [15].

2.3. Galerkin semidiscretization using Nédélec elements

The use of finite elements in computational electromagnetics is a broad topic with a long history, and we do not attempt to give a full account here. We refer the reader to the texts by Monk [28] and Jin [21], as well as the excellent survey article by Hiptmair [19], which relates these methods to the more recent theory of finite element spaces of differential forms. In this section, we briefly review the semidiscretization of Maxwell’s equations using the elements of Nédélec [29, 30], an approach that was subsequently analyzed in a series of papers by Monk [25, 26, 27].

Galerkin semidiscretization of the variational problem (3a) restricts the trial and test functions to some finite-dimensional subspace $V_{h}^{1}\subset\mathring{H}(\operatorname{curl};\Omega)$, resulting in a finite-dimensional system of ODEs. That is, we seek a $C^{1}$ curve $A_{h}\colon t\mapsto A_{h}(t)\in V_{h}^{1}$ such that

[TABLE]

where $E_{h}\coloneqq-\dot{A}_{h}$, $B_{h}\coloneqq\operatorname{curl}A_{h}$, ${D}_{h}\coloneqq\epsilon E_{h}$, and $H_{h}\coloneqq\mu^{-1}B_{h}$. The discrete versions of (1a) and (2a),

[TABLE]

follow immediately. In fact, both hold strongly in $L^{2}$, by the same argument as in Section 2.2.1, since $E_{h}\in V_{h}^{1}\subset\mathring{H}(\operatorname{curl};\Omega)$ and $B_{h}\in\operatorname{curl}V_{h}^{1}\subset\mathring{H}(\operatorname{div};\Omega)$. On the other hand, we cannot conclude that ${D}_{h}$ is in $H(\operatorname{div};\Omega)$, nor that $H_{h}$ is in $H(\operatorname{curl};\Omega)$, since (4) only holds for test functions in $V_{h}^{1}$ and not all of $\mathring{H}(\operatorname{curl};\Omega)$.

Consequently, the charge-conservation constraint (2b) is only preserved in the following, much weaker sense. Let $V_{h}^{0}\subset\mathring{H}^{1}(\Omega)$ be a finite-dimensional subspace such that $\operatorname{grad}V_{h}^{0}\subset V_{h}^{1}$. Then, for all $\xi_{h}\in V_{h}^{0}$, taking $A^{\prime}_{h}=\operatorname{grad}\xi_{h}$ in (4) and applying $\dot{\rho}+\operatorname{div}J=0$ gives

[TABLE]

Hence, if the initial conditions satisfy $\int_{\Omega}(\operatorname{grad}\xi_{h}\cdot{D}_{h}+\xi_{h}\rho)=0$, for all $\xi_{h}\in V_{h}^{0}$, then this is preserved by the flow of (4).

In particular, suppose now that $\Omega$ is polyhedral, and that $\mathcal{T}_{h}$ is a triangulation of $\Omega$ by $3$-simplices (i.e., tetrahedra) $K\in\mathcal{T}_{h}$. We may take $V_{h}^{0}$ to be the space of continuous degree-$r$ piecewise polynomials on $\mathcal{T}_{h}$ vanishing on $\partial\Omega$, corresponding to standard Lagrange finite elements. For $V_{h}^{1}$, we may take either degree-$r$ Nédélec edge elements of the first kind [29] or degree-$(r-1)$ Nédélec edge elements of the second kind [30] with vanishing degrees of freedom on $\partial\Omega$. These are spaces of piecewise-polynomial vector fields in $\mathbb{R}^{3}$ with tangential (but not necessarily normal) continuity between neighboring simplices. These choices ensure that $\operatorname{grad}V_{h}^{0}\subset V_{h}^{1}$, so the weak charge-conservation argument above holds.

Note, however, that $\int_{\Omega}(\operatorname{grad}\xi_{h}\cdot{D}_{h}+\xi_{h}\rho)=0$ only says that $\operatorname{div}{D}_{h}=\rho$ holds in an “averaged” sense, since (unlike in the infinite-dimensional case) nonzero $\xi_{h}\in V_{h}^{0}$ cannot be taken to have arbitrarily small support. We cannot even conclude that the constraint holds in the sense that $\int_{\partial K}{D}_{h}\cdot\mathbf{n}=\int_{K}\rho$, since the indicator function $\mathbf{1}_{K}$ is discontinuous and therefore not an admissible test function. (Christiansen and Winther [15] give a compactness argument for why this weak form of the constraint “might be just as good” as the strong form, in the limit as $h\rightarrow 0$; see also Christiansen [14].) This motivates our proposed hybrid approach, based on domain decomposition, for which piecewise-constants are admissible test functions.

Remark \theremark.

The method above describes the evolution of $A_{h}\in V_{h}^{1}$. Equivalently, one may evolve $E_{h}\in V_{h}^{1}$ and $B_{h}\in\operatorname{curl}V_{h}^{1}\subset V_{h}^{2}\subset\mathring{H}(\operatorname{div};\Omega)$ by augmenting (4) with $\dot{B}_{h}=-\operatorname{curl}E_{h}$. This is the original approach described by Nédélec [29], where $V_{h}^{2}$ is given by face elements on $\mathcal{T}_{h}$.

3. Domain decomposition

In this section, we introduce an alternative variational formulation for Maxwell’s equations, based on domain decomposition. Specifically, we decompose the problem on $\Omega$ into a collection of problems on $K\in\mathcal{T}_{h}$, weakly enforcing internal continuity and external boundary conditions using Lagrange multipliers. This is similar in spirit to the standard approach to domain decomposition for Poisson’s equation, cf. Brezzi and Fortin [12]. We show that the Lagrange multipliers enforcing these conditions on $A$ and $\varphi$ correspond to the traces of $H$ and $D$, respectively, and we show that the latter satisfies an appropriate version of the charge-conservation constraint.

3.1. Function spaces

We begin by introducing the following discontinuous function spaces, which are larger than the spaces used in the previous variational formulation:

[TABLE]

Brezzi and Fortin [12, Proposition III.1.1] show that

[TABLE]

That is, $\mathring{H}^{1}(\Omega)$ is the subspace of $H^{1}(\mathcal{T}_{h})$ where internal continuity and external boundary conditions are enforced by Lagrange multipliers $\lambda\in H(\operatorname{div};\Omega)$. Likewise, [12, Proposition III.1.2] shows that

[TABLE]

Using a similar argument, we now prove the corresponding result for the $H(\operatorname{curl})$ spaces.

Proposition \theproposition.

\mathring{H}(\operatorname{curl};\Omega)=\bigl{\{}u\in H(\operatorname{curl};\mathcal{T}_{h}):\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}(u\times\lambda)\cdot\mathbf{n}=0,\text{ for all }\lambda\in H(\operatorname{curl};\Omega)\bigr{\}}.

Proof.

If $u\in\mathring{H}(\operatorname{curl};\Omega)\subset H(\operatorname{curl};\mathcal{T}_{h})$, then for any $\lambda\in H(\operatorname{curl};\Omega)$, we have

[TABLE]

so the forward inclusion $(\subset)$ holds. To get the reverse inclusion $(\supset)$, suppose that $u\in H(\operatorname{curl};\mathcal{T}_{h})$ satisfies the condition above, and let $\lambda\in C^{\infty}_{c}(\Omega,\mathbb{R}^{3})$. Then, integrating by parts, we have

[TABLE]

where the last line uses the triangle and Cauchy–Schwarz inequalities. It follows that $\operatorname{curl}u\in L^{2}(\Omega,\mathbb{R}^{3})$, so $u\in H(\operatorname{curl};\Omega)$. This implies that $\int_{\partial\Omega}(u\times\lambda)\cdot\mathbf{n}=\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}(u\times\lambda)\cdot\mathbf{n}=0$ for all $\lambda\in H(\operatorname{curl};\Omega)$. Hence, $u\times\mathbf{n}\rvert_{\partial\Omega}=0$ in the trace sense, which completes the proof. ∎

Remark \theremark.

A variant of this result is stated in Boffi et al. [7, Proposition 2.1.3], where $\lambda$ is taken to be in $H^{1}(\Omega,\mathbb{R}^{3})$ rather than $H(\operatorname{curl};\Omega)$. However, the version given here is more natural for the purposes of the hybrid methods discussed in Section 4.

3.2. Domain decomposition of the Lagrangian variational principle

We now introduce a new Lagrangian for Maxwell’s equations, which allows the potentials to live in the discontinuous function spaces defined in the previous section, enforcing continuity and boundary conditions using Lagrange multipliers.

Let $A(t)\in H(\operatorname{curl};\mathcal{T}_{h})$ and $\varphi(t)\in H^{1}(\mathcal{T}_{h})$, and introduce the Lagrange multipliers $\widehat{{H}}(t)\in H(\operatorname{curl};\Omega)$ and $\widehat{{D}}(t)\in H(\operatorname{div};\Omega)$. We adopt the notation, often seen in the literature on discontinuous Galerkin and hybrid methods, of placing hats over variables that act like weak traces/fluxes. As before, suppose that $t\mapsto A(t)$ is $C^{1}$ and that $t\mapsto\varphi(t)$ is $C^{0}$, such that $t\mapsto{D}(t)\in L^{2}(\Omega,\mathbb{R}^{3})$ is $C^{1}$. Furthermore, suppose that $t\mapsto\widehat{H}(t)$ and $t\mapsto\widehat{{D}}(t)$ are both $C^{0}$. Define the Lagrangian

[TABLE]

The Euler–Lagrange equations are then

[TABLE]

where (5a) and (5b) hold for all $K\in\mathcal{T}_{h}$. We now relate this to the classical variational form of Maxwell’s equations, stated in (3).

Proposition \theproposition.

$(A,\varphi,\widehat{{H}},\widehat{{D}})$* is a solution to (5) if and only if $(A,\varphi)$ is a solution to (3) with $\widehat{H}\times\mathbf{n}\rvert_{\partial K}=H\times\mathbf{n}\rvert_{\partial K}$ and $\widehat{{D}}\cdot\mathbf{n}\rvert_{\partial K}={D}\cdot\mathbf{n}\rvert_{\partial K}$. In particular, if $(A,\varphi)$ is a solution to (3), then $(A,\varphi,H,{D})$ is a solution to (5).*

Proof.

Suppose $(A,\varphi,\widehat{{H}},\widehat{{D}})$ is a solution to (5). By Section 3.1, (5c) implies $A(t)\in\mathring{H}(\operatorname{curl};\Omega)$, so taking $A^{\prime}\in\mathring{H}(\operatorname{curl};\Omega)$ and summing (5a) over $K\in\mathcal{T}_{h}$, the integrals over $\partial K$ cancel, yielding (3a). As previously stated, (3a) implies $\operatorname{curl}H=\dot{{D}}+J$, so substituting this into (5a) gives

[TABLE]

so $\widehat{H}\times\mathbf{n}\rvert_{\partial K}=H\times\mathbf{n}\rvert_{\partial K}$. Similarly, (5d) implies $\varphi(t)\in\mathring{H}^{1}(\Omega)$, so taking $\varphi^{\prime}\in\mathring{H}^{1}(\Omega)$ and summing (5b) over $K\in\mathcal{T}_{h}$ yields (3b). This implies $\operatorname{div}{D}=\rho$, and substituting into (5b) gives $\widehat{{D}}\cdot\mathbf{n}\rvert_{\partial K}={D}\cdot\mathbf{n}\rvert_{\partial K}$.

Conversely, suppose $(A,\varphi)$ is a solution to (3). Since $A(t)\in\mathring{H}(\operatorname{curl};\Omega)$ and $\varphi(t)\in\mathring{H}^{1}(\Omega)$, it follows that (5c) and (5d) hold. Furthermore, (3) implies that $\dot{{D}}+J=\operatorname{curl}H$ and $\operatorname{div}{D}=\rho$, so (5a) and (5b) hold with $\widehat{{H}}\times\mathbf{n}\rvert_{\partial K}=H\times\mathbf{n}\rvert_{\partial K}$ and $\widehat{{D}}\cdot\mathbf{n}\rvert_{\partial K}={D}\cdot\mathbf{n}\rvert_{\partial K}$. In particular, we could take $\widehat{H}=H$ and $\widehat{{D}}={D}$. ∎

Remark \theremark.

Note that, in addition to (5b) implying that $\operatorname{div}{D}=\rho$, we also see by taking $\varphi^{\prime}=\mathbf{1}_{K}$ that $\widehat{{D}}$ satisfies the conservation law $\int_{\partial K}\widehat{{D}}\cdot\mathbf{n}=\int_{K}\rho$, for all $K\in\mathcal{T}_{h}$.

3.3. Temporal gauge fixing and the charge-conservation

constraint

As in Section 2.2, if $(A,\varphi,\widehat{H},\widehat{{D}})$ is a solution to (5), then so is $(A+\operatorname{grad}\xi,\varphi-\dot{\xi},\widehat{H},\widehat{{D}})$ for any $C^{1}$ curve $t\mapsto\xi(t)\in\mathring{H}^{1}(\Omega)$. Therefore, we perform temporal gauge fixing by taking $\varphi=0$. This yields the gauge-fixed Lagrangian

[TABLE]

whose Euler–Lagrange equations are simply (5a) and (5c). Of course, (5d) is satisfied trivially, since $\varphi=0$. The next result shows that the charge-conservation constraint (5b) is automatically preserved, for an appropriately-defined $\widehat{{D}}$.

Proposition \theproposition.

Let $(A,\widehat{H})$ be a solution to (5a) and (5c). Suppose initial values for ${D},\widehat{{D}}$ satisfy (5b), and let $\widehat{{D}}$ be the solution to $\dot{\widehat{{D}}}+J=\operatorname{curl}\widehat{H}$. Then $(A,0,\widehat{H},\widehat{{D}})$ is a solution to (5).

Proof.

As we have already mentioned, $\varphi=0$ trivially satisfies (5d), so it suffices to show that (5b) holds. Let $\varphi^{\prime}\in H^{1}(K)$ be arbitrary. Taking $A^{\prime}=\operatorname{grad}\varphi^{\prime}$ in (5a) and integrating by parts gives

[TABLE]

so if (5b) holds at the initial time, then it holds for all time. ∎

Remark \theremark.

As in Section 3.2, taking $\varphi^{\prime}=\mathbf{1}_{K}$ implies $\int_{\partial K}\widehat{{D}}\cdot\mathbf{n}=\int_{K}\rho$. Furthermore, if the initial conditions also satisfy $\operatorname{div}\widehat{{D}}=\rho$, then we have $\operatorname{div}\widehat{{D}}=\rho$ for all time, since $\operatorname{div}\dot{\widehat{{D}}}=\operatorname{div}\operatorname{curl}\widehat{H}-\operatorname{div}J=0+\dot{\rho}$. Finally, if $\widehat{H}=H$, and if the initial conditions for $\widehat{{D}}$ equal those for ${D}$, then we recover $\widehat{{D}}={D}$.

Finally, we express this variational problem in the standard notation used for mixed and hybrid finite element methods, in terms of a pair of bilinear forms [12, Chapter II]. We will make use of this notation throughout the subsequent sections. Defining

[TABLE]

we seek $t\mapsto A(t)\in H(\operatorname{curl};\mathcal{T}_{h})$ and $t\mapsto\widehat{H}(t)\in H(\operatorname{curl};\Omega)$ such that

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the $L^{2}(\Omega,\mathbb{R}^{3})$ inner product. Defining the map $\mathcal{B}\colon H(\operatorname{curl};\mathcal{T}_{h})\rightarrow H(\operatorname{curl};\Omega)^{\ast}$, $A\mapsto b(A,\cdot)$, we see that (6) is equivalent to evolving $A(t)\in\ker\mathcal{B}$ by

[TABLE]

and subsequently solving for $\widehat{H}$ satisfying (6a). Since $\ker\mathcal{B}=\mathring{H}(\operatorname{curl};\Omega)$ by Section 3.1, it follows that $A$ solves the non-domain-decomposed problem (3a).

4. Hybrid semidiscretization

We now perform Galerkin semidiscretization of the domain-decomposed variational problem with temporal gauge fixing, as introduced in the previous section. This results in a hybrid method for Maxwell’s equations, where “hybrid” means that the Lagrange multipliers $\widehat{H}_{h}$ and their test functions $\widehat{H}_{h}^{\prime}$ are both restricted to a subspace of $H(\operatorname{curl};\Omega)$. We then show that a suitably-defined $\widehat{{D}}_{h}$ satisfies the charge-conservation constraint in a strong sense, as opposed to the much weaker sense in which ${D}_{h}$ was seen to satisfy this constraint in Section 2.3. Finally, we discuss how certain choices of elements yield a hybridized version of Nédélec’s method, while others give nonconforming methods, and we remark on how this framework also applies to the 2-D Maxwell equations.

4.1. Semidiscretization of the variational problem

For each $K\in\mathcal{T}_{h}$, let $V_{h}^{1}(K)\subset H(\operatorname{curl};K)$ be a finite-dimensional subspace, so $V_{h}^{1}\coloneqq\prod_{K\in\mathcal{T}_{h}}V_{h}^{1}(K)\subset H(\operatorname{curl};\mathcal{T}_{h})$, and let $\widehat{V}_{h}^{1}\subset H(\operatorname{curl};\Omega)$. We seek $A_{h}\colon t\mapsto A_{h}(t)\in V_{h}^{1}$ and $\widehat{H}_{h}\colon t\mapsto\widehat{H}_{h}(t)\in\widehat{V}_{h}^{1}$ such that

[TABLE]

where (7a) holds for all $K\in\mathcal{T}_{h}$. These are the semidiscretized versions of (5a) and (5c).

Remark \theremark.

Since (7b) only holds for test functions in $\widehat{V}_{h}^{1}$, but not necessarily an arbitrary test function in $H(\operatorname{curl};\Omega)$, in general a solution will have $A_{h}(t)\notin\mathring{H}(\operatorname{curl};\Omega)$. Hence, this method is generally not curl-conforming and is distinct from the conforming methods discussed in Section 2.3.

In terms of the bilinear forms $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$, this method may be written as

[TABLE]

Defining the operator $\mathcal{B}_{h}\colon V_{h}^{1}\rightarrow(\widehat{V}_{h}^{1})^{\ast}$, $A_{h}\mapsto b(A_{h},\cdot)\rvert_{\widehat{V}_{h}^{1}}$, we see that (8) is equivalent to evolving $A_{h}(t)\in\ker\mathcal{B}_{h}$ by

[TABLE]

and subsequently solving for $\widehat{H}_{h}$ satisfying (8a).

Since $V_{h}^{1}$ is finite-dimensional, we may apply Banach’s closed range theorem to deduce that $\langle\dot{{D}}_{h}+J,\cdot\rangle-a(A_{h},\cdot)\in(\ker\mathcal{B}_{h})^{\perp}$ is in the range of $\mathcal{B}_{h}^{\ast}$, so a solution $\widehat{H}_{h}$ exists, although generally not uniquely. A natural choice is to find the solution $\widehat{H}_{h}$ minimizing $\lVert H_{h}-\widehat{H}_{h}\rVert^{2}+\lVert\dot{{D}}_{h}+J-\operatorname{curl}\widehat{H}_{h}\rVert^{2}$, which in a weak sense minimizes the $H(\operatorname{curl};\Omega)$ distance between $H_{h}$ and $\widehat{H}_{h}$. This existence-without-uniqueness is typical of hybrid methods, and one may formally resolve this by replacing $\widehat{V}_{h}^{1}$ by the quotient space $\widehat{V}_{h}^{1}/\ker\mathcal{B}_{h}^{\ast}$ (cf. Brezzi and Fortin [12, IV.1.3]). In practice, the evolution on $\ker\mathcal{B}_{h}$ specified by (9) is the essence of the method, and solving for $\widehat{H}_{h}$ may be seen as an optional post-processing step.

4.2. Preservation of the charge-conservation constraint

In order to discuss the charge-conservation constraint, we first suppose that $V_{h}^{0}(K)\subset H^{1}(K)$ are such that $\mathbf{1}_{K}\in V_{h}^{0}(K)$ and $\operatorname{grad}V_{h}^{0}(K)\subset V_{h}^{1}(K)$ for all $K\in\mathcal{T}_{h}$. We consider whether the following discretization of (5b) is preserved,

[TABLE]

for $\widehat{{D}}_{h}\colon t\mapsto\widehat{{D}}_{h}(t)\in H(\operatorname{div};\Omega)$ suitably defined.

Theorem 4.1**.**

Let $(A_{h},\widehat{H}_{h})$ be a solution to (7). Suppose initial values for ${D}_{h},\widehat{{D}}_{h}$ satisfy (10), and let $\widehat{{D}}_{h}$ be the solution to $\dot{\widehat{{D}}}_{h}+J=\operatorname{curl}\widehat{H}_{h}$. Then (10) holds for all time. In particular, $\int_{\partial K}\widehat{{D}}_{h}\cdot\mathbf{n}=\int_{K}\rho$. Moreover, if $\operatorname{div}\widehat{{D}}_{h}=\rho$ holds at the initial time, then it holds for all time.

Proof.

The proof is essentially similar to that of Section 3.3. Given $\varphi_{h}^{\prime}\in V_{h}^{0}(K)$, taking $A_{h}^{\prime}=\operatorname{grad}\varphi_{h}^{\prime}\in V_{h}^{1}(K)$ in (7a) and integrating by parts,

[TABLE]

so if (10) holds at the initial time, then it holds for all time. The conclusion that $\int_{\partial K}\widehat{{D}}_{h}\cdot\mathbf{n}=\int_{K}\rho$ follows by taking $\varphi_{h}^{\prime}=\mathbf{1}_{K}$, and $\operatorname{div}\dot{\widehat{{D}}}_{h}=\operatorname{div}\operatorname{curl}\widehat{H}_{h}-\operatorname{div}J=0+\dot{\rho}$ implies that if $\operatorname{div}\widehat{{D}}_{h}=\rho$ holds at the initial time, then it holds for all time. ∎

Remark \theremark.

Preservation of $\operatorname{div}\widehat{{D}}_{h}=\rho$ is immediate from $\dot{\widehat{{D}}}_{h}+J=\operatorname{curl}\widehat{H}_{h}$, without appealing to (10). However, it is only a meaningful statement about solutions to (7) when (10) holds. By contrast, if $\widehat{{D}}_{h}$ were instead to satisfy $\dot{\widehat{{D}}}_{h}+J=0$, then $\operatorname{div}\widehat{{D}}_{h}=\rho$ would still be preserved, but this would not say anything about the numerical solution $(A_{h},\widehat{H}_{h})$.

The next result addresses the existence of initial conditions for $\widehat{{D}}_{h}$ satisfying the hypotheses of the previous theorem. Let $V_{h}^{0}\coloneqq\prod_{K\in\mathcal{T}_{h}}V_{h}^{0}(K)\subset H^{1}(\mathcal{T}_{h})$.

Proposition \theproposition.

Suppose that the initial value of ${D}_{h}$ satisfies

[TABLE]

Then there exists an initial value for $\widehat{{D}}_{h}$ such that (10) holds for all $K\in\mathcal{T}_{h}$ and $\operatorname{div}\widehat{{D}}_{h}=\rho$.

Proof.

The first part of the argument is similar to the one we used for the existence of $\widehat{H}_{h}$. Define the bilinear form

[TABLE]

Since $V_{h}^{0}\cap\mathring{H}^{1}(\Omega)$ is the kernel of $\varphi_{h}^{\prime}\mapsto\beta_{h}(\varphi_{h}^{\prime},\cdot)$, the closed range theorem implies that $\varphi_{h}^{\prime}\mapsto\sum_{K\in\mathcal{T}_{h}}\int_{K}\operatorname{grad}\varphi_{h}^{\prime}\cdot{D}_{h}+\int_{\Omega}\varphi_{h}^{\prime}\rho$ is in the range of $\widehat{{D}}_{h}\mapsto\beta_{h}(\cdot,\widehat{{D}}_{h})$. Hence, there exists an initial value for $\widehat{{D}}_{h}$ satisfying (10) for all $K\in\mathcal{T}_{h}$.

Next, suppose $\widehat{{D}}_{h}$ satisfies (10) but not necessarily $\operatorname{div}\widehat{{D}}_{h}=\rho$. Then, on each $K\in\mathcal{T}_{h}$, replace $\widehat{{D}}_{h}$ by $\widehat{{D}}_{h}+\operatorname{grad}u$, where $u$ is the solution to $-\Delta u=\operatorname{div}\widehat{{D}}_{h}-\rho$ with Neumann boundary conditions $\operatorname{grad}u\cdot\mathbf{n}=0$ on $\partial K$. This leaves the normal traces of $\widehat{{D}}_{h}$ unchanged, so the result is still in $H(\operatorname{div};\Omega)$ and satisfies (10), as desired. ∎

Remark \theremark.

The computation of $\widehat{{D}}_{h}$, like that of $\widehat{H}_{h}$, can be seen as an optional post-processing step after computing the solution $A_{h}$ to (9). The key point of Theorem 4.1 is that the evolution of $A_{h}$ is conservative, in the sense that it is consistent with a charge-conserving numerical flux $\widehat{{D}}_{h}$, whether or not one chooses to actually compute $\widehat{{D}}_{h}$.

4.3. Hybridization of Nédélec’s method and nonconforming

methods

As in Section 2.3, let $\Omega$ be polyhedral and $\mathcal{T}_{h}$ be a simplicial triangulation. Let $V_{h}^{0}(K)$ be the space of degree-$r$ polynomials on $K$ and $V_{h}^{1}(K)$ be either degree-$r$ Nédélec edge elements of the first kind or degree-$(r-1)$ Nédélec edge elements of the second kind on $K$. Then $V_{h}^{0}\subset H^{1}(\mathcal{T}_{h})$ and $V_{h}^{1}\subset H(\operatorname{curl};\mathcal{T}_{h})$ correspond to discontinuous Lagrange and Nédélec elements, respectively. Note that discontinuous Nédélec elements of the second kind are just discontinuous piecewise polynomial vector fields.

Now, taking $\widehat{V}_{h}^{1}=H(\operatorname{curl};\Omega)$, it follows that $\ker\mathcal{B}_{h}=V_{h}^{1}\cap\ker\mathcal{B}\subset\mathring{H}(\operatorname{curl};\Omega)$, which corresponds precisely to curl-conforming Nédélec elements with tangential inter-element continuity and boundary conditions. It follows that (9) agrees precisely with Nédélec’s method (4). In fact, it is not necessary to take $\widehat{V}_{h}^{1}$ infinite-dimensional: it suffices to take a large enough finite-dimensional subspace (e.g., Nédélec elements of sufficiently high degree) such that (7b) imposes all the inter-element continuity and boundary conditions on degrees of freedom of $V_{h}^{1}$. (Having $\widehat{V}_{h}^{1}$ infinite-dimensional is not a problem if one is only interested in $A_{h}$, but a finite-dimensional subspace is required if one wishes to compute $\widehat{H}_{h}$.) From these observations, we obtain the following corollary of Theorem 4.1 and Section 4.2

Corollary \thecorollary.

Given $V_{h}^{0}$ and $V_{h}^{1}$ as above, there exists $\widehat{V}_{h}^{1}$ such that solutions $A_{h}$ to Nédélec’s method (4) are equivalent to solutions $(A_{h},\widehat{H}_{h})$ to the hybrid method (7). Consequently, given a solution to Nédélec’s method, there exists $\widehat{{D}}_{h}$ satisfying $\dot{\widehat{{D}}}_{h}+J=\operatorname{curl}\widehat{H}_{h}$, which preserves the charge-conservation constraints (10) and $\operatorname{div}\widehat{{D}}_{h}=\rho$.

In contrast, if $\widehat{V}_{h}^{1}$ is not sufficiently large, we will have $\ker\mathcal{B}_{h}\not\subset\ker\mathcal{B}=\mathring{H}(\operatorname{curl};\Omega)$, so (9) is a nonconforming finite element method for Maxwell’s equations.

4.4. Remarks on the two-dimensional case

This framework may also be adapted to two-dimensional electromagnetics with minor modifications.

For the non-domain-decomposed problem on $\Omega\subset\mathbb{R}^{2}$, the potential $A\in\mathring{H}(\operatorname{curl};\Omega)$ remains a vector field, although $\operatorname{curl}A\in L^{2}(\Omega)$ becomes a scalar field. Consequently, $E$ and $D$ remain vector fields (and $\epsilon$ remains a tensor), while $B$ and $H$ become scalar fields (and $\mu$ becomes scalar). The two-dimensional version of the weak problem (3a) is nearly identical, except the dot product $\operatorname{curl}A^{\prime}\cdot H$ is replaced by the ordinary product $(\operatorname{curl}A^{\prime})H$. For the Galerkin semidiscretization discussed in Section 2.3, one simply replaces the Nédélec edge elements of the first and second kind with Raviart–Thomas (RT) [31] and Brezzi–Douglas–Marini (BDM) [11] edge elements, respectively. These two-dimensional $H(\operatorname{curl})$ elements are just the RT and BDM $H(\operatorname{div})$ elements rotated by 90 degrees, so that tangential traces of the former correspond to normal traces of the latter.

For domain decomposition, Section 3.1 is easily modified to show that

[TABLE]

Alternatively, this can be seen to follow from the corresponding result for $\mathring{H}(\operatorname{div};\Omega)$, where the vector fields are rotated by 90 degrees. Hence, the domain decomposed variational problem in temporal gauge is to find $t\mapsto A(t)\in H(\operatorname{curl};\mathcal{T}_{h})$ and $t\mapsto\widehat{H}(t)\in H^{1}(\Omega)$ such that

[TABLE]

for all $K\in\mathcal{T}_{h}$. Hybrid methods may then be obtained by restricting this variational problem to subspaces $V_{h}^{1}=\prod_{K\in\mathcal{T}_{h}}V_{h}^{1}(K)\subset H(\operatorname{curl};\mathcal{T}_{h})$ and $\widehat{V}_{h}^{0}\subset H^{1}(\Omega)$. As in Section 4.2, one obtains $\widehat{{D}}_{h}(t)\in H(\operatorname{div};\Omega)$ by solving $\dot{\widehat{{D}}}_{h}+J=\operatorname{curl}\widehat{H}_{h}$ (where the curl of a scalar field is its gradient rotated by 90 degrees, i.e., $v\cdot\operatorname{curl}\widehat{H}_{h}\coloneqq v\times\operatorname{grad}\widehat{H}_{h}$ for $v\in\mathbb{R}^{2}$), and the charge-conserving properties follow in the same manner.

For the finite element spaces, one may take $V_{h}^{0}$ to be discontinuous degree-$r$ Lagrange elements and $V_{h}^{1}$ to be discontinuous degree-$r$ RT edge elements or discontinuous degree-$(r-1)$ BDM edge elements. (Note that discontinuous BDM elements are just discontinuous piecewise polynomial vector fields.) In this case, it is much easier to see which $\widehat{V}_{h}^{0}\subset H^{1}(\Omega)$ yield conforming methods, since each edge degree of freedom is either shared by exactly two triangles or lies on the boundary. Both the degree-$r$ RT and degree-$(r-1)$ BDM elements have $r$ degrees of freedom per edge, which match up precisely with those for degree-$(r+1)$ Lagrange elements. Hence, taking $\widehat{V}_{h}^{0}$ corresponding to degree-$(r+1)$ or higher Lagrange elements yields a conforming method. On the other hand, a straightforward counting argument shows that degree-$r$ Lagrange elements have fewer than $r\times\text{\#edges}$ degrees of freedom on element boundaries (unless $\mathcal{T}_{h}$ consists of a single triangle). Since it is impossible to enforce all of the inter-element and boundary conditions in this case, the resulting method is nonconforming.

5. Numerical examples

This section gives numerical illustrations for the simple test problem

[TABLE]

which corresponds to the case where $\epsilon$ and $\mu$ are positive constants with $\epsilon\mu=1$ and $J=0$, as discussed at the end of Section 2.1. As before, $A$ is taken to have vanishing tangential component on the boundary. Preservation of the charge-conservation constraint is equivalent to the condition $\operatorname{div}\ddot{A}=0$.

In the frequency domain, denoting angular frequency by $\omega$, time differentiation becomes multiplication by $i\omega$, so (11) becomes the eigenvalue problem

[TABLE]

In this setting, preservation of the charge-conservation constraint becomes $\omega^{2}\operatorname{div}A=0$, i.e., eigenfunctions with nonzero eigenvalue are divergence-free.

The examples below demonstrate the constraint-preserving properties of the curl-conforming hybridized Nédélec method from Section 4, both in the time domain and in the frequency domain. In the frequency domain, we also observe superconvergence of $\widehat{H}_{h}\rightarrow H$. All finite element computations were performed using FEniCS [23, 1]. For the post-processing step of computing $\widehat{H}_{h}$, whose solution is not unique, we find the solution $\widehat{H}_{h}$ minimizing $\lVert H_{h}-\widehat{H}_{h}\rVert^{2}+\lVert\dot{{D}}_{h}+J-\operatorname{curl}\widehat{H}_{h}\rVert^{2}$, as previously discussed in Section 4.1.

5.1. Time domain

Before turning our attention to the test problem (11), we first describe a discrete time-stepping scheme for the general case of Maxwell’s equations. After semidiscretizing using the hybridized Nédélec method of Section 4, we discretize in time using the following explicit “leapfrog” scheme:

• •

$A_{n+1/2}=A_{n}-\frac{1}{2}\Delta t\epsilon^{-1}{D}_{n}$.

• •

${D}_{n+1}={D}_{n}+\Delta t\dot{{D}}_{n+1/2}$, where $\dot{D}_{n+1/2}\in\ker\mathcal{B}_{h}$ is the solution to

[TABLE]

• •

$\widehat{{D}}_{n+1}=\widehat{{D}}_{n}+\Delta t(\operatorname{curl}\widehat{H}_{n+1/2}-J_{n+1/2})$, where $\widehat{H}_{n+1/2}$ is the solution to

[TABLE]

minimizing $\lVert H_{n+1/2}-\widehat{H}_{n+1/2}\rVert^{2}+\lVert\dot{{D}}_{n+1/2}+J_{n+1/2}-\operatorname{curl}\widehat{H}_{n+1/2}\rVert^{2}$.

• •

$A_{n+1}=A_{n+1/2}-\frac{1}{2}\Delta t\epsilon^{-1}{D}_{n+1}$.

Here, $A_{n}$ denotes the approximation to $A_{h}(t_{n})$, where $t_{n}$ is the $n$th time step and $\Delta t$ is the time step size; similar notation is used for the other variables. This is essentially the Störmer/Verlet method for the semidiscretized system of ODEs (9), augmented by a hybrid post-processing step for $\widehat{H}_{h}$ and $\widehat{{D}}_{h}$. Except for the hybrid post-processing step, which is novel, such leapfrog schemes are widely used for both finite element and finite difference time domain methods in computational electromagnetics (see Yee [32] and Monk [25, Section 5]). The Störmer/Verlet method also has particularly desirable properties when applied to Lagrangian and Hamiltonian dynamics (cf. Hairer et al. [17, 18]).

Figure 1 shows the results of applying this method to the test problem (11) on the 2-D square $\Omega=(0,\pi)^{2}$ and 3-D cube $\Omega=(0,\pi)^{3}$, taking $\epsilon=\mu=1$. For both the 2-D and 3-D problems, we simulate over $t\in[0,2\pi]$ for $1024$ time steps of size $\Delta t=\pi/512$.

For the 2-D problem, the initial conditions are taken to be ${D}_{0}=\widehat{{D}}_{0}=0$ and

[TABLE]

A uniform triangular mesh is taken on a $16\times 16$ grid, with $2\cdot 16^{2}=512$ cells. The space $V_{h}^{1}$ consists of discontinuous piecewise linear vector fields, while $\widehat{V}_{h}^{0}$ consists of cubic Lagrange elements, so that $\ker\mathcal{B}_{h}\subset V_{h}^{1}$ are linear BDM edge elements, as described in Section 4.4 with $r=2$.

For the 3-D problem, the initial conditions are taken to be ${D}_{0}=\widehat{{D}}_{0}=0$ and

[TABLE]

A uniform tetrahedral mesh is taken on an $8\times 8\times 8$ grid, with $6\cdot 8^{3}=3072$ cells. The space $V_{h}^{1}$ consists of discontinuous piecewise linear vector fields, while $\widehat{V}_{h}^{1}$ consists of cubic Nédélec edge elements of the second kind, so that $\ker\mathcal{B}_{h}\subset V_{h}^{1}$ are linear Nédélec edge elements of the second kind, as described in Section 4.3 with $r=2$.

Although the exact solution satisfies $\operatorname{div}{D}=0$, the numerical solution ${D}_{h}$ drifts away from this constraint, as measured by the $H(\operatorname{div};\mathcal{T}_{h})$ seminorm,

[TABLE]

However, $\operatorname{div}\widehat{{D}}_{h}=0$ holds to machine precision, as explained by Theorem 4.1. Looking at ${D}_{h}$ alone, one might think that this method fails to preserve the charge-conservation constraint strongly. In fact, we have illustrated that it actually does preserve this constraint, when expressed in terms of the numerical flux $\widehat{{D}}_{h}$ rather than ${D}_{h}$.

Remark \theremark.

The constraint behavior of ${D}_{h}$ and $\widehat{{D}}_{h}$, observed in Figure 1, is due to the finite element semidiscretization, not the time discretization. Indeed, the charge-conservation constraint is linear, so if it holds for the semidiscretized system of ODEs, then any Runge–Kutta or partitioned Runge–Kutta method preserves it (Hairer et al. [18, Theorem IV.1.2]).

5.2. Frequency domain

We next apply the hybrid approach to the frequency domain, again assuming that $\epsilon$ and $\mu$ are positive constants with $\epsilon\mu=1$ and $J=0$. This is done by first approximating the Maxwell eigenvalue problem (12) on $\ker\mathcal{B}_{h}$ and then applying hybrid post-processing, as follows:

• •

Find eigenpairs $(\omega_{h}^{2},A_{h})\in\mathbb{R}^{+}\times\ker\mathcal{B}_{h}$ satisfying

[TABLE]

and let $H_{h}\coloneqq\mu^{-1}\operatorname{curl}A_{h}$ and ${D}_{h}\coloneqq\epsilon(-i\omega_{h}A_{h})$.

• •

Find $\widehat{H}_{h}$ minimizing $\lVert H_{h}-\widehat{H}_{h}\rVert^{2}+\lVert i\omega_{h}{D}_{h}-\operatorname{curl}\widehat{H}_{h}\rVert^{2}$ such that

[TABLE]

and let $\widehat{{D}}_{h}\coloneqq-i\omega_{h}^{-1}\operatorname{curl}\widehat{H}_{h}$.

Note that this last step is equivalent to $i\omega_{h}\widehat{{D}}_{h}=\operatorname{curl}\widehat{H}_{h}$, so $\widehat{H}_{h}$ can be seen as minimizing $\lVert H_{h}-\widehat{H}_{h}\rVert^{2}+\omega_{h}^{2}\lVert{D}_{h}-\widehat{{D}}_{h}\rVert^{2}$.

We consider the 2-D square $\Omega=(0,\pi)^{2}$, where the exact eigenvalues are sums of squares ($\omega^{2}=1,1,2,4,4,\ldots$). For simplicity, we look at the approximation of the following analytical solution with simple eigenvalue $\omega^{2}=2$, assuming $\epsilon=\mu=1$:

[TABLE]

We take a uniform triangle mesh on an $N\times N$ grid, which has $2N^{2}$ cells. As described in Section 4.4, we take $V_{h}^{1}$ to consist of discontinuous piecewise degree-$(r-1)$ vector fields and $\widehat{V}_{h}^{0}$ to consist of degree-$(r+1)$ Lagrange elements, so that $\ker\mathcal{B}_{h}\subset V_{h}^{1}$ are degree-$(r-1)$ BDM edge elements.

Figure 2 shows ${D}_{h}$ and $\widehat{{D}}_{h}$, along with $\operatorname{div}{D}_{h}$ and $\operatorname{div}\widehat{{D}}_{h}$, for the case $N=16$, $r=2$. Here, by $\operatorname{div}{D}_{h}\in L^{2}(\Omega)$, we mean the element-wise divergence $(\operatorname{div}{D}_{h})\rvert_{K}\coloneqq\operatorname{div}({D}_{h}\rvert_{K})$ for each $K\in\mathcal{T}_{h}$, since ${D}_{h}$ is in $H(\operatorname{div};\mathcal{T}_{h})$ but not in $H(\operatorname{div};\Omega)$. Although the vector fields ${D}_{h}$ and $\widehat{{D}}_{h}$ appear very similar, they behave very differently with respect to the charge-conservation constraint: $\operatorname{div}{D}_{h}\neq 0$, while $\operatorname{div}\widehat{{D}}_{h}=0$ to machine precision. Note that these are purely imaginary when $A_{h}$ is real, so the imaginary parts are plotted.

Table 1 illustrates the convergence behavior of $H_{h}$, $\widehat{H}_{h}$, ${D}_{h}$, and $\widehat{{D}}_{h}$ as the mesh parameter $h\rightarrow 0$, for elements of various degrees. Since $A_{h}$ is simply obtained by using degree-$(r-1)$ BDM edge elements for the Maxwell eigenvalue problem, previous analyses of this problem (e.g., [22, 19, 5, 6, 3] and references therein) show that $\lVert A_{h}-A\rVert=\mathcal{O}(h^{r})$ and $\lVert\operatorname{curl}A_{h}-\operatorname{curl}A\rVert=\mathcal{O}(h^{r-1})$, which imply the observed rates $\lVert{D}_{h}-{D}\rVert=\mathcal{O}(h^{r})$ and $\lVert H_{h}-H\rVert=\mathcal{O}(h^{r-1})$. Interestingly, for $\widehat{H}_{h}$ obtained by hybrid post-processing, we observe the superconvergent rates $\lVert\widehat{H}_{h}-H\rVert=\mathcal{O}(h^{r})$ for $r=2$ and $\mathcal{O}(h^{r+1})$ for $r>2$. For $\widehat{{D}}_{h}$, we observe errors comparable to those for ${D}_{h}$ and the same convergence rate, $\lVert\widehat{{D}}_{h}-{D}\rVert=\mathcal{O}(h^{r})$.

We note that the observed rates of superconvergence, including the reduced rate in the lowest-degree case, are the same as those obtained for scalar elliptic problems by Brezzi, Douglas, and Marini [11] in the original paper on the hybridized BDM method.

6. Conclusion

We have constructed a family of primal hybrid finite element methods for Maxwell’s equations, where the Lagrange multipliers enforcing inter-element continuity and boundary conditions correspond to a numerical trace $\widehat{H}_{h}$ of the magnetic field and a numerical flux $\widehat{{D}}_{h}$ of the electric flux density. These methods strongly preserve the constraints $\operatorname{div}B_{h}=0$ and $\operatorname{div}\widehat{{D}}_{h}=\rho$, the latter of which corresponds to conservation of charge. As a special case, these methods include hybridized versions of standard methods using curl-conforming edge elements, which had previously been thought only to be charge-conserving in a much weaker sense. We emphasize that these conservative properties hold even if the methods are not implemented in a hybrid fashion: if desired, $\widehat{H}_{h}$ and $\widehat{{D}}_{h}$ may be recovered by an optional post-processing step.

There are several natural directions for future work. First, the numerical experiments in Section 5 focused on hybridized curl-conforming methods, due to the fact that their stability and error analysis is already well established. However, as mentioned in Section 4.3, this framework also includes constraint-preserving nonconforming methods, which would be interesting to investigate. Second, we do not yet have a complete explanation of the hybrid superconvergence phenomenon for $\widehat{H}_{h}\rightarrow H$; this is the subject of ongoing work. Finally, the techniques developed here might be applied to study constraint preservation in other families of hybrid methods, particularly hybridizable discontinuous Galerkin (HDG) methods.

Acknowledgments

Ari Stern acknowledges the support of the National Science Foundation (DMS-1913272) and the Simons Foundation (#279968). Yakov Berchenko–Kogan was supported by an AMS–Simons Travel Grant.