INTEGRAL EQUATIONS IN COMPUTATIONAL ELECTROMAGNETICS:

FORMULATIONS, PROPERTIES AND ISOGEOMETRIC ANALYSIS

By

Jie Li

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Electrical Engineering – Doctor of Philosophy

2018

ABSTRACT

INTEGRAL EQUATIONS IN COMPUTATIONAL ELECTROMAGNETICS:

FORMULATIONS, PROPERTIES AND ISOGEOMETRIC ANALYSIS

By

Jie Li

Computational electromagnetics (CEM) provides numerical methods to simulate electro-

magnetic waves interacting with its environment. Boundary integral equation (BIE) based

methods, that solve the Maxwell’s equations in the homogeneous or piecewise homogeneous

medium, are both efficient and accurate, especially for scattering and radiation problems.

Development and analysis electromagnetic BIEs has been a very active topic in CEM re-

search.

Indeed, there are still many open problems that need to be addressed or further

studied. A short and important list includes (1) closed-form or quasi-analytical solutions

to time-domain integral equations, (2) catastrophic cancellations at low frequencies, (3) ill-

conditioning due to high mesh density, multi-scale discretization, and growing electrical size,

and (4) lack of flexibility due to re-meshing when increasing number of forward numerical

simulations are involved in the electromagnetic design process. This dissertation will address

those several aspects of boundary integral equations in computational electromagnetics.

The first contribution of the dissertation is to construct quasi-analytical solutions to

time-dependent boundary integral equations using a direct approach. Direct inverse Fourier

transform of the time-harmonic solutions is not stable due to the non-existence of the in-

verse Fourier transform of spherical Hankel functions. Using new addition theorems for the

time-domain Green’s function and dyadic Green’s functions, time-domain integral equations

governing transient scattering problems of spherical objects are solved directly and stably for

the first time. Additional, the direct time-dependent solutions, together with the newly pro-

posed time-domain dyadic Green’s functions, can enrich the time-domain spherical multipole

theory.

The second contribution is to create a novel method of moments (MoM) framework to

solve electromagnetic boundary integral equation on subdivision surfaces. The aim is to

avoid the meshing and re-meshing stages to accelerate the design process when the geometry

needs to be updated. Two schemes to construct basis functions on the subdivision surface

have been explored. One is to use the div-conforming basis function, and the other one is to

create a rigorous iso-geometric approach based on the subdivision basis function with better

smoothness properties. This new framework provides us better accuracy, more stability and

high flexibility.

The third contribution is a new stable integral equation formulation to avoid catastrophic

cancellations due to low-frequency breakdown or dense-mesh breakdown. Many of the con-

ventional integral equations and their associated post-processing operations suffer from nu-

merical catastrophic cancellations, which can lead to ill-conditioning of the linear systems

or serious accuracy problems. Examples includes low-frequency breakdown and dense mesh

breakdown. Another instability may come from nontrivial null spaces of involving integral

operators that might be related with spurious resonance or topology breakdown. This disser-

tation presents several sets of new boundary integral equations and studies their analytical

properties. The first proposed formulation leads to the scalar boundary integral equations

where only scalar unknowns are involved. Besides the requirements of gaining more stability

and better conditioning in the resulting linear systems, multi-physics simulation is another

driving force for new formulations. Scalar and vector potentials (rather than electromagnetic

field) based formulation have been studied for this purpose.

Those new contributions focus on different stages of boundary integral equations in an

almost independent manner, e.g.

isogeometric analysis framework can be used to solve

different boundary integral equations, and the time-dependent solutions to integral equations

from different formulations can be achieved through the same methodology proposed.

Copyright by
JIE LI
2018

To my family

v

ACKNOWLEDGEMENTS

First, I would like to give sincere thanks to my advisor Prof. Shanker Balasubramaniam for

the continuous guidance, encouragement and financial support. Working with such a great

advisor and researcher has always been amazing and rewarding. Without those many fruitful

discussions with him, the work in the thesis would not be possible.

Next, I thank Profs. Rothwell, Piermarocchi and Albrecht for serving in my committee.

I must also give my thanks to the faculty members in the electomagnetics research group

at MSU for developing great EM courses and creating such a flexible research environment.

Special thanks should also go to Profs. Piermarocchi and Tong for collaborations on projects

during the PhD study.

I owe thanks to group members. I want to thank Dr. Dault for the collaboration and

exciting discussions on many topics and ideas. I also appreciate the opportunity to work with

Drs. Nair, Meierbachtol, Baczewski and Pray in the same lab. With a special mention of my

current graduate fellows Steve Hughey, Connor Glosser, Zane Crawford, Scott O’Connor,

Luke Baumann and Abdel Alsnayyan, I’m so thankful to them for support and help in many

aspects. Those delightful conversations and fruitful discussions have been always valuable

to my graduate life and research. I also thank Xin Fu, Beibei Liu and Rundong Zhao for

collaborations on topics in the thesis.

I would like to give my appreciation to lots of people who have supported and helped me

in the past 5 years. They are Yiqun Yang, Pedro Nariyoshi, Xiaofeng Zhao, Bin Fan, Shichen

Zhang, Qianwei Jiang and Sichao Wang, followed by many others in the MSU community

and a long list of family friends in the Greater Lansing area.

Finally, I thank my family for understanding and love. Without the support and encour-

agement from my parents and sister, I wouldn’t be where I am today. Special thanks go to

my wife for her sacrifices and continuous love. I feel so grateful to have her with me for these

years and to have our remarkable daughters Sophia and Mila.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 1

1.1 Maxwell’s Equations and Boundary Conditions
1.2 Overview of Computational Electromagnetics

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
1.2.1 Time-domain and Frequency-domain Methods . . . . . . . . . . . . .
1.2.2 Differential and Integral Equation-based Methods . . . . . . . . . . .
1.3 Topics that are Studied in this Dissertation . . . . . . . . . . . . . . . . . . .

x

xi

1
1
2
3
4
6

8
CHAPTER 2 ANALYTICAL SOLUTIONS TO INTEGRAL EQUATIONS . . . . .
8
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.3 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.3.1 Natural Spatial basis: Vector spherical harmonics . . . . . . . . . . .
2.3.2
13
Spherical expansion of Time-domain Dyadic Green’s function . . . . .
2.3.3 Volterra Integral kernels and Reduced Time-domain Integral Equations 17
19
21
23
27

2.4 Time-dependent Debye-Mie series solution . . . . . . . . . . . . . . . . . . .
2.5 Eigen-spectrum analysis of the MOT system . . . . . . . . . . . . . . . . . .
2.6 Numerical Examples
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 3 NUMERICAL SOLUTIONS ON HIGHER ORDER GEOMETRY

REPRESENTATIONS . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1
3.2 Subdivision Surfaces
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Formulation and Discretization . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1
Integral Equation for PEC . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Div-conforming Basis on Smooth Subdivision Surface . . . . . . . . .
3.4 Numerical Examples
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 4

ISO-GEOMETRIC ANALYSIS ENHANCED METHOD OF MO-
MENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1
4.2
Integral Equations for Electromagnetic Scattering . . . . . . . . . . . . . . .
4.3 Current Representation and Field Solvers on Subdivision Surfaces . . . . . .
4.3.1
. . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Field Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Low-frequency Stable EFIE . . . . . . . . . . . . . . . . . . . . . . .
4.3.4 Calderón Preconditioner . . . . . . . . . . . . . . . . . . . . . . . . .

Iso-geometric Basis Sets

32
32
33
37
37
38
38
39

41
41
44
46
46
52
54
55

vii

4.4 Numerical Examples

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Accuracy of IGA-MoM . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2
Scattering from Structures at Regular Frequency . . . . . . . . . . .
4.4.3 Conditioning of iso-geometric System at Low-frequency . . . . . . . .
4.4.4 Examples with Multi-scale Mesh . . . . . . . . . . . . . . . . . . . . .
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 5 SCALAR INTEGRAL EQUATIONS FOR ELECTROMAGNETICS .
5.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Charge based Representation . . . . . . . . . . . . . . . . . . . . . .
Scalar Electric Field Integral Equation . . . . . . . . . . . . . . . . .
5.2.3
5.2.4
Scalar Magnetic Field Integral Equation . . . . . . . . . . . . . . . .
5.2.5 Combined Integral Equation . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
Iso-geometric Analysis on Subdivision Surface . . . . . . . . . . . . .
Scalar Basis Function and Linear System Setup . . . . . . . . . . . .
System Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1
Spectral Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2 Numerical Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.3 Method of Moments Solution Procedure

5.4 Numerical Examples

5.3.1
5.3.2
5.3.3

57
57
58
60
61
61

73
73
76
76
76
78
79
80
81
81
82
85
85
86
89
91

6.3.1
6.3.2 Representation of Scattering Potential

CHAPTER 6 DECOUPLED POTENTIAL BASED FORMULATIONS . . . . . . .
6.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Representations of the Decoupled Potential . . . . . . . . . . . . . . . . . . .
Scalar and Vector Potentials . . . . . . . . . . . . . . . . . . . . . . .

94
94
96
98
98
. . . . . . . . . . . . . . . . . 100
6.4 Decoupled Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.4.1 Decoupled Boundary Conditions for the PEC case . . . . . . . . . . . 101
6.4.2 Decoupled Boundary Conditions for Dielectric Case . . . . . . . . . . 102
6.5 Formulation of Decoupled Potential Integral Equations . . . . . . . . . . . . 105
6.5.1
. . . . . . . . . . . . . . . 105
6.5.2 Vector Potential Integral Equation (VPIE) . . . . . . . . . . . . . . . 107
6.6 Analytical Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Stability Properties of SPIE . . . . . . . . . . . . . . . . . . . . . . . 112
Stability Properties of VPIE . . . . . . . . . . . . . . . . . . . . . . . 112
6.7 Perfectly Electrical Conductor Case . . . . . . . . . . . . . . . . . . . . . . . 114
6.8 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.9 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

Scalar Potential Integral Equation (SPIE)

6.6.1
6.6.2

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
APPENDIX A SPHERICAL EXPANSIONS . . . . . . . . . . . . . . . . . . . 123

viii

APPENDIX B BOUNDARY INTEGRAL OPERATORS . . . . . . . . . . . . 126

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

ix

LIST OF TABLES

Table 4.1: Rel. L2 error in the reconstructed surface currents density

. . . . . . . .

58

x

LIST OF FIGURES

Figure 2.1: Time-domain coefficients for mode Ψ1
3

. . . . . . . . . . . . . . . . . . .

Figure 2.2: Time-domain coefficients for mode Φ1
3

. . . . . . . . . . . . . . . . . . .

Figure 2.3: Frequency-domain coefficients for modes Ψ1

3,compared to two
curves (denoted by FD-Mie-psi and FD-Mie-phi respectively) with fre-
quency domain Mie seires approach.

. . . . . . . . . . . . . . . . . . . .

3 and Φ1

Figure 2.4: Time-domain coefficients for mode Ψ1

30 . . . . . . . . . . . . . . . . . . .

Figure 2.5: Time-domain coefficients for mode Φ1

30 . . . . . . . . . . . . . . . . . . .

Figure 2.6: Frequency-domain coefficients for mode Ψ1

30„compared
to two curves (denoted by FD-Mie-psi and FD-Mie-phi respectively)
with frequency domain Mie seires approach.

30 and mode Φ1

. . . . . . . . . . . . . . . .

Figure 2.7: Late-time stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.8: Eigenvalues of MOT system of (a) kernel 1, (b) kernel 2, (c) kernel 3

and (d) kernel 4.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.1:

Irregular vertex and triangle (a) An irregular triangle (vertex 1 is valence-
7), (b)subdivision once

. . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.2: Local averaging during refinement.

. . . . . . . . . . . . . . . . . . . . .

Figure 3.3: Subdivision surface (a) primal mesh , (b) one subdivision, (c) two sub-

divisions and (d) three subdivisions. . . . . . . . . . . . . . . . . . . . . .

Figure 3.4: The parameter domain within an irregular triangle is recursively parti-

tioned into an infinite number of subdomains.

. . . . . . . . . . . . . . .

Figure 3.5: Convergence of the relative error in far field for the sphere . . . . . . . .

Figure 3.6: Convergence of the relative error in far field for the torus . . . . . . . . .

Figure 4.1: Scalar basis function associated with a vertex with valence of 8. (a)

Basis function, (b) Surface Laplacian of the basis function . . . . . . . .

xi

24

25

26

27

28

29

30

31

34

35

36

37

39

40

48

Figure 4.2: Basis function associated with a vertex with valence of 8.

(a) non-

solenoidal type, (b) solenoidal type . . . . . . . . . . . . . . . . . . . . .

Figure 4.3: Magnitude of surface currents density on the sphere: (a) real part and

(b) imaginary part

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.4: Radar cross section of the sphere (φ = 0 cut) . . . . . . . . . . . . . . . .

Figure 4.5: Pointwise relative error (a) real part, (b) imaginary part with 20484 DoFs

Figure 4.6: Magnitude of surface current density on the truncated cone: (a) real

part and (b) imaginary part . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.7: Radar cross section of the truncated cone (φ = 0 cut) . . . . . . . . . . .

Figure 4.8: Magnitude of surface current density on the composite structure: (a)

real part and (b) imaginary part

. . . . . . . . . . . . . . . . . . . . . .

Figure 4.9: Radar cross section of the composite structure (φ = 90 cut)

. . . . . . .

Figure 4.10: Magnitude of surface current density on the air plane model: (a) real

part and (b) imaginary part . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.11: Radar cross section of the air plane model (φ = 90 cut) . . . . . . . . . .

Figure 4.12: Condition number at different frequencies
Figure 4.13: Number of GMRES iterations to converge to a tolerance of 1.0 × 10−6
.
Figure 4.14: Number of GMRES iterations to converge to a tolerance of 1.0 × 10−10 .

. . . . . . . . . . . . . . . . .

Figure 4.15: Comparison of RCS data obtained at 1Hz with analytical solutions

. . .

51

63

64

65

66

66

67

67

68

69

69

70

70

71

Figure 4.16: Multiscale control mesh: (a) the whole mesh and (b) locally refined region 71

Figure 4.17: Convergence of the iterative solver

. . . . . . . . . . . . . . . . . . . . .

Figure 5.1: The support of the scalar basis: two ring of a valence-5 vertex . . . . . .

Figure 5.2: Scalar basis associated with a vertex . . . . . . . . . . . . . . . . . . . .

Figure 5.3: Condition number versus frequency . . . . . . . . . . . . . . . . . . . . .

Figure 5.4: Eigen-condition number versus frequency . . . . . . . . . . . . . . . . . .

72

83

84

87

87

xii

Figure 5.5: Eigenvalues of the system of (a) scalar EFIE, (b) scalar MFIE and (c)

scalar CFIE with α = 0.5/η . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.6: Singular values of the system of the scalar EFIE, the scalar MFIE and

the scalar CFIE with α = 0.5/η . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.7: Radar cross sections of the sphere . . . . . . . . . . . . . . . . . . . . . .

Figure 5.8: Residual of GMRES Iterations, 600Mhz

. . . . . . . . . . . . . . . . . .

Figure 5.9: A plane model: real part of the magnitude of the surface current density

Figure 5.10: Residual of GMRES iterations at different frequencies . . . . . . . . . . .

88

89

90

91

92

93

Figure 6.1: Real and imaginary parts of eigenvalues at 1 Hz . . . . . . . . . . . . . . 117
Figure 6.2: Real and imaginary parts of eigenvalues at 1 × 107 Hz

. . . . . . . . . . 117

Figure 6.3: Condition number of VPIE formulation versus frequency . . . . . . . . . 118

Figure 6.4: Condition number of Müller formulation versus frequency . . . . . . . . . 118

Figure 6.5: Real part of the coefficients for Ψ1

30 mode in the tangential component

of the electric field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

Figure 6.6: Real part of the coefficients for Φ1

30 mode in the tangential component

of the electric field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

xiii

CHAPTER 1

INTRODUCTION

1.1 Maxwell’s Equations and Boundary Conditions

Electromagnetic waves have been very useful in electrical engineering, playing a signif-

icant role in improving people’s lives, industry technologies and even scientific research.

Applications of electromagnetics include communication, radar detection, remote sensing,

bio-medical imaging, lithography and so on, covering a wide range of civil and military ac-

tivities. Engineering with electromagnetic fields involves design and evaluation of complex

systems, where modeling interactions between the system and electromagnetic waves is an

inevitable step. The electromagnetic theory has developed significantly after the publica-

tion of Maxwell’s equations. To model the electromagnetic behaviors of devices or systems,

Maxwell equations have to solved.

In the three-dimensional space, time-dependent electromagnetic fields satisfy Maxwell’s

equations,

∇ × E = −∂B
∂t
∇ × H = J +
∇ · D = ρ
∇ · B = 0

∂D
∂t

(1.1)

where E and H denote electric and magnetic fields respectively; D and B denote electric

flux and magnetic flux respectively. In Maxwell’s equations, electric current density J and

electric charge density J can be considered as sources distributed in the whole space.

In the presence of material discontinuity, the following boundary conditions, regarding

1

the field and flux on both sides of the interface, have to be in play.

n21 × (E1 − E2) = 0
n21 × (H1 − H2) = Js
n21 · (D1 − D2) = ρs
n21 · (B1 − B2) = 0

(1.2)

where n21 represents the unit vector normal to the interface, pointing from medium 2 to
medium 1, and Js and ρs denote the surface current density and charge density, representing
source terms on the interface of discontinuity. It’s noted that they are not the only way to

introduce the exciting source, and another possible source could be the incident wave that

is distributed in the whole space. The latter is very common in scattering analysis.

After taking Fourier transform on both sides of (1.1), Maxwell’s equations in frequency

domain can be written as follows.

∇ × E = −jωB
∇ × H = J + jωD
∇ · D = ρ
∇ · B = 0

(1.3)

where the time-harmonic factor ejωt is assumed.

For linear medium, especially in narrow-band problems, the frequency-domain or time-

harmonic description allows easier discretization, with only spatial variation to be considered

at each frequency. This dissertation will mainly work on time-harmonic problems except

Chapter 2 where time-dependent problems are solved.

1.2 Overview of Computational Electromagnetics

In modern electrical engineering where the electromagnetic system is involved, compu-

tational modeling is almost one inevitable step simply because it’s usually more efficient

and cost-effective to carry out the evaluation task with the help of high-speed processors.

2

Electromagnetic modeling is becoming one important component in the computer-aided engi-

neering (CAE) for designing and manufacturing either consumer or industry products. That

drives the research of solving Maxwell’s equation in a fast computational manner, rather

than the experimental measurement approach. The subject of computational electromag-

netics is to study how to solve Maxwell’s equations, given the characterization of materials,

boundary conditions and the excitation information. Depending on the characteristics (ge-

ometry, material and excitation) of the problem, there are various approaches available to

solve Maxwell’s equations.

1.2.1 Time-domain and Frequency-domain Methods

In most real-life electromagnetic problems, solutions to Maxwell’s equations are time-dependent.

Therefore it seems very natural to solve the equations directly in the time domain. However,

for large amount of linear problems, the time-dependent descriptions can be transformed

into the frequency domain, where the time-harmonic response is to be sought. That’s the

reason why two versions of Maxwell’s equations given earlier are commonly used.

Time-domain methods discretize the governing equations both temporally and spatially,

in a more first-principle sense. Therefore the unknowns are both time and spatial depen-

dent, usually resulting more unknowns than their frequency-domain counterparts. In linear

problems, if the studied system’s parameters are well characterized for its time-dependence,

time-domain methods can produce the transient response which usually contains wide-band

information. Time-dependent simulations can be employed for multi-physics problems as

well, where non-linearity might show up when governing equations for two physics laws are

coupled.

On the other hand, the frequency-domain simulation is the preferable choice for charac-

terizing the frequency-domain response in a narrow-band scenario. One major assumption

for this type of methods is that the material parameters can be determined for the inter-

ested frequency. In frequency-domain methods, the time-harmonic factor is also implicitly

3

assumed on both sides of the equation, therefore no temporal discretization is involved at

all in the solving process.

It is straightforward that all the full-wave (including both analytical and numerical ap-

proaches) and asymptotic methods have both frequency-domain and time-domain versions.

Similarly, both differential and integral equation-based methods can be designed in either

frequency or time domain.

1.2.2 Differential and Integral Equation-based Methods

In full-wave electromagnetic simulations, one can choose to solve either the original set

of coupled Maxwell’s equations or the second order curl-curl equation. Typical numerical

methods include finite-difference method (FDM) and finite element method (FEM). Methods

like FDM and FEM, no matter in time or frequency domain, can be classified as differential

equation(DE)-based methods. These methods have the following advantages.

1. DE-based methods are more suitable to model very generalized problems, like those

with non-uniform materials.

2. The resulting discretized linear system is sparse, and the numerical implementation is

much easier.

3. Significant progress has been made on many theoretical aspects of DE-based methods,

such as linear system solver, convergence and error analysis.

Those features usually enable DE-based methods the first choice when modelling tasks are

planned. That also explains the fact that commercial codes implementing DE-based methods

are available.

On the other hand, disadvantages of DE based methods are listed as follows.

1. Artificial absorbing boundary condition has to be applied to truncate the open domain

problem.

4

2. The differential operator is numerically approximated which caused the phenomenon

of the numerical dispersion.

3. Volume mesh is involved, which explains why they’re not very efficient for some scat-

tering problems.

Compared to the fact that DE-based methods discretize the equation in the volume, an-

other type of methods based on boundary integral equations (BIE) only involve discretization

on the interfaces between two types of materials. The reduction in the dimensionality usually

decreases the number of unknowns for problems have different materials. More importantly,

BIE-based methods don’t suffer from the numerical dispersion error or need artificial ab-

sorbing boundaries. In electromagnetics, classical formulations for BIE include electric field

integral equation (EFIE), magnetic field integral equation (MFIE) and their linear combina-

tion (combined field integral equation, CFIE). The numerical method to solve those boundary

integral equations is called method of moments (MoM), which can be implemented in both

time and frequency domain. Surface integral equation has lots of applications in scattering

and radiation analysis, where the target domain is piece-wise homogeneous in terms of mate-

rial properties. Additionally, integral equation can be also used in hybrid methods to model

even more complicated systems, or to achieve better performance. Though some people

would argue otherwise, the integral equation is more involved than its differential equation

counterpart in modeling electromagnetics. Because of the aforementioned reasons, study on

integral equations by the community has been always intensive in the past several decades.

It’s worth noting that this section by no means serves the goal of giving detailed review on

computational electromagnetics. More information on the progress made in the community

can be found in the journal of IEEE transactions on Antennas and Propagation and similar

periodicals. This dissertation will focus on the boundary integral equations. In the next

section, main topics and contributions of this dissertation will be briefly introduced.

5

1.3 Topics that are Studied in this Dissertation

This dissertation will focus on most of the research results during my PhD study. Several

important aspects of boundary integral equations are covered, including (1) quasi-analytical

time-dependent solutions of scattering and propagation in spherical domains, (2) a novel dis-

cretization framework for boundary integral equations that supports efficient design-through-

analysis process, and (3) new formulations for boundary integral equations and their prop-

erties.

Chapter 2 presents the progress on the direct time-dependent spherical harmonics anal-

ysis in electromagnetics. More than one century ago, analytical time-harmonic solution

to scattering problem of spherical objects was first studied. Similar procedure leading to

closed-form full wave results in time-domain hasn’t been successful for a long time even the

frequency-domain Mie series approach was widely used. A detailed review on this topic in

the history was given in [1]. With the newly derived closed spherical expansion form of time-

domain dyadic Green’s functions, direct time-domain modeling of electromagnetic scattering

or radiation in the spherical domain is made possible in a quasi-analytical manner. The pro-

cedure will be introduced in this chapter, which was published in [2]. The method has been

applied for both acoustics and electromagnetics involving different boundary conditions [1–3].

Chapters 3 and 4 are devoted to the development of new numerical discretization schemes

of boundary integral equations. In Chapter 3, the subdivision surface, a very good geometry

representation scheme in computer graphics, is presented and directly used to construct the

method of moments. Div-conforming basis functions are used, which can achieve super-

convergence in the far field. Chapter 4 presents a new framework to discretize the electro-

magnetic integral equations. The iso-geometric analysis framework, which was first proposed

to flexibly aid the computer-aided engineering process in structure mechanics, is fully ex-

tended and implemented for electromagnetic scattering analysis with the help of subdivision

surfaces. This study on iso-geometric analysis on electromagnetic integral equations has not

been done before. Materials in Chapter 3 and 4 have been published in [4, 5].

6

In Chapters 5 and 6, two types of new boundary integral equations are formulated to

improve the stability. Chapter 5 studies scalar electromagnetic integral equations based on

the idea of applying the Helmholtz decomposition and introducing charge densities as the

only unknowns. Iso-geometric discretization is implemented without introducing any low-

frequency breakdown. Furthermore, the scalar magnetic field integral equation (sMFIE) is

even free from the dense-mesh breakdown due to the fact that it is the second kind integral

equation. Then in Chapter 6, a totally different formulation of boundary integral equations

is presented. The new formulation is based on decoupled scalar and vector potentials, and

in another words, the new formulation gives one integral equation for the scalar potential

and another one for the vector potential. The work in Chapter 6 extends the existing

decoupled potential formulation for perfect electric conductors to dielectric objects, using

a direct approach (vs.

indirect approach) where the unknowns have the physical meaning.

Analytical properties of the resulting integral equations are studied for the spherical object.

Numerical examples are given to illustrate the stability for low frequency and high resolution

modes (corresponding to dense discretization).

It is worth noting that each subsequent chapter is very independent from each other,

because each chapter corresponds to one paper and it has its own detailed literature review,

formulation and section of numerical examples. The notation in each chapter also follows

its own convention, with quite a few cross references (if any) without even breaking its

completeness.

7

CHAPTER 2

ANALYTICAL SOLUTIONS TO INTEGRAL EQUATIONS

2.1

Introduction

Debye-Mie series solution [6, 7] is one of the most useful tools in time-harmonic analy-

sis of electromagnetic (EM) scattering, and has found in extensive use in both optics and

electromagnetics [8, 9]. Their applications range from fields as diverse as light scattering

from small particles [10], to combustion [11], to analysis of antennas [12,13], to weather [14],

to validation and calibration computational and measurement setups, to photonics [15], to

biological applications [16].

However, analytical solutions to scattering from a sphere only exist in the Fourier domain.

Equivalent literature for transient analysis is very sparse. While it is indeed possible to

obtain transient scattering response via Fourier transforms, the challenge of obtaining a

stable direct transient solution has remained unsolved for several decades. While attempts

to arrive at a solution have been made via explicit Fourier transforms [17] of Bessel and

Hankel functions, these approaches are highly oscillatory and therefore, unstable. Aside

from an intellectual challenge, the development of methods to solve for scattering from

sphere can find applications in a number of areas; from scattering from collection of particles

to study of non-linear response of gain particles [18], to imaging [19], to dielectric resonator

antennas [20], and so on.

While direct time-domain Mie series like solutions for scattering do not exist, there has

been extensive interest in developing multipole theory for representations of radiated fields

due to spatially bounded sources. Time domain multipole theory has found applications

in both acoustics as well as electromagnetics [21–26], with applications that range from

radiation due to volumetric sources [27], to near field scanning [28], to near-to-far field

transforms [29]. However, multipole representations require higher order derivatives in time.

8

As a result, there is a paucity of numerical data and the treatment is largely theoretical.

This comment does not apply to fast solvers that use a Whittaker representation of the

radiated field [30], and there are substantial numerical results in both acoustics [31] and

electromagentics [32].

However, compared with the volume of literature on representation of transient fields

due to quiescent sources, research on transient scattering analysis using spherical harmonics

is somewhat limited. The methods investigated thus far fall into two categories; (i) using

inverse Fourier/Laplace transforms of Bessel and Hankel functions within a mode matching

framework, and (ii) using a novel representation of the retarded potential Green’s function.

The first time-domain Mie algorithms were presented for the scalar Helmholtz equation

in [17].

It is evident from the work that obtaining a stable algorithm is difficult as the

inverse Fourier transform of the spherical Hankel function is a Gagenbauer polynomial that

grows as a function of time. The inverse Fourier transform in the work has to be understood

in the convolution sense because the actual inverse Fourier transform doesn’t exist. A more

recent paper examined the causes of instability, errors in the extent literature, and presented

a recursive convolution approach that is stable for both late time and high order harmonics

for scalar fields [33]. Extension of classical Fourier methods to electromagnetic fields is

significantly more complicated [34], and given the instability for scalar fields, it is not clear

that this will be stable for vector fields. It is possible that the method developed in [33]

may be extended to the vector case. The second method for solving scattering from spheres

was developed before [33], but took an entirely different approach [3]. The underlying thesis

of this approach was to utilize the fact that if the trace of quantities of interest were to be

represented in terms of tesseral harmonics, then the solution of an integral equation composed

of these trace quantities would yield the Mie solution. This statement can be readily proven

in the frequency domain [35], and the proofs arise from a spherical harmonic representation of

the Green’s function. In time domain, this would translate to a representation of the retarded

potential Green’s functions, which would bring us back the original set of difficulties. In [3],

9

we proposed a novel representation of the Green’s function that obviates these difficulties,

and produces stable and convergent results.

The work presented herein would follow in the vein of [3] and leverage time domain

integral equation (TDIE) methods to develop a quasi-analytical solution to time dependent

EM scattering from a sphere. The approach presented is based on vector spherical harmonics

expansion for the time-dependent dyadic Green’s function (to be more specific, its tangential

trace) of both electric and magnetic field types, and integrated the results presented in [3].

The use of tesseral harmonics for representing the spatial variation of the trace quantities

on the surface of the sphere permits (i) a mesh free solution and (ii) analytical evaluation

of the integral required for the TDIE, leaving only a Volterra integral in time. The specific

contributions of the chapter are as follows: We

• present a stable spherical expansion for time domain dyadic Green’s functions of both

electric and magnetic types.

• present a spatially mesh-less and singularity-free integral equations for EM scattering

from spherical objects

• develop a set of reduced 1-D Volterra integral kernels that are solved with a higher

order Galerkin scheme

• present several results demonstrating the accuracy, convergence and stability of the

proposed methods.

The remainder of the chapter is organized as follows. Section 2.2 describes the problem to

be solved and gives the integral equation with time-domain dyadic Green’s function. Section

2.3 derives the vector spherical harmonic expansion for the tangential trace of the time-

dependent dyadic Green’s function. Section 2.4 gives a brief approach to solve the reduced

integral equation to obtain the time-dependent spherical multipoles. Section 2.6 presents

a number of results that validate the accuracy, convergence and stability of the presented

10

method. Finally, we summarize our contribution and outline future research directions in

Section 2.7.

2.2 Problem Statement

Consider a perfectly electrically conducting spherical object that occupies a volume D1 ⊂
R3 residing in a homogeneous background medium occupying D0 = R3\D1, characterized
by permittivity 
0 and permeability µ0. The boundary of the scatterer is denoted using
Ω1 = ∂D1, and is equipped with an outward pointing normal ˆn(r) and a radius of r = a.
On the interface Ω1, the tangential electric field vanishes.

The electric field integral equation (EFIE) and magnetic field integral equation (MFIE),

their linear combination, can be used to solve the electric current on the surface Ω1 (usually
with spatial meshing). The EFIE and MFIE are usually written as follows, respectively.

n × n × L(J(r, t)) = n × n × Ei(r, t)

(2.1)

J(r, t) + ˆn × K(J(r, t)) = ˆn × Hi(r, t)

(2.2)
where operators L(·) and K(·) are spatial-temporal (four-dimensional) operators associated
with dyadic Green’s function.

(cid:90)
S(cid:48) ˜Ge0(r, r(cid:48)) ⊗ X(r, t)dS(cid:48)
L(X(r, t)) = µ∂t
(cid:90)
S(cid:48) ˜Gm0(r, r(cid:48)) ⊗ X(r, t)dS(cid:48)

K(X(r, t)) = −

(2.3)

(2.4)

In the above definitions, ˜Ge0 and ˜Gm0, respectively, denote dyadic Green’s functions of
electric type and magnetic type in free space.The solutions to the above equations are typi-

cally effected using a discrete representation of the scatterer, representing the current using

appropriate spatial and temporal basis on this discrete geometric representation and finally,

creating a set of linear equations using testing procedures in both space and time. But if the

scatterer is spherical/canonical, one can take advantage of the geometry. This is done easily

in the frequency domain. In what follows, its transient analogue is presented. First, the

11

time-dependent Dyadic Green’s functions (to be more specific, their tangential traces) are

first expanded using vector spherical harmonics. And then the TD-EFIE and TD-MFIE are

reduced into a simple Volterra integral equations without any spatial meshing and singularity

in the integral kernels. It permits mode-by-mode solutions for coefficients of time-domain

Debye-Mie series.

2.3 Formulation

2.3.1 Natural Spatial basis: Vector spherical harmonics

On a sphere, vector spherical harmonics (VSH) form a natural basis set to represent the

spatial variation of the current [36, 37]. Their relations to spherical harmonics are given as

Ψm

n (ˆr) = Ψm

n (θ, φ) =

Φm

n (ˆr) = Φm

n (θ, φ) =

r(cid:112)n(n + 1)
1(cid:112)n(n + 1)

∇tY m

n (θ, φ)

¯r × ∇tY m

n (θ, φ)

where

Y m
n (θ, φ) =

2n + 1

4π

(n − m)!
(n + m)!

n (cos θ)ejmφ
P m

(2.5)

(2.6)

(2.7)

(cid:115)

(cid:88)

In the above equations, P m

n is the associated Legendre function of degree n and order m of
the mode. VSHs comprise a complete and orthogonal tangential basis set. Therefore, the

current on the sphere can be written as

J(r, t) =

J 1
nm(t)Ψm

n (θ, φ) + J 2

nm(t)Φm

n (θ, φ)

(2.8)

Given the coefficients of both modes, one can get the scattered electric and magnetic fields

n,m

using

Es(r, t) = −L(J(r(cid:48), t))
Hs(r, t) = −K(J(r(cid:48), t))

(2.9)

(2.10)

The basis here used is different from the basis used as in time-harmonic Debye-Mie series,

where the field, rather than the current, is represented using the vector spherical wave

12

functions (VSWF) with radial dependence [12] . However, we use VSWF next to derive

the expansions for the time-dependent dyadic Green’s functions in time-domain integral

equations.

2.3.2 Spherical expansion of Time-domain Dyadic Green’s function

The time-domain dyadic Green’s function of electric type is the solution to the following

wave equation

∇ × ∇ × ˜Ge0(r, r(cid:48), t) − ∂2
t
c2

˜Ge0(r, r(cid:48), t) = ˜Iδ(r − r(cid:48))δ(t)

and its dyadic form is notionally written as

˜Ge0(r, r(cid:48), t) = ( ˜I − c2∂−2

t ∇∇)

δ(t − r−r(cid:48)

c
4πR

)

(2.11)

(2.12)

where ∂−2
t
δ(t − R
c )
4πR

stands for a two-fold integral with respect to time, ˜I is the idempotent, and
is retarded potential. As shown later, ∂−2

is not evaluated explicitly. To glean

t

more insight into the analysis, we start with the frequency domain counterpart of the dyadic

Green’s function:

˜Ge0(r, r(cid:48), ω) =

(cid:104)˜I +

k2∇∇(cid:105)

1

G0(r, r(cid:48), ω)

(2.13)

where k and ω denote the wavenumber and the angular frequency, respectively. The rep-

resentation of the dyadic Green’s function in spherical coordinates is well known, has been

used extensively [35], and reads as

˜Ge0(r, r(cid:48), ω) = jk

(cid:88)

(cid:104)

n,m

N

(4)
nm(r, k)N

(1)∗
nm (r(cid:48), k)
(1)∗
nm (r(cid:48), k)

+ M

(4)
nm(r, k)M

(cid:105)

(2.14)

(p)
nm(k, r) and M

for |r| > |r(cid:48)|. In the above equations, N
(p)
nm(k, r) are the TMr and TEr vector
spherical wave functions (VSWF) of degree n and order m, respectively. In what follows,
we use r = |r| and r(cid:48) = |r(cid:48)|. The explicit expressions for N
found in the appendix in terms of spherical Bessel z

(p)
nm(k, r) can be
n (kr(cid:48)) or spherical Hankel functions

(m)
nm (k, r) and M

(1)

13

(4)
n (kr). As one can see from (2.14), the expansion is in terms of angular dependent wave
z
function (angular dependence is integrated with radial dependence at this moment), rather

than the vector spherical harmonics. Furthermore, the frequency dependence is sprinkled

throughout, and as a result taking an inverse Fourier transform is difficult. In order to fully

exploit the orthogonality of VSHs as well as make the expressions more amenable to taking

an inverse Fourier transform, it is necessary to recast the Green’s function in terms of its

tangential trace. This can be effected via the following: Taking the tangential trace of the

dyadic Green’s function via

e0(r, r(cid:48), ω) = ˆr × ˆr × ˜Ge0(r, r(cid:48), ω) × ˆr(cid:48) × ˆr(cid:48)
˜Gtt

(2.15)

yields

˜Gtt

e0(r, r(cid:48), ω) = jkˆr × ˆr ×(cid:88)
(cid:2)kr(cid:48)z
(cid:104)(cid:2)krz
(cid:88)

n (kr)(cid:3)(cid:48)

+M

n,m

(4)

(4)
nm(k, r)M

(1)∗
n
kr(cid:48)
(1)∗
(4)
n (kr)z
n

+z

N

(cid:104)
(kr(cid:48))(cid:3)(cid:48)

(1)∗
nm (k, r(cid:48))

(4)
nm(k, r)N

(1)∗
nm (k, r(cid:48))

(cid:105) × ˆr(cid:48) × ˆr(cid:48)
Ψnm(ˆr)Ψnm(ˆr(cid:48))
(cid:105)
(kr(cid:48))Φnm(ˆr)Φnm(ˆr(cid:48))

= jk

n,m

kr

(2.16)

(2.17a)

It can be shown that the trace of the Green’s function has the following properties when

operating on the tangential field X that is defined on a spherical surface:

e0(r, r(cid:48), ω) · X(r(cid:48)) = − ˜Ge0(r, r(cid:48), ω) · X(r(cid:48))
˜Gtt
e0(r, r(cid:48), ω) = −X(r) · ˜Ge0(r, r(cid:48), ω)
X(r) · ˜Gtt

As a result of these properties, one can use ˜Gtt

(2.17b)
e0(r, r(cid:48), ω) instead of ˜Ge0(r, r(cid:48), ω). This prop-
e0(r, r(cid:48), ω) is
expressed purely in terms of VSHs, and consequently one can exploit the orthogonality of

erty is valid in both frequency and time domain. As is evident from (2.16), ˜Gtt

these functions as well as take its inverse Fourier transform. The resulting expressions for

14

the dyadic Green’s function can be written as

(cid:88)

n,m

(4)

(cid:2)krz
n (kr)(cid:3)(cid:48)
e0(r, r(cid:48), t) =
(cid:104)F−1(cid:2)jk
˜Gtt
+ F−1(cid:2)jkz
(cid:104)
(cid:88)

kr

=

(1)∗
(4)
n (kr)z
n

nΨnm(ˆr)Ψnm(ˆr(cid:48)) + I2
I1
rr(cid:48)F−1(cid:2)jkz
(cid:105)
(cid:104)
rr(cid:48)I2

n in (2.18) can be rewritten as
c2∂−2
rr(cid:48) ∂r∂r(cid:48)
c2∂−2
rr(cid:48) ∂r∂r(cid:48)

(cid:104)

=

n

t

t

I1
n =

(1)∗
n
kr(cid:48)

(kr(cid:48))(cid:3)(cid:48)
(cid:2)kr(cid:48)z
(cid:3)Ψnm(ˆr)Ψnm(ˆr(cid:48))
(cid:105)
(kr(cid:48))(cid:3)Φnm(ˆr)Φnm(ˆr(cid:48))
(cid:105)
nΦnm(ˆr)Φnm(ˆr(cid:48))
(kr(cid:48))(cid:3)(cid:105)

(1)∗
(4)
n (kr)z
n

n,m
It can be shown that I1

(2.18)

(2.19)

As is evident from (2.18) and (2.19), one needs to obtain the explicit expression for only I2
n
in order to evaluate an expression of transient electric dyadic Green’s function. As is well

known, the approach of using the inverse Fourier transforms of spherical Bessel and Hankel

functions and then evaluating their convolution is unstable, because the inverse Fourier

transform of spherical Hankel function doesn’t exist.

In the next section, we present an

overview of a method that was recently developed to address this very problem [3], and then

use this to recover expressions for both the trace of the electric and magnetic dyadic Green’s

functions.

The starting point for deriving an expression that approximates I2 or provides an alter-
native representation arises from examining the representation for the Green’s function for

the Helmholtz equation, viz.,
exp [−jkR]

(cid:88)
F−1[−jkz

n,m

4πR

δ(t − R/c)

4πR

=

= −jk

(cid:88)

n,m

(4)
n (kr)z
z

(1)

n (kr(cid:48))Y m

n (ˆr)Y m
n

∗(ˆr(cid:48))

(4)
n (kr)z

∗(ˆr(cid:48))

(2.20)

n (ˆr)Y m
n

While this is the well known classical representation, an alternate representation may be

derived using

f (x) =

(cid:88)

n

anPn(x); an =

15

(1)

n (kr(cid:48))]Y m
(cid:90) 1

2n + 1

2

−1

dx f (x)Pn(x)

(2.21)

where Pn(x) is the Legendre function of nth order. Using x = cos γ = (r2 + r(cid:48)2 − R2)/(2rr(cid:48))
where γ denotes the angle between the points r and r(cid:48) or equivalently R =
r2 + r(cid:48)2 − 2rr(cid:48)x,
and f (x) = exp [−jkR]/(4πR), it is trivial to show that

√

The inverse Fourier transform to the above expression yields

n=0

∞(cid:88)
∞(cid:88)
(cid:18) ξ

n=0

(cid:90) r+r(cid:48)
|r−r(cid:48)| exp[−jkR]Pn (x) dR
(cid:18) ξ
2rr(cid:48)
n (ˆr)Y ∗m

Pαβ(t)Pn(x)
n (ˆr(cid:48))

(cid:19)

Pαβ(t)Y m

exp [−jkR]

4πR

=

2n + 1
8πrr(cid:48)

δ(t − R/c)

=

4πR

(cid:88)

n,m

=

c

2πrr(cid:48) Pn

2rr(cid:48)

(2n + 1)c
8πrr(cid:48) Pn

(cid:19)

where ξ = r2 + r(cid:48)2 − c2t2, Pαβ(t) is a pulse function (the value is 1 when t ∈ [α, β], and
zero otherwise) with α = |r − r(cid:48)|/c and β = (r + r(cid:48))/c. In obtaining these expressions, we
have used the addition theorem for Legendre polynomials. Comparing (2.23) to (2.20), it is

apparent that I2

n takes the form

F−1(cid:2) − jkz
(cid:16) ξ

c
2rr(cid:48) Pn

=

2rr(cid:48)

(cid:17)

(1)∗
(4)
n (kr)z
n

(kr(cid:48))(cid:3)

Pαβ(t) ≡ K

(0)

n (r, r(cid:48), t)

(2.24)

Using (2.24) together with (2.19), one can obtain the dyadic Green’s function of the electric

type as follows:

(2.22)

(2.23)

(2.25)

(2.26)

K

(0)

n (r, r(cid:48), t)Φnm(ˆr)Φnm(ˆr(cid:48))
(cid:105)
Ψnm(ˆr)Ψnm(ˆr(cid:48))

(cid:105)
n (r, r(cid:48), t)

(0)

˜Gtt

(cid:104)

e0(r, r(cid:48), t) = −(cid:88)
(cid:104)

n,m

rr(cid:48)K

+

t

c2∂−2
rr(cid:48) ∂r∂r(cid:48)
(cid:88)

Likewise, a similar procedure leads to the representation of the dyadic Green’s function of

the magnetic field.

(cid:104) ∂r(cid:48)

(cid:2)r(cid:48)K
n (r, r(cid:48), t)(cid:3)Φnm(ˆr)Φ∗

r(cid:48)

n,m

m0(r, r(cid:48), t) =
˜Gtt
(cid:2)rK

(0)

− ∂r
r

(cid:105)

nm(ˆr(cid:48))

n (r, r(cid:48), t)(cid:3)Ψnm(ˆr)Ψ∗

(0)

nm(ˆr(cid:48))

Whereas the above expressions are different from those that are used conventionally, their

scalar analogues have proven to be both robust (in terms of stability) and accurate.

16

2.3.3 Volterra Integral kernels and Reduced Time-domain Integral Equations

The expressions for the Green’s functions (or its tangential trace) depend only on VSHs, and

thanks to (2.16) they can be used instead of the dyadic Green’s functions in the expressions

for the EFIE and the MFIE. Specifically, it can be shown that by replacing the dyadic

Green’s function with its tangential trace, one can get the following integral operators;

(cid:90)
e0(r, r(cid:48)) ⊗ X(r, t)dS(cid:48)
Ltt(X(r, t)) = −µ∂t
(cid:90)
S(cid:48) ˜Gtt
m0(r, r(cid:48)) ⊗ X(r, t)dS(cid:48)
S(cid:48) ˜Gtt

Ktt(X(r, t)) = −

The above two operators are effectively equivalent to those in (2.1) and (2.2), but permit the

choice of vector spherical harmonics as basis functions for the unknown current densities.

Using the above equations together with Galerkin testing (that exploits the orthogonality

of the VSHs) results in a set of one-dimensional Volterra equations. For the TD-EFIE, we

obtain the following.

(2.27)

(2.28)

(2.29a)

(2.29b)

nm(ˆr),Ltt(Ψn(cid:48)m(cid:48)(ˆr(cid:48))) >
(cid:105)
n (r, r(cid:48), t)

δnn(cid:48)δmm(cid:48)

(0)

nm(ˆr),Ltt(Φn(cid:48)m(cid:48)(ˆr(cid:48))) >

< Ψ∗

< Φ∗

t

(1)

∂r∂r(cid:48)

nm(ˆr),L(Ψn(cid:48)m(cid:48)(ˆr(cid:48))) >=< Ψ∗
(cid:104)
= −µc2∂−1
rr(cid:48)K
rr(cid:48)
n (r, r(cid:48), t)δnn(cid:48)δmm(cid:48)
≡ K
(cid:104)
(cid:105)
nm(ˆr),L(Φn(cid:48)m(cid:48)(ˆr(cid:48))) >=< Φ∗
n (r, r(cid:48), t)
= −µ∂t
c
n (r, r(cid:48), t)δnn(cid:48)δmm(cid:48)
≡ K
(2)

(0)

K

δnn(cid:48)δmm(cid:48)

where the functions K

and Φnm, respectively. The off-diagonal elements(cid:10)Ψ∗

(1)
n and K

(2)
n are the one-dimensional Volterra integral kernels for Ψnm

nm,Ltt(Ψn(cid:48)m(cid:48))(cid:11) and(cid:10)Φ∗

nm,Ltt(Φn(cid:48)m(cid:48))(cid:11)

are identically zero. These equations can then be used to obtain the requisite Volterra inte-

gral equations for each of the two modes as

(1)

n (r, r(cid:48), t) ⊗ J 1

K

nm(ˆr), n × n × Ei(r, t) >

nm(t) =< Ψ∗
≡ f 1
nm(t)

(2.30a)

17

(2)

n (r, r(cid:48), t) ⊗ J 2

K

nm(ˆr), n × n × Ei(r, t) >

nm(t) =< Φ∗
≡ f 2
nm(t)

(2.30b)

In the above equations, the right-hand-sides are the projection of ˆn × n × Ei(r, t) onto the
two VSH basis sets.

A similar procedure can be followed to derive the necessary equations for the TD-MFIE.

The Volterra integral kernels K

(3)
n

and K

(4)
n

can be defined as follows.

< Ψ∗

< Φ∗

(3)

(0)

nm(ˆr),K(Ψn(cid:48)m(cid:48)(ˆr(cid:48))) >=< Ψ∗
n (r, r(cid:48), t)(cid:3)
(cid:2)r(cid:48)K
= −∂r(cid:48)
r(cid:48)
n (r, r(cid:48), t)δnn(cid:48)δmm(cid:48)
≡ K
n (r, r(cid:48), t)(cid:3)
(cid:2)rK
∂r
=
r
n (r, r(cid:48), t)δnn(cid:48)δmm(cid:48)
≡ K

nm(ˆr),K(Φn(cid:48)m(cid:48)(ˆr(cid:48))) >=< Φ∗

(0)

(4)

nm(ˆr),Ktt(Ψn(cid:48)m(cid:48)(ˆr(cid:48))) >

nm(ˆr),Ktt(Φn(cid:48)m(cid:48)(ˆr(cid:48))) >

Using these kernels, the resulting Volterra equations for the TD-MFIE are

J 1
nm(t) + K

(3)

n (r, r(cid:48), t) ⊗ J 1
=< Ψ∗

nm(t)

nm(ˆr), n × Hi(r, t) >≡ f 3

nm(t)

J 2
nm(t) + K

(4)

n (r, r(cid:48), t) ⊗ J 2
=< Φ∗

nm(t)

nm(ˆr), n × Hi(r, t) >≡ f 4

nm(t)

(2.31a)

(2.31b)

(2.32a)

(2.32b)

Equations (2.30a), (2.30b), (2.32a) and (2.32b) constitute four equations that yield in-

dependent coefficients for the current on the surface of the sphere. In some sense, they can

be considered a VSH transform of the integral operator. If necessary, one can combine the

requisite integrals to form an effective time domain combined field solution as well. That

said, it should be noted that the solution to scattering from a sphere has been reduced to

solving these two equations that are uncoupled. In what follows, we shall provide a method

for solving these equations numerically.

18

2.4 Time-dependent Debye-Mie series solution

Analytic solution to the above Volterra equations do not exist. As a result, numerical
nm(t). Note that,
in this section, subscript (nm) for currents, coefficient vectors, right-hand-side vectors and

solution to these equations rely on discrete representation of J 1

nm(t) and J 2

system matrix representations is suppressed. The approach that we follow is to partition the

time axis uniformly, and represent basis and test function whose support reside on each time

step. The basis functions can be any set of polynomials [38,39], however, we choose Legendre

polynomials with support over the time step. This is in contrast to traditional methods [32]

wherein basis functions were based on backward looking Lagrange polynomials that span

multiple time steps. However, as was shown in [38], using higher order polynomials within

a time step results in higher order convergence. Specifically, our representation

J ξ(ti + τ ∆t) =

ijΠiPj(2τ − 1);
J ξ

τ ∈ [0, 1)

(2.33)

j=0

where Πi = 1 for t ∈ [ti, ti+1), Pj(·) is a Legendre polynomial, ξ = 1, 2, and i = 1··· Nt
where Nt is the total number of time steps. Using Galerkin testing results in the marching-
on-in-time system. In stating the requisite equations, we use a block matrix representation

that is indexed by time step as opposed to unknowns, viz.,

Np(cid:88)

q−1(cid:88)

j=j0

Zν
0Iξ

q = Vν

q −

q−jIξ
Zν
j ;

ν = 1, 2, 3, or 4

(2.34)

(cid:104)

(cid:105)T

where j0 = max(0, q − Nk), q is the time step index ranging from 0 to total number of steps
Nt, and Nk is the length of the discrete kernel. Here, Iξ
are current
coefficients for those polynomials associated with time step j. In a similar fashion, Zν
j and
Vν
j (the one-index subscript denotes the discretization in time) are ,respectively, discrete
system and right-hand-side term at jth time step. At any given time t, the right-hand-side

j1,··· , J ξ
J ξ

j =

jNp

vector can be written as

Vν
q =

(cid:104)(cid:10)ΠqP1(·), f ν(·)(cid:11) ,··· ,

(cid:68)

(cid:69)(cid:105)T
ΠqPNp(·), f ν(·)

(2.35)

19

and the element of the system matrix can be written as

Zν
pi,qj =< ΠpPi(·), K

(ν)

n (·) ⊗ ΠqPj(·) >

(2.36)

In the expressions for vector and matrix elements, the two-index and four-index subscripts

denote the discretization in space.

In the above equations, the temporal dependencies of the function follows from their

explicit definitions earlier in the Section. It should be noted that K

(1)
n has an infinite tail,
very similar to that encountered while evaluating the contribution of the scalar potential in

regular TD-EFIE [38]. As in [32, 40], introduction of an auxiliary charge ameliorates the
computational complexity from O(N 2
t ) to O(Nt). Introducing an auxiliary charge the MOT
system for K

(1)
n

can be written as
ZI
0I1

q = V1

q −

ZI
q−jI1

j − Cj−1

(2.37)

q−1(cid:88)

j=j0

j(cid:88)

where j0 = max(0, q − Nk), Nk is the length of the non-constant part of the discrete kernel
(1)
n . The auxiliary charge coefficient is chosen as
K

Cj = Cj−1 +

ZI
j−kIk

(2.38)

where k0 = min(0, j − Nk). The elements of the modified system (ZI) are defined as follows.

k=k0

ZI
k = Z1

k − Z1

k−1

(2.39)

Equations (2.34) and (2.37) can be solved for the coefficients in every time step, and of

course for each harmonic. On another brief note regarding discretization of these equation;

since both testing and source basis functions are pulse functions that are multiplied by a

higher order polynomial, it can be shown that the inner product in (2.36) can be equivalently

(ν)
n with an interpolatory polynomial that is higher
written as a convolution of the kernel K
order with support from (−∆t, ∆t) and tested by a delta function. The resulting matrix
elements are smooth, well behaved, and causal. More importantly, they do not require

anterpolation unlike that used in [3, 41].

In what follows, we briefly describe the eigen-

spectrum analysis [40, 42] of the kernels presented for sample harmonics.

20

2.5 Eigen-spectrum analysis of the MOT system

Analysis of instability of direct Fourier inversion can be attributed to rapid growth of co-

efficients as elucidated in [33], which also provides a method to overcome this rapid growth.

By contrast, this approach results in Volterra integral equations whose kernels, when con-

volved and tested with temporal basis function, are sufficiently smooth (other than at isolated

points that are handled analytically). As alluded to in a series of papers by Brunner (see [39],

and refs therein) on using DG like methods for Volterra equations, the smoothness of the

resulting kernels assures stability. Next, we examine the stability of the MOT system devel-

oped earlier by analyzing the eigen-spectrum of the marching system for different harmonics.

We note that this section is provided purely for completeness; the equations provided below

are similar to those presented in [40,42] by comparing the similarity of the kernels analyzed

herein to those presented in these papers. Specifically, for equations involving K

MOT can be rewritten in the following matrix system form

 A11 0

I

0

 Ij+1 = F j+1 −

 B11 B12

B21 B21

 Ij,

(1)
n , the

(2.40)

j , 0,··· , 0, 0,··· , 0]T (here the
where Ij = [IT
superscript T denotes transpose) and I is the identity matrix. The elements in other matrices

j−1,··· , CT

j ,··· , J T

]T , F j = [IT

j−1−Nk

j−Nk

, CT

21

can be computed as follows.





A11 =

B11 =

B12 =

B21 =

B22 =



(cid:105)

(cid:104)ZI

0

0

0

0
[I]

0

. . .

0
. . .
[I] 0 . . .

.

(cid:104)ZI

. . .

(cid:105) (cid:104)ZI

(cid:105)

. . .

2

[I]

(cid:104)ZI

0
[I]

[I] 0

0

0

.

[I]

. . .

. . .

1
[I]

0

.

.

. . .





(cid:105)

(cid:105)

Nk
. . .

. . .

0

 ,

Nk
. . .

,





,

,

(2.41a)

(2.41b)

(2.41c)

(2.41d)

(2.41e)

.

(cid:104)ZI

1

0

0

(cid:105) (cid:104)ZI

(cid:105)

. . .

2

. . .

(cid:104)ZI

0

0

.

. . .

.



0

 ,

0

. . .

. . .

[I]

0

[I] 0
[I] 0

.

.

. . .

22

For systems involving other kernels, the corresponding set of equations are

A Ij = F j − B Ij−1

(2.42)

j , 0,··· , 0]T . Compared to the updating equations
where Ij = [IT
for the first kernel, only current vectors are considered, and the elements of the matrices A

j ,··· , J T

j−Nk

]T , F j = [IT

and B can be similarly obtained by changing the corresponding kernels in A11 and B11.

To obtain stability of the marching system, eigen-spectrum analysis should be done for
matrix A−1B. And in next section, the results of this analysis are presented for different
kernels and harmonics.

2.6 Numerical Examples

The contribution of this work lies in deriving quasi-analytic transient Debye-Mie series

via time-domain integral equation and its marching on in time solution. Hence, we provide

several examples of using TD-EFIE and TD-MFIE to test the accuracy of the method and to

demonstrate the stability of the MOT algorithm. However, to understand the computational

complexity of the method, assume that the maximum wavenumber to be modeled is κmax

and the radius of the sphere is a. Then the number of degrees of freedom for classical
tesselation scales as Ns = O(κ2
maxa2), and the order of harmonics required Nh = O(κmaxa).
It follows that the cost of traditional solution scales as O(NtN 2
s ) and that using the method
presented here is O(NtN 2
h), where Nt is the number of time steps used in the simulation.
In what follows, comparison for a single mode is presented as opposed to their sum. The

problem studied here is scattering by a perfectly electrically conducting (PEC) sphere of

radius a = 1m. The incident field is a plane wave with a modulated Gaussian profile

Ei(r, t) = ˆx cos(2πf0t)e−(t−¯r·ˆk/c−tp)2/2σ2

(2.43)

where ˆk = ˆz is the direction of propagation, f0 is center frequency.
σ = 3/(2πB) and tp = 40σ, where B denotes the bandwidth. In all the numerical examples
of this section, center frequency of f0 = 0.4GHz and bandwidth of B = 0.3GHz are chosen.

In the expression ,

23

Figure 2.1: Time-domain coefficients for mode Ψ1
3

With this configuration, the frequency band we are interested in is between 0.1GHz and

0.7GHz. Due to the fact that the power spectral density value at fmax = 0.7GHz is 39dB
lower than that at the center frequency, the number of spherical harmonics to (spatially)
represent the plane wave is bounded by Nm ≈ [2ka] = 30. In this work, unless specified
and the width of the support of temporal
otherwise, the time step size is chosen as ∆t = 1
basis function is ∆t.

20fm

The first test is to illustrate the accuracy of the proposed method by comparing its

frequency-domain counterpart. In this test, the temporal basis function are chosen up to
3 are
plotted in Fig. 2.1 and Fig. 2.2. Each plot contains two curves that are obtained using both

first order Legendre polynomials. The time-dependent coefficients for modes Ψ1

3 and Φ1

TD-EFIE and TD-MFIE with corresponding kernels. The comparison in frequency domain

within the frequency band is demonstrated in Fig.2.3, and the results from frequency domain

24

-0.0015-0.001-0.0005 0 0.0005 0.001 0.0015 0.002 70 80 90 100 110 120 130time-dependent cofficientstime/nsTD-EFIE with kernel (1)TD-MFIE with kenrel (3)Figure 2.2: Time-domain coefficients for mode Φ1
3

Mie (FD-Mie) series are also provided to verify the frequency responses. One can observe

that the agreement is excellent between time and frequency domain. Similarly conclusion

can be drawn for higher order modes. Time-domain and frequency-domain results for modes
30 are given in Figs. 2.4-2.5 and Fig. 2.6, respectively. One can also observe
that coefficients for modes of degree 30 are orders of magnitude smaller that those of degree

30 and Φ1

Φ1

3. From these results, one can recognize the feasibility of extracting transient response of

spherical object with the presented mode by mode MOT approach. In the simulations we
conducted, relative error convergence in frequency domain (down to 10−12) versus the order
of the temporal basis function (up to 6) has been observed.

The second test is to demonstrate the stability of the marching algorithm. Simulations

with two different time steps are carried out for the four reduced Volterra kernels, where the

temporal basis functions are chosen up to first order. The solutions are transient coefficients

25

-0.0015-0.001-0.0005 0 0.0005 0.001 0.0015 0.002 70 80 90 100 110 120 130time-dependent cofficientstime/nsTD-EFIE with kernel (2)TD-MFIE with kenrel (4)Figure 2.3: Frequency-domain coefficients for modes Ψ1
3,compared to two curves
(denoted by FD-Mie-psi and FD-Mie-phi respectively) with frequency domain Mie seires
approach.

3 and Φ1

for the modes Ψ1

30 (solutions to VIE of kernels 1 and 3) and Φ1

30 (solutions to VIE of kernels
are chosen (oversampling of 5 and 25 times

2 and 4). Time steps ∆t = 1
the Nyquist’s sampling rate), and results with 50,000 and 250,000 steps, respectively, are

and ∆t = 1

100fm

20fm

plotted in Fig. 2.7. It’s obvious that no late-time instability is observed in these simulations.

To further show the stability of the resulting MOT scheme, the eigenvalues analysis of the

discretized marching system outlined earlier is given. The real and imaginary parts of the

eigenvalues are plotted in Figs. 2.8a-2.8d, with all the values are within the unit circle

except one lying on the unit circle which is associated with kernel 1. Similar stability and

eigen-analysis results can be also found for simulations with smaller time steps, higher order

temporal basis functions and Volterra kernels of different degrees.

26

 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0.005 0.1 0.2 0.3 0.4 0.5 0.6 0.7frequency-dependent cofficientsFrequency/GHzFD-Mie-psiTD-EFIE with kernel (1)TD-MFIE with kernel (3)FD-Mie-phiTD-EFIE with kernel (2)TD-MFIE with kernel (4)Figure 2.4: Time-domain coefficients for mode Ψ1
30

2.7 Conclusion

In this work, a TDIE-based direct approach is proposed to calculate the time-dependent

spherical multipoles due to the scattering by a PEC sphere. By using spherical harmonics

as basis functions and expanding the tangential trace of time-domain dyadic Green’s func-

tion with VSH, the proposed method analytically reduces the original time-domain integral

into several novel one-dimensional Volterra integral equations. Those reductions allow very

efficient computation for those time-domain spherical multipoles in a marching manner. As

shown in the numerical examples, the accuracy and convergence of discretization in time are

verified. The resulting kernels could also be considered as good testbeds for various tempo-

ral basis function choices and testing schemes. Stability issue arising in conventional TDIE

could be studied with those reduced integral equation kernels, which is a current research

27

-8e-12-6e-12-4e-12-2e-12 0 2e-12 4e-12 6e-12 8e-12 70 80 90 100 110 120 130time-dependent cofficientstime/nsTD-EFIE with kernel (1)TD-MFIE with kenrel (3)Figure 2.5: Time-domain coefficients for mode Φ1
30

topic.

While the solution presented solves a long standing intellectual problem, it may seem

sparse from an application perspective. However, interesting applications exist. For interest,

there is a significant interest in developing methods for scattering from clusters of spherical

particles that could be accelerated by fast methods [32]. Furthermore, it has many practical

applications that its frequency domain counterpart cannot easily handle. First, the method

can be used as a local but non-reflecting boundary conditions for time-dependent differential

equations based solvers. Second, the proposed time-dependent approach, coupled with other

physics-based solvers, can help to model interactions between waves and nonlinear particles.

These ideas are currently active topics of research in our group.

This chapter, c(cid:13)2015 IEEE, is reprinted, with permission, from Li, J. and Shanker,B.,
"Time-Dependent Debye-Mie Series Solutions for Electromagnetic Scattering," Antennas

28

-1.5e-11-1e-11-5e-12 0 5e-12 1e-11 1.5e-11 70 80 90 100 110 120 130time-dependent cofficientstime/nsTD-EFIE with kernel (1)TD-MFIE with kenrel (3)Figure 2.6: Frequency-domain coefficients for mode Ψ1
30„compared to two
curves (denoted by FD-Mie-psi and FD-Mie-phi respectively) with frequency domain Mie
seires approach.

30 and mode Φ1

and Propagation, IEEE Transactions on, August 2015.

29

 0 2e-10 4e-10 6e-10 8e-10 1e-09 1.2e-09 1.4e-09 0.1 0.2 0.3 0.4 0.5 0.6 0.7frequency-dependent cofficientsFrequency/GHzFD-Mie-psiTD-EFIE with kernel (1)TD-MFIE with kernel (3)FD-Mie-phiTD-EFIE with kernel (2)TD-MFIE with kernel (4)Figure 2.7: Late-time stability

30

10−2010−10100Coefficientskernel 1 oversampled 5 timeskernel 1 oversampled 20 times10−2010−10100CoefficientsKernel 2 oversampled 5 timesKernel 2 oversampled 20 times10−2010−10100CoefficientsKernel 3 oversampled 5 timesKernel 3 oversampled 20 times10−2010−10100       0     500    1000    1500    2000    2500    3000    3500Coefficientstime/nsKernel 4 oversampled 5 timesKernel 4 oversampled 20 times(a)

(c)

(b)

(d)

Figure 2.8: Eigenvalues of MOT system of (a) kernel 1, (b) kernel 2, (c) kernel 3 and (d)
kernel 4.

31

−1−0.500.51−1−0.500.51real partimaginary part−1−0.500.51−1−0.500.51real partimaginary part−1−0.500.51−1−0.500.51real partimaginary part−1−0.500.51−1−0.500.51real partimaginary partCHAPTER 3

NUMERICAL SOLUTIONS ON HIGHER ORDER GEOMETRY

REPRESENTATIONS

3.1

Introduction

Surface integral equation has been a very successful tool to analyze the electromagnetic

(EM) scattering from homogeneous scatterers. Among them, electric field integral equation

(EFIE) and its related forms are widely used. To solve the integral equation on the surface,

a geometry representation and a physical quantity (such as the commonly used surface

current density) representation are needed. The choice of physical function representation is

usually strongly related to the choice of geometric representation. Choosing a good geometry

representation is necessary to reduce the numerical error and also helps the construction of

basis functions. Typical geometry representation directly employed in method of moments

(MoM) is surface mesh, suffering from geometry error because the mesh is reduced from the

highly accurate geometry natively defined in computer-aided design (CAD) tools. There are

mesh-free methods and more recently isogeometric analysis (IGA) based methods available to

at least partly solve the problem [5, 43–46]. Furthermore, avoiding geometry representation

reduction explicitly can also improve the efficiency of the design-through-analysis process.

Isogeometric analysis starts with the original geometry representation format, such as non-

uniform rational B-spline (NURBS) surface, subdivision and their modifications, and embeds

the native geometry information in the basis design.

Two aspects, at least, of IGA should be developed further for future industry applications;

(1) generalization to objects of arbitrary shape and topology and (2) convergence properties

of the numerical scheme. Regarding the former, div-conforming basis function defined on

curved surfaces, such as curvilinear RWG basis functions should be able to capture all the

components in the function space of the current density. [47] shows that optimal convergence

32

of solutions can be obtained for high order basis functions if geometry representation satisfies

certain requirement on the curvature of the cell. In this chapter, the subdivision surface,

a powerful geometry representation technique that features high continuity and arbitrary

topology, is explored to construct boundary integral equation solvers using the method of

moments. Subdivision surface, mostly used in the animation industry, is a smooth geometry

representation that can be considered as the limit surface after infinite refinements of primal

control mesh. Loop subdivision [48] is one of the commonly used subdivision surface, which

will be used in this work. The limit surface of Loop subdivision features C2 continuity almost

everywhere except finite number of isolated irregular vertices. Explicit refinement is not

required in practical implementation thanks to the technique proposed by Stam [49]. With

that approach, one can evaluate the limit position and up to second order surface differential

operators at almost arbitrary points in the parameter domain. Generally, advantages of

subdivision surface over other parametric surface like NURBS surface include (1) that it

supports global continuity up to C2 and (2) that it supports control mesh of arbitrary

topology and genus (watertight and no gap issue).

This chapter covers the following,

• Isogeometric analysis on smoothly represented structures of arbitrary topology,

• Constructing Div-conforming basis of various orders on subdivision surface,

• Study of super-convergence properties of IGA-MoM.

3.2 Subdivision Surfaces

In what follows, we provide a brief overview of shape description as effected by subdivision

surfaces. This section is provided purely for completeness and omits details that can be

found in several references, e.g., [48–50]. The traditional workhorse of geometric description

is NURBS. As a result, current iso-geometric methods have been largely developed using

NURBS as basis functions. However, a singular feature stands out that poses to be a

33

bottleneck–while NURBS descriptions are C2 in the interior of a patch, the continuity may

be C0 or even worse (not watertight) at the boundary between patches or at intersection

curves. As a result, one has to define physical basis functions that are div- or curl-conforming

for vector problems. Specifically, to solve for currents through electric field integral equation,

the basis function has to be div-conforming [51–53]. This stands in contrast to subdivision

surfaces that are constructed as limit surfaces obtained by subdivision processes. This leads

to meshes that are C2 almost everywhere except at isolated points where they are C1.

As triangular tessellation is ubiquitous in computational electromagnetics, the choice

of the presentation below is based on the Loop subdivision scheme for closed triangular

surface mesh. Let M0 denote the primal base mesh or control mesh that leads to the limit
surface. As in any tessellation, M0 comprises a set of vertices V and a connectivity map.
The 1-ring (union of incident triangles) of a vertex produced by the connectivity map can be

characterized by the valence of the vertex; specifically, (i) a valence-6 node/vertex is deemed

regular and (ii) any other vertex is called extraordinary or irregular. A triangle is regular if

all its vertices are regular, and irregular otherwise. Fig. 3.1 illustrates an irregular triangle

and one subdivision around it.

(a)

(b)

Figure 3.1: Irregular vertex and triangle (a) An irregular triangle (vertex 1 is valence-7),
(b)subdivision once

In what follows, we present a high level prescription of Loop subdivision. Loop sub-

34

Figure 3.2: Local averaging during refinement.

division prescribes a smooth surface by defining an infinite sequence of meshes that are

determined through repeatedly refining the k-th mesh into the (k + 1)-th mesh by divid-

ing each triangle into four subtriangles. Each subdivision step involves inserting one vertex

per edge and adjusting vertex locations. The location of the newly inserted vertex is an

affine combination of the end vertices of the edge, each with weight 3/8, and the other

two vertices of the two incident triangles, each with weight 1/8. The location of the exist-

ing vertex on the refined mesh is an affine combination of its one-ring vertices with weight
αn = (5/8 − (3/8 + cos(2π/n)/4)2)/n, and itself with weight 1 − nαn, where n is the va-
lence. The limit of the sequence of piecewise flat surface is the Loop subdivision surface

with C2 smoothness almost everywhere, except at the extraordinary vertices, where it has

C1 smoothness. Fig. 3.3 gives an example of subdivision surfaces.

In practice, the above prescription is not followed, i.e., there is not infinite refinement of

meshes. On can alternatively view the surface in terms of basis functions defined at vertices.

Consider the a limit surface S(u, v) =(cid:80) xiξi(u, v) that is defined by coordinates xi of the

vertices Vi; ξi(u, v) is an effective basis function that is associated with quantities associated
with vertex Vi, and (u, v) are the pairwise coordinates on a parameterization chart.

In the case of regular vertex, functions ξi(u, v) are box-splines. In the case of irregular
vertex, ξi(u, v) can be constructed by the subdivision procedure starting with 1 assigned
to Vi and 0 assigned to all other vertices. The process appears to be, again, not efficient
if the evaluation is required at arbitrary location. But one can subdivde k times, instead

35

888nn1nnnnnn8/18/18/38/3(a)

(c)

(b)

(d)

Figure 3.3: Subdivision surface (a) primal mesh , (b) one subdivision, (c) two subdivisions
and (d) three subdivisions.

of unlimited times, when the interested location (location referenced by parametric values

(u, v)) lies in a sub-triangle (Fig.3.4) with three regular vertices. Then the basis function can

vertices of the interested sub-triangle. The relevant values for each subtriangle Ωk

be evaluated through linear combinations of box splines with values on the relevant one-ring
m can be
stored in a vector Ek,m containing 12 scalars, one for each vertex whose box spline basis
function has the triangle in its support. Ek,m can be calculated as

Ek,m = Pm ¯AAkE0,

where E0 is the initial assignment of the n + 6 vertices in Fig. 3.1a, which is 1 for Vl
when evaluating ξl, and 0 for any other vertex; A is the subdivision operator matrix, an
(n + 6)× (n + 6)-matrix denoting how the values influencing the irregular triangle in the next
level are calculated from the values at the current level (Fig. 3.1b shows the vertex at the next

level); ¯A is an extended version of A including how the vertices n + 7 through n + 12 of next

36

level is computed form the current level; Pm is the picking matrix selecting the relevant 12
m within the n+12 values. Note that, in practical implementation, eigenanalysis
values for Ωk
of A is performed so the evaluation involves taking powers of the eigenvalues as an efficient

way of handling the calculation of Ak. Furthermore, the subdominant eigenvectors provide

a means for direct evaluation of derivatives even in irregular triangles [49].

Ω1
3

Ω2
3

Ω3
3
Ω3
2

Ω2
2
Ω3
1

Ω1
2

Ω2
1

Ω1
1

Figure 3.4: The parameter domain within an irregular triangle is recursively partitioned into
an infinite number of subdomains.

3.3 Formulation and Discretization

This section will formulate the discrete scheme to solve the integral equation for EM

scattering from PEC which is represented by the subdivision surface. Without loss of gener-

ality, only EFIE is solved, and extension to uniquely solvable IEs is straightforward. In this

chapter, time harmonic factor ejωt is implicitly assumed and suppressed.

3.3.1

Integral Equation for PEC

From the equivalence theorem, surface current density J(r) can be used as the unknown of

the EFIE. On the interface, from the boundary condition, one can build the EFIE as follows:

Ei(r) + Es(r)

= 0

(3.1)

ˆn × ˆn ×(cid:16)

(cid:17)

37

The scattered field can be represented as

(cid:90)

∂Ω

Es(r) = −jωµ0

dr(cid:48)(cid:20)

(cid:21)

I +

1

k2∇∇

g(r, r(cid:48)) · J(r(cid:48))

(3.2)

The scalar Green’s function in free space g(r, r(cid:48)) = e−jk|r−r(cid:48)|
in free space.

4π|r−r(cid:48)| where k is the wave number

3.3.2 Div-conforming Basis on Smooth Subdivision Surface

In [5], C2 continuous subdivision scalar basis was used to represent the scalar potentials and

the current is expanded in the form of the surface gradient and surface vector curl of the

potentials. Therefore no global loop component is represented when that scheme is applied

to multiply connected structures. To capture that component, one solution is to use edge

oriented basis as in RWG basis [52]. Using the Piola transformation, one can construct high

order div-conforming basis on the subdivision surface:

fi(r) =

DF
Js

¯f (r)

(3.3)

where ¯f (r) is the div-conforming basis of any order in the reference domain, DF is the

jacobian matrix of the transformation from reference domain to subdivision surface, and J

is the jacobian or the determinant of the jacobian matrix. In (3.3), the basis function is

embedded with jacobian information of the subdivision surface.

After the basis function is chosen, Galerkin testing is used to set up the linear system

which can be solved either directly or through an iterative solver.

3.4 Numerical Examples

In this section, numerical examples are given to demonstrate the convergence and effec-

tiveness of the presented numerical scheme.

The first example is the analysis of scattering from a unit sphere illuminated by the an

incidence plane at 100 MHz, where first order basis is used. The convergence in far field

38

solutions versus the number of degrees of freedom (DoFs) is plotted in Fig. 3.5. The order of

the convergence versus the number of DoFs is 1.5, or 3 versus the average size of the control

element. This super-convergence behavior in the far-field zone is expected as in [47].

Figure 3.5: Convergence of the relative error in far field for the sphere

The second example is on the scattering from a genus-1 (with one hole) torus. The inner

and outer diameters are 1.47λ and 2.52λ respectively, and the incidence is along the axis of

the torus. The observed convergence rate shown in Fig.3.6 for first order basis is of the same

order as expected.

3.5 Conclusion

This work applies the IGA to multiply-connected structure by using div-conforming ba-

sis constructed with the help of jacobians of subdivision surfaces. Numerical examples are

presented for validation of the efficacy and excellent convergence performance as well. De-

tails of the implementation, such as function evaluation and singular integral calculation on

high order smooth surfaces, together with more numerical results will be presented at the

conference.

39

103104105number of DoFs10-610-510-410-310-2relative error in far fieldDiv-conforming basisreference curve N-1.5Figure 3.6: Convergence of the relative error in far field for the torus

40

103104105number of DoFs10-510-410-3relative error in far fieldDiv-conforming basisreference curve N-1.5ISO-GEOMETRIC ANALYSIS ENHANCED METHOD OF MOMENTS

CHAPTER 4

4.1

Introduction

Computational methods have become the mainstay of scientific investigation in numer-

ous disciplines, and electromagnetics is no exception. Research in both integral equation

and differential equation based methods has grown by leaps and bounds over the past few

decades. This period has witnessed development of both higher order basis functions [53–57],

and higher order representations of geometry (at least locally) [56–58] amongst many other

equally important advances. These approaches have been applied to a wide range of realistic

problems spanning several wavelengths.

However, despite advances in these areas, there is a fundamental disconnect between the

geometry processing and analysis based on this geometry. Traditional analysis proceeds by

defining a discrete representation of the geometry typically comprising piecewise continu-

ous tessellations. Ironically, this discrete representation of the geometry is obtained using

software or a computer aided design (CAD) tool that contains a higher order differentiable

representation of the geometry. As eloquently elucidated in [44], the rationale for this dis-

connect can be attributed to the different periods in time that CAD tools and analysis tools

developed. As the latter is older, the computational foundation is older as well. As a result,

one is left with awkward communication with the CAD software for refining and remeshing.

This is especially true insofar as accuracy is concerned; lack of higher order continuity in

geometry can cause artifacts if the underlying spaces for field representations are not prop-

erly defined. Indeed, the need to define div/curl conforming spaces on tessellations that are

only C0 led to development of novel basis sets that meet this criterion [59]. An alternate
approach that has recently been espoused is iso-geometric analysis (IGA). In this approach,

the basis functions used to represent the geometry are the same as those used to represent

41

the physics on this geometry. As a result, the features of geometry representation such as

higher order continuity, adaptivity, etc., carry over to function representation as well. Since

the appearance of this approach in 2005 [44], it has been applied to a number of applica-

tions that range from structure mechanics [60] to fluid-structure interactions (FSI) [61] to

contact problems [62] to flow [63] to shell analysis [64, 65] to acoustics [66] and electromag-

netics [67].

In addition to analysis techniques, the power of IGA has been harnessed for

design-through-analysis phase in several practical applications [68–70].

Next, we briefly review some of the existing methods. Most CAD tools use bi/tri-variate

spline based patches/solids like those based upon Bezier, B-splines, and non-uniform rational

B-splines (NURBS). As a result, these basis functions are the most often used as IGA basis,

with the most popular being NURBS. The latter choice is determined by the fact that

NURBS is the industry standard for modern CAD systems. Properties such as non-negativity

and the fact that it provides a partition of unity make it an excellent candidate for defining

function spaces. Finite element methods based on NURBS basis functions that exhibit h- and

p-adaptivity have been demonstrated [44]. Unfortunately, the challenge with using NURBS

arises from the fact that the resulting shapes are topologically either a disk, a tube or a torus.

As a result, stitching together these patches can result in surfaces that are not watertight.

These complexities are exacerbated when the object being meshed is topologically complex or

has multiple scales [71,72]. Two other geometry processing methodologies gaining popularity

for handling shapes that are complex are T-splines and subdivision surface. The former, an

extension to NURBS, can handle T-junctions and hence and greatly reduce the number

of the control points in the control mesh. T-splines, especially analysis-ready T-splines,

comprise a good candidate for constructing iso-geometric analysis. More detailed work on

T-splines and its application in IGA can be found in [72, 73] and references therein. As

opposed to T-splines, subdivision surfaces have played a significant role in the computer

animation industry. Among its many advantages are the ease with which one can represent

complex topologies, scalability, inherently multiresolution features, efficiency and ease of

42

implementation. Furthermore, it converges to a smooth limit surface that is C2 almost

everywhere except at isolated points where it is C1 continuous. There are several subdivision

schemes; Loop, Catmull-Clark, Doo-Sabin to name a few. Generally speaking, all of the

three of these schemes alluded above can be used to construct IGA method. To date, iso-

geometric analysis based on subdivision surface is less well studied. Some work on IGA

based on Catmull-Clark can be found in [74], where IGA is used to solve PDEs defined on

a surface.

While the literature on IGA for differential equations is reasonably widespread across

multiple fields, IGA for integral equations is still at a nascent stage. As a result, it has

recently become the focus of significant attention. Recently, two dimensional iso-geometric

boundary element method (IGBEM) was proposed [75] to study elastostatic problem with

NURBS interpolating basis to represent the geometry, displacement and tractions. In [76],

IGBEM based on unstructured T-splines was developed for a three dimensional linear elas-

tostatic problem. This approach was extended to address IEs associated with hydrodynamic

interactions [68]. Likewise, IGBEM methods have been developed to study acoustic scat-

tering from rigid bodies [66] as well as two-dimensional electromagnetic analysis [77]. To

our knowledge, IGA has not been used for solution to three dimensional IEs associated with

vector electromagnetic fields, and this serves to motivate this chapter.

The focus of this chapter will be on the construction of a well formulated low-frequency

stable IGA solver for the electric field integral equation (EFIE) that is based on subdivision

surfaces specifically, the Loop subdivision scheme. As will be evident, the choice of Loop

subdivision scheme is only incidental; the presented method can be applied to most subdi-

vision surface description. In developing a solver that is robust, several challenges need to

be addressed; these range from definition of basis functions that correctly map the trace of

fields on the surface to formulations that render the resulting system frequency stable to

formulation of effective preconditioners. To set the stage for introduction of this formula-

tion, we will assume that the surfaces are simply connected and have C2 smoothness almost

43

everywhere. As will be evident, the assumption of sufficient “local” smoothness permits sig-

nificant freedom in terms of defining function spaces. Thus, the principal contribution of

this work is four-fold: We will

• present construction of a basis to correctly represent the tangent trace fields on simply-

connected surface,

• demonstrate convergence of EFIE-IGA solver for canonical geometries as well as present

applications for complex targets,

• demonstrate stabilization of EFIE-IGA solvers using the proposed basis sets together

with frequency scaling,

• demonstrate construction of the Calderón preconditioner using the proposed basis sets

(without the need for auxiliary barycentric meshes),

• and, present scattering data from multiscale obstacles.

The chapter is organized as follows: In Section 4.2 outlines the problem, whereas Section

3.2 presents a brief summary of extant literature on subdivision surfaces. Section 4.3 details

the crux of this chapter: (i) defines basis functions on the subdivision surface, (ii) develops

methods to address low frequency breakdown and multiplicative Calderón preconditioners

for the EFIE. Section 4.4 presents several numerical results that verify the accuracy and

efficacy of this approach. Finally, Section 4.5 summarizes the contribution of the work as

well as outlines directions for future research.

4.2

Integral Equations for Electromagnetic Scattering

The model problem for analysis will be a perfect electrically conducting object that

occupies a volume Ω whose surface is denoted by ∂Ω.

It is assumed that this surface is

equipped with an outward pointing normal denoted by ˆn(r). The region external to this
volume (R3\Ω) is occupied by free space. It is assumed that a plane wave characterized by

44

{Ei(r), Hi(r)} is incident upon this object. The scattered field for r ∈ R3\Ω can be obtained
using equivalence theorems as follows:

where

T ◦ J(r)

ˆn(r) × Es(r) = T ◦ J(r)
ˆn(r) × Hs(r) = K ◦ J(r)

(cid:90)

dr(cid:48)(cid:20)

(cid:90)
= −jωµ0ˆn(r) ×
.
= ˆn(r) × ∇ ×
.

I +

1

κ2∇∇

∂Ω

dr(cid:48)g(r, r(cid:48))J(r)

(4.1a)

(cid:21)

g(r, r(cid:48)) · J(r(cid:48))

(4.1b)

K ◦ J(r)

where g(r, r(cid:48)) = exp(cid:2)−jκ|r − r(cid:48)|(cid:3) /(4π|r − r(cid:48)|), κ is the wave number in free space, J(r) is

∂Ω

the equivalent current that is induced on surface, and I is the idempotent. In the above
expression and in everything that follows, we have implicitly assumed (and suppressed)

an exp [jωt] time dependence. Using the above equations, one may prescribe the requisite

integral equations as

EFIE:

ˆn(r) × ˆn(r) ×(cid:16)
ˆn(r) ×(cid:16)

Ei(r) + Es(r)

(cid:17)
= 0 ∀r ∈ ∂Ω−

= 0 ∀r ∈ ∂Ω

(cid:17)

MFIE:

Hi(r) + Hs(r)

(4.2a)

(4.2b)

that are known as the electric/magnetic field integral equations (EFIE/MFIE). In (4.2b),
∂Ω− denotes a surface that is conformal to but just inside ∂Ω. These equations do not have
unique solutions at the so-called irregular frequencies, but may be combined to yield the

combined field integral equation [78] that gives unique solution at all frequencies. It should

be noted that while these are the most popular formulations that are used in practice, they

are not the only ones. Other formulations such as the combined source integral equation [79],

augmented EFIE, augmented MFIE [80] and the charge current integral equation [81] exist

and have seen recent development. However, this chapter’s focus will be on the EFIE and

its discretization. The choice is largely motivated by the numerous challenges that exists in

solving these equations, both at the mid and low frequency regimes. To this end, in what

follows, we will develop basis functions that are well formulated and present modifications

45

to the formulation to account for both low frequency behavior and to impose Calderón

preconditioners.

Methods for solving these equations follow the standard prescriptions: (i) represent the

surface of the scatterer using some tessellation, (ii) define basis function on this approxi-

mation to the geometry, (iii) demonstrate desirable properties of these basis sets, and (iv)

validate solutions to integral equations solved using this procedure. As elucidated earlier,

as opposed to classical tessellation, we use subdivision surfaces for geometric representation.

This is elucidated next.

4.3 Current Representation and Field Solvers on Subdivision Sur-

faces

4.3.1

Iso-geometric Basis Sets

Section 3.2 developed a framework wherein one dealt with the representation of the limit

geometry surface. Using the same representation as presented earlier, we introduce the

concept of scalar functions that are defined on the limit surface. As with geometry, the

“limit” function can be represented in terms of 1-ring control weights of a regular triangle

and the box splines as

f (r) = f (r(u, v)) =

12(cid:88)

i=1

Ni(u, v)wi; for r ∈ T

(4.3)

where Ni(u, v) denotes the box spline function [49]. Several properties of underlying basis
functions for the geometry are worth noting as they are critical to represent both the ge-

ometry as well as the physics defined on the geometry. Note, some of these properties are

available for a NURBS description as well in the interior of the element, but not necessarily

across domains. These include the following:

1. compact support,

2. non-negativity,

46

3. convexity preserving,

4. and C2 continuity across patch boundaries.

The above properties make this subdivision basis functions a good candidate for modeling

physics on both regular and irregular (with the help of subdivision scheme) triangular meshes.

In the above, we have assumed that the mapping r(u, v) exists, where (u, v) are the

barycentric coordinates of a triangle. In what follows, this assumption will be implicit, and

only r will be used. It can be readily verified that if only one vertex has a weight of unity

and other zeros, one would immediately get a smooth (up to 2nd order continuity globally)

scalar function. As a result, one can associate an effective basis function with every vertex

such that the limit function

f (r) =

NV(cid:88)

anξn(r)

(4.4)

n=1

where NV is the number of vertex and ξn(r) is an effective basis function that describes the
influence of the scalar quantities associated with a vertex. Fig. 4.1a gives an example of

the scalar basis function used for formulating iso-geometric analysis on top of subdivision-

described surface. The basis function is associated with an irregular vertex of valence 8. To

demonstrate 2nd order continuity, the surface Laplacian of the scalar function is plotted in
Fig. 4.1b. The behavior of the subdivision basis and its derivatives follow: ξn(r) ≈ O(1),
∇sξn(r) ≈ O(1/h) and ∇2
sξn(r) ≈ O(1/h2) where h is an approximate dimension of the
patch. Two salient properties of the scalar subdivision basis used for this definition are:

(i) as they rely on approximating subdivision they are C2 almost everywhere; and (ii) if
Ω = ∪nΩn, where Ωn is the domain associated with ξn(r), then the function ξn(r) vanishes
on ∂Ωn. In other words, the basis functions ξn(r) ∈ C2

0 almost everywhere.

Next, to use this basis set to represent the current, we note that current on any surface

can be represented using a Helmholtz decomposition as follows:

J(r) = ∇sφ(r) + ∇ × (ˆnψ(r)) + 
(r)

(4.5)

47

(a)

(b)

Figure 4.1: Scalar basis function associated with a vertex with valence of 8.
function, (b) Surface Laplacian of the basis function

(a) Basis

where 
(r) is the harmonic field and has a trivial solution (
(r) = 0) for simply-connected

structures, and φ(r) can be associated with the charge space. Following the subdivision

48

(cid:88)
(cid:88)

n

φ(r) ≈ ˜φ(r) =

ψ(r) ≈ ˜ψ(r) =

a1,nξn(r)

a2,nξn(r)

n

(cid:104)

a1,nJ1

(cid:88)
n = ∇sξn(r)
J1
n = ˆn × ∇sξn(r)
J2

n

(cid:105)

(4.6)

(4.7)

representation in (4.4), we will assume that ˜φ(r) and ˜ψ(r) represent approximations to φ(r)

and ψ(r), respectively, and they can be represented in a manner similar to (4.4); viz.,

Using (4.6), it is now possible to define the approximations to the current as

J(r) ≈ ˜J(r) =

n(r) + a2,nJ2

n(r)

The physical interpretation of this representation is akin to a standard subdivision represen-

tation of a limit surface; we represent the the “limit” current via the “limit” scalar functions.

We note, that the basis for the currents (and the auxiliary potentials) are represented via

operations on the subdivision basis which can be effected numerically rather trivially. In this

work, the standard definition of the inner product (cid:104)X, Y(cid:105) =(cid:82)

Ω drX · Y is used.

Several properties of the basis functions make this definition appealing. These are as

follows:

Continuity Due to the reliance on functions ξn(r) that are C2

0 almost everywhere, the
resulting basis functions have C1 continuity almost everywhere and C0 continuity

at isolated points. This stands in stark contrast with classical Rao-Wilton-Glisson

(RWG) [52] or their higher order counterparts [57] that are div-conforming but not C0

at the boundary.

Orthogonality The inner product of (cid:104)J1

n(cid:105) = 0. The proof for this assertion can be
trivially obtained using Green’s theorems together with properties of the surface curl
and the fact that ξn(r) = 0 for (r) ∈ ∂Ω. Specifically,
n(r) −

drξn(r)∇s · J2

drξn(r)ˆu · J2

n(r) · J2

n(r) =

n, J2

drJ1

n(r)

(cid:73)

(cid:90)

Ωn

∂Ωn

(cid:90)

Ωn

(4.8)

= 0

49

where ˆu is the normal of ∂Ωn tangent to the surface.

Representation and Convergence While the use of these basis for geometry exists [49],

the convergence of the subdivision-enhanced basis function for scalar fields was pre-

sented in the framwork of scalar finite element method [82]. The basis functions used

to represent the currents are orthogonal, with the nonsolenoidal component of the cur-
rent, ∇sφ(r), is only represented through J1
n(r), and the solenoidal part is represented
n(r). The Hemlholtz decomposition used here is complete. As a result, the
through J2

tangent spaces show properties presented in [51, 83–85].

Charge neutrality The basis functions maintain charge neutrality. This is proved using

(cid:90)

Ω

dr∇s · ˜J(r) =

=

(cid:90)
(cid:73)

a1,n

a1,n

(cid:88)
(cid:88)

n

n

dr∇s · ∇sξn(r)

Ωn

∂Ωn

drˆu · ∇sξn(r)

(4.9)

= 0 as ∇sξn(r) = 0 ∀r ∈ ∂Ωn

Compactness The definition of basis function is local.

Refinement The basis functions are subdivision based. As a result, they inherit properties

of subdivision representation including adaptivity.

These properties ensure a complete representation of the currents on the surface of an

simply connected object, and forms a rigorous Helmholtz decomposition of currents on the

surface. Again, as opposed to classical basis functions defined on tessellation, Helmholtz

decomposition is inherent in the definition of basis functions where as the common approach

is to design “quasi”-Helmholtz decomposition using weighted sum of basis functions. It should

be noted that, in the above representation, the potential functions are the actual unknown

functions and not the currents. As a result, extra steps have to be taken to ensure their

uniqueness. These are elucidated below when using these basis functions within a Galerkin

framework.

50

Using these definition of basis functions, each vertex is associated with two degrees of

freedom (DoF). Therefore the number of DoFs is twice the number of vertices. This situation

is quite different from the hierarchical higher order basis functions that involve polynomials

of different degrees (the polynomial degree used is fixed in the proposed IGA.). Fig. 4.2

shows the behavior of the two basis functions associated with a node of valence 8. A point

to note is that the number of solenoidal and non-solenoidal basis functions are equal to each

other. As we will see, this helps creating Calderón preconditioners. It is also worth noting

that the proposed tangent basis set only works for closed surface, since extra requirements,

such as no boundary charge, need to be imposed on open edges. The extension to open

structures that is nontrivial will be studied in the future work.

(a)

(b)

Figure 4.2: Basis function associated with a vertex with valence of 8. (a) nonsolenoidal type,
(b) solenoidal type

51

4.3.2 Field Solvers

Given the prescription of basis functions, the discretized version of (4.2a) can be obtained

using Galerkin testing. The resulting matrix equation can be written as

where

 Z11 Z12
(cid:90)

Z21 Z22



 I1
(cid:90)

I2

Zlk
nm = jωµ0
− jδl1δk1
ωε0

(cid:90)

Ωn

Ωn

drJl

n(r) ·
dr∇s · J1



m(r(cid:48))

V2

 V1

 =
(cid:90)
dr(cid:48)g(r, r(cid:48))Jk

Ωm

Ωm

n(r)

(cid:90)

Ωn

dr(cid:48)g(r, r(cid:48))∇s · J1

m(r(cid:48))

(4.10a)

(4.10b)

(4.10c)

Ik
n = ak,n;Vk

n =

drJk

n(r) · Ei(r)

In the above equations, δij denotes a Kronecker’s delta. Interesting features of the above
expression are apparent; (i) Z11
nm is the only term that has the charge contribution. As will be
shown later, this “decoupling” is an essential component for the construction of low frequency

stable solvers, and (ii) since the system of equations are constructed using conditions on

currents that rely on derivatives of the potential ˜φ(r) and ˜ψ(r), it follows that the potentials

are determined upto a constant, and the number of degrees of freedom is less by one for

each potential. While one can impose this via different means, we have chosen to essentially

constrain the system by choosing a1,N and a2,N to be zero. This implies a trivial change to
the system of equations (4.10). The evaluation of the inner products in the above equations

is effected via higher order quadrature and Duffy integration rules. While the approach

presented thus far is valid for all frequencies, it suffers from low frequency breakdown. In

what follows, we demonstrate that (4.10) can be modified trivially to alleviate low-frequency

breakdown. Likewise, Calderón preconditioners can be easily implemented using this space

of basis functions.

But before we proceed, interesting insight may be gained by examining the inner products

that arise using the Galerkin procedure. Matrices arise from testing the electric field, i.e.,

52

(cid:111)

(cid:110)

Jk
m

.
= E

by Jl

n. If

Zlk
nm corresponds to measuring the radiated electric field due to Em(r)
follows that

(cid:90)

(cid:90)
(cid:73)

Ωn

drJ1

n(r) · Em(r) =

=

dr∇sξn(r) · Em(r)
drξn(r)ˆu · Em(r) −
drξn(r)∇s · Em(r)

Ωn

(cid:90)

∂Ωn

= −

Ωn

(cid:90)

Ωn

drξn(r)∇s · Em(r)

(4.11a)

(cid:90)

Ωn

drJ2

n(r) · Em(r) =

=

=

=

=

(cid:90)
(cid:90)
(cid:90)
(cid:90)
(cid:73)

Ωn

Ωn

Ωn

Ωn

drˆn × ∇sξn(r) · Em(r)
drˆn × ∇ξn(r) · Em(r)
drˆn · (∇ξn(r) × Em(r))
drˆn · ∇ × (ξn(r)Em(r)) −
(cid:90)
drξn(r)Em(r) · ˆt + jωµ0

drξn(r)ˆn · Hm(r)

(cid:90)
(cid:90)

Ωn

Ωn

∂Ωn

= jωµ0

Ωn

(4.11b)

drξn(r)ˆn · ∇ × Em(r)
drξn(r)ˆn · Hm(r)

In the above equation, ˆt is a unit vector tangential to the boundary ∂Ωn, and Hn(r) is the
magnetic field due to Jk
m. From the above equations, it is apparent that (4.11a) tests the
tangential component of the electric field. However, (4.11b) yields equations that test in

normal component of the magnetic field. So the two basis functions used in the analysis

impose tangential continuity of the electric field as well as the normal component of the

magnetic field. As a result, using a basis that satisfies the Helmholtz decomposition (and
0), results in equations that naturally fit into the frame-
work of the Current-Charge Integral Equations (CCIE) [81]. It should also be noted that

relies on scalar function that are C2

while eqns. (4.11a) and (4.11b) were specified for the scattered field, they are equally valid

for the incident field. Indeed, using the final expression in (4.11b) to evaluate the integrals

is more accurate and requires fewer quadrature points.

53

4.3.3 Low-frequency Stable EFIE

It is well known that the EFIE suffers from low-frequency breakdown. The rationale for this

breakdown is readily apparent by examining the components of the impedance matrix in

(4.10). Assuming that the average largest linear dimension of the support of the patch is h.

It follows that the entries corresponding to

Zlk
nm = O(κh2) + δl1δk1O(1/κ)

(4.12)

These results follow from the fact that the functions ξ(r) = O(1). The above scaling indicates
that a portion of the elements associated with source and test basis functions that are
n(r) scale as O(1/κ), whereas all others scale as
associated with irrotational functions J1
O(κh2). This implies that as κ −→ 0, the portion of Z11
nm that corresponds to the charge
contribution dwarfs the rest. This situation is similar to those encountered for by classical

Nedelec elements. Here, one takes recourse to loop-star/loop-tree decompositions [86–88]

that effect an approximate Helmholtz decomposition of the currents. As has been shown,

the resulting decompositions contain a portion that is exactly divergence free and one that

is approximately curl free. Whereas the support of the divergence-free portion is local, the

same is not true of curl-free portion. This is in contrast to the method presented here wherein

both components have local support. Using these basis results in the charge being modeled

correctly, and a matrix whose scaling looks like (4.12). Rescaling these equations has been

shown to render the solution stable. Following a similar procedure, it can be shown that

rescaling both the matrix elements, the coefficients and the right hand side results in

(cid:90)

(cid:90)

Ωm

Zlk
nm = jωβlkµ0
− jδl1δk1

Ωn

drJl

n(r) ·
dr∇s · J1

ε0

dr(cid:48)g(r, r(cid:48))Jk

m(r(cid:48))

dr(cid:48)g(r, r(cid:48))∇s · J1

m(r(cid:48))

n(r)

Ωm

where

(cid:90)
(cid:90)

Ωn

βlk =



ω
ω−1

1

if l = 1, k = 1

if l = 2, k = 2

otherwise

54

(4.13)

(4.14)

n = ak,nω−δk1;Vk
Ik

n = ω−δk2

(cid:90)

Ωn

drJk

n(r) · Ei(r)

(4.15)

As prescribed, these equations achieve the desired stability. As can be trivially shown, these
equations decouple as ω −→ 0. This is akin to similar prescriptions for classical loop-star/tree
algorithms [86,87,89]. However, a salient feature of the IGA-basis is that constructing a well

behaved system is tantamount to using a diagonal preconditioning sans constructing the

complementary system that is required for a loop-tree/star (Hodge) decomposition.

It should be noted that Helmholtz decomposition is not the only way to solve low fre-

quency breakdown, and other techniques [81,90–93] exist. The iso-geometric basis functions

presented herein could help construct those methods.

4.3.4 Calderón Preconditioner

It is well known that the standard EFIE presented in the earlier sections is a first kind

equation, and as with all first kind equations, can be ill-conditioned especially when the

spatial scales in the problem are widely separated. As has been shown by several others (see

Refs. [51,94,95], and references therein), the rationale for the ill-conditioning is the spectral

separation of the eigenvalues of the operator with increase in discretization density, with a set

that clusters around zero and others that cluster around infinity. Given this separation, the

resulting matrix systems become rapidly ill-conditioned. The remedy to this problem exploits
the Calderón identities wherein the EFIE operator (T {·}) preconditions itself resulting in a
second kind integral operator whose eigenvalues accumulate around −1/4 [51, 96, 97]. The
challenge in using these identities was the lack of well behaved Gram matrices that link
the domain and range of the T operator [97, 98]. For Thomas-Raviart/Rao-Wilton-Glisson
basis sets and C0 geometries, basis functions developed by Buffa and Christiansen [99] have

been exploited in a sequence of papers to thoroughly understand and solve this problem. In

what follows, we show that basis functions defined herein result in a well conditioned Gram

matrix, thereby permitting a natural discretization of the Calderón operator. To begin, we

55

T ◦ J(r) = Ta ◦ J(r) + Tφ ◦ J(r)

(cid:90)
(cid:90)

Ta ◦ J(r)
Tφ ◦ J(r)

= −jωµ0ˆn(r) ×
.
= −j
ˆn(r) ×
.

1
ωε0

∂Ω

dr(cid:48)g(r, r(cid:48)) · J(r(cid:48))
dr(cid:48)∇∇g(r, r(cid:48)) · J(r(cid:48))

(4.16a)

(4.16b)

define operators

where

∂Ω
The Calderón projector is defined as T 2 ◦ J(r)
= T ◦ T ◦ J(r).
.
It has been shown that
T 2 ◦ J(r) = (−1/4 + K2) ◦ J(r), where K2 ◦ J(r)
= K ◦ K ◦ J(r) is compact, and the
.
eigenvalues are clustering around −1/4. Alternatively, this operator may also be written

as T 2 ◦ J =(cid:0)Ta ◦ Ta + Ta ◦ Tφ + Tφ ◦ Ta + Tφ ◦ Tφ

(cid:1) ◦ J. Since direct discretization of T 2 is

impossible, typical implementation follows the multiplier approach; i.e., define intermediate
mapping from the range space of the T operator to the domain of the T . This is typically
effected via a Gram matrix and has been extensively explored. Given the usage of Helmholtz

decomposition (4.7) and the orthogonality (4.8) between the two components, it can be
shown that Tφ ◦ Tφ ◦ J = 0. Enforcing the latter condition has been difficult to accomplish
in traditional discretization and representation schemes.

n, Jk

The principal advantage of using (4.7) is that it enforces an exact global Helmholtz
0 almost everywhere.
decomposition in terms of a set of local scalar functions that are C2
m(cid:105) = γnmδlk, where γnm is the result of evaluating the inner
Recall the following: (cid:104)Jl
product over the support of basis functions Ωn ∩ Ωm and is defined in (4.18). As a result,
the Gram matrix G is block diagonal with entries Glk
nm = γnmδlk. Given that, for a pair

of spaces, (J1, J2), it can be shown that the following sequence is satisfied: (cid:0)J1, J2(cid:1) T−→
(cid:0)J2, J1(cid:1) T−→(cid:0)J1, J2(cid:1). Specifically,(cid:0)J1, J2(cid:1) Tφ−−→(cid:0)J2, 0(cid:1) and(cid:0)J1, J2(cid:1) Ta−−→(cid:0)J2, J1(cid:1). It follows
from the above that(cid:0)J1, J2(cid:1) Tφ−−→(cid:0)J2, 0(cid:1) Tφ−−→ (0, 0). In matrix form, the above sequence can

be rewritten as

T G−1T =



a

T 11
T 21
a + T 21

φ

a

T 21
T 22

a



 G11

0

−1

0
G22

56



a

T 11
T 21
a + T 21

φ

a

T 21
T 22

a

(4.17)

Finally, a critical component in discretizing the T 2 operator is the Gram matrix. The
benefits of using local Helmholtz decomposition becomes apparent in its construction with

the off-diagonal blocks being zero and the element of the diagonal blocks being

dr∇sξn(r) · ∇sξm(r)

(4.18)

(cid:90)

γnm =

Ωn∩Ωm

that are a variational form of the Laplace-Beltrami operator on a scalar function. The

resulting system is positive definite leading to well behaved inverse of the Gram system with

diagonal preconditioner or simply.

4.4 Numerical Examples

This section presents several numerical examples to demonstrate the efficacy of the pro-

posed approach to electromagnetic analysis. In order to do so, we shall present data that

demonstrates the following; (i) accuracy of the approach to when compared against analytical

data; (ii) examples illustrating low frequency stability, (iii) analysis of multiscale structures

to illustrate the efficiency of the Calderón preconditioner, and (iv) application to complex

targets,. Unless stated otherwise, the data presented compares in radar scattering cross-

section computed using the IGA-MoM solver and those computed using either analytical

methods (if available) or a validated method of moments code that is based on RWG basis

functions.

4.4.1 Accuracy of IGA-MoM

To validate and demonstrate the accuracy of the proposed approach, we consider scattering

from a sphere discretized at multiple resolutions. To this end, consider a sphere with radius

of 1.0 meter that is modeled using an initial control mesh comprising 642 vertices. In all the
data presented below, a plane wave field propagating in ˆk = −ˆz and polarized along ˆx axis
is incident on the sphere.

First, the performance of the proposed vector basis functions is tested by reconstructing

the surface currents due to the incidence plane waves at two different frequencies, 200MHz

57

Table 4.1: Rel. L2 error in the reconstructed surface currents density

subdivision times
No. of elements

IGA
RWG
IGA
RWG

200MHz

300MHz

0

1280

2.781e-2
9.239e-2
1.072e-1
1.464e-2

1

5120

2.292e-3
4.456e-2
7.027e-3
7.137e-2

2

20480
2.365e-4
2.218e-2
6.633e-4
3.556e-2

3

81920
3.082e-5
1.113e-3
7.915e-5
1.783e-2

4

327680
5.037e-6
5.604e-3
1.086e-5
8.946e-3

and 300MHz. At each frequency, the surface is subdivided four times. That would result

in five sets of control meshes with face number of surface elements ranging from 1280 to

32, 7680. Note, the limit surface for all three mesh densities is identical. As a result, what

changes is the support of the basis functions, and therefore, approximations to the discrete

operators. The convergence of the relative L2 errors versus the number of surface elements
in the control mesh is given in the table 4.1. The same test was also implemented using the

regularly used div-conforming RWG basis. Next, the surface current excited by the incidence

waves is solved at 300MHz. Likewise, starting from the control mesh and performing one and

two Loop subdivision results in two other meshes with 2562 and 10242 vertices, respectively.

The number of degrees of freedom (DoF) for IGA-MoM is twice the number of vertices;

consequently, the three tests have 1284, 5124, and 20484 DOFs, respectively. By comparing

those obtained using IGA-MoM with Mie series approach, the relative errors in surface

currents are 0.0257, 0.0067 and 0.0030, respectively, for these three different discretizations.

Fig. 4.3 plots the current distribution, with its radar cross section (RCS) result given in

Fig. 4.4. Figs. 4.5a and 4.5b present pointwise errors in absolute values of the real and

imaginary parts of the current for the case of two subdivisions. As is evident, the accuracy

in the current both in the L2 and L1 norms is excellent.

4.4.2 Scattering from Structures at Regular Frequency

In what follows, we analyze the performance of the proposed approach and compare these

with a conventional MoM EFIE solution technique that relies on RWG basis functions. Un-

58

less specified otherwise, the metric for comparison is RCS data in the φ = 0 plane. Further-

more, the geometry of all scatterers that are analyzed in this chapter is smooth; extension

to non-smooth structures is non-trivial but realizable from a subdivision perspective [100].

IGA solvers that include in open and sharp surfaces will form the basis of a future work.

The first example is a truncated cone, with the height 3.0λ and the radii of the top

and bottom circular cross-sections being 0.2λ and 1λ, respectively. A plane wave that is
polarized along ˆx and propagating in the −ˆz direction is incident on the object. The object
is represented using 12378 nodes; this number is increased at the top and bottom surfaces

so as to maintain a sufficient degree of sharpness. The number of DoFs for IGA-MoM and

conventional MoM are 24756 and 37128, respectively. The surface current density obtained

using IGA-MoM is depicted in Fig. 4.6.

It is evident from these plots that the currents

are smooth (without unphysical aberrations anywhere, especially near the crease). Fig. 4.7

compares the radar cross section between the proposed method and the conventional method

of moments, and it is apparent that the agreement between the two is excellent.

The second example is a structure composed of a block (3.33λ × 10.0λ × 1.33λ) and five
cylinders (radius is 0.5λ and the height is 0.67λ) uniformly distributed on the top surface
of the block. An electromagnetic field that is propagating along −ˆy and polarized along the
ˆz incident on this object. The number of DoF for the IGA-MoM and conventional MoM

are 29028 and 43536, respectively. Figure4.8 presents the real and imaginary parts of the

current distribution on the surface of the object, and again, it is evident that the results are

smooth without artifacts. Further, excellent agreement can be seen between IGA-MoM and

conventional MoM in the RCS data.

As a final example in this set, we analyze scattering from a model airplane whose di-
mensions are 5.49λ × 5.48λ × 1.52λ. The plane wave incident on the object is propagating
along −ˆy and polarized along ˆz. As before, data is obtained using both IGA-MoM and
conventional MoM using 23988 and 35976 DoFs, respectively. The current distribution using

IGA-MoM is depicted in Fig.4.10. Again, it is evident that the current is well captured as

59

is the scattering cross-section (see Fig. 4.11).

4.4.3 Conditioning of iso-geometric System at Low-frequency

Next, we discuss the stability of the resulting system (and its preconditioned form) at low

frequency using IGA basis and the aforementioned wavenumber scaling. The metrics used

for the discussion are both the condition number of the resulting system and the number

of iterations necessary to achieve a desired error tolerance when solving this system using

an iterative solver.

In all the examples presented, we use a generalized minimal residual

(GMRES) iterative solver and the object being analyzed is a sphere with radius 1m that is

discretized such that the average distance between vertices is around 0.15m.

Fig. 4.12 depicts the condition number of IGA-MoM systems as a function of frequency

from 1Hz to 100MHz. The three curves are, respectively, for the original IGA system, the

diagonal preconditioned system, the Calderon preconditioned system. From the plots, it is

evident that the condition number is almost constant over the whole band for all the systems.

Both diagonal preconditioner and application of the Calderón precondition offer significant

improvements over the original IGA-MoM. Note, the behavior of the original IGA-MoM is

very much unlike conventional MoM, thanks to the properties of the basis function used.

Next, we examine the number of GMRES iterations required to converge to an error
tolerance of 1.0× 10−6 and 1.0× 10−10 at various frequencies. These results are depicted in
Figs. 4.13 and Fig. 4.14, respectively, both sampled at 1Hz, 100Hz, 104Hz, 106Hz and 108Hz.

Since the iterative numbers for the original IGA system are much higher than that of the

diagonal preconditioned and Calderon preconditioned system, only the iteration numbers for

the latter two systems are given. As seen in the two plots, the number of iterations required

for both systems are relatively constant across a range of frequencies (from 1Hz to 100MHz).

These two tests serve to illustrate some of the salient features of the proposed approach.

(i) Low frequency stability is achieved simply by diagonal preconditioning. This is in contrast

60

with conventional approaches where it is explicitly necessary to define auxiliary unknowns

albeit at linear cost. (ii) Imposition of the Calderón preconditioner involves an inversion of

the Gramm matrix as an addition operation. The original matrix system is retained which

makes integration with fast methods trivial. This is different from existing methods that

require Buffa-Christiansen basis functions to be defined on barycentric meshes. Finally, as

is evident from Fig. 4.15, the results obtained from the Calderón preconditioned solver at

1Hz agrees very well with analytical data obtained.

4.4.4 Examples with Multi-scale Mesh

In all the examples analyzed thus far, the initial control mesh is such that the resulting

tesselation is almost uniform in that the ratio of the maximum edge length to the minimum
edge length is O(1). However, a more intellectually interesting and practically applicable
problem is when this ratio is significantly higher. It is challenging to design stable methods

for these problems due to two effects that act in concert with each other; wavenumber and

element size scaling. In this example, scattering from a sphere with locally refined meshes is

simulated, and verified by comparing with the one with almost uniform initial control mesh.

Fig. 4.16 demonstrates the multiscale control mesh for the unit sphere. The radius of the

sphere is 1λ, and the incident plane wave is polarized along the ˆx direction and propagating
in the −ˆz direction. The resulting system of equations is solved using GMRES. Figure 4.17
compares the convergence history for both the Calderón preconditioned system as well as

a diagonal preconditioned system. The efficacy of the Calderón preconditioned system is

evident by the rapid convergence.

4.5 Conclusion

In this chapter, we have developed a novel iso-geometric analysis technique for solving

the electric field integral equation that is encountered in electromagnetic field analysis. The

fundamental thesis of iso-geometric methods is that they inherit the rather significant ad-

61

vantages of modern CAD tools–smoothness of geometry representation, morphing, dynamic

meshes (if necessary), etc–by using the same basis sets for both geometry representation as

well as representation of physics on this geometry. The choice of the electric field integral

equation is predicated upon the fact that it is one of the more challenging equations to solve.

The approach presented here leverages existing geometric construction techniques that use

subdivision to define basis functions that result in well behaved integral operators. Thanks

to these operators, it is possible to trivially modify these to impose Calderón preconditioners,

and construct systems that are low frequency stable and can handle multiscale geometric

features. A number of results demonstrate the convergence and accuracy of the technique,

applicability to computation of scattering at both regular and low frequencies, as well as

to structures with multiscale features. It is very simple to apply the proposed basis set to

other types of integral equations, such as magnetic field integral equation and combined

field integral equation. However, as with the introduction of any new technique several open

problems remain: these include addition of features to handle edges and open domains, in-

tegration with fast solvers, development of techniques for multiply connected objects, using

this framework to include Debye sources, etc. Several of these topics are active areas of

research within the group and will be presented in different forums soon.

This chapter, c(cid:13)2016 Elesvier, is reprinted, with permission, from Li, J.; Dault, D.; Liu,
B.; Tong, Y. and Shanker, B., "Subdivision based isogeometric analysis technique for electric

field integral equations for simply connected structures," Journal of Computational Physics,

August 2016, with slight modifications to fit in this dissertation.

62

(a)

(b)

Figure 4.3: Magnitude of surface currents density on the sphere: (a) real part and (b)
imaginary part

63

Figure 4.4: Radar cross section of the sphere (φ = 0 cut)

64

-20-10 0 10 20 30-150-100-50 0 50 100 150RCS/dBsmScattering angle/degIG-MoMRWG-MoMMie(a)

(b)

Figure 4.5: Pointwise relative error (a) real part, (b) imaginary part with 20484 DoFs

65

(a)

(b)

Figure 4.6: Magnitude of surface current density on the truncated cone: (a) real part and
(b) imaginary part

Figure 4.7: Radar cross section of the truncated cone (φ = 0 cut)

66

-30-20-10 0 10 20 30 0 50 100 150 200 250 300 350RCS/dBsmScattering angle/degIG-MoMRWG-MoM(a)

(b)

Figure 4.8: Magnitude of surface current density on the composite structure: (a) real part
and (b) imaginary part

Figure 4.9: Radar cross section of the composite structure (φ = 90 cut)

67

 20 30 40 50 60 70 80 90 100 0 50 100 150 200 250 300 350RCS/dBsmScattering angle/degIG-MoMRWG-MoM(a)

(b)

Figure 4.10: Magnitude of surface current density on the air plane model: (a) real part and
(b) imaginary part

68

Figure 4.11: Radar cross section of the air plane model (φ = 90 cut)

Figure 4.12: Condition number at different frequencies

69

 0 10 20 30 40 50 0 50 100 150 200 250 300 350RCS/dBsmScattering angle/degIG-MoMRWG-MoM1e+021e+041e+061e+081e+001e+021e+041e+061e+08Condition NumberFrequency/Hzoriginal IGA systemDiagonal Precontioned systemCalderon Preconditioned systemFigure 4.13: Number of GMRES iterations to converge to a tolerance of 1.0 × 10−6

Figure 4.14: Number of GMRES iterations to converge to a tolerance of 1.0 × 10−10

70

 0 200 400 600 800 10001e+001e+021e+041e+061e+08IterationsFrequency/HzCalderon PreconditionerDiagonal Preconditioner 0 200 400 600 800 10001e+001e+021e+041e+061e+08IterationsFrequency/HzCalderon PreconditionerDiagonal PreconditionerFigure 4.15: Comparison of RCS data obtained at 1Hz with analytical solutions

(a)

(b)

Figure 4.16: Multiscale control mesh: (a) the whole mesh and (b) locally refined region

71

-350-340-330-320-310-300-290-280-150-100-50 0 50 100 150RCS/dBsmScattering angle/degIG-MoMMieFigure 4.17: Convergence of the iterative solver

72

Iteration Number05001000150020002500300035004000Relative Residual10-710-610-510-410-310-210-1100CalderondiagonalCHAPTER 5

SCALAR INTEGRAL EQUATIONS FOR ELECTROMAGNETICS

5.1

Introduction

Surface integral equation based method have long been considered an effective and effi-

cient approach to study electromagnetic scattering and radiation [101]. The integral equa-

tion enables the use of surface discretization (hence fewer degree-of-freedom) and its natural

marriage with fast divide-and-conquer methods makes it even more popular for engineering

analysis. Recent research efforts made by the computational electromagnetics community

in surface integral equations can be grouped into categories including fast methods [102],

remedies for low-frequency breakdown [87, 91, 103–105], basis functions and higher order

methods [52, 56, 99, 106], new formulations [79, 80, 91, 103, 107, 108], flexible geometry rep-

resentation for design purpose [3, 5, 45], preconditioning techniques [96, 97, 109] and so on.

Among them, research in stable formulations, geometry representation and basis functions is

of significant importance since it can affect many aspects of the integral equation solvers. For

examples, well-conditioned system would definitely enhance the iteration based fast meth-

ods; good geometry representation can help the development of high order basis functions

and speed up the design and optimization processes. Discussions in this chapter will focus

on these aspects: scalar formulation, geometry representation and iso-geometric basis that

is inspired by the representation of the geometry.

Most commonly used formulations, like electric field integral equation (EFIE), use the

field representations derived from equivalent surface current densities (this is also called direct

approach in mathematics community). In these integral equations, two problems arise: (1)

unbalanced frequency dependence of the contributions from scalar and vector potentials (the

low-frequency breakdown) and (2) div-conforming requirement on basis functions for those

currents. Even in some stable integral equation formulations, the post-processing stage

73

usually involves low-frequency breakdown. Studies focusing on low frequency breakdown

include quasi-Helmholtz decomposition based methods [87, 105, 110] and new formulations

such as current-charge integral equations (CCIE) [108] and augmented EFIE (AEFIE) [90].

Another approach [5] is based on the use of scalar potential as the unknown to implement

the rigorous Helmholtz decomposition (this is cast within an iso-geometric framework). Both

CCIE and AEFIE work with extra continuity equation after introducing auxiliary charges

in the system. In all of the aforementioned research, save for [5], explicit unknown currents

are still involved. In these direct approaches [111], tangent and/or normal components of

the field are considered as the observation field. So a natural question is whether there

is another set of integral equations that can be derived. Other remedies, using indirect

approach with unknowns that don’t usually have physical meaning, include generalized Debye

sources approach [91, 112, 113] and the more recent work on decoupled potential integral

equations [103,114]. The former leads to a second kind integral equation that is stable over

a broadband, which only involves scalar unknowns. The latter solves the problem in a totally

different perspective by re-writing the original Maxwell’s boundary value problem into one in

terms of vector and scalar potentials. The approach presented is similar to the first one, but

with different formulations (they are still categorized as direct formulations) and hence the

unknowns. Similar effort on the conventional EFIE with scalar basis function is presented

in [115], where the tangential component of the electric field is used as the observation field.

In addition to the formulation, discretization of the integral equation also matters in

aspects like accuracy and convergence. The discretization involves geometry representation

and the construction of basis functions on the geometry. Different formulation might impose

certain requirement on geometry and/or basis functions. Studies in these directions include

[5, 45, 56, 57, 116]. Most of the studies use piecewise continuous patch representation that

might again introduce further difficulties in defining basis or test functions. Furthermore, the

mesh on which physical quantities are sampled and represented is usually reduced from high-

fidelity geometry representations in computer-aided-design (CAD) software. In contrast, the

74

iso-geometric analysis [5, 44, 117] embeds more information about the geometry to define

the basis function on it, hence no re-meshing process will be needed. Without further

geometry accuracy loss, simulation on the native geometry can be more accurate and flexible

[5, 44, 117–119]. In [5], scalar potentials are used as the unknown quantity defined on the

subdivision surface and are used to solve the conventional EFIE. However, naive extension to

magnetic field integral equation (MFIE) will not recover the identity plus compact operator.

The scalar integral equation formulations presented in this chapter allow straightforward

application of the iso-geometric to both integral equations. The almost C2 continuous scalar

basis function generated in the Loop subdivision surface [48,49] can be used to represent the

both the unknown scalar charges and the development of test functions.

The principal contribution of this chapter is theoretical development of scalar integral

equations and analytical proofs of stability over a wide frequency range. The results pre-

sented serve to validate this theme as opposed to presenting large scale analysis. To this

end, the contributions of the chapter are the following: We present

• a rigorous formulation of several stable scalar electromagnetic integral equations,

• analytical stability and spectral analysis over wide frequency band,

• iso-geometric implementation of the well-conditioned formulations on simply connected

surfaces,

• and numerical examples demonstrating the accuracy and stability of the iso-geometric

analysis of the scalar formulation.

The remainder of the chapter is organized as follows. Section 5.2 derives the scalar formu-

lation from different forms of boundary conditions. Section 5.3 discusses the iso-geometric

analysis implementation of the scalar formulations. In section 5.4, spectral properties of the

new formulations and performances of the IGA implementations are demonstrated through

numerical examples. Conclusion and related remarks are given in section 5.5.

75

5.2 Formulation

This chapter will focus on electromagnetic (EM) scattering from a simply-connected

perfect electric conductor (PEC) or a collection of simply-connected PEC objects. Without

losing generality, in what follows only one scatterer is considered in formulating the scalar

integral equations.

5.2.1 Problem Statement

Assume the scatterer enclosed by the surface S is illuminated by an incident EM waves
denoted by (Ei, Hi). Traditional approaches to solving for the scattered fields {Es, Hs} rely
on the surface equivalence theorem wherein one posits an equivalent current J that exists

on surface S. Then the total fields satisfy boundary conditions on the surface S. Typically,

these take the form of twisted tangential field components on the boundary, viz.,

n × Et = n ×(cid:16)
n × Ht = n ×(cid:16)

Ei + Es {J}(cid:17)
Hi + Hs {J}(cid:17)

= 0

= J.

(5.1)

(5.2)

and

The properties of these equations have been explored extensively [51,87,101,104,105,108,114]

with regard to their behavior at high and low frequencies, condition number, etc. These

equations are used to develop integral equations (that may or may not be well-conditioned).

In what follows, we explore an alternative approach.

5.2.2 Charge based Representation

To start, consider the standard mixed potential formulation for the scattered electric field

due to surface current densities J(r),

Es(r) = −jωµSk[J] +

∇Sk[∇(cid:48)

s · J],

1

jω


(5.3)

76

where the single layer potential operator is defined as

(cid:90)

S

Sk[J](r) =

G(r, r(cid:48))J(r(cid:48))dS(cid:48)

(5.4)

and G(r, r(cid:48)) is the free space Green’s function for the Helmholtz equation in free space.

As in [5], Helmholtz decomposition of the surface current can help one introduce scalar

sources into the system. Instead of using the potentials, one can introduce charges through

the Poisson equation on the surface. Therefore, given two charge sources, one can represent

the two generalized Debye potentials [91] as follows:

ψ = −jω∆−1

s ρ, φ = −jω∆−1

s ρm,

(5.5)

where ∆−1
the Helmholtz decomposition to the surface currents, one can write the currents with the

is the inverse surface Laplacian or Laplace-Beltrami operator. By applying

s

following representation.

J = −jω∇s∆−1

s ρ − jωn × ∇s∆−1

s ρm.

(5.6)

After replacing the currents in (5.3) with the above one, one can get the following repre-

sentation of the scattered field;

(cid:48)−1
Es = − ω2µSk[∇(cid:48)
s ρ]
s∆
− ω2µSk[n × ∇(cid:48)
s∆

(cid:48)−1
s ρm] − 1



∇Sk[ρ].

Similarly for scattered magnetic field, one obtains

and

Hs = ∇ × Sk[J]

(cid:48)−1
Hs = − jω∇ × Sk[∇(cid:48)
s ρ]
s∆
− jω∇ × Sk[n × ∇(cid:48)
s∆

(cid:48)−1
s ρm].

(5.7)

(5.8)

(5.9)

It’s worth noting that the new representations of both eletromagnetic fields are special cases

of the conventional equivalence theorem with native support of the Helmholtz decomposition.

77

These scalar source based representations, as well as formulations derived later, are different

from those in the indirect approach used in [91] where combined current sources are used.

Given these scalar sources, the question is whether one can identify the corresponding

scalar field quantities that can be used to construct integral equations. The starting point of
the analysis is inspired by [91], wherein divergence of a tangent electric field (∇· (n× n× E))
is chosen. It can be trivially shown that this implies conditions on the divergence of rotated

magnetic current. Here, we propose additional conditions on both the trace and the twisted

trace of both the electric and magnetic fields; viz:

of the electric field and

∇ · (n × n × E),
∇ · (n × E)

∇ · (n × H),
∇ · (n × n × H)

(5.10a)

(5.10b)

(5.11a)

(5.11b)

of the magnetic field.

The main goal of choosing the specific source-observable pair is to allow the use of scalar

basis functions and testing functions. Next, several independent sets of IEs derived from

boundary conditions associated with the above scalar observables will be introduced.

5.2.3 Scalar Electric Field Integral Equation

In PEC case, one can use (5.10) to construct a set of integral equations as in the EFIE.

While one can think of (5.10b) as a condition on surface divergence of the equivalent magnetic

current and (5.10a)a as a condition on the surface divergence of the rotated magnetic current,

it is also possible to show that (5.10a) is tantamount to imposing conditions on the normal

78

component of the electric field and its normal derivative via

∇ · (n × n × E) = −∇ · (E − nEn) = ∇ · (nEn)

= ∇ · nEn + n · ∇En
= JEn + n · ∇En,

(5.12)

where J is the first curvature of the surface.

By using the scalar representation developed earlier, one can get a set of integral equations

for the scattering problem.

(5.13a)

(5.13b)

(cid:48)−1
s ρ] − 1


(cid:48)−1
− ω2µ∇ · n × n × Sk[n × ∇(cid:48)
s∆
s ρm] = g1,

s∆

∇ · n × n ×(cid:2) − ω2µSk[∇(cid:48)
∇ · n ×(cid:2) − ω2µSk[∇(cid:48)

(cid:48)−1
s ρ] − 1


− ω2µ∇ · n × n × Sk[n × ∇(cid:48)
s∆
where g1 = −∇ · (n × n × Ei) and g2 = −∇ · (n × Ei).

s∆

∇Sk[ρ](cid:3)
∇Sk[ρ](cid:3)

(cid:48)−1
s ρm] = g2,

The above equations are the scalar equivalent of the electric field integral equations

(sEFIE) for solving for scattering from a PEC surface as opposed to the usual tangential

components of electric field. It should be pointed out that in this framework we would need

to explicitly impose a constraint on charge neutrality on the scalar sources. This is not done

in conventional framework that uses div-conforming basis.

5.2.4 Scalar Magnetic Field Integral Equation

A similar procedure can be used for the magnetic fields as well. Using the jump condition

on the trace and twisted trace one can obtain two equations

∇ · (n × Hs) + ∇ · (n × Hi) = ∇ · J,

∇ · (n × n × Hs) + ∇ · (n × n × Hi) = ∇ · (n × J)

(5.14)

(5.15)

79

These equations can be used to obtain the following scalar magnetic field integral equation

(sMFIE).

(cid:48)−1
ρ − ∇ · (n × ∇ × Sk[∇(cid:48)
s ρ])
s∆
(cid:48)−1
+ ∇ · (n × ∇ × Sk[n × ∇(cid:48)
s ρm]) =
s∆
(cid:48)−1
ρm − ∇ · (n × n × ∇ × Sk[∇(cid:48)
s ρ])
s∆
(cid:48)−1
+ ∇ · (n × n × ∇ × Sk[n × ∇(cid:48)
s∆
s ρm]) =
where g3 = ∇ · (n × Hi) and g4 = ∇ · (n × n × Hi).

g3−jω

,

g4−jω

,

(5.16a)

(5.16b)

Both of these integral equations are associated with an identity plus compact operators,

which means the eigenvalues of the system cluster around a non-zero value if the frequency

doesn’t correspond to a resonance frequency for the interior domain. Thus, they have better

spectral properties than the one derived using electric fields. Direct implementation of the

iso-geometric magnetic field integral equation using potential as the unknowns, as in [5], will

not necessarily lead to such well-conditioned system.

5.2.5 Combined Integral Equation

It is apparent that both the integral formulations in (5.13) and (5.16) can be used to solve for

scattering from PECs and don’t involve non-physical solutions at low frequencies. Neither

of them has the low-frequency breakdown as a rigorous Helmholtz decomposition is used

(indirectly on the current). However they still suffer resonance problem at higher frequencies.

To remove it, one can use the approach as in the combined field integral equation (CFIE),

which is to use the linear combination of the two integral equation sets. The combined

version or the scalar CFIE can be written as

sMFIE +

sEFIE,

α
η0

(5.17)

where η0 is impedance in free space and α is the weighting factor.

80

5.3 Method of Moments Solution Procedure

In what follows, we prescribe a method of moments solution procedure that is based on

a smooth representation of the geometry. The motivation behind this representation is two

fold: (i) the preceding formulation makes use of a Laplace-Beltrami operator (and other

operators) that rely on a smooth representation of the geometry, and (ii) our using both

scalar unknowns and observables opens the door to using the same basis set to represent

both the geometry and the unknowns on the geometry. These are elucidated in the next

subsection.

5.3.1

Iso-geometric Analysis on Subdivision Surface

To numerically implement the above integral equations, one has to choose suitable geome-

try representation and corresponding scalar basis functions defined on it. The subdivision

surface, unlike polygon or splines based geometric representations, can lead to global smooth-

ness of at least C1 everywhere for arbitrary shaped surfaces. One of the popular subdivision

schemes is called Loop subdivision [48, 49] and it will be used for the implementation pur-

posed in this work. An equivalent basis function that is (a) smooth and (b) of compact

support will be generated through the process of subdivision as is detailed in [5, 49]. This

leads to the idea of iso-geometric analysis.

Generally speaking, iso-geometric analysis on subdivision surface is chosen to implement

the new formulation due to the considerations in two aspects. Firstly, as shown in [5],

subdivision enhanced Galerkin method has natural choice of scalar basis functions defined

on a globally smooth representations (up to C2 almost everywhere or at least C1 everywhere).

The same basis used for both physical simulation and geometry representation allows the

design-through-analysis in the engineering practice. Secondly, smooth geometry (smooth

normal directions in this case) and smooth basis allow the reduction of the singularities, for

example, by transferring derivatives from the integral kernels to the test functions.

81

5.3.2 Scalar Basis Function and Linear System Setup

A scalar function (including each component of the position vector) defined on the the

subdivision surface can be evaluated at local coordinate (u, v) by the following expansion

Nv(cid:88)

f (r(u, v)) =

aiξi(r(u, v))

(5.18)

i=1

where Nv is the number of vertices in the primal control mesh and i is the global index for
each vertex. Since each basis function ξi(r(u, v)) associated with vertex i has finite (two-
neighborhood) support [5,49], the contribution in last expansion can be reduced into a local

one,

Nval+6(cid:88)

j=1

f (r(u, v)) =

ajξj(r(u, v)),

(5.19)

where Nval is the valence (equal to 6 if there is no irregular vertex) of the irregular vertex
in the current element of the control mesh. The expansion (5.19) is quite similar to two

dimensional scalar finite element methods, where unknown function can be represented by

several basis or interpolatory functions associated with the interested element. The only

difference is that, in subdivision basis, the contribution of each vertex is covering triangles in

its two ring around it [5, 49] (Figure 5.1), rather than 1-ring as in hat functions case, hence

there will be more but smooth basis functions contributing to each element.

The baisis function can be evaluated in the same way as explained in chapter 3. An

example of the basis function associated with a vertex that has a valence of 5 is given in

Figure 5.2.

Using Gauss and singular quadrature rules, one can evaluate the elements of impedance

matrix. Two parts of the discrete system are involved, where the first one is from the integral

operator and the the second one is from the surface Laplacian operator. In the first system,

surface Laplacian is not considered, as its inverse operator is not explicit available for matrix

assembly for generalized surfaces. One can consider the Laplacian as a right preconditioner

for general case, though in the sEFIE case the inverse Laplacian operator can be canceled

82

written as the 2 by 2 block matrix,

Z =

 ,

(5.20)

Figure 5.1: The support of the scalar basis: two ring of a valence-5 vertex

in the ∇s · J term. Here only the system assembly for the sMFIE is given, since other two
formulations can be discretized in the similar way. The first linear system for sMFIE can be

A B

C D

83

(cid:90)

(cid:90)
Aij = Gij −

Bij =

Si
(cid:90)
Cij = −

Si

fi(r)P1[∇(cid:48)
fi(r)P1[n(cid:48) × ∇(cid:48)
(cid:90)
fi(r)P2[∇(cid:48)
fi(r)P2[n(cid:48) × ∇(cid:48)

Si

sfj]dS,

sfj]dS,

sfj]dS,

(5.21a)

(5.21b)

(5.21c)

Figure 5.2: Scalar basis associated with a vertex

Each element of the matrix block corresponds to the following integrals

Dij = −Gij +

(5.21d)
where Gij is the matrix element (with index {i, j}) of the Gram system and P1[j(r(cid:48))] =

G(r, r(cid:48))j(r(cid:48))dS(cid:48)] and P2[j(r(cid:48))] = ∇ · [n × n × ∇ ×(cid:82)

∇ · [n × ∇ ×(cid:82)

G(r, r(cid:48))j(r(cid:48))dS(cid:48)].

sfj]dS,

Si

Sj

Sj

The second system is the inverse of the surface Laplacian. One has to take the inverse

operator through the inverse of the discrete Laplacian. The discrete Laplace-Beltrami oper-

ator can be also set up through a Galerkin scheme, the matrix entry of the Laplacian system

84

can be computed through

(cid:90)

S

Lij =

∇sfi · ∇sfjdS

(5.22)

where the numerical integration for each matrix element could only happen on a compact

region due to the compactness of the basis functions. The nonsingular right hand side vector

can be evaluated similarly using Gauss quadrature as in regular methods of moments.

5.3.3 System Solution

The solution to the system of equations is slightly more challenging than a conventional

system due to the inverse of the surface Laplacian that is involved in the solution procedure.

As will become apparent, direct solutions are difficult. As a result we take recourse to an

iterative procedure. The two steps in the solution procedure are as follows: (1) invert the

discrete surface Laplacian operator that operates on charge data and (2) solve the discrete

pseudo-impedance system that operates on the potential data. The first step, in our work, is

implemented through an inner iterative solver, conjugate gradient based solvers are effective

in solving the surface Laplacian system [115]. Special care has to be taken to carry out

the first step since the surface Laplacian operator is exactly rank-deficient by one. One can

compute the inverse by either reducing the rank of the system accordingly as in [5] or using

more rigorous methods as discussed in [115]. The second step is identical to that used in

regular method of moments implementations, which can be accelerated by fast algorithms

[102].

5.4 Numerical Examples

In this section, the proposed formulation will be studied in several aspects through nu-

merical examples. First, the spectral properties and stabilities will be tested on the unit

spherical surface. Next, validation of the numerical implementation will be given, together

with convergence study of the iterative solver. Finally, more generalized examples are stud-

ied.

85

5.4.1 Spectral Properties

In this test, the analytical properties of system of equations derived earlier are studied

using spherical harmonic functions on a unit sphere. Since only scalar basis functions are

required, the scalar tesseral harmonics are used in a Galerkin framework to derive a modal

representation of the system of equations; as is expected, this is block diagonal matrix, with

each block corresponding to one mode. The procedure is very similar to the one that models

quasi-analytical transient acoustic scattering in [3]. If no surface Laplacian is involved, the

modal system is quite similar to the their vector counterparts in [51]. Here we examine the

properties of the scalar formulations, especially the scalar MFIE, that has not been reported

elsewhere.

First, in the low frequency range where ka < 1, the conditioning for three systems

are studied. The condition number and the eigen-condition number are used to examine

the behavior of these equations across the frequency spectrum, where the eigen-condition

number denotes the ratio between largest and smallest absolute values of eigenvalues.

It

is noted that neither sEFIE nor sMFIE will suffer from the spurious resonance problem

within the interested low frequency range that is below the first resonance frequency of the

corresponding spherical cavity. Only scalar EFIE, with its scaled version (as in [2] or other

loop-tree or loop-star based methods), and MFIE are studied here. Figure 5.3 plots condition

number curves versus the frequency for three systems, and the corresponding eigen-condition

number results are shown in Figure 5.4. From both figures, one can see the scaled scalar

EFIE and scalar MFIE have stable condition numbers over the wide low-frequency band.

In scaled scalar EFIE case, the fundamental reason is that Helmholtz decomposition and

frequency-rescaling resolve the low-frequency breakdown. However, in scalar MFIE case, the

existence of second kind integral operators ensures the good conditioning in the interested

frequency range.

Second, at high frequencies, the spectral behaviors are tested by plotting the eigenvalues

of the interested system. At high frequency, neither sEFIE nor sMFIE is immune from

86

Figure 5.3: Condition number versus frequency

Figure 5.4: Eigen-condition number versus frequency

87

100102104106108frequency (Hz)100105101010151020condition numberscalar EFIEscaled scalar EFIEscalar MFIE100102104106108frequency (Hz)100105101010151020eigen-condition numberScalar EFIEScaled Scalar EFIEScalar MFIEinternal resonance, hence, combined version of scalar formulations has to be used. Figure

5.5 and Figure 5.6, respectively, demonstrate the eigenvalues and singular values of three

systems (sEFIE, sMFIE and the combined formulation), where the frequency is f = 10GHz

and the highest mode is at the degree of n = 210. The singular value distributions are

very similar across the formulations while the performances in eigenvalue spectra are very

different. Since the sMFIE is a second kind integral equation, it can keep the compact

spectrum as the regular MFIE, much better than that of scalar EFIE. Due to the spectral

properties of the operators in sEFIE, the combined formulation is not second kind anymore,

however, as is expected its spectrum is better than that of the scaled sEFIE.

(a)

(b)

Figure 5.5: Eigenvalues of the system of (a) scalar EFIE, (b) scalar MFIE and (c) scalar
CFIE with α = 0.5/η

(c)

88

-2-101234real part-5-4-3-2-101imaginary part-0.500.511.5real part-0.6-0.4-0.200.20.40.6imaginary part-0.500.511.522.5real part-2.5-2-1.5-1-0.500.51imaginary partFigure 5.6: Singular values of the system of the scalar EFIE, the scalar MFIE and the scalar
CFIE with α = 0.5/η

5.4.2 Numerical Test

Next, we test the accuracy of the formulation numerically by computing the radar cross

section of the unit sphere. An ˆx-polarized plane wave propagating in the ˆz direction at

frequency f = 600MHz is used as the incidence field.

In the geometry representation,

2562 control vertices are involved, therefore the total number of degree-of-freedom is 5124

for all the three cases. The linear solver chosen is GMRes (generalized minimal residual

method) with restart set to 30, without using any preconditioner. The radar cross section

data of φ = 0 cut in the far field are computed using the presented formulation and given

in Figure 5.7, and analytical result from Mie series approach are also plotted for reference.

Note that all the error in all three formulations at this frequency give almost the same
relative error in far field values (∼ 2 × 10−3); therefore only one result is actually plotted.

89

02004006008001000singular value index10-310-210-1100101singular valuessEFIEsMFIEsCFIEThe performances of the iterative solver for the three formulations are given in Figure 5.8.

At this specific frequency, the sCFIE shows the better performance than sMFIE simply

because of the oscillating nature of eigen-spectrum of the MFIE formulation even though

the frequency chosen does not correspond to a spurious resonance frequency. As shown in

the the result, the sEFIE, like in regular EFIE case doesn’t convergence fast, as is expected

from the spectrum of the sEFIE operator.

Figure 5.7: Radar cross sections of the sphere

Finally, a more complicated structure, the toy plane model with dimension around 11m×
11m×3m in Figure 5.9, to show the stability of the sMFIE formulation at low frequency range.
The model is presented by the subdivision surface with 11, 994 vertices, that corresponds to

23, 988 unknowns when the integral equation is discretized. The incidence field is ˆx polarized

and propagating along ˆy. The scattering problem is solved at different frequencies, mainly

at low frequency range. The current density distribution at 10MHz is also demonstrated in

90

-150-100-50050100150Angle (deg)-20-100102030Radar Cross Section (dBsm)AnalyticalNumericalFigure 5.8: Residual of GMRES Iterations, 600Mhz

Figure 5.9. Actual changes in the residual of the iterative solver (as mentioned previously,

no preconditioner is used for the GMRes solver) at different frequencies are given in Figure

5.10. Wide frequency stability in the iterative solvers can be observed, thanks to the spectral

properties discussed earlier.

5.5 Conclusion

In this chapter, a new well-conditioned scalar charge based formulation has been pre-

sented to model the EM scattering from simply-connected structures. It is noted that fun-

damental unknowns are chosen similarly as in [91]. The the current representation J and

electrical charge representation ρ are physical in this work, whereas they are not physical

in [91] because the indirect approach is used therein. The formulation features wide-band

stability in terms of spectral properties and conditioning of the resulting system. The ad-

91

020040060080010001200Iterations10-1010-810-610-410-2100ResidualsEFIEsMFIEsCFIEFigure 5.9: A plane model: real part of the magnitude of the surface current density

vantage of this formulation in terms of implementation and potential integration with higher

order CAD descriptions is that it involves scalar basis functions; as shown, this enables easy

implementation within an iso-geometric analysis framework via subdivision surfaces or spline

surfaces. We anticipate that the methods developed herein wherein the description of the

geometry and the physics on the geometry are intimately related would be of substantial use

to the engineering community; this is the direction of our future research. Another impor-

tant and non-trivial extension of the current formulation to multiply connected shapes will

be discussed in future communications.

This chapter, c(cid:13)2018 IEEE, is reprinted, with permission, from Li, J., Fu, X. and Shanker,
B., "Formulation and Iso-Geometric Analysis of Scalar Integral Equations for Electromag-

netic Scattering," Antennas and Propagation, IEEE Transactions on, April 2018.

92

Figure 5.10: Residual of GMRES iterations at different frequencies

93

020040060080010001200Iterations10-1010-810-610-410-2100Residual10Hz100Hz1KHz10KHz100KHz1MHz10MHzCHAPTER 6

DECOUPLED POTENTIAL BASED FORMULATIONS

6.1

Introduction

Surface integral equation methods have been widely used for the analysis of electro-

magnetic (EM) scattering and radiation [101, 120–122]. Commonly used integral equations

for perfectly electrical conductors (PECs) include electric field integral equation (EFIE),

magnetic integral equation (MFIE) and combined field integral equation (CFIE) and their

modified forms [101]. For the problem of electromagnetic scattering from dielectric objects

(transmission problem), the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) [121]

and Müller formulations [122] are most popular ones. Issues such as low-frequency break-

down [87, 104], dense mesh breakdown [123] and topology breakdown [124] have been ob-

served in numerical implementations of method of moments when solving these integral

equations. Much of the breakdown phenomenon arises from either catastrophic cancel-

lation, choosing different types of boundary conditions or badly imposed scalings. Di-

rect consequence of the breakdown is ill-conditioning (hence poor convergence of iterative

solvers) in the resulting linear system or lost accuracy in post-processing. Stabilizing ex-

isting integral equations solvers or designing new stable formulation has been extensively

studied by the computational electromagnetics and applied mathematics communities; see

Refs. [5,97,98,101,105,108,110,112,125–128] and references therein for a complete analysis.

Approaches for stabilizing the EFIE or its related incarnations range from loop-tree/star

decomposition [87, 98, 105] (approximate Helmholtz decomposition) to constrained [129]

and rigorous [5] Helmholtz decompositions. Remedies for dense mesh breakdown includes

Calderon preconditioning [97, 123, 125] and quasi-Helmholtz projector-based methods [110]

. All of the aforementioned methods work directly on the ill-conditioned integral equations.

More recently, there has also been an effort to develop new or modify existing formulations.

94

Augmented EFIE (AEFIE) [126,127] is used to fix the low-frequency breakdown by introduc-

ing auxiliary charge terms and continuity constraints. The current-charge integral equation

(CCIE) [108] is very similar to AEFIE but it can be used to develop a second kind integral

equation for analyzing scattering from dielectric objects. Another example of recent work

in this area is the scalar formulations including generalized Debye sources [91,112] or scalar

charge EFIE and MFIE for simply connected structures [128].

The decoupled potential-based approach [103,130,131] is a very recent effort to solve the

low-frequency breakdown problem. Besides its application to addressing breakdown associ-

ated with the low-frequencies, the scalar and vector potential approach can potentially be

applied to simulations that use vector potential directly [103, 130] as opposed to the elec-

tric/magnetic fields. To date, new boundary conditions on the vector and scalar potentials

have been developed so as to describe scattering from perfect electrically conducting (PEC)

bodies [103,130,131]. In [103], a second kind integral equation was constructed based on the

indirect approach, and the formulation presented is well-conditioned in that it does not suffer

from either the low-frequency or dense mesh breakdown. Additionally, it does not have any

spurious resonance issues or suffer from topology breakdown. The integral equation solved

in [131] is not the second kind, but one that behaves like the AEFIE. An effort to address

well-conditioned equations for dielectric objects using the DPIE framework is more recent,

and focuses on further developing and fleshing out ideas presented in [4, 132]. Analysis of

dielectric objects is more complex simply because more unknowns are involved. Another

challenge is choosing suitable unknowns and observables. These are the challenges that we

will address in this chapter; the formulation presented is well-conditioned, is not susceptible

to non-uniqueness due to resonances or breakdown due to either low-frequencies or dense

meshes. Specifically, in this chapter, we will present:

• Decoupled boundary conditions in terms of scalar and vector potential for transmission

problems,

95

• well-conditioned scalar and vector potential integral equations for EM scattering from

homogeneous dielectric objects,

• reduced decoupled integral equations for PEC problems,

• and, a study of the analytical properties of the resulting system to demonstrate features

of the proposed integral equation framework.

The remainder of the paper is organized as follows. Section 6.2 gives preliminary infor-

mation on the classical boundary value problem for EM scattering by dielectric objects. A

new description for the scattering problem is introduced in Section 6.3 involving decoupled

boundary value problems. Sections 6.4 presents the decoupled boundary conditions. Section

6.5 formulates the scalar and vector potential integral equations for the transmission prob-

lem. In Section 6.6, the analytical properties of the presented integral equations are studied

with the help of asymptotic analysis. Section 6.8 provides some numerical results on the

resulting linear system for a sphere. Finally, conclusions and related remarks are given in

Section 6.9.

6.2 preliminaries

Consider a homogeneous dielectric object occupying a volume Ω2 that is immersed in
a homogeneous background Ω1. Let the surface enclosing the domain Ω2 be denoted by
S and equipped with a normal n that points into Ω1. An electromagnetic plane wave

characterized by(cid:8)Ei(r), Hi(r)(cid:9) is incident on the object. Each domain is characterized by

a set of constitutive parameters, permittivity 
i, permeability µi and wavenumber ki for
i = 1, 2. Henceforth, field quantities in domain Ωi will be denoted using the subscript i. For
the purpose of normalization, the permittivity, permeability, and wavenumber in free space

are denoted by 
0, µ0 and k0. The problem to be solved can be posed as follows: Given
the scatterer configuration as described, find the total field (Et(r) and/or Ht(r)) in each
i (r)).

region Ωi that comprises incident field (Ei

i(r) or Hi

i(r)) and scattered field (Es

i (r) or Hs

96

The solution will be given through an integral equation approach which exploits equivalent

sources on the surface that can be used to express the scattered field in Ωi. In this chapter,
a time-harmonic factor ejωt will be used and suppressed, and spatial dependence on r will

be assumed.

To derive the surface integral equation for the scattering problem, the Stratton-Chu

representation for Es

i can be used. Specifically,

i and Hs
i = −jωµSki
Es
i = −jω
Ski
Hs

[Ji] +

[Mi] +

1

jω


1

jωµ

[∇(cid:48)

∇Ski
∇Ski

s · Ji] − ∇ × Ski
[∇(cid:48)
s · Mi] + ∇ × Ski

[Mi]

[Ji]

(6.1)

(6.2)

In the previous representations, the notation for single layer potential operator is used, which

is defined as

Sk[x](r) =

Gk(r, r(cid:48))x(r(cid:48))dS(cid:48),

(6.3)
with the integral path defined as enclosing the object. Gk(r, r(cid:48)) denotes the Green’s function
for Helmholtz equation in free space (with k being the appropriate value corresponding to

S

(cid:90)

the host medium), and

Gk(r, r(cid:48)) =

e−jk|r−r(cid:48)|
4π|r − r(cid:48)|

(6.4)
The integral operators involve J1 = n × H1 ( J2 = −n × H2) and M1 = E1 × n (M2 =
−E2 × n ) as the sources for the radiation field in the exterior (interior) domain. J1 = −J2
and M1 = −M2 can be considered as the boundary conditions for tangential components
required to formulate classical integral equations as discussed later.

The two types of sources are the equivalent electric and magnetic current densities, respec-

tively. Based on the equivalence theorem, formulations like PMCHWT, Müller or combined

field formulations can be derived that involve constructing integral equations associated with

each domain and then imposing the requisite boundary conditions. The manner in which

boundary conditions are used dictates the eventual formulation and its behavior in different

frequency and discretization regimes. For the classical PMCHWT and Müller formulations,

two unknown equivalent current sources are used, two boundary conditions that relate these

97

across boundaries are chosen, which then produce two different integral equations. This

approach is called direct approach, in contrast to the indirect approach that starts from

the boundary condition and then prepares well-chosen integral representations that usually

involve quantities with nonphysical meaning [111, 133].

Another type of integral equation for dielectric objects is based on current and charge

unknowns [108], which introduces charge density in place of surface divergence of current

densities.

In that case, integrands in those operators include electric current, magnetic

current, electric charge and magnetic charge, the unknown surface sources to be solved. As

there are additional unknowns, one needs additional constraints on the system (in this case,

the total number of constraints is four). To obtain these four equations, additional boundary

conditions are imposed on the normal components of electromagnetic fields. Therefore, the

following set of boundary conditions, will be used to set up the four integral equation,

together with extra continuity and charge neutrality constraints:

n × E1 = n × E2
n × H1 = n × H2

1n · E1 = 
2n · E2
µ1n · H1 = µ2n · H2

(6.5)

In the next Section, we present a modified description of the problem based on the two

commonly used potentials, one scalar and the other vector. The new boundary value problem

will comprise a set of two decoupled boundary value problems.

6.3 Representations of the Decoupled Potential

6.3.1 Scalar and Vector Potentials

The starting point of our discussion is the well-known representation of the electric and

magnetic fields in terms of the vector and scalar potentials, viz.,

Ei = −jωAi − ∇Φi

98

(6.6)

Hi =

∇ × Ai

1
µi

(6.7)

The governing partial differential equation (PDE) for the scalar potential Φi is the scalar
Helmholtz equation

and the PDE for the vector potential Ai is the vector Helmholtz equation

(cid:17)
(cid:16)∇2 + k2
(cid:17)
(cid:16)∇2 + k2

i

Φi(r) = 0,

Ai(r) = 0.

i

(6.8)

(6.9)

In the above equations, we have implicitly assumed that the Lorenz gauge is used. Using the

above expressions, one obtains the boundary conditions in terms of the two potentials. Be-

sides the radiation condition at infinity, the coupled description for the boundary conditions

is as follows:

n × (−jωA1 − ∇Φ1) = n × (−jωA2 − ∇Φ2)

n × ∇ × A1 =

1
µ 1

n × ∇ × A2

1
µ 2


1n · (−jωA1 − ∇Φ1) = 
2n · (−jωA2 − ∇Φ2)

n · ∇ × A1 = n · ∇ × A2

(6.10)

It’s worth noting that the A-Φ representation will lead to the same description of the original

problem thanks to the Lorenz gauge. Though another similar pair of potentials, the electric

vector potential (anti-potential) F and magnetic scalar potential Ψ, are used when deriving

formulations for dielectric objects, they have been used mainly for notational purpose only–

to express integral operators involving magnetic currents and their divergence. Therefore,

while we note the existence of a dual potential representation, neither this dual integral

representation of fields no the use of the dual sources associated with these representations

is used here. Such an alternative representation can be constructed using duality principles

and eventually arrives at equations dual to (6.10). As an aside, one can use the pair F-Ψ

instead of A-Φ to represent electric and magnetic fields. The steps for deriving decoupled

potential integral equations for F-Ψ are then dual to those for A-Φ.

As shown later,

different components (trace information) of A and Φ will be used as unknown sources in the

99

representations for them, and all the EM field quantities can be obtained through (6.6) and

(6.7).

6.3.2 Representation of Scattering Potential

Using the Green’s identity and the governing Helmholtz equation, the scattered scalar po-

tential (denoted with superscript "s") representation could be written as

Φs(r) = −Sk

where the double layer potential operator is defined as
∂G(r, r(cid:48))

Dk[x](r) =

(cid:21)

∂n(cid:48)

(cid:20)∂Φ(r(cid:48))
(cid:90)

∂n(cid:48)

S

(cid:2)Φ(r(cid:48))(cid:3)

+ Dk

x(r(cid:48))dS(cid:48)

(6.11)

(6.12)

(6.13)

In the vector potential case, the dyadic Green’s function can be written as

¯G(r, r(cid:48)) = ¯I

e−jk|r−r(cid:48)|
4π|r − r(cid:48)| = ¯IG(r, r(cid:48))

With the Helmholtz equation and the Green’s identity, As(r) can be represented as [130]

As =Sk[n(cid:48) × ∇(cid:48) × A(r(cid:48))] + ∇ × Sk[n(cid:48) × A(r(cid:48))]

− ∇Sk[n(cid:48) · A(r(cid:48))] − Sk[n(cid:48)∇(cid:48) · A(r(cid:48))]

(6.14)

where there are four types of surface sources, associated with the total vector potential

A = Ai + As, to express the scattered vector potential. For notational simplicity, the

following are used to denote surface sources that will be the unknown quantities in the

integral equations derived later;

a(r) = n × ∇ × A(r)
b(r) = n × n × A(r)
γ(r) = n · A(r)
σ(r) = ∇ · A(r)

100

(6.15a)

(6.15b)

(6.15c)

(6.15d)

where ∇ is the 3D volume gradient operator. Using this notation, the integral representation
of vector potential is rewritten as

As = Sk[a] − ∇ × Sk[n(cid:48) × b] − ∇Sk[γ] − Sk[n(cid:48)σ]

(6.16)

6.4 Decoupled Boundary Conditions

In this section, boundary conditions for the scalar potential and the vector potential

will be derived. The development of these models will imitate success in framework used in

both the Poggio-Miller-Chu-Harrington-Wu-Tsai (PMCHWT) and the Müller formulations

wherein one poses the problems in terms of boundary conditions based on surface sources in

associated operators in (6.11) and (6.14). To develop similar equations for the potentials it

follows that one needs to set up boundary conditions on the scalar potential and its normal

derivatives across the interface, so that one can set up two integral equations with two

unknowns. Likewise, for vector potentials, one needs to impose four boundary conditions

(two vector ones and two scalar ones) to be able to formulate four integral equations to solve

the four unknown sources.

Next, although both [103] and [130] cover the PEC case, we will include a brief review

and discussion for completeness.

6.4.1 Decoupled Boundary Conditions for the PEC case

For the PEC, one should have

n × (−jωA1 − ∇Φ1) = 0
∇ × A1 = 0

n · 1
µ

(6.17a)

(6.17b)

From (6.17a), a decoupled potential description of the boundary condition can be derived as

n × A1 = 0
n × ∇Φ1 = 0

101

(6.18a)

(6.18b)

which is stronger than the original boundary condition.

To satisfy (6.18b), ∇sΦ1 = 0 has to be satisfied, which means the surface gradient of
total scalar potential should vanish. Surface gradient data is not commonly used as a source,

so another condition is used.

Φ1 = V0

(6.19)

where V0 denotes a reference potential value. Both (6.18a) and (6.19) are used in both [103]
and [130]. The difference is that [103] allows an extra set of DoFs to deal with the reference

potential, whereas [130] does not deal with this (setting V1 to zero).

One last but very interesting and important point about the PEC case is that the condi-

tion imposed by (6.17b) can be satisfied if (6.18a) holds. The proof is very straightforward,

if the following manipulation is used:

n · (∇ × A) = −∇ · (n × A)

(6.20)

with ∇ × n = 0 being applied.

6.4.2 Decoupled Boundary Conditions for Dielectric Case

In the PEC case, finding the new boundary conditions involving both vector and scalar

potentials is relatively straightforward. However, for the dielectric case, conditions on nor-

mal component quantities have to be satisfied, together with the requirement on tangential

components of electromagnetic fields. A stronger boundary condition set (involving two po-
tentials) has to be derived from the boundary condition on n · E (rather than n · H). This
anti-duality comes from the asymmetric nature of representations of E and H in (6.6) and

102

(6.7). Therefore, a new boundary condition set derived from (6.10) is as follows [130]:

n × A1 = n × A2
n × ∇Φ1 = n × ∇Φ2

n × ∇ × A1 =

1
µ1

n × ∇ × A2

1
µ2


1n · A1 = 
2n · A2

1n · ∇Φ1 = 
2n · ∇Φ2

(6.21a)

(6.21b)

(6.21c)

(6.21d)

(6.21e)

where the first two conditions are derived from the condition on n × E as in PEC case,
the third one corresponds to the requirement on n × H, and the last two are related with
the conditions of normal components of the electric field. Similar to the PEC case, the

requirement on the normal component of H can be satisfied by considering the fact that
n · (∇ × X) = −∇ · (n × X).

By separating vector potential A and scalar potential Φ, and also introducing a constant

jump term in Φ, one obtains the modified boundary conditions

Φ1 = Φ2 + V1

1n · ∇Φ1 = 
2n · ∇Φ2

(6.22a)

(6.22b)

for the scalar potential Φ, where V1 denotes a reference potential value, and boundary
conditions

n × A1 = n × A2

n × ∇ × A1 =

1
µ1

n × ∇ × A2

1
µ2


1n · A1 = 
2n · A2

(6.23a)

(6.23b)

(6.23c)

for the vector potential.

At this point, only three boundary conditions are recovered for the vector potential,

so another scalar boundary condition is required. The last one is the implicitly imposed

103

condition that ∇· E vanishes; this condition on E provides the necessary information on the
divergence of A and Φ, as

∇ · E = −jω∇ · A − ∇2Φ = −jω∇ · A + k2Φ = 0

(6.24)

from the Lorenz gauge. Due to (6.22a), the additional boundary condition associated with
the vector potential is ∇· A1 = ∇· A2 + V2, where V2 is a reference scalar (potential) that is
constant over the surface of a isolated object. The two surface reference potentials (reference

voltages) can be chosen as unknown quantities to be solved.

For the electromagnetic problem, the charge neutrality constraint has to be imposed as

well. This requirement can be satisfied by setting up stronger conditions on both scalar and
vector potential. The idea is to impose zero-mean constraints on both n · ∇Φ and n · A.
These two requirements come directly from examining the physics of the problem and match

a mathematical approach used in [103]. Using the afore-developed conditions, the boundary

condition set for scalar potential would be written as follows:

Φ1 = Φ2 + V1
(cid:90)

1n · ∇Φ1 = 
2n · ∇Φ2
∂n(cid:48) dS(cid:48) = 0
∂Φ1

S
For vector potential, one would get

n × A1 = n × A2

n × ∇ × A1 =

1
µ1

n × ∇ × A2

1
µ2


1n · A1 = 
2n · A2
∇ · A1 = ∇ · A2 + V2

n · A1dS(cid:48) = 0

(cid:90)

S

(6.25a)

(6.25b)

(6.25c)

(6.26a)

(6.26b)

(6.26c)

(6.26d)

(6.26e)

Due to the gauge freedom, it is known that E-H cannot uniquely determine A-Φ, whereas

the reverse is true. Following the same philosophy (choosing specific boundary value descrip-

tions involving A-Φ to represent the Maxwell’s boundary value problem), one can also use

104

an even stronger constraint by setting the reference terms V1 and V2 to zero.
shown that for this choice the solutions to (6.25) and (6.26) are unique and do not admit

It can be

any gauge nullspaces. In the rest of the paper, the formulation and discussions presented

later are under this assumption. Though solutions satisfy the Maxwell’s equation with the

scalar and potentials being auxiliary quantities, effects of the assumption on the true vector

or scalar potential problem are not known and worth being studied to answer a more funda-

mental question–whether a description using scalar-vector potential A-Φ rather than E-H

is possible [134, 135]. If a potential-only description is possible, then there is another open

problem: can we derive an additional physics-based constraint to augment (6.25) and (6.26)

such that they admit unique solutions without setting the reference terms V1 and V2 to zero.
It’s worth noting that several pairs of the decoupled potential boundary conditions such

as (6.22a) and (6.23a) are much stronger than their electric and magnetic fields counterparts

(6.10). This is a fundamental assumption in all of the existing decoupled potential-based

formulations. If the solution satisfies the decoupled boundary value problem, then the solu-

tion is also the solution to the original Maxwell’s equations. The existence and uniqueness of

the decoupled potential-based boundary value problems are essential to set up the decoupled

potential integral equations; see discussions in [103, 132].

6.5 Formulation of Decoupled Potential Integral Equations

In this section, decoupled potential integral equations will be derived based on the rep-

resentation theorems and corresponding decoupled boundary condition sets.

6.5.1 Scalar Potential Integral Equation (SPIE)

derivative ∂Φ

From the representation theorem for scalar potentials, one can choose Φ and its normal
∂n as the sources and observables to construct the integral equations. Therefore,
one needs the information about the two corresponding incident fields, denoted by Φi(r) and
∂Φi(r)

respectively.

∂n

105

On the surface, two integral equations corresponding to the exterior and interior domain

can be written as follows.

(cid:21)

(cid:20)∂Φ1(r(cid:48))
(cid:21)
∂n(cid:48)
− Dk2

(cid:20) ∂Φ2(r(cid:48))
Φ1 = Φi − Sk1
Φ2 = Sk2

∂n(cid:48)

(cid:2)Φ1(r(cid:48))(cid:3)

+ Dk1
(cid:2)Φ2(r(cid:48))(cid:3)

(6.27a)

(6.27b)

In the integrals of the first (second) IE, the observation point approaches the surface from

the exterior (interior) domain. Usually, Cauchy principal values will be taken when working

R.
with operators with singularity order higher than 1

One needs another two singular integral equations to have the same number of equations

as the unknowns.

∂Φ1
∂n
∂Φ2
∂n

∂Φi
=
∂n
= D(cid:48)
k2

(cid:21)

(cid:20)∂Φ1(r(cid:48))
(cid:21)
∂n(cid:48)
− Nk2

− D(cid:48)
k1

(cid:20)∂Φ2(r(cid:48))

∂n(cid:48)

(cid:2)Φ1(r(cid:48))(cid:3)

+ Nk1
(cid:2)Φ1(r(cid:48))(cid:3)

(6.28a)

(6.28b)

where the normal derivatives of operators Sk and Dk are denoted by D(cid:48)

k and Nk respectively.
By linearly combining the two equations in (6.27) and the two equations in (6.28) and

applying the boundary conditions in (6.25), one obtains the following scalar potential integral

equation (SPIE): α1+α2

2

I + C11
C21



 Φ


1
k0
0

∂Φ
∂n

 =

 Φi


1
k0
0

∂Φi
∂n



(6.29)

C12
I + C22

β1+β2

2

where the scale factor


1
k0
0

on the ∂Φ

∂n is used to get the same dimensionality as in Φ and

,

+ α2 ˜Dk2
C11 = −α1 ˜Dk1
− α2k0
0
Sk1
Sk2
C12 =
α1k0
0

1
C21 = −β1
1
Nk1
Nk2
k0
0
− β2 ˜D(cid:48)
C22 = β1 ˜D(cid:48)
k2
k1


2
β2
2
k0
0

+

.

,

,

(6.30a)

(6.30b)

(6.30c)

(6.30d)

106

In (6.30) and equations later, operators with tildes mean the integral is taken in the Cauchy

principal value sense.

Following a similar manner as in the Müller system of equations,

the constraint β1
1 = β2
2 has to be imposed in order to remove the hyper-singularity in
C21. The impact of choices of other parameters on invertibility and actual conditioning is
worth being studied theoretically and numerically, especially when lossy medium is involved.

With that choice, all the operators are bounded and compact. Since the two operators in
the diagonal in (6.30) are in the form of an identity operator plus compact operators (C11
and C22) when the surface is smooth, the integral equation (6.30) is of the second kind.

6.5.2 Vector Potential Integral Equation (VPIE)

The vector potential integral equation corresponding to the vector potential boundary value

problem can be derived in a manner similar to that used for the scalar potential, but it

involves choosing suitable observables and different scalings.

Since the goal is to construct a well-conditioned formulation, it is very natural to choose

the same set of trace information of the vector potential as the observables. As in the scalar
potential case, incident field information including n×∇× Ai, n× n× Ai, n· Ai and ∇· Ai
must be available.

From (6.16), one can write the representations for the four types of observables in the

following two sets (one for the exterior and the other for the interior) of integral equations.

The first set is obtained by allowing the observation point to approach the surface from the

exterior, yielding the following exterior vector potential integral equation (VPIE) set:

 =





a1

b1

γ1

σ1

ai

bi

γi

σi

 +



Kk1
St
k1
Sr
k1
∇ · Sk1

−Tk1
−K(cid:48)
k1
−M3
k1

0

0
−P2
k1
−D(cid:48)
k1
1Sk1
k2







a1

b1

γ1

σ1

−Q1
k1
−Q2
k1
−Q3
k1
Dk1

(6.31)

The operators T , K and K(cid:48), respectively, denote the electromagnetic hyper-singular operator,
the MFIE operator and its adjoint operator. D and D(cid:48) are scalar double layer potential

107

operator and its adjoint, respectively. Their explicit definitions and properties are given in

the Appendix for completeness. The other operators in (6.31) are less well-known and are

defined as follows:

k[b] = n · ∇ × Sk[n(cid:48) × b]
M3
P2
k [γ] = n × n × ∇Sk[γ]
k[σ] = n × ∇ × Sk[n(cid:48)σ]
Q1
k[σ] = n × n × Sk[n(cid:48)σ]
Q2
k[σ] = n · Sk[n(cid:48)σ]
Q3

(6.32)

Among these, all the operators in the diagonal can be written in the form of identity plus a
compact operator (as shown in the appendix). Tk is a hyper-singular operator, and has the
same properties as that of the EFIE. Each operator in the skew-diagonal is bounded but not

compact.

A second, interior VPIE set is obtained by allowing the observation point to approach

the surface from the interior. Using the same operator definitions as in (6.31) and (6.32), but

with subscripts k2 to denote use of the interior region wavenumber in the Green’s function,
the same notation can be used to express the observables as

 = −





a2

b2

γ2

σ2

Kk2
St
k2
Sr
k2
∇ · Sk2

−Tk2
−K(cid:48)
k2
−M3
k2

0

0
−P2
k2
−D(cid:48)
k2
2Sk2
k2

−Q1
k2
−Q2
k2
−Q3
k2
Dk2







a2

b2

γ2

σ2

(6.33)

For convenience, Z1 and Z2 are used to denote the operator matrices in (6.31) and (6.33)
respectively.

In order to simplify application of the boundary conditions and further improve the

conditioning of the system, the following scaled quantities (as in scalar potential case) are

108

used. That is,

√

√

√

µ0
µ

a(cid:48) =
µ0
a =
µ
√
b(cid:48) = −jω

0b = −jω
γ(cid:48) = − jω
√
γ = − jω
√

0

0
σ(cid:48) =
∇ · A

σ =

1√
µ0

1√
µ0

n × ∇ × A


0n × n × A

n · A

(6.34a)

(6.34b)

(6.34c)

(6.34d)

A careful dimensional analysis shows that all the above quantities (a(cid:48), b(cid:48), γ(cid:48), σ(cid:48)) have the
same units. It is not difficult to realize that the first three terms represent scaled electric

current density, scaled tangential electric field and scaled electric charge density, whereas

the last one is scaled potential. Hence the unknown quantities presented in this chapter

can be related to physical quantities, which is an advantage of the direct approach over the

indirect approach in setting up integral equations [101, 111]. In light of this interpretation,

the boundary conditions are tantamount to all the four fields being continuous across the

interface. It is also worth noting that the first two unknowns are part of the scaled electric

current and scaled tangential electric field. The presence of the other two quantities (one

of them is a scaled charge density term), together with the suitable linear combination of

equations for both domains, leads to stability properties as discussed later.

To reflect the changes in the scaling while keeping the identity operator unchanged, one

can define the following block diagonal left and right preconditioners. Since later analysis

will involve objects with electrical size kd, all of these scaling factors are written in terms of

ki, 
i and/or µi, with i = 0, 1, 2 denoting free space, exterior medium and interior medium,
respectively.


0

,− jk0
i
√
µ0
√

,−
0

µ0
jk0
i

1√
µ0

)

√

µ0).

,

,

(6.35)

(6.36)

,− jk0√
µ0
√
µ0
jk0

,−

109

Pl,i = diag(

and

√

µ0
µi

Pr,i = diag(

µi√
µ0

where Z(cid:48)
1 denotes the new operator matrix after introducing Cauchy principal value integrals
for the operators in the diagonal of Pl,1Z1Pr,1. Similarly, the scaled form of the interior VPIE
(6.33) becomes



(6.37)

ai

bi

γi

σi

(cid:18)I

2

(cid:19)

− Z(cid:48)

1



1

1

a(cid:48)
b(cid:48)
γ(cid:48)
σ(cid:48)

1

1

 = Pl,1


 = 0

a(cid:48)
b(cid:48)
γ(cid:48)
σ(cid:48)

2

2

2

2

(cid:18)I

2

(cid:19)

+ Z(cid:48)

2

Thus the scaled form of the exterior VPIE (6.31) becomes

The explicit form for Z(cid:48)



i (i = 1, 2) in (6.37) and (6.38) is
˜Kki

Tki

1

jk0µr
− ˜K(cid:48)
ki
−
rM3
ki

−jk0µrSt
ki
− jk0
Sr
cr
ki
µr∇ · Ski

0
− 1
P2

r
ki
− ˜D(cid:48)
ki
− k2
jk0
r

i

− 1
Q1
µr
ki
jk0Q2
ki
jk0
rQ3
ki
˜Dki

Ski
where the relative light speed cr is defined as cr = 1√
VPIEs for exterior and interior domains, one obtains the final VPIE,


rµr

0

. By linearly combining the two

(cid:18) Q1 + Q2

2

(cid:19)

− Q1Z(cid:48)

1 + Q2Z(cid:48)

2

 = Q1Pl,1





1

1

a(cid:48)
b(cid:48)
γ(cid:48)
σ(cid:48)

1

1

where the linear factors are defined by

Q1 = diag

Q2 = diag

(cid:18) 1
(cid:18) 1

µr2

µr1

,

,

1

r2
1

r1

, 
r2, µr2

, 
r1, µr1

(cid:19)
(cid:19)

110




ai

bi

γi

σi

(6.38)

(6.39)

(6.40)

(6.41a)

(6.41b)

The choice for linear factors is determined by requiring the cancellation of singularities be-
R, that is, to cancel the singularity in those non-compact operators as done in Müller
formulation [122]. This particular linear combination not only has the advantage of eliminat-

yond 1

ing non-compact operators for increased numerical stability, but of also allowing an expanded

choice of basis and testing functions with considerably relaxed continuity/conformability re-

quirements.

6.6 Analytical Properties

In this section, we investigate the analytical properties of the presented SPIE and VPIE

by studying scattering from a sphere of radius a. The resulting linear system from both

integrals can be analytically evaluated if scalar and vector spherical harmonics are used to

represent the unknowns associated in each integral equation. Due to the orthogonality be-

tween basis functions, a block diagonal system will be generated, making it possible to study

the spectral properties of the discrete system. This approach has been used widely as an

efficient tool to understand/capture some of the essential signatures of integral operators

associated with three-dimensional time harmonic or time-dependent acoustic or electromag-

netic problems [2, 3, 51, 91, 103].

Unknown scalar (u) and vector (u) quantities are represented using scalar and vector

spherical harmonics [2] such that

where

Y m
n (θ, Φ) =

(n − m)!
(n + m)!

n (cos θ)ejmφ
P m

Ψm

n (θ, φ) =

∇tY m

n (θ, φ) =

Φm

n (θ, φ) =

ˆr × ∇tY m

n (θ, φ) =

˜Ψm

n (θ, φ)

(cid:112)n(n + 1)
(cid:112)n(n + 1)

n (θ, φ)

˜Φm

(cid:88)

unmY m
n

u1
nmΨm

n + u2

nmΦm
n

u =

u =

2n + 1

(cid:88)
(cid:115)
r(cid:112)n(n + 1)
r(cid:112)n(n + 1)

4π

111

(6.42)

(6.43)

(6.44)

(6.45)

(6.46)

n (θ, Φ) = −ˆr × Φm

monics satisfy the relations Ψm

where r = a and ˆr is the unit vector pointing in the radial direction. The vector spherical har-
n (θ, Φ). The
linear system elements are evaluated as integrals of several types: (1) scalar-scalar integral <
Y m(cid:48)
n(cid:48) ,OP[Y m
n(cid:48))] >, (3) vector-scalar inte-
n(cid:48) (Ψm(cid:48)
gral < Ψm(cid:48)
n(cid:48))] >,
where OP denotes an appropriate operator.

n ] >, (2) scalar-vector integral < Y m(cid:48)
n(cid:48) (Ψm(cid:48)

n (Ψm
n ] > and (4) vector-vector integral < Ψm(cid:48)

n(cid:48) ,OP[Ψm

n(cid:48) ),OP[Ψm

n (Ψm

n (θ, Φ) and Φm

n (θ, Φ) = ˆr × Ψm

n(cid:48) ),OP[Y m

6.6.1 Stability Properties of SPIE

The scalar potential integral equation is well-conditioned and doesn’t suffer from the dense

mesh breakdown problem. The formulation for the scalar potential integral equation is

similar to that for the transmission problem in acoustics [111], hence it is immune to spurious

resonances. The detailed analysis will be not be presented; the analysis procedure is very

similar to that for VPIE, which will be given next.

6.6.2 Stability Properties of VPIE

Except for the identity operators along the diagonal, the vector potential integral equation

involves only compact operators. At high frequency, one should avoid allowing the system
element to grow as ka. As ka → ∞, the asymptotic behavior for the system elements in Z(cid:48)
is of the following form:

i

(6.47)



a2

√

µr
r
µr
1
2

r√
µr
r

0

1
2
µr√
µr
r

O((ka)−1)

µr√
µr
r


r

0
√
1
µr
r
1
2
µr
r

r

√

√
1
µr

rµr
O((ka)−2)

µr

√
1

rµr
1
2



It should be noted that the asymptotic analysis focuses on ka, assuming the mesh resolution

is proportional to the wavenumber set by the spatial Nyquist sampling rate (ka and mesh

density approaching infinity at the same time is not practical in numerical analysis). As seen

112

from (6.47), in each row, operators in the off-diagonal have different scalings in terms of

material properties, therefore it is not possible to have all the off-diagonal elements approach

zero at the same time. Actually, numerical examples show that at high frequencies, due to

the oscillatory nature of each operator, the spectral properties (such as eigenvalues and

eigen-radius) for the system are also oscillatory.

At low frequencies, when ka → 0, one should avoid terms such as O( 1
in Z(cid:48)

ka) that would
i has the same
problem as regular EFIE. The situation can be easily fixed thanks to the fact that at very

lead to catastrophic cancellations.

It is easy to find that

Tki

jkoµr

1

low frequencies, the electric and magnetic fields become decoupled. The frequency scaling

in (6.35) and (6.36) should be removed be removed at low frequencies (i.e., sub-domains

become electrically small), with the resulting system matrix modified to



˜Kki
−jµrSt
ki
− j
Sr
cr
ki
µr∇ · Ski

Tki
1
jµr
− ˜K(cid:48)
ki
−
rM3
ki

0

0
P2
− 1

r
ki
− ˜D(cid:48)
ki
− k2
Ski
i
j
r

− 1
Q1
µr
ki
jQ2
ki
j
rQ3
ki
˜Dki

Low-frequency stability properties can be studied by looking at the asymptotic behav-
ior as ka → 0 for fixed resolution (indicated by fixing the highest mode degree n). The
asymptotic scaling of Z(cid:48)



O(1)

i in (6.48) behaves like
O(1)
1
0
µr
O(1)
O(1)

rO(1) O(1)
0

µrO(1)

1

r

a2

0

0

0

O(1)

1
µr

0

0
O(1)

where each of the O(1) terms is only in terms of spatial resolution parameter n, unit imag-
inary number j and possibly a, but shows no dependence on constitutive parameters. Ap-

parently, after the frequency scaling, all the terms are bounded and no serious cancellation

will occur. Furthermore, by choosing the linear combination factors as in (6.41), one obtains

113





(6.48)

(6.49)



a2

O(1/n)

0

0

1
µr

O(n)
1
0
µr
O(1/n)
O(1)

rO(1) O(1/n)

1

r

0

0

O(1)

µrO(1)

0

0

O(1/n)



vanishing off-diagonal elements in (6.40). It is worth noting that choosing correct boundary
conditions for b(cid:48) and σ(cid:48) is essential to achieving this goal, because it leads to the same
scaling factor in front of the second and fourth operators of the first row in (6.49). At low

frequencies, the system is diagonally dominant.

Another important issue in EFIE or EFIE-like formulations is dense mesh breakdown

when element size h is close to zero or the spatial resolution (mode degree) is very high. For

fixed ka, the dependence of the system elements on spatial resolution n (proportional to 1
h

in piecewise discretization) can be derived as

(6.50)

Again, all the O(·) terms are independent of the material constitutive coefficients that ex-
plicitly appear in front of the operators. From the asymptotic result, one sees that the only

term that can possibly cause density breakdown is the hyper-singular operator. However,

after combining two equations for both interior and exterior domain with the help of (6.41),
all the O(1) and O(n) terms are canceled exactly at low frequencies, owing to the fact that
the resulting operators in the off-diagonal are compact.

6.7 Perfectly Electrical Conductor Case

The presented formulation can be reduced to a simpler one for scattering analysis of PECs.

In this section, several integral equations for PECs will be given briefly, and comparison will

be made to results presented in [103, 130, 131].

Using the condition Φ = 0, one can reduce the SPIE to an equation that has only Φ (and

V1 if necessary) as the unknown quantity,

Sk1

[

∂Φ1(r(cid:48))
∂n(cid:48)

] = Φi

114

(6.51)

where charge neutrality (cid:82) ∂Φ1(r(cid:48))

∂n(cid:48) dS(cid:48) = 0 is imposed. The resulting formulation will suf-
fer from spurious resonances; by combining the equation with its normal derivative (as in

the Burton-Miller [136] approach), one can make the SPIE for PECs immune to spurious

resonances. The indirect approach used in [101, 103] can be also used.

By setting n × n × A and ∇ · A to zero, the VPIE in (6.31) is reduced to



Kk1
St
k1
Sr
k1
∇ · Sk1

0

ai

γi

γ1

a1

bi


 =

 +
Kk1
 −
a1

σi

0

γ1



a1

γ1

0
−D(cid:48)
k1

Sr
k1

a1



γ1

0
−P2
k1
−D(cid:48)
k1
1Sk1
k2


ai
 =

γi



with two unknowns (a1 and γ1) governed by four integral equations. There are several
options for choosing two equations from them. One choice is to formulate a second kind

integral equation,

(6.52)

(6.53)

(6.54)

which suffers neither the low-frequency breakdown nor the dense mesh breakdown, but has

the spurious resonance problems at high frequencies.

Another choice is to choose the second and fourth equations,

 St

k1

∇ · Sk1

−



a1

 =

γ1

−P2
k1
1Sk1
k2

bi



σi
1Sk1

This formulation also suffers from the resonance problem. k2

will vanish as the frequency

approaches zero , and the system leads to a saddle point problem that needs special care [131].

At high frequency, a frequency scaling has to be made to fix the ka dependence of the same

operator.

A linear combination of the above two sets of VPIEs (6.53) and (6.54) should be helpful

in avoiding resonances. Besides those two formulations, other options of choosing two dif-

ferent equations from the four in (6.52) are possible, and they might have different spectral

properties and different immunities against spurious resonances.

115

6.8 Numerical Results

In this section, numerical examples are given to demonstrate the stability properties of

the presented integral equations. For simplicity, the scalar potential Φ is assumed to be zero,

and only the vector potential integral equation is left to be solved. Thus the DPIE solution

only involves the VPIE. This is done for two reasons; (a) the operators involved in the SPIE

are well studied in the context of acoustics [101, 111] and (b) since the operators involved

in the VPIE are are more complex, this approach highlights their behavior. In the general

case, the SPIE and VPIE must be solved simultaneously to recover both the electric and

magnetic fields.

In the first test, the linear system for a spherical object is studied for the condition

numbers and eigenvalue distributions. The transmission problem has the following setup:


2 = 2
1 = 2
0, and µ2 = µ1 = µ0. Modes up to 30th degree are used and the system size
is about 1800 × 1800. The condition numbers of the new formulation for a 1m sphere are
very low across the low-frequency band, and remain near 1.25 for frequencies from 1 Hz to

107 Hz. The stable condition numbers are due to the second kind nature of the integral

equation.

Figs. 6.1 and 6.2 show the eigenvalues of the system at 1.0 Hz and 107 Hz, respectively.

All of the eigenvalues are clustered around 0.5 in the complex plane, both figures illustrating

the nice spectral properties. As a result, an iterative solver will be very efficient for systems

such as this one.

The following test is to study the behavior of conditioning at high frequencies where

the spatial resolution must grow proportional to the frequency.

In the implementation,

the highest degree of the basis functions is set to [2ka] + 1. The condition number of

the VPIE versus frequency is demonstrated in Fig. 6.3. For comparison, the condition

number for the case of the Müller formulation is given in Fig. 6.4. For reference, the dashed

curves in both figures indicate the trend that a linearly varying condition number would

follow.

It’s observed that both formulations lead to an oscillatory increase in condition

116

Figure 6.1: Real and imaginary parts of eigenvalues at 1 Hz

Figure 6.2: Real and imaginary parts of eigenvalues at 1 × 107 Hz

number with frequency. Though the high frequency behavior may not be as ideal as at low

frequency, a growing condition number proportional to the electrical size does not necessarily

lead to the same situation for the iteration count. As in other extant approaches, the

convergence of iterative solvers in the high frequency regime can be accelerated using effective

preconditioning techniques [137–139]. This could be a future topic worth more study and

discussions.

117

0.440.460.480.50.520.540.56real part-505imaginary part10-150.440.460.480.50.520.540.56real part-0.015-0.01-0.00500.0050.01imaginary partFigure 6.3: Condition number of VPIE formulation versus frequency

Figure 6.4: Condition number of Müller formulation versus frequency

In order to show the validity of the formulation, comparison is made between the solution

of the VPIE and that using the Mie series approach. Figs. 6.5 and Fig. 6.6, respectively,
30 in the tangential components of
the electric field (twisted magnetic current density). The error between the Mie series and

give the real part of the coefficients of mode Ψ1

30 and Φ1

the DPIE solution is close to machine precision, thanks to the fact that the basis functions

used are eigenfunctions of the vector Laplace-Beltrami operator. From each of the plots, one

118

1071081091010freq(Hz)100102104106condcond of DPIEreference1071081091010freq(Hz)100101102103104105condcond of MullerreferenceFigure 6.5: Real part of the coefficients for Ψ1
electric field

30 mode in the tangential component of the

Figure 6.6: Real part of the coefficients for Φ1
electric field

30 mode in the tangential component of the

observes that the frequency response of the transmission problem is very oscillatory. As the

frequency decreases, tending to the static limit, the response can still be recovered accurately

by the new formulation.

119

2468freq (Hz)109-1-0.500.511.5CoefficientMieDPIE2468freq (Hz)109-0.3-0.2-0.100.10.20.3CoefficientMieDPIE6.9 Conclusion and Future Work

Decoupled potential integral equations for electromagnetic scattering from homogeneous

dielectric objects have been proposed. The resulting formulations are well-conditioned second

kind integral equations without low frequency or dense mesh breakdown problems. When

reducing the formulation to solve PEC problems, several options are available. Observables

and integral equations must be carefully chosen from (6.52) to obtain formulations avoiding

resonance, low-frequency breakdown or saddle point phenomena.

When setting the scalar potential Φ to zero, the vector potential boundary value problem

becomes an exact (scaled) electric field-based description of the original Maxwell’s transmis-

sion problem. Interestingly, this special case of our formulation is also a direct approach and

is dual to that presented by Vico et al. [132] almost at the same time as this manuscript was

submitted. Their work starts from an indirect approach with rigorous mathematical proof

linking the solution to the resulting integral equation with that of the original transmission

problem. With slight changes, the two formulations can be considered as adjoint to one

another. Therefore, [132] provides mathematical insights on interpreting the uniqueness of

the new vector potential integral equation. From the physical point of view, it’s worthing

asking the question as to whether a constraint on the potential exists that guarantees the

uniqueness of the solutions to the potential-based integral equation. In this chapter, two

scalar potential jumps (V1 and V2) in the decoupled boundary conditions (6.25) and (6.26)
are set to zero.

In terms of using the new set of unknowns (two tangential vectors and two scalars), the

VPIE in the new formulation is also similar to the current-charge integral equation. The

difference lies in that, in our formulation, (1) no continuity constraint is needed and (2)

one charge term and one scalar potential term (rather than two charge terms) are used in

addition to the two tangential terms.

In addition, the SPIE involves another two scalar

unknowns (one is the potential and other one is its normal derivative), again very similar to

the acoustics case.

120

It should be pointed out that both the number of integral operators to be discretized and

the number of unknowns involved in the DPIE equations exceed that of conventional integral

equations. However, the stability properties of the new integral equations and the availability

of fast methods (e.g.

fast multipole, adaptive integral, and algebraic-based compression

methods) make the investigation of these new formulations attractive. Discretization issues,

numerical implementations, and performance, especially at high frequencies, will be studied

and presented in an upcoming communication.

121

APPENDICES

122

APPENDIX A

SPHERICAL EXPANSIONS

A.1 Vector Spherical Harmonics (VSH) and Wave Functions (VSWF)

The two vector spherical wave functions used in conventional expansion of dyadic Green’s

function are defined as follows.

M

(i)

nm(k, ¯r) = −z

(i)
n (kr)Φnm(θ, φ)

(A.1)

(A.2)

N

(i)
nm(k, ¯r) =

[krz

Ψnm(θ, φ)

(i)

n (kr)](cid:48)
kr

(cid:112)n(n + 1)

kr

(i)
n (kr)Ym
z

n (θ, φ)

where Ym

n (θ, φ) = Y m

n (θ, φ)ˆr.

+

In order to derive the spherical harmonics expansion of the tangential trace of the dyadic

Green’s function, the mapping relations between VSH and VSWF should be used.

ˆr × M

(i)
nm(k, ¯r) = z

(i)
n (kr)Ψnm(θ, φ)

ˆr × ˆr × M

(i)
nm(k, ¯r) = z

(i)
n (kr)Φnm(θ, φ)

ˆr × N

(i)
nm(k, ¯r) =

[krz

(i)

n (kr)](cid:48)
kr

Φnm(θ, φ)

ˆr × ˆr × N

(i)

nm(k, ¯r) = −[krz

(i)

n (kr)](cid:48)
kr

Ψnm(θ, φ)

123

(A.3)

(A.4)

(A.5)

(A.6)

A.2 Spherical expansion of the time domain magnetic dyadic Green’s

function

Frequency-domain dyadic Green’s function of magnetic field is the solution to the wave

equation.

∇ × ∇ × ˜Gm0 − k2 ˜Gm0 = ∇ × [˜I δ(¯r − ¯r(cid:48))],

(A.7)

The dyadic Green’s function can be written in terms of the scalar Green’s function in

free space.

˜Gm0(¯r, ¯r(cid:48), ω) = ∇ × [˜I Go(¯r, ¯r(cid:48))] = ∇G0(¯r, ¯r(cid:48)) × ˜I

(A.8)

As in electric field case, the dyadic Green’s function can be expanded using a series of

vector spherical wave functions.

˜Gm0(¯r, ¯r(cid:48), ω) = jk2(cid:88)

(cid:104)

M

(4)
nm(k, ¯r)N

n,m

+ N

(4)
nm(k, ¯r)M

(1)∗
nm (k, ¯r(cid:48))
(cid:105)
(1)∗
nm (k, ¯r(cid:48))

(A.9)

By taking the tangential trace of this dyadic Green’s function, one can get a modified

Green’s function expanded purely by vector spherical harmonics, separating the explicit

radial dependence from non-radial dependence.

m0(¯r, ¯r(cid:48), ω) = ˆr × ˜Gm0(¯r, ¯r(cid:48), ω) × ˆr(cid:48) × ˆr(cid:48)
= jk2(cid:88)
˜Gtt

(4)

(cid:104)(cid:2)krz
n (kr)(cid:3)(cid:48)
(cid:2)kr(cid:48)z

kr

(1)∗
n
kr(cid:48)

n,m

− z

(4)
n (kr)

(1)∗
z
n

(kr(cid:48))(cid:3)(cid:48)

(kr(cid:48))Φnm(ˆr)Φnm(ˆr(cid:48))
(cid:105)
Ψnm(ˆr)Ψnm(ˆr(cid:48))

(A.10)

The reason that ˆr× instead of ˆr × ˆr× is chosen in the left side is due to the fact that there
is already one tangential operator ˆr× in the K operator. Similar as in electric field case, the
modified dyadic Green’s function has the following properties.

m0(¯r, ¯r(cid:48), ω) · X = − ˜Gm0(¯r, ¯r(cid:48), ω) · X
˜Gtt

(A.11a)

124

X · ˆr × ˜Gtt

m0(¯r, ¯r(cid:48), ω) = −X · ˆr × ˜Gm0(¯r, ¯r(cid:48), ω)

(A.11b)

In order to get the time-domain expansion, inverse Fourier transform of the above equa-

n,m

tion gives the following formula.
(cid:104)F−1(cid:8)jk
(cid:88)
m0(¯r, ¯r(cid:48), t) =
˜Gtt
− F−1(cid:8)jk
F−1(cid:8)jk
F−1(cid:8)jk

∂
r(cid:48)∂r(cid:48)

∂
r∂r

(1)

(4)

(4)

(1)

∂
r(cid:48)∂r(cid:48)

∂
r∂r

n (kr(cid:48))(cid:9)Φnm(ˆr)Φnm(ˆr(cid:48))
n (kr)(cid:3)z
(cid:2)rz
(cid:105)
(cid:2)r(cid:48)z
n (kr(cid:48))(cid:3)z
n (kr)(cid:9)Ψnm(ˆr)Ψnm(ˆr(cid:48))
n (r, r(cid:48), t)(cid:3)
(cid:2)r(cid:48)K
n (kr)(cid:9) = −∂r(cid:48)
(cid:2)r(cid:48)z
n (kr(cid:48))(cid:3)z
n (r, r(cid:48), t)(cid:3)
(cid:2)rK
n (kr(cid:48))(cid:9) = −∂r
n (kr)(cid:3)z
(cid:2)rz

r(cid:48)

(4)

(1)

(1)

(4)

(0)

(0)

r

(A.12)

(A.13)

(A.14)

Due to the linearity of (inverse) Fourier transform, one can get the following results.

By using above two kernels in (A.12), one can get the spherical expansion for the time-

domain dyadic Green’s function as in (2.26).

125

APPENDIX B

BOUNDARY INTEGRAL OPERATORS

The following are integral operators commonly used in integral equations for Helmholtz and

Maxwell’s equations. Limiting cases for some of them are also given to help the analysis of

properties of the presented integral equation based formulation.

S[σ] =

(cid:90)
(cid:90) ∂G(r, r(cid:48))
(cid:90) ∂G(r, r(cid:48))

∂n(cid:48)

D[σ] =

D(cid:48)[σ] =

N [σ] =

∂
∂n
K[a] = n × ∇ ×

K(cid:48)[a] = n × n × ∇ ×

σ + ˜D[σ]

σ + ˜D(cid:48)[σ]

=== ±1
2
=== ∓1
2
σ(r(cid:48))dS(cid:48)

G(r, r(cid:48))σ(r(cid:48))dS(cid:48)
σ(r(cid:48))dS(cid:48) r→r±
σ(r(cid:48))dS(cid:48) r→r±
(cid:90) ∂G(r, r)
∂n
D[σ] =
G(r, r(cid:48))a(r(cid:48))dS(cid:48) r→r±
=== ±1
2
G(r, r(cid:48))n × a(r(cid:48))dS(cid:48) r→r±

∂
∂n

∂n(cid:48)

(cid:90)
(cid:90)

(cid:90)

=== ∓1
2
G(r, r(cid:48))n × a(r(cid:48))dS(cid:48)

a + ˜K[a]

a + ˜K(cid:48)[a]

(B.1)

(B.2)

(B.3)

(B.4)

(B.5)

(B.6)

T [a] = n × ∇ × ∇ ×

(B.7)
In the above, r → r± means that the position vector is approaching the surface from the
exterior (+) or interior (−) domain. T and N are hypersingular and unbounded integral
operators, both of which are self-adjoint operators. S, ˜K, ˜K(cid:48), ˜D and ˜D(cid:48) are compact (also
bounded) operators, with the adjoint operators of ˜K and ˜D being ˜K(cid:48) and ˜D(cid:48) respectively [101].
It is straightforward that D, D(cid:48), K and K(cid:48) can be used to construct integral equations of the
second type.

Also we have the following convention for denoting different traces of one operator: Sn =

n × S , St = n × n × S and Sr = n · S.

126

BIBLIOGRAPHY

127

BIBLIOGRAPHY

[1]

J. Li and B. Shanker, “Time-dependent lorentz-mie-debye formulation for electromag-
netic scattering from dielectric spheres,” Progress In Electromagnetics Research, vol.
154, pp. 195–208, 2015.

[2] ——, “Time-dependent debye-mie series solutions for electromagnetic scattering,”
IEEE Transactions on Antennas and Propagation, vol. 63, no. 8, pp. 3644–3653, Aug
2015.

[3]

[4]

[5]

[6]

[7]

J. Li, D. Dault, and B. Shanker, “A quasianalytical time domain solution for scattering
from a homogeneous sphere,” The Journal of the Acoustical Society of America, vol.
135, no. 4, pp. 1676–1685, 2014.

J. Li and B. Shanker, “Isogeometric analysis of em scattering on multiply-connected
subdivision surfaces,” in 2017 IEEE International Symposium on Antennas and Prop-
agation USNC/URSI National Radio Science Meeting, July 2017, pp. 1557–1558.

J. Li, D. Dault, B. Liu, Y. Tong, and B. Shanker, “Subdivision based Isogeometric
Analysis Technique for Electric Field Integral Equations for Simply Connected Struc-
tures,” Journal of Computational Physics, vol. 319, pp. 145–162, 2016.

J. A. Stratton, Electromagnetic theory. John Wiley and Sons, 1941.

A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering. Prentice
Hall, Dec 1990.

[8] M. Kerker, The scattering of light and other electromagnetic radiation. Academic

Press, 1 1969.

[9] W. Hergert and T. Wriedt, Eds., The Mie Theory: Basics and Applications (Springer

Series in Optical Sciences), 2012th ed. Springer, 7 2012.

[10] M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission

of Light by Small Particles, 1st ed. Cambridge University Press, 6 2002.

[11] O. Gulder, G. Smallwood, R. Wong, D. Snelling, R. Smith, B. Deschamps, and J.-C.
Sautet, “Flame front surface characteristics in turbulent premixed propane/air com-
bustion,” Combustion and Flame, vol. 120, no. 4, pp. 407 – 416, 2000.

[12] J. Hansen, Spherical Near-Field Antenna Measurement. London: Peter Peregrinus,

1988.

[13] J. A. Schuller and M. L. Brongersma, “General properties of dielectric optical anten-

nas,” Opt. Express, vol. 17, no. 26, pp. 24 084–24 095, Dec 2009.

[14] D. Deirmendjian, “Scattering and polarization properties of water clouds and hazes in

the visible and infrared,” Appl. Opt., vol. 3, no. 2, pp. 187–196, Feb 1964.

128

[15] M. V. Rybin, A. B. Khanikaev, M. Inoue, K. B. Samusev, M. J. Steel, G. Yushin,
and M. F. Limonov, “Fano resonance between mie and bragg scattering in photonic
crystals,” Phys. Rev. Lett., vol. 103, p. 023901, Jul 2009.

[16] P. Bassan, A. Kohler, H. Martens, J. Lee, H. J. Byrne, P. Dumas, E. Gazi, M. Brown,
N. Clarke, and P. Gardner, “Resonant mie scattering (rmies) correction of infrared
spectra from highly scattering biological samples,” Analyst, vol. 135, pp. 268–277,
2010.

[17] O. M. Buyukdura and S. S. Koc, “Two alternative expressions for the spherical wave
expansion of the time domain scalar free-space Green’s function and an application:
Scattering by a soft sphere,” Journal of the Acoustical Society of America, vol. 101,
no. 1, pp. 87–91, 1997.

[18] B. Liu, A. Yamilov, Y. Ling, J. Y. Xu, and H. Cao, “Dynamic nonlinear effect on lasing

in a random medium,” Phys. Rev. Lett., vol. 91, p. 063903, Aug 2003.

[19] P. K. Jain, K. S. Lee, I. H. El-Sayed, and M. A. El-Sayed, “Calculated absorption
and scattering properties of gold nanoparticles of different size, shape, and composi-
tion: Applications in biological imaging and biomedicine,” The Journal of Physical
Chemistry B, vol. 110, no. 14, pp. 7238–7248, 2006, pMID: 16599493.

[20] K. Leung, K. Luk, K. Lai, and D. Lin, “Theory and experiment of a coaxial probe fed
hemispherical dielectric resonator antenna,” Antennas and Propagation, IEEE Trans-
actions on, vol. 41, no. 10, pp. 1390–1398, Oct 1993.

[21] A. Shlivinski, E. Heyman, and A. J. Devaney, “Time domain radiation by scalar
sources: Plane wave to multipole transform,” Journal of Mathematical Physics, vol. 42,
no. 12, p. 5915, Dec. 2001.

[22] E. Heyman, “Time-dependent plane-wave spectrum representations for radiation from
volume source distributions,” Journal of Mathematical Physics, vol. 37, no. 2, p. 658,
1996.

[23] K. D. Granzow, “Multipole Theory in the Time Domain,” Journal of Mathematical

Physics, vol. 7, no. 4, p. 634, Dec. 1966.

[24] W. C. Davidon, “Time-dependent multipole analysis,” Journal of Physics A: Mathe-

matical, Nuclear and General, vol. 6, no. 11, pp. 1635–1646, Nov. 1973.

[25] E. A. Marengo and A. J. Devaney, “Time-dependent plane wave and multipole expan-
sions of the electromagnetic field,” Journal of Mathematical Physics, vol. 39, no. 7, p.
3643, 1998.

[26] A. J. Devaney, “Multipole expansions and plane wave representations of the electro-

magnetic field,” Journal of Mathematical Physics, vol. 15, no. 2, p. 234, Nov. 1974.

[27] E. Heyman and A. J. Devaney, “Time-dependent multipoles and their application for
radiation from volume source distributions,” Journal of Mathematical Physics, vol. 37,
no. 2, pp. 682–692, 1996.

129

[28] A. Shlivinski and E. Heyman, “Time-domain near-field analysis of short-pulse anten-
nas .I. Spherical wave (multipole) expansion,” IEEE Transactions on Antennas and
Propagation, vol. 47, no. 2, pp. 271–279, 1999.

[29] C.-C. Oetting and L. Klinkenbusch, “Near-to-far-field transformation by a time-
domain spherical-multipole analysis,” IEEE Transactions on Antennas and Propaga-
tion, vol. 53, no. 6, pp. 2054–2063, Jun. 2005.

[30] E. Whittaker, “On the partial differential equations of mathematical physics,” Mathe-

matische Annalen, vol. 57, no. 3, pp. 333–355, 1903.

[31] A. A. Ergin, B. Shanker, and E. Michielssen, “Fast analysis of transient acoustic wave
scattering from rigid bodies using the multilevel plane wave time domain algorithm,”
The Journal of the Acoustical Society of America, vol. 107, no. 3, pp. 1168–1178, 2000.

[32] B. Shanker, A. Ergin, and E. Michielssen, “Fast analysis of transient electromagnetic
scattering phenomena using the multilevel plane wave time domain algorithm,” IEEE
Transactions on Antennas and Propagation, vol. 51, no. 3, pp. 628–641, Mar. 2003.

[33] L. Greengard, T. Hagstrom, and S. Jiang, “The solution of the scalar wave equation
in the exterior of a sphere,” Journal of Computational Physics, vol. 274, pp. 191–207,
Oct. 2014.

[34] S. Azizoglu, S. Koc, and O. Buyukdura, “Spherical Wave Expansion of the Time-
Domain Free-Space Dyadic Green’s Function,” IEEE Transactions on Antennas and
Propagation, vol. 52, no. 3, pp. 677–683, Mar. 2004.

[35] C.-T. Tai, Dyadic green functions in electromagnetic theory, 2nd ed.

IEEE Press, 1994.

Piscataway:

[36] E. L. Hill, “The Theory of Vector Spherical Harmonics,” American Journal of Physics,

vol. 22, no. 4, p. 211, Apr. 1954.

[37] B. Carrascal, G. A. Estevez, P. Lee, and V. Lorenzo, “Vector spherical harmonics and
their application to classical electrodynamics,” European Journal of Physics, vol. 12,
no. 4, pp. 184–191, Jul. 1991.

[38] A. J. Pray, Y. Beghein, N. V. Nair, K. Cools, H. Bağcı, and B. Shanker, “A higher order
space-time galerkin scheme for time domain integral equations,” IEEE Transactions
on Antennas and Propagation, vol. 62, no. 12, pp. 6183–6191, 2014.

[39] H. Brunner, P. J. Davies, and D. B. Duncan, “Global convergence and local super-
convergence of first-kind volterra integral equation approximations,” IMA Journal of
Numerical Analysis, 2011.

[40] A. Pray, N. Nair, and B. Shanker, “Stability Properties of the Time Domain Electric
Field Integral Equation Using a Separable Approximation for the Convolution With
the Retarded Potential,” IEEE Transactions on Antennas and Propagation, vol. 60,
no. 8, pp. 3772–3781, Aug. 2012.

130

[41] D. Weile, G. Pisharody, N.-W. Chen, B. Shanker, and E. Michielssen, “A Novel Scheme
for the Solution of the Time-Domain Integral Equations of Electromagnetics,” IEEE
Transactions on Antennas and Propagation, vol. 52, no. 1, pp. 283–295, Jan. 2004.

[42] S. P. Walker, M. J. Bluck, and I. Chatzis, “The stability of integral equation time-
domain scattering computations for three-dimensional scattering; similarities and
differences between electrodynamic and elastodynamic computations,” International
Journal of Numerical Modelling: Electronic Networks, Devices and Fields, vol. 15, no.
5-6, pp. 459–474, 2002.

[43] L. Valle, F. Rivas, and M. F. Catedra, “Combining the moment method with geometri-
cal modelling by nurbs surfaces and bezier patches,” IEEE Transactions on Antennas
and Propagation, vol. 42, no. 3, pp. 373–381, Mar 1994.

[44] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: CAD, finite
elements, NURBS, exact geometry and mesh refinement,” Comput. Methods Appl.
Mech. Eng., vol. 194, no. 39-41, pp. 4135–4195, 2005.

[45] D. L. Dault, N. V. Nair, J. Li, and B. Shanker, “The generalized method of moments
for electromagnetic boundary integral equations,” IEEE Transactions on Antennas and
Propagation, vol. 62, no. 6, pp. 3174–3188, June 2014.

[46] M. Feischl, G. Gantner, and D. Praetorius, “Reliable and efficient a posteriori error
estimation for adaptive {IGA} boundary element methods for weakly-singular integral
equations,” Computer Methods in Applied Mechanics and Engineering, vol. 290, pp.
362 – 386, 2015.

[47] A. F. Peterson and A. F., “Cell Curvature and Far-Field Superconvergence in Numeri-
cal Solutions of Electromagnetic Integral Equations,” International Journal of Anten-
nas and Propagation, vol. 2016, pp. 1–7, 2016.

[48] C. Loop, “Smooth Subdivision Surfaces Based on Triangles,” M.S. thesis, University

of Utah, 1987.

[49] J. Stam, “Evaluation of loop subdivision surfaces,” in Computer Graphics Proceedings

ACM SIGGRAPH, 1998.

[50] F. Cirak, M. Ortiz, and P. Schröder, “Subdivision surfaces: a new paradigm for thin-
shell finite-element analysis,” International Journal for Numerical Methods in Engi-
neering, vol. 47, no. 12, pp. 2039–2072, Apr. 2000.

[51] G. C. Hsiao and R. E. Kleinman, “Mathematical foundations for error estimation in
numerical solutions of integral equations in electromagnetics,” IEEE Transactions on
Antennas and Propagation, vol. 45, no. 3, pp. 316–328, Mar 1997.

[52] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary

shape,” IEEE Trans. Antennas Propag., vol. 30, no. 3, pp. 409–418, May 1982.

131

[53] W. Cai, T. Yu, H. Wang, and Y. Yu, “High-order mixed RWG basis functions for elec-
tromagnetic applications,” IEEE Transactions on Microwave Theory and Techniques,
vol. 49, no. 7, pp. 1295–1303, Jul. 2001.

[54] A. D. Hellicar, J. S. Kot, G. James, and G. K. Cambrell, “A Comparison of Higher
Order Nodal- and Edge-Basis Functions in the MFIE on Rational BÉzier Geometries,”
IEEE Transactions on Antennas and Propagation, vol. 56, no. 6, pp. 1812–1818, Jun.
2008.

[55] R. Rieben, G. Rodrigue, and D. White, “A high order mixed vector finite element
method for solving the time dependent Maxwell equations on unstructured grids,”
Journal of Computational Physics, vol. 204, no. 2, pp. 490–519, Apr. 2005.

[56] E. Jorgensen, J. L. Volakis, P. Meincke, and O. Breinbjerg, “Higher order hierarchical
legendre basis functions for electromagnetic modeling,” IEEE Transactions on Anten-
nas and Propagation, vol. 52, no. 11, pp. 2985–2995, Nov 2004.

[57] R. D. Graglia, D. R. Wilton, and A. F. Peterson, “Higher order interpolatory vec-
tor bases for computational electromagnetics,” IEEE Transactions on Antennas and
Propagation, vol. 45, no. 3, pp. 329–342, Mar 1997.

[58] R. Sevilla, S. Fernández-Méndez, and A. Huerta, “3D NURBS-enhanced finite ele-
ment method (NEFEM),” International Journal for Numerical Methods in Engineer-
ing, vol. 88, no. 2, pp. 103–125, Oct. 2011.

[59] J. Nedelec, “Mixed finite elements in R3,” Numerische Mathematik, vol. 35, no. 3, pp.

315–341, 1980.

[60] J. Cottrell, A. Reali, Y. Bazilevs, and T. Hughes, “Isogeometric analysis of structural
vibrations,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no.
41-43, pp. 5257–5296, Aug. 2006.

[61] Y. Bazilevs, M.-C. Hsu, and M. A. Scott, “Isogeometric fluid–structure interaction
analysis with emphasis on non-matching discretizations, and with application to wind
turbines,” Computer Methods in Applied Mechanics and Engineering, vol. 249–252,
no. 0, pp. 28–41, 2012.

[62] L. De Lorenzis, J. Evans, T. Hughes, and A. Reali, “Isogeometric collocation: Neu-
mann boundary conditions and contact,” Computer Methods in Applied Mechanics and
Engineering, vol. 284, pp. 21–54, Feb. 2015.

[63] K. Nordanger, R. Holdahl, A. M. Kvarving, A. Rasheed, and T. Kvamsdal, “Implemen-
tation and comparison of three isogeometric Navier–Stokes solvers applied to simula-
tion of flow past a fixed 2D NACA0012 airfoil at high Reynolds number,” Computer
Methods in Applied Mechanics and Engineering, vol. 284, pp. 664–688, Feb. 2015.

[64] P. Kang and S.-K. Youn, “Isogeometric analysis of topologically complex shell struc-

tures,” Finite Elements in Analysis and Design, vol. 99, pp. 68–81, Jul. 2015.

132

[65] R. Bouclier, T. Elguedj, and A. Combescure, “An isogeometric locking-free NURBS-
based solid-shell element for geometrically nonlinear analysis,” International Journal
for Numerical Methods in Engineering, vol. 101, no. 10, pp. 774–808, Mar. 2015.

[66] R. Simpson, M. Scott, M. Taus, D. Thomas, and H. Lian, “Acoustic isogeometric
boundary element analysis,” Computer Methods in Applied Mechanics and Engineering,
vol. 269, pp. 265–290, Feb. 2014.

[67] A. Buffa, G. Sangalli, and R. Vázquez, “Isogeometric analysis in electromagnetics:
B-splines approximation,” Computer Methods in Applied Mechanics and Engineering,
vol. 199, no. 17-20, pp. 1143–1152, Mar. 2010.

[68] K. Kostas, A. Ginnis, C. Politis, and P. Kaklis, “Ship-hull shape optimization with
a T-spline based BEM-isogeometric solver,” Computer Methods in Applied Mechanics
and Engineering, Oct. 2014.

[69] G. Kuru, C. V. Verhoosel, K. G. van der Zee, and E. H. van Brummelen, “Goal-
adaptive Isogeometric Analysis with hierarchical splines,” Computer Methods in Ap-
plied Mechanics and Engineering, vol. 270, no. 0, pp. 270–292, 2014.

[70] J. Kiendl, R. Schmidt, R. Wüchner, and K.-U. Bletzinger, “Isogeometric shape opti-
mization of shells using semi-analytical sensitivity analysis and sensitivity weighting,”
Computer Methods in Applied Mechanics and Engineering, vol. 274, no. 0, pp. 148–167,
2014.

[71] T. DeRose, M. Kass, and T. Truong, “Subdivision surfaces in character animation,”
in Proceedings of the 25th Annual Conference on Computer Graphics and Interactive
Techniques, ser. SIGGRAPH ’98. New York, NY, USA: ACM, 1998, pp. 85–94.

[72] Y. Bazilevs, V. M. Calo, J. A. Cottrell, J. A. Evans, T. J. R. Hughes, S. Lipton, M. A.
Scott, and T. W. Sederberg, “Isogeometric analysis using T-splines,” Comput. Methods
Appl. Mech. Eng., vol. 199, no. 5-8, pp. 229–263, 2010.

[73] M. R. Dörfel, B. Jüttler, and B. Simeon, “Adaptive isogeometric analysis by local h-
refinement with T-splines,” Computer Methods in Applied Mechanics and Engineering,
vol. 199, no. 5-8, pp. 264–275, Jan. 2010.

[74] P. J. Barendrecht, “IsoGeometric Analysis with Subdivision Surfaces,” Ph.D. disserta-

tion, Eindhoven University of Technology, 2013.

[75] R. Simpson, S. Bordas, J. Trevelyan, and T. Rabczuk, “A two-dimensional Isogeometric
Boundary Element Method for elastostatic analysis,” Computer Methods in Applied
Mechanics and Engineering, vol. 209-212, pp. 87–100, Feb. 2012.

[76] M. Scott, R. Simpson, J. Evans, S. Lipton, S. Bordas, T. Hughes, and T. Seder-
berg, “Isogeometric boundary element analysis using unstructured T-splines,” Com-
puter Methods in Applied Mechanics and Engineering, vol. 254, pp. 197–221, Feb.
2013.

133

[77] R. Vazquez, A. Buffa, and L. Di Rienzo, “NURBS-Based BEM Implementation of High-
Order Surface Impedance Boundary Conditions,” IEEE Transactions on Magnetics,
vol. 48, no. 12, pp. 4757–4766, Dec. 2012.

[78] J. Mautz and R. F. Harrington, “H-field, e-field and combined field solution for con-
ducting bodies of revolution,” Arch. Elektron. Übertragung.tech., vol. 32, no. 4, pp.
157–164, 1978.

[79] J. Mautz and R. Harrington, “A combined-source solution for radiation and scattering
from a perfectly conducting body,” IEEE Transactions on Antennas and Propagation,
vol. 27, no. 4, pp. 445–454, Jul 1979.

[80] A. D. Yaghjian, “Augmented electric- and magnetic-field integral equations,” Radio

Science, vol. 16, no. 6, pp. 987–1001, 1981.

[81] M. Taskinen and P. Yla-Oijala, “Current and charge integral equation formulation,”
Antennas and Propagation, IEEE Transactions on, vol. 54, no. 1, pp. 58–67, Jan 2006.

[82] Q. Pan, G. Xu, G. Xu, and Y. Zhang, “Isogeometric analysis based on extended loop’s

subdivision,” Journal of Computational Physics, vol. 299, pp. 731 – 746, 2015.

[83] A. Bendali, “Numerical analysis of the exterior boundary value problem for the time-
harmonic maxwell equations by a boundary finite element method part 1: The contin-
uous problem,” Mathematics of Computation, vol. 43, no. 167, pp. 29–46, 1984.

[84] ——, “Numerical analysis of the exterior boundary value problem for the time- har-
monic maxwell equations by a boundary finite element method part 2: The discrete
problem,” Mathematics of Computation, vol. 43, no. 167, pp. 47–68, 1984.

[85] A. Bespalov and N. Heuer, “The hp-bem with quasi-uniform meshes for the electric field
integral equation on polyhedral surfaces: A priori error analysis,” Applied Numerical
Mathematics, vol. 60, no. 7, pp. 705 – 718, 2010.

[86] W.-L. Wu, A. W. Glisson, and D. Kajfez, “A study of two numerical solution proce-
dures for the electric field integral equation at low frequency,” Applied Computational
Electromagnetics Society Journal, vol. 10, no. 3, pp. 69–80, 1995.

[87] G. Vecchi, “Loop-star decomposition of basis functions in the discretization of the
EFIE,” IEEE Transactions on Antennas and Propagation, vol. 47, no. 2, pp. 339–346,
1999.

[88] F. P. Andriulli, “Loop-Star and Loop-Tree Decompositions: Analysis and Efficient
Algorithms,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 5, pp.
2347–2356, May 2012.

[89] J.-F. Lee and R. Burkholder, “Loop star basis functions and a robust preconditioner for
efie scattering problems,” Antennas and Propagation, IEEE Transactions on, vol. 51,
no. 8, pp. 1855–1863, Aug 2003.

134

[90] Z.-G. Qian and W. C. Chew, “Enhanced A-EFIE With Perturbation Method,” IEEE
Transactions on Antennas and Propagation, vol. 58, no. 10, pp. 3256–3264, Oct. 2010.

[91] C. L. Epstein and L. Greengard, “Debye Sources and the Numerical Solution of the
Time Harmonic Maxwell Equations,” COMMUNICATIONS ON PURE AND AP-
PLIED MATHEMATICS, vol. 63, no. 4, pp. 413–463, APR 2010.

[92] C. L. Epstein, Z. Gimbutas, L. Greengard, A. Klockner, and M. O’Neil, “A consistency
condition for the vector potential in multiply-connected domains,” IEEE Transactions
on Magnetics, vol. 49, no. 3, pp. 1072–1076, March 2013.

[93] X. Y. Z. Xiong, L. J. Jiang, W. E. I. Sha, and Y.-H. Lo, “A New Efie Method Based on
Coulomb Gauge for the Low-Frequency Electromagnetic Analysis,” Prog. Electromagn.
Res., vol. 140, no. April, pp. 613–631, 2013.

[94] F. Vipiana, P. Pirinoli, and G. Vecchi, “Spectral properties of the efie-mom matrix
for dense meshes with different types of bases,” Antennas and Propagation, IEEE
Transactions on, vol. 55, no. 11, pp. 3229–3238, Nov 2007.

[95] F. Andriulli, “Well-posed boundary element formulations in electromagnetics,” Ph.D.
dissertation, University of Michigan, 2008, copyright - Copyright ProQuest, UMI Dis-
sertations Publishing 2008; Last updated - 2014-01-08; First page - n/a.

[96] S. H. Christiansen and J.-C. Nédélec, “A Preconditioner for the Electric Field Integral
Equation Based on Calderon Formulas,” SIAM Journal on Numerical Analysis, vol. 40,
no. 3, pp. 1100–1135, Jan. 2002.

[97] F. P. Andriulli, K. Cools, H. Bagci, F. Olyslager, A. Buffa, S. Christiansen, and
E. Michielssen, “A Multiplicative Calderon Preconditioner for the Electric Field In-
tegral Equation,” IEEE Transactions on Antennas and Propagation, vol. 56, no. 8, pp.
2398–2412, Aug. 2008.

[98] S. Yan, J.-M. Jin, and Z. Nie, “EFIE Analysis of Low-Frequency Problems With Loop-
Star Decomposition and Calderón Multiplicative Preconditioner,” IEEE Transactions
on Antennas and Propagation, vol. 58, no. 3, pp. 857–867, Mar. 2010.

[99] A. Buffa and S. H. Christiansen, “A dual finite element complex on the barycentric
refinement,” MATHEMATICS OF COMPUTATION, vol. 76, no. 260, pp. 1743–1769,
2007.

[100] R. Ling, W. Wang, and D. Yan, “Fitting sharp features with loop subdivision surfaces,”

Computer Graphics Forum, vol. 27, no. 5, pp. 1383–1391, 2008.

[101] D. Colton and R. Kress, Integral Equation Methods in Scattering Theory. John Wiley

& Sons, 1983.

[102] J. Song, C.-C. Lu, and W. C. Chew, “Multilevel fast multipole algorithm for electro-
magnetic scattering by large complex objects,” IEEE Transactions on Antennas and
Propagation, vol. 45, no. 10, pp. 1488–1493, Oct 1997.

135

[103] F. Vico, M. Ferrando, L. Greengard, and Z. Gimbutas, “The decoupled potential inte-
gral equation for time-harmonic electromagnetic scattering,” Communications on Pure
and Applied Mathematics, vol. 69, no. 4, pp. 771–812, 2016.

[104] Z. G. Qian and W. C. Chew, “A quantitative study on the low frequency breakdown of
efie,” Microwave and Optical Technology Letters, vol. 50, no. 5, pp. 1159–1162, 2008.

[105] J. S. Zhao and W. C. Chew, “Integral equation solution of Maxwell’s equations from
zero frequency to microwave frequencies,” IEEE Trans. Antennas Propag., vol. 48,
no. 10, pp. 1635–1645, 2000.

[106] N. V. Nair and B. Shanker, “Generalized method of moments: A novel discretization
technique for integral equations,” IEEE Trans. Antennas Propag., vol. 59, no. 6 PART
2, pp. 2280–2293, 2011.

[107] J. R. Mautz and R. F. Harrington, “H-field, e-field, and combined field solutions for

bodies of revolution,” DTIC Document, Tech. Rep., 1977.

[108] M. Taskinen and P. Ylä-Oijala, “Current and charge integral equation formulation,”

IEEE Trans. Antennas Propag., vol. 54, no. 1, pp. 58–67, 2006.

[109] H. Bagci, F. P. Andriulli, K. Cools, F. Olyslager, and E. Michielssen, “A calderon mul-
tiplicative preconditioner for the combined field integral equation,” IEEE Transactions
on Antennas and Propagation, vol. 57, no. 10, pp. 3387–3392, Oct 2009.

[110] F. P. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, “On a well-conditioned elec-
tric field integral operator for multiply connected geometries,” IEEE Transactions on
Antennas and Propagation, vol. 61, no. 4, pp. 2077–2087, April 2013.

[111] R. E. Kleinman and P. A. Martin, “On single integral equations for the transmission
problem of acoustics,” SIAM Journal on Applied Mathematics, vol. 48, no. 2, pp. 307–
325, 1988.

[112] C. L. Epstein, L. Greengard, and M. O’Neil, “Debye sources and the numerical solution
of the time harmonic maxwell equations ii,” Communications on Pure and Applied
Mathematics, vol. 66, no. 5, pp. 753–789, 2013.

[113] E. V. Chernokozhin and A. Boag, “Method of generalized debye sources for the analysis
of electromagnetic scattering by perfectly conducting bodies with piecewise smooth
boundaries,” IEEE Trans. Antennas Propag., vol. 61, no. 4, pp. 2108–2115, 2013.

[114] F. Vico, M. Ferrando-Bataller, T. Jiménez, and A. Berenguer, “A high order locally
corrected Nyström implementation of the Decoupled Potential Integral Equation,” in
2015 9th European Conference on Antennas and Propagation (EuCAP), May 2015, pp.
1–4.

[115] X. Fu, J. Li, L. J. Jiang, and B. Shanker, “Generalized debye sources-based efie solver
on subdivision surfaces,” IEEE Transactions on Antennas and Propagation, vol. 65,
no. 10, pp. 5376–5386, Oct 2017.

136

[116] Z. Peng, K. H. Lim, and J. F. Lee, “A discontinuous galerkin surface integral equa-
tion method for electromagnetic wave scattering from nonpenetrable targets,” IEEE
Transactions on Antennas and Propagation, vol. 61, no. 7, pp. 3617–3628, July 2013.

[117] A. Buffa, G. Sangalli, and R. Vázquez, “Isogeometric methods for computational elec-
tromagnetics: B-spline and t-spline discretizations,” Journal of Computational Physics,
vol. 257, Part B, pp. 1291 – 1320, 2014.

[118] Y. BAZILEVS, L. BEIRÃO DA VEIGA, J. A. COTTRELL, T. J. R. HUGHES, and
G. SANGALLI, “Isogeometric analysis: Approximation, stability and error estimates
for h-refined meshes,” Mathematical Models and Methods in Applied Sciences, vol. 16,
no. 07, pp. 1031–1090, 2006.

[119] B. Juttler, A. Mantzaflaris, R. Perl, and M. Rumpf, “On numerical integration in
isogeometric subdivision methods for pde on surfaces,” Computer Methods in Applied
Mechanics and Engineering, vol. 302, pp. 131 – 146, 2016.

[120] M. Costabel, E. Darrigrand, and E. Kone, “Volume and surface integral equations for
electromagnetic scattering by a dielectric body,” Journal of Computational and Applied
Mathematics, vol. 234, no. 6, pp. 1817–1825, 2010.

[121] A. J. Poggio and E. K. Miller, Integral equation solutions of three-dimensional scatter-

ing problems. MB Assoc., 1970.

[122] C. Müller, Foundations of

the Mathematical Theory of Electromagnetic Waves.

Springer-Verlag Berlin Heidelberg, 1969.

[123] K. Cools, F. P. Andriulli, F. Olyslager, and E. Michielssen, “Time domain calderon
identities and their application to the integral equation analysis of scattering by pec
objects part i: Preconditioning,” IEEE Transactions on Antennas and Propagation,
vol. 57, no. 8, pp. 2352–2364, Aug 2009.

[124] ——, “Nullspaces of mfie and calderon preconditioned efie operators applied to toroidal
surfaces,” IEEE Transactions on Antennas and Propagation, vol. 57, no. 10, pp. 3205–
3215, Oct 2009.

[125] R. J. Adams, “Physical and analytical properties of a stabilized electric field integral
equation,” IEEE Transactions on Antennas and Propagation, vol. 52, no. 2, pp. 362–
372, Feb 2004.

[126] Z. G. Qian and W. C. Chew, “An augmented electric field integral equation for high-
speed interconnect analysis,” Microwave and Optical Technology Letters, vol. 50, no. 10,
pp. 2658–2662, 2008.

[127] J. Cheng, R. J. Adams, J. C. Young, and M. A. Khayat, “Augmented efie with nor-
mally constrained magnetic field and static charge extraction,” IEEE Transactions on
Antennas and Propagation, vol. 63, no. 11, pp. 4952–4963, 2015.

137

[128] J. Li, X. Fu, and B. Shanker, “Formulation and Iso-geometric Analysis of Scalar Inte-
gral Equations for Electromagnetic Scattering,” IEEE Transactions on Antennas and
Propagation, submitted.

[129] J. Cheng and R. J. Adams, “Electric Field-Based Surface Integral Constraints for
Helmholtz Decompositions of the Current on a Conductor,” IEEE Transactions on
Antennas and Propagation, vol. 61, no. 9, pp. 4632–4640, Sep. 2013.

[130] W. C. Chew, “Vector Potential Electromagnetics with Generalized Gauge for Inho-
mogeneous Media: Formulation,” Progress In Electromagnetics Research, vol. 149, pp.
69–84, 2014.

[131] Q. S. Liu, S. Sun, and W. C. Chew, “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, July 2015, pp. 744–
745.

[132] F. Vico, L. Greengard, and M. Ferrando, “Decoupled field integral equations for elec-

tromagnetic scattering from homogeneous penetrable obstacles,” arxiv, Apr 2017.

[133] R. F. Harrington, “Boundary integral formulations for homogeneous material bodies,”

Journal of Electromagnetic Waves and Applications, vol. 3, no. 1, pp. 1–15, 1989.

[134] A. Stewart, “Role of the nonlocality of the vector potential in the aharonov–bohm

effect,” Canadian Journal of Physics, vol. 91, no. 5, pp. 373–377, 2013.

[135] G. Trammel, “Aharonov-bohm paradox,” Physical Review, vol. 134, no. 5B, p. B1183,

1964.

[136] A. Burton and G. Miller, “The application of integral equation methods to the numer-
ical solution of some exterior boundary-value problems,” in Proceedings of the Royal
Society of London A: Mathematical, Physical and Engineering Sciences, vol. 323, no.
1553. The Royal Society, 1971, pp. 201–210.

[137] K. Chen, “An analysis of sparse approximate inverse preconditioners for boundary
integral equations,” SIAM Journal on Matrix Analysis and Applications, vol. 22, no. 4,
pp. 1058–1078, 2001.

[138] A. Bunse-Gerstner and I. Gutiérrez-Cañas, “A preconditioned gmres for complex dense
linear systems from electromagnetic wave scattering problems,” Linear algebra and its
applications, vol. 416, no. 1, pp. 135–147, 2006.

[139] B. Carpentieri, I. S. Duff, L. Giraud, and G. Sylvand, “Combining fast multipole
techniques and an approximate inverse preconditioner for large electromagnetism cal-
culations,” SIAM Journal on Scientific Computing, vol. 27, no. 3, pp. 774–792, 2005.

138