INTEGRAL EQUATIONS IN COMPUTATIONAL ELECTROMAGNETICS AND THEIR
IMPLEMENTATIONS

By

Luke Baumann

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Electrical and Computer Engineering—Doctor of Philosophy
Computational Mathematics, Science, and Engineering—Dual Major

2023

ABSTRACT

Integral equations in Computational Electromagnetics (CEM) are one branch of a diverse

field. There are many methods to solve for electromagnetic scattering and transmission, with

boundary integral equations being one of the most efficient. This is due to only needing to

discretize the object’s surface, leading to smaller, dense systems as opposed to the larger,

sparse systems encountered with Finite Element Method (FEM). Combining the boundary

integral method with FEM leads to the creatively named Finite Element Boundary Integral

(FEBI) method. It can use the more appropriate method as needed for a given region of

space. We turn our focus to boundary integral methods and their implementations.

The subfield of boundary integral equations comprises many subparts, including for-

mulations, representations, testing, singularity treatment, acceleration techniques, solvers,

preconditioning, and others. In this thesis, I will present several new and existing formulations

using the same formulation framework, demonstrate how to perform the integrals for analytic

and piecewise basis and testing functions, modify acceleration techniques for various integral

equations, and present supporting results.

The new formulations are well-conditioned, free from traditional breakdowns, and compa-

rable to state-of-the-art formulations. Most of the implementation of all the formulations

presented is shared to limit unintended comparisons.

Copyright by
LUKE BAUMANN
2023

To my loving wife
and kids, Ocean and Bode,
you are everything

iv

ACKNOWLEDGEMENTS

I want to acknowledge my advisor first and foremost for his continued support through the

years. He is incredibly passionate about excellence and has shaped how I approach life, work,

and research. It has been my joy and great privilege to study under him.

A special thanks to my committee members, Drs. Aktulga, Qian, Rothwell, Van, and

Zayernouri. I made it a point to attend official and unofficial office hours as much as possible

for each of them and have a broader and deeper understanding because of it.

I owe so much thanks to my fellow group members, both past and present. I have many

memories of that first semester in the group with Drs. Li, Hughey, O’Connor, Crawford,

Alsnayyan, Glisson. I want to thank the current members of the group, Jacob Hawkins, Chad

Moorman, Omkar Ramachandranharishankar, and Elliot Lu, for making this experience more

fulfilling and collaborative.

The Electromagnetics Research Group at Michigan State University is incredible. The

faculty have been extremely accessible for advice and eager to help. The students are brilliant

and motivated.

Most of all, I would like to thank my family. My wife, for all of her sacrifices, encouragement,

and love. I would not have finished this journey without her. My wonderful kids Ocean and

Bode who are probably halfway to a degree from all of the online classes and meetings they

have been present for. My parents and brother who have always supported and loved me.

v

TABLE OF CONTENTS

CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 What is a Formulation? . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 2

FORMULATIONS . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1
2.2 Common Formulation Framework . . . . . . . . . . . . . . . . . . . . . .
2.3 Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 DFIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 DPIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 LC-CFIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 3

DISCRETIZATION . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Analytic Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Piecewise Analysis
3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 4

ACCELERATION . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Multilevel Fast Multipole Algorithm (MLFMA) Background . . . . . . .
4.3
Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 5

CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
11

12
12
14
17
18
22
31
33

34
34
35
52
63

64
64
64
65
67
71

72
72

74

vi

CHAPTER 1

INTRODUCTION

This thesis builds on decades of research in the field of CEM, a vast and varied field that

spans physics, mathematics, engineering, and computer science. This conflux of disciplines

makes for a fertile research environment where there is rarely a lack of topics to expand

upon. This thesis is razor-focused and aimed at developing new surface integral equation

formulations capable of solving problems over a broad band of frequencies. Background

on surface equivalence theorem, surface integral equations, direct methods, and current

formulations is necessary to fully appreciate this body of work. As an introduction, a brief

overview of the key concepts will be covered.

1.1 What is a Formulation?

In our context, a formulation is a boundary value problem that can be used to model

electromagnetic physics. A formulation can be separated into two parts: the differential

equation that governs the behavior of fields interior to a region and boundary conditions

relating fields between touching regions.

1.1.1 Differential Equations

Maxwell’s equations are the cornerstone of all electromagnetic formulations because they

are the differential equations that govern all electromagnetic phenomena. For the scenarios

that we are interested in, we can specialize Maxwell’s equations for homogeneous, source-free

regions with linear materials

∇ × E = − jκηH

∇ × H =

jκ
η

E

∇ · E =0

∇ · H =0

1

(1.1)

(1.2)

(1.3)

(1.4)

with electric field E, magnetic field H, imaginary unit j =
wave impedance η = (cid:112) µ

ϵ , permittivity ϵ, permeability µ. Other material constants used

√

−1, wavenumber κ = ω

√

µϵ,

throughout include relative permittivity ϵr and permeability µr defined with respect to the

permittivity and permeability of free space ϵ = ϵrϵ0, µ = µrµ0 and the refractive index

√

n =

µrϵr. These equations describe how the fields interact with each other and can be

combined to see that each field satisfies the vector Helmholtz differential equation

∇ × E = − jκηH

∇ × ∇ × E = − jκη∇ × H

∇ × ∇ × E = − jκη

jκ
η

E

∇ (∇ · E) − ∇2E =κ2E

∇2E + κ2E =0

∇ × H =

∇ × ∇ × H =

∇ × ∇ × H =

jκ
η
jκ
η
jκ
η

E

∇ × E

(−jκηH)

∇ (∇ · H) − ∇2H =κ2H

∇2H + κ2H =0

(1.5)

(1.6)

(1.7)

(1.8)

(1.9)

(1.10)

(1.11)

(1.12)

(1.13)

(1.14)

1.1.2 Boundary Conditions

The boundary conditions can be formed by examining a pill-box and Stokesian loop across

an interface separating two regions of different materials as is typically done [Balanis, 2012].

Boundary conditions for more complex junctions are easily developed as in [Yl¨a-Oijala et al.,

2005]. Still, for simplicity, we restrict ourselves to a single object, dividing space into two

regions with a single interface. Once we denote values associated with each region with the

proper subscript p and equip the interface with two normal vectors that point into their

2

respective regions, we can write the boundary conditions as

E1 × ˆn1 + E2 × ˆn2 =Ms

ϵr1ˆn1 · E1 + ϵr2ˆn2 · E2 =

ˆn1 × H1 + ˆn2 × H2 =Js
ρs
ϵ0
ρms
µ0

µr1ˆn1 · H1 + µr2ˆn2 · H2 =

(1.15)

(1.16)

(1.17)

(1.18)

with electric surface current and charge densities Js and ρs along with fictitious magnetic

surface current and charge densities Ms and ρms as is typically done.

1.1.3 Surface Integral Formulations

Surface integral formulations reformulate Maxwell’s equations into an integrodifferential

equation, discretize the boundary, and introduce unknown surface sources responsible for

the scattered fields within each region. There are two classes of surface integral equations:

direct and indirect methods. For direct methods, the sources are directly related to physical

quantities. This means the boundary conditions can be imposed on the surface sources.

For indirect methods, the sources have no direct physical meaning but can still represent

the scattered field within each region. Because there is no physical meaning to the sources,

the boundary conditions are imposed on the measurement of the field or charge on the

boundary. This thesis will exclusively utilize the direct method, but equivalent indirect

method corollaries exist.

For all of the formulations, we decompose the fields into the incident and scattered

parts xt

p = xi

p + xs

p with x being any vector that satisfies the vector Helmholtz equation

∇2x + κ2x = 0 and superscripts t, i, and s denoting the total, incident, and scattered part.

The incident field is given due to a primary source. This incident field induces a reactionary

scattered field such that the boundary conditions on the total fields are satisfied at the

interface. In addition, the scattered field must also satisfy a radiating condition to generate a

practical and unique solution. The scattered field is written such that it is due to unknown

sources on the surface of the boundary. The integral equation used to represent the scattered

3

field on the surface is derived from Green’s vector identity and the vector Helmholtz equation.

We let xs

p be the scattered field, and Gp be the Green’s function for the inhomogeneous

Helmholtz equation

Gp (r, r′) =

exp (−jκpR)
4πR

∇2 Gp (r, r′) + κ2

p Gp (r, r′) = −δ (R)

(1.19a)

(1.19b)

with r as the observation point, r′ as the source point, R = r − r′, and R = ∥R∥2. After

much manipulation, we can arrive at what is known as Green’s representation theorem

p = Sp ◦ (cid:0)ˆnp × ∇ × xt
xs

p

(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × xt

p

(cid:1) − ∇Sp ◦ (cid:0)ˆnp · xt

p

(cid:1) − Sp ◦ (cid:0)ˆnp∇ · xt

p

(cid:1) (1.20)

where

(cid:90)

(cid:90)

Sp ◦ fp =

Sp ◦ f p =

Gp (r, r′) fp (r′) dS′

Gp (r, r′) f p (r′) dS′.

(1.21a)

(1.21b)

is the Single-Layer potential for vectors and scalars. We identify four surface sources per

region ˆnp × ∇ × xt

p, ˆnp × xt

p, ˆnp · xt

p, and ∇ · xt

p. The complete steps from Green’s identity to

(1.20) can be found in [Wilcox, 1956].

We can immediately substitute E and H in for x because they satisfy the vector Helmholtz

equation and arrive at

p =Sp ◦ (cid:0)ˆnp × ∇ × Et
Es

p

(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × Et

p

(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Et

p

(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Et

p

(cid:1)

Hs

p =Sp ◦ (cid:0)ˆnp × ∇ × Ht

p

(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Ht

p

(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Ht(cid:1).

(1.22b)

(1.22a)

Because we are using electromagnetic fields and have Maxwell’s equations, we can simplify

by using Gauss’ Laws

p =Sp ◦ (cid:0)ˆnp × ∇ × Et
Es
p =Sp ◦ (cid:0)ˆnp × ∇ × Ht

Hs

p

(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × Et(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Et(cid:1)
(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Ht
(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

p

p

(1.23a)

(1.23b)

(cid:1).

4

It is also possible to use Faraday’s and Ampere’s Laws to create two additional representations

by taking the curl of both sides and substituting for ∇ × E and ∇ × H

p =∇ × Sp ◦ (cid:0)ˆnp × Et
Es

p

(cid:1) +

Hs

p =∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:1) −

ηp
jκp
1
jκpηp

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:1)

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Et

p

(cid:1).

(1.24a)

(1.24b)

We have three representation integrals for each field and will use them to develop several

formulations with different properties.

To introduce the formulation process, we will step through it for the Electric Field Inte-

gral Equation (EFIE), Magnetic Field Integral Equation (MFIE), Combined Field Integral

Equation (CFIE), Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT), and M¨uller formu-

lations. The process for all direct method formulations is similar: use the surface equivalence

principle to break the problem into a homogeneous problem for each region, represent the

field in each region with an integral equation, enforce boundary conditions on the surface,

and solve for the unknowns.

1.1.4 EFIE

We focus on a single Perfect Electrical Conductor (PEC) object immersed in free space

with a plane wave incident upon it. For the incident electric and magnetic field plane wave,

we will utilize

Ei

1 (r) =E0 exp (−jκ1 · r)

Hi

1 (r) =H0 exp (−jκ1 · r)

(1.25a)

(1.25b)

where E0 and H0 = 1
η1

ˆκ×E0 are polarization vectors that are perpendicular to the propagation

direction ˆκ and each other. There is no incident field produced within the interior of the

object Ei

2 = Hi

2 = 0.

We use the surface equivalence principle to break our problem into two problems. For the

first problem, we replace the PEC object with free space, prescribe a null field interior to the

object boundary, and place equivalent currents J1 and M1 on its surface. For the second

5

problem, we replace free space with PEC (such that everything is PEC), the field exterior to

the object boundary is the null field, and J2 and M2 are placed on its surface. This second

problem is already solved because no field is supported inside a PEC, so the interior fields are

also the null field, and consequently J2 = M2 = 0. We use the scattered field representation

integral (1.24a), add the incident field, and take −ˆnp × ˆnp× on both sides

−ˆnp × ˆnp × Et

p = −ˆnp × ˆnp ×

(cid:18)

p + ∇ × Sp ◦ (cid:0)ˆnp × Et
Ei

p

(cid:1) +

ηp
jκp

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:1)

(cid:19)

.

(1.26)

At this point, we can apply the boundary conditions (1.15) and (1.16), recall that Et

2 = Ht

2 =

Ms = 0 because the PEC supports no fields interior to it and no magnetic currents on the

surface, and substitute

Et

1 × ˆn1 =0

−ˆn1 × ˆn1 × Ei

ˆn1 × ˆn1 × ∇ × ∇ × S1 ◦ ˆn1 × Ht

1

ˆn1 × Ht

1 =

1 =Js
η1
jκ1
1
κ2
p

Lp ◦ x =

∇ × ∇ × Sp ◦ x

−ˆn1 × ˆn1 × Ei

1 = − jκ1η1ˆn1 × ˆn1 × L1 ◦ Js.

(1.27a)

(1.27b)

(1.27c)

(1.27d)

(1.27e)

1.1.5 MFIE

We reuse some of the work done for the EFIE to formulate the MFIE. This time we use

(1.24b) to represent the scattered field. Following a similar procedure, we add the incident

field and take ˆnp× on either side

ˆnp × Ht

p = ˆnp ×

(cid:18)

Hi

p + ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:1) −

1
jκpηp

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Et

p

(cid:1)

(cid:19)

.

(1.28)

6

We apply the boundary conditions (1.15) and (1.16), recall that Et

2 = Ht

2 = Ms = 0, and

substitute

Et

1 × ˆn1 =0

ˆn1 × Ht

1 =Js

Js =ˆn1 × Hi

1 + ˆnp × ∇ × S1 ◦ ˆn1 × Ht

1

Kp ◦ x =ˆnp × ∇ × Sp ◦ x

ˆn1 × Hi

1 =Js − K1 ◦ Js.

(1.29a)

(1.29b)

(1.29c)

(1.29d)

(1.29e)

1.1.6 CFIE

The EFIE and MFIE operators each have a null space, meaning there exist corresponding

sources Js that produce no scattered field that can be added to any solution. To eliminate

the null space, we add the two formulations together with a scaling factor for each equation






(cid:19)
1

(cid:18)

α
η0

−ˆn1 × ˆn1 × Ei
ˆn1 × Hi

1

1




 =

(cid:18)

α
η0






(cid:19)
1

−jκ1η1ˆn1 × ˆn1 × L1

I − K1




 ◦ Js.

(1.30)

The purpose of dividing the EFIE formulation by η0 is so that both integral equations and

incident fields are of the same scale, which leads to a better-conditioned system. The purpose

of α is to give the ability to tune the formulation with different relative weights for each

formulation as desired.

1.1.7 PMCHWT

We enter the world of dielectrics with the PMCHWT formulation [Poggio and Miller,

1973]. We will make use of two representation integrals (1.24b) and (1.24a), add the incident

7

fields, and take ˆnp × ˆnp× on (1.24b) and −ˆnp × ˆnp× on (1.24a) to either side

−ˆnp × ˆnp × Et

p = − ˆnp × ˆnp ×

(cid:18)

p + ∇ × Sp ◦ (cid:0)ˆnp × Et
Ei

p

(cid:1) +

ηp
jκp

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:1)

(cid:19)

ˆnp × ˆnp × Ht

p =ˆnp × ˆnp ×

(cid:18)

Hi

p + ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:1) −

1
jκpηp

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Et

p

(1.31a)
(cid:19)
(cid:1)

.

We apply the boundary conditions (1.15) and (1.16), but this time, neither electric nor

magnetic currents are supported on the boundary Js = Ms = 0. However, we can relate

the two traces that satisfy the boundary conditions to each other by creating two common

“source-like” functions ¯J and ¯M

(1.31b)

ˆn1 × Ht

1 = − ˆn2 × Ht

2

Et

1 × ˆn1 = − Et




2 × ˆn2






P 1

1

ˆn1 × Ht
Et

1 × ˆn1


 =P 2




P p = diag

εp,






 =




¯J

¯M

2

ˆn2 × Ht
Et
(cid:18)

2 × ˆn2
(cid:19)

εp
η0






(1.32a)

(1.32b)

(1.32c)

(1.32d)

where εp is + for one region and − for the other (e.g., εp = 3 − 2p for the two region case

and p ∈ {1, 2}) and P p contains the boundary condition relationship between the traces of
the fields on either side of the interface and a scaling parameter such that ¯J and ¯M are of

the same scale. We can use a scaling matrix, W p (along with P p), to combine the equations

8

from each region



Zp =




−jκpηpˆnp × ˆnp × Lp

ˆnp × (I − Kp)

ˆnp × (I − Kp)

jκp
ηp

ˆnp × ˆnp × Lp











ˆnp × Ei
p × ˆnp
ˆnp × ˆnp × Hi

p

2
(cid:88)

p=1

P pW p






ˆnp × Ei
p × ˆnp
ˆnp × ˆnp × Hi

p




 =Zp






p

ˆnp × Ht
Et

p × ˆnp









 =

(cid:32) 2

(cid:88)

p=1

P pW pZp

(cid:1)−1

(cid:0)P p

(cid:33)






¯J

¯M




 .

(1.33a)

(1.33b)

(1.33c)

The ˆnp × I operator can result in degenerate matrices if not correctly discretized. With

the right choice of W p, we can eliminate it, analytically, canceling out the operator, and

preconditioning the system

1.1.8 M¨uller

W p = diag

(cid:19)

, η0

.

(cid:18) 1
η0

(1.34)

The M¨uller formulation [M¨uller, 1969] is the final formulation covered in the introduction.

We employ the same representation integrals (1.24b) and (1.24a), add the incident fields, but

this time take ˆnp× on (1.24b) and −ˆnp× on (1.24a) to either side resulting in

ˆnp × Ht

p =ˆnp ×

(cid:18)

−ˆnp × Et

p = − ˆnp ×

Hi

p + ∇ × Sp ◦ (cid:0)ˆnp × Ht
(cid:18)

p + ∇ × Sp ◦ (cid:0)ˆnp × Et
Ei

p

p

(cid:1) −

(cid:1) +

1
jκpηp
ηp
jκp

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Et

p

(cid:19)
(cid:1)

(1.35a)

∇ × ∇ × Sp ◦ (cid:0)ˆnp × Ht

p

(cid:19)
(cid:1)

.

(1.35b)

As before, we apply the boundary conditions (1.15) and (1.16), eliminate sources Js = Ms = 0,

and relate the two boundary condition traces by creating two common “source-like” functions

9

¯J and ¯M

P 1

ˆn1 × Ht

1 = − ˆn2 × Ht

2

Et

1 × ˆn1 = − Et




2 × ˆn2






1

ˆn1 × Ht
Et

1 × ˆn1


 =P 2









 =




¯J

¯M

2

ˆn2 × Ht
Et
(cid:18)

2 × ˆn2
(cid:19)

εp
η0

P p = diag

εp,

.






(1.36a)

(1.36b)

(1.36c)

(1.36d)

We use a scaling matrix, W p (along with P p), to combine the equations from each region






jκp
ηp

ˆnp × Lp

I − Kp






Zp =




 =Zp




I − Kp

−jκpηpˆnp × Lp




ˆnp × Ht
Et × ˆnp







 =

(cid:32) 2

(cid:88)

p=1

P pW pZp

(cid:1)−1

(cid:0)P p

(cid:33)






¯J

¯M




 .






p

ˆnp × Hi
Ei

p × ˆnp






p

ˆnp × Hi
Ei

p × ˆnp

2
(cid:88)

p=1

P pW p

(1.37a)

(1.37b)

(1.37c)

For the M¨uller formulation, we can make a second-kind formulation by using the linear factors

W p = diag (µrp, ϵrp) .

(1.38)

1.1.9 Need and Outline of this Thesis

Well-conditioned formulations free of low-frequency, dense-mesh, and topology break-

downs have recently been developed. These new well-conditioned surface integral equation

formulations apply to the same problems as current formulations, can be solved with a few

iterations across a broad band of frequencies, and are free from typical breakdowns. Like the

ones presented above, existing formulations suffer from one breakdown or another. For that

reason, I will present several new formulations and the necessary details needed to implement

the operators involved in these formulations on spheres and tessellated objects.

10

Much of what is covered will apply tried and true methods and techniques to new

formulations. Here is an incomplete list of existing tools and techniques that will be enlisted

to take the new formulations to market.

1. A surface equivalence framework to represent scattered fields with surface sources

2. Using a direct method formulation and enforcing boundary conditions to link surface

unknowns between regions

3. Spherical harmonics, Rao-Wilton-Glisson (RWG), and pulse basis functions to represent

unknowns and measure functions

4. Method of Moments to create the system of equations

5. Iterative solvers to find the solution for the unknown sources

6. Using MLFMA to accelerate the matrix-vector multiplications needed for an iterative

solver.

1.2 Organization

This thesis will be organized into three main sections. Chapter 2 will cover the new

formulations and discretization process. Chapter 3 will cover implementing these new

formulations with piecewise basis sets. Chapter 4 will cover the acceleration details needed

to implement the new formulations with MLFMA. Results will be presented throughout as

appropriate.

11

CHAPTER 2

FORMULATIONS

2.1 Introduction

A formulation serves as the foundation for a numerical method. The discretized system of

equations to be solved depends on the equations, first and foremost. Many main formulations

utilized today were introduced decades ago and have been under constant development.

In addition to the formulations introduced in Chapter 1, other formulations include the

single integral equation [Marx, 1984, Glisson, 1984] as well as various combinations of

surface equivalence theorems; see [Yla-Oijala and Taskinen, 2005], [Li et al., 2014], and

references therein for some varieties. The necessity of analyzing composite objects has driven

vast advances in the machinery and techniques necessary to compute scattering accurately.

The primary challenges include dense mesh breakdown [Cools et al., 2009b], low-frequency

breakdown [Vecchi, 1999, Qian and Chew, 2008], and topology breakdown [Cools et al.,

2009a]. These breakdowns find their cause in the formulation due to either catastrophic

cancellation, bad constraints, or improper scaling. When representing and measuring these

integral equations on tessellated representations, they manifest as ill-conditioned and poorly

convergent discrete systems. As a result, there has been a concerted effort to develop well-

conditioned formulations in both the electromagnetics and applied mathematics communities

for a while [Kress and Roach, 1978, Costabel and Stephan, 1988, Wilde, 1987, Costabel,

1991, Dely, 2020].

Indeed, early work recognizing these challenges and efforts toward amelioration date

back four decades [Wilton and Glisson, 1981, Wu et al., 1995, Burton and Kashyap, 1995,

Vecchi, 1999] with the introduction of loop-star and loop-tree decompositions. Developing

preconditioners, formulations, and appropriate basis sets to work around these breakdowns

numerically has been a topic of intense research over the past decade [Andriulli et al., 2008].

These methods rely on the representation integrals (1.24b) and (1.24a). More recently,

research has focused on alternative formulations which use (1.22a), (1.22b), and in general

12

(2.17) and (1.20). The two main classes of new formulations are the Decoupled Field Integral

Equation (DFIE) [Vico et al., 2018] and Decoupled Potential Integral Equation (DPIE)

[Chew, 2014, Liu et al., 2015, Vico et al., 2016, Li et al., 2017]. Other formulations and

techniques to overcome some of the bottlenecks of the classical integral equations have been

the reliance on Debye potentials [Epstein and Greengard, 2010, Fu et al., 2017], projectors

[Andriulli et al., 2008, Adrian et al., 2019], and Dirac formulations [Helsing et al., 2020].

The DPIE is a solution to the transmission problem posed in terms of potentials, whereas

the DFIE utilizes fields directly. It has been shown that the DFIE and DPIE are robust to

breakdowns associated with low frequency, mesh discretization, and topology [Vico et al.,

2018, Vico et al., 2016]. The DPIE, as is usually the case, was first developed to analyze

scattering from PECs [Chew, 2014, Liu et al., 2015, Vico et al., 2016]. A corrected Nystr¨om’s

method implementation was presented in [Vico et al., 2015], and a time domain potential

based integral equation solvers [Roth and Chew, 2018, Roth and Chew, 2020] currently

exist. As an aside, it can be shown that [Liu et al., 2015], while based on potential, is

akin to an Augmented Electric Field Integral Equation (A-EFIE) [Qian and Chew, 2009].

As expected, developing equations for analyzing scattering from dielectric objects is more

complicated. It has its genesis in [Vico et al., 2016, Li et al., 2017, Li et al., 2019]. It is

well-conditioned and not susceptible to non-uniqueness due to resonances (under assumptions

on constitutive parameters [Vico et al., 2018, Helsing et al., 2020]) nor breakdown due to

either low frequencies or topology. All these features were demonstrated using analytic basis

sets on spheres or Nystr¨om’s method on canonical geometries. The analysis using potentials

was first presented by [Vico et al., 2016], who then continued their work by developing the

DFIE formulation for dielectrics using an indirect method [Vico et al., 2018]. It has been

noted that the DFIE is non-unique for specific pairs of constitutive parameters when loss

is present. The formulations presented here are presumed to have the same problem; only

lossless materials are considered. The Dirac formulation has recently been introduced with

similar properties for a broader range of constitutive properties [Helsing et al., 2020]. Since

13

the publication of [Vico et al., 2016, Li et al., 2019, Vico et al., 2018], there has yet to be

transition to scattering analysis on tessellated geometries. The complexity of the operators is

deceptively daunting. In what follows, we present details necessary to implement all of the

operators used by the formulations in this chapter in Chapter 3 on spheres using spherical

harmonics basis functions and on tessellated surfaces using piecewise basis functions and

show that the properties demonstrated using analytic basis sets are preserved.

This chapter uses a general formulation process to describe each formulation presented

concisely. Several DFIE and DPIE formulations are developed for dielectric objects, and the

Local Calder´on-Combined Field Integral Equation (LC-CFIE) formulation for PEC objects

is presented using a Calder´on identity that was encountered while developing the DFIE and

DPIE.

2.2 Common Formulation Framework

It is not difficult to see from looking at Chapter 1 that there is a fairly general way of

writing all of these formulations

bχ
p =Z χ

(cid:16)(cid:0)P χ

(cid:1)−1 xχ(cid:17)

p

Np
(cid:88)

p=1

P χ

p W χ

p bχ

p =

p ◦
(cid:32) Np
(cid:88)

p=1

P χ

p W χ

p Z χ
p

(cid:33)

◦ xχ

(cid:1)−1

(cid:0)P χ

p

(2.1a)

(2.1b)

with χ being a placeholder for a specific formulation, Np being the number of regions,

bχ
p being a collection of incident functions, xχ

p being a collection of source or “source-like”

unknown functions, P χ

p enforcing boundary conditions (for dielectric formulations) and scaling

unknowns and measurements, W χ

p being another scaling matrix, and Z χ

p being the block

operator matrix. While this common formulation framework is not entirely general, it applies

to all the formulations discussed in this thesis. We can further incorporate the sum over p by
creating a stacked P χ, (P χ)◦−1T

, and bχ and block W χ and Z χ

P χW χbχ = P χW χZ χ(P χ)◦−1T

◦ xχ

(2.2)

where xχ is now the number of independent unknowns and A = B◦−1, Aij = B−1

ij denotes the

element-wise inverse.

14

For the remainder of this thesis, we will concern ourselves with only the two region cases

with either a PEC or dielectric closed object immersed in free space. However, it is important

to note that our discussion is mainly extendable to composite and open objects.

We will recast the formulations from the introduction within this framework for complete-

ness. For EFIE, we have

For MFIE, we have

P EFIE
1

=1

W EFIE
1

=

1
η0

Z EFIE
1

= − jκ1η1ˆn1 × ˆn1 × L1

bEFIE
1

= − ˆn1 × ˆn1 × Ei

1

xEFIE =Js.

P MFIE
1

=1

W MFIE
1

=1

Z MFIE
1

=I − K1

bMFIE
1

=ˆn1 × Hi

1

xMFIE =Js.

For CFIE, we scale and combined the EFIE and MFIE formulations

P CFIE
1

=






P EFIE
1

P MFIE
1



T




W CFIE
1

= diag (cid:0)αW EFIE

1

, W MFIE
1

(cid:1)

Z CFIE
1

= diag (cid:0)Z EFIE
1




, Z MFIE
1

(cid:1)

bCFIE
1

=




bEFIE
1

bMFIE
1




xCFIE =Js.

15

(2.3a)

(2.3b)

(2.3c)

(2.3d)

(2.3e)

(2.4a)

(2.4b)

(2.4c)

(2.4d)

(2.4e)

(2.5a)

(2.5b)

(2.5c)

(2.5d)

(2.5e)






(2.6a)

(2.6b)

(2.6c)

(2.6d)

(2.6e)

(2.7a)

(2.7b)

(2.7c)

(2.7d)

(2.7e)

For the PMCHWT, we have

P PMCHWT

p

= diag

W PMCHWT

p

= diag

εp
η0

(cid:19)

(cid:19)

, η0

(cid:18)

εp,
(cid:18) 1
η0

Z PMCHWT

p

=






−jκpηpˆnp × ˆnp × Lp

ˆnp × (I − Kp)

ˆnp × (I − Kp)

jκp
ηp

ˆnp × ˆnp × Lp

bPMCHWT

1

=






ˆnp × Ei
p × ˆnp
ˆnp × ˆnp × Hi

p






xPMCHWT =






¯J

¯M




 .

For M¨uller, we have

P M¨uller
p

(cid:18)

= diag

εp,

(cid:19)

εp
η0

W M¨uller
p

= diag (µrp, ϵrp)



I − Kp

−jκpηpˆnp × Lp






jκp
ηp

ˆnp × Lp

I − Kp

Z M¨uller
p

=

bM¨uller
1

=

xM¨uller =



















p

ˆnp × Hi
Ei

p × ˆnp


¯J

¯M


 .

16

2.3 Operators

Before presenting additional formulations, all of the operators used will be defined. In

total, fifteen different integral operators are defined as follows:





T




x

x




Sp ◦

(cid:90)

=

Gp (r, r′)






x (r′)

x (r′)


T




dS′

Dp ◦ x = − ∇ · Sp ◦ (cid:0)ˆnpx(cid:1)
Np ◦ x = − ˆnp · ∇∇ · Sp ◦ (cid:0)ˆnpx(cid:1)

D′

p ◦ x =ˆnp · ∇Sp ◦ x
p ◦ x = − ∇ × Sp ◦ (cid:0)ˆnp × x(cid:1)
K′
p ◦ x =ˆnp × Sp ◦ (cid:0)ˆnp × x(cid:1)
J (2)
p ◦ x =ˆnp · Sp ◦ (cid:0)ˆnp × x(cid:1)
J (3)
J (4) ◦ x =∇ · Sp ◦ (cid:0)ˆnp × x(cid:1)

Lp ◦ x =

1
κ2
p

∇ × ∇ × Sp ◦ x

Kp ◦ x =ˆnp × ∇ × Sp ◦ x

M(3)

p ◦ x =ˆnp · ∇ × Sp ◦ x

P (2)

Q(1)

p ◦ x =ˆnp × ∇Sp ◦ x
p ◦ x =ˆnp × ˆnp × ∇ × Sp ◦ (cid:0)ˆnpx(cid:1)
p ◦ x =ˆnp × Sp ◦ (cid:0)ˆnpx(cid:1)
p ◦ x =ˆnp · Sp ◦ (cid:0)ˆnpx(cid:1).

Q(2)

Q(3)

(2.8a)

(2.8b)

(2.8c)

(2.8d)

(2.8e)

(2.8f)

(2.8g)

(2.8h)

(2.8i)

(2.8j)

(2.8k)

(2.8l)

(2.8m)

(2.8n)

(2.8o)

We define the Single Layer Potential operator (2.8a) as accepting either a scalar or vector-

valued function, denote adjoint operators with a prime, and introduce the shorthand for

operators Ot

p = ˆnp × ˆnp × Op.

17

2.4 DFIE

The DFIE formulations represent and enforce boundary conditions on fields on either side

of an interface as was done for the PMCHWT and M¨uller formulations. There are two forms

of the DFIE, one for E and one for H that are very closely related to each other.

2.4.1 E-Decoupled Field Integral Equation (E-DFIE)

For the E-DFIE, we employ the representation integral (1.22a), add the incident fields,

and take ˆnp × ˆnp × ∇×, ˆnp×, ˆnp·, ∇· on either side of the equation. But before that, we

chose to add an extra ˆnp× to the first equation and unknown such that the representation

integral reads

Es = − Sp ◦ (cid:0)ˆnp × ˆnp × ˆnp × ∇ × Et

(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × Et(cid:1)
− ∇Sp ◦ (cid:0)ˆnp · Et(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Et

p

p

(2.9)

(cid:1).

This is not strictly required but is done so that Galerkin testing with the piecewise RWG

basis functions can easily be used with both vector unknowns being represented and measured

with RWG functions as will be done in Section 3.3.

ˆnp × ˆnp × ∇ × Et

p =ˆnp × ˆnp × ∇ × (cid:0)Ei − Sp ◦ (cid:0)ˆnp × ˆnp × ˆnp × ∇ × Et

p

(cid:1)

(2.10a)

+∇ × Sp ◦ (cid:0)ˆnp × Et(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Et(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Et

p

(cid:1)(cid:1)

ˆnp × Et

p =ˆnp × (cid:0)Ei − Sp ◦ (cid:0)ˆnp × ˆnp × ˆnp × ∇ × Et

p

(cid:1)

+∇ × Sp ◦ (cid:0)ˆnp × Et(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Et(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Et

p

(cid:1)(cid:1)

ˆnp · Et

p =ˆnp · (cid:0)Ei − Sp ◦ (cid:0)ˆnp × ˆnp × ˆnp × ∇ × Et

p

(cid:1)

+∇ × Sp ◦ (cid:0)ˆnp × Et(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Et(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Et

p

(cid:1)(cid:1)

∇ · Et

p =∇ · (cid:0)Ei − Sp ◦ (cid:0)ˆnp × ˆnp × ˆnp × ∇ × Et

p

(cid:1)

(2.10b)

(2.10c)

(2.10d)

+∇ × Sp ◦ (cid:0)ˆnp × Et(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Et(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Et

p

(cid:1)(cid:1)

We apply the boundary conditions (1.15), (1.16), (1.17), and create a fourth boundary

condition on ∇ · Et

p (such that it is continuous across the boundary). As for PMCHWT

and M¨uller, there are no currents or charges supported on the boundary Js = Ms = 0 and

18

ρs == ρms = 0 but we can create four “source-like” functions, ¯a, ¯b, ¯γ, and ¯σ to use

1
µr1

ˆn1 × ˆn1 × ∇ × Et

1 =

1
µr2

ˆn2 × ˆn2 × ∇ × Et

2

ˆn1 × Et

1 = − ˆn2 × Et

2

ϵr1ˆn1 · Et

1 = − ϵr2ˆn2 · Et

2

∇ · Et
1 =∇ · Et


2












1

ˆn1 × ˆn1 × ∇ × Et
ˆn1 × Et
ˆn1 · Et
∇ · Et
1

1

1










=P EDFIE
2












2

ˆn2 × ˆn2 × ∇ × Et
ˆn2 × Et
ˆn2 · Et
∇ · Et
2

2

2












=























¯a

¯b

¯γ

¯σ

P EDFIE

1

We can define the matrices, operators, and vectors for this formulation as

P EDFIE

p

= diag

, −jκ0εp, −jκ0εpϵrp, 1

(cid:19)

(cid:18) 1
µrp

(cid:18)

(cid:19)

1
ϵrp

,

1
µrp

W EDFIE
p

= diag

µrp, ϵrp,

Z EDFIE

p

=












I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

J (4)
p

I − Kp

−M(3)
p

0












p

p

ˆnp × ˆnp × ∇ × Ei
ˆnp × Ei
ˆnp · Ei
∇ · Ei
p
(cid:19)T

p

¯a ¯b ¯γ ¯σ

.

bEDFIE

p

=

(cid:18)

xEDFIE =












0

P (2)
p

Q(1)
p

Q(2)
p

I + D′

p Q(3)

p

−κ2Sp I − Dp











(2.11a)

(2.11b)

(2.11c)

(2.11d)

(2.11e)

(2.12a)

(2.12b)

(2.12c)

(2.12d)

(2.12e)

The linear factors of W EDFIE

p

are chosen such that the E-DFIE formulation results in a

second-kind integral equation.

19

2.4.2 E-Reduced Decoupled Field Integral Equation (E-RDFIE)

Rather than using (1.22a), we could explicitly enforce Gauss’s law (1.3) by using (1.23a).

A second-kind integral equation can be made by removing the fourth row and column of the

E-DFIE.

By reusing many of the same definitions from the E-DFIE, we can jump straight to

P ERDFIE

p

= diag

, −jκ0εp, −jκ0εpϵrp

(cid:19)

(cid:18) 1
µrp

(cid:18)

W ERDFIE

p

= diag

µrp, ϵrp,

(cid:19)

1
ϵrp

Z ERDFIE

p

=









I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

I − Kp

−M(3)
p









(cid:18)

p

ˆnp × ˆnp × ∇ × Ei
ˆnp × Ei
ˆnp · Ei
(cid:19)T

p

p

¯a ¯b ¯γ

bERDFIE

p

=

xERDFIE =









0

P (2)
p

I + D′
p








(2.13a)

(2.13b)

(2.13c)

(2.13d)

(2.13e)

2.4.3 H-Decoupled Field Integral Equation (H-DFIE)

For the H-DFIE, we employ the same tactics as with the E-DFIE, and because of duality,

we can replace E with H, µ with ϵ, and ϵ with µ and arrive at the formulation. Boundary

conditions are expressed as

P HDFIE

1












1

ˆn1 × ˆn1 × ∇ × Ht
ˆn1 × Ht
ˆn1 · Ht
∇ · Ht
1

1

1












2

ˆn2 × ˆn2 × ∇ × Ht
ˆn2 × Ht
ˆn2 · Ht
∇ · Ht
2

2

2












=












= P HDFIE
2

20












(2.14a)



¯a



¯b







¯γ

¯σ

and the remaining matrices, operators, and vectors are

P HDFIE

p

= diag

, −jκ0εp, −jκ0εpµrp, 1

(cid:19)

(cid:18) 1
ϵrp

(cid:18)

(cid:19)

1
µrp

,

1
ϵrp

W HDFIE
p

= diag

ϵrp, µrp,

Z HDFIE

p

=












I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

J (4)
p

I − Kp

−M(3)
p

0












p

p

ˆnp × ˆnp × ∇ × Hi
ˆnp × Hi
ˆnp · Hi
∇ · Hi
p
(cid:19)T

p

¯a ¯b ¯γ ¯σ

bHDFIE

p

=

(cid:18)

xHDFIE =

0

P (2)
p

Q(1)
p

Q(2)
p

I + D′

p Q(3)

p

−κ2Sp I − Dp






















(2.14b)

(2.14c)

(2.14d)

(2.14e)

(2.14f)

21

2.4.4 H-Reduced Decoupled Field Integral Equation (H-RDFIE)

Similarly, we could enforce Gauss’s law (1.4) explicitly by using (1.23b). A second-kind

integral equation can be made by removing the fourth row and column of the H-DFIE.

P HRDFIE

p

= diag

, −jκ0εp, −jκ0εpµrp

(cid:19)

(cid:18) 1
ϵrp

(cid:18)









0

P (2)
p

I + D′
p


W HRDFIE

p

= diag

ϵrp, µrp,

(cid:19)

1
µrp

















(cid:18)

Z HRDFIE

p

=

bHRDFIE

p

=

xHRDFIE =

I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

I − Kp

−M(3)
p

p

ˆnp × ˆnp × ∇ × Hi
ˆnp × Hi
ˆnp · Hi
(cid:19)T

p

p







¯a ¯b ¯γ

(2.15a)

(2.15b)

(2.15c)

(2.15d)

(2.15e)

2.5 DPIE

All other formulations discussed thus far have used representations of fields, either E or

H, as a starting point. The novelty of the DPIE formulation is that it uses potentials.

It is well-known that fields can be represented with potentials. For the vector magnetic

and scalar electric potential A − ϕ decomposition and the vector electric and scalar magnetic

potential F − ψ decomposition, we have

E = − jωA − ∇ϕ

H =

1
µ

∇ × A

H = − jωF − ∇ψ

E =

−1
ϵ

∇ × F.

22

(2.16a)

(2.16b)

(2.16c)

(2.16d)

Furthermore, under the Lorenz gauge,

∇ · A =

∇ · F =

κ2
jω
κ2
jω

ϕ

ψ

these potentials satisfy the Helmholtz equation

∇2A + κ2A =0

∇2ϕ + κ2ϕ =0

∇2F + κ2F =0

∇2ψ + κ2ψ =0.

(2.16e)

(2.16f)

(2.16g)

(2.16h)

(2.16i)

(2.16j)

We start diverging from the previous formulations by representing the potentials A and F

rather than the fields with (1.20). The scalar counterpart to Green’s Vector Representation

Theorem is derived from Green’s scalar identities. With x being any scalar that satisfies the

scalar Helmholtz equation ∇2x + κ2x = 0, we have

p = −∇ · Sp ◦ (cid:0)ˆnpxt
xs

p

(cid:1) − Sp ◦ (cid:0)ˆnp · ∇xt

p

(cid:1).

(2.17)

We use Green’s Scalar Representation Theorem to represent ϕ and ψ.

Explicitly, we have four representation equations

As

p =Sp ◦ (cid:0)ˆnp × ∇ × At

p

(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × At

p

(cid:1) − ∇Sp ◦ (cid:0)ˆnp · At

p

(cid:1) − Sp ◦ (cid:0)ˆnp∇ · At

p

(cid:1)

p = − ∇ · Sp ◦ (cid:0)ˆnpϕt
ϕs
p =Sp ◦ (cid:0)ˆnp × ∇ × Ft
Fs

(cid:1) − Sp ◦ (cid:0)ˆnp · ∇ϕt
(cid:1) + ∇ × Sp ◦ (cid:0)ˆnp × Ft

(cid:1)

p

p

p

p

p = − ∇ · Sp ◦ (cid:0)ˆnpψt
ψs

p

(cid:1) − Sp ◦ (cid:0)ˆnp · ∇ψt

p

(cid:1).

(2.18a)

(2.18b)

(cid:1) − ∇Sp ◦ (cid:0)ˆnp · Ft

p

(cid:1) − Sp ◦ (cid:0)ˆnp∇ · Ft

p

(cid:1)

(2.18c)

(2.18d)

At this point, we can impose boundary conditions; however, the boundary conditions we

have available are on the fields, not potentials. Because there is a many-to-one relationship

23

from potentials to fields, many options exist for deriving potential-based boundary conditions

from the field boundary conditions. One method in particular results in decoupling the vector

and scalar potentials [Li et al., 2019].

As for the DFIE, there are several variants to the DPIE, and a collection of them will be

covered in the following sections.

2.5.1 Incident Wave

The incident waves ϕi

1 and Ai

1 decompose Ei

1 (1.25a) while ψi

1 and Fi

1 decompose Hi

1

(1.25b) and are in the form found in [Vico et al., 2015]:

ϕi
1 = − (r · E0) exp (−jκ1 · r)

Ai

1 = −

κ1
ω

(r · E0) exp (−jκ1 · r)

ψi

1 = − (r · H0) exp (−jκ1 · r)

Fi

1 = −

κ1
ω

(r · H0) exp (−jκ1 · r) .

Alternate decompositions, such as

ϕi
1 =ϕ0 exp (−jκ1 · r)

Ai

1 =A0 exp (−jκ1 · r)

A0 =

−1
jω

(E0 − jκ1ϕ0)

ψi

1 =ψ0 exp (−jκ1 · r)

Fi

1 =F0 exp (−jκ1 · r)

F0 =

−1
jω

(H0 − jκ1ψ0) .

(2.19a)

(2.19b)

(2.19c)

(2.19d)

(2.20a)

(2.20b)

(2.20c)

(2.20d)

(2.20e)

(2.20f)

with ϕ0 and ψ0 being arbitrary constants, have limitations at low frequencies but can be

useful for analysis with spherical harmonics as done in Section 3.2.4.2.

2.5.2 A − ϕ-Decoupled Potential Integral Equation (A − ϕ-DPIE)

We use boundary conditions derived from the boundary conditions on fields (1.15), (1.16),

(1.17), and (1.18) that decouples A from ϕ. We create six “source-like” functions, ¯a, ¯b, ¯γ, ¯σ,

24

(2.21a)

(2.21b)

(2.21c)

(2.21d)

(2.21e)

(2.21f)

(2.21g)





































=



¯a

¯γ



¯b














¯σ

¯α

¯β

¯α, and ¯β to use

1
µr1

ˆn1 × ˆn1 × ∇ × At

1 =

1
µr2

ˆn2 × ˆn2 × ∇ × At

2

1 = − ˆn2 × At

2

ˆn1 · At

ˆn1 × At
1
ϵr1
1
µr1ϵr1

∇ · At

1 =

1 =

−1
ϵr2
1
µr2ϵr2

ˆn2 · At

2

∇ · At
2

ϕ1 =ϕ2

ϵr1ˆn1 · ∇ϕ1 = − ϵr2ˆn2 · ∇ϕ2



















2

ˆn2 × ˆn2 × ∇ × At
ˆn2 × At
ˆn2 · At
∇ · At
2

2

2

ϕ2

ˆn2 · ∇ϕ2

P A−ϕDPIE

1



















1

ˆn1 × ˆn1 × ∇ × At
ˆn1 × At
ˆn1 · At
∇ · At
1

1

1

ϕ1

ˆn1 · ∇ϕ1



















=P A−ϕDPIE

2

25

We can define the matrices, operators, and vectors for this formulation as

P A−ϕDPIE

p

= diag

, −jκ0εp, −jκ0εpϵrp,

(cid:19)

1
n2
p

,

−jκ0
c0

,

εpϵrp
c0

(cid:18) 1
µrp

(cid:18)

1
ϵrp

, ϵrp, 1,

(cid:19)

1
ϵrp

W A−ϕDPIE

p

= diag

µrp, ϵrp,

Z A−ϕDPIE

p

=



















I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

J (4)
p

0

0

I − Kp

−M(5)
p

0

0

0

bA−ϕDPIE

p

=



















p

ˆnp × ˆnp × ∇ × Ai
ˆnp × Ai
ˆnp · Ai
∇ · Ai
p

p

p

ϕi
p

ˆnp · ∇ϕi
p

0

P (2)
p

Q(1)
p

Q(2)
p

I + D′

p Q(3)

p

−κ2

pSp I − Dp

0

0

0

0

0

0

0

0

0

0

0

0

I − Dp

Sp

−Np

I + D′
p



















xA−ϕDPIE =

(cid:18)

¯a ¯b ¯γ ¯σ ¯α ¯β

(cid:19)T

.



















(2.22a)

(2.22b)

(2.22c)

(2.22d)

(2.22e)

The linear factors of W p are chosen such that the DPIE formulation only has one operator

whose strong singularity is not canceled.

2.5.3 A − ϕ-Reduced Decoupled Potential Integral Equation (A − ϕ-RDPIE)

Alternate formulations are possible due to the Lorenz gauge, which relates two unknowns

to each other ¯σ = ¯α. We can immediately create several Reduced-Decoupled Potential

Integral Equation (R-DPIE) formulations with the six equations and five unknowns and show

26

two as examples. The first is created by eliminating ¯σ

P A−ϕRDPIE 1

p

= diag

, −jκ0εp, −jκ0εpϵrp,

(cid:18) 1
µrp

(cid:18)

1
ϵrp

, 1,

(cid:19)

1
ϵrp

(cid:19)

−jκ0
c0

,

εpϵrp
c0

p

p

κ2
jω Q(1)
p
κ2
jω Q(2)
p
κ2
jω Q(3)
p
I − Dp

p

0

0

0

Sp

−Np

I + D′
p

0

P (2)
p

I + D′
p

0

0





























(2.23a)

(2.23b)

(2.23c)

(2.23d)

(2.23e)

W A−ϕRDPIE 1

p

= diag

µrp, ϵrp,

Z A−ϕRDPIE 1

p

=















I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

0

0

I − Kp

−M(5)
p

0

0

bA−ϕRDPIE 1

p

=















p

ˆnp × ˆnp × ∇ × Ai
ˆnp × Ai
ˆnp · Ai

p

p

ϕi
p

ˆnp · ∇ϕi
p

xA−ϕRDPIE 1 =

(cid:18)

¯a ¯b ¯γ

¯α ¯β

(cid:19)T

.

27

The second is create by eliminating ¯α

P A−ϕRDPIE 2

p

= diag

(cid:18) 1
µrp

(cid:18)

W p = diag

µrp, ϵrp,

, −jκ0εp, −jκ0εpϵrp,

1
ϵrp

, ϵrp,

(cid:19)

1
ϵrp

0

P (2)
p

(cid:19)

1
n2
p

,

εpϵrp
c0

Q(1)
p

Q(2)
p

0

0

0

0

I + D′

p Q(3)

p

−κ2

pSp I − Dp

0

−jω
κ2
p

Np I + D′
p





























(2.24a)

(2.24b)

(2.24c)

(2.24d)

(2.24e)

Z A−ϕRDPIE 2

p

=















I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

J (4)
p

0

I − Kp

−M(5)
p

0

0

bA−ϕRDPIE 2 =















p

ˆnp × ˆnp × ∇ × Ai
ˆnp × Ai
ˆnp · Ai
∇ · Ai
p

p

p

ˆnp · ∇ϕi
p

xA−ϕRDPIE 2 =

(cid:18)

¯a ¯b ¯γ ¯σ ¯β

(cid:19)T

.

Additional A − ϕ-RDPIEs not shown here result from scaling and combining the fourth

and fifth rows of the A − ϕ-DPIE after combining the fourth and fifth unknown functions.

2.5.4 F − ψ-Decoupled Potential Integral Equation (F − ψ-DPIE)

We can make use of the duality between the F − ψ-DPIE and A − ϕ-DPIE formulations

by replacing A with F, ϕ with ψ, µ with ϵ, and ϵ with µ and arrive at the formulation.

28

Boundary conditions are expressed as

P F−ψDPIE

1



















1

ˆn1 × ˆn1 × ∇ × Ft
ˆn1 × Ft
ˆn1 · Ft
∇ · Ft
1

1

1

ψ1

ˆn1 · ∇ψ1



















= P F−ψDPIE
2



















2

ˆn2 × ˆn2 × ∇ × Ft
ˆn2 × Ft
ˆn2 · Ft
∇ · Ft
2

2

2

ψ2

ˆn2 · ∇ψ2



















=





































¯a

¯b

¯γ

¯σ

¯α

¯β

and the remaining matrices, operators, and vectors are

P F−ψDPIE

p

= diag

, −jκ0εp, −jκ0εpµrp,

(cid:19)

,

−jκ0
c0

,

εpµrp
c0

1
n2
p

(cid:18) 1
ϵrp

(cid:18)

(cid:19)

1
µrp

, µrp, 1,

1
µrp

W F−ψDPIE

p

= diag

ϵrp, µrp,

Z F−ψDPIE

p

=



















I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

J (4)
p

0

0

I − Kp

−M(5)
p

0

0

0

bF−ψDPIE

p

=



















p

ˆnp × ˆnp × ∇ × Fi
ˆnp × Fi
ˆnp · Fi
∇ · Fi
p

p

p

ψi
p

ˆnp · ∇ψi
p



















0

P (2)
p

Q(1)
p

Q(2)
p

I + D′

p Q(3)

p

−κ2

pSp I − Dp

0

0

0

0

0

0

0

0

0

0

0

0

I − Dp

Sp

−Np

I + D′
p



















(cid:18)

xF−ψDPIE =

¯a ¯b ¯γ ¯σ ¯α ¯β

(cid:19)T

.

29

(2.25a)

(2.25b)

(2.25c)

(2.25d)

(2.25e)

(2.25f)

2.5.5 F − ψ-Reduced Decoupled Potential Integral Equation (F − ψ-RDPIE)

The same observation for the A − ϕ-RDPIE is applicable for the F − ψ-RDPIE, and

alternate formulations are possible due to the Lorenz gauge because ¯σ = ¯α. The first of the

two examples are created by eliminating ¯σ

P F−ψRDPIE 1

p

= diag

, −jκ0εp, −jκ0εpµrp,

(cid:19)

−jκ0
c0

,

εpµrp
c0

(cid:18) 1
ϵrp

(cid:18)

(cid:19)

1
µrp

, 1,

1
µrp

p

p

κ2
jω Q(1)
p
κ2
jω Q(2)
p
κ2
jω Q(3)
p
I − Dp

p

0

0

0

Sp

−Np

I + D′
p

0

P (2)
p

I + D′
p

0

0





























(2.26a)

(2.26b)

(2.26c)

(2.26d)

(2.26e)

W F−ψRDPIE 1

p

= diag

ϵrp, µrp,

Z F−ψRDPIE 1

p

=















I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

0

0

I − Kp

−M(5)
p

0

0

bF−ψRDPIE 1

p

=















p

ˆnp × ˆnp × ∇ × Fi
ˆnp × Fi
ˆnp · Fi

p

p

ψi
p

ˆnp · ∇ψi
p

xF−ψRDPIE 1 =

(cid:18)

¯a ¯b ¯γ

(cid:19)T

.

¯α ¯β

30

The second is create by eliminating ¯α

P F−ψRDPIE 2

p

= diag

(cid:18) 1
ϵrp

(cid:18)

W p = diag

ϵrp, µrp,

1
µrp

, µrp,

, −jκ0εp, −jκ0εpµrp,

(cid:19)

1
n2
p

,

εpµrp
c0

Q(1)
p

Q(2)
p

(cid:19)

1
µrp

0

P (2)
p

I + D′

p Q(3)

p

−κ2

pSp I − Dp

0

0

0

0

0

−jω
κ2
p

Np I + D′
p





























(2.27a)

(2.27b)

(2.27c)

(2.27d)

(2.27e)

Z F−ψRDPIE 2

p

=















I − K′t

p −κ2

pLt
p

J (2)
p

J (3)
p

J (4)
p

0

I − Kp

−M(5)
p

0

0

bF−ψRDPIE 2 =















p

ˆnp × ˆnp × ∇ × Fi
ˆnp × Fi
ˆnp · Fi
∇ · Fi
p

p

p

ˆnp · ∇ψi
p

xF−ψRDPIE 2 =

(cid:18)

¯a ¯b ¯γ ¯σ ¯β

(cid:19)T

.

Again, combining the fourth and fifth rows of F − ψ-DPIE, we can create further F − ψ-

RDPIEs formulations.

2.6 LC-CFIE

Thus far, the new formulations have been for dielectric objects. Indeed, my research goal

was to advance dielectric formulations, and a subtle observation led to a new preconditioned

formulation for PEC objects.

The system matrices that we have been working with are all Calder´on type projectors.

Specifically, we can view the system matrices Zp of some of the dielectric formulations (M¨uller

(2.7c), DFIE (2.12c) (2.14d), DPIE (2.22c) (2.25d)) as projecting the traces of total fields or

potentials on the surface to the incident traces on the surface. The complementary projector

I − Zp projects the traces of the total fields or potentials on the surface to the scattered

31

traces on the surface. From the perspective of projectors, it is obvious that the result of

applying a projector and its complementary projector will be a block null operator matrix.

Zp ◦ (cid:0)I − Zp

(cid:1) = 0.

(2.28)

Each entry in the operator matrix above will be the null operator, and Calder´on type identities

can be derived. One of these Calder´on identities that results from using any of the new Z χ
p

in this chapter can be used to precondition the EFIE formulation

p ◦ (cid:0)κ2
J (2)

pLp

(cid:1) = Kp ◦ (I − Kp) − P (2)

p ◦ M(3)
p .

(2.29)

With this identity, we can arrive at a Calder´on preconditioned EFIE formulation

1 ◦ (cid:0)−jκ1ˆn1 × ˆn1 × Ei
J (2)

1

(cid:1) = J (2)

1 ◦ (cid:0)κ2

1η1Lt

1 ◦ Js

(cid:1).

(2.30)

This formulation can be combined with the MFIE to create a CFIE formulation.






(cid:19)
1

(cid:18)

α
η0

1 ◦ (cid:0)−jκ1ˆn1 × ˆn1 × Ei
J (2)
ˆn1 × Hi

1

1


(cid:1)

 =

(cid:18)

α
η0

(cid:19)
1






J (2)

1 ◦ (κ2

1η1Lt
1)

I − K1




 ◦ Js.

(2.31)

This formulation has a null space because the MFIE and this Calder´on EFIE share a null

space. Fortunately, the stabilizing properties of J (2)

p

are local, and we can use a lossy ˜κ1

which not only shifts the null space but also improves high-frequency stability

p ◦ x =ˆnp × ˜Sp ◦ (ˆn × x)
˜J (2)

(cid:90)

˜Sp ◦ x =

˜Gp (r, r′) x (r′) dS′

˜Gp (r, r′) =

exp (−j˜κp |r − r′|)
4π |r − r′|

.

arriving at






(cid:19)
1

(cid:18)

α
η0

1 ◦ (cid:0)−jκ1ˆn1 × ˆn1 × Ei
˜J (2)
ˆn1 × Hi

1

1


(cid:1)

 =

(cid:18)

α
η0






(cid:19)
1

˜J (2)

1 ◦ (κ2

1η1Lt
1)

I − K1




 ◦ Js

(2.32a)

(2.32b)

(2.32c)

(2.32d)

which is unique and stable across a broad range of frequencies with the proper choice of

˜κ = κ − j0.4H 2/3κ1/3 [Antoine et al., 2006] with H being the max mean curvature of the

object.

32

Writing this in the common formulation framework results in

P LC-CFIE

1

=

(cid:18)

(cid:18)

(cid:19)

1 1

W 1 =

α
η0

Z LC-CFIE

1

= diag


bLC-CFIE

1

=




(cid:19)
1
(cid:16) ˜J (2)

1 ◦ (cid:0)κ2

1η1Lt
1

(cid:17)

(cid:1), I − K1


1 ◦ (cid:0)−jκ1ˆn1 × ˆn1 × Ei
˜J (2)
ˆn1 × Hi

1

1

(cid:1)




xLC-CFIE =Js.

2.7 Conclusion

(2.33a)

(2.33b)

(2.33c)

(2.33d)

(2.33e)

The several formulations presented in this chapter cover most used in CEM. The direct

method framework developed here is general and applicable to most formulations. The

following chapter will build on this framework by discretizing these formulations.

33

CHAPTER 3

DISCRETIZATION

3.1 Introduction

The integral equations presented thus far are agnostic to the basis functions used to

discretize these systems. For the formulations in Chapter 1, only vector unknowns need to be

represented, but in formulations presented in Chapter 2, both vector and scalar unknowns are

needed. We introduce some additional notation for ease of presentation (and implementation).

Let the space of scalar basis functions be denoted using Bs

n, for n ∈ [1, Ns]. Likewise, the

space of vector basis functions is denoted using Bv

n, for n ∈ [1, Nv]. In the above, we note that

Ns and Nv are the number of scalar and vector basis functions. For the single, closed object

scatterer, the number of degrees of freedom will be the Nb = nsNs + nvNv where ns and nv

are the numbers of scalar and vector unknown source functions respectively. It follows that

one can represent the collection of basis functions needed to represent all of the unknowns

using F χ = diag (Bs, . . . , Bv, . . . ), where “diag” here is used to mean a block-diagonal matrix
and a Bs for each of the ns scalar unknown functions and a Bv for each of the nv vector

unknown functions.

The discretized system for all formulations in this thesis is constructed through a Galerkin

framework using inner products defined as

(cid:90)

(cid:90)

⟨g, f ⟩ =

⟨g, f ⟩ =

g∗ (r) f (r) dS

g∗ (r) · f (r) dS

(3.1a)

(3.1b)

with ∗ indicating complex conjugate and g, f , g, and f denoting arbitrary scalar and vector

functions of r. For each region, we test (2.1b) over the limiting surface as r approaches the

interface from within Ωp and add the systems together, resulting in a single system, coupled

through the boundary conditions, that can be written as

Zχy = bχ

34

(3.2a)

where the elements of this system are defined as

(cid:42)

bχ
k =

F χ
k ,

Np
(cid:88)

P χ

p W χ

p bχ

p

(cid:43)

(cid:42)

Zχ

kn =

F χ
k ,

p=1
(cid:32) Np
(cid:88)

p=1

P χ

p W χ

p Z χ
p

(cid:1)−1

(cid:0)P χ

p

(cid:33)

(cid:43)

◦ F χ
n

.

(3.2b)

(3.2c)

As the boundary conditions and scale factors manifest themselves as diagonal matrices, it is

trivial to show that

bχ =

Zχ =

2
(cid:88)

p=1

2
(cid:88)

p=1

¯P χ
p

¯W χ
p

¯P χ
p

¯W χ
p

(cid:10)F χ

p , bχ

p

(cid:11)

(cid:10) ¯F χ

p , Z χ

p ◦ ¯F χ

p

(cid:11) (cid:0) ¯P χ

p

(cid:1)−1

(3.3a)

(3.3b)

where for the two region case ¯W χ

p and ¯P χ

p are built with either Ns or Nv repeated elements

and ¯F χ

p is used to represent an individual region’s unknown source and “source-like” functions.
p and ¯P χ

p we test and measure each region individually and with ¯W χ

Effectively, with ¯F χ

p

we scale and combine inner-products and enforce boundary conditions and scaling for the

unknown source and “source-like” functions on their coefficients. As an aside for composite

objects, there is an added step of selecting elements and functions only related to the particular

region but conceptually the same.

All that is now required is to select the basis sets for Bs and Bv, which will be covered

in the remainder of this chapter. For analytic analysis on a sphere, we will use spherical

harmonics in Section 3.2. For numerical analysis on arbitrarily shaped tessellated objects, we

will use piecewise RWG functions in Section 3.3.

3.2 Analytic Analysis

The functions used to extract characteristics for all formulations will be vector and scalar

spherical harmonics basis and testing functions. These functions will represent the source

and “source-like” unknowns and test the integral equations on the surface of a spherical

scatterer with radius a. Because of the spherical harmonic expansion of Green’s function

and orthogonality properties of the spherical harmonics, we can compute all of the integrals

35

analytically. This is as in [Li et al., 2019, Yuan and Shanker, 2005, Hsiao and Kleinman,

1997].

We will use spherical harmonic indices l and m and understand that special care must

be taken when l = 0 so that normalization terms do not introduce any division by 0. For

simplicity of presentation, we will ignore this detail and set a = 1 m.

3.2.1 Scalar Spherical Harmonics for Bs

To construct a basis for the space of scalar functions on the surface of a sphere, we utilize

the orthonormalized spherical harmonics defined as

Ym

l (ˆr) =

(cid:115)

2l + 1
4π

(l − m)!
(l + m)!

Pm

l (cos θ) exp (jmϕ)

(3.4a)

with Pm

l (cos θ) denoting associated Legendre polynomials, l ≥ 0, and |m| ≤ l.

With these basis functions, we can define Bs as

(cid:18)

Bs =

Y0
0

. . . Ym
l

. . . YNh
Nh

(cid:19)

(3.4b)

with Nh being the highest order of harmonic examined and Ns = (Nh + 1)2. For the remainder

of this chapter, l = Nh = 2 ⌈maxp ℜ (κpa)⌉ + 2 where ℜ (κpa) takes the real component of κpa

and maxp takes the max over the single exterior region for the PEC case or the two regions

for the dielectric case.

3.2.2 Vector Spherical Harmonics for Bv

For the space of vector functions, we utilize orthonormalized vector spherical harmonics

defined with respect to the scalar spherical harmonics

Y m

l (ˆr) =ˆr Ym

l (θ, ϕ)

Ψm

l (ˆr) = − ˆr × Φm

l (ˆr) =

Φm

l (ˆr) =ˆr × Ψm

l (ˆr) =

r∇ Ym
l (ˆr)
(cid:112)l (l + 1)
r × ∇ Ym
l (ˆr)
(cid:112)l (l + 1)

(3.5a)

(3.5b)

(3.5c)

The space of functions we want to represent is tangential to the surface. Because of this, we

only require Ψ and Φ type functions to represent them. With these basis functions, we can

36

define Bv as

(cid:18)

(cid:18)

(cid:18)

BΨ =

BΨ =

Bv =

(cid:19)

Ψ0
0

. . . Ψm
l

. . . ΨNh
Nh

Φ0
0

. . . Φm
l

. . . ΦNh
Nh

(cid:19)

(cid:19)

BΨ BΦ

.

(3.6a)

(3.6b)

(3.6c)

with Nv = 2 (Nh + 1)2.

3.2.3 Zero-Mean Constraint

It is important to note that for some unknown functions, there are additional constraints.

For example, when representing γp = ˆnp · Et

p on the surface of a closed region p, it must

satisfy a zero-mean constraint

(cid:90)

γp (r) dS = 0.

(3.7)

When using spherical harmonics to represent these unknowns, this can be performed by

zeroing the coefficient associated with l = m = 0.

3.2.4 Analytic Integration

Evaluation of all of the operators in the formulations in Chapter 2 can be performed

exactly with the spherical harmonic basis sets by exploiting their orthogonality properties.

Rather than stepping through each operator individually, the necessary tools will be given to

perform all of the integrals, and only a single operator will be shown in complete detail. At

the end of this section, all operators will be listed for completeness.

3.2.4.1 Additional Spherical Functions

Spherical Bessel functions are denoted using

b(α)
l

(z) =






jl (z)

α = 1

yl (z)

α = 2

(z) α = 3

(z) α = 4

h(1)
l

h(2)
l

37

(3.8)

where the index α is used to denote the kind of spherical Bessel function.

Using spherical Bessel and harmonic functions, the scalar and vector spherical wave

functions that satisfy the wave equation are

lm (κ, r) = b(α)
φ(α)

l

(κr)

Ym
l (ˆr)
(cid:112)l (l + 1)

L(α)
lm (κ, r) =∇ φ(α)

lm (κ, r)

M(α)

lm (κ, r) =

N(α)

lm (κ, r) =

1
κ
1
κ

∇ × N(α)

lm (κ, r) = −r × ∇ φ(α)

lm (κ, r)

∇ × M(α)

lm (κ, r)

(3.9a)

(3.9b)

(3.9c)

(3.9d)

Note, in typical CEM spherical harmonic analysis, E and H are represented with M(α)

lm and

N(α)

lm only and L(α)

lm is not required because the fields do not radiate radially. However, that

is not the case for representing the vector potentials A and F [Stratton, 2015].

3.2.4.2 Plane Wave Expansion

For completeness, we include plane wave expansion in terms of spherical wave functions.

For ˆκ = ˆz, E0 = ˆx and ϕ0 = 1, we have ϕi

1, Ei

1, and Ai

(3.10a)

(3.10b)

(3.10c)

(3.10d)

ϕi
1 (r) = exp (−jκ1ˆz · r) =

∞
(cid:88)

l
(cid:88)

l=0

m=−l

alm φ(1)

lm (κ1, r)

Ei

1 (r) =ˆx exp (−jκ1ˆz · r) =

∞
(cid:88)

l
(cid:88)

l=0

m=−l

blm M(1)

lm (κ1, r) + clm N(1)

lm (κ1, r)

Ai

1 (r) =

=

−1
jω

−1
jω

(cid:0)Ei

1 + ∇ϕi
1

(cid:1)

∞
(cid:88)

l
(cid:88)

l=0

m=−l

alm L(1)

lm (κ1, r) + blm M(1)

lm (κ1, r) + clm N(1)

lm (κ1, r) .

38

For ˆκ = ˆz, H0 = 1
η1

ˆy and ψ0 = j
η1

, we have ψi

1, Hi

1, and Fi

1

ψi

1 (r) =

Hi

1 (r) =

=

j
η1

1
η1

j
η1

j
η1

∞
(cid:88)

l
(cid:88)

l=0

m=−l

alm φ(1)

lm (κ1, r)

exp (−jκ1ˆz · r) =

ˆy exp (−jκ1ˆz · r)

∞
(cid:88)

l
(cid:88)

l=0

m=−l

clm M(1)

lm (κ1, r) + blm N(1)

lm (κ1, r)

Fi

1 (r) =

−1
jω

(cid:0)Hi

1 + ∇ψi
1

(cid:1)

=

−1
ωη1

∞
(cid:88)

l
(cid:88)

l=0

m=−l

alm L(1)

lm (κ1, r) + clm M(1)

lm (κ1, r) + blm N(1)

lm (κ1, r) .

Defining in this way allows the coefficients alm, blm, and clm to be shared



j−l(cid:113) 4π(2l+1)

l(l+1) m = 0

alm =

blm =

clm =


0



j−l+|m|

√

π(2l+1)
l(l+1)






0

j−(l+m)

√

π(2l+1)
l(l+1)

0

m ̸= 0

|m| = 1

|m| ̸= 1

|m| = 1

.

|m| ̸= 1

(3.10e)

(3.10f)

(3.10g)

(3.10h)

(3.10i)

(3.11a)

(3.11b)

(3.11c)

3.2.4.3 Green’s Function Expansion

The final ingredients necessary for the analytic evaluation of the integrals in all of the

operators used in this thesis are the expansions of Green’s functions. The scalar and dyadic

Green’s function expansions are

Gp (r, r′) = − jκp

∞
(cid:88)

l=0

l (l + 1)

l
(cid:88)

m=−l

lm (κp, r) φ(βp)∗
φ(αp)

lm

(cid:0)κ∗

p, r′(cid:1)

G
p

(r, r′) = − jκp

∞
(cid:88)

l
(cid:88)

l=0

m=−l

M(αp)

lm (κp, r) M(βp)∗

lm

(cid:0)κ∗

p, r′(cid:1)

(3.12a)

(3.12b)

+ N(αp)

lm (κp, r) N(βp)∗

lm

(cid:0)κ∗

p, r′(cid:1)

39

with

αp =

βp =









4 p = 1

1 p = 2

1 p = 1

.

3 p = 2

(3.13a)

(3.13b)

These expressions and the spherical harmonic and wave functions constitute the underpinnings

of our analytical analysis.

3.2.4.4 Candidate Inner Product Evaluation

As there are many different inner products to be evaluated, it is cumbersome to provide a

detailed prescription for each. Indeed, there are 17 unique operators for all of the formulations

presented in Chapters 1 and 2. As most of these follow the same template, we provide one

illustrative example while leaving the reader with the necessary tools to follow the same steps

for the others.

Our example is an operator that is less commonly encountered, P (2)

p . For illustration, we

will focus on a single region and omit the scaling matrices, as these are diagonal and do not

add to the integral. Specifically,

P(2)
p

(cid:12)
(cid:12)kn

= (cid:10)Bv

k, P (2)

p ◦ Bs
n

(cid:11) = (cid:10)Bv

k, ˆnp × ∇S ◦ Bs
n

(cid:11)

(3.14)

where k ∈ [1, Nv] and n ∈ [1, Ns] are the testing and basis indices for this matrix.

Substituting (3.12a) in leads to

(cid:90)

(cid:90)

P(2)
p

(cid:12)
(cid:12)kn

=

=

Bv

k (r) · ˆnp × ∇

Bv

k (r) · ˆnp × ∇

(cid:90)

(cid:90)

G (r, r′) Bs

n (r′) dS′dS

−jκp

∞
(cid:88)

l=0

l (l + 1)

l
(cid:88)

m=−l

lm (κp, r) φ(βp)∗
φ(αp)

lm

(3.15a)

(cid:0)κ∗

p, r′(cid:1) Bs

n (r′) dS′dS

We use (3.4b) and (3.6c) for analytic analysis. If we “unroll” the testing index k into

corresponding harmonic indices l and m for Ψ when k ≤ (Nh + 1)2 and for Φ when k >

(3.15b)

40

(Nh + 1)2, and “unroll” basis index n into corresponding to harmonic indices l′ and m′, it

is easy to see that we can take advantage of the orthogonality of the spherical harmonics.

When l ̸= l′ or m ̸= m′, the system element is zero. This allows us to focus on a specific

harmonic specified by l and m. We will use superscripts Y, Ψ, and Φ to denote the basis

and testing functions as in |Ψ Y

lm to denote testing function Ψm

l and basis function Ym

l . We

analyze the two cases when k selects a testing function from BΨ

P(2)
p

(cid:12)
(cid:12)

Ψ Y

lm

= (cid:10)Ψm
(cid:90)

l , ˆnp × ∇S ◦ Ym
l
(cid:90)

(cid:11)

=

Ψm

l (ˆr) · ˆnp × ∇

−jκpl (l + 1) φ(αp)

lm (κp, r) φ(βp)∗

lm

(cid:0)κ∗

p, r′(cid:1) Ym

l (ˆr′) dS′dS

(3.16a)

= − εpjκpl (l + 1)

(cid:90)

=0

Ψm

l (ˆr) × ˆr · L(αp)

lm (κp, r) dS

(cid:68)

Ym

l , φ(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

and BΦ after skipping some shared steps

P(2)
p

(cid:12)
(cid:12)

Φ Y

lm

= (cid:10)Φm

l , ˆnp × ∇S ◦ Ym
l

(cid:11)

= − εpjκpl (l + 1)

= − εpjκpl (l + 1)

(cid:90)

(cid:68)

Φm

l (ˆr) × ˆr · L(αp)
(cid:1)(cid:69) (cid:68)

Ψm

l , L(αp)

lm

(cid:0)κp

lm (κp, r) dS

Ym

l , φ(βp)

lm

(cid:0)κ∗

p

(cid:1)(cid:69)∗

(cid:0)κ∗

p

(cid:68)

lm

Ym

l , φ(βp)
(cid:1)(cid:69)∗

.

(3.16b)

(3.16c)

(3.16d)

(3.16e)

(3.16f)

(3.16g)

3.2.4.5 All Operators

Below are all the operators used in Chapters 1 and 2 implemented with spherical harmonic

basis functions. All pairs of test and basis functions not shown are analytically 0. The

41

vector-vector operators are

K′
p

(cid:12)
(cid:12)

Ψ Ψ

lm

K′
p

(cid:12)
(cid:12)

Φ Φ

lm

= − K′t
p

= − K′t
p

(cid:12)
(cid:12)

(cid:12)
(cid:12)

Ψ Ψ

lm

Φ Φ

lm

(cid:68)

= εpjκ2
p

Ψm

l , N(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Φm

l , M(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

= −εpjκ2
p
(cid:68)

(cid:68)

Φm

l , M(αp)

lm

Ψm

l , N(βp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp
(cid:1)(cid:69) (cid:68)

(cid:1)(cid:69)∗

(cid:0)κ∗
p
(cid:1)(cid:69)∗

Lp|Ψ Ψ

lm = − Lp|Ψ Ψ

lm = −jκp

Ψm

l , N(αp)

lm

(cid:0)κp

Ψm

l , N(βp)

lm

(cid:0)κ∗

p

(cid:12)
(cid:12)

Φ Φ

lm

= −jκp

(cid:68)

Φm

l , M(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Φm

l , M(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

=jκp

Φm

l , M(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Φm

l , M(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

Lp|Φ Φ

lm = − Lt
p
(cid:68)

Ψ Ψ

J(2)
p

(cid:12)
(cid:12)

lm

J(2)
p

(cid:12)
(cid:12)

Φ Φ

lm

=jκp

(cid:16)(cid:68)

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ψm

lm

l , N(αp)
(cid:68)

Ψm

lm

l , N(βp)
(cid:1)(cid:69) (cid:68)

Ψm

l , L(αp)

lm

(cid:0)κp

Ψm

l , L(βp)

lm

(cid:1)(cid:69)∗(cid:19)

(cid:0)κ∗

p

(cid:1)(cid:69)∗

(cid:0)κ∗

p

l (l + 1)
κ2
p
(cid:68)

Kp|Ψ Ψ

lm =εpjκ2

p

Φm

l , M(αp)
(cid:68)

lm

Kp|Φ Φ

lm = − εpjκ2

p

Ψm

l , N(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Φm

l , M(βp)

lm

(cid:1)(cid:69)∗

.

(cid:0)κ∗

p

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ψm

l , N(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

The vector-scalar operators are

Q(1)
p

P(2)
p

Q(2)
p

(cid:12)
(cid:12)

(cid:12)
(cid:12)

(cid:12)
(cid:12)

Φ Y

lm

Φ Y

lm

Φ Y

lm

(cid:68)

=εpjκ2
p

Φm

l , M(αp)
(cid:68)

lm

= − εpjκpl (l + 1)

(cid:16)(cid:68)

= − jκp

Ψm

lm

l , N(αp)
(cid:68)

+

l (l + 1)
κ2
p

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Y m

(cid:1)(cid:69)∗

(cid:0)κ∗

p

l , N(βp)
lm
(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ψm

l , L(αp)
lm
(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ym

l , φ(βp)
lm
(cid:1)(cid:69)∗

(cid:0)κ∗

p

Y m

l , N(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

Y m

l , L(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ψm

l , L(βp)

lm

The scalar-vector operators are

J(3)
p

(cid:12)
(cid:12)

Y Φ

lm

=jκp

(cid:16)(cid:68)

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ψm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

lm

l , N(βp)
(cid:1)(cid:69) (cid:68)

Y m

l , L(αp)

lm

(cid:0)κp

Ψm

l , L(βp)

lm

Y m

lm

l , N(αp)
(cid:68)

+

l (l + 1)
κ2
p
(cid:68)

M(3)
p

J(4)
p

(cid:12)
(cid:12)

(cid:12)
(cid:12)

Y Φ

lm

Y Φ

lm

= − εpjκ2
p

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Y m

l , N(αp)
lm
(cid:68)

Φm

lm

l , M(βp)
(cid:1)(cid:69) (cid:68)

(cid:1)(cid:69)∗

(cid:0)κ∗

p

Ψm

l , L(βp)

lm

(cid:1)(cid:69)∗

.

(cid:0)κ∗

p

= − εpjκpl (l + 1)

Ym

l , φ(αp)

lm

(cid:0)κp

(cid:1)(cid:69)∗(cid:19)

.

(cid:0)κ∗

p

(cid:1)(cid:69)∗(cid:19)

(cid:0)κ∗

p

(3.17a)

(3.17b)

(3.17c)

(3.17d)

(3.17e)

(3.17f)

(3.17g)

(3.17h)

(3.17i)

(3.17j)

(3.17k)

(3.17l)

(3.17m)

(3.17n)

42

The scalar-scalar operators are

D′
p

Q(3)
p

(cid:12)
(cid:12)

(cid:12)
(cid:12)

Y Y

lm

Y Y

lm

= − εpjκpl (l + 1)

(cid:68)

Y m

(cid:16)(cid:68)

= − jκp

Y m

lm

l , N(αp)
(cid:68)

+

l (l + 1)
κ2
p

l , L(αp)
lm
(cid:1)(cid:69) (cid:68)

(cid:0)κp

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Y m

l , N(βp)

lm

Ym

l , φ(βp)
lm
(cid:1)(cid:69)∗

(cid:0)κ∗

p

(cid:1)(cid:69)∗

(cid:0)κ∗

p

Y m

l , L(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Y m

l , L(βp)

lm

Sp|Y Y

lm = − jκpl (l + 1)

(cid:68)

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ym

l , φ(βp)

lm

Ym

l , φ(αp)
(cid:68)

lm

Ym

l , φ(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Y m

l , L(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

Dp|Y Y

lm = − εpjκpl (l + 1)
(cid:68)

lm = − jκpl (l + 1)

Np|Y Y

Y m

l , L(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Y m

l , L(βp)

lm

(cid:0)κ∗

p

.

(3.17s)

(3.17o)

(3.17p)

(3.17q)

(3.17r)

(cid:1)(cid:69)∗(cid:19)

(cid:0)κ∗

p

(cid:1)(cid:69)∗

(cid:0)κ∗

p

p
(cid:1)(cid:69)∗

For completeness, additional vector-vector operators defined in Chapter 1 are provided

here

ˆnp × Lp

ˆnp × Lp

ˆnp × Kp

ˆnp × Kp

(cid:12)
(cid:12)

(cid:12)
(cid:12)

(cid:12)
(cid:12)

(cid:12)
(cid:12)

Ψ Φ

lm

Φ Ψ

lm

Ψ Φ

lm

Φ Ψ

lm

(cid:68)

=εpjκp

= − εpjκp

Φm

l , M(αp)
(cid:68)

lm

Ψm

l , N(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Φm

l , M(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κ∗

p

(cid:1)(cid:69) (cid:68)

Ψm

l , N(βp)

lm

(cid:1)(cid:69)∗

(cid:0)κp
(cid:1)(cid:69) (cid:68)

(cid:68)

(cid:68)

=jκ2
p

=jκ2
p

Ψm

l , N(αp)

lm

(cid:0)κp

Φm

l , M(βp)

lm

Φm

l , M(αp)

lm

(cid:1)(cid:69) (cid:68)

(cid:0)κp

Ψm

l , N(βp)

lm

p

(cid:0)κ∗
(cid:1)(cid:69)∗

(cid:1)(cid:69)∗

.

(cid:0)κ∗

p

(cid:0)κ∗

p

(3.17t)

(3.17u)

(3.17v)

(3.17w)

3.2.5 Results

With these basis sets, we can analyze the eigenvalues and condition numbers of band-

limited versions of the operators, and compute Radar Cross Section (RCS) scattering from

spheres.We do this by isolating the block of operators associated with a given harmonic

specified by l and m. Because of the orthogonality of the basis functions, the system matrix

becomes block diagonal with block size (2Nv + Ns) × (2Nv + Ns). The eigenvalues of the

entire system comprise the 2Nv + Ns eigenvalues of each block. Likewise, the singular values

of the entire system comprise the 2Nv + Ns singular values of each block. Computing each

block’s eigenvalues and singular values is significantly more efficient, and separating the

eigenvalues by harmonic enables additional insight into the system’s behavior. Also, the

system is not dependent on harmonic index m.

43

For this analysis, the first harmonic is excluded where l = m = 0 and the highest

harmonic is for l = Nh. The finite-dimensional expansion of the operator, in effect, creates a

band-limited version. If the operator is compact, then this series converges, and the condition

number represents the system for which Nh → ∞. As is well known, the EFIE and PMCHWT

integral equations have hyper-singular operators, so the operator’s actual condition number

is unbounded.

The eigenvalues are from the generalized eigenvalue problem, Av = λBv with Akn = Zχ
kn

and Bkn = ⟨F χ

k , F χ

n⟩ (or equivalently the regular eigenvalue problem B−1Av = λv). Condition

numbers are computed from the matrix B−1A.

Rather than present results for each formulation, for dielectric formulations E-DFIE,

A − ϕ-DPIE will be shown with PMCHWT and M¨uller as a reference, and for the PEC

formulations LC-CFIE will be compared with EFIE and MFIE. For these comparisons, we

will show the following:

1. the condition number over frequency

2. the eigenvalue spectrum for several fixed orders l as frequency varies

3. the eigenvalue spectrum for all orders l < Nh for a specific frequency

3.2.5.1 Dielectric Formulations

For the first case, we examine a low contrast scenario µr2 = 1 and ϵr2 = 1.5 in Figure 3.1.

We see in Figure 3.1a that the condition number stays constant at low frequencies for M¨uller,

DFIE, and DPIE but increases for PMCHWT. At high frequencies, all formulations grow

at roughly the same rate. Next, we focus on the harmonics l ∈ [1, 3, 5] and m = 0, as

representatives of the whole system. We see in Figure 3.1b and 3.1c that at low frequency,

the eigenvalues collect away from zero and are bounded at high frequency for M¨uller, DFIE,

DPIE but collect at zero and infinity for PMCHWT. If the frequency is held constant at

10 MHz, 100 MHz, and 1 GHz, we notice that as frequency increases, we consider additional

harmonics as designed. From Figures 3.1d and 3.1e, we can justify extrapolating the findings

from examining l = 1, 3, 5 and see that the eigenvalues even for higher l are bounded and do

44

(a) Condition numbers across frequencies

(b) Eigenvalues for fixed harmonic indices for
M¨uller, E-DFIE, and A − ϕ-DPIE

(c) Eigenvalues for fixed harmonic indices for
PMCHWT.

(d) Eigenvalues for fixed frequencies for M¨uller,
E-DFIE, and A − ϕ-DPIE

(e) Eigenvalues for fixed frequencies for PM-
CHWT

Figure 3.1 Condition numbers and eigenvalues for different band-limited, dielectric integral
equations using analytic basis sets for a dielectric sphere with µr2 = 1 and ϵr2 = 1.5. Real
and imaginary parts of eigenvalues on the x- and y-axes, respectively.

45

103102101100101Electrical Size []101100101102103104105106Condition NumberMüllerPMCHWTE-DFIEA-DPIE0.50.00.5Müllerl=1l=3l=50.50.00.5E-DFIE120.50.00.5A-DPIE1212101101103Electircal Size []101101103Electircal Size []101101103Electircal Size []2.50.02.52.50.02.5PMCHWTl=12.50.02.5l=32.50.02.5l=5101101103Electircal Size []01Müller10MHz100MHz1GHz01E-DFIE0201A-DPIE0202100101l100101l100101l2.50.02.5420246PMCHWT10MHz2.50.02.5100MHz2.50.02.51GHz100101lnot collect at zero for M¨uller, DFIE, DPIE and collect closer and closer to zero and infinity

as frequency goes down for PMCHWT.

For the second case, we examine the dual scenario where µr2 = 1.5 and ϵr2 = 1 in

Figure 3.2. This is equivalent to examining the first case with the dual formulations, H-DFIE

and F − ψ-DPIE. As evidenced by Figures 3.2, there is little noticeable difference between

the first two scenarios.

For the third scenario, we consider a case where n2

1 = n2

2 with µr2 = 1/1.5 and ϵr2 = 1.5

in Figure 3.3. This scenario is unique because the DPIE formulations are second-kind only

under this condition. We see in Figure 3.3a that the condition number is much more flat

for all formulations In Figure 3.3b and 3.3c, the eigenvalues seem to follow some additional

structure but otherwise follow the same behavior as the first two cases. From Figures 3.3d

and 3.3e, we again notice minimal distinction between the M¨uller, E-DFIE, and A − ϕ-DPIE

formulations and similar takeaways for the PMCHWT formulation.

For the fourth and final case, we consider a high contrast scenario with µr2 = 1 and

ϵr2 = 20 in Figure 3.4. We see in Figure 3.4a that the condition number is worse for all

formulations. We also observe several “near-resonances” where the condition number spikes

but not to infinity. Figure 3.4b and 3.4c show that the eigenvalues are still bounded but

appear to be near and around zero. For certain frequencies, an eigenvalue will get relatively

close to the origin, which will cause what we call a “near-resonance”. From Figures 3.4d and

3.4e, we notice little significant difference between the three well-conditioned formulations

and no new takeaways for the PMCHWT formulation.

3.2.5.2 PEC Formulations

Before we examine the analytic results for the LC-CFIE formulation, we briefly describe

the well-studied Calder´on formulation [Andriulli et al., 2008] using the common framework.

We will preface the existing Calder´on formulation with ⋆ to differentiate between the two

46

(a) Condition numbers across frequencies

(b) Eigenvalues for fixed harmonic indices for
M¨uller, E-DFIE, and A − ϕ-DPIE

(c) Eigenvalues for fixed harmonic indices for
PMCHWT.

(d) Eigenvalues for fixed frequencies for M¨uller,
E-DFIE, and A − ϕ-DPIE

(e) Eigenvalues for fixed frequencies for PM-
CHWT

Figure 3.2 Condition numbers and eigenvalues for different band-limited, dielectric integral
equations using analytic basis sets for a dielectric sphere with µr2 = 1.5 and ϵr2 = 1. Real
and imaginary parts of eigenvalues on the x- and y-axes, respectively.

47

103102101100101Electrical Size []101100101102103104105106Condition NumberMüllerPMCHWTE-DFIEA-DPIE0.50.00.5Müllerl=1l=3l=50.50.00.5E-DFIE0.51.01.50.50.00.5A-DPIE0.51.01.50.51.01.5101101103Electircal Size []101101103Electircal Size []101101103Electircal Size []2.50.02.52.50.02.5PMCHWTl=12.50.02.5l=32.50.02.5l=5101101103Electircal Size []0.50.00.5Müller10MHz100MHz1GHz0.50.00.5E-DFIE020.50.00.5A-DPIE0202100101l100101l100101l2.50.02.5642024PMCHWT10MHz2.50.02.5100MHz2.50.02.51GHz100101l(a) Condition numbers across frequencies

(b) Eigenvalues for fixed harmonic indices for
M¨uller, E-DFIE, and A − ϕ-DPIE

(c) Eigenvalues for fixed harmonic indices for
PMCHWT.

(d) Eigenvalues for fixed frequencies for M¨uller,
E-DFIE, and A − ϕ-DPIE

(e) Eigenvalues for fixed frequencies for PM-
CHWT

Figure 3.3 Condition numbers and eigenvalues for different band-limited, dielectric integral
equations using analytic basis sets for a dielectric sphere with µr2 = 1/1.5 and ϵr2 = 1.5. Real
and imaginary parts of eigenvalues on the x- and y-axes, respectively.

48

103102101100101Electrical Size []101100101102103104105106Condition NumberMüllerPMCHWTE-DFIEA-DPIE0.250.000.25Müllerl=1l=3l=50.250.000.25E-DFIE1.01.50.250.000.25A-DPIE1.01.51.01.5101101103Electircal Size []101101103Electircal Size []101101103Electircal Size []2.50.02.52.50.02.5PMCHWTl=12.50.02.5l=32.50.02.5l=5101101103Electircal Size []0.250.000.25Müller10MHz100MHz1GHz0.250.000.25E-DFIE1.01.50.250.000.25A-DPIE1.01.51.01.5100101l100101l100101l504202468PMCHWT10MHz50100MHz501GHz100101l(a) Condition numbers across frequencies

(b) Eigenvalues for fixed harmonic indices for
M¨uller, E-DFIE, and A − ϕ-DPIE

(c) Eigenvalues for fixed harmonic indices for
PMCHWT.

(d) Eigenvalues for fixed frequencies for M¨uller,
E-DFIE, and A − ϕ-DPIE

(e) Eigenvalues for fixed frequencies for PM-
CHWT

Figure 3.4 Condition numbers and eigenvalues for different band-limited, dielectric integral
equations using analytic basis sets for a dielectric sphere with µr2 = 1 and ϵr2 = 20. Real
and imaginary parts of eigenvalues on the x- and y-axes, respectively.

49

103102101100101102Electrical Size []101100101102103104105106Condition NumberMüllerPMCHWTE-DFIEA-DPIE10010Müllerl=1l=3l=510010E-DFIE02010010A-DPIE020020100102104Electircal Size []100102104Electircal Size []100102104Electircal Size []2.50.02.52.50.02.5PMCHWTl=12.50.02.5l=32.50.02.5l=5100102104Electircal Size []10010Müller10MHz100MHz1GHz10010E-DFIE02010010A-DPIE020020100101102l100101102l100101102l1001050510PMCHWT10MHz100100MHz1001GHz100101102lformulations.

P ⋆LC-CFIE

1

=

W 1 =

(cid:18)

(cid:19)

1 1

(cid:18)

α
η0

(cid:19)
1

(cid:16)
ˆn1 × ˜L1 ◦ (cid:0)κ2

1η1ˆn1 × L1

(cid:1), I − K1

(cid:17)

Z ⋆LC-CFIE

1

= diag


b⋆LC-CFIE

1

=




ˆn1 × ˜L1 ◦ (cid:0)jκ1Ei
ˆn1 × Hi

1

1 × ˆn1

(cid:1)






x⋆LC-CFIE =Js.

(3.18a)

(3.18b)

(3.18c)

(3.18d)

(3.18e)

We examine the results in Figure 3.5 where we show the standard EFIE, MFIE, and

CFIE behavior for reference and the related LC-CFIE and ⋆ LC-CFIE formulations. What

is immediately noticeable from Figure 3.5a is that neither the EFIE, MFIE, nor CFIE

formulation are well-conditioned at either low or high frequencies. In contrast, in Figure 3.5b,

we see excellent conditioning for ⋆ LC-CFIE and excellent low-frequency and comparable

high-frequency conditioning from the LC-CFIE. From Figure 3.5c, we observe that the EFIE

and CFIE are unbounded and how the MFIE and CFIE do not collect at the origin at low

frequencies. That is not to say that the MFIE does not have resonances. The CFIE has no

resonances as all the eigenvalues are away from the origin. In contrast, both the LC-CFIE and

⋆ LC-CFIE are bounded, do not collect at the origin, and are far from the origin. The main

difference between the two Calder´on formulations is that the new LC-CFIE formulation has a

“bubble” that grows with harmonic index l. This is directly responsible for the difference in

high-frequency behavior from the starred counterpart. Finally, in Figures 3.5e and 3.5f, we see

eigenvalues always clustering away from zero for the LC-CFIEand ⋆ LC-CFIE formulations

but not so with the EFIE, MFIE, or CFIE formulation, and for the LC-CFIE formulation, the

“bubble” manifests itself at higher frequencies and is the cause of the increased conditioning

compared to ⋆ LC-CFIE.

50

(a) Condition numbers across frequencies for
EFIE, MFIE, and CFIE

(b) Condition numbers across frequencies for LC-
CFIE and ⋆ LC-CFIE

(c) Eigenvalues for fixed harmonic indices for
EFIE, MFIE, and CFIE

(d) Eigenvalues for fixed harmonic indices for LC-
CFIE and ⋆ LC-CFIE.

(e) Eigenvalues for fixed frequencies for EFIE,
MFIE, and CFIE

(f) Eigenvalues for fixed frequencies for LC-CFIE
and ⋆ LC-CFIE.

Figure 3.5 Condition numbers and eigenvalues for different band-limited, integral equations
using analytic basis sets for a PEC sphere. Real and imaginary parts of eigenvalues on the x-
and y-axes, respectively.

51

103102101100101Electrical Size []101100101102103104105106Condition NumberEFIEMFIECFIE103102101100101Electrical Size []101100101102103104105106Condition NumberLC-CFIE LC-CFIE101EFIEl=1l=3l=5101MFIE101101CFIE101101101101103Electircal Size []101101103Electircal Size []101101103Electircal Size []0.00.5LC-CFIEl=1l=3l=50.51.00.00.5 LC-CFIE0.51.00.51.0101101103Electircal Size []101101103Electircal Size []0.02.5EFIE10MHz100MHz1GHz0.02.5MFIE200.02.5CFIE2020100101l100101l100101l02LC-CFIE10MHz100MHz1GHz0202 LC-CFIE0202100101l100101l3.3 Piecewise Analysis

For analysis on non-spherical objects, we can utilize a triangular mesh with Nf flat-

triangles, Ne edges, and Nn nodes to approximate the true surface as is typically done.

Piecewise basis and testing functions will not have the helpful analytic integration properties

of the spherical harmonics. However, they will enable the analysis of arbitrarily shaped

objects, which is necessary for practical utility.

3.3.1 RWG functions for Bv

RWG functions are a heavily utilized set of basis functions for representing vector functions

in CEM. They are excellent for demonstrating how to implement the formulations presented

in Chapter 1 and 2. These functions are defined for each edge with a domain of the two

connected triangles as

where

f ne (r) =




lne
2A±
ne



0

ϱ±
ne (r)

r ∈ T ±
ne

otherwise

ϱ±
ne (r) =




± (cid:0)r − p±
ne

(cid:1)

r ∈ T ±
ne



0

otherwise

(3.19a)

.

(3.19b)

In the above equations, lne is the length of edge ne, A±

ne is the area associated with triangle

ne is the node of the triangle T ±

T ±
ne, and p±
these basis functions, we can define Bv as

ne opposite the edge ne [Rao et al., 1982]. With

Bv =

(cid:18)

(cid:19)

f 1

. . .

f ne

. . .

f Ne

(3.20)

with Nv = Ne.

3.3.2 Pulse functions for Bs

Pulse functions are the simplest basis functions to represent scalar functions over a

triangular mesh. These functions are constant over each face of the mesh

pnf (r) =

1 r ∈ Tnf

0 otherwise





52

(3.21)

where Tnf is the nf triangle.

With these basis functions, we can define Bs as

Bs =

(cid:18)

(cid:19)

p1

. . . pnf

. . . pNf

(3.22)

with Ns = Nf .

3.3.3 Hat functions for Bs

While hat functions have been implemented for these new formulations, they are unneces-

sary for the primary purpose of this chapter. Rather than showing two sets of results, one

with pulse functions and one with hat functions, pulse functions are used for the remainder

of this thesis. This section serves as a reminder of how additional basis sets can be used to

construct Bv and Bs without any change to the underlying formulation.

The second most straightforward scalar functions on a triangular mesh are hat functions

(pyramidal functions). These functions are associated with each node and linearly decrease

from 1 at the node to 0 at the edges, forming a ring around the node. We can define these

functions fairly simply by reusing RWG functions

hnn =




1 − ˆui

nn · f i
nn

r ∈ T i
nn



0

otherwise

(3.23)

where each node has Nnn triangles connected to it with index i ∈ [1, Nnn]. Each triangle

nn has an edge that does not include node nn. Those edges have unit normal vectors ˆui
T i
nn
perpendicular to them pointing out of the triangle and an RWG function f i

nn oriented such

that T i

nn is its positive triangle.

With these basis functions, we can define Bs as

Bs =

(cid:18)

(cid:19)

h1

. . . hnn

. . . hNn

(3.24)

with Ns = Nn.

These functions typically result in fewer unknowns for Bs than the pulse functions because

Nn ≤ Nf for closed objects. However, more quadrature points will be needed for each integral

53

because the domain of support covers more than one triangle. If we count the number of

integrals, each triangle is integrated over once for pulse functions but three times for hat

functions.

3.3.4 Zero-Mean Constraint

In order to effect the zero-mean constraint for a closed tessellated object through the

basis function coefficients, we use a Lagrange multiplier as in [Dault and Shanker, 2015]. This

adds a row and column to the matrix for each Lagrange multiplier while not significantly

modifying the convergence properties of the formulation as confirmed experimentally.

As an alternative, a rank-1 update will avoid the matrix with the added cost of an

additional matrix solution as done in [Hawkins et al., 2022].

3.3.5 Numeric Integration

The key challenge of working with piecewise tessellations is evaluating these operators

when the k testing domain is close to or shared with the n basis domain. While several

operators are familiar, some are not. We have taken a straightforward approach to evaluating

these integrals—singularity subtraction. There are other methods of evaluating these integrals

[Tihon and Craeye, 2018]. However, they are not necessary to demonstrate the crux of this

chapter—the demonstration of how to implement the new formulations for discrete piecewise

tessellations and their properties. To that end, we note that even though there are several

unfamiliar operators, the singular integrals necessary are the same as those encountered in

the EFIE and MFIE; all other singular integrals can be evaluated via a combination of four

base singular integrals. These singular integrals, along with an example of the complete

treatment of one of the operators, are presented in the next section. Standard quadrature

techniques are applicable for numeric integration when the basis and testing domains are

sufficiently separated and not discussed here.

54

r′

d

p0

R+

R0 R− R

ρ0 ρ−

p0

p−

ρ+

ρ

ρ±

p+

r

ϱ±

p±

Figure 3.6 Vector and point definitions for singular integration over a testing triangle.

3.3.5.1 Singular Integrals

We subtract the first two terms from the Taylor series expansion of the exponential

function. As will be evident, the singular integrals needed are identical to the ones needed

for the EFIE and MFIE and can be found in [Wilton et al., 1984, Hodges and Rahmat-

Samii, 1997]. It follows that, if necessary, one can use better rules [Fink et al., 2008, Botha,

2013, Tihon and Craeye, 2018] if so desired. We use Figure 3.6 to set the stage for definitions,

and it is here purely for completeness. The definitions here follow those in [Wilton et al.,

1984].

It is well known that the following integrals can be evaluated analytically as r′ approaches

a testing triangle.

I1/R =

Iρ/R =

I1/R3 =

Iρ/R3 =

dS

dS

(cid:90) 1
R
(cid:90) ρ
R
(cid:90) 1

(cid:90) ρ

R3 dS

R3 dS

55

(3.25a)

(3.25b)

(3.25c)

(3.25d)

where the integrals are defined according to Figure 3.6. Similar integrals can be defined and

evaluated as r approaches a source triangle with the only difference being that R points away

from the triangle, so the integrals involving R have an extra factor of −1. Singular integrals

over the source domain will be denoted with a prime.

We used a seven-point Gauss-Legendre rule for the integrals over a triangle to effect

the integration. For the hypersingular integral Np, we treat it as in [Terai, 1980] and use a

fourteen-point Gauss-Legendre rule for the line integral.

As an aside, two issues naturally crop up when surfaces are modeled using higher-order

geometric representations: (a) evaluation of singular integrals and (b) cost of a higher-order

quadrature rule. The integrals in (3.25) have been dealt with in the community using either

a mapping of the contour to a flat patch, which allows the use of singularity subtraction or

singularity cancellation techniques. Amelioration of the costs associated with the evaluation of

these integrals for higher order geometries has been dealt with in [Alsnayyan et al., 2020, Dault

and Shanker, 2015] by integrating with a wideband Fast Multipole Method (FMM), which

uses a combination of an adaptive quadrature rule and singularity cancellation around a

small neighborhood of the singularity.

The techniques used below rely on the fact that we use a flat tessellation and are specialized

in using RWG and pulse functions for basis and testing functions.

3.3.5.2 Candidate Inner Product Evaluation

To illustrate a general approach, consider the evaluation of P(2)
p

The integration is over the testing domain T ±

(cid:12)
(cid:12)
(cid:12)kn
k and the basis domain Tn. To simplify the

as defined in (3.14).

presentation, only consider T +

k (the other follows trivially). Given that ˆnp is constant over

T +
k , manipulation of the integral results in

P(2)
p

(cid:12)
(cid:12)kn

=

=

(cid:90)

T +
k

(cid:90)

Tn

f k (r) · ˆnp × ∇

(cid:90)

Tn

Gp (r, r′) pn (r′) dS′dS

pn (r′) ˆnp ·

(cid:90)

T +
k

∇ Gp (r, r′) × f k (r) dSdS′.

(3.26a)

(3.26b)

56

As is usually done for singularity subtraction, we use a Taylor series expansion for Gp or

∇ Gp depending on the integral. We then decompose the integral into two parts: one with a

removable singularity and another non-singular integral that can be evaluated numerically. As

mentioned earlier, we subtract the first two terms of the expansions to construct the singular

part (although the second terms are not singular), and the rest of the series is captured in

the non-singular part

Gp (r, r′) =

∇ Gp (r, r′) =

=

1
4π
1
4π
(cid:124)

1
4π
1
4π
(cid:124)

(cid:18) 1
R
(cid:18) 1
R

(cid:19)

− jκp

+ O (R)

(cid:19)

(cid:18)

− jκp

+

G −

(cid:123)(cid:122)
Singular

(cid:125)

(cid:124)

1
4π

(cid:18) 1
R
(cid:123)(cid:122)
Non-singular

− jκp

(cid:19)(cid:19)

(cid:125)

(cid:19)

(cid:18)

−

(cid:18)

−

R
R3 −
R
R3 −
(cid:123)(cid:122)
Singular

κ2
pR
2R
κ2
pR
2R

+ RO (1)

(cid:19)

(cid:18)

+

∇ G −

(cid:125)

(cid:124)

(cid:18)

1
4π

−

R
R3 −
(cid:123)(cid:122)
Non-singular

(3.27a)

(3.27b)

(3.27c)

(3.27d)

(cid:19)(cid:19)

κ2
pR
2R

(cid:125)

We only discuss the singular portion of the integral denoted below using superscript S as in

|S. For (3.26), this can be written as

P(2)
p

(cid:12)
(cid:12)

S

kn

=

−1
4π

(cid:90)

Tn

ˆnp ·

(cid:32)(cid:90)

T +
k

R
R3 × f k (r) dS +

κ2
p
2

(cid:90)

T +
k

It follows that to evaluate this integral we need to evaluate (cid:82) R

(cid:33)

dS′.

× f k (r) dS

R
R
R3 ×f k (r) dS and (cid:82) R

(3.28)

R ×f k (r) dS.

In this example, the integral over the source domain can be evaluated numerically.

From (3.25), one can derive the following analytic integrals:

(cid:90) R
R

(cid:90)

T ±
k

(cid:90)

T ±
k
(cid:90) R

R
R
f k (r)
R

IR/R =

IR/R×f =

If/R =

IR/R3 =

IR/R3×f =

dS = Iρ/R + dI1/R

× f k (r) dS = ±

lk
2A±
k

(cid:0)d + ρ±(cid:1) × IR/R

dS = ±

lk
2A±
k

(cid:0)Iρ/R − ρ±I1/R

(cid:1)

R3 dS = Iρ/R3 + dI1/R3
R
lk
R3 × f k (r) dS = ±
2A±
k

(cid:90)

T ±
k

(cid:0)d + ρ±(cid:1) × IR/R3.

(3.29a)

(3.29b)

(3.29c)

(3.29d)

(3.29e)

57

The integral If k/R is not necessary for this example but is for the other operators in the

following section.

We finally arrive at

P(2)
p

(cid:12)
(cid:12)

S

kn

=

−1
4π

(cid:90)

Tn

(cid:18)

ˆnp ·

IR/R3×f +

(cid:19)

dS′.

IR/R×f

κ2
p
2

(3.30)

3.3.5.3 All Operators

The process for addressing the rest of the singular integrals is similar. Again, the integrals

are defined according to Figure 3.6 where the integration is over a testing triangle, and it is

important to remember if integrating over the source domain to add an extra factor of −1

for singular integrals involving R. Singular integrals over the source domain will be denoted

with a prime.

Note, the singular integrals are either due to Gp or ∇ Gp, are over either testing or source

domains, and can all be formed as combinations of (3.25). Just as the normal over each

triangle is constant due to the flat-tessellation, ˆn · R and ˆn′ · R are also constant over the

testing and source domains respectively. In addition, ∇s · f and ∇′

s · f ′ are constant.

Below, the expressions we use to treat the singular integrals are cataloged. This catalog

is one of several equivalent sets of expressions that can be used for singularity subtraction.

The singular integrals for the vector-vector operators are

K′
p

(cid:12)
(cid:12)

S

kn

= − K′t
p

(cid:12)
(cid:12)

S

kn

=

Lp|S

kn = − Lt

p

(cid:12)
(cid:12)

S

kn

=

(cid:0)ˆn′

n

p × f ′
(cid:18)

T +
n

(cid:90)

−1
4π
1
4π

(cid:90)

T +
k

f k ·

I ′
f/R − jκp

f ′
ndS′

(cid:90)

T +
n

κ2
p
2
(cid:19)

(cid:18)

IR/R3×f +

(cid:1) ·

(cid:19)

dS′

IR/R×f

J(2)
p

(cid:12)
(cid:12)

S

kn

=

Kp|S

kn =

1
4π

1
4π

(cid:90)

T +
k

(cid:90)

T +
k

−

1
κ2
p
(cid:0)ˆnp × f k

(cid:1) ·

∇s · f k∇′
(cid:18)

s · f ′
(cid:18)

n

ˆn′

p ×

I ′
f/R − jκp

(cid:0)I ′

1/R − jκpA+
n

(cid:1) dS

(cid:90)

T +
n
(cid:19)

(cid:19)(cid:19)

dS

f ′
ndS′

dS.

I ′
R/R×f

(cid:0)ˆnp × f k

(cid:1) ·

(cid:18)

I ′
R/R3×f +

κ2
p
2

58

(3.31a)

(3.31b)

(3.31c)

(3.31d)

The singular integrals for the vector-scalar operators are

The singular integrals for the scalar-vector operators are

Q(1)
p

P(2)
p

Q(2)
p

(cid:12)
(cid:12)

(cid:12)
(cid:12)

(cid:12)
(cid:12)

S

kn

S

kn

S

kn

=

=

=

−1
4π

−1
4π

−1
4π

(cid:90)

Tn

(cid:90)

Tn

(cid:90)

Tn

ˆn′

p ·

(cid:18)

(cid:18)

IR/R3×f +

ˆnp ·

IR/R3×f +

(cid:19)

(cid:19)

dS′

dS′

IR/R×f

IR/R×f

κ2
p
2
κ2
p
2

(cid:32)

ˆn′

p · ˆnp ×

If/R − jκp

J(3)
p

(cid:12)
(cid:12)

S

kn

=

M(3)
p

(cid:12)
(cid:12)

S

kn

=

J(4)
p

(cid:12)
(cid:12)

S

kn

=

1
4π

−1
4π

1
4π

(cid:90)

Tk
(cid:90)

Tk

(cid:90)

Tk

ˆnp · ˆn′

p ×

(cid:18)

I ′
f/R − jκp

(cid:18)

I ′
R/R3×f +

(cid:18)

I ′
R/R3×f +

ˆnp ·

ˆn′

p ·

κ2
p
2
κ2
p
2

I ′
R/R×f

(cid:19)

dS.

I ′
R/R×f

(cid:90)

T +
k

(cid:90)

T +
n

(cid:33)

f kdS

dS′.

(cid:19)

dS

f ′
ndS′

(cid:19)

dS

The singular integrals for the scalar-scalar operators are

D′
p

(cid:12)
(cid:12)

S

kn

=

Q(3)
p

(cid:12)
(cid:12)

S

kn

=

Sp|S

kn =

Dp|S

kn =

−1
4π
1
4π
1
4π

1
4π

(cid:90)

Tn

(cid:90)

Tk

(cid:90)

Tk

(cid:90)

Tk

(cid:18)

I1/R3 +

(cid:0)ˆnp · R(cid:1)

κ2
p
2
ˆnp · ˆn′ (cid:0)I1/R − jκpA+

n

(cid:1) dS

(cid:19)

dS′

I1/R

I1/R − jκpA+

n dS

(cid:0)ˆn′

p · R(cid:1)

(cid:18)

I1/R3 +

(cid:19)

I1/R

dS.

κ2
p
2

(3.31e)

(3.31f)

(3.31g)

(3.31h)

(3.31i)

(3.31j)

(3.31k)

(3.31l)

(3.31m)

(3.31n)

As stated earlier in this chapter, the Np operator is hyper-singular and is dealt with as in

[Terai, 1980] rather than with singularity subtraction.

For completeness, additional vector-vector operators defined in Chapter 1 are provided

here

ˆnp × Lp

(cid:12)
(cid:12)

S

kn

=

1
4π

(cid:90)

T +
k

(cid:18)(cid:18)

(cid:0)f k × ˆnp

(cid:1) ·

I ′
f/R − jκp

(cid:90)

(cid:19)

f ′
ndS′

(3.31o)

(cid:19)(cid:19)

dS

I ′
R/R

T +
n
κ2
p
2
(cid:19)

(cid:18)

I ′
R/R3 +

κ2
p
2

I ′
R/R×f

dS.

(3.31p)

−

1
κ2
p

∇′

s · f ′

n

ˆnp × Kp

(cid:12)
(cid:12)

S

kn

=

1
4π

(cid:90)

T +
k

f k ·

(cid:18)

I ′
R/R3×f +

59

(a) 100 MHz

(b) 1 Hz

(c) 1 µHz

Figure 3.7 RCS of the dielectric sphere with µr2 = 1 and ϵr2 = 1.5.

The integrals (cid:82)

T +
k

f kdS and (cid:82)

T +
n

f ′
ndS′ are not singular and can be evaluated analytically

or numerically.

Piecewise geometric descriptions come with the challenge of accurately describing the

scattering physics from objects with singular features, e.g., cones, tips, edges, These are

well-known challenges (and the formulations from Chapter 2 are not immune to these as they

arise from the geometry as opposed to the formulation), and one could use singular basis

functions [Graglia et al., 2013] or adaptive basis functions [Nair and Shanker, 2011, Dault

and Shanker, 2015] to overcome some of the bottlenecks partially.

Utilizing hat instead of pulse functions changes the expressions above but not in a way

that requires additional singular integrals and is relatively straightforward due to how hat

functions defined in (3.23) are related to RWG functions.

3.3.6 Results

We begin the piecewise analysis by considering the scattering from the dielectric sphere

with radius 1 m. The sphere is meshed into 1280 patches and 1920 edges with an average edge

length at 100 MHz of about λ1/19.9 with λ1 being the wavelength in the exterior region. In

Figures 3.7, we examine RCS data at three different frequencies, specifically at 100 MHz, 1 Hz,

and 1 µHz due to an incident wave with an arrival vector of ˆκ = ˆz, polarized along E0 = ˆx,

and measured at points along θ ∈ [−π, π] and ϕ = 0. In Figure 3.7a, we see that there is

good agreement among all formulations, including the analytic formulation at 100 MHz. In

Figure 3.7b, we see that the E-DFIE, M¨uller, and PMCHWT formulations have already

60

/20/23020100RCS [dBm²]MüllerE-DFIEA-DPIEMie/20/2600550500450400350300RCS [dBm²]MüllerE-DFIEA-DPIEMie/20/2800700600500400300RCS [dBm²]MüllerE-DFIEA-DPIEMie(a) Matvec count using TFQMR method for dif-
ferent integral equations across frequencies.

(b) Convergence plot of the relative residual
∥residual∥ / ∥RHS∥ using the TFQMR method
for different integral equations. The thickness of
the line indicates the frequency, while the style
(dotted, dashed, solid) indicates the formulation.

Figure 3.8 Results from the dielectric sphere discretized with RWG and pulse basis sets with
µr2 = 1 and ϵr2 = 1.5.

diverged from the analytic solution at 1 Hz whereas the A − ϕ-DPIE and Mie agree very

well with each other. Finally, in Figure 3.7c, we see that the DPIE still matches the analytic

solution while the E-DFIE and M¨uller solutions are off the chart at 1 µHz.

For this next experiment, we used the same piecewise discretized sphere and collected data

from 1 µHz to 100 MHz. For this demonstration, we use the Transpose-Free Quasi-Minimal

Residual (TFQMR), but we have verified that other iterative solvers show similar behavior. A

relative tolerance of 10−12 is used for two reasons. First, low-frequency convergence, especially

to −500 dB m2, necessitated the need for a smaller than typical tolerance and a relative as

opposed to absolute tolerance

∥residual∥ ≤ tol × ∥RHS∥ .

(3.32)

Second, as we examine convergence across frequencies, we wanted to keep the relative tolerance

constant so that the iteration count comparison would be fair. The TFQMR solver is stopped

after 500 iterations unless it stops due to a TFQMR related breakdown or the tolerance has

been met. For TFQMR, there are two matrix-vector multiplications (matvecs) per iteration.

61

104101102105108Frequency [Hz]101102103Matvec CountMüllerPMCHWTE-DFIEA-DPIE100101102Iteration Count1011109107105103101101Relative Residual NormMüllerPMCHWTE-DFIEA-DPIE100 MHz10 kHz10 Hz(a) Tessellated geometry with in-
cident and observation angles.

(b) RCS at 90 MHz.

(c) Matvec count using TFQMR
method for different integral
equations across frequencies

Figure 3.9 The dielectric National Aeronautics and Space Administration (NASA) Almond
using piecewise basis sets with µr2 = 1 and ϵr2 = 1.5.

From Figure 3.8a, it is apparent that the matvec and iteration count remains low as the

frequency goes toward 0 Hz. It is also evident that PMCHWT requires significantly more

iterations than the threshold except at 1 µHz, where the iterative solver does converge within

the tolerance but not to the correct solution. There is a moderate increase for the E-DFIE

and A − ϕ-DPIE as one tends to higher frequencies. This trend mirrors that seen in the

M¨uller system [Li et al., 2019]. We note that methods to mollify the behavior of M¨uller

systems have been addressed by using intermediate Buffa-Christianssen basis sets [Yan et al.,

2010], but investigating these phenomena lies outside the main goals of this thesis and can

be investigated in the future.

We turn our attention to the convergence of the iterative solver. In Figure 3.8b, we see

that the E-DFIE, A − ϕ-DPIE, and M¨uller formulations converge very quickly in a few

iterations to a solution within the relative tolerance specified. The E-DFIE and M¨uller

formulation are well-conditioned but do not recover the RCS at a low enough frequency. In

other words, even though the solver finds a solution within a few iterations, it is not helpful

in the low-frequency regime.

The results of Figure 3.8 demonstrate fast convergence for the E-DFIE and A − ϕ-DPIE

62

/20/280706050403020RCS [dBm²]MüllerE-DFIEA-DPIE104105106107108Frequency [Hz]101102103Matvec CountMüllerPMCHWTE-DFIEA-DPIE3.3.6.1 Dielectric Almond

Next, we analyze scattering from a dielectric NASA almond, as seen in Figure 3.9a. The

NASA almond meshes into 896 patches and 1344 edges. The average edge length at 90 MHz

is about λ1/28.5.

In Figure 3.9b, the RCS is shown at 90 MHz with an incident wave arriving from ˆκ = ˆz,

polarization E0 = ˆy, and measured at points along θ ∈ [−π, π] and ϕ = π

2 . We see that there

is excellent agreement among the E-DFIE, A − ϕ-DPIE, and M¨uller formulations at 90 MHz.

Next, the number of TFQMR iterations needed to converge to an error of 10−12 is shown

in Figure 3.9c. The data shown in this figure is for a range of frequencies from 9 kHz to

90 MHz. Again, at the lower frequencies, only A − ϕ-DPIE is the only formulation that is

accurate. The number of iterations is relatively constant for the E-DFIE, A − ϕ-DPIE, and

M¨uller.

Overall, as is evident from Figure 3.9, the favorable properties of the E-DFIE and

A − ϕ-DPIE hold for non-canonical geometries.

3.4 Conclusion

Discretization with spherical harmonics enables analysis that provides insight into the

formulations, even more than what is presented here. With piecewise basis functions, we can

use these formulations to analyze arbitrary shapes and model actual use cases.

Implementing both discretizations provides a very practical advantage because they can

be used to cross-validate each other. A simple change of basis from piecewise functions to

spherical harmonics of a particular order can help identify implementation issues and give

valuable insight into discretization limitations.

At this point in the thesis, a framework for formulations was presented in Chapter 2,

and all of the formulations are implementable using the framework and integrals provided in

this chapter. In addition, indirect formulations and formulations that are developed under

different frameworks but involve the same operators are implementable with what is covered

above.

63

CHAPTER 4

ACCELERATION

4.1 Introduction

As noted earlier, the formulations introduced in Chapter 2 involve more operators than

encountered in the classical formulations and result in large systems. As shown in Chapter 3,

it is possible to implement these formulations using piecewise basis sets and construct a dense

system of equations to solve for the unknown coefficients of (3.3b). However, an increased

system size typically translates to a significant increase in solve time even if the system is

well-conditioned due to the necessity of evaluating each entry in the system. Fortunately,

this burden can be alleviated by using the MLFMA for computing far interactions, which

constitute most of the system. The advantages are realized by reconstructing the potential or

fields before they (or their derivatives or traces) are measured. While it is well understood, we

briefly overview the MLFMA to introduce the vocabulary necessary to discuss modifications

to evaluate the new operators.

4.2 MLFMA Background

The derivation of FMM and its multilevel variant MLFMA are well known [Coifman et al.,

1993, Song and Chew, 1995, Song et al., 1997, Sarvas, 2003, Vikram et al., 2009, Hughey

et al., 2018]. Without any loss of generality, we will present the MLFMA through the 2-level

description and drop the region index p from the notation. FMM starts with the integral

representation of the Green’s function

G (X + d) ≈ −

jκ
(4π)2

(cid:90)

S2

exp (−jκ · d) T (κ, X) d2ˆκ,

(4.1)

where the translation operator T is given by

T (κ, X) =

∞
(cid:88)

n=0

(−j)n (2n + 1) h(2)

n (κX) Pn

(cid:16)

ˆκ · ˆX

(cid:17)

,

(4.2)

with κ = κˆκ. Here, S2 denotes the unit ˆκ-sphere, parameterized by (θ, ϕ) ∈ [0, π] × [0, 2π].

We note that ˆκ = ˆκ (θ, ϕ), and use these notations interchangeably. As in any MLFMA

64

scheme, the computational domain is embedded into a bounding box that is used to construct

a uniform oct-tree, the smallest of whose boxes we refer to as leaf boxes. As usual, one designs

operators to map from sources in the leaf boxes to observers in other leaf boxes, provided the

boxes satisfy an interaction criterion.

4.3 Implementation Details

As is evident from how all of the formulations in this thesis are constructed, the repre-

sentation equations relate the scattered potential or field to traces and traces of derivatives

acting on the total potential or field collectively referred to as x in (2.1b). Likewise, one

measures traces or traces of derivatives acting on the scattered potential or field. Within an

MLFMA scheme, all these operators can be evaluated spectrally before and after translation.

This implies that one can send scattered potentials or fields up and down the tree and then

evaluate the necessary derivatives (and their traces) just before testing.

Consider a source box s whose domain is denoted by Ωs and centroid by rc

s. The potential

or field radiated by these sources is observed in an observation box with domain Ωo and

centroid rc

o. Assume that in the source box, there exist N s

f functions and N s

e RWG basis

functions. Each representation integral is evaluated using a quadrature rule. Each pulse

basis function is evaluated at N f

q quadrature points located at r′

nf i with quadrature weights

wnf i. Each RWG basis function is evaluated at N e

q quadrature points located at r′

nei with

quadrature weights wnei and unit normal ˆn′

nei (which changes based on which triangle T ±

ne of

the RWG function the quadrature point resides on).

Several formulations are presented in this thesis, each with slightly different representation

integrals. The following presents the most general representation integrals (1.20) and (2.17).

We will utilize placeholders x and x to denote any vector and scalar that satisfies the vector

and scalar Helmholtz equation, respectively. In that way, the following prescription will cover

all of the DFIE and DPIE formulations. For clarity, the coefficients in y will be split into

coefficients αnf , βnf , ane, bne, γnf , and σnf as needed for the formulations and for simplicity

P and W will be omitted without any loss in generality. Then Charge-to-Multipole (C2M)

65

expansions are as follows

fs (r, κ) = exp (−jκ · (rc

s − r))

Vx

s (κ) =

N s
f
(cid:88)

N f
q
(cid:88)

wnf ifs

(cid:16)

r′
nf i, κ

(cid:17)

i=1

nf =1
(cid:16)

jκ · ˆn′
nf

αnf − βnf

(cid:17)

(cid:17)

(cid:16)

r′
nf i

pnf

Vx

s (κ) =

N s
f
(cid:88)

N f
q
(cid:88)

wnf ifns

(cid:16)

r′
nf i, κ

(cid:17)

(4.3a)

(4.3b)

(4.3c)

i=1

nf =1
(cid:16)

jκγnf − ˆnnf
N s
e(cid:88)

N e
q
(cid:88)

+

(cid:17)

(cid:17)

(cid:16)

r′
nf i

pnf

σnf

wneifs

(cid:0)r′

nei, κ(cid:1)

ne=1

i=1

(cid:0)−ˆnneiane − jκbne

(cid:1) × f nei

(cid:0)r′

nei

(cid:1) .

Multipole-to-Local (M2L) transformations from Ωs multipoles to local expansion Ωo are

unchanged

Ux

o (κ) =

Ux

o (κ) =

(cid:88)

s
(cid:88)

s

T (κ, rc

o − rc

s) Vx

s (κ)

T (κ, rc

o − rc

s) Vx

s (κ) .

(4.4a)

(4.4b)

As an aside for the multilevel case, Multipole-to-Multipole (M2M) and Local-to-Local (L2L)

expansions and M2L translations on all levels are unchanged.

For Local-to-Observer (L2O) measurements, dipole and monopole observations are located

66

at r ∈ box Ωo can be evaluated using requisite inner products with

go (r, κ) =

xF F (r) =

∇xF F (r) =

xF F (r) =

∇ × xF F (r) =

∇ · xF F (r) =

S2

(cid:90)

S2

(cid:90)

S2

(cid:90)

S2

(cid:90)

S2

−jκ
(4π)2 exp (−jκ · (r − rc
o))
(cid:90)
o (κ) d2ˆκ

go (r, κ) Ux

−jκgo (r, κ) Ux

o (κ) d2ˆκ

go (r, κ) Ux

o (κ) d2ˆκ

−jκ × go (r, κ) Ux

o (κ) d2ˆκ

−jκ · go (r, κ) Ux

o (κ) d2ˆκ.

(4.5a)

(4.5b)

(4.5c)

(4.5d)

(4.5e)

(4.5f)

Here, superscript F F denotes the far-field scattered quantity. Minor modifications are

necessary to use this for the reduced formulations, for example, setting σ = 0 for the Reduced-

Decoupled Field Integral Equation (R-DFIE) and utilizing σ = κ2

jω α for R-DPIE. For the

classical formulations, utilize their respective representation integrals to develop the proper

C2M operator.

We note the following: (a) only changes to an existing MLFMA implementation to

accommodate the new formulations are at the leaf levels in the C2M and L2O operators. (b)

For the DPIE, we chose to have a four-tree configuration to represent the three Cartesian

components of As and ϕs or Fs and ψs and a three-tree configuration for the DFIE to

represent the Cartesian components of Es or Hs. We chose to send the Cartesian (rather than

θ − ϕ) components primarily to enable simpler future integration with wideband MLFMA

[Vikram et al., 2009]. Furthermore, (c) because we take all derivatives spectrally, we must

carefully account for the overall spectral content to design integration rules.

4.4 Results

The main results of this chapter are twofold: the efficacy of DFIE and DPIE for analysis

of electrically large objects and the effectiveness of the MLFMA. To that end, we will show

several RCS plots for a variety of objects and timings of the MLFMA accelerated (i.e. far-field)

portion of the matrix-vector multiplications.

67

(a) 300 MHz

(b) 500 MHz

(c) 1 GHz

Figure 4.1 RCS plots from MLFMA accelerated solutions for a 1 m radius, dielectric sphere
with µr2 = 1 and ϵr2 = 1.5.

Like before, rather than showing results for each form of the DFIE and DPIE, we choose

to focus on the E-DFIE and A − ϕ-DPIE as representatives. The main reference point will

be the three-tree MLFMA accelerated M¨uller formulation.

4.4.1 Accuracy

We first examine the accuracy of the MLFMA by comparing the RCS to the Mie series

RCS. After that, we show agreement between the different formulations of arbitrarily shaped

objects.

4.4.1.1 Spheres

We begin by considering the scattering from different refinements of a dielectric sphere

with a radius of 1 m. We present three RCS measurements due to an incident plane wave

arriving from ˆκ = ˆz, polarized along E0 = ˆx, and measured at points along θ ∈ [−π, π] and

ϕ = 0. These are compared against the Mie series solution.

The first result is the RCS of a sphere meshed into 20,480 patches and 30,720 edges with

an average edge length of about λ1/26.6 with an incident plane wave at 300 MHz. A three-level

MLFMA tree is used for this result. Figure 4.1a shows excellent agreement between the

formulations with almost no noticeable disagreement at any θ.

The second result is the RCS of a sphere meshed into 103,680 patches and 155,520 edges

with an average edge length of about λ1/35.7 with an incident wave at 500 MHz. A four-level

MLFMA tree is used for this result. We see in Figure 4.1b excellent agreement between all of

68

/20/21001020RCS [dBm²]MüllerE-DFIEA-DPIEMie/20/220100102030RCS [dBm²]MüllerE-DFIEA-DPIEMie/20/2100102030RCS [dBm²]MüllerE-DFIEA-DPIEMie(a) The NASA almond geometry with incident
and observation angles.

(b) The RCS at 900 MHz.

Figure 4.2 The dielectric NASA almond with µr2 = 1 and ϵr2 = 1.5.

the formulations with only slight disagreement with the DPIE formulation for certain RCS

valleys.

The third result is the RCS of a sphere meshed into 184,320 patches and 276,480 edges

with an average edge length of about λ1/23.8 with an incident wave at 1 GHz. A five-level

MLFMA tree is used for this result. We see in Figure 4.1c excellent agreement between all of

the formulations except the DPIE formulation, which still captures the general shape but

misses the peaks of some side lobes by up to 5 dB m2 for |θ| > π/8.

In Figure 4.1, we have shown that the accuracy is not impacted by the MLFMA acceleration

and have identified that the RCS from the A − ϕ-DPIE starts to deviate at higher frequencies.

4.4.1.2 Arbitrary Objects

Next, we analyze scattering from a dielectric NASA almond, as seen in Figure 4.2a. The

NASA almond is meshed into 72,576 patches and 108,864 edges with an average edge length

of about λ1/25.6 with an incident 900 MHz plane wave with ˆκ = ˆz, polarization E0 = ˆy,

and measured at points along θ ∈ [−π, π] and ϕ = π

2 . A five-level MLFMA tree is used

for this result. Figure 4.2b shows excellent agreement between the formulations with slight

disagreement with the DPIE formulation for certain RCS valleys.

Finally, we analyze scattering from a dielectric arrowhead, as seen in Figure 4.3a. This is

a challenging geometry in that it has sharp tips and edges. The arrowhead is meshed into

69

/20/230201001020RCS [dBm²]MüllerE-DFIEA-DPIE(a) The arrowhead geometry with incident and
observation angles.

(b) The RCS of the dielectric arrowhead at
300 MHz.

Figure 4.3 The dielectric arrowhead with µr2 = 1 and ϵr2 = 1.5.

23,024 patches and 34,536 edges with an average edge length of about λ1/15.6 with an incident

300 MHz plane wave directed along ˆκ = ˆz, polarization E0 = ˆy, and measured at points

along θ ∈ [−π, π] and ϕ = π

2 . A five-level MLFMA tree is used for this result.

We see in Figure 4.3b excellent agreement between all of the formulations with almost no

noticeable disagreement at any θ.

The results from the analysis of the sphere in Figure 4.1, almond in Figure 4.2, and arrow

in Figure 4.3 are for geometries that result in systems that are much larger than those from

Chapter 3. The results presented in this chapter would be too large to fit in memory on the

machines we were utilizing if solved without MLFMA acceleration. In addition, they were

solved significantly faster than some of the smaller problems from Chapter 3. The following

section will explore a timing study for our implementation.

4.4.2 Timing

For the timing study, we use a 1 m × 1 m × 1 m cube with a variety of meshes discretized

with different edge lengths. The frequency is chosen for each mesh so that the edge lengths

are approximately λ1/15. Only the outer region’s discretized system ⟨F , Z1 ◦ F ⟩ y is timed

for the various formulations because all other matrices involved are sparse. The number of

unknowns N for a given mesh depends on the formulation: 2Nv, 2Ns + 2Ns, and 2Nv + 4Ns

70

/20/220100102030RCS [dBm²]MüllerE-DFIEA-DPIEFigure 4.4 Timings for direct and MLFMA accelerated matvecs.

for the M¨uller, DFIE, and DPIE formulations respectively. As is seen in Figure 4.4, the direct

matrix vector computation scales with O (N 2) and we confirm that the MLFMA accelerated

computation scales with O (N log N ). The cost of MLFMA can start paying dividends for

relatively small meshes because the cross-over point occurs on the order of thousands of

unknowns. However, the exact cross-over point of around 3 × 104 unknowns is not meaningful.

It is somewhat of an upper bound because many standard MLFMA optimizations have

been omitted in our implementation. The grand takeaway is that the additional operators

and basis function types of the DFIE and DPIE formulations do not change the scaling

properties nor significantly increase the matvec computation time compared to standard

MLFMA accelerated formulations.

4.5 Conclusion

This chapter presents how to accelerate both the DFIE and DPIE with MLFMA. As

is observed in canonical and non-canonical geometries, the discrete piecewise systems are

accurate, and MLFMA successfully accelerates the formulations with only minor modifications

to existing implementations.

71

103104105106Number of Unknowns103102101100101Time [s]MüllerDFIEDPIEFMMFullCHAPTER 5

CONCLUSION

To briefly summarize, in this thesis, a generalized formulation process for the direct method

is developed and applied to both classic and several new formulations in Chapter 2, details

required for discretizing the operators involved in all formulations with spherical harmonics

and piecewise basis functions with results detailing the characteristics of each formulation

were provided in Chapter 3. The necessary modifications needed to implement all of the

formulations with MLFMA with RCS and timing results are given in Chapter 4.

This thesis pushes the state-of-the-art forward by (1) generalizing the direct method

formulation process and utilizing it to present old and new formulations, (2) explicitly

describing how to implement all of the formulations with analytic basis sets, enabling

eigenvalue and condition number analysis in a unified manner, (3) developing the required

integrals for implementing all operators used in many formulations on piecewise basis and

testing functions, and (4) demonstrating how to accelerate these formulations using a common

MLFMA framework.

5.1 Future Work

Some areas of further development include but are not limited to

• additional excitations such as point sources

• impedance boundary conditions

• extending the framework to support multiply connected and open geometries

• more advanced singularity treatment for the piecewise basis sets

• alternate basis functions and non-Galerkin testing such that the testing functions lie in

the dual space of the range of each operator

• higher order piecewise basis functions

• hybrid Multilevel Accelerated Cartesian Expansion (MLACE)-MLFMA acceleration for

wide-band performance

• FEBI with DFIE and DPIE formulations

72

Much of this work can be accomplished abstractly, as was done in this thesis, such that

comparisons to existing formulations and methods are always apples-to-apples and share as

much logic and implementation as possible. Future novel formulations will also benefit from

this and can be more quickly matured with less effort.

73

BIBLIOGRAPHY

[Adrian et al., 2019] Adrian, S. B., Andriulli, F. P., and Eibert, T. F. (2019). On a refinement-
free calder´on multiplicative preconditioner for the electric field integral equation. Journal
of Computational Physics, 376:1232–1252.

[Alsnayyan et al., 2020] Alsnayyan, A., Li, J., Hughey, S., Diaz, A., and Shanker, B. (2020).
Efficient isogeometric boundary element method for analysis of acoustic scattering from
rigid bodies. The Journal of the Acoustical Society of America, 147(5):3275–3284.

[Andriulli et al., 2008] Andriulli, F. P., Cools, K., Bagci, H., Olyslager, F., Buffa, A., Chris-
tiansen, S., and Michielssen, E. (2008). A multiplicative calderon preconditioner for
the electric field integral equation. IEEE Transactions on Antennas and Propagation,
56(8):2398–2412.

[Antoine et al., 2006] Antoine, X., Darbas, M., and Lu, Y. Y. (2006). Computer methods in

applied mechanics and engineering. Arch. Elektr. Ubertrag, 195(33-36):4060–4074.

[Balanis, 2012] Balanis, C. A. (2012). Advanced engineering electromagnetics. John Wiley &

Sons.

[Botha, 2013] Botha, M. M. (2013). A family of augmented duffy transformations for near-
singularity cancellation quadrature. IEEE Transactions on Antennas and Propagation,
61(6):3123–3134.

[Burton and Kashyap, 1995] Burton, M. and Kashyap, S. (1995). A study of a recent mo-
ment method algorithm that is accurate to very low frequencies. Applied Computational
Electromagnetics Society Journal, 10(3):58–68.

[Chew, 2014] Chew, W. C. (2014). Vector potential electromagnetics with generalized gauge
for inhomogeneous media: Formulation. Progress In Electromagnetics Research, 149:69–84.

[Coifman et al., 1993] Coifman, R., Rokhlin, V., and Wandzura, S. (1993). The fast multipole
method for the wave equation: A pedestrian prescription. IEEE Antennas and Propagation
magazine, 35(3):7–12.

[Cools et al., 2009a] Cools, K., Andriulli, F. P., Olyslager, F., and Michielssen, E. (2009a).
Nullspaces of mfie and calder ´On preconditioned efie operators applied to toroidal surfaces.
IEEE Transactions on Antennas and Propagation, 57(10):3205–3215.

[Cools et al., 2009b] Cools, K., Andriulli, F. P., Olyslager, F., and Michielssen, E. (2009b).
Time domain calder ´On identities and their application to the integral equation analysis of
scattering by pec objects part i: Preconditioning. IEEE Transactions on Antennas and
Propagation, 57(8):2352–2364.

74

[Costabel, 1991] Costabel, M. (1991). A coercive bilinear form for maxwell’s equations.

Journal of Mathematical Analysis and Applications, 157(2):527–541.

[Costabel and Stephan, 1988] Costabel, M. and Stephan, E. P. (1988). Strongly elliptic
boundary integral equations for electromagnetic transmission problems. Proceedings of the
Royal Society of Edinburgh: Section A Mathematics, 109(3-4):271–296.

[Dault and Shanker, 2015] Dault, D. and Shanker, B. (2015). An interior penalty method
for the generalized method of moments. IEEE Transactions on Antennas and Propagation,
63(8):3561–3568.

[Dault and Shanker, 2015] Dault, D. and Shanker, B. (2015). A mixed potential mlfma for
higher order moment methods with application to the generalized method of moments.
IEEE Transactions on Antennas and Propagation, 64(2):650–662.

[Dely, 2020] Dely, A. (2020). Computational strategies for impedance boundary condition
integral equations in frequency and time domains. PhD thesis, Ecole nationale sup´erieure
Mines-T´el´ecom Atlantique.

[Epstein and Greengard, 2010] Epstein, C. L. and Greengard, L. (2010). Debye sources and
the numerical solution of the time harmonic maxwell equations. Communications on Pure
and Applied Mathematics, 63(4):413–463.

[Fink et al., 2008] Fink, P. W., Wilton, D. R., and Khayat, M. A. (2008). Simple and
efficient numerical evaluation of near-hypersingular integrals. IEEE Antennas and Wireless
Propagation Letters, 7:469–472.

[Fu et al., 2017] Fu, X., Li, J., Jiang, L. J., and Shanker, B. (2017). Generalized debye
sources-based efie solver on subdivision surfaces. IEEE Transactions on Antennas and
Propagation, 65(10):5376–5386.

[Glisson, 1984] Glisson, A. (1984). An integral equation for electromagnetic scattering
from homogeneous dielectric bodies. IEEE Transactions on Antennas and Propagation,
32(2):173–175.

[Graglia et al., 2013] Graglia, R. D., Peterson, A. F., and Matekovits, L. (2013). Singular,
hierarchical scalar basis functions for triangular cells. IEEE Transactions on Antennas
and Propagation, 61(7):3674–3692.

[Hawkins et al., 2022] Hawkins, J., Baumann, L., Aktulga, H., Dault, D., and Shanker, B.
(2022). Analytic preconditioners for decoupled potential integral equations for wideband
analysis of scattering from pec objects. arXiv preprint arXiv:2211.10560.

[Helsing et al., 2020] Helsing, J., Karlsson, A., and Ros´en, A. (2020). Comparison of integral

equations for the maxwell transmission problem with general permittivities.

75

[Hodges and Rahmat-Samii, 1997] Hodges, R. E. and Rahmat-Samii, Y. (1997). The eval-
uation of mfie integrals with the use of vector triangle basis functions. Microwave and
Optical Technology Letters, 14(1):9–14.

[Hsiao and Kleinman, 1997] Hsiao, G. C. and Kleinman, R. E. (1997). Mathematical founda-
tions for error estimation in numerical solutions of integral equations in electromagnetics.
IEEE Transactions on Antennas and Propagation, 45(3):316–328.

[Hughey et al., 2018] Hughey, S., Aktulga, H., Vikram, M., Lu, M., Shanker, B., and
Michielssen, E. (2018). Parallel wideband mlfma for analysis of electrically large, nonuniform,
multiscale structures. IEEE Transactions on Antennas and Propagation, 67(2):1094–1107.

[Kress and Roach, 1978] Kress, R. and Roach, G. F. (1978). Transmission problems for the

helmholtz equation. Journal of Mathematical Physics, 19(6):1433–1437.

[Li et al., 2014] Li, J., Dault, D., Nair, N., and Shanker, B. (2014). Analysis of scattering
from complex dielectric objects using the generalized method of moments. JOSA A,
31(11):2346–2355.

[Li et al., 2017] Li, J., Fu, X., and Shanker, B. (2017). Potential integral equations in
electromagnetics. In 2017 IEEE International Symposium on Antennas and Propagation
USNC/URSI National Radio Science Meeting.

[Li et al., 2019] Li, J., Fu, X., and Shanker, B. (2019). Decoupled potential integral equations
for electromagnetic scattering from dielectric objects. IEEE Transactions on Antennas
and Propagation, 67(3):1729–1739.

[Liu et al., 2015] Liu, Q. S., Sun, S., and Chew, W. C. (2015). An integral equation method
based on vector and scalar potential formulations. In 2015 IEEE International Symposium
on Antennas and Propagation USNC/URSI National Radio Science Meeting, pages 744–745.

[Marx, 1984] Marx, E. (1984). Integral equation for scattering by a dielectric. IEEE Trans-

actions on Antennas and Propagation, 32(2):166–172.

[M¨uller, 1969] M¨uller, C. (1969). Electromagnetic Waves in a Homogeneous Medium. Springer

Berlin Heidelberg, Berlin, Heidelberg.

[Nair and Shanker, 2011] Nair, N. V. and Shanker, B. (2011). Generalized method of mo-
ments: A novel discretization technique for integral equations. IEEE transactions on
antennas and propagation, 59(6):2280–2293.

[Poggio and Miller, 1973] Poggio, A. and Miller, E. (1973). Chapter 4 - integral equation
solutions of three-dimensional scattering problems. In Mittra, R., editor, Computer Tech-
niques for Electromagnetics, International Series of Monographs in Electrical Engineering,
pages 159 – 264. Pergamon.

76

[Qian and Chew, 2009] Qian, Z. and Chew, W. C. (2009). Fast full-wave surface integral
equation solver for multiscale structure modeling. IEEE Transactions on Antennas and
Propagation, 57(11):3594–3601.

[Qian and Chew, 2008] Qian, Z. G. and Chew, W. C. (2008). A quantitative study on the low
frequency breakdown of efie. Microwave and Optical Technology Letters, 50(5):1159–1162.

[Rao et al., 1982] Rao, S., Wilton, D., and Glisson, A. (1982). Electromagnetic scattering by
surfaces of arbitrary shape. IEEE Transactions on Antennas and Propagation, 30(3):409–
418.

[Roth and Chew, 2018] Roth, T. E. and Chew, W. C. (2018). Development of stable A-ϕ
time-domain integral equations for multiscale electromagnetics. IEEE Journal on Multiscale
and Multiphysics Computational Techniques, 3:255–265.

[Roth and Chew, 2020] Roth, T. E. and Chew, W. C. (2020). Radiation gauge potential-
based time domain integral equations for penetrable regions. Progress In Electromagnetics
Research, 168:73–86.

[Sarvas, 2003] Sarvas, J. (2003). Performing interpolation and anterpolation entirely by
fast fourier transform in the 3-d multilevel fast multipole algorithm. SIAM Journal on
Numerical Analysis, 41(6):2180–2196.

[Song et al., 1997] Song, J., Lu, C.-C., and Chew, W. C. (1997). Multilevel fast multipole
algorithm for electromagnetic scattering by large complex objects. IEEE transactions on
antennas and propagation, 45(10):1488–1493.

[Song and Chew, 1995] Song, J. M. and Chew, W. C. (1995). Multilevel fast-multipole
algorithm for solving combined field integral equations of electromagnetic scattering.
Microwave and optical technology letters, 10(1):14–19.

[Stratton, 2015] Stratton, J. A. (2015). Spherical Waves, chapter VII, pages 392–423. John

Wiley & Sons, Ltd.

[Terai, 1980] Terai, T. (1980). On calculation of sound fields around three dimensional objects

by integral equation methods. Journal of Sound and Vibration, 69(1):71 – 100.

[Tihon and Craeye, 2018] Tihon, D. and Craeye, C. (2018). All-analytical evaluation of the
singular integrals involved in the method of moments. IEEE Transactions on Antennas
and Propagation, 66(4):1925–1936.

[Vecchi, 1999] Vecchi, G. (1999). Loop-star decomposition of basis functions in the discretiza-

tion of the efie. IEEE Transactions on Antennas and Propagation, 47(2):339–346.

77

[Vico et al., 2016] Vico, F., Ferrando, M., Greengard, L., and Gimbutas, Z. (2016). The
decoupled potential integral equation for time-harmonic electromagnetic scattering. Com-
munications on Pure and Applied Mathematics, 69(4):771–812.

[Vico et al., 2015] Vico, F., Ferrando-Bataller, M., Jim´enez, T. B., and Berenguer, A. (2015).
A high order locally corrected nystr¨om implementation of the decoupled potential integral
equation. In 2015 9th European Conference on Antennas and Propagation (EuCAP), pages
1–4.

[Vico et al., 2018] Vico, F., Greengard, L., and Ferrando, M. (2018). Decoupled field in-
tegral equations for electromagnetic scattering from homogeneous penetrable obstacles.
Communications in Partial Differential Equations, 43(2):159–184.

[Vikram et al., 2009] Vikram, M., Huang, H., Shanker, B., and Van, T. (2009). A novel
wideband fmm for fast integral equation solution of multiscale problems in electromagnetics.
IEEE Transactions on Antennas and Propagation, 57(7):2094–2104.

[Wilcox, 1956] Wilcox, C. H. (1956). An expansion theorem for electromagnetic fields.

Communications on Pure and Applied Mathematics, 9(2):115–134.

[Wilde, 1987] Wilde, P. (1987). Transmission problems for the vector helmholtz equation.
Proceedings of the Royal Society of Edinburgh: Section A Mathematics, 105(1):61–76.

[Wilton et al., 1984] Wilton, D., Rao, S., Glisson, A., Schaubert, D., Al-Bundak, O., and
Butler, C. (1984). Potential integrals for uniform and linear source distributions on
polygonal and polyhedral domains. IEEE Transactions on Antennas and Propagation,
32(3):276–281.

[Wilton and Glisson, 1981] Wilton, D. R. and Glisson, A. W. (1981). On improving the
electric field integral equation at low frequencies. In URSI National Radio Science Meeting
Digest, page 24.

[Wu et al., 1995] Wu, W., Glisson, A. W., and Kajfez, D. (1995). A study of two numerial
solution procedures for the electric field integral equation at low frequency. Applied
Computational Electromagnetics Society Journal, 10(3):69–80.

[Yan et al., 2010] Yan, S., Jin, J., and Nie, Z. (2010). Calder´on preconditioner: From
efie and mfie to n-m¨uller equations. IEEE Transactions on Antennas and Propagation,
58(12):4105–4110.

[Yla-Oijala and Taskinen, 2005] Yla-Oijala, P. and Taskinen, M. (2005). Well-conditioned
muller formulation for electromagnetic scattering by dielectric objects. IEEE Transactions
on Antennas and Propagation, 53(10):3316–3323.

78

[Yl¨a-Oijala et al., 2005] Yl¨a-Oijala, P., Taskinen, M., and Sarvas, J. (2005). Surface integral
equation method for general composite metallic and dielectric structures with junctions.
Progress In Electromagnetics Research, 52:81–108.

[Yuan and Shanker, 2005] Yuan, J. and Shanker, B. (2005). Analysis of the spectrum of the
single integral equation for scattering from dielectric objects. In 2005 IEEE Antennas and
Propagation Society International Symposium, volume 2B, pages 93–96 vol. 2B.

79