Computation of Maxwell's equations on Manifold using DEC

TL;DR

This paper presents a numerical method using discrete exterior calculus to solve Maxwell's equations on manifolds, analyzing stability and errors, with implementation for simulating electromagnetic waves in vacuum.

Contribution

Introduces a discrete exterior calculus approach for Maxwell's equations on manifolds, including stability analysis and implementation details.

Findings

Stable numerical scheme for Maxwell's equations on manifolds

Error analysis of the discretization method

Simulation of TE/M waves in vacuum

Abstract

In this paper, the method of discrete exterior calculus for numerically solving Maxwell's equations in space manifold and the time is discussed, which is a kind of lattice gauge theory. The analysis of its stable condition and error is also accomplished. This algorithm has been implemented on C++ plateform for simulating TE/M waves in vacuum.

Figures

Equations

$$\begin{array}[]{lll}dF=0&&d\ast F=\ast J\end{array}$$

$$d\ast J=0.$$

$$D({\sigma }_0):=\bigcup _{{\sigma }_0\in {\sigma }_1\in \cdots \in {\sigma }_r}\mathrm{Int}(c({\sigma }_0)c({\sigma }_1)\cdots c({\sigma }_r)),$$

$$A:E\to \mathbb{R}.$$

$$F=dA:P\to \mathbb{R}.$$

$$dF=0.$$

$$\begin{array}[]{lll}L(A,J)&=&-\frac{1}{2}\langle dA,dA\rangle+\langle A,J\rangle,\end{array}$$

$$\begin{array}[]{lll}\langle dA,dA\rangle&:=&(A)_{1\times|E|}(d)_{|E|\times|F|}(\ast)_{|F|\times|F|}(d)^{T}_{|F|\times|E|}(A)^{T}_{|E|\times 1}\\ \langle A,J\rangle&:=&(A)_{1\times|E|}(\ast)_{|E|\times|E|}(J^{T})_{|E|\times 1}\end{array}$$

$$\begin{array}[]{lll}\partial_{A_{i}}L(A,J)&=&\partial_{A_{i}}(-\frac{1}{2}\langle dA,dA\rangle+\langle A,J\rangle)\\ &=&\partial_{A_{i}}(-\frac{1}{2}(A)_{1\times|E|}(d)_{|E|\times|F|}(*)_{|F|\times|F|}(d)^{T}_{|F|\times|E|}(A)^{T}_{|E|\times 1}\\ &&+(A)_{1\times|E|}(*)_{|E|\times|E|}(J^{T})_{|E|\times 1})\\ &=&-\frac{1}{2}(0,...,\underbrace{1}_{i},...,0)_{1\times|E|}(d)_{|E|\times|F|}(*)_{|F|\times|F|}(d)^{T}_{|F|\times|E|}(A)^{T}_{|E|\times 1}\\ &&-\frac{1}{2}(A)_{1\times|E|}(d)_{|E|\times|F|}(*)_{|F|\times|F|}(d)^{T}_{|F|\times|E|}(0,...,\underbrace{1}_{i},...,0)^{T}_{|E|\times 1}\\ &&+(0,...,\underbrace{1}_{i},...,0)_{1\times|E|}(*)_{|E|\times|E|}(J^{T})_{|E|\times 1}\\ &=&-(0,...,\underbrace{1}_{i},...,0)_{1\times|E|}(d)_{|E|\times|F|}(*)_{|F|\times|F|}(d)^{T}_{|F|\times|E|}(A)^{T}_{|E|\times 1}\\ &&+(0,...,\underbrace{1}_{i},...,0)_{1\times|E|}(*)_{|E|\times|E|}(J^{T})_{|E|\times 1}\end{array}$$

$$-(d{)}_{|E|\times |F|}(\ast {)}_{|F|\times |F|}(d{)}_{|F|\times |E|}^T(A{)}_{|E|\times 1}^T+(\ast {)}_{|E|\times |E|}(J^T{)}_{|E|\times 1}=0,$$

$$\delta F=J,$$

$$d^T\ast J=0.$$

$$A\to A+df$$

$$F\to F+d^{2}f=F,$$

$$⟨df,J⟩=(f{)}_{1\times |V|}(d^T{)}_{|V|\times |E|}(\ast {)}_{|E|\times |E|}(J^T{)}_{|E|\times 1}=0,$$

$$\begin{array}[]{lll}L(A+df,J)&=&-\frac{1}{2}\langle d(A+df),d(A+df)\rangle+\langle A+df,J\rangle=L(A,J).\end{array}$$

$$J=(-{\rho }_{e}dt,J_e)F^{n+\frac{1}{2}}=E^{n+\frac{1}{2}}\wedge dt+B^n\ast F^n=H^n\wedge dt-D^{n-\frac{1}{2}},$$

$$\begin{array}[]{lll}d_{s}B^{n}&=&0\\ d_{s}E^{n+\frac{1}{2}}\wedge d{t}&=&-d_{{t}}B^{n}\\ d^{T}_{s}D^{n-\frac{1}{2}}&=&\ast(\rho_{e}dt)^{n-\frac{1}{2}}\\ d^{T}_{s}H^{n}\wedge d{t}&=&d^{T}_{{t}}D^{n-\frac{1}{2}}+\ast{J}^{n}_{e},\end{array}$$

$$\begin{array}[]{lll}d_{{t}}B^{n}:=\dfrac{B^{n+1}-B^{n}}{\Delta t}\wedge{d{t}}&&d^{T}_{{t}}D^{n-\frac{1}{2}}:=\dfrac{D^{n+\frac{1}{2}}-D^{n-\frac{1}{2}}}{\Delta t}\wedge d{t}.\end{array}$$

$$\bar{J}=({\rho }_m,-J_m\wedge dt),$$

$$\begin{array}[]{lll}d_{s}B^{n}&=&\rho^{n}_{m}\\ d_{s}E^{n+\frac{1}{2}}\wedge d{t}&=&-d_{{t}}B^{n}-{J}^{n+\frac{1}{2}}_{m}\wedge d{t}\\ d^{T}_{s}D^{n-\frac{1}{2}}&=&\ast(\rho_{e}dt)^{n-\frac{1}{2}}\\ d^{T}_{s}H^{n}\wedge d{t}&=&d^{T}_{{t}}D^{n-\frac{1}{2}}+\ast{J}^{n}_{e},\end{array}$$

$$\begin{array}[]{lll}dF=\bar{J}&&d^{T}\ast F=\ast J,\end{array}$$

$$\begin{array}[]{lll}d\bar{J}=0&&d^{T}\ast J=0.\end{array}$$

$$d_s^T\ast ({\rho }_{e}dt)=0d_t^T\ast J_e=0d_s{\rho }_m=0d_t{J}_m\wedge dt=0,$$

$$d_t^T(\ast {\rho }_{e}d{t}^{n-\frac{1}{2}})-d_s^T\ast J_e^n=0-d_s{J}_m^{n+\frac{1}{2}}\wedge dt+d_t{\rho }_m^n=0.$$

$$\begin{array}[]{lll}d^{T}_{t}d^{T}_{s}D^{n-\frac{1}{2}}-d^{T}_{t}\ast(\rho_{e}dt)^{n-\frac{1}{2}}&=&-d^{T}_{t}\ast(\rho_{e}dt)^{n-\frac{1}{2}}-d^{T}_{s}(d^{T}_{s}H^{n}\wedge d{t}-\ast J^{n}_{e})\\ &=&0\\ d_{t}d_{s}B^{n}-d_{t}\rho^{n}_{m}&=&-d_{t}\rho^{n}_{m}+d_{s}(d_{s}E^{n+\frac{1}{2}}\wedge d{t}+{J}^{n+\frac{1}{2}}_{m}\wedge d{t})\\ &=&0.\end{array}$$

$$\begin{array}{rlr}& \frac{D_1^{n+\frac{1}{2}}-D_1^{n-\frac{1}{2}}}{\Delta t}+{J_e^n}_1=\frac{H_1^n-H_2^n}{|\ast e_1|}& \\ & -\frac{B_1^{n+1}-B_1^n}{\Delta t}=\frac{E_1^{n+\frac{1}{2}}|e_1|+E_2^{n+\frac{1}{2}}|e_2|+E_3^{n+\frac{1}{2}}|e_3|}{|P_1|}.& \end{array}\}$$

$$\begin{array}[]{lll}D_{i}=\varepsilon_{0}E_{i}&&B_{i}=\mu_{0}H_{i},\end{array}$$

$$\begin{array}{rlr}& {ϵ}_0\frac{E_1^{n+\frac{1}{2}}-E_1^{n-\frac{1}{2}}}{\Delta t}+J_{e1}^n=\frac{H_1^n-H_2^n}{|\ast e_1|}& \\ & {\mu }_0\frac{H_1^{n+1}-H_1^n}{\Delta t}=-\frac{E_1^{n+\frac{1}{2}}|e_1|+E_2^{n+\frac{1}{2}}|e_2|+E_3^{n+\frac{1}{2}}|e_3|}{|P_1|}& \end{array}\}\mathrm{TE}$$

$$\begin{array}{rlr}& {ϵ}_0\frac{E_1^{n+\frac{1}{2}}-E_1^{n-\frac{1}{2}}}{\Delta t}+J_{e1}^n=\frac{H_1^n-H_2^n}{|\ast e_1|}& \\ & {\mu }_0\frac{H_1^{n+1}-H_1^n}{\Delta t}+J_{m1}^{n+\frac{1}{2}}=-\frac{E_1^{n+\frac{1}{2}}|e_1|+E_2^{n+\frac{1}{2}}|e_2|+E_3^{n+\frac{1}{2}}|e_3|}{|P_1|}& \end{array}\}\mathrm{TE}$$

Full text

**Keywords: **Discrete exterior calculus, Discrete variation, Maxwell’s equations, Lattice gauge theory.

PACS(2010): 41.20.Jb, 02.30.Jr, 02.40.Sf, 02.60.Cb.

1 Introduction

Computational electromagnetism is concerned with the numerical study of Maxwell’s equations. The Yee scheme is known as finite difference time domain and is one of the most successful numerical methods, particularly in the area of microwave problems[1]. It preserves important structural features of Maxwell’s equations[2, 3, 4, 5, 6]. Bossavit et al present the Yee-like scheme and extend Yee scheme to unstructured grids. This scheme combines the best attributes of the finite element method (unstructured grids) and Yee scheme (preserving geometric structure)[7, 8]. Stern et al [9] generalize the Yee scheme to unstructured grids not just in space, but in 4-dimensional spacetime by discrete exterior calculus(DEC)[10, 11, 14, 15, 16, 20, 17, 18, 19, 12, 13]. This relaxes the need to take uniform time steps.

In this paper, we generalize the Yee scheme to the discrete space manifold and the time. The spacetime manifold used here is split as a product of $1$D time and $2$D or $3$D space manifold. The space manifold can be approximated by triangular and tetrahedrons depending on dimension, and the time by segments. So the spacetime manifold is approximated by prism lattice, on which the discrete Lorentz metric can be defined.

With the technique of discrete exterior calculus, the $\mathbb{R}$ value discrete connection, curvature and Bianchi identity are defined on prim lattice. With discrete variation of an inner product of discrete $1-$forms and their dual forms, the discrete source equation and continuity equation are derived.

• 2.

Those equations compose the discrete Maxwell’s equations in vacuum case, which just need the local information of triangulated manifold such as length and area.

The discrete Maxwell’s equations here can be re-grouped into two sets of explicit iterative schemes for TE and TM waves, respectively. Those schemes can directly use acute triangular, rectangular, regular polygon and their combination, which has been implemented on C++ plateform to simulate the electromagnetic waves propagation and interference on manifold.

2 DEC for Maxwell’s equations

Maxwell’s equations can be simply expressed once the language of exterior differential forms is used. The electric and magnetic fields are jointly described by a curvature $2-$form $F$ in a 4-D spacetime manifold. The Maxwell’s equations reduce to the Bianchi identity and the source equation

[TABLE]

where $d$ denotes the exterior differential operator, $\ast$ denotes the Hodge star operator, and 1-form $J$ is called the electric current form satisfying the continuity equation

[TABLE]

As the exterior derivative is defined on any manifold, the differential form version of the Bianchi identity makes sense for any 3D or 4D spacetime manifold, whereas the source equation is defined if the manifold is oriented and has a Lorentz metric. Now, we introduce the discrete counterpart of those differential geometric objects to derive the numerical computational schemes for Maxwell’s equations.

Discrete Lorentz metric

The spacetime manifold used here is split as a product of $1$D time and $2$D or $3$D space manifold. The $2$D or $3$D space manifold can be approximated by triangular or tetrahedrons, and the time by segments. The length of edge and area of triangular and volume of tetrahedrons gives the discrete Riemann metric on space grids. The metric on time grid is the minus of length square. The spacetime manifold is approximated by prism lattice, on which the discrete Lorentz metric can be defined as the product of discrete metric on space and time.

Discrete exterior calculus

A discrete differential $k$-form, $k\in\mathbb{Z}$, is the evaluation of the differential $k$-form on all $k$-simplices. Dual forms, i.e., forms that we evaluate on the dual cell. Suppose each simplex contains its circumcenter. The circumcentric dual cell $D(\sigma_{0})$ of simplex $\sigma_{0}$ is

[TABLE]

where $\sigma_{i}$ is all the simplices which contains $\sigma_{0}$,…, $\sigma_{i-1}$, and $c(\sigma_{i})$ is the circumcenter of $\sigma_{i}$.

The two operators in Eqs.(1) can be discretized as follows:

Discrete exterior differential operator $d$, this operator is the transpose of the incidence matrix of $k$-cells on $k+1$-cells.

• 2.

Discrete Hodge Star $\ast$, the operator scales the cells by the volumes of the corresponding dual and primal cells.

Discrete connection and curvature

Discrete connection $1-$form or gauge field $A$ assigns to each element in the set of edges $E$ an element of the gauge group $\mathbb{R}$:

[TABLE]

Discrete curvature $2-$form is the discrete exterior derivative of the discrete connection $1-$form

[TABLE]

The value of $F$ on each element in the set of triangular $P$ is the coefficient of Holonomy group of this face. The $2-$form $F$ automatically satisfies the discrete Bianchi identity

[TABLE]

Note that since the gauge group $\mathbb{R}$ used here is abelian, we need not pick a starting vertex for the loop. We may traverse the edges in any order, so long as taking orientations into account.

Discrete Maxwell’s equations

For source case, we need discrete current $1-$form $J$. Let $A=\sum\limits_{E}A_{i}$ and the Lagrangian functional be

[TABLE]

where

[TABLE]

Supposing that there is a variation of $A_{i}$, vanishing on the boundary, we have

[TABLE]

The Hamilton’s principle of stationary action states that this variation must equal zero for any such vary of $A_{i}$, implying the Euler-Lagrange equations

[TABLE]

which is the discrete source equation

[TABLE]

where $\delta=\ast^{-1}d^{T}\ast$. Since $(d^{T})^{2}=0$, the discrete continuity equation can express as:

[TABLE]

The equations of discrete Bianchi identity (2), source equation (3), and continuity equation (4) are called discrete Maxwell’s equations.

Discrete Gauge transformations

Discrete gauge transformations are maps

[TABLE]

for any $0-$form or scalar function $f$ on vertex. Since the discrete exterior derivative maps

[TABLE]

the discrete Maxwell’s equations(2-4) are invariant under discrete gauge transformations. Since the discrete continuity equation (4) ensures

[TABLE]

so we have

[TABLE]

That is to say the Lagrangian function is also invariant under discrete gauge transformations.

3 Explicit schemes

Schemes for TE wave

The discrete current $1-$form, discrete curvature $2-$form and its dual can be written as

[TABLE]

where $n$ and $n+\frac{1}{2}$ denote the coordinate of the time, $E=\sum\limits_{E}{E}_{i}e^{i}$ (electric field) is discrete $1-$form on space, $B=\sum\limits_{P}{B}_{i}P^{i}$ (magnetic field) is discrete $2-$form on space, $H=\sum\limits_{P}{H}_{i}*P^{i}$ (magnetizing field) is the dual of $B$ on space, $D=\sum\limits_{E}{D}_{i}*e^{i}$(electric displacement field) is the dual of $E$ on space, $\rho_{e}dt$ (charge density) is the discrete $1-$form on time, $J_{e}=\sum\limits_{E}J_{ei}e^{i}$ (electric current density) is the discrete $1-$form on space. The discrete Maxwell’s equations can be rewritten as

[TABLE]

where $d_{s}$, $d^{T}_{{s}}$ are the restriction of $d$, $d^{T}$ on space, and

[TABLE]

If allowing for the possibility of magnetic charges and current discrete $3-$form

[TABLE]

the symmetric scheme can be written as

[TABLE]

where

• $\rho_{m}=\sum\limits_{Tet}\rho_{mi}T^{i}$(magnetic charges) is discrete 3-form on space,
• ${J}_{m}=\sum\limits_{P}J_{mi}P^{i}$ (current) is discrete $2-$form on space.

The compact form of Eqs.(6) can be written as

[TABLE]

with discrete continuity equations or integrability conditions

[TABLE]

Proposition 3.1

If the initial condition satisfies the first and third equations in Eqs.(6), the solution of the second and fourth equations in Eqs.(6) automatically satisfy Eqs.(6).

Proof. Because the dimension of spacetime is $3+1$, therefore

[TABLE]

and the continuity equations can be reduced to

[TABLE]

So we have

[TABLE]

$\Box$

Now we show the scheme (5) on the product of 2D discrete space manifold and time. The second and fourth equations in Eqs.(5) based on Fig.1 are

[TABLE]

where $|~{}|$ denotes the measure of forms and dual. The summation on the right is orient, that is to say, inverse the orientation of $e_{i}$, then multiply $-1$ with $\bar{E}_{i}$. Eqs.(7) can be implemented on 2D discrete manifold directly(see Fig.1). Eq.(7) on rectangular gird is just the Yee scheme.

[TABLE]

In the absence of magnetic or dielectric materials, the relations are simple:

[TABLE]

where $\varepsilon_{0}$ and $\mu_{0}$ are two universal constants, called the permittivity of free space and permeability of free space, respectively. With relations (8), Eqs.(7) can be rewritten into an explicit iterative scheme for TE wave.

[TABLE]

The symmetric TE wave scheme induced from Eqs.(6) can be written as follows.

[TABLE]

Schemes for TM wave

If writing

[TABLE]

where $H=\sum\limits_{E}{H}_{i}e^{i}$ is the discrete $1-$form on space, $D=\sum\limits_{P}{H}_{i}P^{i}$ is the discrete $2-$form on space, $E=\sum\limits_{*P}{D}_{i}*P^{i}$ is the dual of $D$ on space, $B=\sum\limits_{*E}{B}_{i}*e^{i}$ is the dual of $H$ on space, $\rho_{e}=\sum\limits_{Tet}\rho_{ei}T^{i}$ is the discrete $3-$form on space, ${J}_{e}=\sum\limits_{P}J_{ei}P^{i}$ is the discrete $2-$form on space, $\rho_{m}dt$ is the discrete $1-$form on time, ${J}_{m}=\sum\limits_{E}J_{mi}E^{i}$ is the discrete $1-$form on space, the discrete Maxwell’s equations can be rewritten as

[TABLE]

Proposition 3.2

If the initial condition satisfies the first and third equations in Eqs.(11), the solution of the second and fourth equations in Eqs.(11) automatically satisfy Eqs.(11).

Proof.Because the dimension of spacetime is $3+1$ or $2+1$, therefore

[TABLE]

and the continuity equations can be reduced to

[TABLE]

So we have

[TABLE]

$\Box$

Now we show the scheme (11) on the product of 2D discrete space manifold and time. The second and fourth equations in Eqs.(11) based on Fig.1 are

[TABLE]

With relations (8), Eqs.(12) can be rewritten into an explicit iterative scheme for TM wave.

[TABLE]

General schemes

For real world materials, the constitutive relations are not simple proportionalities, except approximately. The relations can usually still be written:

[TABLE]

but $\varepsilon$ and $\mu$ are not, in general, simple constants, but rather functions. With Ohm’s law

[TABLE]

where $\sigma$ is the electrical conductivity and $\sigma_{m}$ is magnetic conductivity. The DEC schemes can be written as

[TABLE]

[TABLE]

4 Stability, convergence and accuracy

Stability

The Courant-Friedrichs-Lewy condition is a necessary condition for convergence while solving certain partial differential equations numerically. Now, we find this condition for scheme (10). Condition for scheme (13) can be induced in the same way. First, we decompose DEC algorithm into temporal and spacial eigenvalue problems.

The temporal eigenvalue problem:

[TABLE]

It can approximated by difference equation

[TABLE]

Supposing

[TABLE]

and

[TABLE]

and substituting those into Eq.(14), we obtain

[TABLE]

[TABLE]

therefore

[TABLE]

This is the stabile condition for the temporal eigenvalue problem.

The spacial eigenvalue problem:

[TABLE]

It can be approximated by difference equation (15) based on Fig.2.

[TABLE]

[TABLE]

Let

[TABLE]

and substitute into Eq.(15) to obtain

[TABLE]

[TABLE]

So we have

[TABLE]

In order to keep the stability of scheme (13), we need

[TABLE]

or

[TABLE]

Convergence

By the definition of truncation error, the exact solution of Maxwell’s equations satisfy the same relation as DEC scheme except for an additional term $O((\Delta t)^{2}+\Delta t|\ast e|)$. This expresses the consistency, and so convergence for DEC scheme by Lax equivalence theorem (consistency $+$ stability $=$ convergence).

Accuracy

The derivative of Maxwell’s equations is approximated by first order difference in schemes (10) and (13). Equivalently, $H$ and $E$ are approximated by linear interpolation functions. Consulting the definition about accuracy of finite volume method, we can also say that schemes (10) and (13) have first order spacial and temporal accuracy, and have second order spacial and temporal accuracy on rectangular grid with equivalent space and time steps.

5 Implementation

The DEC algorithm of Maxwell’s equations was implemented in C++ platform. The Fig.3 shows the flowchart of DEC schemes for Maxwell’s equations.

[TABLE]

In the common practice, not every simulation step needs to be visualized, especially when the time step size is too small. Fig.4 exhibit Gaussian pluses’ waveforms simulated by DEC.

[TABLE]

Fig.5 exhibits two sources Gaussian pluses interference simulated by DEC. Our algorithm can simulate more complex situation on surface and 3D-space manifold.

[TABLE]