1. Introduction
Maxwell’s equations consist of two vector evolution equations, together with two scalar constraint equations,
and , where is magnetic flux density, is electric flux density, and 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 curlconforming edge elements for the electric field and divergenceconforming face elements for , then the semidiscretized equations preserve strongly. On the other hand, holds only in the Galerkin sense (i.e., when both sides are integrated against certain continuous, piecewisepolynomial 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 , another idea is to use divergencefree 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 and scalar potential , using Lagrange multipliers and to enforce interelement continuity and boundary conditions. These Lagrange multipliers are shown to correspond to boundary traces of the magnetic field and electric flux density . After using gauge symmetry to fix , we show that the evolution of automatically preserves the constraints and . Finally, we semidiscretize this domaindecomposed 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 (nondomaindecomposed) formulation of Maxwell’s equations, and study the domaindecomposed version of the constraints and their preservation.

In Section 4, we consider primal hybrid finite element methods for semidiscretizing the domaindecomposed 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 postprocessed magnetic field , 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,
(1a)  
(1b) 
together with the scalar constraint equations,
(2a)  
(2b) 
Here, and denote the electric field and magnetic field, and denote the electric flux density and magnetic flux density, and
are the electric permittivity and magnetic permeability tensors, and
and are current density and charge density, respectively. We use the “dot” notation to denote partial differentiation with respect to time.The evolution equations (1) automatically preserve the constraints (2). Indeed, taking the divergence of (1a) implies , so (2a) is preserved. Similarly, taking the divergence of (1b) implies , so (2b) is preserved if and only if and satisfy , which is the law of conservation of charge. We refer to (2b) as the chargeconservation 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 , called the vector potential, and a scalar field , called the scalar potential. Given and , we define the electric field and magnetic flux density by
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 . Observe that if is any timedependent scalar field, then the transformation leaves , , , 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 may be transformed into one of the form by taking to be a solution of . Therefore, we may restrict our attention to solutions with .
Remark .
This procedure of restricting to particular solutions, which are related to a general solution by some gauge transformation, is called gauge fixing. The choice , 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 for any constant in time.
After temporal gauge fixing, we can write (1b) as either a firstorder system in , ,
or as a secondorder equation in alone,
In the special case where and are simply positive constants with (as in vacuum, with units chosen so that the speed of light is ) and , the latter equation just becomes
Taking the Fourier transform with respect to time (the socalled
frequency domain or timeharmonicapproach), this latter equation transforms into the eigenvalue problem for the
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 and , and then fixing the temporal gauge to arrive at a weak formulation in terms of alone.
2.2.1. Function spaces and regularity
Let be a bounded Lipschitz domain, and define the function spaces
We also define the following subspaces, with boundary conditions imposed:
Here, denotes the outer unit normal to , and restrictions to are interpreted in the trace sense.
Let be a curve in and be a curve in . It follows that is a curve in , that is a curve in , and that (1a) and (2a) hold strongly in . We also assume that both and are , symmetric, and uniformly elliptic. In particular, this implies that and are both curves in . Henceforth, we restrict our attention to such that is in fact a curve in .
Finally, let the current density be a given curve in and the charge density be a given curve in , satisfying the charge conservation condition .
2.2.2. The Lagrangian and Euler–Lagrange equations
For as above, define the Lagrangian
The Euler–Lagrange equations are
(3a)  
(3b) 
These Euler–Lagrange equations imply that solutions have additional regularity properties. Since is in , we have that is in . Likewise, since is in , we have that is in . Hence, solutions to this weak problem are in fact strong solutions of Maxwell’s equations.
Remark .
When and are constant in time, the electric and magnetic fields have precisely the same regularity assumed by Monk [27, eqs. (7)–(8)], namely: is in and in , while is in and in .
As in Section 2.1, this formulation is symmetric with respect to gauge transformations , where is now an arbitrary curve in . Fixing the temporal gauge , the Lagrangian becomes
and the Euler–Lagrange equations are just (3a). This again implies that is in , so (1b) holds strongly. By the same argument as in Section 2.1, this automatically preserves the chargeconservation constraint.
Remark .
Preservation of the chargeconservation constraint may also be seen as a consequence of the remaining gauge symmetry , mentioned in Section 2.1.2, where 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 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 finitedimensional subspace , resulting in a finitedimensional system of ODEs. That is, we seek a curve such that
(4) 
where , , , and . The discrete versions of (1a) and (2a),
follow immediately. In fact, both hold strongly in , by the same argument as in Section 2.2.1, since and . On the other hand, we cannot conclude that is in , nor that is in , since (4) only holds for test functions in and not all of .
Consequently, the chargeconservation constraint (2b) is only preserved in the following, much weaker sense. Let be a finitedimensional subspace such that . Then, for all , taking in (4) and applying gives
Hence, if the initial conditions satisfy , for all , then this is preserved by the flow of (4).
In particular, suppose now that is polyhedral, and that is a triangulation of by simplices (i.e., tetrahedra) . We may take to be the space of continuous degree piecewise polynomials on vanishing on , corresponding to standard Lagrange finite elements. For , we may take either degree Nédélec edge elements of the first kind [29] or degree Nédélec edge elements of the second kind [30]
with vanishing degrees of freedom on
. These are spaces of piecewisepolynomial vector fields in with tangential (but not necessarily normal) continuity between neighboring simplices. These choices ensure that , so the weak chargeconservation argument above holds.Note, however, that only says that holds in an “averaged” sense, since (unlike in the infinitedimensional case) nonzero cannot be taken to have arbitrarily small support. We cannot even conclude that the constraint holds in the sense that , since the indicator function 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 ; see also Christiansen [14].) This motivates our proposed hybrid approach, based on domain decomposition, for which piecewiseconstants are admissible test functions.
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 into a collection of problems on , 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 and correspond to the traces of and , respectively, and we show that the latter satisfies an appropriate version of the chargeconservation 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:
Brezzi and Fortin [12, Proposition III.1.1] show that
That is, is the subspace of where internal continuity and external boundary conditions are enforced by Lagrange multipliers . Likewise, [12, Proposition III.1.2] shows that
Using a similar argument, we now prove the corresponding result for the spaces.
Proposition .
.
Proof.
If , then for any , we have
so the forward inclusion holds. To get the reverse inclusion , suppose that satisfies the condition above, and let . Then, integrating by parts, we have
where the last line uses the triangle and Cauchy–Schwarz inequalities. It follows that , so . This implies that for all . Hence, in the trace sense, which completes the proof. ∎
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 and , and introduce the Lagrange multipliers and . 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 is and that is , such that is . Furthermore, suppose that and are both . Define the Lagrangian
The Euler–Lagrange equations are then
(5a)  
(5b)  
(5c)  
(5d) 
where (5a) and (5b) hold for all . We now relate this to the classical variational form of Maxwell’s equations, stated in (3).
Proposition .
Proof.
Suppose is a solution to (5). By Section 3.1, (5c) implies , so taking and summing (5a) over , the integrals over cancel, yielding (3a). As previously stated, (3a) implies , so substituting this into (5a) gives
so . Similarly, (5d) implies , so taking and summing (5b) over yields (3b). This implies , and substituting into (5b) gives .
Remark .
Note that, in addition to (5b) implying that , we also see by taking that satisfies the conservation law , for all .
3.3. Temporal gauge fixing and the chargeconservation constraint
As in Section 2.2, if is a solution to (5), then so is for any curve . Therefore, we perform temporal gauge fixing by taking . This yields the gaugefixed Lagrangian
whose Euler–Lagrange equations are simply (5a) and (5c). Of course, (5d) is satisfied trivially, since . The next result shows that the chargeconservation constraint (5b) is automatically preserved, for an appropriatelydefined .
Proposition .
Proof.
Remark .
As in Section 3.2, taking implies . Furthermore, if the initial conditions also satisfy , then we have for all time, since . Finally, if , and if the initial conditions for equal those for , then we recover .
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
we seek and such that
(6a)  
(6b) 
where is the inner product. Defining the map , , we see that (6) is equivalent to evolving by
and subsequently solving for satisfying (6a). Since by Section 3.1, it follows that solves the nondomaindecomposed problem (3a).
4. Hybrid semidiscretization
We now perform Galerkin semidiscretization of the domaindecomposed 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 and their test functions are both restricted to a subspace of . We then show that a suitablydefined satisfies the chargeconservation constraint in a strong sense, as opposed to the much weaker sense in which 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 2D Maxwell equations.
4.1. Semidiscretization of the variational problem
For each , let be a finitedimensional subspace, so , and let . We seek and such that
(7a)  
(7b) 
where (7a) holds for all . These are the semidiscretized versions of (5a) and (5c).
Remark .
Since (7b) only holds for test functions in , but not necessarily an arbitrary test function in , in general a solution will have . Hence, this method is generally not curlconforming and is distinct from the conforming methods discussed in Section 2.3.
In terms of the bilinear forms and , this method may be written as
(8a)  
(8b) 
Defining the operator , , we see that (8) is equivalent to evolving by
(9) 
and subsequently solving for satisfying (8a).
Since is finitedimensional, we may apply Banach’s closed range theorem to deduce that is in the range of , so a solution exists, although generally not uniquely. A natural choice is to find the solution minimizing , which in a weak sense minimizes the distance between and . This existencewithoutuniqueness is typical of hybrid methods, and one may formally resolve this by replacing by the quotient space (cf. Brezzi and Fortin [12, IV.1.3]). In practice, the evolution on specified by (9) is the essence of the method, and solving for may be seen as an optional postprocessing step.
4.2. Preservation of the chargeconservation constraint
In order to discuss the chargeconservation constraint, we first suppose that are such that and for all . We consider whether the following discretization of (5b) is preserved,
(10) 
for suitably defined.
Theorem 4.1.
Proof.
The proof is essentially similar to that of Section 3.3. Given , taking in (7a) and integrating by parts,
so if (10) holds at the initial time, then it holds for all time. The conclusion that follows by taking , and implies that if holds at the initial time, then it holds for all time. ∎