# Generalized Helmholtz Decomposition for Modal Analysis of Electromagnetic
Problems in Inhomogeneous Media

Jie Zhu, Thomas E. Roth,
Dong-Yeop Na, and Weng Cho Chew

###### Abstract

Potential-based formulation with generalized Lorenz gauge can be used in the
quantization of electromagnetic fields in inhomogeneous media [1]. However,
one often faces the redundancy of modes when finding eigenmodes from
potential-based formulation. In free space, this can be explained by the
connection to the well-known Helmholtz decomposition. In this work, we
generalize the Helmholtz decomposition to its generalized form, echoing the
use of generalized Lorenz gauge in inhomogeneous media [2]. We formulate
electromagnetics eigenvalue problems using vector potential formulation which
is often used in numerical quantization [1]. The properties of the
differential operators are mathematically analyzed. Orthogonality relations
between the two classes of modes are proved in both continuous and discrete
space. Completeness of two sets of modes and the orthogonality relations are
numerically validated in inhomogeneous anisotropic media. This work serves as
a foundation for numerical quantization of electromagnetic fields in
inhomogeneous media with potential-based formulation.

##  I Introduction

The formulation of electromagnetic theory based on E, H, D, and B offers
physical insight that has resulted in the development of many electromagnetic-
related technologies. However, classical electromagnetic equations formulated
in terms of E-H have low-frequency or long-wavelength breakdown [2]. Many
numerical methods based on E-H formulation are unstable at low frequencies.
Therefore, E-H formulation is not truly multi-scale, but exhibits catastrophic
breakdown when the dimensions of objects become much smaller than the local
wavelength, or when the frequency is low.

Much previous work has been devoted to eliminating low-frequency breakdown of
E-H formulation. For differential equation solvers, Manges and Cendes [3]
proposed a generalized tree-cotree gauge to eliminate the matrix null space of
the magnetostatic curl-curl equation. In another work [4], hierarchical basis
functions and tree-cotree splitting are used to enhance the stability of the
finite-element procedure. For time domain solvers, method based on changing
curl-curl operator to Laplacian is proposed to handle the ill-conditioning of
the system matrix [5]. For integral equation solvers, loop-tree decomposition
was first proposed which separates the electrostatic and magnetostatic physics
[6, 7, 8, 9, 10, 11]. Augmented electric field integral equation was then
proposed which eliminated the need to search for loops [12].

Electromagnetic formulation based on potential theory has received increasing
attention recently since it eliminates the low-frequency breakdown issue
elegantly. It removes the null space solution by converting singular matrices
to full-rank Laplacian-like matrices. In a homogeneous medium, it is easy to
find its connection to the Helmholtz decomposition. A number of previous works
have dealt with potential-based (or $\mathbf{A}$-$\Phi$) formulation with both
differential equations [13, 14, 15, 16], and integral equations [17, 18, 19,
20, 21, 22, 23, 24]. Moreover, potential-based formulation is more compatible
with quantum theory. Aharonov-Bohm effect is an example in quantum theory
where E-H are zero, but vector potential A is not zero, and yet, the effect of
A is felt [25]. Besides, in many quantum optics studies, both vector potential
$\mathbf{A}$ and scalar potential $\Phi$ are used when the electromagnetic
fields are incorporated into the Schrödinger equation [15].

To perform numerical quantization of electromagnetic fields in the mode
picture with potential-based formulation, it is often needed to find the
electromagnetic eigenmodes of the system with inhomogeneous media by numerical
methods [1, 26, 27]. When solving potential-based formulation as
electromagnetics eigenvalue problems, null space solution does not exist, but
one often obtains some “extra solutions” that feature zero $\mathbf{E}$ and
$\mathbf{H}$ even though $\mathbf{A}$ and $\Phi$ are nonzero. In a homogeneous
medium, this can be explained by making connection to the Helmholtz
decomposition, as will be shown in Section II.

Due to the need of numerical quantization in complex media in modern quantum
technology [1, 26], an understanding of the solution space of vector potential
wave equation in general inhomogeneous media is called for. The connection
between the extra modes and the solution space should be made. A general
scheme is needed to explain it in inhomogeneous media.

In this work, we investigate how to generalize Helmholtz decomposition to
inhomogeneous media. We formulate electromagnetics eigenvalue problems using
the vector potential formulation, and investigate the modes from solving the
system equations. The formulation can be readily used in the numerical
quantization scheme in quantum electromagnetics [1, 26, 27]. We first review
Helmholtz theorem and motivate the need for its generalization. We then begin
the proof by analyzing the differential operator associated with the wave
equation for vector potential in Section II. In Section III, we define
div-$\varepsilon$-free modes in inhomogeneous media, which is the counterpart
of div-free modes in homogeneous case. We prove analytically that the
div-$\varepsilon$-free modes are orthogonal to the curl-free modes in source-
less inhomogeneous regions, by both continuous-space and discrete-space
calculus. No such proof has been found in the literature. The proof can be
readily generalized to anisotropic media. We provide numerical results in
Section IV: We first show numerically that Helmholtz decomposition can be
extended to inhomogeneous case, which we call generalized Helmholtz
decomposition (GHD). We then show numerically the orthogonality condition. We
conclude the work in Section V.

##  II A Review of Helmholtz Theorem and Motivation of Generalization

###  II-A A Traditional Proof

Two important theorems establish the conditions for the existence and
uniqueness of solutions to time-independent problems in electromagnetic
theory. The first one is [28]:

> A vector field is uniquely specified by giving its divergence and its curl
> within a simply connected region and its normal component on the boundary.

It can be proved by contradiction. Assume that two different vector fields
$\mathbf{A}$ and $\mathbf{A^{\prime}}$ have the same divergence and curl

$\displaystyle\nabla\cdot\mathbf{A}$
$\displaystyle=\nabla\cdot\mathbf{A^{\prime}}=s$ (1)
$\displaystyle\nabla\times\mathbf{A}$
$\displaystyle=\nabla\times\mathbf{A^{\prime}}=\mathbf{c}$ (2)

and their normal component $\mathbf{A}_{n}$ and $\mathbf{A}^{\prime}_{n}$ on
the boundary are also the same. Let
$\mathbf{F}=\mathbf{A}-\mathbf{A^{\prime}}$. Then $\nabla\cdot\mathbf{F}$,
$\nabla\times\mathbf{F}$, and $\mathbf{F}_{n}$ are all zero. Thus there exists
a $\phi$ such that $\mathbf{F}=-\nabla\phi$ and it satisfies

$\displaystyle\nabla^{2}\phi=0.$ (3)

Due to Earnshaw’s theorem, the solution to (3) must have its maximum and
minimum on the boundary. Together with the condition $\mathbf{F}_{n}=0$, or
$\dfrac{\partial\phi}{\partial n}=0$, it must have that $\phi$ is a constant
everywhere. Hence $\mathbf{F}=0$ everywhere. Thus $\mathbf{A}$ and
$\mathbf{A^{\prime}}$ are identical.

Given the proof of uniqueness, we then review the Helmholtz theorem. The
Helmholtz theorem states [28, 29, 30]:

> A vector $\mathbf{A}$ with both source and circulation densities vanishing
> at infinity may be written as the sum of two parts, one of which is
> irrotational, curl-free, or longitudinal, the other of which is solenoidal,
> div-free, or transverse.

That is,

$\displaystyle\mathbf{A}=\mathbf{A}_{\parallel}+\mathbf{A}_{\perp}=-\nabla\phi+\nabla\times\mathbf{P}$
(4)

where $\mathbf{A}_{\parallel}=-\nabla\phi$ and
$\mathbf{A}_{\perp}=\nabla\times\mathbf{P}$. Here we use
$\mathbf{A}_{\parallel}$ and $\mathbf{A}_{\perp}$ to denote longitudinal
fields and transverse fields, respectively.111It should be noted that the term
longitudinal and transverse are first used in the context of plane wave, where
$\nabla\times\mathbf{A}=0\rightarrow\mathbf{k}\times\mathbf{A}=0$ and
$\nabla\cdot\mathbf{A}=0\rightarrow\mathbf{k}\cdot\mathbf{A}=0$, suggesting
that the longitudinal (transverse) field is parallel (perpendicular) to the
propagation direction in the Fourier space. Here we port the use of these
terms to curl-free and divergence free in general cases, with a caveat that
they do not have the same physical meaning as in plane wave. Apparently
$\nabla\times\mathbf{A}_{\parallel}=0$ and $\nabla\cdot\mathbf{A}_{\perp}=0$.
Since a vector can be uniquely defined by its divergence
$\nabla\cdot\mathbf{A}=s$ and curl $\nabla\times\mathbf{A}=\mathbf{c}$, taking
the divergence and curl of the above, we have

$\displaystyle\nabla\cdot\mathbf{A}$ $\displaystyle=-\nabla^{2}\phi=s$ (5)
$\displaystyle\nabla\times\mathbf{A}$
$\displaystyle=\nabla\times\nabla\times\mathbf{P}=\mathbf{c}$ (6)

where $s$ acts as the source for Poisson equation, and $\mathbf{c}$ the source
for $\mathbf{P}$. We proceed the proof of Helmholtz theorem by showing the
expression for $\phi$ and $\mathbf{A}$ that recover $s$ and $\mathbf{c}$.
Since the solution of $\mathbf{P}$ in (6) is not unique, we set
$\nabla\cdot\mathbf{P}=0$, and (6) becomes vector Poisson equation with a
unique solution. The solution to (5) and (6) are then well known as [30]

$\displaystyle\phi(\mathbf{r})$
$\displaystyle=\frac{1}{4\pi}\int\frac{s(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}d\mathbf{r}^{\prime}$
(7) $\displaystyle\mathbf{P}(\mathbf{r})$
$\displaystyle=\frac{1}{4\pi}\int\frac{\mathbf{c}(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}d\mathbf{r}^{\prime}.$
(8)

The above can be proved by direct substitution of (7) and (8) into (5) and
(6).

###  II-B Solution Space of Wave Equations

Helmholtz theorem can be used to explain the null-space solution of electric
field wave equation. We consider the electric field wave equations in a
homogeneous medium without sources

$\displaystyle\nabla\times\nabla\times\mathbf{E}=\omega^{2}\varepsilon\mu\mathbf{E}.$
(9)

It can be formulated as an eigenvalue problem. According to the Helmholtz
theorem, the solution to (9) must be either div-free, or curl-free. Apparently
curl-free solutions are the null-space solutions to (9), where the eigenvalues
are zero. On the other hand, div-free solutions are the non-null-space
solutions. This is because when $\nabla\cdot\mathbf{E}=0$, we can add
$-\nabla\nabla\cdot\mathbf{E}$ to the left-hand side of (9) without changing
the equation. Then the left-hand side becomes
$\nabla\times\nabla\times\mathbf{E}-\nabla\nabla\cdot\mathbf{E}=-\nabla^{2}\mathbf{E}$
and the Laplacian $\nabla^{2}$ is a negative definite operator. Thus the div-
free solutions must have nonzero eigenvalues.

Helmholtz theorem also connects the solution to the vector potential wave
equation with the solution to the electric field wave equation. We consider
the vector potential wave equations with Lorenz gauge
$\nabla\cdot\mathbf{A}=-\varepsilon\mu\partial_{t}\Phi$ in the same
homogeneous medium without sources

$\displaystyle\nabla\times\nabla\times\mathbf{A}-\nabla\nabla\cdot\mathbf{A}=\omega^{2}\varepsilon\mu\mathbf{A}.$
(10)

Both (9) and (10) have two families of solutions following the Helmholtz
theorem. The div-free solution to (9) must also be the div-free solution to
(10), and vice versa. This is because the second term of (10) can be dropped
in the case of div-free solution. However, since null-space solution does not
exist in (10) due to the negative-definiteness of Laplacian, the space spanned
by the null-space solutions to (9) must be filled by the same number of non-
null-space solutions to (10). Again due to the Helmholtz theorem, these
solutions are curl-free solutions in (10) as well. The mapping between the
solution space of (9) and (10) is illustrated in Fig. 1. Apparently curl-free
vector potential in electrodynamics leads to zero $\mathbf{E}$ and
$\mathbf{H}$ even though $\mathbf{A}$ and $\Phi$ are nonzero. This also
justifies the use of $\Phi=0$ gauge due to Lorenz gauge to remove the
redundancy of the curl-free solutions in homogeneous case for many practical
problems.

Figure 1: Mapping between solution space of the electric field wave equation
and vector potential wave equation. The double-headed arrow represents the
same space spanned by the eigenvectors of both equations.

###  II-C A New Look of Helmholtz Theorem

We can also prove Helmholtz theorem from a different perspective by looking at
(9). We first prove that the curl operator is a self-adjoint operator [31]
using inner product. We begin with

$\displaystyle\langle\mathbf{A}^{*},\overline{\bm{\mathcal{C}}}\mathbf{B}\rangle$
$\displaystyle=\int_{V}d\mathbf{r}\,\mathbf{A}^{*}\cdot\left(\nabla\times\mathbf{B}\right)$
$\displaystyle=\int_{V}d\mathbf{r}\,\mathbf{B}\cdot\left(\nabla\times\mathbf{A}^{*}\right)-\int_{S}d\mathbf{r}\,\hat{n}\cdot(\mathbf{A}^{*}\times\mathbf{B})$
(11)

where $\overline{\bm{\mathcal{C}}}$ denotes the curl operator. With the
boundary condition that $\hat{n}\times\mathbf{A}=0$ and
$\hat{n}\times\mathbf{B}=0$,222The rationale of this boundary condition is
explained in the next section. and the surface integral term vanishes, and
hence the curl operator is self-adjoint for fields of the prescribed boundary
conditions. With the same approach, it can also be proved that
$\nabla\times\nabla\times$ is a self-adjoint operator [16]. Thus the solution
vectors to (9) form a complete basis.

The proceeding proof also suggests that the curl operator has the same null-
space as its power. In general when an operator maps a domain space function
(or a row space vector, if the operator is a finite-sized matrix) into its
range space (or column space), the new function may contain null-space
component. However, the self-adjointness and symmetry of the curl operators
ensure that its domain space and range space are the same, and its
(right-)null-space and left-null-space are the same. Thus the curl operator
maps a domain space function into the same space. Hence it will not generate a
null-space component.

Now we are ready to prove the Helmholtz theorem by analyzing (9). We only
consider the non-trivial solution to (9) where $\mathbf{E}\neq 0$ almost
everywhere. Taking the divergence of (9), we have
$\omega^{2}\nabla\cdot\mathbf{E}=0$. It has three possible classes of
solutions: (i) $\omega=0$ and $\nabla\cdot\mathbf{E}\neq 0$; (ii) $\omega\neq
0$ and $\nabla\cdot\mathbf{E}=0$; and (iii) $\omega=0$ and
$\nabla\cdot\mathbf{E}=0$. We can rule out (iii) because it would make (9) the
same as (10), and $\omega=0$ solution cannot exist in (10). This is because
$\nabla\times\nabla\times\mathbf{A}-\nabla\nabla\cdot\mathbf{A}=-\nabla^{2}\mathbf{A}$
and the Laplacian $\nabla^{2}$ is a negative definite operator. For (i),
apparently $\nabla\times\nabla\times\mathbf{E}=0$, and in turn
$\nabla\times\mathbf{E}=0$. For (ii), we have $\nabla\times\mathbf{E}\neq 0$.
Concluding the above, we can clearly see two classes of solution to (9)
similar to the Helmholtz theorem, where

  1. 1.

$\nabla\times\mathbf{E}=0$, $\nabla\cdot\mathbf{E}\neq 0$;

  2. 2.

$\nabla\cdot\mathbf{E}=0$, $\nabla\times\mathbf{E}\neq 0$.

Due to the self-adjointness of $\nabla\times\nabla\times$ operator, the
solution to (9) form a complete basis to the entire space. Hence any function
satisfying the proposed boundary condition can be separated into two classes.
This completes the proof of Helmholtz theorem. Note that analyzing (10) also
leads to the same two families of solutions described above. This will be
shown in greater details similar to the inhomogeneous media case next.

###  II-D Motivation of Generalization

As shown above, the Helmholtz theorem works well in explaining the solution
space of the electric field wave equation and vector potential wave equation
with Lorenz gauge in homogeneous medium. However, in inhomogeneous
dispersionless media, it is suggested that generalized Lorenz gauge should be
used instead [2]. Thus, to better explain the eigenmodes found by potential-
based formulation, it would be ideal if Helmholtz theorem can be generalized
in a way that
$\nabla\cdot\varepsilon(\mathbf{r})\mathbf{A}_{\varepsilon\perp}(\mathbf{r})=0$,333Position-
dependent $\varepsilon(\mathbf{r})$ is assumed from now on. For simplicity, we
only write $\varepsilon$.  where $\mathbf{A}_{\varepsilon\perp}$ denotes the
div-$\varepsilon$-free (or generalized transverse) fields. Considering this,
we seek to prove the following generalized Helmholtz decomposition:

> The vector field $\mathbf{A}$ given by the Helmholtz wave equation in
> source-free inhomogeneous media can be decomposed into two components
> $\mathbf{A}_{\varepsilon\perp}$ and $\mathbf{A}_{\parallel}$, where the
> div-$\varepsilon$-free (or generalized transverse) fields satisfy
> $\nabla\cdot\varepsilon\mathbf{A}_{\varepsilon\perp}=0$ and the curl-free
> (or longitudinal) fields satisfy $\nabla\times\mathbf{A}_{\parallel}=0$.

In the remainder of this work, we use solution of the wave equation for vector
potential $\mathbf{A}$ to arrive at the generalized Helmholtz decomposition. A
modified generalized Lorenz gauge
$\nabla\cdot\varepsilon\mathbf{A}=-\varepsilon_{0}^{2}\mu_{0}\partial_{t}\Phi$
is used. We follow matrix theory generally in the proof. It will be shown that
the vector potential wave equation has two classes of solutions. They
degenerate into the Helmholtz decomposition in free space. Thus the vector
potential wave equation for inhomogeneous media can be used to motivate the
generalized Helmholtz decomposition. We also note that there is no null space
in the resulting generalized eigenvalue problem.

##  III Proof of Orthogonality

We start from proving the orthogonality between the two families of fields
proposed in the generalized Helmholtz decomposition, assuming the
decomposition does exist. We leave it to the next section to demonstrate the
existence of the decomposition through numerical examples.

We consider the space terminated by continuous perfect electric conductor
(PEC) boundary, filled with inhomogeneous media (Fig. 2). Because of the gauge
$\nabla\cdot\varepsilon\mathbf{A}=-\varepsilon_{0}^{2}\mu_{0}\partial_{t}\Phi$,
the boundary condition for $\nabla\cdot\varepsilon\mathbf{A}$ is actually that
for $\Phi$ on a PEC, which is $\Phi$ equals a constant voltage [18, 16]. If
there is only one PEC, we can choose this as our reference. Thus

$\displaystyle\hat{n}\times\mathbf{A}=0$ (12)

which is consistent with $\hat{n}\times\mathbf{E}=0$, and

$\displaystyle\nabla\cdot\varepsilon\mathbf{A}=0.$ (13)

The latter is equivalent to $\Phi=0$ on the boundary when the PEC has a
constant voltage of zero, i.e. “grounded”.

If there are additional PECs inside the domain, they may not be grounded. In
this case, $\nabla\cdot\varepsilon\mathbf{A}=\Phi$ is equal to a constant
voltage on each conductor. Additionally, charge neutrality constraint can be
enforced on each conductor, which requires [16, 18]

$\displaystyle\int_{S}d\mathbf{r}\,\hat{n}\cdot\varepsilon\mathbf{A}\;=0.$
(14)

where the surface integral is over the surface of each PEC.

Figure 2: Schematic of the region of interest.

Next, we proceed to show the proof of orthogonality first in continuous space,
followed by a similar procedure in discrete space. Different definitions of
orthogonality will be discussed. It will also be shown that discrete
approximation of Maxwell’s equations is homomorphic to the continuum case if
done properly.

###  III-A Continuous Space

Wave equation is usually the starting point of many $\mathbf{A}$-$\Phi$ based
solver. We start the mathematical proof of mode orthogonality by defining
proper operator and vector space. We consider the wave equation for vector
potential $\mathbf{A}$ in lossy dispersionless inhomogeneous isotropic
media444To satisfy the Kramers–Kronig relations, this assumption is only valid
over a narrow bandwidth. with modified generalized Lorenz gauge [2]

$\displaystyle\nabla\times\frac{1}{\mu}\nabla\times\mathbf{A}-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\varepsilon\nabla\nabla\cdot\varepsilon\mathbf{A}=\omega^{2}\varepsilon\mathbf{A}.$
(15)

Let

$\displaystyle\overline{\bm{\mathcal{W}}}_{1}$
$\displaystyle=\nabla\times\frac{1}{\mu}\nabla\times$ (16)
$\displaystyle\overline{\bm{\mathcal{W}}}_{2}$
$\displaystyle=-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\varepsilon\nabla\nabla\cdot\varepsilon$
(17) $\displaystyle\overline{\bm{\mathcal{L}}}$
$\displaystyle=\overline{\bm{\mathcal{W}}}_{1}+\overline{\bm{\mathcal{W}}}_{2}$
(18) $\displaystyle\overline{\bm{\mathcal{U}}}$ $\displaystyle=\varepsilon$
(19) $\displaystyle\lambda$ $\displaystyle=\omega^{2}.$ (20)

Then (15) can be written abstractly and concisely as

$\displaystyle(\overline{\bm{\mathcal{W}}}_{1}+\overline{\bm{\mathcal{W}}}_{2})\mathbf{A}=\lambda\bm{\overline{\mathcal{U}}}\mathbf{A}$
(21)

where the operator notation is defined as a differential operator acting on a
function $\mathbf{A}$. Note that we use boldface calligraphic font with an
overbar to denote a Hilbert space operator. The function $\mathbf{A}$ in the
abovementioned space is in turn square integrable [32].

It has been proved that the $\overline{\bm{\mathcal{L}}}$ operator is self-
adjoint in lossless media given the boundary conditions defined in the
beginning of this section [16]. Here we investigate the general case with
lossy media. Since the $\overline{\bm{\mathcal{L}}}$ operator consists of two
terms $\overline{\bm{\mathcal{W}}}_{1}$ and $\overline{\bm{\mathcal{W}}}_{2}$,
it is clearer to look at the two terms separately [16].

We first define $\mathbf{A}_{\varepsilon\perp}$ to be the solution to the
eigenvalue problem
$\overline{\bm{\mathcal{W}}}_{1}\mathbf{A}_{\varepsilon\perp}=\lambda_{\varepsilon\perp}\overline{\bm{\mathcal{U}}}\mathbf{A}_{\varepsilon\perp}$
with eigenvalue $\lambda_{\varepsilon\perp}$. Here we only study the solutions
with nonzero eigenvalues. This gives us, from the definition above,

$\displaystyle\nabla\times\frac{1}{\mu}\nabla\times\mathbf{A}_{\varepsilon\perp}=\lambda_{\varepsilon\perp}\varepsilon\mathbf{A}_{\varepsilon\perp}.$
(22)

Taking divergence of the above equation gives

$\displaystyle
0=\lambda_{\varepsilon\perp}\nabla\cdot\varepsilon\mathbf{A}_{\varepsilon\perp}.$
(23)

Thus $\mathbf{A}_{\varepsilon\perp}$ is in the null space of the
$\overline{\bm{\mathcal{W}}}_{2}$ operator defined in (17), and it is also a
solution to the original eigenvalue problem (21) with eigenvalue
$\lambda_{\varepsilon\perp}$.

Similarly, we define $\mathbf{A}_{\parallel}$ to be the solution to
$\overline{\bm{\mathcal{W}}}_{2}\mathbf{A}_{\parallel}=\lambda_{\parallel}\overline{\bm{\mathcal{U}}}\mathbf{A}_{\parallel}$
with nonzero eigenvalue $\lambda_{\parallel}$, or more explicitly,

$\displaystyle-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\varepsilon\nabla\nabla\cdot\varepsilon\mathbf{A}_{\parallel}=\lambda_{\parallel}\varepsilon\mathbf{A}_{\parallel}.$
(24)

Dividing the above by $\varepsilon$ and taking curl of the resulting equation
yields

$\displaystyle 0=\lambda_{\parallel}\nabla\times\mathbf{A}_{\parallel}.$ (25)

Thus $\mathbf{A}_{\parallel}$ is in the null space of the
$\overline{\bm{\mathcal{W}}}_{1}$ operator defined in (16), and it is also a
solution to the original eigenvalue problem (21) with eigenvalue
$\lambda_{\parallel}$.

Next, we will manipulate the above equations, and use (23) and (25) to prove
orthogonality of two sets of modes.

####  III-A1 Proof 1

To begin, we take the left dot product of (24) by
$\mathbf{A}_{\varepsilon\perp}$, and integrate it over the space filled with
dielectrics to yield

$\displaystyle-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\int_{V}d\mathbf{r}\left(\mathbf{A}_{\varepsilon\perp}\cdot\varepsilon\nabla\nabla\cdot\varepsilon\mathbf{A}_{\parallel}\right)=\int_{V}d\mathbf{r}\left(\mathbf{A}_{\varepsilon\perp}\cdot\lambda_{\parallel}\varepsilon\mathbf{A}_{\parallel}\right).$
(26)

Using integration by parts, the left-hand side becomes

$\displaystyle\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\int_{V}d\mathbf{r}\,(\nabla\cdot\varepsilon\mathbf{A}_{\varepsilon\perp})(\nabla\cdot\varepsilon\mathbf{A}_{\parallel})$
$\displaystyle-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\int_{S}d\mathbf{r}\,\hat{n}\cdot$
$\displaystyle\left[(\nabla\cdot\varepsilon\mathbf{A}_{\parallel})\varepsilon\mathbf{A}_{\varepsilon\perp}\right]$
(27)

where the surface integral is over the PEC surfaces. From (23), the first term
vanishes. At the grounded PEC boundary,
$\nabla\cdot\varepsilon\mathbf{A}_{\parallel}=0$, and then the second term
also vanishes. At the PEC not connected to the ground, the second term becomes

$\displaystyle-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}(\nabla\cdot\varepsilon\mathbf{A}_{\parallel})\int_{S}d\mathbf{r}\,\hat{n}\cdot\varepsilon\mathbf{A}_{\varepsilon\perp}=0$
(28)

where the divergence term is taken outside the integral because $\Phi$ is
equal to a constant voltage on the conductor, and the integral vanishes
because of the charge neutrality constraint (14). Thus

$\displaystyle\langle\mathbf{A}_{\varepsilon\perp},\varepsilon\mathbf{A}_{\parallel}\rangle=0$
(29)

as long as $\lambda_{\parallel}\neq 0$, where
$\langle\mathbf{F},\mathbf{G}\rangle=\int\mathbf{F}\cdot\mathbf{G}\;d\mathbf{r}$
is defined as the reaction inner product by Rumsey [33]. In other words, the
two sets of modes $\mathbf{A}_{\varepsilon\perp}$ and $\mathbf{A}_{\parallel}$
are $\varepsilon$-orthogonal to each other.

####  III-A2 Proof 2

Now, taking the left dot product of (22) with $\mathbf{A}_{\parallel}$, and
integrating it over space yields

$\displaystyle\int_{V}d\mathbf{r}\left(\mathbf{A}_{\parallel}\cdot\nabla\times\frac{1}{\mu}\nabla\times\mathbf{A}_{\varepsilon\perp}\right)=\int_{V}d\mathbf{r}\left(\mathbf{A}_{\parallel}\cdot\lambda_{\varepsilon\perp}\varepsilon\mathbf{A}_{\varepsilon\perp}\right).$
(30)

The right-hand side is denoted
$\lambda_{\varepsilon\perp}\langle\mathbf{A}_{\parallel},\varepsilon\mathbf{A}_{\varepsilon\perp}\rangle$.
Using integration by parts, the left-hand side becomes

$\displaystyle\int_{V}d\mathbf{r}\,\left(\frac{1}{\mu}\nabla\times\mathbf{A}_{\varepsilon\perp}\right)\cdot\left(\nabla\times\mathbf{A}_{\parallel}\right)$
$\displaystyle+\int_{S}d\mathbf{r}\,\hat{n}\cdot$
$\displaystyle\left[\left(\frac{1}{\mu}\nabla\times\mathbf{A}_{\varepsilon\perp}\right)\times\mathbf{A}_{\parallel}\right].$
(31)

From (25), the first term vanishes. Implementing the boundary condition that
$\hat{n}\times\mathbf{A}_{\parallel}=0$, the second term also vanishes. Thus

$\displaystyle\langle\mathbf{A}_{\parallel},\varepsilon\mathbf{A}_{\varepsilon\perp}\rangle=0$
(32)

as long as $\lambda_{\varepsilon\perp}\neq 0$. In other words, the two sets of
modes $\mathbf{A}_{\parallel}$ and $\mathbf{A}_{\varepsilon\perp}$ are
$\varepsilon$-orthogonal to each other.

In the above derivation, if we replace $\mathbf{A}_{\parallel}$ by
$\mathbf{A}_{\parallel}^{*}$, the derivation is still valid. That is

$\displaystyle\langle\mathbf{A}_{\parallel}^{*},\varepsilon\mathbf{A}_{\varepsilon\perp}\rangle=0.$
(33)

However, the same token does not apply to Proof 1 in III-A. This may motivate
new definition of inner product in the problems concerning potential-based
formulation.

###  III-B Discrete Space

The problem in the previous subsection is homomorphic to numerical linear
algebra if we find the matrix representation of the operators using subspace
projection methods [34], including finite difference method [35], finite
element method [36], or discrete exterior calculus [37]. These methods should
be chosen properly to preserve certain properties of the continuum calculus,
i.e. $\nabla\cdot(\nabla\times\mathbf{A})=0$ and $\nabla\times\nabla f=0$. In
this work we use finite difference method. After discretization, we have

$\displaystyle\nabla\times\mathbf{A}$
$\displaystyle\Rightarrow\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{1}$ (34)
$\displaystyle\nabla\times\nabla\times\mathbf{A}$
$\displaystyle\Rightarrow\overline{\mathbf{C}}_{2}\cdot\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{1}$
(35) $\displaystyle\nabla\cdot\varepsilon\mathbf{A}$
$\displaystyle\Rightarrow\overline{\mathbf{D}}_{1}\cdot\mathsf{A}_{1}$ (36)
$\displaystyle\varepsilon\nabla f$
$\displaystyle\Rightarrow\overline{\mathbf{E}}_{1}\cdot\mathsf{f}_{1}$ (37)
$\displaystyle\varepsilon$ $\displaystyle\Rightarrow\overline{\mathbf{U}}$
(38)

where the operators on the LHS are in the continuous Hilbert space, while the
RHS are the linear algebra approximation of the continuum space. Here we use
boldface font with overbar to denote finite size matrices, and sans serif
typestyle to denote finite length vectors stored as 1D arrays in computers.

Due to the discretization, each operator may have different representations
depending on the space it operates on. To use finite difference method, Yee
grid is constructed with both primal and dual grids [38, 35]. On each grid, we
place a 3D vector on the face of the cube, and a scalar on the center of the
cube. Details of the discretization scheme can be found in the Appendix A.

In (34) to (38), the discrete operators with subscript $1$ denotes operations
on the primal mesh, while subscript $2$ denotes operations on the dual mesh.
The discrete curl operator $\overline{\mathbf{C}}_{1}$ acts on a vector on the
primal mesh, resulting in a vector on the dual mesh, while
$\overline{\mathbf{C}}_{2}$ acts on a vector on the dual mesh and results in a
vector on the primal mesh. It can be shown that these two discrete curl
operators are transpose to each other
$\overline{\mathbf{C}}_{1}=\overline{\mathbf{C}}_{2}^{T}$, where T denotes
matrix transpose. Also, it is easy to prove that the discrete divergence
operator and the negative of discrete gradient operator are transpose to each
other $\overline{\mathbf{D}}_{1,2}^{T}=-\overline{\mathbf{E}}_{1,2}$. The
$\overline{\mathbf{U}}$ operator is a square matrix whose dimension depends on
the vector it operates on. The following vector identities are preserved after
discretization using finite difference method [39, 40, 41]

$\displaystyle\nabla\cdot(\nabla\times\mathbf{A})=0$
$\displaystyle\Rightarrow\overline{\mathbf{D}}_{2}\cdot\overline{\mathbf{U}}^{-1}_{2}\cdot\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{1}=0$
$\displaystyle\mathrm{and}\;\overline{\mathbf{D}}_{1}\cdot\overline{\mathbf{U}}^{-1}_{1}\cdot\overline{\mathbf{C}}_{2}\cdot\mathsf{A}_{2}=0$
(39) $\displaystyle\nabla\times\nabla f=0$
$\displaystyle\Rightarrow\overline{\mathbf{C}}_{1}\cdot\overline{\mathbf{U}}^{-1}_{1}\cdot\overline{\mathbf{E}}_{1}\cdot\mathsf{f}_{1}$
$\displaystyle=-\overline{\mathbf{C}}_{1}\cdot\overline{\mathbf{U}}^{-1}_{1}\cdot\overline{\mathbf{D}}_{1}^{T}\cdot\mathsf{f}_{1}=0.$
(40)

Without loss of generality, we let $\mu=1$.555It can be seen that the proof is
still valid when $\mu$ is inhomogeneous, anisotropic, or complex. Also we
denote $\alpha=1/\varepsilon_{0}^{2}\mu_{0}$. We assume the solution vector is
on the primal mesh, and thus the subscript $1$ is omitted. Then (21) becomes

$\displaystyle\overline{\mathbf{C}}_{2}\cdot\overline{\mathbf{C}}_{1}\cdot\mathsf{A}-\alpha\overline{\mathbf{E}}_{1}\cdot\overline{\mathbf{D}}_{1}\cdot\mathsf{A}=\lambda\overline{\mathbf{U}}_{1}\cdot\mathsf{A}.$
(41)

We break the above into two eigenvalue problems corresponding to the two
operators on the left-hand side

$\displaystyle\overline{\mathbf{C}}_{2}\cdot\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}$
$\displaystyle=\lambda_{\varepsilon\perp}\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}$
(42)
$\displaystyle-\alpha\overline{\mathbf{E}}_{1}\cdot\overline{\mathbf{D}}_{1}\cdot\mathsf{A}_{\parallel}$
$\displaystyle=\lambda_{\parallel}\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\parallel}.$
(43)

Taking the left dot product of (42) with
$\overline{\mathbf{D}}_{1}\cdot\overline{\mathbf{U}}_{1}^{-1}$, we get

$\displaystyle\overline{\mathbf{D}}_{1}\cdot\overline{\mathbf{U}}_{1}^{-1}\cdot\overline{\mathbf{C}}_{2}\cdot\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}$
$\displaystyle=\lambda_{\varepsilon\perp}\overline{\mathbf{D}}_{1}\cdot\overline{\mathbf{U}}_{1}^{-1}\cdot\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}.$
(44)

Due to (39), LHS of the above vanishes. Thus

$\displaystyle\overline{\mathbf{D}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}=0.$
(45)

Taking the left dot product of (43) with
$\overline{\mathbf{C}}_{1}\cdot\overline{\mathbf{U}}_{1}^{-1}$, we get

$\displaystyle-\alpha\overline{\mathbf{C}}_{1}\cdot\overline{\mathbf{U}}_{1}^{-1}\cdot\overline{\mathbf{E}}_{1}\cdot\overline{\mathbf{D}}_{1}\cdot\mathsf{A}_{\parallel}$
$\displaystyle=\lambda_{\parallel}\overline{\mathbf{C}}_{1}\cdot\overline{\mathbf{U}}_{1}^{-1}\cdot\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\parallel}.$
(46)

Due to (40), the LHS of the above vanishes. Thus

$\displaystyle\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{\parallel}=0.$ (47)

####  III-B1 Proof 1

Taking the left dot product of (43) with $\mathsf{A}_{\varepsilon\perp}^{T}$
and using the symmetry property between divergence and gradient operators, we
have

$\displaystyle\alpha\mathsf{A}_{\varepsilon\perp}^{T}\cdot\overline{\mathbf{D}}_{1}^{T}\cdot\overline{\mathbf{D}}_{1}\cdot\mathsf{A}_{\parallel}$
$\displaystyle=\lambda_{\parallel}\mathsf{A}_{\varepsilon\perp}^{T}\cdot\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\parallel}.$
(48)

Due to (45) and its transpose, LHS of the above vanishes. Thus

$\displaystyle\mathsf{A}_{\varepsilon\perp}^{T}\cdot\overline{\mathbf{U}}\cdot\mathsf{A}_{\parallel}=0.$
(49)

which is analogous to (29).

####  III-B2 Proof 2

Now, taking the left dot product of (42) with $\mathsf{A}_{\parallel}^{T}$,
where † denotes conjugate transpose, and using the symmetry property of curl
operator, we have

$\displaystyle\mathsf{A}_{\parallel}^{T}\cdot\overline{\mathbf{C}}_{1}^{T}\cdot\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}$
$\displaystyle=\lambda_{\varepsilon\perp}\mathsf{A}_{\parallel}^{T}\cdot\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}.$
(50)

Due to (47) and its transpose, LHS of the above vanishes. Thus

$\displaystyle\mathsf{A}_{\parallel}^{T}\cdot\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}=0$
(51)

which is homomorphic to (32). However, we can also left dot product (42) by
$\mathsf{A}_{\parallel}^{\dagger}$, and get

$\displaystyle\mathsf{A}_{\parallel}^{\dagger}\cdot\overline{\mathbf{C}}_{1}^{T}\cdot\overline{\mathbf{C}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}$
$\displaystyle=\lambda_{\varepsilon\perp}\mathsf{A}_{\parallel}^{\dagger}\cdot\overline{\mathbf{U}}_{1}\cdot\mathsf{A}_{\varepsilon\perp}.$
(52)

Due to (47) and its Hermitian conjugate, and the realness of the curl operator
under finite difference method, LHS of the above vanishes. Thus

$\displaystyle\mathsf{A}_{\parallel}^{\dagger}\cdot\overline{\mathbf{U}}\cdot\mathsf{A}_{\varepsilon\perp}=0$
(53)

which is analogous to (33). The above shows that the finite difference
approximation of Maxwell’s equations is homomorphic to the continuum case if
done properly.

###  III-C Anisotropic Case

The mathematical proof in the previous two subsections still applies in
anisotropic media. Consider the wave equation for vector potential
$\mathbf{A}$ in dispersionless inhomogeneous anisotropic media with modified
generalized Lorenz gauge in frequency domain [2]

$\displaystyle\nabla\times\overline{\bm{\mu}}^{-1}\nabla\times\mathbf{A}-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\overline{\bm{\varepsilon}}\cdot\nabla\left(\nabla\cdot\overline{\bm{\varepsilon}}\cdot\mathbf{A}\right)=\omega^{2}\overline{\bm{\varepsilon}}\cdot\mathbf{A}.$
(54)

The mathematical proof of mode orthogonality follows exactly as Section III-A
for continuous space and as Section III-B for discrete space, with the
requirement that

$\displaystyle\overline{\bm{\varepsilon}}=\overline{\bm{\varepsilon}}^{T}$
(55)

which implies that the material is reciprocal.

##  IV Numerical Results

###  IV-A Numerical Demonstration of Completeness

In this section, we first demonstrate that the two sets of modes given in the
previous section are complete in the discrete space. We solve (15)
numerically. Finite difference method is used corresponding to the discrete
space proof in Section III. Uniform discretization is used, and the
discretization scheme follows [39, 15], as described in Appendix A.

Figure 3: Geometries of the numerical examples with dimensions.

For an empty or homogeneous structure, the proposed generalized Helmholtz
decomposition degenerates into the usual Helmholtz decomposition. To
demonstrate it, we consider a rectangular $7\times 8\times 6$ cavity
homogeneously filled with a $\varepsilon_{r}=3$ material, as shown in Fig.
3(a). Discretization step is $0.5$, resulting in $6938\times 6938$ matrices.
We use natural unit so that $\varepsilon_{0}=\mu_{0}=1$. After solving the
generalized eigenvalue problem, we take divergence and curl on each mode to
identify the two classes of modes. The results are displayed in Fig. 4. Note
that the indices of curl-free modes and div-$\varepsilon$-free modes are
skipped in (b) and (c), respectively. The total number of
div-$\varepsilon$-free and curl-free modes equals to the dimension of the
eigenvalue problem. Thus it clearly shows two distinct groups of modes
corresponding to the Helmholtz decomposition, which complete the discrete
space, and are consistent across field-based and potential-based formulations.

Figure 4: Two set of modes found in the geometry described by Fig. 3(a) can be
clearly distinguished. (a) The absolute value of eigenvalues are shown for
both potential-based formulation and $\mathbf{E}$-$\mathbf{H}$ formulation. It
can be seen clearly that the potential-based formulation does not have null
space. For the potential-based formulation, only the div-$\varepsilon$-free
(or generalized transverse) modes are shown in (b), and only the curl-free (or
longitudinal) modes are shown in (c). Totally there are 4793
div-$\varepsilon$-free modes in (b) and 2145 curl-free modes in (c). The
latter is consistent with the number of null-space eigenvalues with
$\mathbf{E}$-$\mathbf{H}$ formulation in (a). The total number of both
families of modes is equal to the dimension of the eigenvalue problem. The
indices of curl-free modes and div-$\varepsilon$-free modes are skipped in (b)
and (c), respectively. It can be seen that the div-$\varepsilon$-free modes
satisfy $\nabla\cdot\varepsilon\mathbf{A}=0$, while the curl-free modes
satisfy $\nabla\times\mathbf{A}=0$. Blue open circles: for each mode, take
$\nabla\cdot\varepsilon\mathbf{A}$, and record the maximum value of
abs($\nabla\cdot\varepsilon\mathbf{A}$) in the entire grid. Red stars: for
each mode, take $\nabla\times\mathbf{A}$, and record the maximum magnitude of
all components in the entire grid.

We then consider a more general example with anisotropic material. As shown in
3(b), we construct a cylindrical cavity with radius $r=2$ and height $h=2$. It
is filled with a block of anisotropic material centered at $(1,0.4,0.2)$ with
dimension $0.6,0.4,0.8$ in $x,y,z$ directions, respectively. Discretization
step is $0.2$, resulting in $5822\times 5822$ matrices. The permittivity
tensor of the anisotropic material is

$\displaystyle\overline{\bm{\varepsilon}}=\begin{pmatrix}2&0&1+3i\\\
0&3+0.2i&0\\\ 1+3i&0&2\end{pmatrix}$ (56)

which satisfies reciprocity. The matrix of anisotropic permittivity is treated
the same way as [42]. The discretized operators satisfy the discrete vector
identities (39) and (40).

Mode degeneracy is observed after solving the generalized eigenvalue problem.
We apply the analysis described in Appendix B to separate degenerate modes.
The results are then displayed in Fig. 5. We find that the modes can be
categorized into two distinct groups:

  1. 1.

div-$\varepsilon$-free (or generalized transverse) modes (Fig. 5(a)):
$\nabla\cdot\varepsilon\mathbf{A}=0$ for the entire grid;

  2. 2.

Curl-free (or longitudinal) modes (Fig. 5(b)): $\nabla\times\mathbf{A}=0$ for
the entire grid.666The curl-free modes can be eliminated by $\Phi=0$ gauge.
This will be explained in future work.

It can be seen that each mode is either div-$\varepsilon$-free, or curl-free.
The div-$\varepsilon$-free modes are consistent with the modes found by
$\mathbf{E}$-$\mathbf{H}$ formulation since (15) returns to the wave equation
of $\mathbf{E}$ field [2, 43]. The curl-free modes have zero $\mathbf{E}$ and
$\mathbf{H}$ field by the fact that $\mathbf{B}=\nabla\times\mathbf{A}=0$.

In both numerical examples, the two classes of modes form a complete basis.
Although the completeness of basis expanded by the eigenvectors in the
discrete space from the numerical examples does not necessarily guarantee the
completeness of the two sets of modes in the continuous Hilbert space, one
could use convergence and existence of the numerical solution to a partial
differential equation with finite difference approximation [44] to imply
completeness of these modes in the continuous Hilbert space. Thus the
numerical proof could motivate further study in the continuous Hilbert space
as future work.

Figure 5: Two set of modes found in the geometry described by Fig. 3(b) can be
clearly distinguished. Totally there are 4049 div-$\varepsilon$-free modes and
1773 curl-free modes. The description of the figure follows the caption of
Fig. 4(b) and (c).

###  IV-B Numerical Validation of Orthogonality

We then proceed to validate the orthogonality conditions. For the numerical
example described by Fig. 3(b), we calculate the four sets of inner products
between the div-$\varepsilon$-free modes and curl-free modes defined in
Section III-B after separating the degenerate modes. The colormaps of the
resulting matrices are shown in Fig. 6, and the maximum magnitude of the
matrix elements are shown in Table I. For readability, only the first
$100\times 100$ elements of each matrix are displayed in Fig. 6, but these
results are representative of the full matrices. It can be seen from Fig. 6
and Table I that

$\displaystyle\langle\mathbf{A}_{\varepsilon\perp},\varepsilon\mathbf{A}_{\parallel}\rangle=0$
(57)
$\displaystyle\langle\mathbf{A}_{\parallel},\varepsilon\mathbf{A}_{\varepsilon\perp}\rangle=0$
(58)
$\displaystyle\langle\mathbf{A}_{\varepsilon\perp}^{*},\varepsilon\mathbf{A}_{\parallel}\rangle\neq
0$ (59)
$\displaystyle\langle\mathbf{A}_{\parallel}^{*},\varepsilon\mathbf{A}_{\varepsilon\perp}\rangle=0$
(60)

which are consistent with the mathematical proofs in Section III. This result,
supplemented by the numerical proof of completeness in the previous
subsection, completes our demonstration of the generalized Helmholtz
decomposition.

Figure 6: Color map of four sets of inner products from the numerical example described by Fig. 3(b). (a) $\overline{\mathbf{V}}_{\varepsilon\perp}^{T}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\parallel}$; (b) $\overline{\mathbf{V}}_{\parallel}^{T}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\varepsilon\perp}$; (c) $\overline{\mathbf{V}}_{\varepsilon\perp}^{\dagger}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\parallel}$; (d) $\overline{\mathbf{V}}_{\parallel}^{\dagger}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\varepsilon\perp}$. Here $\overline{\mathbf{V}}_{\varepsilon\perp}$($\overline{\mathbf{V}}_{\parallel}$) is a matrix where each column vector is a div-$\varepsilon$-free (curl-free) mode. Each pixel represents the magnitude of the corresponding matrix element. The indices are the indices of matrix elements. Note the different scales of the colorbars. The values of the resulting matrices are essentially 0 in (a), (b), and (d).  TABLE I:  Orthogonality results from the second numerical example. The maximum magnitude of matrix elements in the corresponding matrices are shown in the table. Here $|\cdot|$ denotes element-wise magnitude. max$(|\overline{\mathbf{V}}_{\varepsilon\perp}^{T}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\parallel}|)$ | 5.1694e-10
---|---
max$(|\overline{\mathbf{V}}_{\parallel}^{T}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\varepsilon\perp}|)$ | 5.1694e-10
max$(|\overline{\mathbf{V}}_{\varepsilon\perp}^{\dagger}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\parallel}|)$ | 0.5779
max$(|\overline{\mathbf{V}}_{\parallel}^{\dagger}\cdot\overline{\bm{\varepsilon}}\cdot\overline{\mathbf{V}}_{\varepsilon\perp}|)$ | 7.2430e-10

##  V Conclusion

In this work, we combine mathematical proof and numerical results to
demonstrate the generalization of Helmholtz decomposition to inhomogeneous
media. The two families of fields in the generalized Helmholtz decomposition,
div-$\varepsilon$-free and curl-free fields, are connected to the two classes
of solutions to the vector potential wave equation. The div-$\varepsilon$-free
and curl-free solutions form a complete set of the basis, as demonstrated in
the numerical results and its discussion. They are also orthogonal to each
other. The div-$\varepsilon$-free field manifests charge-free condition, while
the curl-free field is associated with charges. In the literature, the curl-
free fields are often eliminated by setting $\Phi=0$. However, when there are
charges present, we need both families of solutions. It also provides a
theoretical background to the numerical quantization [1] based on the
generalized Lorenz gauge. Future work will be devoted to the application of
$\Phi=0$ gauge and rigorously determining suitable conditions for eliminating
the curl-free modes in general inhomogeneous media while preserving full rank
of the resulting system matrix for practical applications.

##  Appendix A Discrete Vector Calculus on a Grid

Discrete vector calculus on a Yee lattice [38, 35] is the basis for the proof
of orthogonality in discrete space. We give a brief overview of the
discretization. The details can be found in [39].

The staggered grid with both primal (solid line) and dual (dashed line) grids
are shown in Fig. 7. The dual grid is shifted half a grid width in all three
directions from the primal grid. On each grid, we place a 3D vector on the
face of the cube, and a scalar in the center of the cube.

Figure 7: Schematic of the grid. The primal and dual grids are shifted by half
grid step. The solution vector $\mathbf{A}$ is placed on the face of the
primal grid.

We index the grid point at the cube center of the primal grid by integers
$(i,j,k)$. Below we explicitly formulate the discrete differential operations,
where central difference is used.

A discrete gradient $\mathbf{g}=\nabla f$ on the primal grid is given by

$\displaystyle g_{x}^{i+0.5,j,k}$
$\displaystyle=\frac{f^{i+1,j,k}-f^{i,j,k}}{\Delta x},$ (61) $\displaystyle
g_{y}^{i,j+0.5,k}$ $\displaystyle=\frac{f^{i,j+1,k}-f^{i,j,k}}{\Delta y},$
(62) $\displaystyle g_{z}^{i,j,k+0.5}$
$\displaystyle=\frac{f^{i,j,k+1}-f^{i,j,k}}{\Delta z}.$ (63)

A discrete divergence $f=\nabla\cdot\mathbf{g}$ on the primal grid is given by

$\displaystyle f^{i,j,k}$
$\displaystyle=\frac{g_{x}^{i+0.5,j,k}-g_{x}^{i-0.5,j,k}}{\Delta x}$
$\displaystyle+\frac{g_{y}^{i,j+0.5,k}-g_{y}^{i,j-0.5,k}}{\Delta y}$
$\displaystyle+\frac{g_{z}^{i,j,k+0.5}-g_{z}^{i,j,k-0.5}}{\Delta z}.$ (64)

A discrete curl $\mathbf{h}=\nabla\times\mathbf{g}$ on a vector on the primal
grid resulting in a vector on the dual grid is given by

$\displaystyle h_{x}^{i,j+0.5,k+0.5}$
$\displaystyle=\frac{g_{z}^{i,j+1,k+0.5}-g_{z}^{i,j,k+0.5}}{\Delta y}$
$\displaystyle-\frac{g_{y}^{i,j+0.5,k+1}-g_{y}^{i,j+0.5,k}}{\Delta z},$ (65)
$\displaystyle h_{y}^{i+0.5,j,k+0.5}$
$\displaystyle=\frac{g_{x}^{i+0.5,j,k+1}-g_{x}^{i+0.5,j,k}}{\Delta z}$
$\displaystyle-\frac{g_{z}^{i+1,j,k+0.5}-g_{z}^{i,j,k+0.5}}{\Delta x},$ (66)
$\displaystyle h_{z}^{i+0.5,j+0.5,k}$
$\displaystyle=\frac{g_{y}^{i+1,j+0.5,k}-g_{y}^{i,j+0.5,k}}{\Delta x}$
$\displaystyle-\frac{g_{x}^{i+0.5,j+1,k}-g_{x}^{i+0.5,j,k}}{\Delta y}.$ (67)

A discrete curl $\mathbf{g}=\nabla\times\mathbf{h}$ on a vector on the dual
grid resulting in a vector on the primal grid is given by

$\displaystyle g_{x}^{i+0.5,j,k}$
$\displaystyle=\frac{h_{z}^{i+0.5,j+0.5,k}-h_{z}^{i+0.5,j-0.5,k}}{\Delta y}$
$\displaystyle-\frac{h_{y}^{i+0.5,j,k+0.5}-h_{y}^{i+0.5,j,k-0.5}}{\Delta z},$
(68) $\displaystyle g_{y}^{i,j+0.5,k}$
$\displaystyle=\frac{h_{x}^{i,j+0.5,k+0.5}-h_{x}^{i,j+0.5,k-0.5}}{\Delta z}$
$\displaystyle-\frac{h_{z}^{i+0.5,j+0.5,k}-h_{z}^{i-0.5,j+0.5,k}}{\Delta x},$
(69) $\displaystyle g_{z}^{i,j,k+0.5}$
$\displaystyle=\frac{h_{y}^{i+0.5,j,k+0.5}-h_{y}^{i-0.5,j,k+0.5}}{\Delta x}$
$\displaystyle-\frac{h_{x}^{i,j+0.5,k+0.5}-h_{x}^{i,j-0.5,k+0.5}}{\Delta y}.$
(70)

##  Appendix B Analysis of Degenerate Modes

Given a set of eigenmodes $\mathbf{A}_{1},\mathbf{A}_{2},\dots,\mathbf{A}_{n}$
with corresponding eigenvalues $\lambda_{1},\lambda_{2},\dots,\lambda_{n}$
obtained from the generalized eigenvalue problem, the space expanded by

$\displaystyle\frac{1}{\lambda_{i}}\varepsilon^{-1}\nabla\times\frac{1}{\mu}\nabla\times\mathbf{A}_{i}$
(71)

has the same dimension as the number of independent div-$\varepsilon$-free
modes. This is because the curl-free modes are eliminated by the operation due
to (25) $\lambda_{\parallel}\nabla\times\mathbf{A}_{\parallel}=0$, and the
div-$\varepsilon$-free modes are recovered due to (22)
$\nabla\times\frac{1}{\mu}\nabla\times\mathbf{A}_{\varepsilon\perp}=\lambda_{\varepsilon\perp}\varepsilon\mathbf{A}_{\varepsilon\perp}$.
Besides, the space expanded by

$\displaystyle-\frac{1}{\lambda_{i}}\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\nabla\nabla\cdot\varepsilon\mathbf{A}_{i}$
(72)

has the same dimension as the number of independent curl-free modes. This is
because the div-$\varepsilon$-free modes are eliminated by the operation due
to (23)
$\lambda_{\varepsilon\perp}\nabla\cdot\varepsilon\mathbf{A}_{\varepsilon\perp}=0$,
and the curl-free modes are recovered due to (24)
$-\frac{1}{\varepsilon_{0}^{2}\mu_{0}}\varepsilon\nabla\nabla\cdot\varepsilon\mathbf{A}_{\parallel}=\lambda_{\parallel}\varepsilon\mathbf{A}_{\parallel}$.

The dimension of the two spaces from the numerical example described by Fig.
3(b) is shown in Fig. 8. The number of nonzero singular values indicate the
dimension of the space spanned by a set of vectors. It can be seen that the
sum of the dimension of the space of (71) and (72) equals the dimension of the
problem. Thus it proves that the div-$\varepsilon$-free and curl-free modes as
defined in the main text form a complete basis of the discrete space.
Otherwise, if there exists an eigenmode that is nonzero after both (71) and
(72) operations, it must contribute one “extra dimension”.

Degenerate modes with the same eigenvalue can also be separated based on this
analysis. Given $n$ degenerate modes
$\mathbf{A}_{1},\mathbf{A}_{2},\dots,\mathbf{A}_{n}$ with the same eigenvalue,
we first learn the number of independent div-$\varepsilon$-free and curl-free
modes by evaluating the dimension of the subspace expanded by (71) and (72),
as shown in Fig. 9. The div-$\varepsilon$-free modes are isolated by (71)
because of (22), and the curl-free modes are isolated by (72) because of (24).
Gram-Schmidt orthogonalization is then performed within each subspace to
extract the independent div-$\varepsilon$-free and curl-free modes,
respectively.

Figure 8: The eigenvalues and the resulting solution space of the geometry
described by Fig. 3(b) are shown. (a) Eigenvalues are shown in blue curve;
Orange upward spikes denote a cluster of degenerate modes. (b) The singular
values of the matrix constructed by lining up
$\overline{\bm{\varepsilon}}^{-1}\cdot\nabla\times\frac{1}{\mu}\nabla\times\mathbf{A}_{i}$
are shown. It can be seen that there are 4049 nonzero singular values. (c) The
singular values of the matrix constructed by lining up
$-\nabla\nabla\cdot\overline{\bm{\varepsilon}}\cdot\mathbf{A}_{i}$ are shown.
It can be seen that there are 1773 nonzero singular values. The number of
nonzero singular values indicate the dimension of the space spanned by a set
of vectors. The total number of nonzero singular values from (b) and (c) is
5822, equal to the dimension of the problem.  Figure 9: A closer look at the
first cluster of degenerate modes (index 1428 to 1451) in Fig. 8(a). Similar
results are observed for the other clusters of degenerate modes. Singular
values of (a) and (b) are obtained in the same way as Fig. 8(b) and (c). The
total number of nonzero singular values from (a) and (b) is 24, equal to the
dimension of the degenerate subspace.

## References

  * [1] D.-Y. Na, J. Zhu, W. C. Chew, and F. L. Teixeira, “Quantum information preserving computational electromagnetics,” _Physical Review A_ , vol. 102, no. 1, p. 013711, 2020.
  * [2] W. C. Chew, “Vector potential electromagnetics with generalized gauge for inhomogeneous media: Formulation,” _Progress In Electromagnetics Research_ , vol. 149, pp. 69–84, 2014.
  * [3] J. B. Manges and Z. J. Cendes, “A generalized tree-cotree gauge for magnetic field computation,” _IEEE Transactions on Magnetics_ , vol. 31, no. 3, pp. 1342–1347, 1995.
  * [4] S.-C. Lee, J.-F. Lee, and R. Lee, “Hierarchical vector finite elements for analyzing waveguiding structures,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 51, no. 8, pp. 1897–1905, 2003.
  * [5] L. Xue and D. Jiao, “Fast method for accelerating time-domain solutions of ill-conditioned electromagnetic problems by changing curl–curl operator to Laplacian,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 69, no. 11, pp. 4613–4622, 2021.
  * [6] D. Wilton, J. Lim, and S. Rao, “A novel technique to calculate the electromagnetic scattering by surfaces of arbitrary shape,” in _URSI Radio Science Meeting_ , 1993, p. 322.
  * [7] M. Burton and S. Kashyap, “A study of a recent, moment-method algorithm that is accurate to very low frequencies,” _Applied Computational Electromagnetics Society Journal_ , vol. 10, pp. 58–68, 1995.
  * [8] W.-L. Wu, A. W. Glisson, and D. Kajfez, “A study of two numerical solution procedures for the electric field integral equation at low frequency,” _Applied Computational Electromagnetics Society Journal_ , vol. 10, no. 3, pp. 69–80, 1995.
  * [9] G. Vecchi, “Loop-star decomposition of basis functions in the discretization of the EFIE,” _IEEE Transactions on Antennas and Propagation_ , vol. 47, no. 2, pp. 339–346, 1999.
  * [10] J.-S. Zhao and W. C. Chew, “Integral equation solution of Maxwell’s equations from zero frequency to microwave frequencies,” _IEEE Transactions on Antennas and Propagation_ , vol. 48, no. 10, pp. 1635–1645, 2000\.
  * [11] J. Mautz and R. Harrington, “An e-field solution for a conducting surface small or comparable to the wavelength,” _IEEE Transactions on Antennas and Propagation_ , vol. 32, no. 4, pp. 330–339, 1984.
  * [12] Z. G. Qian and W. C. Chew, “An augmented electric field integral equation for high-speed interconnect analysis,” _Microwave and Optical Technology Letters_ , vol. 50, no. 10, pp. 2658–2662, 2008.
  * [13] Y.-L. Li, S. Sun, Q. I. Dai, and W. C. Chew, “Finite element implementation of the generalized-Lorenz gauged A-$\Phi$ formulation for low-frequency circuit modeling,” _IEEE Transactions on Antennas and Propagation_ , vol. 64, no. 10, pp. 4355–4364, 2016.
  * [14] R. Dyczij-Edlinger, G. Peng, and J.-F. Lee, “A fast vector-potential method using tangentially continuous vector finite elements,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 46, no. 6, pp. 863–868, 1998.
  * [15] C. J. Ryu, A. Y. Liu, W. E. I. Sha, and W. C. Chew, “Finite-difference time-domain simulation of the Maxwell–Schrödinger system,” _IEEE Journal on Multiscale and Multiphysics Computational Techniques_ , vol. 1, pp. 40–47, 2016.
  * [16] T. E. Roth and C. A. Braun, “Theoretical and numerical study of wave port boundary conditions for potential-based finite element methods,” 2022, manuscript submitted.
  * [17] P. De Doncker, “A volume/surface potential formulation of the method of moments applied to electromagnetic scattering,” _Engineering Analysis with Boundary Elements_ , vol. 27, no. 4, pp. 325–331, 2003.
  * [18] F. Vico, M. Ferrando, L. Greengard, and Z. Gimbutas, “The decoupled potential integral equation for time-harmonic electromagnetic scattering,” _Communications on Pure and Applied Mathematics_ , vol. 69, no. 4, pp. 771–812, 2016.
  * [19] Q. S. Liu, S. Sun, and W. C. Chew, “A potential-based integral equation method for low-frequency electromagnetic problems,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 3, pp. 1413–1426, 2018.
  * [20] J. Li, X. Fu, and B. Shanker, “Decoupled potential integral equations for electromagnetic scattering from dielectric objects,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 3, pp. 1729–1739, 2018.
  * [21] B. Khalichi, Ö. Ergül, and V. B. Ertürk, “Broadband solutions of potential integral equations with NSPWMLFMA,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 6, pp. 4307–4312, 2019.
  * [22] T. E. Roth and W. C. Chew, “Development of stable A-$\Phi$ time-domain integral equations for multiscale electromagnetics,” _IEEE Journal on Multiscale and Multiphysics Computational Techniques_ , vol. 3, pp. 255–265, 2018\.
  * [23] ——, “Lorenz gauge potential-based time domain integral equations for analyzing subwavelength penetrable regions,” _IEEE Journal on Multiscale and Multiphysics Computational Techniques_ , vol. 6, pp. 24–34, 2021\.
  * [24] R. Chen, H. A. Ulku, F. P. Andriulli, and H. Bagci, “On the low-frequency behavior of vector potential integral equations for perfect electrically conducting scatterers,” _IEEE Transactions on Antennas and Propagation_ , 2022.
  * [25] Y. Aharonov and D. Bohm, “Significance of electromagnetic potentials in the quantum theory,” _Physical Review_ , vol. 115, no. 3, p. 485, 1959.
  * [26] D.-Y. Na, J. Zhu, and W. C. Chew, “Diagonalization of the Hamiltonian for finite-sized dispersive media: Canonical quantization with numerical mode decomposition,” _Physical Review A_ , vol. 103, no. 6, p. 063707, 2021.
  * [27] W. C. Chew, A. Y. Liu, C. Salazar-Lazaro, and W. E. I. Sha, “Quantum electromagnetics: A new look—Part II,” _IEEE Journal on Multiscale and Multiphysics Computational Techniques_ , vol. 1, pp. 85–97, 2016.
  * [28] G. B. Arfken, H.-J. Weber, and F. E. Harris, _Mathematical Methods for Physicists (Waltham, MA_. Elsevier, 2013\.
  * [29] R. E. Collin, _Foundations for Microwave Engineering_. John Wiley & Sons, 2007.
  * [30] J. D. Jackson, _Classical Electrodynamics_. American Association of Physics Teachers, 1999.
  * [31] W. C. Chew, _Waves and Fields in Inhomogeneous Media_. IEEE press, 1995.
  * [32] O. Axelsson and V. A. Barker, _Finite Element Solution of Boundary Value Problems: Theory and Computation_. SIAM, 2001.
  * [33] V. Rumsey, “Reaction concept in electromagnetic theory,” _Physical Review_ , vol. 94, no. 6, p. 1483, 1954.
  * [34] W. C. Chew, “Lectures on electromagnetic field theory,” 2022. [Online]. Available: https://engineering.purdue.edu/wcchew/ece604f22/EMFTEDX103122R.pdf
  * [35] A. Taflove and S. C. Hagness, _Computational Electrodynamics: The Finite-Difference Time-Domain Method_. Artech House, 2005.
  * [36] J.-M. Jin, _The Finite Element Method in Electromagnetics_. John Wiley & Sons, 2015.
  * [37] S. C. Chen and W. C. Chew, “Electromagnetic theory with discrete exterior calculus,” _Progress In Electromagnetics Research_ , vol. 159, pp. 59–78, 2017.
  * [38] K. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” _IEEE Transactions on antennas and propagation_ , vol. 14, no. 3, pp. 302–307, 1966.
  * [39] W. C. Chew, “Electromagnetic theory on a lattice,” _Journal of Applied Physics_ , vol. 75, no. 10, pp. 4843–4850, 1994.
  * [40] F. L. Teixeira and W. Chew, “Lattice electromagnetic theory from a topological viewpoint,” _Journal of mathematical physics_ , vol. 40, no. 1, pp. 169–187, 1999.
  * [41] T. Weiland, “Time domain electromagnetic field computation with finite difference methods,” _International Journal of Numerical Modelling: Electronic Networks, Devices and Fields_ , vol. 9, no. 4, pp. 295–319, 1996.
  * [42] R. C. Rumpf, C. R. Garcia, E. A. Berry, and J. H. Barton, “Finite-difference frequency-domain algorithm for modeling electromagnetic scattering from general anisotropic objects,” _Progress In Electromagnetics Research B_ , vol. 61, pp. 55–67, 2014.
  * [43] A. Raman and S. Fan, “Photonic band structure of dispersive metamaterials formulated as a Hermitian eigenvalue problem,” _Physical Review Letters_ , vol. 104, no. 8, p. 087401, 2010.
  * [44] R. J. LeVeque, _Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems_. SIAM, 2007.