# Differential Forms in Lattice Field Theories: An Overview

Authors: F. L. Teixeira
Journal: ISRN Mathematical Physics
 (2013)

---

## Abstract

We provide an overview on the application of the exterior calculus of differential forms to the ab initio formulation of lattice field theories, with a focus on irregular or “random” lattices.


---

## Body

## 1. Introduction

The need to formulate field theories on a lattice (mesh, grid) arises from two main reasons, which may occur simultaneously or not. First, the lattice provides a natural “regularization” of divergences in lieu of renormalization techniques [1]. Such regularization does not need to be viewed as an ad hoc step, but instead as a natural consequence of assuming the field theory to be, at some fundamental level, an effective (“low’’-energy) description [2]. Second, the lattice provides a direct route to compute, in a nonperturbative fashion, quantities of interest by numerical simulations. Nontrivial domains and complex boundary conditions can then be easily treated as well [3–6]. For these, the use of irregular (“random”) lattices are often of interest to gain geometrical flexibility. Irregular lattices are also of interest as a means to provide a potentially faster convergence to the continuum limit, near-isotropic lattice dispersion properties, and better “conservation” of some (e.g., long-range translational and rotational) symmetries [7, 8]. In some cases, irregular lattices are useful for universality tests as well [9, 10].Lattice theories are typically developed by taking the counterpart continuum theory as starting point and then applying discretization techniques whereby derivatives are approximated by finite differences or some constraints are enforced on the functional space of admissible solutions to be spanned by a finite set of “basis” functions (e.g., “Galerkin methods” such as spectral elements and finite elements). These discretization strategies have proved very useful in many settings; however, they often produce difficulties in the case of irregular (“random”) lattices. Among such difficulties are (i) numerical instabilities in marching-on-time algorithms (regardless of the time integration method used), (ii) convergence problems in algorithms relying on iterative linear solvers, and (iii) spurious (“ghost”) modes and/or extraneous degrees of freedom. These problems often (but not always) appear associated with highly skewed or obtuse lattice elements, or at the boundary between heterogeneous (hybrid) lattices subcomponent, comprising overlapped domains or “mesh-stitching” interfaces, for example. Clearly, such difficulties put a constraint on the geometric flexibility that irregular lattices are intended for, and may require stringent (and computationally demanding) mesh quality controls. These difficulties also impact the ability to utilize “mesh refinement” strategies based on a priori error estimates. The reasons behind these difficulties can be traced to an inconsistent rendering of the differential calculus and degrees of freedom on the lattice. A rough classification of those inconsistencies is provided in AppendixB.The objective of this paper is to present an overview on the application of the exterior calculus of differential forms to the ab initio formulation of field theories on irregular lattices [11–29]. In the exterior calculus framework, the lattice is treated as a cell complex (in the parlance of algebraic topology [30]) instead of simply a collection of discrete points, and dynamic fields are represented by means of discrete differential forms (cochains) of various degrees [28, 31, 32]. This prescription provides a basis for developing a consistent “discrete calculus” on irregular lattices, and discrete analogues to partial differential equations that better adhere to the underlying physics.This topic intersects many disparate application areas. For concreteness, we employ classical electrodynamics as the standard example here; however, whenever particularly relevant to do so, we provide brief pointers to other field theories as well. Although some familiarity with the exterior calculus of differential forms is assumed [18, 19, 33–39], the discussion is mostly kept at a tutorial level. Finally, we stress that this is a review paper and no claim of originality is intended.


## 2. Premetric Lattice Equations

Let us denote the space of differentialp-forms on a smooth connected manifold Ω as Λp(Ω). From a geometric perspective, a differential p-form αp∈Λp(Ω) can be viewed as an oriented p-dimensional density, or an object naturally associated with p-dimensional domains of integration Up such that the lattice contraction (“pairing”) below:
(1)〈Up,αp〉≐∫Upαp
gives a real number (in our context) for each choice ofUp [22]. On a lattice 𝒦,Up is restricted to be a union of elements from the finite set of p-dimensional Np-oriented lattice elements, which we denote by Γp(𝒦)={σp,i,i=1,…,Np}. These are collectively called “p-chains.” In four dimensions for example, they correspond to the possible unions of elements from the set of vertices (nodes) σ0, edges (links) σ1, facets (plaquettes) σ2, volume cells (voxels) σ3, and hypervolume cells σ4, for p=1,…,4, respectively. In the discrete setting, the degrees of freedom are reduced to the set of pairings (1) on each one of the lattice elements.On the lattice, the pairing above can be understood as a mapℛp:Λp(Ω)→Γp(𝒦) such that
(2)ℛp(αp)=〈σp,i,αp〉≐∫σp,iαp
defines its action on the basis of p-chains. Note that we use Γp(𝒦) to denote the space dual to Γp(𝒦), that is, the space p-cochains. The latter can be viewed as the space of “discrete differential forms.” Because of this, and with some abuse of language, we use the terminology “differential forms” and “cochains” interchangeably to denote the same objects in what follows. The map ℛp is called the de Rham map [22].The basic differential operator of exterior calculus is the exterior derivatived, applicable to any number of dimensions. The discretization of d on a general irregular lattice can be effected by a straightforward application of the generalized Stokes' theorem [22]:
(3)∫σp+1dαp=∫∂σp+1αp
with p=0,…,3 in n=4. In the above, ∂ is the boundary operator, which simply maps a p-dimensional lattice element to the set of (p-1)-dimensional lattice elements that comprise its boundary, preserving orientation. This theorem sets ∂ as the formal adjoint of d in terms of the pairing given in (1), that is, 〈σp+1,dαp〉=〈∂σp+1,αp〉. Computationally, the boundary operator can be implemented by means of incidence matrices [22, 29, 40] such that
(4)∂σp+1,i=∑jCijpσp,j,
where the indices i and j run over all (p+1)- and p-dimensional lattice elements, respectively. The incidence matrix entries are such that Cijp∈{-1,0,1} for all p, with sign determined by the relative orientation of lattice elements i and j. The restriction to this set of integer values reflects the “metric-free” nature of the exterior derivative: only information about element connectivity, that is, the combinatorial aspects of the lattice, is involved here. It turns out that the metric is fully encoded by Hodge star operators, the discretization of which will be discussed further down below.Using (3) and (4), one can write
(5)∫σp+1,idαp=∑jCijp∫σp,jαp
for all i, so that the derivative operation is replaced by a proper sum over j. On the lattice, the nilpotency of the operators ∂∘∂=d∘d=0 [41] is recovered by the constraint [22]:
(6)∑kCikp+1Ckjp=0
for all i and j.


## 3. Example: Lattice Electrodynamics

We write Maxwell's equations in a four-dimensional Lorentzian manifoldΩ as [34]:
(7)dF=0,dG=*𝒥,
where d is the four-dimensional exterior derivative, F and G are the so-called Faraday and Maxwell 2-forms, respectively, and *𝒥 is the charge-current density 3-form. The Hodge star operator * is an isomorphism that maps p-forms to (4-p)-forms, and more generally p forms to (n-p) forms in a n-dimensional manifold, and, as mentioned before, depends on the metric of Ω [22, 23, 34, 35, 42–45]. The above equations are complemented by the relation G=*F, which indicates that F and G are “Hodge duals” of each other.


### 3.1. Primal and Dual Lattices

SinceF and G are 2-forms, they should be discretized as 2-cochains residing on plaquettes (2-chains) of the 4-dimensional lattice; however, it is important to recognize that these two forms are of different types: F is an “ordinary” (or “nontwisted”) differential form, whereas G (as well as *𝒥) is a “twisted” (or “odd”) differential form [46]. The basic difference here has to do with orientation: ordinary forms have internal orientation whereas twisted forms have external orientation [20, 22, 46–48]. These two types of orientations exhibit different symmetries under reflection, a distinction akin to that between proper (or polar) tensors and pseudo (or axial) tensors. Only twisted forms admit integration in nonorientable manifolds. These two types of forms are associated with two distinct “cell complexes” (lattices), each one inheriting the corresponding orientation: the ordinary form F is associated with the set of plaquettes Γ2 on the “ordinary cell complex” 𝒦, thus belonging to Γ2(𝒦), while the twisted forms G and *𝒥 are associated with the set of plaquettes Γ~2 on the “twisted cell complex” 𝒦~ [22, 27, 48, 49], thus belonging to Γ2(𝒦~). Consequently, we also have two sets of incidence matrices Cijp and C~ijp, one for each lattice. It is convenient to denote 𝒦 as the “primal lattice” and 𝒦~ as the “dual lattice” [22].As detailed further below, these two lattices become intertwined by the Hodge dualityF=*G. The need for dual lattices can also be motivated from a purely combinatorial standpoint (as a means to preserve key topological properties from the continuum theory) [24] or from a strictly computational standpoint (e.g., to provide higher-order convergence to the continuum) [50–52].


### 3.2.3+1 Theory

At this point, it is suitable to degeometrize time and treat it simply as a parameter. This corresponds to the majority of low-energy applications involving Maxwell's equations, in which one is interested in predicting the field evolution along different spatial slices for a given set of initial and boundary conditions. In this case, we still use the symbols𝒦 and 𝒦~ for the primal and dual lattices, but they now refer to three-dimensional spatial lattices. Similarly, Ω now refers to a three-dimensional Euclidean manifold. In such a 3+1 setting, one can decompose F and G as
(8)F=E∧dt+B,G=D-H∧dt,
and the source density as
(9)*𝒥=-J∧dt+ρ,
where ∧ is the wedge product, E and H are the electric intensity and magnetic intensity 1-forms on Γ1 and Γ~1, respectively, D and B are the electric flux and magnetic flux 2-forms on Γ~2 and Γ2, respectively, J is the electric current density 2-form on Γ~2, and ρ is the electric charge density 3-form on Γ~3 (corresponding assignments for the 2+1 and 1+1 cases are provided in [32]). As a result, Maxwell's equations reduce to
(10)dE=-∂tB,(11)dH=∂tD+J,
representing Faraday's and Ampere's laws, respectively. Here, d stands for the 3-dimensional spatial exterior derivative. Note that both (10) and (11) are metric-free. They are supplemented by Hodge star relations given by
(12)D=⋆ϵE,H=⋆μ-1B
now involving two Hodge star maps in three-dimensional space: ⋆ϵ:Λ1(Ω)→Λ2(Ω) and ⋆μ-1:Λ2(Ω)→Λ1(Ω). On the lattice, we have the corresponding discrete counterparts: [⋆ϵ]:Γ1(𝒦)→Γ2(𝒦~) and [⋆μ-1]:Γ2(𝒦)→Γ1(𝒦~). The subscripts ϵ and μ in ⋆ϵ and  ⋆μ-1 serve to indicate that these operators also incorporate macroscopic constitutive material properties through the local permittivity and permeability values [53] (we assume dispersionless media for simplicity). In Riemannian manifolds (and in particular, Euclidean space) and reciprocal media, these two Hodge star operators are symmetric and positive-definite [54].In what follows, we employ the following short-hand notation for cochains:〈σ1,i,E〉=Ei, 〈σ~1,i,H〉=Hi, 〈σ~2,i,D〉=Di, 〈σ2,i,B〉=Bi, 〈σ~2,i,J〉=Ji, and 〈σ~3,i,ρ〉=ρi, where the indices run over the respective basis of p-chains in either 𝒦 or 𝒦~, p=1,2,3. With the exception of Appendix A, we restrict ourselves to the 3+1 setting throughout the remainder of this paper.


## 3.1. Primal and Dual Lattices

SinceF and G are 2-forms, they should be discretized as 2-cochains residing on plaquettes (2-chains) of the 4-dimensional lattice; however, it is important to recognize that these two forms are of different types: F is an “ordinary” (or “nontwisted”) differential form, whereas G (as well as *𝒥) is a “twisted” (or “odd”) differential form [46]. The basic difference here has to do with orientation: ordinary forms have internal orientation whereas twisted forms have external orientation [20, 22, 46–48]. These two types of orientations exhibit different symmetries under reflection, a distinction akin to that between proper (or polar) tensors and pseudo (or axial) tensors. Only twisted forms admit integration in nonorientable manifolds. These two types of forms are associated with two distinct “cell complexes” (lattices), each one inheriting the corresponding orientation: the ordinary form F is associated with the set of plaquettes Γ2 on the “ordinary cell complex” 𝒦, thus belonging to Γ2(𝒦), while the twisted forms G and *𝒥 are associated with the set of plaquettes Γ~2 on the “twisted cell complex” 𝒦~ [22, 27, 48, 49], thus belonging to Γ2(𝒦~). Consequently, we also have two sets of incidence matrices Cijp and C~ijp, one for each lattice. It is convenient to denote 𝒦 as the “primal lattice” and 𝒦~ as the “dual lattice” [22].As detailed further below, these two lattices become intertwined by the Hodge dualityF=*G. The need for dual lattices can also be motivated from a purely combinatorial standpoint (as a means to preserve key topological properties from the continuum theory) [24] or from a strictly computational standpoint (e.g., to provide higher-order convergence to the continuum) [50–52].


## 3.2.3+1 Theory

At this point, it is suitable to degeometrize time and treat it simply as a parameter. This corresponds to the majority of low-energy applications involving Maxwell's equations, in which one is interested in predicting the field evolution along different spatial slices for a given set of initial and boundary conditions. In this case, we still use the symbols𝒦 and 𝒦~ for the primal and dual lattices, but they now refer to three-dimensional spatial lattices. Similarly, Ω now refers to a three-dimensional Euclidean manifold. In such a 3+1 setting, one can decompose F and G as
(8)F=E∧dt+B,G=D-H∧dt,
and the source density as
(9)*𝒥=-J∧dt+ρ,
where ∧ is the wedge product, E and H are the electric intensity and magnetic intensity 1-forms on Γ1 and Γ~1, respectively, D and B are the electric flux and magnetic flux 2-forms on Γ~2 and Γ2, respectively, J is the electric current density 2-form on Γ~2, and ρ is the electric charge density 3-form on Γ~3 (corresponding assignments for the 2+1 and 1+1 cases are provided in [32]). As a result, Maxwell's equations reduce to
(10)dE=-∂tB,(11)dH=∂tD+J,
representing Faraday's and Ampere's laws, respectively. Here, d stands for the 3-dimensional spatial exterior derivative. Note that both (10) and (11) are metric-free. They are supplemented by Hodge star relations given by
(12)D=⋆ϵE,H=⋆μ-1B
now involving two Hodge star maps in three-dimensional space: ⋆ϵ:Λ1(Ω)→Λ2(Ω) and ⋆μ-1:Λ2(Ω)→Λ1(Ω). On the lattice, we have the corresponding discrete counterparts: [⋆ϵ]:Γ1(𝒦)→Γ2(𝒦~) and [⋆μ-1]:Γ2(𝒦)→Γ1(𝒦~). The subscripts ϵ and μ in ⋆ϵ and  ⋆μ-1 serve to indicate that these operators also incorporate macroscopic constitutive material properties through the local permittivity and permeability values [53] (we assume dispersionless media for simplicity). In Riemannian manifolds (and in particular, Euclidean space) and reciprocal media, these two Hodge star operators are symmetric and positive-definite [54].In what follows, we employ the following short-hand notation for cochains:〈σ1,i,E〉=Ei, 〈σ~1,i,H〉=Hi, 〈σ~2,i,D〉=Di, 〈σ2,i,B〉=Bi, 〈σ~2,i,J〉=Ji, and 〈σ~3,i,ρ〉=ρi, where the indices run over the respective basis of p-chains in either 𝒦 or 𝒦~, p=1,2,3. With the exception of Appendix A, we restrict ourselves to the 3+1 setting throughout the remainder of this paper.


## 4. Casting the Metric on a Lattice

### 4.1. Whitney Forms

The Whitney map𝒲:Γp(𝒦)→Λp(Ω) is the right-inverse of the de Rham map (2), that is, ℛ∘𝒲=ℐ, where ℐ is the identity operator. In simplicial lattices, this morphism can be constructed using the so-called Whitney forms [15, 22, 36, 43, 55–61] which are basic interpolants from cochains to differential forms [33] (other interpolants are also possible [62, 63]). By definition, all cell elements of a simplicial lattice are simplices, that is, cells whose boundaries are the union of a minimal number of lower-dimensional cells. In other words, 0-simplices are nodes, 1-simplices are links, 2-simplices are triangles, 3-simplices are tetrahedra, and so on. Note that if the primal lattice is simplicial, the dual lattice is not [31]. For a p-simplex σp,i, the (lowest-order) Whitney form is given by
(13)ωp[σp,i]≐p!∑j=0p(-1)iλi,jdλi,0∧dλi,1⋯dλi,j-1∧dλi,j+1⋯dλi,p,
where λi,j, j=0,…,p, are the barycentric coordinates associated with σp,i. In the case of a 0-simplex (node), (13) reduces to ω0[σ0,i]=λi.From its definition, it is clear that Whitney forms have compact support. Among its important structural properties are(14)〈σp,i,ωp[σp,j]〉=∫σp,iωp[σp,j]=δij,
where δij is the Kronecker delta, which is simply a restatement of ℛ∘𝒲=ℐ, and
(15)ωp[∂Tσp-1,i]=d(ωp-1[σp-1,i]),
where ∂T is the coboundary operator [56], consistent with the generalized Stokes' theorem. Further structural properties are provided in [57, 58]. Higher-order version of Whitney forms also exist [59, 60]. The key result 𝒲∘ℛ→ℐ holds in the limit of zero lattice spacing. This is discussed, together with other related convergence results in various contexts, in [15, 33, 64–68].Using the short-handωp[σp,i]=ωip, we can write the following expansions for E and B in a irregular simplicial lattice, in terms of its cochain representations:
(16)E=∑iEiωi1,B=∑iBiωi2,
where the sums run over all primal lattice edges and faces, respectively.One could argue that Whitney forms are continuum objects that should have no fundamental place on a truly discrete theory. In our view, this is only partially true. In many applications (see, e.g., the discussion on space-charge effects below), it is less natural to consider the lattice as endowed with some a priori discrete metric structure than it is to consider it instead as embedded in an underlying continuum (say, Euclidean) manifold with metric and hence inheriting all metric properties from it. In the latter case, Whitney forms provide the standard route to incorporate metric information into the discrete Hodge star operators, as described next.


### 4.2. Discrete Hodge Star Operator

In a source-free media, we can write the Hamiltonian as(17)ℋ=12∫Ω‍(E∧D+H∧B)=∫Ω(E∧⋆ϵE+⋆μ-1B∧B).
Using (16), the lattice Hamiltonian assumes the expected quadratic form:
(18)ℋ=∑i∑jEi[⋆ϵ]ijEj+∑i∑jBi[⋆μ-1]ijBj,
where we immediately identify the symmetric positive definite matrices:
(19)[⋆ϵ]ij=∫Ωωi1∧⋆ϵωj1,[⋆μ-1]ij=∫Ω(⋆μ-1ωi2)∧ωj2
as the discrete realization of the Hodge star operator(s) on a simplicial lattice [23, 69] so that
(20)Di=∑j[⋆ϵ]ijEj,Hi=∑j[⋆μ-1]ijBj.
From the above, the Hamiltonian can be also expressed as
(21)ℋ=∑iEiDi+∑iHiBi.


### 4.3. Symplectic Structure and Dynamic Degrees of Freedom

The Hodge star matrices[⋆ϵ] and [⋆μ-1] have different sizes. The number of elements in [⋆ϵ] is equal to N1×N1, whereas the number of elements in [⋆μ-1] is equal to N2×N2. In other words, Θ(E)=Θ(D)≠Θ(B)=Θ(H), where Θ denotes the number of (discrete) degrees of freedom in the corresponding field.One important property of a Hamiltonian system is its symplectic character, associated with area preservation in phase space. The symplectic character of the Hamiltonian in principle would require a canonical pair such asE,B to have identical number of degrees of freedom. This apparent contradiction can be explained by the fact that Maxwell's equations (10) and (11) can be thought as aconstrained dynamic system (by the divergence conditions) so that, even though Θ(E)≠Θ(B), we still have Θd(E)=Θd(B), where Θd denotes the number of dynamic degrees of freedom. This is discussed further below in Section 6, in connection with the discrete Hodge decomposition on a lattice.


## 4.1. Whitney Forms

The Whitney map𝒲:Γp(𝒦)→Λp(Ω) is the right-inverse of the de Rham map (2), that is, ℛ∘𝒲=ℐ, where ℐ is the identity operator. In simplicial lattices, this morphism can be constructed using the so-called Whitney forms [15, 22, 36, 43, 55–61] which are basic interpolants from cochains to differential forms [33] (other interpolants are also possible [62, 63]). By definition, all cell elements of a simplicial lattice are simplices, that is, cells whose boundaries are the union of a minimal number of lower-dimensional cells. In other words, 0-simplices are nodes, 1-simplices are links, 2-simplices are triangles, 3-simplices are tetrahedra, and so on. Note that if the primal lattice is simplicial, the dual lattice is not [31]. For a p-simplex σp,i, the (lowest-order) Whitney form is given by
(13)ωp[σp,i]≐p!∑j=0p(-1)iλi,jdλi,0∧dλi,1⋯dλi,j-1∧dλi,j+1⋯dλi,p,
where λi,j, j=0,…,p, are the barycentric coordinates associated with σp,i. In the case of a 0-simplex (node), (13) reduces to ω0[σ0,i]=λi.From its definition, it is clear that Whitney forms have compact support. Among its important structural properties are(14)〈σp,i,ωp[σp,j]〉=∫σp,iωp[σp,j]=δij,
where δij is the Kronecker delta, which is simply a restatement of ℛ∘𝒲=ℐ, and
(15)ωp[∂Tσp-1,i]=d(ωp-1[σp-1,i]),
where ∂T is the coboundary operator [56], consistent with the generalized Stokes' theorem. Further structural properties are provided in [57, 58]. Higher-order version of Whitney forms also exist [59, 60]. The key result 𝒲∘ℛ→ℐ holds in the limit of zero lattice spacing. This is discussed, together with other related convergence results in various contexts, in [15, 33, 64–68].Using the short-handωp[σp,i]=ωip, we can write the following expansions for E and B in a irregular simplicial lattice, in terms of its cochain representations:
(16)E=∑iEiωi1,B=∑iBiωi2,
where the sums run over all primal lattice edges and faces, respectively.One could argue that Whitney forms are continuum objects that should have no fundamental place on a truly discrete theory. In our view, this is only partially true. In many applications (see, e.g., the discussion on space-charge effects below), it is less natural to consider the lattice as endowed with some a priori discrete metric structure than it is to consider it instead as embedded in an underlying continuum (say, Euclidean) manifold with metric and hence inheriting all metric properties from it. In the latter case, Whitney forms provide the standard route to incorporate metric information into the discrete Hodge star operators, as described next.


## 4.2. Discrete Hodge Star Operator

In a source-free media, we can write the Hamiltonian as(17)ℋ=12∫Ω‍(E∧D+H∧B)=∫Ω(E∧⋆ϵE+⋆μ-1B∧B).
Using (16), the lattice Hamiltonian assumes the expected quadratic form:
(18)ℋ=∑i∑jEi[⋆ϵ]ijEj+∑i∑jBi[⋆μ-1]ijBj,
where we immediately identify the symmetric positive definite matrices:
(19)[⋆ϵ]ij=∫Ωωi1∧⋆ϵωj1,[⋆μ-1]ij=∫Ω(⋆μ-1ωi2)∧ωj2
as the discrete realization of the Hodge star operator(s) on a simplicial lattice [23, 69] so that
(20)Di=∑j[⋆ϵ]ijEj,Hi=∑j[⋆μ-1]ijBj.
From the above, the Hamiltonian can be also expressed as
(21)ℋ=∑iEiDi+∑iHiBi.


## 4.3. Symplectic Structure and Dynamic Degrees of Freedom

The Hodge star matrices[⋆ϵ] and [⋆μ-1] have different sizes. The number of elements in [⋆ϵ] is equal to N1×N1, whereas the number of elements in [⋆μ-1] is equal to N2×N2. In other words, Θ(E)=Θ(D)≠Θ(B)=Θ(H), where Θ denotes the number of (discrete) degrees of freedom in the corresponding field.One important property of a Hamiltonian system is its symplectic character, associated with area preservation in phase space. The symplectic character of the Hamiltonian in principle would require a canonical pair such asE,B to have identical number of degrees of freedom. This apparent contradiction can be explained by the fact that Maxwell's equations (10) and (11) can be thought as aconstrained dynamic system (by the divergence conditions) so that, even though Θ(E)≠Θ(B), we still have Θd(E)=Θd(B), where Θd denotes the number of dynamic degrees of freedom. This is discussed further below in Section 6, in connection with the discrete Hodge decomposition on a lattice.


## 5. Semidiscrete Equations

### 5.1. Local and Ultralocal Lattice Coupling

By using a contraction in the form of (2) on both sides of (10) with every face σ2,j of 𝒦, and using the fact that 〈σ2,j,ωi2〉=〈σ1,j,ωi1〉=δij from (14), we get
(22)〈σ2,j,∂tB〉=∂t∑iBi〈σ2,j,ωi2〉=∂tBj,〈σ2,j,dE〉=〈∂σ2,j,E〉=∑iEi∑kCjk1〈σ1,k,ωi1〉=∑iCji1Ei
so that
(23)-∂tBi=∑jCij1Ej,
where the index i runs over all faces of the primal lattice. On the dual lattice 𝒦~, we can similarly contract both sides of (11) with every dual face σ~2,j to get
(24)∂tDi=∑jC~ij1Hj,
where now the index i runs over all faces of the dual lattice. Using (20) and the fact that, in three-dimensions C~ij1=Cji1 [22] (up to possible boundary terms ignored here), we can write the last equation in terms of primal lattice quantities as
(25)∂t∑j[⋆ϵ]ijEj=∑jCji1∑k[⋆μ-1]jkBk
or, by using the inverse Hodge star matrix [⋆ϵ]ij-1, as
(26)∂tEi=∑jΥijBj
with
(27)Υij≐∑k∑l[⋆ϵ]ik-1Clk1[⋆μ-1]lj.
The matrix [Υ] can be viewed as the discrete realization, for p=2, of the codifferential operator δ=(-1)p*-1d* that maps p-forms to (n-p)-forms [35].Since the continuum operators⋆ϵ and ⋆μ-1 are local [46] and, as seen, Whitney forms (13) have local support, it follows that the matrices [⋆ϵ] and [⋆μ-1] are sparse, indicative of an ultralocal coupling (in the terminology of [70]). In contrast, the numerical inverse [⋆ϵ]-1 used in (27) is, in general, not sparse so that the field coupling between distant elements is nonzero. The lack of sparsity is a potential bottleneck in practical simulations. However, because the coupling strength in this case decays exponentially [29, 44], we can still say (using again the terminology of [70]) that the resulting discrete operator encoded by the matrix in (27) is local. In practical terms, the exponential decay allows one to set a cutoff on the nonzero elements of [⋆ϵ], based on element magnitudes or on the sparsity pattern of the original matrix [⋆ϵ], to build a sparse approximate inverse for [⋆ϵ] and hence recover back an ultralocal representation for ⋆ϵ-1 [29, 71]. The sparsity pattern of [⋆ϵ] encodes the nearest-neighbor edge information of the mesh and, consequently, the sparsity pattern of [⋆ϵ]k likewise encodes successive “k-level” neighbors. The latter sparsity patterns can be used to build, quite efficiently, sparse approximations for [⋆ϵ]-1, as detailed in [29]. Once such sparse representations are obtained, (23) and (26) can be used in tandem to construct a marching-on-time algorithm (e.g., see Section 9.1 ahead) with a sparse structure and hence amenable for large-scale problems.


### 5.2. Barycentric Dual and Barycentric Decomposition Lattices

An alternative approach, aimed at constructing a sparse discrete Hodge star for⋆ϵ-1 directly from the dual lattice geometry is described in [27], based on earlier ideas exposed in [24, 72]. This approach is based on the fact that both primal 𝒦 and dual 𝒦~ lattices can be decomposed into a third (underlying) lattice 𝒦^ by means of a barycentric decomposition, see [24]. The dual lattice 𝒦~ in this case is called the barycentric dual lattice [27, 72] and the underlying lattice 𝒦^ is called the barycentric decomposition lattice. Importantly, 𝒦^ is simplicial and hence admits Whitney forms built on it using (13). Whitney forms on 𝒦^ can be used as building blocks to construct (dual) Whitney forms on the (nonsimplicial) 𝒦~, and from that, a sparse inverse discrete Hodge star [⋆ϵ-1] using integrals akin to (19). An explicit derivation of such dual lattice Whitney forms is provided in [73]. Furthermore, a recent comprehensive survey of this and other approaches based on dual lattices to construct discrete sparse inverse Hodge stars is provided in [74].The barycentric dual lattice has the important property below associated with Whitney forms:(28)〈σ~(n-p),i,⋆ωp[σp,j]〉=∫σ~(n-p),i⋆ωp[σp,j]=δij,
where ⋆ stands for the spatial Hodge star operator (distilled from constitutive material properties), and σ~(n-p),i is the dual element to σp,i on the barycentric dual lattice. The operator ⋆ is such that
(29)∫Ωωp∧⋆ωp=∫Ω|ω|2dv,
where |ω|2 is the two-norm of ωp, and dv is the volume element.The identity (28) plays the role of structural property (14), on the dual lattice side. We stress that identity (28) is a distinctively characteristic feature of the barycentric dual lattice not shared by other geometrical constructions for the dual lattice. In other words, compatibility with Whitney forms via (28) naturally forces one to choose the dual lattice to be the barycentric dual.From the above, one can also define a (Hodge) duality operator directly on the space of chains, that is,⋆K:Γp(𝒦)↦Γn-p(𝒦~) with ⋆K(σp,i)=σ~(n-p),i and ⋆K~:Γp(𝒦~)↦Γn-p(𝒦) with ⋆K(σ~p,i)=σ~(n-p),i, so that ⋆K⋆K~=⋆K~⋆K=1. This construction is detailed in [24].


### 5.3. Galerkin Duality

Even though we have chosen to assignE and B to the primal (simplicial) lattice, and consequently D, H, J, and ρ to the dual (nonsimplicial) lattice, the reverse is equally possible. In this case, the fields D, H become associated to a simplicial lattice and hence can be expressed in terms of Whitney forms; the expressions dual to (16) are now
(30)H=∑iHiωi1,D=∑iDiωi2
with sums running over primal edges and primal faces, respectively, and where
(31)Ei=∑j[⋆ϵ-1]ijDj,Bi=∑j[⋆μ]ijHj
with
(32)[⋆ϵ-1]ij=∫Ω(⋆ϵ-1ωi2)∧ωj2,[⋆μ]ij=∫Ωωi1∧⋆μωj1,
and the two Hodge star maps now used are such that, in the continuum, ⋆ϵ-1:Λ2(Ω)→Λ1(Ω) and ⋆μ:Λ1(Ω)→Λ2(Ω), and, on the lattice, [⋆ϵ-1]:Γ2(𝒦)→Γ1(𝒦~) and [⋆μ]:Γ1(𝒦)→Γ2(𝒦~). This alternate choice entails a duality between these two formulations, dubbed “Galerkin duality.” This is explored in more detail in [44].


## 5.1. Local and Ultralocal Lattice Coupling

By using a contraction in the form of (2) on both sides of (10) with every face σ2,j of 𝒦, and using the fact that 〈σ2,j,ωi2〉=〈σ1,j,ωi1〉=δij from (14), we get
(22)〈σ2,j,∂tB〉=∂t∑iBi〈σ2,j,ωi2〉=∂tBj,〈σ2,j,dE〉=〈∂σ2,j,E〉=∑iEi∑kCjk1〈σ1,k,ωi1〉=∑iCji1Ei
so that
(23)-∂tBi=∑jCij1Ej,
where the index i runs over all faces of the primal lattice. On the dual lattice 𝒦~, we can similarly contract both sides of (11) with every dual face σ~2,j to get
(24)∂tDi=∑jC~ij1Hj,
where now the index i runs over all faces of the dual lattice. Using (20) and the fact that, in three-dimensions C~ij1=Cji1 [22] (up to possible boundary terms ignored here), we can write the last equation in terms of primal lattice quantities as
(25)∂t∑j[⋆ϵ]ijEj=∑jCji1∑k[⋆μ-1]jkBk
or, by using the inverse Hodge star matrix [⋆ϵ]ij-1, as
(26)∂tEi=∑jΥijBj
with
(27)Υij≐∑k∑l[⋆ϵ]ik-1Clk1[⋆μ-1]lj.
The matrix [Υ] can be viewed as the discrete realization, for p=2, of the codifferential operator δ=(-1)p*-1d* that maps p-forms to (n-p)-forms [35].Since the continuum operators⋆ϵ and ⋆μ-1 are local [46] and, as seen, Whitney forms (13) have local support, it follows that the matrices [⋆ϵ] and [⋆μ-1] are sparse, indicative of an ultralocal coupling (in the terminology of [70]). In contrast, the numerical inverse [⋆ϵ]-1 used in (27) is, in general, not sparse so that the field coupling between distant elements is nonzero. The lack of sparsity is a potential bottleneck in practical simulations. However, because the coupling strength in this case decays exponentially [29, 44], we can still say (using again the terminology of [70]) that the resulting discrete operator encoded by the matrix in (27) is local. In practical terms, the exponential decay allows one to set a cutoff on the nonzero elements of [⋆ϵ], based on element magnitudes or on the sparsity pattern of the original matrix [⋆ϵ], to build a sparse approximate inverse for [⋆ϵ] and hence recover back an ultralocal representation for ⋆ϵ-1 [29, 71]. The sparsity pattern of [⋆ϵ] encodes the nearest-neighbor edge information of the mesh and, consequently, the sparsity pattern of [⋆ϵ]k likewise encodes successive “k-level” neighbors. The latter sparsity patterns can be used to build, quite efficiently, sparse approximations for [⋆ϵ]-1, as detailed in [29]. Once such sparse representations are obtained, (23) and (26) can be used in tandem to construct a marching-on-time algorithm (e.g., see Section 9.1 ahead) with a sparse structure and hence amenable for large-scale problems.


## 5.2. Barycentric Dual and Barycentric Decomposition Lattices

An alternative approach, aimed at constructing a sparse discrete Hodge star for⋆ϵ-1 directly from the dual lattice geometry is described in [27], based on earlier ideas exposed in [24, 72]. This approach is based on the fact that both primal 𝒦 and dual 𝒦~ lattices can be decomposed into a third (underlying) lattice 𝒦^ by means of a barycentric decomposition, see [24]. The dual lattice 𝒦~ in this case is called the barycentric dual lattice [27, 72] and the underlying lattice 𝒦^ is called the barycentric decomposition lattice. Importantly, 𝒦^ is simplicial and hence admits Whitney forms built on it using (13). Whitney forms on 𝒦^ can be used as building blocks to construct (dual) Whitney forms on the (nonsimplicial) 𝒦~, and from that, a sparse inverse discrete Hodge star [⋆ϵ-1] using integrals akin to (19). An explicit derivation of such dual lattice Whitney forms is provided in [73]. Furthermore, a recent comprehensive survey of this and other approaches based on dual lattices to construct discrete sparse inverse Hodge stars is provided in [74].The barycentric dual lattice has the important property below associated with Whitney forms:(28)〈σ~(n-p),i,⋆ωp[σp,j]〉=∫σ~(n-p),i⋆ωp[σp,j]=δij,
where ⋆ stands for the spatial Hodge star operator (distilled from constitutive material properties), and σ~(n-p),i is the dual element to σp,i on the barycentric dual lattice. The operator ⋆ is such that
(29)∫Ωωp∧⋆ωp=∫Ω|ω|2dv,
where |ω|2 is the two-norm of ωp, and dv is the volume element.The identity (28) plays the role of structural property (14), on the dual lattice side. We stress that identity (28) is a distinctively characteristic feature of the barycentric dual lattice not shared by other geometrical constructions for the dual lattice. In other words, compatibility with Whitney forms via (28) naturally forces one to choose the dual lattice to be the barycentric dual.From the above, one can also define a (Hodge) duality operator directly on the space of chains, that is,⋆K:Γp(𝒦)↦Γn-p(𝒦~) with ⋆K(σp,i)=σ~(n-p),i and ⋆K~:Γp(𝒦~)↦Γn-p(𝒦) with ⋆K(σ~p,i)=σ~(n-p),i, so that ⋆K⋆K~=⋆K~⋆K=1. This construction is detailed in [24].


## 5.3. Galerkin Duality

Even though we have chosen to assignE and B to the primal (simplicial) lattice, and consequently D, H, J, and ρ to the dual (nonsimplicial) lattice, the reverse is equally possible. In this case, the fields D, H become associated to a simplicial lattice and hence can be expressed in terms of Whitney forms; the expressions dual to (16) are now
(30)H=∑iHiωi1,D=∑iDiωi2
with sums running over primal edges and primal faces, respectively, and where
(31)Ei=∑j[⋆ϵ-1]ijDj,Bi=∑j[⋆μ]ijHj
with
(32)[⋆ϵ-1]ij=∫Ω(⋆ϵ-1ωi2)∧ωj2,[⋆μ]ij=∫Ωωi1∧⋆μωj1,
and the two Hodge star maps now used are such that, in the continuum, ⋆ϵ-1:Λ2(Ω)→Λ1(Ω) and ⋆μ:Λ1(Ω)→Λ2(Ω), and, on the lattice, [⋆ϵ-1]:Γ2(𝒦)→Γ1(𝒦~) and [⋆μ]:Γ1(𝒦)→Γ2(𝒦~). This alternate choice entails a duality between these two formulations, dubbed “Galerkin duality.” This is explored in more detail in [44].


## 6. Discrete Hodge Decomposition and Euler's Formula

For anyp-form αp, we can write
(33)αp=dζp-1+δβp+1+χp,
where χp is a harmonic form [31]. This Hodge decomposition is unique. In the particular case of the 1-form E, we have
(34)E=dϕ+δA+χ,
where ϕ is a 0-form and A is a 2-form, with dϕ representing the static field, δA the dynamic field, and χ the harmonic field component (if any). In a contractible domain, χ is identically zero and the Hodge decomposition simplifies to
(35)E=dϕ+δA,
more usually known as Helmholtz decomposition in three dimensions.In the discrete setting, the degrees of freedom ofϕ are associated to the nodes of the primal lattice. Likewise, the degrees of freedom of A are associated to the facets of the primal lattice. Consequently, we have from (35) that
(36)Θd(E)=NEh-NVh=(NE-NEb)-(NV-NVb)=NE-NV,
where NV is the number of primal nodes, NE the number of primal edges, and NF the number of primal facets, with superscript b standing for boundary (fixed) elements and h for interior (free) elements.On the other hand, once we identify the lattice as a network of (in general) polyhedra, we can apply Euler's polyhedron formula on the primal lattice to obtain [44]
(37)NV-NE=1-NF+NP,
where NP represents the number of volume cells comprising the primal lattice. A similar Euler's polyhedron formula applies to the (closed, two-dimensional) boundary of the primal lattice:
(38)NVb-NEb=2-NFb.Combining (37) and (38), we have
(39)(NE-NEb)-(NV-NVb)=(NF-NFb)-(NP-1).
From the Hodge decomposition (35), we see that Θd(E) is
(40)Θd(E)=NEin-NVin=(NE-NEb)-(NV-NVb).Note that the divergence free conditiondB=0 produces one constraint on the 2-form B for each volume element. This constraint also spans the whole lattice boundary. The total number of the constrains for B is, therefore, (NP-1). Consequently, we have
(41)Θd(B)=NFin-(NP-1)=(NF-NFb)-(NP-1)
so that
(42)Θd(B)=Θd(E).
This discussion can be generalized to lattices on noncontractible domains with any number of holes (genus), where the identity Θd(B)=Θd(E) is also satisfied [31]. Moreover, from Hodge star isomorphism, we have Θd(D)=Θd(E) and Θd(H)=Θd(B).In general, we can trace a direct correspondence between quantities in the Euler polyhedron formula to the quantities in the Hodge decomposition formula. For example, each term in the two-dimensional Euler's formulaNE=NV+(NF-1)+g is associated to a corresponding term in E=dϕ+δA+χ, that is, the number of edges NE corresponds to the dimension of the space of lattice 1-forms E, which is the sum of the number of nodes NV (dimension of the space of discrete 0-forms ϕ), the number of faces (NF-1) (dimension of the space of discrete 2-forms A), and the number of holes g (dimension of the space of harmonic forms χ). A similar correspondence can be traced on a three-dimensional lattice [31]. This correspondence provides a physical picture to Euler's formula and a geometric interpretation to the Hodge decomposition.


## 7. Absorbing Boundary Conditions

In many wave scattering simulations, the presence of long-range interactions with slow (algebraic) decay, together with practical limitations in computer memory resources, implies that open-space problems necessitate the use of special techniques to suppress finite-volume effects and emulate, for example, the Sommerfeld radiation condition at infinity. Perfectly matched layers (PML) are absorbing boundary conditions commonly used for this purpose [75–78]. In the continuum limit, the PML provides a reflectionless absorption of outgoing waves, in such a way that when the PML is used to truncate a computational lattice, finite-volume effects such as spurious reflections from the outer boundary are exponentially suppressed. When first introduced in the literature [75], the PML relied upon the use of matched artificial electric and magnetic conductivities in Maxwell's equations and of a splitting of each vector field component into two subcomponents. Because of this, the resulting fields inside the PML layer are rendered “non-Maxwellian.” The PML concept was later shown to be equivalent in the Fourier domain (∂t→-iω) to a complex coordinate stretching of the coordinate space (or an analytic continuation to a complex-valued coordinate space) [76–78] and, as such, applicable to any linear wave phenomena.Inside the PML, the (local) spatial coordinateζ along the outward normal direction to each lattice boundary point is complexified as
(43)ζ→ζ~=∫0ζsζ(ζ′)dζ′,
where sζ is the so-called complex stretching variable written as sζ(ζ,ω)=aζ(ζ)+iΩζ(ζ)/ω with aζ≥1 and Ωζ≥0 (profile functions). The first inequality ensures that evanescent waves will have a faster exponential decay in the PML region, and the second inequality ensures that propagating waves will decay exponentially along ζ inside the PML. As opposed to some other lattice truncation techniques, the PML preserves the locality of the underlying differential operators and hence retains the sparsity of the formulation.For Maxwell's equations, the PML can also be affected by means of artificial material tensors (Maxwellian PML) [79]. In three dimensions, the Maxwellian PML can be represented as a media with anisotropic permittivity and permeability tensors exhibiting stratification along the normal to the boundary S that parametrizes the lattice truncation boundary. The PML tensors properties depend on the local geometry via the two principal curvatures of S [80–82]. The boundary surface S is assumed (constructed) as doubly differentiable with non negative radii of curvature, otherwise dynamic instabilities ensue during a marching-on-time evolution [83].From (43), the PML also admits a straightforward interpretation as a complexification of the metric [38, 84]. As a result, the use of differential forms readily unifies the Maxwellian and non-Maxwellian PML formulations because the metric is explicitly factored out into the Hodge star operators—any transformation the metric corresponds, dually, to a transformation on the Hodge star operators that can be mimicked by modified constitutive relations [37]. In the differential forms framework, the PML is obtained by a mapping on the Hodge star operators: ⋆ϵ→⋆~ϵ and ⋆μ-1→⋆~μ-1 induced by the complexification of the metric. The resulting differential forms inside the PML, E~, D~, H~, and B~ therefore obey “modified” Hodge relations D~=⋆~ϵE~ and B~=⋆~μ-1H~, but identical premetric equations (10) and (11). In other words, (10) and (11) are invariant under the transformation (43) [38, 84].


## 8. Implementation of Space Charge Effects

In many applications related to plasma physics or electronic devices, it is necessary to include space charges (uncompensated charge effects) into lattice models of macroscopic Maxwell's equations. This is typically done by representing the charged plasma media using particle-in-cell (PIC) methods that track the individual particles on the lattice [85–87]. The field/charge interaction is then modeled by (i) interpolating lattice fields (cochains) to particle positions (gather step), (ii) advancing particle positions and velocities in time using equations of motion, and (iii) interpolating back charge densities and currents onto the lattice as cochains (scatter step). In general, the “particles” do not need to be actual individual particles, but can be a collection thereof (macroparticles). To put it simply, incorporation of space charges requires two extra steps during the field update in any marching-on-time algorithm, which transfer information from the instantaneous field distribution to the particle kinematic update and vice versa. Conventionally, this information transfer relies on spatial interpolations that often violates the charge continuity equation and, as a result, leads to spurious charge deposition on the lattice nodes. On regular lattices, this problem can be corrected, for example, using approaches that either subtract a static solution (charges) from the electric field solution (Boris/DADI correction) or directly subtract the residual error on the Gauss law (Langdon-Marder correction) at each time step [88]. On irregular lattices, additional degrees of freedom can be added as coupled elliptical constraints to produce an augmented Lagrange multiplier system [89]. All these approaches necessitate changes on the original equations, while still allowing for small violations on charge conservation. In contrast, Whitney forms provide a direct route to construct gather and scatter steps that satisfy charge conservation exactly even on unstructured lattices [90, 91], as explained next. To conform to the vast majority of the plasma and electronic devices literature, we once more restrict ourselves here to the 3+1 setting, even though a four-dimensional analysis in Minkowski space would have provided a more succinct discussion.For the gather step, Whitney forms can be used to directly compute (interpolate) the fields at any location from the knowledge of its cochain values, such as in (16), for example. For the scatter step, charge movement can be modeled as the Hodge-dual of the current 2-form J, that is, as the 1-form ⋆J which can be expanded in terms of Whitney 1-forms on the primal lattice. Here, ⋆ represents again the spatial Hodge star in three dimensions distilled from macroscopic constitutive properties. The Hodge-dual current associated to an individual point charge can be expressed as ⋆J=qv♭, where q is the charge value, v is the associated velocity vector, and ♭ is the “flat” operator or index-lowering canonical isomorphism that maps a vector to a 1-form, given by the Euclidean metric.  Similarly, point charges can be encoded as the Hodge-dual of the charge density 3-form ρ, that is, as the 0-form ⋆ρ, which can be expanded in terms of Whitney 0-forms on the primal lattice. These two Whitney maps are linked in such a way that the rate of change on the value of the 0-cochain representing ⋆ρ at a node is associated to the presence of a 1-cochain representing ⋆J along the edges that touch that particular node, leading to exact charge conservation at the discrete level. To show this, consider for simplicity the two-dimensional case of a point charge q moving from point x(s) to point x(f) during a time interval τ inside a triangular cell with nodes σ0,0, σ0,1, and σ0,2, or simply 0, 1, and 2. At any point x inside this cell, the 0-form ⋆ρ can be scattered to these three adjacent nodes via
(44)⋆ρ=q∑i=13〈x,ωi0〉ωi0,
where we are again using the short-hand ω0[σ0,i]=ωi0, and the brackets represent the pairing expressed by (1). In this case, p=0 and the pairing integral in (1) reduces to a function evaluation at a point. Since Whitney 0-forms are equal to the barycentric coordinates associated of a given node, that is, 〈x,ωi0〉=λi(x), we have the scattered charge qλis≐qλi(x(s)) on node i for a charge q at x(s), and, similarly, the scattered charge qλif on node i for a charge q at x(f). The rate of scattered charge variation on a given node i is, therefore, equal to q˙(λif-λis), where q˙=q/τ.Duringτ, the particle travels through a path ℓ from x(s) to x(f), and the corresponding ⋆J can be expanded as a sum of Whitney 1-forms ωij-1 associated to the three adjacent edges ij¯=01¯,12¯,20¯, that is,
(45)⋆J=q˙∑ij¯〈ℓ,ωij¯1〉ωij¯1.
The coefficients 〈ℓ,ωij-1〉 represent the (oriented) current flow along the associated oriented edge, that is, the cochain representation of ⋆J along edge ij¯. Using (13), the sum of the total current magnitude scattered along edges 01¯ and 20¯ that flows into node 0 is, therefore,
(46)q˙(-〈ℓ,ω01¯1〉+〈ℓ,ω20¯1〉)=q˙∫ℓ(-ω01¯1+ω20¯1).
Using ωij-1=λidλj-λjdλi and λ1+λ2+λ3=1, the above reduces to
(47)q˙∫ℓdλ0=q˙(λ0f-λ0s),
which exactly matches the rate of scattered charge variation on node 0 obtained before. It is clear that similar equalities hold for nodes 1 and 2. More fundamentally, these equalities are a direct consequence of the structural property (15).


## 9. Outline of Related Discretization Methods

We outline below various discretization programs that rely, one way or another, on tenets exposed above. The delineation is informed mostly by applications related to electrodynamics. As expected, this delineation is not too sharp because the programs share much in common.


### 9.1. Finite-Difference Time-Domain Method

In cubical lattices, the (lowest-order) Whitney forms can be represented by means of a product of pulse and “rooftop” functions on the three Cartesian coordinates [92]. This choice, together with the use of low-order quadrature rules to compute the Hodge star integrals in (19), leads to diagonal matrices [⋆ϵ], [⋆μ-1], and, consequently, also diagonal [⋆ϵ]-1, [⋆μ-1]-1 and sparse [Υ] so that an ultralocal equation results for (26). In this fashion, one obtains a “matrix-free” algorithm where no linear algebra is needed during a marching-on-time solution for the fields. This prescription exactly recovers the Yee's scheme [50] that forms the basis for the celebrated finite-difference time-domain (FDTD) method (see [51, 93] and references therein). FDTD adopts the simplest explicit, energy-conserving (symplectic) time-discretization for (23) and (26), which can be constructed by staggering the electric and magnetic fields in time and replacing time derivatives by central differences. This results in the following “leap-frog” marching-on-time scheme:
(48)Bin+1/2=Bin-1/2-Δt(∑jCij1Ejn),Ein+1=Ein+Δt(∑jΥijBjn+1/2),
where the superscript n denotes the time-step index, and Δt is the time increment (assumed uniform for simplicity). The staggering of the fields in both space and time is consistent with the presence of two staggered hypercubical spacetime lattices [48, 94] that can be viewed as prismatic extrusions, along the time coordinate, from the two (dual) staggered spatial lattices. The staggering in time also provides a O(Δt2) truncation error. Higher-order FDTD schemes with faster convergence to the continuum can be constructed by using less local approximations for the spatial derivatives (or equivalently less sparse [⋆ϵ] and [Υ]) and/or for the time derivatives [95–97].


### 9.2. Finite-Integration Technique

The finite-integration technique (FIT) [98–100] is closely related to FDTD, with the main distinction being that in FIT the discretized equations are derived from the integral form of Maxwell's equations applied to every cell. Assuming piecewise constant fields over each cell, the latter is equivalent to applying the (discrete version) of the generalized Stokes' theorem to the cochains in (23) and (24). Another difference is that the incidence matrices and material (Hodge star) matrices are treated separately in FIT. In other words, metric-free and metric-dependent parts of the equations are factorized a priori, in a manner akin to that exposed in Sections 3 and 4. Like FDTD, FIT is based on dual staggered lattices and, for cubical lattices, it turns out that the lowest-order FIT is algorithmically equivalent to the lowest-order FDTD. The spatial operators in FIT can all be viewed as discrete incarnations of the exterior derivative for the various p, and as such, the exact sequence property of the underlying de Rham complex is automatically enforced by construction [55]. Because of this, it could perhaps be claimed that FIT provides a more systematic route for generalizations to irregular lattices than Yee's FDTD. Historically, FIT generalizations to irregular lattices have relied on the use of either projection operators [101] or Whitney forms [102] to construct discrete versions of the Hodge star operators (or their procedural equivalents); however, these generalizations do not necessarily recover the specific form of the discrete Hodge matrix elements expressed in (19).


### 9.3. Cell Method

Another related discretization method, based on general principles originally put forth in [47–49], is the Cell method [103–108]. Even though this method does not rely on Whitney forms for constructing discrete Hodge star operators (other geometrically based constructions are instead used), it is nevertheless still based upon the use of “domain-integrated” discrete variables that conform to the notion of discrete differential forms or cochains of various degrees and, as such, it is naturally suited for irregular lattices. The Cell method also employs metric-free discrete operators that satisfy the exactness property of the de Rham complex and make explicit use of a dual lattice (but not necessarily barycentric) motivated by the notion of inner and outer orientations. The relationships between the various discrete operators and “domain-integrated” field quantities (cochains) in the Cell method are built into general classification diagrams referred to as “Tonti diagrams” that reproduce correct commuting diagram properties of the underlying operators [47, 48].


### 9.4. Mimetic Finite Differences

“Mimetic” finite-difference methods, originally developed for nonorthogonal hexahedral structured lattices (“tensor-product grids”) and later extended for irregular and polyhedral lattices [109–118], also share many of the properties exposed above. The thrust here is towards the construction of discrete versions of the differential operators divergence, gradient, and curl of vector calculus having “compatible” (in the sense of the exactness property of the underlying de Rham complex) domains and ranges and such that the resulting discrete equations exactly satisfy discrete conservation laws. In three dimensions, this naturally leads to the definition of three “natural” operators and three “adjoint” operators that can be associated with exterior derivative d and the codifferential δ, respectively, for p=1,2,3 (although the exterior calculus terminology is often not used explicitly in this context). Metric aspects are not factored out into Hodge star operators because the latter do not appear explicitly in the formulation; instead, their procedural analogues are embedded into the definition of the discrete differential operators themselves through a properly defined set of discrete inner products for discrete scalar and vector fields. In mimetic finite differences, the discrete analogues of the codifferential operator δ are full matrices, and the matrix-free character of FDTD is lacking even on orthogonal lattices. In spite of that, an obvious advantage of mimetic finite differences versus conventional FDTD is that the former methodology provides a more natural extension to nonorthogonal and irregular lattices. Note that higher-order versions of mimetic finite differences also exist [119, 120].


### 9.5. Compatible Discretizations and Finite-Element Exterior Calculus

In recent years, much attention has been devoted to the development of “compatible discretizations,” an umbrella term used to denote spatial discretizations of partial differential equations seeking to provide finite-element spaces that reproduce the exactness of the underlying de Rham complex (or the correct cohomology in topologically nontrivial domains) [121–126]. In this program, Whitney forms play a role of providing “conforming” vector-valued functional (finite-element) spaces of Sobolev type. Specifically, Whitney 1-forms recover the space of “Nedelec edge-elements” or curl-conforming Sobolev space H(curl,Ω) [127] and Whitney 2-forms recover the space of “Raviart-Thomas elements” or div-conforming Sobolev space H(div,Ω) [128]. In this regard, a relatively new advance here has been the development of new finite-element spaces, beyond those provided by Whitney forms, based on the Koszul complex [129]. The latter is key for the stable discretization of elastodynamics, which had been an outstanding problem for many decades. An excellent first-hand summary of these advances is provided in [130]. Another recent comparable approach aimed at the stable discretization of elastodynamics using bundle-valued discrete differential forms is described in [131].We should note that the link between stability conditions of some mixed finite-element methods [127] and the complex of Whitney forms has a long history in the context of electrodynamics. This link was first established in [55, 132], and further explored, for example, in [18, 19, 21, 23, 32, 36, 61, 133–136].


### 9.6. Discrete Exterior Calculus

The “discrete exterior calculus” (DEC) is another discretization program aimed at developing ab initio consistent discrete models to describe field theories [91, 137–141]. The main thrust of this program is not tied to any particular field theory but rather seeks to develop fundamental discrete tools (field variables, operators), amenable to tackle a whole gamut of theories (electrodynamics, fluid dynamics, elastodynamics, etc.). This discretization program recognizes the key role played by discrete differential forms as well as the need to define primal and dual cell complexes. There is a perceived focus on the use of circumcentric dual lattices as opposed to barycentric duals [138, 139] (even though the former does not admit a metric-free construction), and the program does not emphasize the role of Whitney forms (at least on its earlier stages). On the other hand, it recognizes the need to address group-valued differential forms, as well as the mathematical objects that exist on the dual-bundle space together with the associated operators (such as contractions and Lie derivatives), in connection to discrete problems in mechanics, optimal control, and computer vision/graphics [137]. A recent discussion on obstacles associated with some of the DEC underpinnings is provided in [142].


## 9.1. Finite-Difference Time-Domain Method

In cubical lattices, the (lowest-order) Whitney forms can be represented by means of a product of pulse and “rooftop” functions on the three Cartesian coordinates [92]. This choice, together with the use of low-order quadrature rules to compute the Hodge star integrals in (19), leads to diagonal matrices [⋆ϵ], [⋆μ-1], and, consequently, also diagonal [⋆ϵ]-1, [⋆μ-1]-1 and sparse [Υ] so that an ultralocal equation results for (26). In this fashion, one obtains a “matrix-free” algorithm where no linear algebra is needed during a marching-on-time solution for the fields. This prescription exactly recovers the Yee's scheme [50] that forms the basis for the celebrated finite-difference time-domain (FDTD) method (see [51, 93] and references therein). FDTD adopts the simplest explicit, energy-conserving (symplectic) time-discretization for (23) and (26), which can be constructed by staggering the electric and magnetic fields in time and replacing time derivatives by central differences. This results in the following “leap-frog” marching-on-time scheme:
(48)Bin+1/2=Bin-1/2-Δt(∑jCij1Ejn),Ein+1=Ein+Δt(∑jΥijBjn+1/2),
where the superscript n denotes the time-step index, and Δt is the time increment (assumed uniform for simplicity). The staggering of the fields in both space and time is consistent with the presence of two staggered hypercubical spacetime lattices [48, 94] that can be viewed as prismatic extrusions, along the time coordinate, from the two (dual) staggered spatial lattices. The staggering in time also provides a O(Δt2) truncation error. Higher-order FDTD schemes with faster convergence to the continuum can be constructed by using less local approximations for the spatial derivatives (or equivalently less sparse [⋆ϵ] and [Υ]) and/or for the time derivatives [95–97].


## 9.2. Finite-Integration Technique

The finite-integration technique (FIT) [98–100] is closely related to FDTD, with the main distinction being that in FIT the discretized equations are derived from the integral form of Maxwell's equations applied to every cell. Assuming piecewise constant fields over each cell, the latter is equivalent to applying the (discrete version) of the generalized Stokes' theorem to the cochains in (23) and (24). Another difference is that the incidence matrices and material (Hodge star) matrices are treated separately in FIT. In other words, metric-free and metric-dependent parts of the equations are factorized a priori, in a manner akin to that exposed in Sections 3 and 4. Like FDTD, FIT is based on dual staggered lattices and, for cubical lattices, it turns out that the lowest-order FIT is algorithmically equivalent to the lowest-order FDTD. The spatial operators in FIT can all be viewed as discrete incarnations of the exterior derivative for the various p, and as such, the exact sequence property of the underlying de Rham complex is automatically enforced by construction [55]. Because of this, it could perhaps be claimed that FIT provides a more systematic route for generalizations to irregular lattices than Yee's FDTD. Historically, FIT generalizations to irregular lattices have relied on the use of either projection operators [101] or Whitney forms [102] to construct discrete versions of the Hodge star operators (or their procedural equivalents); however, these generalizations do not necessarily recover the specific form of the discrete Hodge matrix elements expressed in (19).


## 9.3. Cell Method

Another related discretization method, based on general principles originally put forth in [47–49], is the Cell method [103–108]. Even though this method does not rely on Whitney forms for constructing discrete Hodge star operators (other geometrically based constructions are instead used), it is nevertheless still based upon the use of “domain-integrated” discrete variables that conform to the notion of discrete differential forms or cochains of various degrees and, as such, it is naturally suited for irregular lattices. The Cell method also employs metric-free discrete operators that satisfy the exactness property of the de Rham complex and make explicit use of a dual lattice (but not necessarily barycentric) motivated by the notion of inner and outer orientations. The relationships between the various discrete operators and “domain-integrated” field quantities (cochains) in the Cell method are built into general classification diagrams referred to as “Tonti diagrams” that reproduce correct commuting diagram properties of the underlying operators [47, 48].


## 9.4. Mimetic Finite Differences

“Mimetic” finite-difference methods, originally developed for nonorthogonal hexahedral structured lattices (“tensor-product grids”) and later extended for irregular and polyhedral lattices [109–118], also share many of the properties exposed above. The thrust here is towards the construction of discrete versions of the differential operators divergence, gradient, and curl of vector calculus having “compatible” (in the sense of the exactness property of the underlying de Rham complex) domains and ranges and such that the resulting discrete equations exactly satisfy discrete conservation laws. In three dimensions, this naturally leads to the definition of three “natural” operators and three “adjoint” operators that can be associated with exterior derivative d and the codifferential δ, respectively, for p=1,2,3 (although the exterior calculus terminology is often not used explicitly in this context). Metric aspects are not factored out into Hodge star operators because the latter do not appear explicitly in the formulation; instead, their procedural analogues are embedded into the definition of the discrete differential operators themselves through a properly defined set of discrete inner products for discrete scalar and vector fields. In mimetic finite differences, the discrete analogues of the codifferential operator δ are full matrices, and the matrix-free character of FDTD is lacking even on orthogonal lattices. In spite of that, an obvious advantage of mimetic finite differences versus conventional FDTD is that the former methodology provides a more natural extension to nonorthogonal and irregular lattices. Note that higher-order versions of mimetic finite differences also exist [119, 120].


## 9.5. Compatible Discretizations and Finite-Element Exterior Calculus

In recent years, much attention has been devoted to the development of “compatible discretizations,” an umbrella term used to denote spatial discretizations of partial differential equations seeking to provide finite-element spaces that reproduce the exactness of the underlying de Rham complex (or the correct cohomology in topologically nontrivial domains) [121–126]. In this program, Whitney forms play a role of providing “conforming” vector-valued functional (finite-element) spaces of Sobolev type. Specifically, Whitney 1-forms recover the space of “Nedelec edge-elements” or curl-conforming Sobolev space H(curl,Ω) [127] and Whitney 2-forms recover the space of “Raviart-Thomas elements” or div-conforming Sobolev space H(div,Ω) [128]. In this regard, a relatively new advance here has been the development of new finite-element spaces, beyond those provided by Whitney forms, based on the Koszul complex [129]. The latter is key for the stable discretization of elastodynamics, which had been an outstanding problem for many decades. An excellent first-hand summary of these advances is provided in [130]. Another recent comparable approach aimed at the stable discretization of elastodynamics using bundle-valued discrete differential forms is described in [131].We should note that the link between stability conditions of some mixed finite-element methods [127] and the complex of Whitney forms has a long history in the context of electrodynamics. This link was first established in [55, 132], and further explored, for example, in [18, 19, 21, 23, 32, 36, 61, 133–136].


## 9.6. Discrete Exterior Calculus

The “discrete exterior calculus” (DEC) is another discretization program aimed at developing ab initio consistent discrete models to describe field theories [91, 137–141]. The main thrust of this program is not tied to any particular field theory but rather seeks to develop fundamental discrete tools (field variables, operators), amenable to tackle a whole gamut of theories (electrodynamics, fluid dynamics, elastodynamics, etc.). This discretization program recognizes the key role played by discrete differential forms as well as the need to define primal and dual cell complexes. There is a perceived focus on the use of circumcentric dual lattices as opposed to barycentric duals [138, 139] (even though the former does not admit a metric-free construction), and the program does not emphasize the role of Whitney forms (at least on its earlier stages). On the other hand, it recognizes the need to address group-valued differential forms, as well as the mathematical objects that exist on the dual-bundle space together with the associated operators (such as contractions and Lie derivatives), in connection to discrete problems in mechanics, optimal control, and computer vision/graphics [137]. A recent discussion on obstacles associated with some of the DEC underpinnings is provided in [142].



---

*Source: 487270-2013-02-03.xml*