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.

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 [36]. 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 Appendix B.

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 [1129]. 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, 3339], 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 differential -forms on a smooth connected manifold as . From a geometric perspective, a differential -form can be viewed as an oriented -dimensional density, or an object naturally associated with -dimensional domains of integration such that the lattice contraction (“pairing”) below: gives a real number (in our context) for each choice of [22]. On a lattice is restricted to be a union of elements from the finite set of -dimensional -oriented lattice elements, which we denote by . These are collectively called “-chains.” In four dimensions for example, they correspond to the possible unions of elements from the set of vertices (nodes) , edges (links) , facets (plaquettes) , volume cells (voxels) , and hypervolume cells , for , 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 such that defines its action on the basis of -chains. Note that we use to denote the space dual to , that is, the space -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 is called the de Rham map [22].

The basic differential operator of exterior calculus is the exterior derivative , applicable to any number of dimensions. The discretization of on a general irregular lattice can be effected by a straightforward application of the generalized Stokes' theorem [22]: with in . In the above, is the boundary operator, which simply maps a -dimensional lattice element to the set of -dimensional lattice elements that comprise its boundary, preserving orientation. This theorem sets as the formal adjoint of in terms of the pairing given in (1), that is, . Computationally, the boundary operator can be implemented by means of incidence matrices [22, 29, 40] such that where the indices and run over all - and -dimensional lattice elements, respectively. The incidence matrix entries are such that for all , with sign determined by the relative orientation of lattice elements and . 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 for all , so that the derivative operation is replaced by a proper sum over . On the lattice, the nilpotency of the operators [41] is recovered by the constraint [22]: for all and .

3. Example: Lattice Electrodynamics

We write Maxwell's equations in a four-dimensional Lorentzian manifold as [34]: where is the four-dimensional exterior derivative, and 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 -forms to -forms, and more generally forms to forms in a -dimensional manifold, and, as mentioned before, depends on the metric of [22, 23, 34, 35, 4245]. The above equations are complemented by the relation , which indicates that and are “Hodge duals” of each other.

3.1. Primal and Dual Lattices

Since and 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: is an “ordinary” (or “nontwisted”) differential form, whereas (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, 4648]. 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 is associated with the set of plaquettes on the “ordinary cell complex” , thus belonging to , while the twisted forms and are associated with the set of plaquettes on the “twisted cell complex” [22, 27, 48, 49], thus belonging to . Consequently, we also have two sets of incidence matrices and , 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 duality . 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) [5052].

3.2. 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 setting, one can decompose and as and the source density as where is the wedge product, and are the electric intensity and magnetic intensity 1-forms on and , respectively, and are the electric flux and magnetic flux 2-forms on and , respectively, is the electric current density 2-form on , and is the electric charge density 3-form on (corresponding assignments for the and cases are provided in [32]). As a result, Maxwell's equations reduce to representing Faraday's and Ampere's laws, respectively. Here, 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 now involving two Hodge star maps in three-dimensional space: and . On the lattice, we have the corresponding discrete counterparts: and . The subscripts and in and 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: , , , , , and , where the indices run over the respective basis of -chains in either or , . With the exception of Appendix A, we restrict ourselves to the setting throughout the remainder of this paper.

4. Casting the Metric on a Lattice

4.1. Whitney Forms

The Whitney map 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, 5561] 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, -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 -simplex , the (lowest-order) Whitney form is given by where , , are the barycentric coordinates associated with . In the case of a -simplex (node), (13) reduces to .

From its definition, it is clear that Whitney forms have compact support. Among its important structural properties are where is the Kronecker delta, which is simply a restatement of , and where 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, 6468].

Using the short-hand , we can write the following expansions for and in a irregular simplicial lattice, in terms of its cochain representations: 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 Using (16), the lattice Hamiltonian assumes the expected quadratic form: where we immediately identify the symmetric positive definite matrices: as the discrete realization of the Hodge star operator(s) on a simplicial lattice [23, 69] so that From the above, the Hamiltonian can be also expressed as

4.3. Symplectic Structure and Dynamic Degrees of Freedom

The Hodge star matrices and have different sizes. The number of elements in is equal to , whereas the number of elements in is equal to . In other words, , 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 as 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 , we still have , where 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 of , and using the fact that from (14), we get so that where the index runs over all faces of the primal lattice. On the dual lattice , we can similarly contract both sides of (11) with every dual face to get where now the index runs over all faces of the dual lattice. Using (20) and the fact that, in three-dimensions [22] (up to possible boundary terms ignored here), we can write the last equation in terms of primal lattice quantities as or, by using the inverse Hodge star matrix , as with The matrix can be viewed as the discrete realization, for , of the codifferential operator that maps -forms to -forms [35].

Since the continuum operators and are local [46] and, as seen, Whitney forms (13) have local support, it follows that the matrices and are sparse, indicative of an ultralocal coupling (in the terminology of [70]). In contrast, the numerical inverse 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 [29, 71]. The sparsity pattern of encodes the nearest-neighbor edge information of the mesh and, consequently, the sparsity pattern of likewise encodes successive “-level” neighbors. The latter sparsity patterns can be used to build, quite efficiently, sparse approximations for , 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 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 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: where stands for the spatial Hodge star operator (distilled from constitutive material properties), and is the dual element to on the barycentric dual lattice. The operator is such that where is the two-norm of , and 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, with and with , so that . This construction is detailed in [24].

5.3. Galerkin Duality

Even though we have chosen to assign and to the primal (simplicial) lattice, and consequently , , , and to the dual (nonsimplicial) lattice, the reverse is equally possible. In this case, the fields , become associated to a simplicial lattice and hence can be expressed in terms of Whitney forms; the expressions dual to (16) are now with sums running over primal edges and primal faces, respectively, and where with and the two Hodge star maps now used are such that, in the continuum, and , and, on the lattice, and . 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 any -form , we can write where is a harmonic form [31]. This Hodge decomposition is unique. In the particular case of the -form , we have where is a -form and is a -form, with representing the static field, the dynamic field, and the harmonic field component (if any). In a contractible domain, is identically zero and the Hodge decomposition simplifies to 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 are associated to the facets of the primal lattice. Consequently, we have from (35) that where is the number of primal nodes, the number of primal edges, and the number of primal facets, with superscript standing for boundary (fixed) elements and 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] where 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:

Combining (37) and (38), we have From the Hodge decomposition (35), we see that is

Note that the divergence free condition produces one constraint on the 2-form for each volume element. This constraint also spans the whole lattice boundary. The total number of the constrains for is, therefore, . Consequently, we have so that This discussion can be generalized to lattices on noncontractible domains with any number of holes (genus), where the identity is also satisfied [31]. Moreover, from Hodge star isomorphism, we have and .

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 formula is associated to a corresponding term in , that is, the number of edges corresponds to the dimension of the space of lattice -forms , which is the sum of the number of nodes (dimension of the space of discrete -forms ), the number of faces (dimension of the space of discrete -forms ), and the number of holes (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 [7578]. 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 () to a complex coordinate stretching of the coordinate space (or an analytic continuation to a complex-valued coordinate space) [7678] 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 where is the so-called complex stretching variable written as with and (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 that parametrizes the lattice truncation boundary. The PML tensors properties depend on the local geometry via the two principal curvatures of [8082]. The boundary surface 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 induced by the complexification of the metric. The resulting differential forms inside the PML, , , , and therefore obey “modified” Hodge relations and , 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 [8587]. 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 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 , that is, as the 1-form 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 , where is the charge value, 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 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 moving from point to point during a time interval inside a triangular cell with nodes , , and , or simply , , and . At any point inside this cell, the 0-form can be scattered to these three adjacent nodes via where we are again using the short-hand , and the brackets represent the pairing expressed by (1). In this case, 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, , we have the scattered charge on node for a charge at , and, similarly, the scattered charge on node for a charge at . The rate of scattered charge variation on a given node is, therefore, equal to , where .

During , the particle travels through a path from to , and the corresponding can be expanded as a sum of Whitney 1-forms associated to the three adjacent edges , that is, The coefficients represent the (oriented) current flow along the associated oriented edge, that is, the cochain representation of along edge . Using (13), the sum of the total current magnitude scattered along edges and that flows into node is, therefore, Using and , the above reduces to which exactly matches the rate of scattered charge variation on node 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).

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 , , and, consequently, also diagonal , 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: where the superscript denotes the time-step index, and 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 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 [9597].

9.2. Finite-Integration Technique

The finite-integration technique (FIT) [98100] 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 , 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 [4749], is the Cell method [103108]. 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 [109118], 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 and the codifferential , respectively, for (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) [121126]. 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 [127] and Whitney 2-forms recover the space of “Raviart-Thomas elements” or div-conforming Sobolev space H [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, 133136].

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, 137141]. 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].


A. Differential Forms and Lattice Fermions

Differential -forms can be viewed as antisymmetric covariant tensor fields on rank . Therefore, the ingredients discussed above are applicable to any antisymmetric tensor field theory, including non-Abelian gauge field theories and even topological field theories such as Chern-Simons theory [72]. However, for (Dirac) fermion fields the situation is different and, at first, it would seem unclear how differential forms could be used to describe spinors. Nevertheless, a useful connection can indeed be established [1, 16, 143]. To briefly address this point, we consider the lattice transcription of the (one-flavor) Dirac equation here.

Needless to say, the topic of lattice fermions is vast and we cannot do much justice to it here; we focus only on aspects that are more germane to main theme of this paper. In accordance to the related literature on lattice fermions, we work on Euclidean spacetime with in this appendix, and adopt the repeated index summation convention with , as coordinate indices, where is a point in four-dimensional space.

It is well known that fermion fields defy a lattice description with local coupling that gives the correct energy spectrum in the limit of zero lattice spacing and the correct chiral invariance [144]. This is formally stated by the no-go theorem of Nielsen-Ninomiya [145] and is associated to the well-known “fermion-doubling” problem [146]. A perhaps less known fact is that it is possible to arrive at a “geometrical” interpretation of the source of this difficulty by considering the “generalization” of the Dirac equation given by the Dirac-Kähler equation: The square of the Dirac-Kähler operator can be viewed as the counterpart of the Dirac operator in the sense that recovers the Laplacian operator in the same fashion as the Dirac operator squared does, that is, , where represents Euclidean gamma matrices.

The Dirac-Kähler equation admits a direct transcription on the lattice because both the exterior derivative and the codifferential can be simply replaced by its lattice analogues, as discussed before. However, for the Dirac equation the analogy has to further involve the relationship between the 4-component spinor field and the object . This relationship was first established in [16, 17] for hypercubic lattices and later extended to nonhypercubic lattices in [10, 147]. The analysis of [16, 17] has shown that can be represented by a 16-component complex-valued inhomogeneous differential form: where is a (1-component) scalar function of position or 0-form, is a (4-component) 1-form, and likewise for representing -, -, and -forms with -, -, and -components, respectively. By employing the following Clifford algebra product: as using the anticommutative property of the exterior product , we have which exactly matches the anticommutator result of the matrices, . This suggests that plays the role of the matrix in the space of inhomogeneous differential forms with Clifford product [148], that is, keeping in mind that while acts on spinors, acts on inhomogeneous differential forms. This analysis leads to a “geometrical” interpretation of the popular Kogut-Susskind staggered lattice fermions [149, 150] because the latter can be made identical to lattice Dirac-Kähler fermions after a simple relabeling of variables [17].

The 16-component object can be viewed as a matrix that produces a fourfold degeneracy with respect to the Dirac equation for . This degeneracy is actually not a problem in the continuum because there is a well-defined procedure to extract the 4-components of from those of [16, 17], whereby the 16 scalar equations encoded by (A.1) all reduce to the same copy of the four equations encoded by the standard Dirac equation. This procedure is performed by a set of “projection operators” that form a group [16, 151]. On the lattice, however, the operators and , as well as (which plays a role on the space of inhomogeneous differential forms analogous to that of on the space of spinors [152]), behave in such a way that their action leads to lattice translations. This is because cochains with different necessarily live on different lattice elements and also because is a map between different lattice elements. As a consequence, the product operation of such “group” is not closed anymore. This nonclosure also stems from the fact that the lattice operators and do not satisfy Leibnitz's rule [148]. Because of this, the degeneracy of the Dirac equation on the lattice is present at a more fundamental level and is harder to extricate using the Dirac-Kähler description than the analogous degeneracy in the continuum. In this regard, a new approach to identify the extraneous degrees of freedom away from the continuum was recently described in [153]. In addition, a split-operator approach to solve Dirac equation based on the methods of characteristics that purports to avoid fermion doubling while maintaining chiral symmetry on the lattice was very recently put forth in [154]. This approach preserves the linearity of the dispersion relation by a splitting of the original problem into a series of one-dimensional problems and the use of a upwind scheme with a Courant-Friedrichs-Lewy (CFL) number equal to one, which provides an exact time evolution (i.e., with no numerical dispersion effects) along each reduced one-dimensional problem. The main (practical) obstacle in this case is the need to use very small lattice elements.

B. Classification of Inconsistencies in Naïve Discretizations

We provide below a rough classification scheme of inconsistencies arising from naïve discretizations of the differential calculus on irregular lattices.

(i) Premetric Inconsistencies of First Kind

We call premetric inconsistencies of the first kind those that are related to the primal or dual lattices taken as separate objects and that occur when the discretization violates one or more properties of the continuum theory that is invariant under homeomorphisms–-for example, conservations laws that relate a quantity on a region with an associated quantity on the boundary of the region, (a topological invariant). Perhaps the most illustrative example is violation of “divergence-free” conditions caused by improper construction of incidence matrices, whereby the nilpotency of the (adjoint) boundary operator, , is not observed. This implies, in a dual fashion, that the identity is violated [22]. Stated in another way, the exact sequence property of the underlying de Rham differential complex is violated [155]. In practical terms, this leads to the appearance spurious charges and/or spurious (“ghost”) modes. As the classification suggests, these properties are not related to metric aspects of the lattice, but only to its “topological aspects,” that is, on how discrete calculus operators are defined vis-à-vis the lattice element connectivity. In more mathematical terms, one can say that the structure of the (co)homology groups of the continuum manifold is not correctly captured by the cell complex (lattice). We stress again that, given any dual lattice construction, premetric inconsistencies of the first kind are associated to the primal or dual lattice taken separately, and not necessarily on how they intertwine.

(ii) Premetric Inconsistencies of Second Kind

The second type of premetric inconsistency is associated to the breaking of some discrete symmetry of the Lagrangian. In mathematical terms, this type of inconsistency can occur when the bijective correspondence between -cells of the primal lattice and -cells of the dual lattice (an expression of Poincaré duality at the level of cellular homology [156], up to boundary terms) is violated. This is typified by “nonreciprocal” constructions of derivative operators, where the boundary operator effecting the spatial derivation on the primal lattice is not the dual adjoint (or the incidence matrix transpose) of the boundary operator on the dual lattice , for example, the identity (up to boundary terms) used to obtain (25) is violated. One basic consequence of this violation is that the resulting discrete equations break time-reversal symmetry. Consequently, the numerical solutions will violate energy conservation and produce either artificial dissipation or late-time instabilities [22]. Many algorithms developed over the years for hyperbolic partial differential equations do indeed violate these properties: they are dissipative and cannot be used for long integration times [157, 158].

It should be noted at this point that lattice field theories invariably break Lorentz covariance and many of the continuum Lagrangian symmetries and, as a result, violate conservation laws (currents) by virtue of Noether's theorem. For example, angular momentum conservation does not hold exactly on the lattice because of the lack of continuous rotational symmetry (note that discrete rotational symmetries can still be present). However, this latter type of symmetry breaking is of a fundamentally different nature because it is “controllable,” that is, their effect on the computed solutions is made arbitrarily small in the continuum limit. More importantly, discrete transcriptions of the Noether's theorem can be constructed for Lagrangian symmetries on a lattice [13, 159], to yield exact conservation laws of (properly defined) quantities such as discrete energy and discrete momentum [3].

(iii) Hodge Star Inconsistencies

In the third type of inconsistency, we include those that arise in connection with metric properties of the lattice. Because the metric is entirely encoded in the Hodge star operators [22, 42, 160], such inconsistencies can be simply understood as inconsistencies on the construction of discrete Hodge star operators (or their procedural analogues). For example, it is not uncommon for naïve discretizations in irregular lattices to yield asymmetric discrete Hodge operators, as noted in [161, 162]. Even if symmetry is observed, nonpositive definiteness might ensue that is often associated with portions of the lattice with highly skewed or obtuse cells [101]. Lack of either of these properties leads to unconditional instabilities that destroy marching-on-time solutions [22]. When very long integration times are needed, asymmetry in the discrete Hodge matrices can be a problem even if produced at the level of machine rounding-off errors.


The author thanks Weng Chew, Burkay Donderici, Bo He, and Joonshik Kim for discussions. The author also thanks the editorial board for the invitation to contribute with this paper.