FullText (40)

8/17/2019 FullText (40)

1/142

TitleSurface integral equation method for analyzing electromagneticscattering in layered medium

Advisor(s) Chew, WC; Jiang, L

Author(s) Chen, Yongpin.;– Œ˜‘

.

Citation

Issued Date 2011

URL http://hdl.handle.net/10722/174463

RightsThe author retains all proprietary rights, (such as patent rights)and the right to use in future works.

8/17/2019 FullText (40)

2/142

Surface Integral Equation Method for AnalyzingElectromagnetic Scattering in Layered Medium

by

Yongpin CHEN (陳涌頻)

A thesis submitted in partial fulfillment of the requirements for the degree of

Doctor of Philosophy at The University of Hong Kong

October 2011

8/17/2019 FullText (40)

3/142

Abstract of thesis entitled

Surface Integral Equation Method for Analyzing

Electromagnetic Scattering in Layered Medium

Submitted by

Yongpin CHEN

for the degree of Doctor of Philosophy

at The University of Hong Kong

in October 2011

Surface integral equation (SIE) method with the kernel of layered medium Green's

function (LMGF) is investigated in details from several fundamental aspects. A

novel implementation of discrete complex image method (DCIM) is developed to

accelerate the evaluation of Sommerfeld integrals and especially improve the far

field accuracy of the conventional one. To achieve a broadband simulation of thin

layered structure such as microstrip antennas, the mixed-form thin-stratified

medium fast-multipole algorithm (MF-TSM-FMA) is developed by applyingcontour deformation and combining the multipole expansion and plane wave

expansion into a single multilevel tree. The low frequency breakdown of the

integral operator is further studied and remedied by using the loop-tree

decomposition and the augmented electric field integral equation (A-EFIE), both

in the context of layered medium integration kernel. All these methods are based

on the EFIE for the perfect electric conductor (PEC) and hence can be applied in

antenna and circuit applications. To model general dielectric or magnetic objects,

the layered medium Green's function based on pilot vector potential approach is

generalized for both electric and magnetic current sources. The matrix

representation is further derived and the corresponding general SIE is setup.

Finally, this SIE is accelerated with the DCIM and applied in quantum optics,

such as the calculation of spontaneous emission enhancement of a quantum

emitter embedded in a layered structure and in the presence of nano scatterers.

8/17/2019 FullText (40)

4/142

Declaration

I declare that the thesis and the research work thereof represents my own work,

except where due acknowledgement is made, and that it has not been previously

included in a thesis, dissertation or report submitted to this University or to any

other institution for a degree, diploma or other qualifications.

Signed

Yongpin CHEN

8/17/2019 FullText (40)

5/142

c 2011 Yongpin CHEN

8/17/2019 FullText (40)

6/142

SURFACE INTEGRAL EQUATION METHOD FOR ANALYZINGELECTROMAGNETIC SCATTERING IN LAYERED MEDIUM

BY

YONGPIN CHEN

DISSERTATION

Submitted in partial fulfillment of the requirements

for the degree of Doctor of Philosophy in Electrical and Electronic Engineering

in the Graduate School of the

University of Hong Kong, 2011

Hong Kong

8/17/2019 FullText (40)

7/142

ABSTRACT

Surface integral equation (SIE) method with the kernel of layered medium

Green’s function (LMGF) is investigated in details from several fundamental

aspects. A novel implementation of discrete complex image method (DCIM)

is developed to accelerate the evaluation of Sommerfeld integrals and espe-

cially improve the far field accuracy of the conventional one. To achieve abroadband simulation of thin layered structure such as microstrip antennas,

the mixed-form thin-stratified medium fast-multipole algorithm (MF-TSM-

FMA) is developed by applying contour deformation and combining the mul-

tipole expansion and plane wave expansion into a single multilevel tree. The

low frequency breakdown of the integral operator is further studied and reme-

died by using the loop-tree decomposition and the augmented electric field

integral equation (A-EFIE), both in the context of layered medium integra-

tion kernel. All these methods are based on the EFIE for the perfect electric

conductor (PEC) and hence can be applied in antenna and circuit applica-

tions. To model general dielectric or magnetic objects, the layered medium

Green’s function based on pilot vector potential approach is generalized for

both electric and magnetic current sources. The matrix representation is

further derived and the corresponding general SIE is setup. Finally, this

SIE is accelerated with the DCIM and applied in quantum optics, such as

the calculation of spontaneous emission enhancement of a quantum emitter

embedded in a layered structure and in the presence of nano scatterers.

ii

8/17/2019 FullText (40)

8/142

To my wife and my parents, for their love and support.

iii

8/17/2019 FullText (40)

9/142

ACKNOWLEDGMENTS

First and foremost, I would like to thank my mentor and thesis adviser,

Professor Weng Cho CHEW, for providing me this opportunity to study

with him. His patient and careful guidance during the last four years has

totally changed my way of thinking, re-shaped my knowledge structure, and

finally built my self-confidence. His enthusiasm, persistence, and diligentwork ethic has inspired me much and will continue to affect my career and

life.

I would also like to thank my co-adviser, Dr. Lijun JIANG, who always

encouraged me when I was frustrated, shared with me his working experience

in industry, and helped me a lot in technical details.

I feel grateful to my former supervisors Professor Zaiping NIE and Pro-

fessor Jun HU at UESTC, who opened me the door of the fascinating CEM

world, encouraged and supported me to pursue study at HKU, and always

pay kind attention to my recent status.

Many thanks to my friends and colleagues with the Electromagnetics and

Optics Laboratory, to Dr. Wallace C. H. CHOY for providing me TA op-

portunities and inviting me to audit his course on organic devices, to Dr.

Sheng SUN for giving me suggestions on possible directions about circuit

simulation, to Dr. Wei SHA for sharing his knowledge on modern physics,

to Dr. Yang LIU for teaching me lots of mathematics, to Dr. Yat Hei LO for

sharing his expertise on computer and optics, to Dr. Shaoying HUANG, Dr.

Min TANG, Dr. Osman GONI, Dr. Bo ZHU, Phillip ATKINS, Peng YANG,Qi DAI, Jun HUANG, Zuhui MA, Shiquan HE, Yumao WU, Yan LI, Nick

HUANG, Ping LI for their help and friendship.

I also wish to thank alumni from both UIUC and UESTC. A special thanks

to Dr. Yuan LIU for helping me a lot when I first came to HKU, though we

have never met with each other, to Dr. Liming XU for sharing his code and

teaching me the transmission line method, to Dr. Jie XIONG for sharing her

iv

8/17/2019 FullText (40)

10/142

code and polishing my understanding on many important concepts, to Dr.

Su YAN for working together to develop the MLFMA code, to Dr. Zhiguo

QIAN for his help on A-EFIE.

I thank my parents and my little sister for their love and support. Though

I cannot go back to my hometown every year, their encouragement through

the telephone line always warms my heart.

Finally, I thank my wife, Min MENG, for her love, her sharing of joy

and sadness, and her shouldering of so much responsibility alone when I was

pursuing my study.

v

8/17/2019 FullText (40)

11/142

TABLE OF CONTENTS

LIST OF TA B LES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v i i i

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

LIST OF ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . xiii

CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Organization of the Thesis . . . . . . . . . . . . . . . . . . . . 3

CHAPTER 2 A NOVEL IMPLEMENTATION OF DISCRETECOMPLEX IMAGE METHOD (DCIM) . . . . . . . . . . . . . . . 42.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Layered Medium Green’s Function (LMGF) . . . . . . . . . . 62.3 Discrete Complex Image Method (DCIM) . . . . . . . . . . . 92.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

CHAPTER 3 MIXED-FORM THIN-STRATIFIED MEDIUM FAST-MULTIPOLE ALGORITHM (MF-TSM-FMA) . . . . . . . . . . . 193.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.2 Review of Thin-Stratified Medium Fast-Multipole Algo-

rithm (TSM-FMA) . . . . . . . . . . . . . . . . . . . . . . . . 213.3 Mixed-Form Thin-Stratified Medium Fast-Multipole Algo-

rithm (MF-TSM-FMA) . . . . . . . . . . . . . . . . . . . . . . 243.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

CHAPTER 4 REMEDIES FOR LOW FREQUENCY BREAK-DOWN OF INTEGRAL OPERATOR IN LAYERED MEDIUM . . 374.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374.2 Loop-Tree Decomposition and Frequency Normalization . . . . 394.3 Basis Rearrangement via Connection Matrix . . . . . . . . . . 424.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . 454.5 Augmented Electric Field Integral Equation (A-EFIE) . . . . 45

vi

8/17/2019 FullText (40)

12/142

4.6 Charge Neutrality Issue . . . . . . . . . . . . . . . . . . . . . 514.7 Electrostatic Limit . . . . . . . . . . . . . . . . . . . . . . . . 564.8 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . 594.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

CHAPTER 5 A NEW LAYERED MEDIUM GREEN’S FUNC-TION (LMGF) FORMULATION FOR GENERAL OBJECTS . . . 635.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.2 Dyadic Form of LMGF . . . . . . . . . . . . . . . . . . . . . . 655.3 Matrix Representation of LMGF . . . . . . . . . . . . . . . . . 735.4 Duality Principle of LMGF . . . . . . . . . . . . . . . . . . . . 825.5 Surface Integral Equation . . . . . . . . . . . . . . . . . . . . 845.6 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . 855.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

CHAPTER 6 DCIM ACCELERATED SURFACE INTEGRALEQUATION (SIE) METHOD FOR NANO-OPTICAL APPLI-CATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 906.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 906.2 SIE with DCIM Acceleration . . . . . . . . . . . . . . . . . . . 916.3 Green’s Function Approach in Spontaneous Emission (SE) . . 956.4 Surface Plasmon Polaritons (SPPs) and Localized Surface

Plasmons (LSP) . . . . . . . . . . . . . . . . . . . . . . . . . . 986.5 Quantum Emitter Coupled to Surface Plasmon Resonance . . 1016.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

CHAPTER 7 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . 106

R E F E R E N C E S . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 0 7

LIST OF PUBLICATIONS . . . . . . . . . . . . . . . . . . . . . . . . 119

AUTHOR’S BIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . 122

vii

8/17/2019 FullText (40)

13/142

LIST OF TABLES

2.1 CPU time comparison (seconds) . . . . . . . . . . . . . . . . . 15

3.1 Layout of radiation and receiving patterns . . . . . . . . . . . 27

5.1 Matrix element of KE: no line integral (×103) . . . . . . . . . 795.2 Matrix element of

KE: testing line integral (

×103) . . . . . . . 79

5.3 Matrix element of KE: testing & basis line integrals (×104) . . 816.1 CPU time for single frequency calculation . . . . . . . . . . . 102

viii

8/17/2019 FullText (40)

14/142

LIST OF FIGURES

2.1 A general layered medium. . . . . . . . . . . . . . . . . . . . . 62.2 A typical microstrip structure (unit: mm). . . . . . . . . . . . 112.3 One branch-cut case: Sommerfeld integration path (SIP),

Sommerfeld branch cut (SBC) and poles in the complex kρplane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.4 Two branch-cut case: Sommerfeld integration path (SIP),Sommerfeld branch cuts (SBC) and poles in the complexkρ plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.5 The magnitude of gss versus k0ρ for the microstrip struc-ture. The asymptotic behavior is 1/ρ2, which agrees withthe theoretical prediction. . . . . . . . . . . . . . . . . . . . . 15

2.6 The magnitude of gφ versus k0ρ for the microstrip struc-ture. The far field is dominated by the pole contribution,which has the asymptotic behavior of 1/

√ ρ. . . . . . . . . . . 16

2.7 A 3-layer model (unit: mm). . . . . . . . . . . . . . . . . . . . 16

2.8 The magnitude of gss versus k0ρ for the 3-layer model. Theasymptotic behavior is 1/ρ2. . . . . . . . . . . . . . . . . . . . 17

2.9 The magnitude of gφ versus k0ρ for the 3-layer model. Theasymptotic behavior is 1/ρ2. . . . . . . . . . . . . . . . . . . . 17

3.1 The worst interaction in applying the addition theorem. . . . . 233.2 Relative error by using plane wave expansion. . . . . . . . . . 233.3 Relative error by using multipole expansion. . . . . . . . . . . 283.4 The MF-TSM-FMA structure. . . . . . . . . . . . . . . . . . . 303.5 The geometrical structure of the 8× 4 microstrip array. . . . . 323.6 The radiation patterns of the microstrip array in E plane

(φ = 0◦) and H plane (φ = 90◦) at 9.42 GHz. . . . . . . . . . . 323.7 The bistatic RCS of the microstrip array at f = 2.5 GHz

and (θi, φi) = (60◦, 0◦). . . . . . . . . . . . . . . . . . . . . . . 33

3.8 A dipole on an EBG structure, top and side view. . . . . . . . 333.9 Radiation pattern of the dipole in E plane, due to the lay-

ered medium assumption, the ground plane is infinitely large. . 343.10 The geometrical model of a 30× 30 microstrip array. . . . . . 343.11 Bistatic RCS of the 30× 30 microstrip array. . . . . . . . . . . 35

ix

8/17/2019 FullText (40)

15/142

3.12 Memory requirement (solid line) and CPU time per itera-tion (dash line) versus number of unknowns. . . . . . . . . . 35

4.1 Local loop and global loop (dashed lines). . . . . . . . . . . . 414.2 RCS of a sphere. . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.3 Geometrical model of the capacitor. . . . . . . . . . . . . . . 444.4 Negative input reactance of the capacitor, ǫr2 = 6.0. . . . . . . 444.5 The geometrical structure of the loop inductor embedded

in a seven-layer medium, unit: mm. The central layer isa magnetic material. A delta gap excitation is applied atthe center of the top arm. . . . . . . . . . . . . . . . . . . . . 53

4.6 The condition number versus frequency for the rectangularloop. The condition number is unbounded when decreas-ing the frequency. Charge neutrality enforcement (CNE)makes the condition number constant. . . . . . . . . . . . . . 54

4.7 The eigenvalue distribution for the rectangular loop at 1Hz. The smallest eigenvalue is removed away from theorigin after the charge neutrality enforcement (CNE). . . . . . 54

4.8 The geometrical model of the half loop embedded in a five-layer medium (including the PEC layer), unit: mm. Adelta gap excitation is applied at the center of the top arm. . . 55

4.9 The condition number versus frequency for the half loop.Since it is connected to the ground plane, charge neutralitycannot be guaranteed. The condition number is boundedwhen decreasing the frequency without any special treatment. 55

4.10 The geometrical model of the circular parallel plate capac-itor, with a dielectric layer (ǫr = 2.65) inserted in between.A delta gap is applied at the edge. The mesh is refined tocapture the fringing effect. . . . . . . . . . . . . . . . . . . . . 57

4.11 The input reactance of the rectangular loop. A-EFIE agreeswell with the loop-tree (LT) decomposition, while the tra-ditional EFIE breaks down quickly when decreasing thefrequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.12 The input reactance of the half loop. A-EFIE maintainsthe scale invariance very well while the traditional EFIEbreaks down quickly when decreasing the frequency. Since

the non-magnetic dielectric is transparent to the inductor,a PEC half space model is applied to validate the results. . . . 61

x

8/17/2019 FullText (40)

16/142

4.13 The capacitance of the circular parallel plate capacitor. A-EFIE I represents the capacitance extracted from current,while A-EFIE Q means the capacitance extracted fromcharge. The A-EFIE current suffers from an inaccuracyproblem while the A-EFIE charge is stable. The result

agrees with the static solver. Both are further validatedby the analytic solution. When the frequency is below 1MHz, the relative error of A-EFIE Q is around 0 .1%. . . . . . 61

5.1 A homogeneous dielectric object is embedded in a layeredmedium. The external excitation is either a plane wave ora Hertzian dipole. Equivalent electric and magnetic cur-rents are induced on the boundary which then generate thescattered field. . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.2 An electric or magnetic dipole is radiating in a seven-layer

medium (unit: m). The layered medium is both dielectricand magnetic, and the layer constants are shown in thefigure. The source point is in Layer 2 and the observationline is in Layer 5. . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.3 The magnitude of electric field generated by an electricdipole. The dipole polarization is (θ = 20o, φ = 30o) andis working at f = 300 MHz. The result is validated by thetransmission line method (solid line). . . . . . . . . . . . . . . 74

5.4 The magnitude of electric field generated by a magneticdipole. The dipole polarization is (θ = 20o, φ = 30o) andis working at f = 300 MHz. The result is validated by thetransmission line method (solid line). . . . . . . . . . . . . . . 74

5.5 Cases where testing line integral exists. The testing func-tion is straddling the interface, all radiation from the RWGor half-RWG basis functions (triangles) in color needs in-voking testing line integral, while radiation from othersdoes not need this line integral. . . . . . . . . . . . . . . . . . 80

5.6 Line integral test. The testing function is at the top in-terface. Radiation from: basis function 1: no line integralactivated; basis function 2: testing line integral activated;basis function 3: both testing and basis line integrals acti-

vated. (unit: m). . . . . . . . . . . . . . . . . . . . . . . . . . 805.7 A homogeneous capsule structure embedded in a 5-layer

medium, where h = 0.6, r = 0.15 (unit: m). . . . . . . . . . . 865.8 The scattered field inside the object with ǫ = 2 and µ = 2.

Since there is no contrast, the scattered field recovers theincident field (solid line). . . . . . . . . . . . . . . . . . . . . . 86

xi

8/17/2019 FullText (40)

17/142

5.9 The scattered field outside the object with ǫ = 4 and µ = 1.The tangential components of the E field are continuous atthe interface. Symbol: field at d+5 , solid line: field at d

−5 . . . . 87

5.10 The scattered field outside the object with ǫ = 4 and µ = 1.The normal component of the E field differs by a factor of

3 (ǫ4/ǫ5) at the interface. . . . . . . . . . . . . . . . . . . . . . 875.11 A homogeneous cuboid embedded in a 4-layer medium

(unit: m). It is penetrating different layers. . . . . . . . . . . . 885.12 The scattered field of the PEC object, calculated from

EFIE and MFIE. . . . . . . . . . . . . . . . . . . . . . . . . . 885.13 The scattered field of the dielectric object. . . . . . . . . . . . 89

6.1 Transition of a two-level atomic system. . . . . . . . . . . . . 966.2 SPPs can be sustained in the interface of dielectric-metal

half space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

6.3 LSP can be excited by a metallic nano sphere under an EMfield excitation. . . . . . . . . . . . . . . . . . . . . . . . . . . 1006.4 A gold nano sphere with radius 20 nm is located above a

gold slab with thickness 30 nm. A z -polarized emitter islocated at the middle of the sphere and slab. . . . . . . . . . . 103

6.5 Normalized SER for the three cases: (a) only slab (SPPsenhanced SE); (b) only nano sphere (LSP enhanced SE);(c) both (both effects). . . . . . . . . . . . . . . . . . . . . . . 103

6.6 A gold nano bowl is located in a layered medium with twogold layers. A z -polarized emitter is located at the centerof the aperture. The dimensions are shown in the figure. . . . 104

6.7 Mesh of the gold nano bowl. . . . . . . . . . . . . . . . . . . . 1046.8 Normalized SER. The first peak corresponds to the peak

from SPPs of the layered medium at 520 nm, the secondpeak is from the nano bowl effect around 590 nm. . . . . . . . 105

xii

8/17/2019 FullText (40)

18/142

LIST OF ABBREVIATIONS

SIE Surface Integral Equation

LMGF Layered Medium Green’s Function

PEC Perfect Electric Conductor

CEM Computational Electromagnetics

FDM Finite Difference Method

FEM Finite Element Method

MoM Method of Moments

GPR Ground-Penetrating Radar

LED Light-Emitting Diode

GMRES Generalized Minimal Residual

DCIM Discrete Complex Image Method

PML Perfectly Matched Layers

RFFM Rational Function Fitting Method

SBC Sommerfeld Branch Cut

SIP Sommerfeld Integration Path

TE Transverse Electric

TM Transverse Magnetic

MF-TSM-FMA Mixed-Form Thin-Stratified Medium Fast-Multipole Algo-rithm

FMM Fast Multipole Method

MLFMA Multilevel Fast Multipole Algorithm

xiii

8/17/2019 FullText (40)

19/142

AIM Adaptive Integral Method

CG-FFT Conjugate Gradient Fast Fourier Transform

MLMDA Multilevel Matrix Decomposition Algorithm

ACA Adaptive Cross Approximation

RCS Radar Cross Section

EBG Electromagnetic Band Gap

CPU Central Processing Unit

PC Personal Computer

A-EFIE Augmented Electric Field Integral Equation

CAD Computer-Aided Design

RWG Rao-Wilton-Glisson

CCIE Current and Charge Integral Equation

SPIE Separated Potential Integral Equation

MPIE Mixed Potential Integral Equation

CNE Charge Neutrality Enforcement

LT Loop-Tree

2D Two-Dimensional

3D Three-Dimensional

VIE Volume Integral Equation

BOR Body of Revolution

TL Transmission Line

PMCHWT Poggio-Miller-Chang-Harrington-Wu-TsaiEFIE Electric Field Integral Equation

MFIE Magnetic Field Integral Equation

SE Spontaneous Emission

SPPs Surface Plasmon Polaritons

LSP Localized Surface Plasmons

xiv

8/17/2019 FullText (40)

20/142

LDOS Local Density of States

EM Electromagnetics

QED Quantum Electrodynamics

PC Photonic Crystal

SER Spontaneous Emission Rate

DOS Density of States

OLED Organic Light-Emitting Diode

MEMS Microelectromechanical Systems

GIBC Generalized Impedance Boundary Condition

xv

8/17/2019 FullText (40)

21/142

CHAPTER 1

INTRODUCTION

1.1 Motivation

With the development of computer technology and scientific computation,

numerical simulation becomes as important as theoretical study and experi-

mental verification in understanding the nature of our world. Computational

electromagnetics (CEM) is one branch of numerical techniques to model all

electromagnetic phenomena governed by the Maxwell’s equations, covering

from circuit engineering, microwave and antenna engineering, to optical en-

gineering.

The methods in CEM can be classified into different categories according to

different criteria, such as different computation domains (time or frequency),

or different equation forms (differential or integral) [1]. The time domain

methods solve the time domain EM fields directly, while the frequency do-main methods solve the harmonic fields based on Fourier transform. On the

other hand, the differential equation methods solve the Maxwell’s equations

in differential form, such as finite difference method (FDM) [2], finite ele-

ment method (FEM) [3], while the integral equation methods first extract

the point response (Green’s function) of a differential equation, and solve the

integral equation formulated from equivalence principle (surface or volume),

such as the method of moments (MoM) [4].

Integral equation method in the frequency domain is one of the most pow-erful methods in analyzing electromagnetic radiation and scattering. The key

of this method is the Green’s function, which describes the EM response of

a point source in a specific environment. It satisfies the boundary condition

automatically and hence does not require any domain truncation or absorp-

tion boundary condition in the computation. Also, if the surface equivalence

principle is applied, the surface integral equation (SIE) can be obtained,

1

8/17/2019 FullText (40)

22/142

where the unknowns are only associated with the boundary of the scatterers

and thus can be made minimum.

Although the SIE can always be derived mathematically, the Green’s func-

tions for arbitrary inhomogeneous medium is not trivial and can only be

obtained numerically, which may be as complex as solving the original prob-

lem. For this reason, most research on surface integral equation methods are

focused on homogeneous environment, such as aerospace applications. How-

ever, if the inhomogeneity is one dimensional piecewise, or so called layered

medium, the Green’s function can be determined analytically in the spec-

tral (Fourier) domain, and the spatial domain counterpart can be obtained

by simply inverse Fourier transforming it [5]. In fact, this scenario con-

sists of a broad class of applications in both microwave and optical regimes,

such as microstrip antennas and microwave circuits, geophysical exploration,ground-penetrating radar (GPR), solar cell, light-emitting diode (LED), and

lithography, etc.

This thesis is focused on the SIE method with layered medium Green’s

function (LMGF), and addresses several fundamental problems of this method.

First, the spatial domain LMGF can only be obtained by numerically in-

tegrating the spectral domain counterpart, expressed as a Sommerfeld in-

tegral [6]. This integral is oscillatory and slowly convergent, which makes

the repeated evaluation in MoM extremely time expensive. An acceleration

technique at the Green’s function level is necessary.

Second, the SIE leads to a full matrix, which requires an O(N 2) memory

and an O(N 3) or O(N 2) computational time if a direct inverse such as LU

decomposition or an iterative method such as the GMRES is applied [7] [8]. A

fast algorithm is indispensable for large scale computation, such as radiation

of an antenna array. The fast algorithm developed shall also cover a broad

frequency range to achieve a broadband simulation.

Third, when the structure under investigation is much smaller then the

wavelength—which is common in circuit simulation—the integral operator

may suffer from a low frequency breakdown. Remedies shall be proposed to

overcome this difficulty for the specific integration kernel of LMGF.

Fourth, though the commonly used metal is a good conductor at the radio

or microwave frequencies and the electric field integral equation (EFIE) [9]

for a perfect electric conductor (PEC) can be implemented in most applica-

tions, the ability of dielectric modeling is important in other situations such

2

8/17/2019 FullText (40)

23/142

as underground exploration or optical applications. For example, metal is

penetrable to the order of wavelength at the optical frequency band, and

hence becomes dielectric. A full sets of LMGFs shall be developed to obtain

a general SIE.

Finally, the power of CEM should not be restricted to classical electrody-

namics [10], but shall also be highlighted in quantum electrodynamics such

as quantum optics [11]. Such attempt will be made by applying the tools

developed in this thesis.

1.2 Organization of the Thesis

The organization of the thesis is as follows: Chapter 2 first reviews one pop-

ular acceleration technique, the discrete complex image method (DCIM), for

fast evaluation of LMGF. The cause of far field inaccuracy of the conven-

tional DCIM is explained and a novel implementation based on Sommerfeld

branch cut is proposed. In Chapter 3, a fast algorithm suitable for mi-

crostrip structures—the quasi-3D thin-stratified medium fast-multipole algo-

rithm (TSM-FMA)—is first introduced, followed by the discussion of low fre-

quency breakdown of it. A mixed-form TSM-FMA (MF-TSM-FMA) comb-

ing the multipole expansion and plane wave expansion is proposed to achieve

a fast and broadband simulation. Next, in Chapter 4, the low frequency

breakdown of the integral operator is further investigated. The loop-tree

decomposition based on quasi-Helmholtz decomposition and the augmented

electric field integral equation (A-EFIE) based on the augmentation tech-

nique in the context of layered medium are both studied. After that, in

Chapter 5, a new LMGF formulation is developed to obtain a general SIE

for dielectric modeling. The matrix representation is derived to reduce the

singularity of LMGF and the line integral issue associated with it is dis-

cussed in details. On top of the SIE solver, in Chapter 6, the DCIM isincorporated and is applied to study the matter-field interaction in a com-

plex environment—spontaneous emission enhancement of a quantum emitter

embedded in a layered medium and in the presence of nano scatterers. Fi-

nally, in Chapter 7, the conclusion is made and several possible research

topics are discussed.

3

8/17/2019 FullText (40)

24/142

CHAPTER 2

A NOVEL IMPLEMENTATION OF

DISCRETE COMPLEX IMAGE METHOD(DCIM)

A novel implementation of discrete complex image method (DCIM) based

on the Sommerfeld branch cut is proposed to accurately capture the far-field

behavior of the layered medium Green’s function (LMGF) as a complement

to the traditional DCIM. By contour deformation, the Green’s function can

be naturally decomposed into branch-cut integration (radiation modes) andpole contributions (guided modes). For branch-cut integration, matrix pencil

method is applied, and the alternative Sommerfeld identity in terms of kz

integration is utilized to get a closed-form solution. The guided modes are

accounted for with a pole-searching algorithm. Both one-branch-cut and

two-branch-cut cases are studied. Several numerical results are presented to

validate this method. For the sake of completeness, the LMGF formulation

will also be briefly summarized in this chapter.

2.1 Introduction

Sommerfeld first investigated a dipole radiating above a half space using

Hertzian potential [6]. This problem was further studied and extended to

general multilayered medium by various researchers [12]–[15], [5]. There are

several approaches to derive the layered medium Green’s function (LMGF),

such as the popular transmission line analog [16], [17], Hertzian potential

approach [6], [18], E z –H z formulation [19], [20], and pilot vector potentialapproach [5], [21], etc. No matter which method is applied, however, the

LMGF can only be expressed as an infinite, oscillatory and slowly-convergent

integral, Sommerfeld integral, making the numerical evaluation process very

inefficient. Several methods have been developed to expedite the evalua-

tion of the LMGF, such as the function approximation approach [22], [23],

the path deformation technique [5], [24], the tabulation and interpolation

4

8/17/2019 FullText (40)

25/142

method [25], [26], the perfectly matched layers (PML) method [27], [28], the

singularity subtraction method [29], and the fast all-modes combined with

numerical modified steepest descent path method [30], [31] etc.

Function approximation in the spectral domain is one of the most popular

methods. In this method, the integration kernel is first approximated by

certain “simple functions”, and the integral is then evaluated in a closed

form by applying relevant integration identities. Though lots of function

approximation techniques are available from a numerical analysis point of

view, those candidates with closed-form identities of the infinite integrals

in our context can finally be utilized. This leads to the following methods:

the complex discrete complex image method (DCIM) [22], [32] (based on the

complex exponential functions), the rational function fitting method (RFFM)

[23], [33] (based on the rational fraction functions), or their combination [34].The popular DCIM has unpredictable errors when the interaction is in

the far field region (ρ ≫ 0). The original sampling path cannot effectivelycapture certain singularities, which correspond to the guided modes or surface

waves and lateral waves physically. Several efforts have been made to remedy

this problem. A two-level approximation [35] was proposed to separate the

sharp-transition region from the smooth-varying region, with higher sampling

rate in the former part to capture the singularities. In [36], the surface-

wave poles are extracted explicitly for a general multilayer medium before

applying the DCIM, which is also suggested in [22]. Other attempts are made

to deform the sampling path to carry more pole singularities. For instance,

in [37], a direct DCIM was developed to push the sampling path closer to

the poles and rely on the matrix pencil method [38] itself to take care of the

singularities. Meanwhile, spatial error control is another big issue in DCIM

and a recent progress can be found in [39].

Though the pole singularities can be captured successfully by the above

methods, the branch-point singularities are still not considered. These singu-

larities contribute to the lateral wave when the source and observation points

are at the interface of the layers [5]. Recently, a three-level DCIM [40] was

proposed to bring the sampling path closer to the branch point to capture

more information of this singularity. In the RFFM, the continuous spectrum

in the far field is also considered based on a vertical path and asymptotic

analysis [41]. In this chapter, we propose an alternative implementation

of the DCIM to capture these singularities. Other than introducing extra

5

8/17/2019 FullText (40)

26/142

Figure 2.1: A general layered medium.

segments of the sampling path to approach the singularities, we simply de-

form the sampling path to the Sommerfeld branch cut (SBC) [5] when ρ is

relatively large. The matrix pencil method is applied to approximate the

function along the SBC, which can be mapped into the real axis in the kz

plane [42]. The pole contributions are accounted for by applying a robust

pole-searching algorithm [43].

For the sake of completeness, the basic formulation of the LMGF such

as the propagation factor and the matrix-friendly formulation will be briefly

reviewed. It can also be found in Chapter 3, Chapter 4 and will be generalized

in Chapter 5.

2.2 Layered Medium Green’s Function (LMGF)

The LMGF is the impulse response of a dipole in this medium. The con-

figuration is shown in Figure 2.1, where each layer is characterized by the

permittivity ǫ and permeability µ (losses are included). The dipole is as-

sumed to be at r′ in layer m and the observation point is at r in layer n.

The LMGF (in dyadic form) can be expressed as [21],

Ḡ(r, r′) = (∇× ẑ )(∇′ × ẑ )gTE(r, r′)

+ 1k2nm

(∇×∇× ẑ )(∇′ ×∇′ × ẑ )gTM(r, r′)(2.1)

where k2nm = ω2ǫnµm. The g

TE/TM(r, r′) is expressed as a Sommerfeld inte-

6

8/17/2019 FullText (40)

27/142

gral,

g(r, r′) = i

8π

+∞−∞

dkρkmzkρ

H (1)0 (kρρ)F (kρ, z , z

′) (2.2)

where F (kρ, z , z ′) is the propagation factor [5], kmz = k2m −

k2ρ, and H (1)0 (kρρ)

is the first kind Hankel function of order 0 (the time convention is assumed

to be e−iωt).

2.2.1 Propagation Factor

The expression of propagation factor depends on the relative positions of the

source point and the observation point. In the following expressions, the

argument kρ is neglected for simplicity.

A. n = mThe source point and the observation point are in the same layer. There

are primary field as well as secondary field. Then

F (z, z ′) = eikmz |z−z′|

+ M̃ m R̃m,m−1eikmz(−2dm+z′+z)

+ M̃ m R̃m,m+1eikmz(2dm+1−z′−z)

+ M̃ m R̃m,m+1 R̃m,m−1eikmz(2dm+1−2dm+z′−z)

+ M̃ m R̃m,m−1 R̃m,m+1eikmz(2dm+1−2dm−z′+z)

(2.3)

where the generalized reflection coefficients and Fresnel reflection coefficients

[5] are

R̃i,i+1 = Ri,i+1 + R̃i+1,i+2e

iki+1,z(2di+2−2di+1)

1 + Ri,i+1 R̃i+1,i+2eiki+1,z(2di+2−2di+1)(2.4)

R̃i,i−1 = R

i,i−1 + R̃

i−1,i−2eiki−1,z(2di−2di−1)

1 + Ri,i−1 R̃i−1,i−2eiki−1,z(2di−2di−1)(2.5)

Rij = p jkiz − pik jz p jkiz + pik jz

(2.6)

p =

µ, TE

ǫ, TM(2.7)

7

8/17/2019 FullText (40)

28/142

and

M̃ m =

1 − R̃m,m−1 R̃m,m+1e2ikmz(dm+1−dm)−1

(2.8)

B. n > m

The observation layer is above the source layer. There is only secondaryfield. Then

F (z, z ′) =

eikmz(dm+1−z′) + R̃m,m−1e

ikmz(dm+1−2dm+z′)

·

eiknz(z−dn) + R̃n,n+1eiknz(2dn+1−dn−z)

M̃ m T̃ mn

(2.9)

where

T̃ mn =n−1

j=m+1 eikjz (dj+1−dj)S j−1,j

S n−1,n (2.10)

S j−1,j = T j−1,j

1 − R j,j−1 R̃ j,j+1eikjz (2dj+1−2dj). (2.11)

C. n < m

The observation layer is below the source layer. There is also only sec-

ondary field. Then

F (z, z ′) =

eikmz(z′−dm) + R̃m,m+1e

ikmz(2dm+1−dm−z′)

· eiknz(dn+1−z) + R̃n,n−1eiknz(z+dn+1−2dn) M̃ m T̃ mn(2.12)

where

T̃ mn =m−1

j=n+1

eikjz (dj+1−dj)S j+1,j

S n+1,n (2.13)

S j+1,j = T j+1,j

1−R j,j+1 R̃ j,j−1eikjz (2dj+1−2dj) (2.14)

2.2.2 Matrix-Friendly Formulation

The matrix-friendly formulation can be derived in the implementation of

method of moments (MoM), where the spatial derivatives can be moved from

the Green’s function to either testing or basis functions. In this manner,

the singularity of the LMGF can be reduced and the Sommerfeld integral

converges more rapidly. Here we only summarize the basic expressions and

8

8/17/2019 FullText (40)

29/142

the detailed derivation can be found in [21]. In the electric field integral

equation (EFIE) formulation, there are five scalar Green’s functions in the

matrix-friendly formulation.

gss(r, r′) = k2ρgTE(r, r′) (2.15)

gzz(r, r′) = k2mng

TM(r, r′)− ∂ z∂ z′gTE(r, r′) (2.16)

gz1(r, r′) =

µnµm

∂ z′gTM(r, r′) + ∂ zg

TE(r, r′) (2.17)

gz2(r, r′) =

ǫmǫn

∂ zgTM(r, r′) + ∂ z′g

TE(r, r′) (2.18)

gφ(r, r′) =

∂ z∂ z′

k2nmgTM(r, r′)− gTE(r, r′) (2.19)

where the partial derivative with respect to z and z ′ can be easily imple-

mented in the spectral domain, namely, ∂ z = ±iknz and ∂ z′ = ±ikmz wherethe signs are determined by the relative positions of the source and observa-

tion points.

2.3 Discrete Complex Image Method (DCIM)

2.3.1 Traditional DCIM

The Green’s functions in (2.15)–(2.19) can be expressed as an infinite integral

of the following type

g(ρ) = i

8π

+∞−∞

dkρkρkz

H (1)0 (kρρ)g̃(kρ). (2.20)

If the integration kernel can be approximated by a series of complex expo-

nentials,

g̃(kρ) =M i=1

aieikzbi (2.21)

by applying the Sommerfeld identity [5],

eikr

r =

i

2

+∞−∞

dkρkρkz

H (1)0 (kρρ)e

ikzz, r =

ρ2 + z 2 (2.22)

9

8/17/2019 FullText (40)

30/142

the infinite integral can be calculated in a closed form,

g(ρ) =M i=1

aieikri

4πri, ri =

ρ2 + b2i . (2.23)

The complex exponential series can be obtained by, for example, the matrix

pencil method [38], which approximates a function with real argument by

y(t) =M i=1

RieS it. (2.24)

The mapping of the the real variable t to kz is given by [35],

kz = k

it +

1− t

T 0

, 0 ≤ t ≤ T 0 (2.25)

where ai and bi in (2.23) can be obtained by

bi = iS iT 0

k(1 − iT 0) , ai = Rie−ikbi (2.26)

The far field prediction of the traditional DCIM is poor, and various remedies

have been proposed [35]–[37], [40].

2.3.2 DCIM Based on the Sommerfeld Branch Cut

When the transverse distance is large, we deform the integration path to the

Sommerfeld branch cut.

One Branch-Cut Case

If the layered medium is backed by a PEC ground plane, there is only one

branch cut associated with the top layer [5]. A typical microstrip structure

falls into this case, as is shown in Figure 2.2. Due to the Cauchy’s theorem

and Jordan’s lemma, the integral defined along the Sommerfeld integration

path (SIP) is equivalent to the path integral along the Sommerfeld branch

cut (SBC) and some pole contributions, as shown in Figure 2.3.

Based on the deformed path, the Green’s function can be expressed as a

10

8/17/2019 FullText (40)

31/142

Figure 2.2: A typical microstrip structure (unit: mm).

Figure 2.3: One branch-cut case: Sommerfeld integration path (SIP),Sommerfeld branch cut (SBC) and poles in the complex kρ plane.

11

8/17/2019 FullText (40)

32/142

superposition of the following two terms,

g = gbranch + gpole (2.27)

wheregbranch =

i

8π

SBC

dkρkρ

kNzH

(1)0 (kρρ)g̃(kρ) (2.28)

gpole = −14

q

kρ,qkNz,q

H (1)0 (kρ,qρ)Res[g̃(kρ,q)] (2.29)

where q is the number of poles and Res [g̃(kρ,q)] is the residue of the kernel.

The locations of the poles and relevant residues can be obtained by a ro-

bust pole-searching algorithm [43]. Pole contributions are represented by the

Hankel function, which has the following asymptotic behavior,

H (1)0 (kρρ) ∼

2

πkρρei(kρρ−

π4 ) (kρρ →∞) (2.30)

If kρ is real, we have

gpole ∼

1/ρ (2.31)

The gbranch represents the radiation modes from the branch cut integration,

which can be obtained in closed form by the new DCIM. By transforming

the variable from kρ to kNz, (2.28) becomes

gbranch = i

8π

+∞−∞

dkNzH (1)0 (kρρ)g̃(kρ) (2.32)

with

dkNz = − kρkNz

dkρ (2.33)

In order to use (2.24) to approximate the kernel of (2.32), we can let

kNz = t − T 02

, 0 ≤ t ≤ T 0. (2.34)

Then, ai and bi can be obtained similarly from S i and Ri,

bi = S i

i , ai = Rie

ibiT 0/2 (2.35)

Once it is approximated by complex exponentials, we can apply the alterna-

12

8/17/2019 FullText (40)

33/142

Figure 2.4: Two branch-cut case: Sommerfeld integration path (SIP),Sommerfeld branch cuts (SBC) and poles in the complex kρ plane.

tive Sommerfeld identity in terms of kz integration [5] to get a closed-form

solution of gbranch,

eikr

r =

i

2

+∞−∞

dkzH (1)0 (kρρ)e

ikzz, r =

ρ2 + z 2 (2.36)

The radiation modes include spatial wave and the lateral wave, which have

the following asymptotic behavior respectively

gspatial ∼ 1/ρ (2.37)

glateral ∼ 1/ρ2 (2.38)From (2.31), (2.37), and (2.38), we can see that usually the surface wave

dominates in the far field. However, at the interface, the spatial waves of the

primary term and the secondary term cancel each other and the lateral wave

can be observed, if there are no pole contribution for certain cases. This

cancelation in DCIM was first analyzed in details in [40].

Two Branch-Cut Case

For a general layered medium, there are two branch cuts associated with the

top and bottom layers, where radiation modes can be supported. In this

case, the path and possible poles are shown in Figure 2.4.

13

8/17/2019 FullText (40)

34/142

In this case, (2.20) becomes

g = gbranch,N + gbranch,1 + gpole (2.39)

where gbranch,N and gpole are similar to those in (2.32) and (2.29), whilegbranch,1 has the form of

gbranch,1 = i

8π

+∞−∞

dk1zH (1)0 (kρρ)

k1zkNz

g̃(kρ). (2.40)

2.4 Numerical Results

Several numerical results are presented in this section. Only gss and gφ are

used as illustrative examples here. The microstrip shown in Figure 2.2 is

first studied. The working frequency is f = 3 GHz, and the source point

and observation point are at the interface between the air and dielectric

substrate. We apply the traditional DCIM for small ρ, and switch it to the

new implementation when ρ is large. The transition region can be set in

100 < k0ρ < 101. In the following examples, we set the transition point at

k0ρ = 100.5. For this microstrip problem, only one real TM pole is found. The

gss and gφ are calculated in Figure 2.5 and Figure 2.6; both agree well with

those from numerical integration. In Figure 2.6, since there is no TE pole,

we can observe that the asymptotic behavior of gss is 1/ρ2, which represents

the lateral wave. It agrees with the results by the three-level DCIM in [40]

except for a constant due to the definition of the Green’s function. It is also

reported in [40] that the popular two-level DCIM cannot correctly capture

the branch point contribution in this case. In Figure 2.6, since gφ contains

both TE and TM waves, the pole contribution dominates in the far field,

which is of 1/√

ρ.

The acceleration of DCIM lies in the fact that multiple source–observationpoints can share the same image coefficients. For the numerical examples

shown in Figure 2.5 and Figure 2.6, we compare the computational time

for calculating 5,000 k0ρ samples. The central processing unit (CPU) time

is listed in Table 2.1. We can observe that the computation can be much

accelerated by using DCIM.

To validate the two-branch cut case, a 3-layer model with lossy material

14

8/17/2019 FullText (40)

35/142

−3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

log10

(k0ρ)

l o g 1 0 | g

s s

|

Numerical

DCIM

Figure 2.5: The magnitude of gss versus k0ρ for the microstrip structure.The asymptotic behavior is 1/ρ2, which agrees with the theoreticalprediction.

Table 2.1: CPU time comparison (seconds)

LMGF Numerical DCIM

gss 56.82 0.59

gφ 74.69 2.44

shown in Figure 2.7 is studied. The working frequency is f = 1.5 GHz.

Both gss and gφ are calculated and compared with the numerical integration

results, as are shown in Figure 2.8 and Figure 2.9. Again, good agreement

can be observed. In this case, the poles are general complex numbers away

from the real axis, so their contributions decay quickly when ρ is large. The

lateral wave dominates in the far field with the asymptotic behavior of 1 /ρ2,

as are shown in Figure 2.8 and Figure 2.9.

2.5 Summary

A novel implementation of the DCIM based on Sommerfeld branch cut is

proposed to improve the far field prediction of the LMGF. By contour defor-

mation, the Green’s function can be naturally decomposed into the radiation

modes and guided modes. The guided modes can be obtained by a robust

pole-searching algorithm and the radiation modes can be calculated in a

15

8/17/2019 FullText (40)

36/142

−3 −2 −1 0 1 2 3 4

−6

−5

−4

−3

−2

−1

0

log10

(k0ρ)

l o g 1 0

| g φ

|

Numerical

DCIM

Figure 2.6: The magnitude of gφ versus k0ρ for the microstrip structure.The far field is dominated by the pole contribution, which has theasymptotic behavior of 1/

√ ρ.

Figure 2.7: A 3-layer model (unit: mm).

16

8/17/2019 FullText (40)

37/142

−3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

log10(k0ρ)

l o g 1 0

| g s s

|

NumericalDCIM

Figure 2.8: The magnitude of gss versus k0ρ for the 3-layer model. Theasymptotic behavior is 1/ρ2.

−3 −2 −1 0 1 2 3 4−12

−10

−8

−6

−4

−2

0

2

log10

(k0ρ)

l o g 1 0

| g φ

|

Numerical

DCIM

Figure 2.9: The magnitude of gφ versus k0ρ for the 3-layer model. Theasymptotic behavior is 1/ρ2.

17

8/17/2019 FullText (40)

38/142

closed form, so that the evaluation can be made efficient compared to the

direct numerical integration. For small ρ interaction, we simply switch back

to the traditional DCIM to capture the near field. One should note that

in this new implementation, when ρ becomes small, the length of sampling

path in spectral domain increases, and it becomes harder to approximate the

kernel. Efforts can be made to improve this DCIM in the near field, such

as extracting the asymptotic behavior analytically. At the same time, for

cases when poles are very close to the branch cut, the accuracy of function

approximation may be affected and more careful treatment of the poles is

necessary. Such efforts shall also be carried out in the future to improve this

DCIM.

18

8/17/2019 FullText (40)

39/142

CHAPTER 3

MIXED-FORM THIN-STRATIFIED

MEDIUM FAST-MULTIPOLEALGORITHM (MF-TSM-FMA)

A mixed-form thin-stratified medium fast-multipole algorithm (MF-TSM-

FMA) is proposed for fast simulation of general microstrip structures at both

low and mid-frequencies. The newly developed matrix-friendly formulation

of layered medium Green’s function (LMGF) is applied in this algorithm.

For well-separated interactions, the contour deformation technique is imple-mented to achieve a smoother and exponentially convergent integral. The

two-dimensional addition theorem is then incorporated into the integrand to

expedite the matrix-vector product. In our approach, multipole expansion

(low-frequency fast-multipole algorithm) as well as the plane wave expansion

(mid-frequency fast-multipole algorithm) of the translational addition theo-

rem are combined into a single multilevel tree to capture quasi-static physics

and wave physics simultaneously. The outgoing wave is represented first in

terms of multipole expansion at leafy levels, and then switched to plane wave

expansion automatically at higher levels. This seamless connection makes the

algorithm applicable in simulations, where subwavelength interaction (circuit

physics) and wave physics both exist.

3.1 Introduction

The simulation of microstrip structures has attracted intensive study for

many years. Since this kind of structure usually involves a thin layeredmedium, the LMGF can be applied to reduce the number of unknowns. All

the acceleration techniques introduced at the beginning in Chapter 2 can be

implemented to expedite the process during the matrix filling stage in the in-

tegral equation method. However, since the integral equation method leads

to a full matrix, the memory consumption is large and the matrix-vector

product is still time-consuming in iterative solvers. Fast algorithms are in-

19

8/17/2019 FullText (40)

40/142

dispensable for large-scale computation. In the last two decades, various

fast algorithms were proposed, first based on the free space Green’s function,

and then extended to other kernels. The incomplete list includes the fast

multipole method (FMM) [44], [45] and the multilevel fast multipole algo-

rithm (MLFMA) [46]–[48], the adaptive integral method (AIM) [49], [50],

the conjugate gradient fast Fourier transform (CG-FFT) [51], [52] and the

precorrected FFT method [53]–[55], the multilevel Green’s function inter-

polation method [56], the kernel-independent or black-box FMM [57], [58],

the multilevel matrix decomposition algorithm (MLMDA) [59], the adaptive

cross approximation (ACA) algorithm [60], [61], the IES3 [62], and the IE-

QR algorithm [63] etc. These algorithms can roughly be categorized into two

classes, the algebraic algorithms, which are kernel independent and are thus

more flexible, and the physics-based algorithms, which utilize the physicalproperty of the concrete kernel and are thus usually more efficient, especially

for dynamic field. In this chapter, we only focus on one of the physical al-

gorithms and the applications in microstrip structures. A detailed survey of

the algebraic algorithms can be found in [60].

For the physical algorithm, both MLFMA- and FFT-based methods have

been extended to layered medium problems, for instance, in the parasitic

extraction and simulation of microstrip antennas and arrays. The former one

factorizes the Green’s function by using the translational addition theorem

in a multilevel manner [47], [48], [64]–[67] and the latter one takes advantage

of the transverse translational invariance property of the interactions and

applies FFT to accelerate the computation in the Fourier domain [50], [52],

[55], [68], [84]. Both these methods can achieve an O(N log N ) computational

complexity in the dynamic regime.

MLFMA is more efficient than FFT if the objects are sparse, however,

MLFMA usually suffers from a low frequency breakdown problem [70], which

means that the plane wave expansion of the addition theorem, attractive at

mid-frequencies, is error uncontrollable due to the finite precision of the com-

puter. The multipole expansion is usually implemented instead in the low

frequency regime. This expansion is suitable for scale-invariant interaction,

but becomes inefficient rapidly when entering the electrodynamic regime. For

many broadband large-scale simulation, both electrodynamic interaction and

sub-wavelength interaction exist in one problem. This requires a broadband

fast algorithm capable of capturing the wave physics and quasi-static physics

20

8/17/2019 FullText (40)

41/142

simultaneously. Motivated by the idea of mixed-form fast multipole algo-

rithm in free space [70], a mixed-form thin-stratified medium fast-multipole

algorithm (MF-TSM-FMA), based on the newly-formulated LMGF [21], is

developed in this chapter. In this MF-TSM-FMA, the multipole expansion

and the plane wave expansion are combined into one multilevel tree, where

different scales of interaction can be separated by the multilevel nature of

the MLFMA [71]. A transition equation can be derived to connect these two

expansions.

3.2 Review of Thin-Stratified Medium Fast-Multipole

Algorithm (TSM-FMA)

3.2.1 TSM-FMA Based on the Plane Wave Expansion

We show the LMGF again here [21],

Ḡ(r, r′) = (∇× ẑ )(∇′ × ẑ )gTE(r, r′)

+ 1k2nm

(∇×∇× ẑ )(∇′ ×∇′ × ẑ )gTM(r, r′)(3.1)

where

gTE/TM(r, r′) = i

8π

+∞−∞

dkρkmzkρ

H (1)0 (kρρ)F

TE/TM(kρ, z , z ′) (3.2)

To obtain a TSM-FMA, the path deformation technique is necessary to im-

prove the convergence of the Sommerfeld integral.

Deformation of the Integration Path

For well separated interactions, exponential convergence can be achieved by

deforming the Sommerfeld integration path into the vertical branch cut for

a thin layered medium. In this structure, the vertical dimension is tightly

confined compared to the horizontal dimension and the source-observation

dependent steepest descent path for a thick layered medium evolves to be

source-observation independent [5], letting us to implement the fast algo-

rithm. Possible pole singularity should be included in such deformation ac-

21

8/17/2019 FullText (40)

42/142

cording to Cauchy’s theorem [47]. The detailed implementation about the

path integration, and the choice of Riemann sheets, can be found in [67] and

will not be repeated here.

Two Dimensional MLFMA

The transverse interaction of the LMGF is the Hankel function, so 2D MLFMA

can be implemented to accelerate the computation. The core equations of

the MLFMA are listed below [72],

H (1)0 (kρρ ji) = 1

2π

2π0

dαβ̃ jJ (α)α̃JI (α)β̃ Ii(α) (3.3)

α̃JI (α) =P

p=−P

H (1) p (kρρJI )e−ip(φJI −α−

π2) (3.4)

β̃ jJ (α) = eikρ·ρjJ β̃ Ii(α) = e

ikρ·ρIi (3.5)

where J and I are box centers for the observation point j and the source

point i, respectively. Here kρ is in general complex. In implementing TSM-

FMA, the z -dependent propagation factor F α(z, z ′) can be easily factorized

and embedded into the radiation and receiving patterns [67].

3.2.2 Low Frequency Breakdown of the TSM-FMA

When the operation frequency decreases, the plane wave expansion becomes

unstable due to the numerical overflow of the outgoing-to-outgoing translator

α̃. It has been discussed for free space MLFMA (spherical harmonics) by

[70]. We will show that this is also the case for layered medium problems

(cylindrical harmonics) for general complex argument kρρ. The asymptotic

expression of the Hankel function is

H (1)0 (kρρ) ∼ Ceik0ρe−u

2ρ (3.6)

where kρ = k0 + iu2, eik0ρ is the propagating wave while e−u

2ρ is the decaying

factor. For a fixed quadrature point u, we can show that numerical insta-

bility of the plane wave expansion exists in the low frequency regime. In

implementing the addition theorem, the worst interaction situation is shown

22

8/17/2019 FullText (40)

43/142

Figure 3.1: The worst interaction in applying the addition theorem.

Figure 3.2: Relative error by using plane wave expansion.

in Figure 3.1. For a given quadrature u, if we fix the frequency, then the

imaginary part of the argument may introduce different decays for different

ρ, making the Hankel function become negligible rapidly for large ρ. To

make a fair comparison, we may fix the transverse distance ρ and change the

frequency. Let the box size to be b = 5 mm and u = 4. The relative errors of

the plane wave expansion with the expansion order P for different box sizes

are shown in Figure 3.2. It is shown that the expansion is unstable when

the box size is much smaller than the wavelength. Since in near region, the

evanescent wave dominates and the field is scale invariant, the plane wave

expansion based on wave physics is no longer suitable.

23

8/17/2019 FullText (40)

44/142

3.3 Mixed-Form Thin-Stratified Medium

Fast-Multipole Algorithm (MF-TSM-FMA)

3.3.1 TSM-FMA Based on the Multipole ExpansionThe multipole expansion, which uses the undiagonalized addition theorem,

is free of low-frequency breakdown since normalization can be made and

different translators (α̃ and β̃ ) are well balanced. However, since the ∇operator cannot be calculated analytically in the multipole expansion, we

will apply the matrix-friendly formulation of the LMGF.

Matrix-Friendly LMGF

The matrix element in the MoM procedure can be expressed in the matrix-

friendly LMGF as [21]

Z ji = Z ss ji + Z

zz ji + Z

z1 ji + Z

z2 ji + Z

φ ji (3.7)

where

Z ss ji = iωµmf s(r j), gss(r j, r′i), f s(r′i) (3.8)

Z zz

ji = iωµmẑ · f (r j), gzz(r j, r′

i), ẑ · f (r′

i) (3.9)Z z1 ji = −iωµmẑ · f (r j), gz1(r j , r′i),∇′ · f (r′i) (3.10)

Z z2 ji = −iωµm∇ · f (r j), gz2(r j , r′i), ẑ · f (r′i) (3.11)Z φ ji = iωµm∇ · f (r j), gφ(r j , r′i),∇′ · f (r′i) (3.12)

gss(r j, r′i) = k

2ρg

TE (r j, r′i) (3.13)

gzz(r j, r′

i) = k2

mngTM

(r j, r′

i)− ∂ z∂ z′gTE

(r j, r′

i) (3.14)

gz1(r j, r′i) =

µnµm

∂ z′gTM (r j, r

′i) + ∂ zg

TE (r j, r′i) (3.15)

gz2(r j , r′i) =

εmεn

∂ zgTM (r j, r

′i) + ∂ z′g

TE (r j , r′i) (3.16)

gφ(r j, r′i) =

∂ z∂ z′

k2nmgTM (r j, r

′i)− gTE (r j , r′i) (3.17)

24

8/17/2019 FullText (40)

45/142

Compared with (3.1), the matrix-friendly Green’s function has lower-order

singularity, since several spatial derivatives have been moved to the basis

function.

Factorization of the Propagation Factor

The factorization of the propagation factor is critical in implementing a gen-

eralized TSM-FMA, because the 2D MLFMA can only factorize the trans-

verse interaction, and the vertical interaction should be factorized manually

to bind to the radiation and receiving patterns. The propagation factor takes

different forms if the positions of source and observation points are different.

According to the expressions shown in Chapter 2, the propagation factor can

be easily factorized as followings,A. n = m

F (z, z ′) = F 1(z, z ′) + f v2(z )I 2(m)f r2(z

′)

+f v3(z )I 3(m)f r3(z ′)

(3.18)

where

F 1(z, z ′) = eikmz|z−z

′| (3.19)

f v2(z ) = eikmz(dm+1−z) (3.20)

f r2(z ′) = [eikmz(dm+1

−z

′

) + R̃m,m−1eikmz(dm+1

−2dm+z

′

)] (3.21)

I 2(m) = M̃ m R̃m,m+1 (3.22)

f v3(z ) = eikmz(z−dm) (3.23)

f r3(z ′) = [eikmz(z

′−dm) + R̃m,m+1eikmz(2dm+1−dm−z′)] (3.24)

I 3(m) = M̃ m R̃m,m−1 (3.25)

B. n > m

F (z, z ′

) = f v(z )I (m, n)f r(z ′

) (3.26)

where

f v(z ) = [eiknz(z−dn) + R̃n,n+1e

iknz(2dn+1−dn−z)] (3.27)

f r(z ′) = [eikmz(dm+1−z

′) + R̃m,m−1eikmz(dm+1−2dm+z′)] (3.28)

I (m, n) = M̃ m T̃ mn (3.29)

25

8/17/2019 FullText (40)

46/142

C. n < m

F (z, z ′) = f v(z )I (m, n)f r(z ′) (3.30)

where

f v(z ) = [e

iknz(dn+1−z)

+ ˜Rn,n−1e

iknz(z+dn+1−2dn)

] (3.31)f r(z

′) = [eikmz(z′−dm) + R̃m,m+1e

ikmz(2dm+1−dm−z′)] (3.32)

I (m, n) = M̃ m T̃ mn (3.33)

Here I (m, n) only depends on layer parameters. When n = m, the modulus

sign in the direct term F 1(z, z ′) = eikmz|z−z

′| is a problem, however, we can

show that this modulus sign can be removed in the context of vertical branch

cut integration [67]. After implementing factorization of the propagation

factor, the addition theorem can be applied to accelerate the computation.However, there are totally ten scalar components that need to be factorized

separately, hence optimization of the patterns are necessary.

Layout of Radiation and Receiving Patterns

For a certain source-observation-layer pair, the partial derivative ∂ z or ∂ z′ in

the matrix-friendly formulation (3.13)–(3.17) is acting directly on the propa-

gation factor. It can be calculated easily according to the concrete expression

of F (z, z ′). For example, if n > m,

∂ z∂ z′F (z, z ′)

= ∂ zf (z )∂ z′f (z ′)I (m, n)

= ikmz

−eikmz(dm+1−z′) + R̃m,m−1eikmz(dm+1−2dm+z′)

·iknz

eiknz(z−dn) − R̃n,n+1eiknz(2dn+1−dn−z)

T̃ mn M̃ m

(3.34)

By investigating the matrix-friendly formulas (3.7)–(3.17), we know that

there are eight scalar components and one 2D vector component requiring

factorization separately. The memory requirement is huge if we store all

the radiation and receiving patterns directly. In our implementation, the

radiation and receiving patterns are filled on the fly in the process of matrix-

vector product. In order to reduce the number of MLFMA processes, we can

set two receiving patterns for one radiation pattern. For example, The TE

26

8/17/2019 FullText (40)

47/142

Table 3.1: Layout of radiation and receiving patterns

Radiation Receiving A Receiving B

gTE (r′)x̂ · f (r′) x̂ · f (r)gTE (r) 0

gTE

(r′

)ŷ · f (r′

) ŷ · f (r)gTE

(r) 0gTE (r′)∇′ · f (r′) ∇ · f (r)gTE (r) ẑ · f (r)∂ zgTE (r)

∂ z′gTE (r′)ẑ · f (r′) ∇ · f (r)gTE (r) ẑ · f (r)∂ zgTE (r)

gTM (r′)ẑ · f (r′) ẑ · f (r)gTM (r) ∇ · f (r)∂ zgTM (r)∂ z′g

TM (r′)∇′ · f (r′) ẑ · f (r)gTM (r) ∇ · f (r)∂ zgTM (r)

part of Z φ ji and Z z1 ji can be accelerated at the same time. If we denote

g

α

(r, r′

) = g

α

(r)g

α

(r′

) (3.35)

then we can set one radiation pattern and two receiving patterns to be

I r = ∇′ · f (r′)gTE (r′) (3.36)

I va = ∇ · f (r)gTE (r) (3.37)

I vb = ẑ · f (r)∂ zgTE (r) (3.38)

By doing this, two interactions can share the same radiation patterns,outgoing waves and incoming waves. However, minor modification is needed

in the last stage, where incoming waves at the leafy level are disaggregated

into observation points by using different receiving patterns. The pattern

layout is listed in Table 3.1. By categorizing the patterns, one may only

need to calculate six kinds of interactions instead of ten.

Two Dimensional Low Frequency MLFMA

In the low frequency regime, multipole expansion is applied to factorize the

Hankel function. The general core equations are listed in the following [72]

αmn(ρ ji) =M

N

β mM (ρ jJ )αMN (ρJI )β Nn(ρIi) (3.39)

αmn(ρ ji) = H (1)m−n(kρρ ji)e

−i(m−n)φji (3.40)

27

8/17/2019 FullText (40)

48/142

Figure 3.3: Relative error by using multipole expansion.

β mn(ρ ji) = J m−n(kρρ ji)e−i(m−n)φji (3.41)

For Hankel function, we have

H (1)0 (kρρ ji) =M N

β 0M (ρ jJ )αMN (ρJI )β N 0(ρIi) (3.42)

Though this multipole expansion is equivalent to the plane wave expan-

sion mathematically, the former is stable in the low frequency regime since

normalization can be easily implemented to balance the α and β , to make

each numerical step error controllable. For the same setup as in Figure 3.1,

the relative error of the multipole expansion is shown in Figure 3.3. It is

obvious that the expansion is stable at low-frequencies. However, when the

frequency increases, the multipole expansion becomes very inefficient since

more and more expansion terms are required to capture the wave physics.So it should be switched back to plane wave representation to describe the

dynamic field.

28

8/17/2019 FullText (40)

49/142

3.3.2 Mixed-Form Expansion

For many applications such as broadband microstrip antennas or arrays, the

interaction may involve both quasi-static field and dynamic field. Neither of

the two expansions can be used independently to capture the two different

fields. Following the idea of [70], one can combine the two expansions into

a multilevel tree and switch between them dynamically. Based on addition

theorem and integral representation of the Bessel function, we can derive the

following mixed-form core equation

αmn(ρ ji) = 1

2π

2π0

dαe−i(α+π2)mβ̃ jJ (α)α̃JI (α)β̃ Ii(α)e

i(α+π2)n (3.43)

This equation combines the multiple expansion and plane wave expansion.

Substituting it into (3.39), we can construct a mixed-form TSM-FMA. At

leafy levels, the box size is much smaller compared to the wavelength, where

interaction is in the sub-wavelength regime, the outgoing wave is calculated

by multipole expansion. With the process of aggregation, the box size be-

comes large enough, where the interaction may enter the dynamic regime, the

outgoing wave is switched into the plane wave representation to capture the

wave physics. Similar idea can be applied in the process of disaggregation.

In matrix notation, the MF-TSM-FMA for a certain kρ is expressed as,

H (1)0 (kρρ ji) = [β jJ 1]1×P 1 · [β J 1J 2]P 1×P 2 · [T ]†P 2×K 3 · [β̃ J 2J 3]K 3×K 3

· [I ]T K 3×K 4 · [β̃ J 3J 4]K 4×K 4 · [α̃J 4I 4]K 4×K 4 · [β̃ I 4I 3]K 4×K 4

· [I ]K 4×K 3 · [β̃ I 3I 2]K 3×K 3 · [T ]K 3×P 2 · [β I 2I 1]P 2×P 1 · [β I 1i]P 1×1(3.44)

where [I ] is the interpolation matrix in the dynamic regime, and [T ] trans-

forms the multipole outgoing wave into plane wave outgoing wave.

T =

ei(α+π2)nK ×P

(3.45)

The implementation scheme is shown in Figure 3.4. We denote the level

where I 2 in the figure resides as the switch level. From our experience, the

switch level can be chosen with box size to be around Re [kρ] b = 0.2π −0.4π, namely 0.1λ0 − 0.2λ0 for Re[kρ] = k0, where λ0 is the wavelength infree space. If the switch box size is smaller than the box size of the leafy

29

8/17/2019 FullText (40)

50/142

Figure 3.4: The MF-TSM-FMA structure.

level, all the interaction are calculated in plane wave expansion and the MF-

TSM-FMA goes back to the mid-frequency TSM-FMA.

For dynamic field, the computational complexity is O(N log N ), while for

quasi-static field, the computational complexity is O(N ). The complexity

of MF-TSM-FMA should stay in between depending on the location of the

switch level. However, the real computational time may depend on other

factors. For example, the translation matrix in multipole expansion is full,

which may take longer time even though the complexity is lower for a given

problem.

3.4 Numerical Results

In this section, several numerical results are demonstrated to validate our

algorithm. We have first simulated an 8 × 4 corporate-fed microstrip arrayavailable in the literature [48], [50], [52], [54]. The geometrical structure of

this array is shown in Figure 3.5. The detailed parameters are l = 10.08 mm,

ω = 11.79 mm, d1 = 1.3 mm, d2 = 3.93 mm, l1 = 12.32 mm, l2 = 18.48

mm, D1 = 23.58 mm, D2 = 22.40 mm. The thickness of the substrate is

30

8/17/2019 FullText (40)

51/142

d = 1.59 mm, and the relative permittivity and permeability are ǫr = 2.2

and µr = 1. The radiation patterns are calculated by feeding at the input

port at the working frequency of f = 9.42 GHz. The patterns in both E-

plane (φ = 0o) and H-plane (φ = 90o) are shown in Figure 3.6. The results

agree well with the results extracted from [36]. Here, the box size at the

leafy level is 0.3λ0, and a five-level FMA is set up. Since the switch level is

below the leafy level, the MF-TSM-FMA applies the plane-wave expansion

for all levels. In order to activate the low-frequency FMA, we decrease the

frequency to f = 2.5 GHz, and calculate the radar cross section (RCS) of

this microstrip array. The results are shown in Figure 3.7. The leafy level

box size is 0.08λ0, and the switch level is 4. The number of unknown is 9,394,

and the memory consumption of the MF-TSM-FMA is 40 MB. Next, to show

the capability of modeling vertical structures, a dipole antenna mounted ona 4× 4 mushroom-type electromagnetic band gap (EBG) is simulated, withits geometrical structure shown in Figure 3.8. The geometrical parameters

are: w = 29 mm, g = 1 mm, d = 1 mm, h = 2 mm, L = 63 mm, W = 1

mm, H = 2 mm. The FR4 is used as the substrate with ǫ = 4.4 and

tan δ = 0.02. The antenna is working at f = 2.2 GHz. The radiation pattern

in the E plane is calculated and shown in Figure 3.9 (in unit of dB). Due

to the layered medium assumption, there is no back lobe and side lobes in

the radiation pattern. A 30×

30 microstrip array [47] shown in Figure 3.10

is then simulated at frequency of 2.2 GHz, 1.1 GHz, and 550 MHz, where

a = b = 6 cm, L = 3.66 cm, W = 2.6 cm, the thickness of the substrate is

d = 0.158 cm, ǫr = 2.17. The number of unknown is 117,000. A seven-level

MF-TSM-FMA is constructed. At the high frequency end (f = 2.2 GHz),

the RCS is compared with the one sampled from [47] and good agreement is

achieved. At the low frequency end (f = 550 MHz), the leafy level box size is

0.03 λ0, and the switch level is 5, which corresponds to a box size of 0.12 λ0.

The memory requirement is 250 MB, the central processing unit (CPU) time

for each iteration is 2.2 min. Finally, we list the memory consumption and

CPU time per iteration versus the number of unknowns in Figure 3.12, where

the proportions (the number of levels) of multipole expansion and plane wave

expansion are set to be similar. All the examples run on a personal computer

(PC) with Intel 2.66 GHz processor. It is obvious that this algorithm can

be used in the simulation of large scale microstrip structures with acceptable

computational cost.

31

8/17/2019 FullText (40)

52/142

8/17/2019 FullText (40)

53/142

0 50 100 150 200 250 300 350−80

−75

−70

−65

−60

−55

−50

−45

−40

−35

φ (°)

B i s t a t i c R C S ( d B s m )

MoM

MF−TSM−FMA

Figure 3.7: The bistatic RCS of the microstrip array at f = 2.5 GHz and(θi, φi) = (60

◦, 0◦).

Figure 3.8: A dipole on an EBG structure, top and side view.

33

8/17/2019 FullText (40)

54/142

Figure 3.9: Radiation pattern of the dipole in E plane, due to the layeredmedium assumption, the ground plane is infinitely large.

Figure 3.10: The geometrical model of a 30 × 30 microstrip array.

34

8/17/2019 FullText (40)

55/142

0 10 20 30 40 50 60 70 80 90

−160

−140

−120

−100

−80

−60

−40

θ (°)

B i s t a t i c R C S ( d B s m )

f=2.2 GHz

Ref.

f=1.1 GHz

f=0.55 GHz

Figure 3.11: Bistatic RCS of the 30 × 30 microstrip array.

0 2 4 6 8 10 12

x 104

0

100

200

300

M e m o r y r e q u i r e m e n t ( M B )

0 2 4 6 8 10 12

x 104

0

1

2

3

T i m e p e r i t e r a t i o n ( m i n )

Figure 3.12: Memory requirement (solid line) and CPU time per iteration(dash line) versus number of unknowns.

35

8/17/2019 FullText (40)

56/142

3.5 Summary

A mixed-form thin-stratified medium fast multipole algorithm is developed in

this chapter. By path deformation, the LMGF can be expressed by several-

term summation of Hankel functions weighted by z -dependent propagationfactors. For each quadrature point, the Hankel function can be accelerated

by 2D MLFMA. The matrix-friendly formulation of the LMGF is applied,

and interaction is categorized into six components if we introduce two re-

ceiving patterns for one radiation pattern. Due to the numerical instability

of the TSM-FMA in the low frequency regime, a mixed-form TSM-FMA is

proposed for broadband simulation. By utilizing the multilevel nature of

the MLFMA, interactions with different scales are accelerated by different

FMA expansions, namely multipole expansion and plane wave expansion. A

transform equation is derived to embed the two expansions into one MLFMA

tree. Numerical results show the accuracy and efficiency of this algorithm. It

should be noted that in the very low frequency regime, the integral equation

in layered medium may suffer from another kind of low frequency breakdown,

where the existence of the null space of the operator plagues the numerical

procedure. Before modeling large problems in extremely low frequency, one

must develop an efficient algorithm accounting for the two breakdowns at

the same time.

36

8/17/2019 FullText (40)

57/142

CHAPTER 4

REMEDIES FOR LOW FREQUENCY

BREAKDOWN OF INTEGRALOPERATOR IN LAYERED MEDIUM

When the working frequency further decreases, the low frequency breakdown

of the integral operator—the electric field integral equation (EFIE) occurs.

Two remedies—the loop-tree decomposition and the augmented electric field

integral equation (A-EFIE)—which were first developed in free space, are

extended to layered medium in this chapter. In the loop-tree decomposition,the current is decomposed into divergence-free part and non-divergence-free

part according to quasi-Helmholtz decomposition when frequency tends to

zero, in order to capture both capacitance and inductance physics. Fre-

quency normalization and basis rearrangement are applied to stabilize the

matrix system. In the A-EFIE, the traditional EFIE can be cast into a

generalized saddle-point system, by separating charge as extra unknown list

and enforcing the current continuity equation. Frequency scaling for the

matrix-friendly layered medium Green’s function (LMGF) is analyzed when

frequency tends to zero. Rank deficiency and the charge neutrality enforce-

ment of the A-EFIE for LMGF is discussed in detail. The electrostatic limit

of the A-EFIE is also analyzed. Without any topological loop-searching algo-

rithm, electrically small conducting structures embedded in a general layered

medium can be simulated by using this new A-EFIE formulation.

4.1 Introduction

Computational electromagnetics becomes indispensable as a computer-aided

design (CAD) methodology in various electrical engineering applications,

such as in integrated circuit and wireless communication device. The op-

erating frequency of the electrical systems keeps on increasing to several

gigahertz, meanwhile fabrication process has achieved nanoscale. Hence, a

broadband simulation tool is badly needed for capturing circuit physics of

37

8/17/2019 FullText (40)

58/142

the tiny structures as well as wave physics for the whole package. Unfor-

tunately, however, in the sub-wavelength regime, or so-called low frequency

regime, the electric field and magnetic field decouple, and the total current

in Maxwell’s equations decomposes into a divergence-free part and a curl-

free part, with different frequency scaling properties. In this situation, the

commonly used EFIE method solved by the method of moments (MoM)

[4] with the Rao-Wilton-Glisson (RWG) basis function [73] suffers from a

“low frequency breakdown” problem, where vector potential gradually looses

its significance compared with the scalar potential part when the frequency

decreases, and the EFIE operator becomes singular [1]. Various approaches

have been proposed to overcome this problem in the last few years. One of the

most popular remedies is the loop-tree or loop-star decomposition [74], [75],

where the solenoidal and irrotational components of the unknown current canbe separated due to the quasi-Helmholtz decomposition (also known as Hodge

decomposition), to capture inductance physics and capacitance physics when

the frequency tends to zero. However, even after frequency normalization,

the matrix is still ill-conditioned. Preconditioning is necessary to improve

the convergence when iterative solvers are applied. Several effective precon-

ditioners have been proposed, either based on the basis-rearrangement, where

the favorable property of electrostatic problems is utilized [76], or based on

the near-field interactions, where the incomplete factorization with a heuris-

tic drop strategy is applied [77]. By using the Calderón identity and the dual

basis or Buffa-Christiansen basis function [78], [79], a more effective precon-

ditioner has been constructed [80]–[82]. The loop-tree or loop-star method

has also been implemented with the LMGF, which is more versatile in the

simulation of printed antenna and planarly integrated circuit [83]–[85].

However, one big issue associated with the loop-tree or loop-star method

is the loop

Date post:	06-Jul-2018
Category:	Documents
Upload:	kan
View:	216 times
Download:	0 times

FullText (40)

Documents