close

Вход

Забыли?

вход по аккаунту

?

Global modeling of microwave circuits

код для вставкиСкачать
INFORMATION TO USERS
This manuscript has been reproduced from the microfilm master. UMI films
the text directly from the original or copy submitted. Thus, some thesis and
dissertation copies are in typewriter face, while others may be from any type of
computer printer.
The quality of this reproduction is dependent upon the quality of the
copy submitted. Broken or indistinct print, colored or poor quality illustrations
and photographs, print bleedthrough, substandard margins, and improper
alignment can adversely affect reproduction.
In the unlikely event that the author did not send UMI a complete manuscript
and there are missing pages, these will be noted. Also, if unauthorized
copyright material had to be removed, a note will indicate the deletion.
Oversize materials (e.g., maps, drawings, charts) are reproduced by
sectioning the original, beginning at the upper left-hand comer and continuing
from left to right in equal sections with small overlaps.
Photographs included in the original manuscript have been reproduced
xerographicaliy in this copy.
Higher quality 6" x 9* black and white
photographic prints are available for any photographs or illustrations appearing
in this copy for an additional charge. Contact UMI directly to order.
ProQuest Information and Learning
300 North Zeeb Road. Ann Arbor. Ml 48106-1346 USA
800-521-0600
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
GLOBAL MODELING OF MICROWAVE CIRCUITS
by
Sebastien Goasguen
A Dissertation Presented in Partial Fulfillment
o f the Requirements for the Degree of
Doctor o f Philosophy
ARIZONA STATE UNIVERSITY
December 2001
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
UMI Number 3023891
UMI*
UMI Microform 3023891
Copyright 2001 by Bell & Howell Information and Learning Company.
All rights reserved. This microform edition is protected against
unauthorized copying under Title 17, United States Code.
Beil & Howell Information and Learning Company
300 North Zeeb Road
P.O. Box 1346
Ann Arbor, Ml 48106-1346
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
GLOBAL MODELING OF MICROWAVE CIRCUITS
by
Sebastien Goasguen
has been approved
August 2001
.
APPROVED:
'6 > .
.Chair
1/Ai jZJJ
wind/?- S
CP.
,
c W t l U I 'm
Supervisory Committee
ACCEPTED:
ipan inentjChair
Dean, Gradi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ABSTRACT
In this dissertation, the global modeling o f microwave transistors and circuits
is investigated. A good understanding o f the tools and methods necessary to perform such
simulations is attained. The main idea o f global modeling is the coupling o f
electromagnetic and semiconductor simulations. This dissertation presents the theory
required to perform such simulations. First, the numerical device simulation is studied
using a hydrodynamic model based on the moments o f Boltzmann’s transport equations.
Several approximations are made to simplify this model. The time and space derivatives
can be neglected, resulting in a simplified hydrodynamic model (SHM) and an energy
model (ENM). A drift diffusion model (DDM) is also implemented. Secondly, a hybrid
technique is developed using a standard finite difference time domain (FDTD) method
and a neural network model o f a field effect transistor (FET). This method achieved good
central processor unit (CPU) time requirements to allow for the small-signal and largesignal simulations o f microwave circuits. Different equivalent circuits can be
implemented inside a FDTD mesh to simulate various types o f transistors simulated using
various types of physics-based models.
Then, extensive work is done to investigate the potential o f using a wavelet
numerical method to solve the global modeling problem. A large bibliography of
wavelets research is presented, then a wavelet method is used to solve the DDM
characterizing a FET. This method is different than the standard wavelet-Galerkin
method and more suitable to deal with semiconductor equations. Good results are
obtained but this method generates a high computational overhead. This method is also
used to solve for a two-dimensional cavity enclosing a cylinder. Good compression
iii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
characteristics are obtained and a trade-off is found between compression and accuracy
the same type o f behavior can be achieved. Finally, fast wavelet transform (FWT) is used
to solve the Poisson’s equation. The work presented in this dissertation represents the
first step toward a numerical tool that can solve the global modeling problem efficiently.
iv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
I dedicate this dissertation to the four plus one people that are always in my heart. My
mother Marie-France, my father Jean-Etienne, my sisters Sabine and Fabienne, and
Caroline, my love, my friend, my partner.
v
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ACKNOWLEDGMENTS
I would like to express my gratitude to my advisor, Prof. Samir El-Ghazaly,
who believed in me and contributed his time and ideas towards the completion o f this
dissertation. The past three years have been filled with academic and personal
experiences that will be part o f me forever. I will always be thankful to Prof. El-Ghazaly
for being the one to make me come to Arizona. I would also like to thank my committee
members, Prof. Robert Grondin, Prof. Steven Goodnick, Prof. Bruno Welfert, and Prof.
Constantine Balanis for their interest and their time, devoted to the consideration o f this
work. All o f them taught me something and I will always be grateful.
Sincere appreciation and respect are due to my friends and colleagues o f the
Telecommunications Research Center for their help and support. In particular, I would
like to thank Kathy Muckenhim for dealing so efficiently with so many problems.
Salvatore Bellofiore will be forever my personal Windows/Linux/Unix FAQ solver,
many thanks goes to him. Yuri Tretyakov is a friend and one o f the most intelligent
people I have ever met, may he live a great life in the land o f the free. To Dr. Samir
Hammadi, Munes Tomeh, Aly Aly, Hazem Alassaly, Yasser Hussein, Mohammed
Mostageer, Kai Liu, Ekram Bhuiyan, Mohammed Waliullah, Marios, Marinos and
Stavros, I give many thanks for the discussions, the jokes and the arguments that took
place in the labs during the long hours o f work.
Finally, I will be in debt to all my friends who love me for who I am. May
they continue to accept my credit cards for the rest o f my life and may “The Great
Wallalah” keep an eye on us.
vi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
TABLE OF CONTENTS
Page
LIST OF TABLES................................................................................................................xi
LIST OF FIGURES.............................................................................................................xu
CHAPTER
1 INTRODUCTION.............................................................................................................. 1
1.1
Commercial CAD packages and their simulation techniques............................... 2
1.1.1
Linear frequency domain analysis................................................................ 5
1.1.2
Time-doraain simulation............................................................................... 6
1.1.3
Harmonic balance simulation........................................................................7
1.1.4
S-parameters from EM simulation................................................................ 7
1.1.5
Active Device Models....................................................................................8
1.2
Numerical techniques in global modeling............................................................. 9
1.2.1
Finite Difference Time Domain (FDTD)...................................................... 9
1.2.2
Finite Element Method (FEM)..................................................................... 10
1.2.3
Hybrid techniques........................................................................................ 10
1.3
New trends.............................................................................................................12
1.3.1
Neural networks modeling........................................................................... 13
1.3.2
Wavelets techniques..................................................................................... 14
1.4
Summary................................................................................................................ 15
2 DEVICE SIMULATION.................................................................................................17
2.1
Introduction............................................................................................................ 17
vii
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
CHAPTER
2.2
Page
The Boltzmann transport equation...................................................................... 18
2.2.1
The relaxation time approximation.............................................................20
2.2.2
Moments o f Boltzmann equation................................................................ 21
2.3
Presentation o f several physics-based model......................................................27
2.3.1
Hydrodynamic model...................................................................................28
2.3.2
Approximations o f the full hydrodynamic model....................................... 35
2.4
Equivalent circuit parameters extraction............................................................37
2.4.1
The small signal equivalent circuit............................................................. 37
2.4.2
Determination o f the lumped elements values............................................ 38
2.5
Hydrodynamic model simulations results........................................................... 40
2.5.1
Preliminary results....................................................................................... 41
2.5.2
Determination o f the lumped elements........................................................ 45
2.5.3
Comparison o f various hydrodynamic models........................................... 47
2.6
Summary.............................................................................................................. 49
3 THE NEURO FDTD METHOD..................................................................................... 50
3.1
Introduction...........................................................................................................50
3.2
The FDTD method............................................................................................... 51
3.2.1
FDTD update equations................................................................................53
3.2.2
Bench mark results........................................................................................ 56
3.2.3
Order o f accuracy.......................................................................................... 59
3.3
3.3.1
The extended FDTD method..............................................................................62
The resistor.................................................................................................... 63
viii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER
Page
3.3.2
Parallel resistor and capacitor..................................................................... 64
3.3.3
Parallel resistor and current source............................................................. 65
3.3.4
Multiple cell formulation.............................................................................66
3.3.5
Simulation results......................................................................................... 67
3.4
Principle o f neural networks................................................................................ 68
3.5
Hybrid ANN-FDTD method................................................................................ 73
3.6
Full-wave simulation results................................................................................ 76
3.7
Summary............................................................................................................... 84
4 ON THE USE OF WAVELETS FOR GLOBAL MODELING SIMULATIONS
4.1
85
Introduction........................................................................................................... 85
1} (9t) and a few definitions.......................................................................86
4.1.2
Multiresolution approximation....................................................................86
4 .1.3
Separable multiresolution............................................................................ 91
4.1.4
Two examples: the Haar and the Battle-Lemarie wavelets.......................94
4.2
A review o f wavelet research................................................................................98
4.2.1
The MRTD algorithm.................................................................................102
4.2.2
The wavelet collocation technique............................................................. 104
4.2.3
Mesh refinement and Daubechie wavelet expansion.................................105
4.2.4
Operators in wavelet bases and wavelet preconditioning......................... 106
4.3
Interpolating wavelets applied to semiconductor simulations........................... 107
4.3.1
Wavelet numerical techniques....................................................................107
ix
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER
Page
4.3.2
Interpolating wavelets................................................................................ 109
4.3.3
Drift diffusion example.............................................................................. 114
4.3.4
Results......................................................................................................... 116
4.4
Interpolating wavelets applied to a 2-D cavity................................................. 124
4.4.1
Sparse point representation.........................................................................125
4.4.2
Compression characteristics........................................................................128
4.5
Wavelet transform for efficient Poisson’s solver.............................................. 129
4.5.1
Poisson’s operator for FET structure......................................................... 129
4.5.2
Iterative solvers........................................................................................... 132
4.5.3
Wavelet transform and diagonal preconditioning...................................... 135
4.5.4
Results..........................................................................................................142
4.6
Summary............................................................................................................. 145
5 CONCLUSIONS............................................................................................................. 146
REFERENCES.................................................................................................................. 149
x
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
LIST OF TABLES
Table
Page
1.1
Yield o f wafer versus chip size...............................................................................2
1.2
Commercial CAD software..................................................................................... 3
2.1
MESFET simulation parameters........................................................................... 40
4.1
Condition number o f A versus wavelet type and resolution levels................... 142
xi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
LIST OF FIGURES
Figure
Page
1.1
Hybrid techniques for simulations o f transistors................................................. 12
2.1
Steady-state simulation flow-chart.......................................................................35
2.2
Small signal equivalent circuit o f a FET transistor.............................................. 37
2.3
I-V characteristics for the 0.3 pm gate-length MESFET..................................... 41
2.4
The potential distribution...................................................................................... 42
2.5
The electron density distribution...........................................................................42
2.6
The total current density J............................................ ........................................43
2.7
The energy distribution......................................................................................... 44
2.8
The temperature distribution................................................................................. 44
2.9
The drain to source capacitance............................................................................ 45
2.10
The gate to drain capacitance................................................................................ 46
2.11
The gate to source capacitance............................................................................. 46
2.12
Cross sections o f the carrier density, potential, andvelocity (x component).
First column is for CHM and second column is for DDM. First row is the
carrier density, second row is the potential and final rowis the velocity............48
2.13
I-V Characteristics of the MESFET. CHM in blue,SHM in black,ENM in
green and DDM in red. For three gate voltages: 0 V, -0.5 V and - I Volts......... 49
3.1
Time domain response.......................................................................................... 57
3.2
FDTD bench mark results......................................................................................58
3.3
Order o f accuracy o f finite different schemes....................................................... 62
3.4
Resistor modeling using Extended FDTD............................................................ 68
xii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure
Page
3.5
Three-layer artificial neural network....................................................................69
3.6
Error curves for training and testing data............................................................. 72
3.7
Basic amplifier with matching networks and intrinsic transistor model............. 73
3.8
Equivalent circuit parameters o f the transistor inserted in the FDTD grid......... 74
3.9
I-V characteristics o f the simulated MESFET, hydrodynamic model and
ANN...................................................................................................................... 77
3.10
Transconductance computed from the ANN results............................................ 78
3.11
Gate to Source charges, Qgs. Hydrodynamic model and ANN results............... 79
3.12
Gate to source capacitance,
3.13
Small-signal 5-parameters comparison between the FDTD method and HP-
C gs computed
from the ANN results...................... 79
Libra for the transistor. (Dashed lines: FDTD, Solid line: Libra)...................... 80
3.14
Small-signal 5-parameters comparison between the FDTD method and HPLibra for the amplifier. (Dashed lines: FDTD, Solid line: Libra)....................... 81
3.15
Time domain response o f the small-signal amplifier........................................... 81
3.16
Output voltage versus distance in the z-axis at three different times................... 82
3.17
Time domain response o f the large-signal simulation......................................... 83
3.18
Harmonics obtained from a large signal simulation using the ANN................... 84
4.1
Haar scaling function.............................................................................................94
4.2
Haar wavelet function............................................................................................95
4.3
Batlle-Lemarie scaling function............................................................................ 96
4.4
Battle-Lemarie wavelet function............................................................................97
4.5
A set o f dyadic grid.............................................................................................. 110
xiii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Page
Figure
4.6
Scaling function for (a) p=2 and (b) p = 4 ......................................................... 110
4.7
Left boundary scaling function forp= 4............................................................ 111
4.8
Wavelet subspaces representation..................................................................... 112
4.9
rV-Characteristics o f the simulated MESFET..................................................116
4.10
Contour plot o f the Electrostatic potential for a drain voltage o f 2 Volts and
an applied external gate voltage of - 1 Volt........................................................................117
4.11
Contour plot o f the carrier density (the labels are normalized by 1017) for 2
Volts drain voltage and - 1 Volt applied gate voltage.................................................... 118
4.12
Sparse Point representation o f the carrier density at 2 Volts drain voltage and
- I Volt applied gate voltage.............................................................................. 119
4.13
Time Domain drain, source and gate currents at 2 Volts drain voltage and - I
Volt applied gate voltage................................................................................... 120
4.14
Mesh Adaptability for three different values of wavelet threshold. ( Solid
lines p = 2 , dashed lines p = 4 )...........................................................................121
4.15
Relative Error on the drain current for three different values o f wavelet
threshold. ( Solid lines p = 2 , dashed lines p = 4 ) ............................................. 123
4.16
Electric field in the y direction, (a) t=520 ps, (b) t=840 ps.............................. 125
4.17
Sparse point representation o f the y component o f the electric field at (a)
t=520 ps and (b) t=840 p s .................................................................................127
4.18
Number o f unknowns remaining in the SPR for three different value o f
wavelet threshold...............................................................................................128
4.19
Conventional FET structure..............................................................................130
xiv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure
Page
4.20
Sparsity o f A using Haar wavelets......................................................................138
4.21
Sparsity o f A using db8 wavelets........................................................................139
4.22
Sparsity o f A using db20 wavelets......................................................................139
4.23
Sparsity o f A for one resolution level.................................................................140
4.24
Sparsity o f A for two resolution levels...............................................................141
4.25
Sparsity o f A for three resolution levels.............................................................141
4.26
Convergence behavior o f SSOR, CG and PCG on A and Am........................... 143
4.27
Convergence behavior of SSOR, CG and PCG on Amand Am......................... 144
xv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 1
INTRODUCTION
The standard modeling tools for computer aided design (CAD) are based on
circuit theory and equivalent circuits o f active and passive components. These equivalent
circuits are obtained after several run o f fabrication-measurements-optimization. This
process is costly and takes a considerable amount o f time and effort. In order to simulate
microwave circuits using new devices and new technology, there is a need for a
simulation tool that would give the performance trend o f a specific circuit before
fabrication. Such a tool is especially needed for high-frequency and high-speed circuits
that are very sensitive to electromagnetic coupling, process variation and propagation
effects. Moreover the density o f microwave circuits is constantly increasing due to new
technology. This reduces the cost of fabrication but increases the design complexity due
to the proximity o f the different components o f the circuit. Table l.l shows that as the
chip size increase you get less chips off the wafer and a lower yield.
Therefore, there has been an effort to analyze active devices based on
numerical techniques that incorporate the electromagnetic effects [1-5]. Such simulations
are called global modeling simulations, because they couple the electromagnetic fields to
the active phenomena.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Table 1.1 Yield of wafer versus chip size
Chip size (mm)
Yield (%)
Working circuits
per 3" wafer
Bare chip cost at
$4000 per wafer
Ix l
80
3600
l.l
2x2
70
800
5
5x5
45
80
50
7x7
30
25
160
10x10
20
9
440
In this chapter we review the standard CAD tools available and present the
simulation techniques that are used. Then, we will introduce the global modeling
techniques and present the previous research done on the subject in the last decade. We
will then review the new trends in CAD techniques and layout the main ideas developed
in this dissertation.
1.1
Commercial CAD packages and their simulation techniques
This section reviews the different CAD packages that are available for
microwave circuit design. At one end o f the scale, it is possible to design a simple MMIC
using a basic PC simulator with simple equivalent circuit models. At the other end o f the
scale, it is possible to use a fully integrated MMIC CAD workstation which has libraries
for layout automation and direct integration o f electromagnetic simulation for non­
standard structures. It is to be noted, however, that with the constant increase in clock
speed o f PC microprocessors, PC based CAD software can compete with Unix
workstations software and are only limited by the amount o f RAM necessary for the
simulations.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3
Table 1.2 Commercial CAD software
PRODUCT
COMPANY
HP-EEsof (HP range)
TYPE
MDS
Integrated package
MLS
Linear simulator
MNS
Harmonic balance
Impulse
Time-domain
Momentum
3D Planar eletromagnetic
HFSS
3D
MW Artwork Generator
electromagnetic
Arbitrary
Layout
HP-EEsof (EEsof range)
Touchstone
Linear
Linecalc
Physical parameters
Libra
Harmonic balance
MicroWave Spice (now in Time-domain
Libra)
Filter synthesis
E-Syn
System simulation
Communications
Design
Suite
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
4
Compact Software
Sonnet Software
EASi
Integrated suite
Serenade
Schematic capture
SuperCompact
Linear
Harmonica
Harmonic balance
SuperSpice
Time-domain
Microwave Explorer
3-D Planar electromagnetic
Em
3-D Planar electromagnetic
Xgeom
Layout entry
Emvu
Current display
The special requirements for microstrip circuit modeling were first met by
programs like Compact™, and later by Super-Compact™ and Touchstone™. The main
limitations of these types o f simulators are the limited accuracy, especially at high
microwave frequencies, and the inability to treat non-linear circuits. As a result, there has
been a rapid development o f other CAD techniques, especially for electromagnetic
simulation and for non-linear circuit modeling. Electromagnetic simulation is important
for microwave circuits because the approximate models used in conventional simulators
assume quasi-TEM behavior and neglect any interaction between individual
discontinuities and transmission lines.
Electromagnetic simulators can look at the entire picture, rather than breaking
the circuit down into building blocks, and can conduct a fill 1-wave analysis by solving
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
5
Maxwell's equations for two-dimensional and three-dimensional transmission-line
problems.
The key players in the major microwave CAD business are listed in Table 1.2.
HP/EEsof has now only one product line. The u n ifie d product is called Advanced Design
System (ADS) and has become the standard. There is a wide range o f software on offer,
and the function o f each o f the different types o f microwave CAD packages is described.
1.1.1
Linear frequency domain analysis
Microwave circuit analysis is largely performed in the frequency domain.
Working in the frequency-domain means that only the steady-state response of
distributed elements (transmission-lines) needs to be considered. This is very important,
because the multiple reflections encountered in a microwave circuit would take a lot of
time to solve in the time-domain. Furthermore, the frequency-dependence of
transmission-line parameters, resulting from dispersion, cannot readily be handled in the
time domain.
Linear frequency-domain analysis gives a time-independent solution for linear
networks with sinusoidal excitation at a range o f frequencies. Non-linear components are
treated with linearised small-signal models; the component’s complex behavior is
simplified by assuming that the ac signal is causing only small perturbations around the
dc operating point. This type of simulator generally uses Y-parameter matrix techniques
to solve for the steady-state frequency response.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.1.2
Time-domain simulation
The basic design o f an amplifier can be carried out entirely using linear
simulations based on the transistor’s measured S-parameters or a small-signal equivalent
circuit model. However, often the designer needs to be able to predict and optimise the
amplifier’s output power and intermodulation performance. Other circuits which require
large-signal simulations are mixers and oscillators. In these cases, non-linear simulation
is necessary. In addition, there are particular difficulties attached to using time-domain
simulation for microwave circuits: firstly, distributed circuits can take a long time to
reach a steady state, and, secondly, microstrip lines are dispersive, with frequencydependant behaviour. As a result, special microwave non-linear simulators have been
developed which use either a convolution technique or the harmonic-balance method.
In order to calculate the true response o f non-linear components such as
transistors it is necessary to work in the time-domain. This is because the equivalent
circuit elements vary with signal amplitude, and are thus functions o f time. Direct timedomain methods are generally based on Spice, with more convenient display options
added to create programmes such as MicroWave-Spice™ (now known as the Time
Domain Test Bench in Series IV). Spice models a circuit by representing the differential
equations whereby each component is modelled in terms o f finite difference equations
(such as i=Cdv/dt for a capacitor), and solving the set o f non-linear algebraic functions at
each time sample in an iterative manner.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
7
1.1.3
Harmonic balance simulation
For steady-state large-signal and repetitive (optimisation) analyses (such as
gain compression o f a power amplifier), the harmonic balance technique [6] is faster than
time-domain methods. This uses a combination o f time- and frequency-domain analysis.
The harmonic balance analysis can be performed in three ways. In the nodal form [7], the
circuit unknowns are the voltages at every node in the circuit. This method has proved to
be an efficient one when the number o f non-linear elements is large compared to the
number o f linear elements.
In the piecewise harmonic balance method, the circuit is segmented into a
linear part and a non-linear part connected to the linear part at a series o f ports. The linear
and non-linear interface defines corresponding nodes. The voltages and currents in the
non-linear components are calculated in the time domain. The linear elements are treated
in the frequency domain. Thus, every port is connected to a non-linear device. Fourier
transformation is then used to compare the currents and voltages o f the non-linear
components with the terminating impedances presented by the linear components at all
the harmonic frequencies used in the analysis. Any error is reduced by successive
iterations.
1.1.4
S-parameters from EM simulation
Electromagnetic simulation o f arbitrary planar structures is now a reality. This
is particularly useful for microstrip circuitry where standard CAD elements are invalid
due the close proximity o f bends, Tee-junctions, etc. So, after the bulk o f the design has
been carried out using CAD elements, the initial layout is carried out. Once the required
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
8
geometry o f the microstrip is known, the layout (o f a meandered line, for example) can
be imported directly into the em analysis package. The S-parameters are calculated and
the new circuit response checked. This process is now transparent with some o f the latest
CAD packages.
1.1.5
Active Device Models
For linear simulations these can be S-parameters, equivalent circuits for a
fixed device and bias point, or scaleable/bias-dependent equivalent circuits. A scaleable
model is one which can model different gate-widths. A bias-dependent model is one in
which the small-signal equivalent circuit elements vary with a bias parameter. For
example, the drain current o f a FET might be a parameter, this allows the designer to see
how the bias affects the gain o f an amplifier, but it is not a large-signal model.
Large-signal device models (transistors and diodes) are ones in which the DC
and RF parameters are all calculated dynamically for the instantaneous voltages applied
to the device. As a result, when using a large-signal model, the DC bias supplies must be
added to the simulation as voltage or current sources. Many different large-signal models
exist for FETs, such as the Curtice Cubic, Triquint's Own, HP Root, etc. Most foundries
provide at least one large-signal model type. The better foundries will offer a model
optimised for your particular simulator. The accuracy o f large-signal models still leaves a
lot to be desired, especially for highly non-linear circuits such as mixers. That is why
there is a need for a new era o f simulations that would enclose semiconductor and
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
electromagnetic models in one global model. One could even imagine adding mechanical
and thermal phenomena in the simulations.
1.2
Numerical techniques in global modeling
Global modeling could be the new simulation tool o f future CAD software.
Active devices would be characterized by solving physical models and passive devices
would be simulated using EM solver, the two characterizations being coupled for
consistency. However, in order to perform one global simulation o f a microwave circuit,
a considerable amount o f memory and CPU time is required. This extensive load on the
computational aspect o f the problem is a considerable drawback. A numerical technique
must be chosen such that the computational cost can be minimized, keeping the
simulation time within reasonable limits.
Nowadays, global simulations are unsuitable for design purposes. The time
domain simulations performed are too expensive to be part o f an iterative process, such
as an optimization algorithm. That is why, new numerical techniques must be invented to
make global modeling a reality in the design process. In the following section we will
briefly present the finite difference method, the finite element method and hybrid
techniques based on those methods.
1.2.1
Finite Difference Time Domain (FDTD)
The FDTD method was invented by K.S Yee in 1966 [8], when he developed
the so called Yee’s cell. This computational cell naturally interleaves electric and
magnetic grid points (grid points where the electric and magnetic fields are evaluated) so
that the rotational o f the electric or magnetic fields can be expressed in a straightforward
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
manner using finite difference approximations. This method remained in the shadow up
until the late 80’s when the simulation o f open structures was made possible using
absorbing boundary conditions (ABC). The Mur’s ABC [9] and the Perfectly Matched
Layer [10] where a big breakthrough that allowed the FDTD method to be used in a
variety o f problems. The basics of FDTD will be explained in chapter three but more
information can be found in the book by Alen Taflove [11] and in this paper referencing
most o f the publications on FDTD [12]. This numerical technique is known for its
versatility and simplicity of implementation. It has been the method o f choice to tackle
the global modeling challenge at an early stage [I].
1.2.2
Finite Element Method (FEM)
The FEM is radically different from the FDTD because it is a frequency
domain method. There is some research being done on finite element time domain FETD
techniques but this has not reached a wide use in the research community [13]. The FEM
is explained in [14]. It is the method used in HFSS. It has been widely used in areas other
than electromagnetics such as hydrodynamics and mechanics. In global modeling it is
used in conjunction with a circuit simulator to model multi-finger FET [IS], This shows
that hybrid techniques can be advantageous when dealing with a dual problem like global
modeling.
1.2.3
Hybrid techniques
Nowadays, commercial CAD software are based on circuit theory and
equivalent circuits. This was extensively discussed in the previous sections. Some o f the
current design software already contain an electromagnetic solver (Momentum, Sonnet)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
11
which generates a scattering parameters matrix that can be used in a circuit simulation. It
represents the first type o f hybrid simulation.
These techniques use the strength o f all the numerical tools used in order to
obtain a result faster and with a better accuracy, hi the context o f global modeling, two
different type o f structure need to be simulated, namely: the passive and the active part.
Therefore, one can imagine using the best possible technique for each part and find a way
to couple the different tool.
As explained before, FDTD is the numerical technique o f choice for
implementing a versatile EM solver. The problem resides in the difficulty to incorporate
active phenomena within the passive structure simulated. A solution is to use the
extended-FDTD method or the Lumped-FDTD. This method is presented and used in
chapter three. Lumped elements are added in the finite difference grid by modifying the
Maxwell’s equations to take into account the conduction current due to the lumped
elements. This technique was developed by Sui et al. [16], expanded by PiketMay [17]
and used to simulate microwave amplifier by Kuo et al. [18-19]. It was also used to
simulate an active antenna [20] and a diode mixer [21] and it was connected with SPICE
models [22] and transmission line matrix (TLM) method in [23],
Larique et al. [IS] used a finite element technique to import a scattering
parameters matrix within a circuit simulator. The active part is modelled by an equivalent
circuit that represents an elementary slice o f the transistor. The propagation along the
width o f the device is simulated by adding several slices. The equivalent circuit can be
based on measurements or on simulation using physical models (SILVACO).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Fig. 1.1 illustrates the different combination that can be used to analyze a
transistor. The more accurate method being the fully numerical one, where an EM solver
and a physical model are used. The method that uses a circuit simulator will be the more
efficient and is the one used in optimisation process. Unfortunately, this is also the one
the least accurate in terms of EM coupling.
Field
Effect
Transistor
Circuit
Simulator
Electromagnetic
Simulator
Physics-based
Simulator
Hybrid
Simulator
Fig. 1.1 Hybrid techniques for simulations of transistors
U
New trends
hi the past years there has been extensive efforts made to reduce the
computation time of global modeling simulations [24-27]. This can be done either by
using new hybrid techniques or by using new and more efficient numerical tools. New
ideas are very often found in a different engineering field than the one we are working in:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
electromagnetism. In signal processing and automatic, we can find very interesting tools
that have been used for a while but remain obscure to the micro-wave engineer. Neural
networks [28] and wavelets [29] represent such kind o f tools. These techniques have been
used for non-linear modeling and for data compression, they started to be used in the last
decade [30-31].
A general view o f the new trends in CAD was presented very nicely in a paper
by K.C Gupta [32]. It presents the advantages o f using neural networks to model
scattering matrices obtained from EM solver. Wavelets are much more mathematically
involved, in order to be implemented efficiently. The strength o f wavelet techniques is to
use multi-resolution in order to simulate details much more accurately and efficiently. It
was first introduced in the method o f moments [33].
1.3.1
Neural networks modeling
Artificial neural network (ANN) will be presented in details in chapter three.
We used ANN to model non-linear parameters obtained from hydrodynamic simulations.
The ANN is then used to update lumped elements values inside a FDTD grid. Large
signal simulations o f microwave circuits can be performed very efficiently [34-37].
ANN have been used in a wide range o f applications in electromagnetism. It
remains a universal non-linear approximator with multiple inputs and multiple outputs, hi
CAD, it was used mainly to create scaleable model that predicts the correct scattering
parameters for dimensions o f the structure that were not simulated [38]. It was used for
spiral inductors and interconnects [39] but also for optimization purposes [40-41].
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
14
Different types o f neural networks exist [42], the perceptron is the more common but the
radial basis function (RBF) is also widely used. Space mapping algorithms are also being
studied for knowledge-based design [43].
1.3.2
Wavelets techniques
Wavelets are very interesting because they implement a multi-resolution
analysis of a problem by filtering the details with a band pass filter [44]. It is then
possible to obtain a low resolution approximation o f an unknown and add the details that
have been filtered to reconstruct the unknown at a higher resolution.
They are widely used in signal processing for de-noising techniques and data
compression such as image compression [29]. Ten years ago, wavelets started to be used
in the method o f moments [33] to obtain sparse matrices. Sparse matrices contain few
non-zero elements, which speeds up iterative methods for linear system. Then, wavelets
started to be used in a more theoretical way [45]. It was shown that FDTD can be seen as
a wavelet based Galerkin technique using Haar wavelets [46].
That discovery launched
a
new
area
o f research
for numerical
electromagnetics: wavelets in EM. Numerous papers have been written [47-51]. We
personally used wavelets in semiconductor simulations to determine its potential in
global modeling techniques.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
1.4
Summary
Commercial CAD softwares are based on circuit theory and incorporate
transistor models obtained from measurements. They are limited in their ability to model
electromagnetic effect and high-frequency effect o f particle wave interaction. Despite the
incorporation o f built-in EM solver that generates scattering matrices for transmission
line structures, there is a need for a global modeling simulation.
In this dissertation we present a hybrid technique, using the extended FDTD
method that incorporates an equivalent circuit of a transistor inside an FDTD mesh.
Therefore the matching networks can be accurately simulated. The effects o f the active
device can be modeled up to a frequency where the numerical errors-due to the presence
o f lumped elements-are negligible. The lumped elements values are updated through a
neural network trained using numerical results from the hydrodynamic model simulation.
This allows us to get a very practical large-signal model and performed large-signal
simulation of the complete amplifier.
The main issue when dealing with a global modeling simulation is the
difference in scale when comparing the mesh used to discretize the passive parts and the
one used to discretize the active part. This problem can be solved by using time domain
diakoptics [S2-S7] that make use o f convolution to connect different part o f a circuit or
by using alternative techniques that speed up FDTD simulations like system
identification [58-61] and Prony’s method [62]. But there might be a way to develop a
numerical tool that deals with these different scales.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A new mathematical tool has started to being used in the last decade in
electromagnetic: a wavelet approach. Ideally such a technique can implement
multiresolution. Therefore in this research we focused on investigating the potential o f
wavelets numerical techniques to deal with the global modeling approach. We present a
wavelet like method that helps us to generate a non-uniform mesh adaptive in time.
The second chapter presents the fundamentals o f device simulations. The
hydrodynamic model and its approximations are presented. Simulations results are also
shown such as two-dimensional distribution o f the potential, the carrier density and the
energy. Lumped elements are extracted and plotted. Then the third chapter presents the
NEURO-FDTD technique. Principles o f neural networks are given and the extended
FDTD method is explained. Full-wave simulation results are also shown. The fourth
chapter presents
an extended bibliographical research on wavelets
used in
electromagnetics and in the problem o f solving PDE’s. Then the interpolating wavelet
technique is introduced and results are shown for device simulation and a twodimensional cavity enclosing a cylinder. The wavelet preconditioning technique is
explained in detail to solve the electrostatic potential o f a FET structure.
Finally, conclusions are drawn in chapter five and future recommendations in
state o f the art global modeling are stated.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 2
DEVICE SIMULATION
2.1
Introduction
Physical models provide a link between physical parameters (doping profile,
gate length, recessed gate depth, etc...) and electrical performance parameters (e.g. DC
characteristics, RF transconductance, junction capacitances, etc...). For these reasons,
physical device models should be essential tools not only in the development of new
devices but also in performance optimization and yield optimization o f microwave
circuits. The major limitation o f such models, however, is that they are computationally
intensive, preventing their direct use in circuit design. In most cases, it is inevitable to
introduce simplifying assumptions in order to make the model computationally tractable.
Most commonly used physics-based models can be classified into two
categories: particle based models and fluid-based or hydrodynamic models. The first
category is represented by the Monte Carlo technique [63-64]. The second is based on a
set o f conservation or continuity equations that can be derived from the Boltzmann
Transport Equation (BTE) [65]. These models usually involve several approximations
and they comprise from the simplest to the more complex: Analytical models, DriftDiffusion models, Energy models and Full-hydrodynamic models.
The semiconductor equations consist o f a set o f partial differential equations,
which must be solved subject to a pre-defined set o f boundary conditions over a specified
domain. Although generalized solutions are not available for all devices there are many
closed-form analytical solutions available for a wide variety o f devices. These closed-
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
18
form solutions are more suited to lower frequency devices, such as vertical diodes and
long gate length FET’s, with predominantly one-dimensional field and carrier profiles.
They are severely limited in their range o f application and accuracy because o f the multi­
dimensional non-linear nature o f most modem devices, such as planar transistors.
A more generalized method o f solution frequently applied to the
semiconductor equations is to solve them using numerical techniques [66-70]. This latter
approach requires considerably more computer time than the closed-form methods, but
usually produces more accurate results and provides greater flexibility.
In this chapter we show how to perform such numerical simulations. First, the
balance equations that characterize the device are derived. These balance equations
represent the moments o f Boltzmann transport equation. Then we show how to obtain
different models using basic approximations based on inertia effects. Finally, numerical
simulations are performed and it is shown how to extract an equivalent circuit.
2.2
The Boltzmann transport equation
The distribution function gives the probability o f finding a carrier at time t
located at position r with momentum p. The Boltzmann transport equation (BTE) relates
any net change of the distribution function/fop,f) to the different physical mechanisms
occurring inside a device. The BTE is derived according to the theory explained in [65].
The differential o f the distribution function can be written as follows:
(2l)
dt
dr d t
dp dt
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
19
dr
do
We set v = — and F = — . The net increase in carrier density due to generation
dt
dt
recombination processes is defined as v
dt
then, the net change o f carrier density in
G -R
time can then be expressed as:
(2.2)
dt
dt
G -R
Which leads to the following expression using equation 2.1
“ +vVr/ + F V / = —
dt
r
p
dt
(2.3)
G -R
The generation-recombination processes can be divided in two principal
mechanisms. First, the actual generation-recombination such as photogeneration. Such a
process is characterized by the function s ( r ,p ,t) . Then, collisions between carriers that
send one carrier from one momentum state to another. Separating these two mechanisms,
we write the continuity equation for the distribution function. This equation is the
Boltzmann transport equation.
(2.4)
$ . + r t rf + FV ' f m s( r , , , l y £ .
colt
When studying unipolar devices generation-recombination can be neglected. In the rest o f
this section s(r,p ,f)w ill be neglected. If one decides to investigate bipolar devices such
as BJT or HBT, this assumption will not be possible.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
20
The probability that carrier will scatter from momentum state p to momentum state p’ is
defined by the transition rate S ( p ,p ) . The probability that a carrier is at momentum
state p is / ( p ) , and the probability that the momentum state p’ is empty is [ l - / ( p ) ] .
Therefore the net rate o f increase in the distribution function due to scattering processes
can be written as:
&
= Z /( p ') [ l - /G O ]s (p '>p )~ %p' /O O P ~ / 0 7')] s fa p ')
dt coll P'
2.2.1
(2.5)
The relaxation time approximation
Equation 2.5 and 2.4 form a complicated differential equation that requires
simplifying assumptions in order to be solved easily. The generation recombination terms
has already been neglected because in this study we are mainly concerned with unipolar
devices. Another widely used approximation is the relaxation time approximation (RTA).
This assumption simplifies dramatically the collision term.
The distribution function can be decomposed in a symmetric and an anti-symmetric part.
f ( r , p ,t ) = f s (r, p ,t ) + f A(r, p ,t)
(2.6)
The first term is the symmetric component and is assumed to be large whereas the second
term is anti-symmetric and is assumed to be small. At equilibrium, f s - f 0 and f A = 0 .
The collision term o f equation 2.4 can be written as:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Under low-field transport, the symmetric component o f the collision term vanishes and
can be set to zero. But under high-field transport this is not the case. In the remainder of
this study we will consider the RTA that is valid for low-field transport and for scattering
mechanism that are elastic, isotropic or both. The RTA is shown in equation 2.8.
(2.8)
Where t is a characteristic time, which describes how the distribution function relaxes
back to its equilibrium value once the driving forces are removed.
2.2.2
Moments of Boltzmann equation
Boltzmann transport equation is extremely difficult to solve directly in the
form presented in equation 2.4 even though assumptions have already been made. A
common approach to simplify the derivation o f the solution o f BTE is to use alternative
equations known as balance equations. They represent the moments o f the BTE. An
infinite number of balance equations can be derived but fortunately only a small number
are required to characterize entirely a device. In this section we will show how to derive
the current continuity equation, the energy conservation equation and the momentum
conservation equation. Those three equations are moments o f the BTE and are much
easier to solve than the standard BTE presented in equation 2.4. Let’s now derive a
generic balance equation o f the BTE.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
22
Given a function o f momentum /> -> $ (/> ), the total value o f the quantity
associated with O (p ) is the weighted sum:
(2.9)
To find the balance equation for n0 , we multiply the BTE by <D(p)/Q . The first three
moments o f the BTE are obtained by choosing d>(p) such that n# represents the carrier
density, the carrier momentum or the carrier energy. Respectively, this means
Q (p ) = l,p ,E ( p ) . Where E represents the energy. Using equation 2.4 and neglecting
the generation recombination term, we get:
Where E is the electric field. This is not to be misunderstood with the scalar quantity E,
the energy. The first term in that equation can be modified because <D(p) is independent
o f time. It reduces to:
(2. 11)
The second term in equation 2.10 can be modified in the same way because 0 ( p ) is
independent o f position. The spatial gradient can then be moved outside o f the
summation.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
i£ « > ( P> .v ,/= w ,
(212)
f . ( w ) - 5 Z * 0 >K
(2 1 2 )
Where
Represents the flux associated with
.
The third term in equation 2.10 is re-written in the following way:
( - , ) £ . j ; < |) ( p ^ , / = ( - « ) £ . 5 ; v , ( 4 ) / ) + « £ . £ / V , 4 .
P
P
(2.13)
P
The first term on the right hand-side o f equation 2.13 is zero because the distribution
function rapidly reaches zero for large p. We can define a generation rate:
( - « ) E .£ < t> ( p ) V = - G .
<2' ,4 >
P
Where
O ,( r ,0 = H ) f . ^ Z / V , ® ( p )
(2 1 5 )
p
We can see that the BTE starts to look like a balance equation where the net change o f a
quantity in a given volume is given by the gradient o f the flux associated with that
quantity plus a generation rate. A recombination rate is missing. Now we will show that
the collision term in the BTE can be expressed as a recombination rate.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
24
(2.16)
Using the RTA, equation 2.16 reduces to:
(2.17)
Finally, these mathematical manipulations give us the final form o f the balance equation
for the moments o f the BTE
To analyze most devices, only the first three moments o f the BTE are necessary. This
means that we have to derive three balance equations based on equation 2.18 using
® { p ) = I P ,E { p ) '
2.2.2.1
The carrier density balance equation
hi this case we take d>(/?) = 1. This means that n<p(r,f) = n (r,r),
Vp<D(p)=0and that — = 0 . Therefore equation 2.18 reduces to:
(2.19)
This equation is most often known as the current continuity equation.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
25
2.2.2.2
The momentum balance equation
To derive the momentum balance equation we let 0 ( p ) = p <?with q=x,y,z.
Three balance equations can be derived, one for each vector component. From equation
2.9, we get:
/ \ I v
<•/
\
•
" 4 r > 0 = o L W / ,-, ) =nm
(2.20)
p
Which reduces equation 2.11 to:
Where
d n ^ t)
d (n m \)
dt
dt
(2*2 l>
is the q component o f the average carrier velocity. The flux associated with
appears as a flux o f momentum and is derived from equation 2.12.
„
I v
,
(2.22)
The carrier velocity can be decomposed in two components, one that represents the
velocity due to the applied field and another one due to the collision process
representative o f the thermal energy. Therefore 2.22 can be written as:
= n m w q + nkaTc
(2.23)
The generation term G0 by noticing that:
VcD(p )|
=j?
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2-24)
26
Where uqis the unit vector in the q direction. This leads to:
= (-q )n E q
(2.25)
Where Eq is the component of the electric field in the q direction. The recombination term
defined in equation 2.17 can be re-expressed as:
R* = — (nni (v^ -
))
(226)
This finally allows us to write the final form o f the momentum balance equation:
™
2.2.2.3
= ~V- (nm\ ) - v r(nkBTc)+ (-q )n E q + j - ( n m ( v llq - v ^ ) )
The energy balance equation
This represents the third moment o f the BTE. Here, we let <t>(p) = E ( p ) .
Where E (p ) represents the kinetic energy o f a carrier at momentum state p . Therefore,
equation 2.11 reduces to:
d("»M ) | „ £ . , 8/_ a (n £ )
a*
dt
r»
Cl*?
y r 'fl#
/ dt
(2.28)
A*
dt
Where E is the average kinetic energy. The associated energy flux expressed by equation
2.12 can be written as:
„
m v
2 r- nm —
• = 2Q ^ - V
_ 2~Vv
ttt
(2.29)
V
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
27
By using the same argument as in the previous section relating the average drift velocity
and the thermal velocity, we get
!?=?(F+frflrc)
(2-3°)
The generation of energy comes from the electric field that accelerates the carriers. Using
equation 2.20, equation 2.15 becomes:
p
£ r [j M - 9 >>£.p
12m
(2.31)
Finally, the recombination term is expressed as:
R ^ = -L--/f
n [ E - E^° )
*E
E
(2.32)
Now, the complete energy balance equation can be written as follows:
(2.33)
= -V r (nEv
2.3
(-q)nE .v +— n{^E - E° ^
Presentation of several physics-based model
Up until the 1970’s, semiconductor devices were well modeled with a drift-
diffusion transport approach. The drifi-diffusion approach includes a drift velocity
controlled by the electric field and a diffusion term in opposite direction from the carrier
density gradients. A typical drift-diffusion model is shown on the following equations
valid for a unipolar device:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
28
(2.34)
Jn =QnMnE+ qD nVn
Where Jn is the electron current density,
(2.35)
the electron mobility and D„ the carrier
diffusion coefficient. For a bipolar device, one would have to consider the same
equations for the holes and add a generation-recombination term.
It assumes that the microscopic distribution o f the momentum and energy over
the charge carriers at any location and time inside the device is equal to the one we would
find in a large sample with a DC field equal to the local instantaneous field.
However, the assumptions in the Drift-Diffusion model break down for sub­
micron devices, where carrier transport is predominantly non-stationary. For small size
devices with gate-length less than 0.5 microns, non-stationary transport effects, such as
hot electron effects and velocity over-shoot, become very important and must be
accounted for in the device model [66,70].
The models that can be used are called hydrodynamic model and are based on
the balance equations derived in the previous section. We will know show how to use the
balance equations to perform a real device simulation.
2.3.1
Hydrodynamic model
The model used during this study is a hydrodynamic model with no
simplifying assumptions on the energy or momentum equations [66]. Assumptions will
be made later on [71].
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
29
It is based on the semi-classical hydrodynamic model that states conservation
o f particles, momentum and energy. It is referred to as semi-classical because the
descriptions used for the band structure and scattering processes are generated by
quantum mechanical calculations.
The first three moments o f the BTE may be summarized for electrons
(neglecting generation and recombination) as:
(2.36)
— + V (nv)= 0
dt
v ’
(2.37)
dt
r.
Where n is the electron density, t is time, v is the average electron velocity, e is the
average electron energy, q is the electronic charge, \lq is the Boltzmann’s constant, T is
the electronic temperature, So is the equilibrium thermal energy, px is the x-component o f
the electron momentum and t m and t e are the momentum and energy relaxation times
respectively. These equations have to be linked with equations 2.19,2.33 and 2.27.
The electric field is obtained from quasi-static simulations using Poisson’s equation:
(2.39)
Similar equations are written for the momentum in the other directions.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
30
The electronic current density distribution J inside the active device at any time is given
by:
J = -q n v
(2-4°)
These equations coupled with Poisson’s equation form a complete system that
can be solved for the carrier density, momentum and energy. They give the DC
characteristics o f the device.
2.3.1.1
Single Gas approximation
The conservation equations are obtained by proper integration o f the
Boltzmann’s transport equation over the momentum space and averaging over the
multivalley conduction band [72]. Indeed, in semiconductors with multi-valley energy
structures, such as GaAs, there exist one electron gas with different distribution function
for each o f its conduction band valleys. Therefore, the BTE and the hydrodynamic
equations are separately valid for each o f the conduction band valleys. To reduce the
number o f partial differential equations that need to be solved, the single electron gas
approximation is almost always used. This consists o f using average quantities over all
main valleys for the carrier density momentum and energy [73]. For example, the new
quantities used in the single gas approximation are defined as follows:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Where the summations are taking over all conduction band valleys.
2.3.1.2
Transport parameters
The solution o f the balance equations requires prior knowledge o f different
transport parameters such as the ensemble energy and momentum relaxation times, the
electron effective mass, mobility and temperature. Two main approaches are commonly
used to evaluate these parameters. One approach is to guess the form o f the distribution
and then use the balance equations to solve for the parameters in this function. The most
commonly used guess is the displaced Maxwellian. Another approach to evaluate these
transport parameters is based on Ensemble Monte Carlo simulations. The parameters are
evaluated under steady-state conditions in bulk GaAs. Several simulations are performed
using different electric field and impurity concentration profiles. The resulting statistical
information can then be empirically fitted in analytical expressions suitable for numerical
simulations.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
32
So, the electron energy can be expressed in terms o f the steady-state applied field E<
0.316
£■= 0.316—
(2.42)
Ess
1+lTj
for Ess ^ 12.5 kV/cm and:
( E s s - 12.5)
£ = 0.308+3.353--------5— 1
103
(2.43)
for Es. >12.5 kV/cm.
The electron saturation velocity is given by:
(2.44)
1.53+0.9
U 04J
vf l = ! 07
104
and the electron drift velocity is given by:
(2.45)
f-M
1,4500J
1+ f - ^ - l
UsooJ
Where p0 is the low-field mobility and is given by [74].
(2.46)
8000
M0 =
1+ yjNd
X
lO’17
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
33
The electron temperature, Te, is expressed as:
(2.47)
*T
■o
1
=— 0.085 I —
1+
'f + 0 0 8 ] 13
^ 0.115 ,
For e < 0.3 eV and as in equation 2.16 for e > 0.3 eV.
s - 0 .3
.1.6
(2.48)
T = — 0.05+0.03 1+
. 0.075
e ka
The electron effective mass is approximated by:
0.0623mo
(2.49)
[0.0623 + 0.00526(£„ -1.944)] m0
m* = • [0.074+0.0216(£w-4 .1 6 7 )] m0
[0.1172 + 0.01053(£„-6 .1 9 4 )] m0
[0.1397+0.005(£Jt -8 .3 3 )] m0
Finally, the electron mobility and the energy and momentum relaxation time are
respectively evaluated from:
Vj
H \ £) = ~z~
(2.50)
And
(2.51)
Ts ( e ) ~ -------
v.E.
d^u
The final solution is obtained in a self-consistent evaluation o f the three
conservation equations in conjunction with Poisson’s equation.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
34
It was found that the order in which these equations are solved is critical to the
stability o f the solution. This is due to the different stiffness o f the momentum and the
energy conservation equations.
As shown in the following flowchart, Poisson’s equation is first used to
update the electric field. Because the momentum relaxation time is about one order o f
magnitude shorter than dielectric and energy relaxation times, the momentum equation is
solved first using the new values o f the electric field. Then, the carrier energy is updated
and used to compute the new transport parameters. Finally, the continuity equation is
solved for a new carrier density distribution. This process is repeated until a stable selfconsistent solution is obtained. This is shown in Fig. 2.2.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
35
Start
J
Initialization
Evaluate Fields
Solve Momentum Eq.
Solve Energy Eq.
Update Transport
Parameters
No
time=tmax
I Yes
Stop
Fig. 2.1 Steady-state simulation flow-chart.
2.3.2
Approximations of the full hydrodynamic model
As explained in the previous sections, the semiconductor models used are
based on the moments o f Boltzmann’s transport equations obtained by integration over
the momentum space. Three equations need to be solved together with Poisson’s
equations in order to get the quasi-static characteristics o f the transistor.
The solution o f this system o f partial differential equations represents the
complete hydrodynamic model (CHM). Simplified models are obtained by neglecting
time and/or space derivatives in the momentum equation [71]. The simplified
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
36
hydrodynamic model (SHM) is given by (2.52), which is the momentum conservation
equation where the space derivatives have been neglected. The energy model (ENM) is
given by (2.53), which is equation (2.27) where space and time derivatives have been
neglected. Finally, a drift diffusion model (DDM) is derived by assuming that the
mobility is a function o f the electric field (2.54).
*1"
8t
T\=-0L P-.-.PA
dx
(2.52)
zm
(2.53)
(2.54)
Vjflf „
These different models were used for quasi-static simulations. It was shown that the
CHM model predicts a higher current compared to DDM due to overshoot phenomena
taking place in the device. Deep level traps, surface charge and magnetic effects could be
included [75-79] as well as thermal effects [80-81].
It is to be noted that these models have different computational cost. The
CHM model is the most costly and the DDM model is the less costly. The CPU time
spent to obtain the quasi-static characteristics o f a two-dimensional FET can vary by 30%
if one chooses to use the DDM or the CHM.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
37
2.4
Equivalent circuit parameters extraction
Although the potential of physical models for nonlinear CAD has been
demonstrated, equivalent circuit models remain the most commonly used approach in
circuit design because of its computational efficiency. There is an increasing trend,
however, in generating physics-based equivalent circuit models. These are circuit models,
whose parameters are extracted from physical model simulations, rather than from
experimental measurements.
2.4.1
The small signal equivalent circuit
Fig. 2.2 shows a typical equivalent circuit model for a FET with its intrinsic
and extrinsic parameters [82].
[ntrinik part
Lg
Rg
Cgd
Rd
Ld
'VW
Drain
Source
Fig. 2.2 Small signal equivalent circuit of a FET transistor.
In our simulations, only the intrinsic elements can be extracted. The extrinsic
elements are found by cold measurements and package analysis.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
38
2.4.2
Determination of the lumped elements values
The values of these elements can be computed from the variations in voltages,
currents and charges due to small changes in the DC bias voltages and/or currents. For
instance, the small-signal gate-to-source capacitance can be computed as:
A Q ,,
(2.54)
*'** " AV lv"=0
p
Where Qg is the charge on the gate electrode while
and Vgd are respectively the gate-
to-source and gate-to-drain voltages. The charges Qg is evaluated from the integral form
of Gauss’s law given as:
V i?
V
rfv=J/0
dv
(2.55)
V
Where theintegral is taken over a volume containing thegate electrode. Assuming that
the longitudinal fields are zero and that the transversal Heldsare uniform along the width
of the device (2D case), the volume integral can be reduced to a contour integral of the
form:
,
G,
[e V £ dl = —
J
W
V
"
(2.56)
Where W is the width of the device and the integral is over a contour enclosing the gate
electrode. Using the divergence integral theorem we obtain the following Hnal expression
for Q, which can be easily evaluated numerically:
Q%= W $ e E n dl
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2*57)
39
Where n is the unit vector normal to the contour path. The equations for computing the
remaining intrinsic small-signal parameter are given below:
r
AC, |
(2.58)
= AVrrfJlv*’=0
r
_Aai
* ~ AV* lv,- =0
. A/„ ,
gm ~ A V
(2-59)
(2.60)
*»
A/„ |
(2.61)
f- " A V > _AVpf
(2-62)
^ " A / ** > ^ °
Gate to drain capacitance (2.58), drain to source capacitance (2.59), transconductance
(2.60), output conductance (2.61), input resistance (2.62) and finally, cut-off frequency
by (2.63).
,
(2.63)
2*C.,
An important figure of merit that can be evaluated from these parameters is the device
unit gain cut-off fiequency, ft, given by (2.63).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40
2.5
Hydrodynamic model simulations results
The metal-semiconductor field-effect transistor (MESFET) is a unipolar
transistor where conduction involves predominantly one type of carriers.
Table 2.1 MESFET simulation parameters.
Drain and source contacts
0.5 pm
Gate-source separation
0.5 pm
Gate-drain separation
2.0 pm
Device thickness
0.4 pm
Active layer thickness
0.1 pm
Active layer doping
2 x 10l / cm'J
Substrate doping
1014 cm'3
Schottky barrier height
0.8 V
The geometry of the MESFET device used in this simulation consists of a 0.1
pm-thick active layer with doping concentration Nd = 2 x 1017 cm'3 on a 0.4 pm-thick
buffer layer. The doping of the buffer layer is taken to be 10u cm'3 and the doping profile
between the active and the buffer layers is assumed abrupt. A gate length of 0.3 pm was
used, with 0.5 pm gate-to-source separation and 1 pm gate-to-drain separation. The builtin potential at the Schottky contact under the gate is set to 0.8V. The source and drain
electrodes are assumed to form perfect ohmic contacts with the semiconductor region.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
41
The lattice temperature is taken as 300K. Table 2.1 summarizes all the important physical
parameters used in the simulation.
2.5.1
Preliminary results
First, the drain-current versus drain-voltage characteristics of the MESFET are
presented Fig. 2.3.
300
ld(mAAnm)
280
200
180
100
-80
Veto
Fig. 2 3 I-V characteristics for the 0 3 |im gate-length MESFET.
Fig. 2.4 to Fig. 2.6 show the electrostatic potential, electron density and current density
distributions for the hydrodynamic transport model at the bias point Vgs=-lV and
Vds=3V. hi Fig. 2.4, we can see the equipotential lines. The most of the voltage drop
takes place in the region between the gate and the drain.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Fig. 2.4 The potential distribution.
Dmity(fcm9)
Fig. 2.5 The electron density distribution.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
43
In Fig. 2.5, the shape of depletion region is presented. The current path for the
hydrodynamic transport model is underlined in the vector current density plot (Fig. 2.6).
It is clear that high values of the current density penetrate deeply into the substrate
region, not only under the gate, but also in the quasi-neutral regions on either side of the
gate depletion region.
r
■ i
1
Ornate d r c o w rt
i
——t
t
i -----
.
A.
At.
A.
~vli
M
. .\i
,.ii
10
10
30
40
50
80
Fig. 2.6 The total current density J.
Because the region between the gate and the drain is a high-field region, the
carriers are accelerated and injected into the buffer layer. As a result, carriers with high
energies are produced. These hot electrons gain energies in excess, especially at the gate
edge on the drain side, as seen with the energy distribution o f the hydrodynamic transport
model on Fig. 2.7.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
44
Fig. 2.7 The energy distribution.
These results are confirmed Fig. 2.8, which shows the electron temperature.
T n v n l« (IQ
10
»
30
40
00
N
70
Fig. 2.8 The temperature distribution.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
45
2.5.2
Determination of the lumped elements.
By using equations 2.59-2.63 one can obtain the following lumped element
values.
First the drain-to-source capacitance (Fig. 2.9), the gate-to-drain capacitance
(Fig. 2.10) and the gate-to-source capacitance (Fig. 2.11) are proposed for two different
bias points.
Cds(pFfrnm)
IS
1.4
1.2
as
0.2
1.8
2.8
3.8
4.5
Vd
Fig. 2.9 The drain to source capacitance
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
46
Cgd(pfAnm)
O lO S
as
Vd
49
Fig. 2.10The gate to drain capacitance
Cga(pFftnm)
ar
a»
at
ass
-is
as
VQI
Fig. 2.11The gate to source capacitance
From these curves, we can develop an intrinsic small-signal model, which
represents the behavior of a transistor at several bias points.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
47
2.5.3
Comparison of various hydrodynamic models
The quasi-static characteristics are obtained by solving the hydrodynamic
equations together with Poisson’s equation on a two-dimensional domain. The FET under
analysis is a conventional MESFET of gate length 0.3 pm with a 0.1 pm active channel
doped at 2.1&7 A/cm3 and a buffer layer of 0.4 pm doped at 2.1014 A/cm3. The source to
gate separation is 0.5 pm and the gate to drain separation is 2.0 pm. The transistor is
biased at a gate voltage
V s •0.5
V and at a drain voltage of V*= 3 V with a 0.8 V
diffusion voltage.
Fig. 2.12. Shows the carrier density, potential and velocity (x-component)
obtained using the CHM and DDM models. It demonstrates that the hydrodynamic model
predicts that more electrons move into the buffer layer and that overshoot phenomena
takes place with the velocity exceeding the saturation velocity.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
48
30
30
20
20
10
10
20
40
60
30
30
20
20
10
10
20
40
60
------------------ :
20
40
60
20
40
60
20
40
60
30
20
10
Fig. 2.12Cross sections of the carrier density, potential, and velocity (x component).
First column is for CHM and second column is for DDM. First row is the carrier
density, second row is the potential and final row is the velocity.
Fig. 2.13 shows the IV characteristics obtained by steady state simulations. As
predicted, the energy models give a much higher current than the DDM. The difference
between SHM and ENM is indistinguishable on the graph. Simulations of different gate
length transistor will be performed to determine when the SHM and ENM differ.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
49
300
250
200
E
E
150
100
0.5
2.5
1.5
3.5
Vda
Fig. 2.13I-V Characteristics of the MESFET. CHM in blue, SHM in black, ENM in
green and DDM in red. For three gate voltages: 0 V, -0.5 V and -1 Volts.
2.6
Summary
We showed how to derive the balance equations of the BTE that form the
hydrodynamic model. Various hydrodynamic models were obtained by neglecting inertia
effects. These models have been used to characterize a standard FET using finite
difference techniques. Finally it was shown how to extract equivalent circuit parameters
from those numerical simulations.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 3
THE NEURO FDTD METHOD
3.1
Introduction
Artificial Neural Network (ANN) is a new trend in mm-wave CAD that
dramatically reduces the computation time o f electromagnetic (EM) models, replacing
the expensive EM model with ANN trained by EM simulation results [32]. This new
approach has been successfully used [39-40] to model passive devices, such as CPW
structures, spiral inductors, and interconnects. In [28] and [41], neural networks are used
to design a microstrip corporate feed, and also to optimize interconnects in high-speed
VLSI circuits. ANNs have also been used to model the high non-linearities o f transistors
obtained by measurements or empirical models [30].
Commercial design softwares are based on circuit approach. In order to
perform a more accurate and efficient simulation o f MMICs, the extended FDTD method
employs equivalent current/voltage sources to substitute the device inside the FDTD grid
[19]. Lumped elements such as resistor, capacitor and inductor can also be directly
included in the FDTD marching time algorithm [16-17]. This technique has been used to
simulate active antenna and amplifier [18,20] and was connected to SPICE simulator in
order to analyze more complicated transistor models or loads [22].
The extended FDTD approach can include non-linear devices in the FDTD
grid as in [21,34]. This method leads to an equivalent circuit that characterizes the wavedevice interactions, which provides a good approximation o f a global modeling
simulation providing that the transistor model comes from a semi-conductor simulation.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
51
In this chapter, we propose to model the steady state solution o f the hydrodynamic
equations by an artificial neural network that accurately, and efficiently predicts the
behavior o f the simulated MESFET. The transistor is implemented in an extended FDTD
code to simulate an amplifier. The non-linearities o f the MESFET are described by the
ANN, which updates the circuit parameters values inside the FDTD mesh according to
the electromagnetic field computed. This new approach provides a very efficient, and
practical first order global modeling o f a MMIC.
The chapter is organized as follows: Section 3.2 gives an introduction about
the FDTD method. In section 3.3 the extended FDTD approach is reviewed. The
following section describes artificial neural network, and gives a very comprehensive
matrix formulation. The interaction between FDTD and the ANN is explained in section
3.5. Finally, section 3.6 presents the results o f the ANN model compared with the
hydrodynamic model. To validate the method, the frequency response o f the small signal
simulation is compared with HP-Libra. The time domain responses o f the small signal,
and large signal simulations are also shown.
3.2
The FDTD method
The finite difference time domain (FDTD) method is a numerical technique
used to solve Maxwell’s equations in the time domain. It approximates the spatial
derivatives using central differencing scheme and a standard backward Euler scheme to
solve the resulting ordinary differential equations (ODE). The Yee’s algorithm [8] based
on a leapfrog scheme allows to take care o f the coupling between the equations. The
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
52
equations to be solved are as presented below. They are known as Faraday’s and
Ampere’s equations.
(3.1)
Vx£ =
dt
m
with
\Jn =P'H
\ j e = oE
Which leads to six scalar PDEs. In the case o f a homogeneous medium with no lumped
elements:
(3.2)
BEt ^ 1
dt
e
dy
dz
^ L =I tJ L J J L -r t,
dx
dt
e dz
dE. ^ 1 ( dHy
dt
dHx
dx dy
(3.4)
-o E .
(3.5)
S ^ = J_
dt
(3-3)
u Kdz
dy
(3.6)
dt
M \d x
dz
(3.7)
dH. ^ 1
dt
e
dy
dx
For a lossless medium, they get simplified to become:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.2.1
FDTD update equations
In order to derive the update equations for the electric and the magnetic field,
we need to approximate the derivatives with some type o f numerical scheme. The FDTD
as said earlier uses a central difference scheme. Section 3.2.3 will present in more details
the finite difference (FD) technique used to approximate the derivatives. As a small
introduction to FD, let’s F represent the electric o f the magnetic field. Then finite
difference can be used to approximate the derivatives, hi time the very simple Backward
Euler scheme is used.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
54
dFx , y , z
^
(3.14)
F
- F* "x , y , i
* x ,y ,z
dt
A/
Where At is the time increment, hi space,
(
dF_ _ v
dx
(3.15)
2
Ac
(3.16)
f
dF_
dy
F , -F
V '-7V
Ay
,
(3.17)
5F
&
F
, -F
^ ,V’^2
A2
,
Where Ax, Ay, Az are the space increments in the three cartesian coordinates. For
simplicity we define the following coefficient:
(3*18)
I—
i.e
£ l
—
>.y.*
2e-‘-J-k
■u —o->,j,t
■..At
2e -t +o- ■tAt
^£ i,j,k
And
(3.19)
At
g».y.*
C1
Hy.* =■
1+
i.e
C2
Hy.* =
2At
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
55
Which finally gives the FDTD update equations for the electric field in the three
components x,y,z.
i
i
,/i+—
,/i ~
H.\ \ - H .\ \
-'u+\,k
i
i
i /i +—
in+—
H \ * , - H y\ 2 ,
'^.*4
Ay
i
I
/5TJ. * , -is rj. * t
i,j,k + —
'i , j , k —
_7
J
J 1
+ Cj . k
,n + —
„l „ |»
Z_s.Jl*+l = C
- .-t Ev\
y\ij.k
y U ,j,k
'J '
//I
£ • .I 1r,7 ‘,* = C *»y»*
‘ £.1"
+C*'»y»#
-ii.y.jfc
i
in*—
*
i
in + —
-//I l
y\i-i,j,k
(3.20)
Az
I
(3.21)
ffj?
-ffj.
i* -,i,k
' i— ,j,k
iJ
■>
Ax
1
I
,n + —
H x\ \ - H x\ \
x|».j4*
x,<j4 *
Ax
(3.22)
Ay
And the update equations for the magnetic field in the three components x,y,z.
(3.23)
1
,n + —
,« — 1
A #t
A
ij *-lx
x
~ Ey\
y
l*
H x\. }. - H x\. 1 +
(3.24)
I
|rt+—
1/t—t A£u.
^ yky,*
U =^ U
'U,y.*+ //----rt,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
56
H •Uj,k
± .\= H ± .\ +
i
At
Mu*
Ay
- E , i~ U
(3.25)
Ax
These update equations represent the core o f the FDTD algorithm. But
absorbing boundary conditions (ABC) have to be used to truncate the computational
domain. Mur’s ABC and PML are explained in [11,9].
3.2.2
Bench mark results
The standard low pass filter simulated in [83-84] was used to calibrate our
code. The time domain results are shown in Fig. 3.1. We can see the gaussian pulse
propagation along the transmission line without any reflections (case (a)) and the
reflections at the input when the low pass filter is simulated (case (b)).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
57
0.5
-1.5
100
200
300
400
100
200
300
400
600
700
800
900
1000
600
500
Time daps
700
800
900
1000
500
0.5
!>
-0.5
-1.5
Fig. 3.1 Time domain response
The results were compared to standard circuit theory simulations using Series
IV. Results are presented in the following Fig. 3.2.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
58
-10
c/T-30
-40
—
—
FDTD
Series IV
-50
Frequency (GHz)
-10
—
—
>20
FDTD
Series IV
-30
-40
Frequency (GHz)
Fig. 3.2 FDTD bench mark results
A reasonable agreement is observed. This shows that the FDTD method gives
meaningful results. The discrepancy can be explained by the EM coupling occurring in
the structure. Circuit theory does not model the EM coupling [85]. It is believed that for
this reason the FDTD results would be more accurate [86].
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
3 JJ
Order of accuracy
This scheme is second order accurate in time and in space. This means that the
approximation of the derivatives is correct up to the second order. Higher order schemes
can be derived for both space and time. These schemes offer the advantage of better
accuracy, less numerical dispersion but have the drawback of requiring more CPU time
and complexity of implementation. In this section we will present the basic methodology
used to derive finite difference schemes based on Taylor’s expansion.
Let’s consider a function x j^ > u ( x ,t) . The Taylor’s expansion at x( +Ax is expressed
as follows:
(Ax)* d*u
24 dx4
+0
2
(3.26)
(Ax)3 d3u
(Ac)2 d2u
du
u fo + A r)! = u| +At.^=v‘
\
Ua
dx
a?l ,
6
dx3
N s]
At x - Ax we have:
«(jcf -A *)| =w|
v 1
\
(Ax)* d4u
24 dx4
( a i )!
du
- A c .^
+0N
,,
2
(3.27)
(Ac)3 d3u
6
**
dx3
J]
By substracting equation 3.26 to equation 3.27 we get:
(Ac)3 d3u
du
X iJ .
3
dx3
(3.28)
+0
N
s]
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
60
Which leads us to the finite difference approximation of the first derivative:
(3.29)
du
+0[(A*)’ ]
2Ax
dx
In order to get a scheme that is accurate at a higher order, we need to use points that are
further. Taylor’s expansion at x;. +2Ax and x, - 2 Ax are used:
4,
. / . , i d 2u
(3.30)
vj83«
+2(4t) a? +3(a,) a?l
2 /a \*d*u
a ? 4iAi
+o[(M)3]
And,
/
4,
v1 d2u
♦2M ^
(3.31)
, i d 3u
a?
+0[(Ax)J]
iiJ.
By substracting 3.30 and 3.31 we get:
(3.32)
du
«(*< +2A*)L - u ( x £+2Ax)|f< =4Ax.—
+0H
•*!**»
]
Multiplying equation 3.28 by 8 we get:
du
8u(xf +Ar)|f_ -8m (x; -A x)|f< =16Ar.—
8/ A
3^
\ 3 ^ 3“
>a^
+0
[<M!]
We see that we can now substract 3.33 and 3.32 and we obtain:
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
(3.33)
61
8«(xf +Ar)|f<-8 u (x f-A r )[ -w (x { +2Ar)|f< + « ( * -2Ax)|f< =
. . . du
12Ajc.-~
dx
(3.34)
+o[(A*)5]
Finally the fourth order accurate central difference scheme can be expressed as:
(3.35)
du
dx
J_
- § « ( * - H , + f “ ( ^ + i t )L
Ax
+ 2H
+0H * ]
Fig. 3.3 shows the relative error versus the number of discretization points
while computing the derivative of the sine function. We see that the slopes are different
which is in perfect agreement with the order of accuracy of the schemes. Backward and
forward up-wind scheme have been added and exhibit first order of accuracy.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
62
e
i
— 2nd Order
• Forward
i
— 4th Order
—
Backward
.-10
.-12
Number of points
Fig. 3 3 Order of accuracy of finite different schemes
33
The extended FDTD method
Conduction current can be added to Ampere’s law. This conduction current
can therefore be formulated to represent the existence of a lumped element within a
FDTD cell. In this chapter we present the theory o f the extended FDTD method [16]. The
derivation of the new update equations for the FDTD algorithm is performed. Study cases
include: A resistor, a capacitor and current source.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
63
3.3.1
The resistor
When a resistor is connected inside a FDTD cell, the FDTD update equation
needs to be modified to account for the current due to the resistor. These changes are
done in the direction in which the resistor is oriented.
The current flowing through a resistor is given by Ohm’s law, the voltage
across the cell in which the resistor is located is obtained by integration of the electric
field. In equation 3.36 the resistor is located in the y direction.
(3.36)
V _ d yE n
R
R
R is the value of the resistor implemented in the cell, the superscript n in the electric field
states the time increment.
The current density flowing through the cell is obtained from 3.36 averaging the electric
field for stability purposes.
i
rf,
dxdz
dxdzR
( r t f )
(3.37)
2
Ampere’s law is then modified to include the current density of equation 3.37.
(3.38)
In order to get an explicit formulation of the electric field at the new time step a few
manipulations are required. It leads to the new FDTD update equations that accounts for
a resistor inside the cell in the y direction.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
64
1—
E"+l =
1+
33.2
(3.39)
Aftfy
At
A/rfy
2Redxdz
En +Aftfy \
-V x ff
e\ l+ IRedxdZj
2Redxdz)
Parallel resistor and capacitor
The methodology is the same than in the previous section. But here the current
is due to a parallel resistor and capacitor.
. V „dV
/ = —+ C ——
R
(3.40)
dt
In terms of electric field and current density, equation 3.40 is expressed as:
_dy(E n + Ea+l) Cdy(En* - E n)
2 Rdxdz
(3-41)
dxdzdt
We set the following parameters to simplify the equations:
dy
-— +
IRdxdz
o_ dy
IRdxdz
a-
—
Cdy
—
dxdzdt
Cdy
dxdzdt
(3.42)
—
This allows us to write Ampere’s law as follows:
En+i = En + — V x / / - — o r r +1~— 0E n
s
e
e
After manipulations of the expression in order to get explicit formulation, we get:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3'43)
Equation 3.44 represents the modified FDTD update equations of the electric field when
a parallel RC circuit is inserted inside one cell.
3JJ
Parallel resistor and current source
To implement an equivalent circuit of a transistor we need to modify the
FDTD update equations to simulate the current source of the drain side. A resistor in
parallel with a current source can be incorporated inside a cell as follows:
The current flowing through the cell is presented in equation 3.4S.
V\
f.
n SnYgs
J = - ---------dxdz
(3.45)
Where Vgs is the voltage at the gate side computed by integration of the electric field.
Ampere’s law is then modified, by including the current density of equation 3.45.
En+l = E n + — V x f f — & d y _ t £ n+l + En) ~ — — Vgs
e
2eRdxdz
£dxdz
(3'46)
In explicit form this lead to equation 3.47.
^
Atdy
"j
Ar
rntl _V 2R£dxdZj pit .
~ fu
tody y
f
t
IRedxdz J
I
£
tody
(3-47)
Hy u
\
2Redxdz)
^ Sm
e d x d z fu
I
Ardy
2 Redxdz.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
66
3.3.4
Multiple cell formulation
In the previous sections, all the update equations represent the implementation
of lumped elements in one cell of the finite difference grid. Due to the physical
dimensions of the device, one may want to distribute these lumped elements over several
cells. The main problem of this formulation is that it leads to implicit schemes difficult to
solve [87-88]. To fully deal with the physical dimensions of the lumped element we
would have to modify the schemes presented in the previous section to account for the
voltage across the number of cells needed to satisfy the physical length of the element.
The current would have to be written as:
(3.48)
Where j spans from m to p to cover the cells enclosing the resistance. This transforms
Ampere’s law to:
(3.49)
edxdzR
Unfortunately this scheme is unstable because the electric field is not averaged in time.
One can decide to average in time only the updated electric field component. This would
result in the following update equation:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
67
'Z E j
edxdzlR
^
edxdzlR J
£
edxdzlRj
(3.50)
e dxdzR
edxdz2R
This scheme is stable for high permittivity as explained in [88]. It is easy to
understand that averaging in time every component of the electric field that are involved
in the update would result into an explicit scheme. This new scheme needs to solve a
system of linear equations to update the electric field [87]. This would be the ideal way of
implementing lumped elements through multiple cells.
3.3.5
Simulation results
In our approach [35], we used a much simpler technique. If the lumped
elements inserted are in the y direction, the actual value o f those parameters per cell is as
follow:
(3.51)
Where nx and ny are the number of cells below the microstrip in the x and y direction.
As an example, a resistor loading a microstrip line is simulated. The results are shown in
Fig. 3.4. We see that the reflection coefficient for the four different values o f resistance is
approximately the one that we obtain theoretically. The frequency response is relatively
flat. The errors are due to numerical implementation and capacitive effects that are
explained in [89] and could be corrected by a technique presented in [90].
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
68
R=10
§ -1 0
-1 2
-14
-16
-18
-20
Frequency (GHz)
Fig. 3.4 Resistor modeling using Extended FDTD
3.4
Principle of neural networks
Artificial neural networks (ANN) are non-linear approximators for multi­
input, multi-output functions. They are based on a training process that mimics our
understanding of how the human brain works. A set of training data is used to optimize
some weights. Then a set of testing data is used to verify the validity o f the non-linear
approximation. In this section we present the basic theory o f neural networks and we
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
69
derive a comprehensive matrix formulation for easy implementation. Fig. 3.S shows a
three-layer neural network.
Cgd
Vdi
ELdi
Cdi
Input Ujur
KiddnUjwr
Oupu»Uj«t
Fig. 3.5 Three-layer artificial neural network
There are two inputs - Gate to source voltage
V*
and drain to source voltage
and five outputs - the drain current Ids, the gate to drain capacitance Cgd, the drain
to source capacitance Cds, the gate to source capacitance Cds and the drain to source
resistance R&.
The input of each neuron in the hidden layer is the sum of all the outputs of
the input layer multiplied by a set of weighting factors. The output is given by a so-called
activation function, which is a non-linear function applied to the input in order to
correctly model our data. The adaptability o f the network is defined as the discrepancy
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
70
between the outputs given by the network, and the desired outputs. The computation of
the outputs can be seen as a simple multiplication of matrices.
Let’s consider a general three layer neural network, with / input, m hidden
neurons and n outputs. Let Xu, be a set of input values, the weighting factors can be stored
in a matrix
.
'
(3.52)
in
*\
X"
•
•
•
* <
<
'
•
r in
*1-1
„HI
<>
. <.
Then, the input of the hidden neurons is a vector
computed as the multiplication of X,„
and Wu,.
(3.53)
The output of the hidden neurons is a vector lm obtained by applying the activation
function to the vector lin. We chose to use an hyperbolic tangent function shown below.
2
.
—1
=uig{lin) w i* tsig{x)=l+ e x p ( -2.x)
(3.54)
Finally, the output of the neural network is obtained in the same way than the input of the
hidden layer. You multiply the hidden layer output vector Iout , and a weighting factor
matrix Wo u .
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
This procedure computes the output of the neural network for a given input vector.
During the learning process this has to be done thousands o f time. For a given input
matrix (representing the set o f input vectors), let’s define the set of data used for training
'
' jT ts t
jT r a in
d\
and the set of data used for testing as DTat =
as DTmin =
I -1
jT r a in
An
{djmin ,d]“‘
. Each
.
is in fact a pair of vector containing the output data of the i h output for
different input vectors.
The network is trained, when the error defined by equation 3.56 is less than a chosen
error bound e-iram.
E rm ^
* i=l
- < ‘ 't
° X)
This training is validated, when the error defined by:
(3.57)
E rro rs = ± £ (< /,r”
1=1
Is less than a chosen error bound Exest-
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
72
A set of data obtained from the steady state solution of the hydrodynamic
model is used to train the network. The data are normalized to speed up the convergence;
and the Levenberg-Marquardt method is used to update the weighting factors. The
number of neurons in the hidden layer is optimized to avoid over training, and to give a
good accuracy for testing data. Fig. 3.6 presents the average error for different number of
neurons in the hidden layer. It can be seen that for 6 neurons in the hidden layer, the RMS
error is the same for training, as for testing.
If much more neurons are needed to accurately model the non-linearities, a
multi-layer neural network can also be used [30] to reduce the number of neurons.
Tatting
ui
IS
Numbar o f Hidden naurana
Fig. 3.6 Error curves for training and testing data.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
73
3.5
Hybrid ANN-FDTD method
The extended FDTD method and an ANN are used to simulate an amplifier.
Physics based simulation of a MESFET provide a set of data to train an ANN. Then
lumped elements representing a transistor are introduced in a FDTD grid. Large signal
simulation can then be performed without the need to stop the FDTD marching time for
simulation of the physics based MESFET.
In this research we restricted ourselves to an intrinsic transistor model due to
the results obtained from the hydrodynamic model and for simplicity of implementation
in the FDTD mesh. Fig. 3.7 shows the basic amplifier that is to be simulated using the
extended FDTD approach.
t opyt Mi fcl i unNwcft
Out put M t e t n o n Nh t wo r t
Bit
n /\H
Ua
TbBiMtEqimlMrtcttcial
Fig. 3.7 Basic amplifier with matching networks and intrinsic transistor model
Parallel and series lumped elements can be inserted in a cell by adding the
current density flowing through this cell, due to the presence of the lumped element.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
74
Single loads have been investigated and the derivation was presented in section 3.3. Here,
we present the implementation of a parallel resistance and capacitance with a current
source. This represents the drain to source resistance and drain to source capacitance plus
the drain current source of the drain of a transistor. Each lumped element is distributed in
each cell, hence the value of the actual lumped element per cell is not the total value of
the lumped element. Fig. 3.8 presents an insight in the meshing of the structure at the
location of the transistor. This represents only a section of the actual structure, it repeats
over the number of cells under the transmission line.
v*
bcbQ*
Fig. 3.8 Equivalent circuit parameters of the transistor inserted in the FDTD grid
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
75
In the case of a current source in the y direction, we can derive the following
formulation to implement the drain of a transistor inside a FDTD grid. The drain current,
and the drain current density are given by:
h = 8 myas
(3.58)
jd = 8-*'Vaf
Ax.Az
For a parallel resistance and capacitance we define the two coefficients Coeffl and Coeffl
as:
C W -i
4 £ flZ _ + _ £ 5 z _
2.erJ lA x A y erA xA z
0 * 2.,♦ _ * & _ + _ £ & _
2.£r.RA xA z £rAx.Az
<159>
© *>
Which leads to the following update equations:
Ar
„ in+i
Coeffl _ i"
EA
= — — *£„
’w
C o effl
Those coefficients
When we
(3.61)
£r
_ rr i«4
+ — -— * V x / / v 1
C o effl
y|'^*
are derived from a semi-implicit scheme used to avoidunstability.
add the current source, it is straightforward to obtainthefollowing update
equation:
AT
.~i = Coeff 1
C o effl
i- + _ £ r _ * v x f f H
y^ - k C o effl
(3.62)
+—
— £ I"
erA x M £ o e ffl y|w ^
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
76
Where
Vgs
is computed as the integral of the electric field from the ground
plane to the strip. In the previous equation, we assumed that the gate, and the drain of the
transistor were separated by p cells. Thus, the gate voltage is calculated at the k-p cell in
the z direction.
The values of R* , Cds» gm at the output, and Cp at the input, are updated
according to the gate, and the drain voltages obtained from the electric field, thanks to the
ANN model. This leads to a first order approximation of the device-wave interaction, and
to a large signal analysis.
3.6
Full-wave simulation results
A typical MESFET was simulated according to the hydrodynamic model
presented in the second chapter. Then, a three layer neural network was trained to model
the non-linearities of Qgs, Rjs, Qd» Qgd and gm. Six neurons were used in the hidden layer.
Fig. 3.9 shows the I-V characteristics of the simulated MESFET and the results obtained
after training of the ANN. The agreement for both sets of data used for training-testing,
and the hydrodynamic model initial data is very good.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
77
i -v a u n c m ta lc i
ffiOr
T
T
T
jjj- # - •<
...
200-
^
i
*
d
*
:
:
o
6
2
A ..I
ISO-
* •
+
<s»
8 *
•
100-
j 'T
;
:
a
*
*
1 - * # ! *
SO
f
K
r
r
t e f8 T
*
« 8 ?
T r i r - - M
-t*
0
I
-SO
OS
1 .6
2.6
;
3.6
1
J
d a ta fra in e d
d a ta te a ta d
m odel
;
4
i
4.6
Vdt
Fig. 3.9 I-V characteristics of the simulated MESFET, hydrodynamic model and
ANN.
The accuracy is better than 1%. At a given Vp, the drain current can be
computed for the all range of
almost instantaneously. This has to be compared with
the 45 minutes necessary to compute 21 points using the hydrodynamic model on the
same computer. Hence, an operating point can be predicted using the ANN without
requiring to solve for the hydrodynamic set of equations.
The ANN can be used to calculate the transconductance, gm. Fig. 3.10 shows
the transconductance for three different V* values. The agreement is very good, which
means that the non-linearities of the MESFET are well modeled using the ANN. Fig. 3.11
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
78
shows the gate to source charge that leads to the capacitance of the equivalent circuit
presented in Fig. 3.7. The neural network predicts very well the non-linearities of the
charges. The same results are obtained for the drain to source charges.
Tranaconductance
200
160
100
E
IE
£at
-50l
-2
-
1.8
-
1.6
- 1.2
-
0.8
Vga
Fig. 3.10 Transconductance computed from the ANN results
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
79
XMncr" i
v
O p of a MESFET from Hydrodynamic m m and ANN
• * •
XA “
x.%-
9
i * *
9
1.1
.....
a - * -
• :
^ . r r .r .*.
• * ......................... - ..................—
i • *
........ .........
-■
• • •
• •
+
♦
0
as -
1
18
2
..............
i - i • * • * •
!
• •
i . • • ?
r • .
as
—
28
3
Fortaining
For toting
ANN
38
4
48
8
VU»
Fig. 3.11 Gate to Source charges, Q gs - Hydrodynamic model and ANN results.
I1 Q - “
Q d tto d n ln o tia s tin c a
U p
848
•
4-
38-
9
3-
28 2T8 -
la
-1 .8
-1 8
-1 .4
-1 8
-1
-OB
-a s
-0 .4
-0 ,2
0
Vs .
Fig. 3.12Gate to source capacitance, Cgs computed from the ANN results.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
We simulated the transistor without any matching network, to validate the
method. The dielectric constant of the substrate is £(=2.2, and the height of the substrate
is
0.794 mm. The amplifier is assumed to be bias at Vgs=-0.6 Volts and
V d s= 3
Volts.
At this bias point, Cs.,=0.5 pF, C<ts=l7fF, Rds-120 Q and gm=I20 mS. These results are
compared with a simulation performed using HP-Libra, the trend of S2 1 , and 5/y is close
to the FDTD results as shown in Fig. 3.13 and Fig. 3.14. Fig. 3.15 and Fig. 3.16 show the
time domain results of a small signal FDTD simulation.
*■ - ►_ ^
S21 :
-10
Libra, HP-EESOF
Extended FDTD
-15
Fig. 3.13Small-signal S-parameters comparison between the FDTD method and HPLibra for the transistor. (Dashed lines: FDTD, Solid line: Libra)
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
81
Libra. HP-EESOF
Extended FDTD
Fig. 3.14 Small-signal S-parameters comparison between the FDTD method and HPLibra for the amplifier. (Dashed lines: FDTD, Solid line: Libra)
O utput
100
200
300
600
600
700
aoo
800
1000
Fig. 3.I5Time domain response of the small-signal amplifier
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
82
0. 2 :
-
0.1
-0.4
-0.5
100
120
Number of call*. Z n i t (fli*0.4333 mm)
Fig. 3.160utput voltage versus distance in the z-axis at three different times.
The structure is excited by a gaussian pulse (dashed line in Fig. 3.1S). The
total input voltage is represented by the solid line, whereas, the dashed-dotted lines
represent the output total voltage at two different positions in the z-axis. We can obsreve
at the output, that the incident wave has been amplified, and that it travels in the +z
direction. It is to be noted that the output is perfectly matched with a PML wall. Fig. 3.16
demonstrates the same behavior, but in this case, we have the total output voltage versus
z. The structure is divided in 111 cells in the z direction, with dzpO.4233 mm. At t=0.88
ps, the incident wave has not reached the transistor yet, but at t-2 .7 6 p s, it goes out of the
transistor. Then, we can notice that the wave is amplified, and also that a discontinuity is
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
83
observed at n-55, which is the location of the current source. At t=2.64 ps, the wave is
out of the active region, and it is widely spread over the z-axis. A reflected wave is also
observed.
Then, the neural network has been used to take into account the non-linearities
of the intrinsic parameters. The structure is excited by a sinusoidal wave of amplitude
, and frequency f- 1 0 GHz. Fig. 3.17 presents the output voltage for three different cases.
As the signal gets larger, the wave starts to saturate, generating harmonics in the
frequency domain. Fig. 3.18 shows the output voltage in the frequency domain.
a is
IBn >0.3
|Ein
|Bn •1.5
0.1
006
-
0.1
- 0.2,
100
200
300
700
900
400
GOO
Tima »«ratiora, d M .44ia-i 2 1
BOO
900
1000
Fig. 3.17Time domain response of the large-signal simulation
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
84
OJ
08
08
08
03
Fig. 3.18Harmonics obtained from a large signal simulation using the ANN.
In each case, the output voltage has been normalized to its maximal value. We
can observe that at small signals, no harmonics are generated, but as expected, second
and third harmonics are obtained for the two other cases.
3.7
Summary
This chapter presented a new technique for global simulation of active
microwave circuits. The neuro-FDTD method uses an artificial neural network to model
physics based simulation results and update lumped elements values inside the FDTD
grid. This technique has been presented in [34-35] and was recognized as a state o f the art
method in [91].
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 4
ON THE USE OF WAVELETS FOR GLOBAL MODELING SIMULATIONS
4.1
Introduction
The principle o f multiresolution is very interesting for global modeling
problems because they exhibit several levels o f resolution. Therefore there is a need to
develop a wavelet-based technique to tackle such kind o f simulations.
In this chapter we will introduce the basic theory o f wavelets in section 4.1,
then in section 4.2 we will present a thorough review o f the current state o f research in
wavelets for electromagnetics problems. Section 4.3 will develop the use of interpolating
wavelets as a non-uniform mesh generator for semi-conductor modeling. Finally in
section 4.S the fast wavelet transform is used to solve efficiently the Poisson’s equation.
Wavelets represent an orthonormal basis o f square integrable functions. It lies
within the multiresolution framework. The lowest resolution can be characterized by the
scaling function that filters smooth variations o f the chosen unknown. Higher resolutions
can be added by filtering details using the wavelets functions, hi this section we present
the basic mathematical formulation o f wavelet theory. This will help the reader to get
familiar with wavelets and understand how it can be used to develop efficient numerical
method.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
86
4.1.1
I 2(ft) and a few definitions
For the one-dimensional case, we look for a solution o f a PDG in I 2(ft)
which is the vector space o f measurable, square-integrable one-dimensional functions.
The inner product associated with this vector space is given by:
V /,g e Z 2(ft)
(4.1)
(s(“ ) / ( 4 =
) / ( “ )*«
The norm o f a function is given by:
V/* 6 Is (ft)
(4.2)
I f f = ( / > / ) = l / ( “ ) 2rfM
-0 0
The convolution o f two functions is given by:
V / , g e I 2(ft)
f*g(x)=
(4-3)
(/M*s(“)XO
/ * s(x)= * }/(« > (* - «V«
-oo
And finally the Fourier transform o f a function is given by:
V /e Z 2(ft)
(43)
/(<»)=
-<o
4.1.2
Multiresolution approximation
We arelooking for an approximate solution o f our problem
projection onto
by orthogonal
a vector space, which is a subspace o f I 2(ft). Given
a specific
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
87
discretization o f a domain, the solution o f a problem on this particular grid represents an
approximation o f the solution at a particular resolution.
This solution at a given resolution represents an orthogonal projection on a
, where V is the resolution. This vector space
vector space that we will call
represents the set o f all possible approximations at the resolution 2j o f functions in
1} (9l). If one uses a finer grid (i.e a higher resolution) the solution on the coarser grid
will be also solution on the finer grid, which means:
V /e Z , Vv a V ^ c ^ . . .
(4.4)
As the resolution increases to infinity, the solution tends to the exact solution
and all details are added. Conversely, as the resolution decreases to zero, fewer and fewer
details are included and the solution converges to zero, hi a mathematical formulation we
have:
lim V,j =
/-♦ -mo
is dense in £2(ft)
jm
(4.5)
-<0
And
Um Vv = f V , , = M
Such a set o f vector spaces {v^ )
^
is called a multiresolution approximation o f I 2(9l).
To characterize a function (i.e the solution o f our problem) at a given resolution, we need
to find an orthonormal basis o f Vv . Such a basis can be defined by dilating and
translating a unique function x -> ^ (x ).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
88
Theorem
Let (F,y)
be a multiresolution approximation o f L2(SR). There exists a unique function
x - > 0 ( x ) e £ 2(9t), called a scaling function, such that if we set <f>v {x) = 2J ifa 1x) for
y 'e Z , t h e n ( j c - 2~Jn )]^ is an orthonormal basis o f Vv . In other terms,
y ¥ ^ x - n ) \ . z is an orthonormal basis o f VzJ
The term
appears in order to normalize the functions with respect
to the £2(9t) norm. It is to be noted that there exist such a family o f scaling
functions for a given multiresolution approximation and that for a different one,
the scaling functions would be different.
Now that we have an orthonormal basis for Vv , we can expand the
orthogonal projection o f a function x -> f{ x ) on this subspace as a summation o f the
basis functions namely:
Let’s define A^ , the orthogonal projection operator.
Then,
V / e L 1fy),A lif(x)= 2 -i
n ) ^ , ( r - 2 ' y/i)
^
rtte-oo
The basis functions are known, thus to determine the functions we just need to determine
the set o f inner products.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
89
The idea behind multiresolution expansion is to characterize a signal at a
given resolution and add the detailed signal that leads to the next higher resolution level,
hi other words, we want to compute the difference between the projection on V2, and on
V,j„. It can be shown that the detailed signal is obtained by orthogonal projection on the
orthogonal complement of P,, in P,,.,: Wv .
(4.8)
To compute the detailed signal we need to find an orthonormal basis o f Wv , like we
needed to find an orthonormal basis o f P,; .
Theorem
Let (P,; )
be a multiresolution approximation o f I 2(SR),
x
-» 0 (x )e L2(9l) the scaling
function associated with this resolution, and H be the corresponding conjugate filter (cf.
Appendix A). Let x -> y{x) be a function whose fourier transform is given by:
with, G(eo) = e~'“ I f ( a + jc)
Let \ffv (x )= 2 yV (2yx) for j e Z , th e n (V ? V 2J ( x - 2 ~JnjjneZ is an orthonormal basis o f
W2j . And, y & y r f e 'x - n j l ^ j ^ i is an orthononnal basis o f
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
90
To compute a wavelet, we define a conjugate filter a -> H{ca), then compute
the corresponding scaling function and finally according to the last theorem we define the
orthogonal wavelet. Depending upon the choice o f the conjugate filter, both scaling and
wavelet functions can have good localization in spatial and Fourier domains. Daubechies
[29] studied the properties o f the scaling and wavelets functions depending on the
conjugate filter. She came up with special wavelets, o f compact support and n times
continuously differentiable. Other wavelets may not have compact support: BattleLemarie wavelets [92], but have been used extensively recently due to their fast
computation compared to Daubechies wavelets.
As explained before in the case of Vl t , to compute the detailed signal we want
to expand the projection o f the function / i n Wv as a sum o f all the basis functions,
namely:
Let’s define P2i, the orthogonal projection operator in Wl t .
Then,
V / € L2(< R )^ ,/(x )= 2~J £ ( / ( « ) {t - 2~Jn)yr2j ( t - 2~Jn )
(4'9)
/!*-«
The inner products to be determined are the wavelet coefficients.
All of the previous discussion represents the basics o f wavelet and
multiresolution theory for the 1-D case. For other dimensions, the extension is straight­
forward but a very important theorem very useful for engineer must be mentioned: This
as to deal with the particular case o f separable multiresolutions.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
4.1.3
Separable multiresolution
In most practical problems, a three dimensional formulation is needed. Let’s
consider a three dimensional case, for a separable multiresolution approximation o f
Z,2(iR3). Each vector space V t can be decomposed as a tensor product o f three identical
subspaces o f L2(9l).
Vv =V[i ^ V \ i ® V\i
(4.10)
*s a nwltiresototi011 approximation o f Z,2(iH3) if and only if (k,1, ) ,eZ is a
multiresolution approximation o f L 2(SR). The scaling function (x ,y ,z)-> 0 (x ,y , z) can
then be written as:
fay*z)=j(x)Mfe)
(4 U )
Where x -» <j>(x) is the one-dimensional scaling function o f the multiresolution
approximation (k,1, )
.
Therefore the projection of a signal/at a resolution 2J is characterized by the set o f inner
products:
K f =
( x - 2 - Jn)t2/( y - 2 - Jm )t2J( * - 2^ ) ) ) ^ ^
(4-12)
Another extension due to the separable property states that an orthonormal
basis o f W2j (The orthogonal complement o f V2t) is obtained by scaling and translating 7
(2n-l where n is the space dimension) wavelet functions:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
92
(x,y,z)-> yrl(x,y,z)
(4.13)
( x , y , r ) - n / 2(x,y,z)
( x ,y ,z ) - » ^ 3(x,y,z)
fa ,y ,z ) ^ \f/* fa ,y ,z )
(x ,y ,z )-» y r5(x,y,z)
(x ,y ,z ) - > ^ 6(x,y,z)
(x ,y ,z)-> < /7(x,y,z)
From this, the following theorem can be said:
Let
(^,y)JgZ be
a separable multiresolution
approximation o f l
Let
(x, y , z) -> ^(x, y , z) = <j>{x)^{y)ift{z) be the associated three dimensional scaling function.
Let x -> y{x) be the one dimensional wavelet function associated with x -> 0(x). Then
the seven “wavelets”
(x ,y ,z)-> </1(x ,y ,r) = ^(x)^(yV (^)
(x,y, z)-> <f/2fa,y, z )= 0 (x )y (y )t(z )
(x ,y ,z)-> y 3(x ,y ,z)= if/(x)fi(y)fi(z)
(x ,y ,z)-> if/* fa, y ,z )= tfa )y (y )ff(z)
^ f a y , * )-> v 5fa ,y ,z )= v fa fy fa ty fa )
(x ,y ,z)-> \}f6fa y ,z ) = y/fatyfatyfa)
(x ,y ,z)-> if/1f a y , z)= tf/fa)f/(y)f/(z)
are such that,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4 l4 )
93
(x- 2~Jn ,y -2 ~ Jm ,z - 2~sl \
(
2~VJ
2~Jn ,y
(4.15)
~2~Jm, z “ 2 ~Jl )
2"Vi> ( ? - 2 'j n ,y - 2 - J m, z - 2 - Jl )
2~J ip*, ^ c - 2 'Jn , y - 2~j m , z - 2~Jl )
2~Jy \, (r - 2 'j n ,y -2 ~ j m, z - 2~Jl )
2~Jy/\, ( r - 2~Jn ,y -2 ~ Jm , z - 2 ~J/ )
2 " V i ( ? - 2 - j n , y - 2 - j m ,z - 2 - Jl ) \ nl>z3
is an orthonormal basis o f WnJ.
And,
( 2
(x -2 ~ Jn ,y -2 ~ Jm ,z -2 ~ Jl \
2- > ; ( x - 2 - Jn , y - 2 - Jm , z - 2 'Jl )
2 - ' r l (? - 2~j n, y - 2 ~ Jm ,z -2 ~ Jl )
2’V i (* - 2~J n, y - 2~J m, z - 2~J/ )
2~ V 25y(x -2 ~ J n ,y -2 ~ J m ,z -2 ~ J/ )
2- > :6,(* :-2~J n ,y -2 ~ J m ,z -2 ~ J l )
2- > 27^ :- 2~J n , y - 2 ~Jm, z - 2~Jl ) )<
is an orthonormal basis o f /.'(w 3).
This theorem represents the basis o f the multiresolution time domain method
[SO] that will be introduced later in this chapter. The same can said about four and more
dimension, so time dependent functions can be expanded using wavelets in both space
and time dimensions [93].
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
94
4.1.4
Two examples: the Haar and the Battle-Lemarie wavelets
Wavelets can be defined by closed form formulations or by iterative relations.
The Haar wavelets are the first of the Daubechie wavelets and represent the pulse
functions. The Battle-Lemarie wavelets have a closed form expression and are therefore
easier to handle than Daubechie wavelets of high order. We will now present those two
basic wavelets to give an idea of what they look like.
The Haar scaling and wavelet functions are pulse functions defined by:
l,/ o r |s |< l / 2
(4 .1 7 )
l /2 ,/o r |i| = l / 2
0 , otherwise
1.3}
i
t
-a t
-a t
—■,
H a n rS o a ln g fu n e tD n
t
,
-i —
i
»
i
t*
as•
0
i«l
-t
■ i
-a.4
i
-02
»
o
i
02
Spn—lorTlwjQoitnim
■>
a*
*■
at
»
at
- .1
i
Fig. 4.1 Haar scaling function
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
95
And,
(4.18)
112, fo r ts = —1/2
l,fo r : - H 2 < s < 0
0 = - l ,fo r : 0 < s < 1 /2
- 1 / 2 , f o r ts = 1/ 2
0, otherwise
1.5
i'"
i
i
-0.6
-0.4
Haar WivtW function
i
i '
l
l
i
r
0.4
0.8
0.8
0.5
-0.5
''- I
-0.8
-0.2
0
0.2
Spatal or Tim* coonMrmm*
1
Fig. 4.2 Haar wavelet function
We see from Fig. 4.1 and Fig. 4.2 that the scaling function and the wavelet are
orthogonal and normal. By translating and dilating these functions we can create a basis.
Unfortunately, these wavelets are not differentiable which make them of little use to
solve PDEs.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
96
The Battle-Lemarie wavelets presented here, are the cubic-spline wavelets,
they are smooth and continuously differentiable. We can find a closed form expression of
the scaling and wavelet function in the frequency domain:
->'*
.
4
1—
sin(*x /2 )
* > : J 5 _ (M 2)
s •i n —
3
4
315
x —sin
2 •
+
U J
.
sin
5
(4.19)
(—
kA *
U ,
12
Fig. 4.3 shows the Battle-Lemarie scaling function, while Fig. 4.4 shows the BattleLemarie wavelet function.
Battle-Lemarie scaling function
1.4
0.8
0.6
0.4
0.2
-10
Fig. 4J Batlle-Lemarie scaling function
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
97
Equation 4.19 represents the closed form expression of the Battle-Lemarie
scaling function and equation 4.20 represents the closed form expression of the BattleLemarie wavelet function. Both of these expressions are in the frequency domain. Fast
fourier transform (FFT) needs to be performed to compute these two functions in the time
domain.
j— Slk +2ft) -/
>
(4.20)
4(^12+ ft)
BaMa-lamariawavaiat function
0.9
0.8
0.7
0.6
0.4
0.3
02
-10
Fig. 4.4 Battle-Lemarie wavelet function
From the two previous graphs, one can see the low-pass property of the
scaling function and the band-pass property of the wavelet function. In most papers, one
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
98
chooses a starting resolution given by particular space mesh densities. This resolution is
said to be the zero resolution. It represents the vector space at which the first
approximation of the solution is made: Vv . Sometimes the fields are so smooth that the
expansion using only the scaling function is enough. Sometimes however, we may need
the next resolution level, this can be done with few effort by expanding the fields on the
wavelet basis, where it is needed. Thus we add the detailed signal, which is a vector of
w2l. By doing so we construct the approximation at the next resolution level: v2,„ .The
scaling functions are orthogonal to each other at a given resolution level, but we cannot
mix scaling functions of different resolution level. On the contrary the wavelet basis is
orthogonal, thus we can expand the fields using the wavelet functions at the same and at
higher levels if necessary. This means that we are looking for an approximation in the
vector space:
vv ®wv ®
4.2
®w2„z ®
-
(4.21)
A review of wavelet research
The main papers of the past ten years are presented. We focus on the
applications of wavelets in the method of moments (MoM) and in FDTD. Then we
broaden the subject with a list of paper that gives the trend toward efficient numerical
solution of PDE’s. The different approaches are presented, wavelet-based Galerkin
method, wavelet collocation method, and post-processing using wavelet preconditioning.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
99
This section gives a good understanding of applied wavelets techniques with multiple
references that will highlight some points developed in the previous section.
In the late 1980’s a new class of functions appeared: The wavelets. These
functions are highly localized both in space and time and lie into the construction of
multiiesolution analysis to create an orthonormal basis of square integrable functions.
The principles of muitiresolution analysis are explained by Mallat in [44], a more recent
reference is [94]. A signal can be approximated by his projection on a subspace
corresponding to a given resolution using the so called scaling function. To obtain the
approximation of the signal at the next higher resolution (scale), you need to add the
detailed signal that can be expressed in terms of wavelets. The two subspaces spanned by
the scaling function and the wavelets are orthogonal. One may want to construct such
scaling and wavelet function to satisfy the orthogonal relations. Different kind of
wavelets have to be considered. The compactly supported wavelets and the noncompactly supported wavelets. In this two categories one can choose different style of
wavelets. The Battle-Lemarie wavelets are constructed from B-splines functions of
different order, they have infinite support and they have closed form expressions in both
the time and frequency domain [95]. The famous Daubechies wavelets [29] have compact
support and are obtained by an iterative algorithm on the dyadic grid. One special case of
Daubechies wavelets are the pulse functions, they are called the Haar wavelets, they have
sharp discontinuities due to their zero vanishing moment property. One can also consider
the Meyer wavelets, that are constructed in the transformed domain, satisfying the
orthogonal relations in the frequency domain with functions of compact support in the
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
frequency domain. These Meyer Wavelets can also have closed form expressions: Rised
cosine filter. Other wavelets can be considered: Coifman wavelets, periodic wavelets and
intervalic wavelets. The sets of functions that we talk about so far, form an orthogonal
basis of square integrable functions. They are defined on the line. So what happens when
one deals with an interval ? The completeness of the basis is broken. This gives just the
idea that wavelet decomposition on an interval has to be treated with care. Periodic
wavelets can be easily constructed on the interval and satisfy periodic boundary
conditions. When non-periodic boundary conditions have to be treated, we must use
intervallic wavelets such has coiflets.
Armed with all these tools Mathematicians tried to tackle engineering
problems. The most fruitful use of wavelets is in image compression and signal analysis
where a wavelet decomposition is seen as a convolution with a low-pas filter and a highpass filter. These filters named filter bank can be found in [29] for Daubechies and
Coifman wavelets. Because of the control on the number of vanishing moments,
electromagnetics engineers realised that using a wavelet expansion in the method of
moments may result in a sparsified admittance matrix [33], thus saving memory
resources and CPU time during the solve of the linear system equation. An approach that
is philosophically opposed to this one is to treat the moment method matrix as we treat an
image in signal processing. By applying different Fast Wavelet Transform (FWT) one
can choose the best wavelet bases to compress that matrix and solve the linear system in a
dramatically reduced computation time [96]. This represents the idea explained in [97]
where it is demonstrated that a wavelet preconditioning can reduce the condition number
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
101
to 0(1). This technique is used in section 4.5 to implement an efficient solver for
Poisson’s equation applied to a standard FET structure.
The biggest achievement in wavelets applied to electromagnetics is with no
doubt the one of Krumpholz [46]. He demonstrated that finite difference time domain
algorithms can be derived from the method of moments applied to Maxwell’s equations
using wavelet functions as expansion and testing functions. The resulting dispersion is
greatly enhanced [98] and overall, time and space adaptivity can be obtained. In [46], it
was demonstrated that using Haar wavelets one can derive the Yee algorithm of the
regular second order accurate FDTD. If one decides to use higher order Daubechie
wavelets it is expected that higher-order schemes will be obtained. This is in agreement
with the enhance linear dispersion. The main problem is that Daubechie wavelets do not
have closed form expressions thus, Battle-Lemarie have been employed. But on the
interval it results in tremendous problems to implement correctly the absorbing boundary
conditions.
Wavelet-based technique seem very promising for linear problems, but in
engineering we often deal with non-linear equations. Mathematicians have been working
on wavelet techniques to solve PDE’s since the beginning of the decade. They all seem to
agree that a pure wavelet expansion is too time consuming when dealing with non-linear
equations [99]. Instead they propose to use the auto-correlation function of Daubechie
wavelets. This function is known as the interpolating wavelet even if it is not really a
wavelet. It comes from the first stage of the lifting scheme that creates second generation
wavelet from the “lazy” wavelets namely Kronecker delta functions defined on a dyadic
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
102
grid [100-101]. Differentiation and multiplication is simplified in this wavelet basis
because wavelet coefficients represent the difference between the exact solution and the
interpolated values from the coarser grid at dyadic grid points.
In the following sections we will present the main papers published in the
literature concerning wavelets applied to EM analysis and the numerical solution of
PDE’s. The first part of this work will deal with the Multiresolution Time Domain
Analysis called MRTD algorithm. In the second part will present some results on wavelet
collocation technique and its different applications to Maxwell’s equations and semi­
conductor modeling. In the third part, we show that Daubechie wavelet expansions
represent a mesh refinement, which demonstrate that the best way to use wavelet is to
generate a non-uniform mesh. Finally we will present the work o f G. Beyikin who
worked on the discretization of operator in wavelet bases and wavelet preconditioning for
efficient iterative solver.
4.2.1
The MRTD algorithm
The MRTD algorithm originally presented in [46] is currently under a lot of
investigation [102]. The new edition of the book from A. Taflove on time domain
electrodynamics includes a chapter on MRTD [92]. It represents a good introduction to
those who would like to get a good start with MRTD.
The algorithm is derived from method of moments applied to Maxwell’s
equations with wavelet expansions of the fields. Problems from 1-D to 3-D have been
investigated. The major restriction so far is in the number of resolution levels. The
formulation is quite advanced and requires a systematic notation used in quantum
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
103
mechanics. The authors considered the S-MRTD algorithm, which is a basic expansion of
the field using the scaling function and no wavelets. But the W-MRTD has of course also
been used. It is to be noted that in the 3-D case the W-MRTD formulation is so big that it
cannot be written in a book. The authors only presents the Wy-MRTD, which is the
expansion of the field in the y-direction using wavelet [92].
The Battle-Lemarie wavelets are used because they have closed form
expressions. But they also have infinite support, which means that they have to be
truncated to satisfy the boundary conditions. The problem of the boundary condition is
not clear in those publications. They must be enforced on the wavelet coefficient and not
on the total field, so symmetric wavelets make it much easier. But what if we use another
type o f wavelets. This problem is not discussed. The derivatives of the Helds need to be
computed so complex integral need to be evaluated, it is not sure that the Battle-Lemarie
functions give the best computational result.
However, it has been demonstrated that these algorithms represents higher
order finite difference schemes and lead to improved linear dispersion and accuracy. It
would be very interesting to investigate MRTD scheme for other wavelets. In case we
would like to use more that one wavelet subspace, the integration might generate extra
computational overhead. That is the reason why the one may want to work in the physical
space instead of the wavelet space. The wavelet coefficients would represent the total
Held in a more direct way. The wavelet collocation method is an alternative to this
algorithm.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
104
4.2.2
The wavelet collocation technique
The main source of information for this method comes from the Uppsala
University in Sweden [103-104]. Walden worked on signal processing oriented algorithm
[10S] using the filter banks. He obtained a very attractive but extremely complicated
method. He also worked on Galerkin methods using wavelets [106], but has demonstrated
by Keiser [99] this is computationally inefficient for non-linear equations. Holmstrom
and Walden applied Galerkin based wavelet method to hyperbolic equations [107] before
Holmstrom used interpolating wavelets to solve PDE’s [103-104].
The interpolating wavelet method or collocation method is based on the sparse
point representation (SPR). Given a dyadic grid we define a “starting” coarse mesh and a
“finishing” fine mesh, i.e we can define as many levels as we want. Then, we go from
one mesh to the finer one by interpolating from the coarse mesh with a given polynomial
usually cubic polynomial. The wavelet coefficients represent the error between the
interpolated value and the exact value. This method allows compression of the unknowns.
Hence, matrix multiplication and differentiation can be done much faster. There is no
expansion using the wavelets it just represents an algorithm for fast differentiation and
multiplication in the physical space.
This scheme is based on the interpolating scheme developed by Donoho
[100]. The wavelet used is nothing else than the auto-correlation function of compactly
supported wavelets has demonstrated by Saito and Beyikin [108]. A thorough
mathematical study can be found in [109-110]. Cai and Wang investigated initial
boundary value problems for non-linear PDE’s [111]. Lippert et al. presented multilevel
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
105
computation using interpolating wavelets in [112]. Bertoluzza and Naldi [113]
demonstrated that Galerkin method using the auto-correlation function leads to the same
scheme and that wavelet preconditioning can reduce the condition number to 0 ( 1).
In the engineering field, Rubinacci [47] applied this scheme to l-D Maxwell’s
equations without absorbing boundary conditions and Toupikov [114] solved the 2-D PN
junction problem.
4.23
Mesh refinement and Daubechie wavelet expansion
Jameson does not work in the field of electromagnetics but has done some
extensive research on wavelets while working at NASA Langley Research Center
ICASE. He demonstrated in [115] that Daubechie wavelet expansion represents a
localized mesh refinement. His so called Wavelet Optimized Finite Difference method is
based on a non-uniform grid generation using Daubechie wavelet. As shown in MRTD
you can obtain higher order schemes and control the order of the method on different
localization of the grid [116]. His wavelet based grid generation algorithm and program
can be found in [117]. In [118], the latest version of this method, 2-D problems are
treated with a multi-domain technique which allow to use different grid and different
schemes on various domains. The computational problems of differentiation in the
wavelet based were explained in [119]. Jameson concludes by saying that wavelet are
best used in finite difference method if used just to generate a non-uniform grid, which
represents a similar idea than the interpolating wavelet scheme. This is specially the case
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
106
for multidimensionals problems where the matrix differentiation becomes really involved
and computational efficiency can be lost. Some applications can be found in [120],
4.2.4
Operators in wavelet bases and wavelet preconditioning
G. Beyikin seems to have produced the first engineering oriented work using
wavelets [108]. His idea was to discretize operators in a wavelet basis to obtain fast
computation. This can linked to the idea used in moment method, where sparse MoM
matrices are expected to reduce the CPU time. He demonstrated that wavelet
preconditioning can reduce the condition number especially for high order vanishing
moment wavelet [97-121]. Standard fast wavelet transform and non-standard FWT are
investigated. This idea is closer to post-processing technique. An operator is descritized
and the resulting matrix is transformed in the wavelet domain. The iterative solver used
afterwards deal with a sparse matrix with reduced condition number. Therefore, it
dramatically reduces the computation time. It is to be noted that FWT is an orthogonal
transformation that does not change the condition number this is really something that
should be discussed in the case of MoM where non preconditioning has been used so far.
His Ph-d student: Keiser worked on non-linear PDE’s [99] that we talked about
previously. Beylkin’s latest work can be found in [96], Jaffard also applied the same idea
in [122]. His main concern as a mathematician was to demonstrate the real advantage of
wavelets as a tool to solve PDE’s. His feelings can be followed through out his
publications. The advantages are not so clear except for linear problems.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
107
43
Interpolating wavelets applied to semiconductor simulations
MRTD has been extensively developed. In Global Modeling one need to solve
Maxwell’s equations consistently with the semiconductor transport equations. These
equations are non-linear and therefore cause a significant problem using standard
wavelet-Galerkin based technique. We use an interpolating wavelet scheme to solve the
non-linear partial differential equations that characterize the behavior of semiconductor
devices. We apply this method to a typical Field Effect Transistor. The I-V characteristics
are obtained and the accuracy is compared with the basic finite difference scheme. A
reduction of 83% reduction in the number of unknowns at steady state is obtained with 3
% error. The performance of the scheme is investigated using different values of
threshold and two types of interpolating wavelet namely the second order and the fourth
order wavelets. Due to the specific problem analyzed here, a tread off appears between
good compression, accuracy and order of the wavelet. This is a significant step toward a
unified numerical technique that uses wavelets to solve Maxwell’s equations and the
semiconductor equations for global modeling of high-frequency circuits.
43.1
Wavelet numerical techniques
The full wave analysis of microwave circuits (Global Modeling) is a
tremendous task that requires involved numerical techniques and algorithms. In [1] the
authors self-consistently solve Maxwell’s equations together with the semi-conductor
equations that characterize a sub-micrometer gate device. Alternative techniques have
shown interesting results. The extended FDTD method [16] that can deal with an
equivalent circuit of the active device has been used to predict nonlinear phenomena and
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
108
EMC effects [18]. A hybrid technique using frequency domain finite element method and
a circuit simulator has also demonstrated good capability to predict the behavior of
different transistor topology using the same process [IS]. Those techniques however
cannot simulate the propagation effects occurring inside the transistor at very high
frequencies. Therefore it is our belief that there is an urgent need to present a new
approach to the computational problem of global modeling to make this technique
practical for circuit design at millimeter-waves.
On one hand, one can decide to use the Finite-Difference Time-Domain
(FDTD) technique to solve for the electromagnetic fields of the passive and active parts,
but soon will face a difficult problem: The cell size of passive parts will be almost as big
as the whole transistor. Therefore computationally expensive techniques like time domain
diakoptics must be employed [52-57]. On the other hand, implementing a technique that
adaptively refines the mesh in domains where the unknown quantities vary rapidly,
would considerably reduce the number of unknowns. Such a technique corresponds to a
multiresolution analysis of the problem. A very attractive way of implementing a
multiresolution analysis is to use wavelets [44]. Wavelets have been used in
electromagnetics for a few years, first in the method of moments [33], and later in FDTD.
It was demonstrated [46] that finite difference schemes can be derived using the method
of moments with wavelet expansions of the fields. The resulting numerical technique has
been called the Multiresolution Time Domain Technique (MRTD)[92]. This method has
been studied extensively and shows very good performance as for the accuracy, memory
requirements and CPU time. Nevertheless, the implementation of wavelets to semi­
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
109
conductor equations is still to be investigated thoroughly in order to determine if a
multiresolution numerical method can be used in the context of global modeling. The
MRTD can be regarded as a wavelet-based Galerkin method. For non-linear equations
such as semiconductor modeling equations, this method can become quite time
consuming [99]. Therefore, a different wavelet approach is investigated, bearing in mind
the possible hybridization of different multiresolution techniques in the future.
We propose to apply a wavelet-interpolating scheme to the semiconductor
equations. The equations are solved on a dyadic grid. The wavelet coefficients are
directly related to the physical domain as they represent the error between the exact
solution on the grid and the interpolated value from the previous coarser mesh. This
scheme allows us to multiply and differentiate very quickly. We will follow the algorithm
explained in [103] that has been used to solve 1-D Maxwell’s equations in [47] and a 2-D
PNjunctionin[114].
4.3.2
Interpolating wavelets
It is based on the interpolating subdivision scheme studied by Deslauriers and
Dubuc [95]. It was later on extended in an interpolated wavelet transform by Donoho
[100]. The idea is to consider a set o f dyadic grids that defines different resolution levels.
A grid contains all points of the coarser level plus points inserted half way between those.
Values at odd grid points are kept unchanged while values at even points are interpolated
by a polynomial. The set of dyadic grids can be represented as shown in Fig. 4.5.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
110
Vo
•
•
•
•
•
V,
.......................................................................
V2
.......................................................................
Fig. 4.5 A set of dyadic grid
If we start on the coarser grid with the Kronecker delta function and
interpolate the values at even grid points, we build the scaling function. According to the
order of the interpolating polynomial the scaling function built is different. Fig. 4.6 shows
the scaling function for two different type of interpolation. Linear interpolation p -2 and
cubic interpolation p=4.
(a)
0.5
-0.5
(b)
0.5
-0.5
Fig. 4.6 Scaling fiinction for (a) p=2 and (b) p=4
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
I ll
This scaling function is defined on the real line so special attention should be
made to the boundaries. In the case of linear interpolation the scheme remains the same,
but for other polynomial ip -4 ) the standard stencil must be modified. Fig. 4.7 shows the
left boundary scaling functions for level 3 when p -4 .
*3.0
*3.1
1.2
0.8
0.5
0.6
0.4
0.2
-0.2.
0.2
0.4
0.6
0.8
-0.5
0.2
*3.2
0.4
0.6
0.8
0.6
0.8
*3.3
1.2
0.8
0.6
0.8
0.4
0.6
0.2
0.4
0.2
-
0.2
-0.4
0.2
0.4
0.6
0.8
-
0.2
0.2
0.4
Fig. 4.7 Left boundary scaling function for p=4
Interestingly enough, this scaling function <pis the autocorrelation function of
the Daubechies wavelet [108]. For a given order p , <p{x) = J$(y)p(y -x jd y where $ is
the scaling function associated to Daubechies wavelets of p/2 vanishing moments. It was
shown in [113] that the expansion of the solution of an elliptic PDE in terms of
interpolating scaling function leads to the same linear system than the one obtained by
using Daubechies wavelets as test and trial function in a Galerkin method. From
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
112
construction this scaling function has a compact support [-p+l.p-1] and verifies the
dilation equation specific to wavelets:
*=q-i
<P(x) = £
(4.22)
gk<p{2x-k)
* = -p + l
This scaling function is used to refine the mesh in other terms move from
subspace Vj to subspace Vj+, . We can introduce subspace Wj to define the difference
between V} and Vj*i . These new subspaces can be represented as shown in Fig. 4.8.
Vo
•
•
Wo
W,
•
•
•
•
•
•
•
........................................ ......
Fig. 4.8 Wavelet subspaces representation.
A possible choice of functions that span the subspaces Wj is the set of scaling
functions itself. It defines therefore a non-orthogonal multi-resolution. In this
representation, the wavelet coefficients are the error between the value at odd grid points
and the interpolated value at a coarser grid. With this wavelet representation, one can
refine or coarsen the mesh as desired.
Now we can create the so-called sparse point representation of a function /
[103]. Starting with the coarsest grid, we create an irregular mesh by computing the error
between the exact value of / o n the coarse grid and the value obtained by interpolation.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
113
The wavelet coefficients cany the detailed signal: the error between the exact value and
the interpolated value. By removing the points that can be interpolated, we greatly
compress the data defining a sparse point representation (SPR) that can be used for
computation. The number of points remaining in the SPR depends on the smoothness of
the function. If the function varies smoothly it will be easily described by a polynomial
interpolation and thus the SPR will contain very few points. If the function varies rapidly
then the SPR will contain more points, as the interpolation will be insufficient to coarsen
the mesh. The compression that one can obtain depends on the threshold value set as the
minimum error for the interpolation procedure.
This extremely versatile technique is expected to be very useful to solve
partial differential equations (PDE), especially non-linear ones. Multiplication and
derivatives can be computed in this SPR, which may considerably reduce the
computation time. The typical way of solving a PDE using the SPR representation is to
transform the PDE in a system of ordinary differential equation (ODE) using the method
of lines. The initial conditions are set and the first SPR is obtained, the spatial derivatives
are approximated using the interpolated wavelet scheme [104] and the new time iteration
is performed using an ODE solver. The type of PDE that are to be solved may introduce
some changes in the algorithm, as the SPR may not need to be computed at each time
step. The ODE solver also offers a great degree of freedom as one can choose between
standard Runge-Kutta methods, multistep methods and so on. In general, the more
expensive the ODE solver the more efficient this scheme will be as the overhead created
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
114
by the SPR manipulation will be overcome by the spatial derivatives approximations. In
summary the algorithm to solve a differential equation using this scheme is as follows:
•
•
•
Set initial conditions
Obtain the SPR of the unknowns
Approximate spatial derivative using SPR interpolation
•
Advance in time using ODE solver
•
Go back to the SPR computation
A similar scheme has been employed to solve the problem of a 2-D PN
junction [114]. However, in this previous work the electron and hole mobilities were
assumed constant. This assumption severely limits the validity of the technique for
simulating modem semiconductor devices used for high speed or high-frequency
applications. Therefore, in this research, we propose to treat the problem of a Held effect
transistor with mobility and diffusion coefficient as functions of the electric field.
4.3.3
Drift diffusion example
Poisson’s equations can be discretized by finite difference schemes. This
results in a linear system that can be solved by any iterative solver provided that the
matrix representation of the discretized operator satisfies the stability and convergence
criteria such as: positive definite matrix. As the Poisson’s equations is to be replaced in a
future work by the Maxwell’s equations, no significant efforts where put to gain
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
115
computational time while solving this equation. However section 4.5 will deal with this
particular problem using Beylkin technique. All the effort was concentrated on the
continuity equation that exhibit non-linearity. However as mentioned in the previous
section, a discretization of Poisson’s equations using Daubechies wavelets leads to the
same system as the one obtained by interpolating scaling functions [113]. Thus a very
effective technique called diagonal pre-conditioning could be used to dramatically reduce
the condition number of the discretized operator matrix and therefore reduce the number
of iterations needed for convergence [97]. The continuity equation is solved using the
interpolating scheme presented in the previous section. Time domain simulations are
performed. One can verify that the gate current is zero as effectively no current must flow
through the schottky contact. A relative error is defined in order to compare the standard
finite difference method that uses a uniform mesh and the interpolating wavelet scheme
that generates the non-uniform mesh. This relative error is defined on the drain current
as:
(4.23)
relative
Where Je is the drain current density computed by the interpolating wavelet scheme using
a threshold of e on the SPR of the carrier density. J fd is the drain current density
computed by a standard finite difference code.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
116
4.3.4
Results
Initial results were introduced in [71]. A MESFET with the following
dimensions is simulated: 0.6 fan gate length, 1 fan long source and drain electrodes, 0.7
fan source-gate gap, 2.5 fan gate-drain separation, 0.2 fan deep channel layer and a 0.8
fan deep buffer layer. The doping of the active layer is 2.2xl017 A/cm3 and the one of the
buffer layer is lxlO 14 cm'3. The zero-field mobility is 800 cm2/V .s , the saturation velocity
Vo, is l ( f cm / s . A 65*65 mesh was used for spatial discretization while the Euler’s
method was used as the ODE solver. The time step was chosen as 8.10'15 s , but due to
the non-uniform meshes used some tuning was required to ensure stability. Fig. 4.9
shows the I-V characteristics of the simulated MESFET.
150
vgrt.
vg*.-1v
too
I
I
0.5
1.5
Z5
3.5
45
V<*
Fig. 4.9 IV-Characteristics of the simulated MESFET
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
117
A typical behavior is observed. The electrostatic potential contour plot is
shown in Fig. 4.10, for a drain voltage of 2 Volts and a gate voltage of -1 Volt.
0.9
0.8
0.7
0.6
0.S
0.4
0.3
0.2
0.1
0.5
1.5
2.5
3.5
4.5
Fig. 4.10Contour plot of the Electrostatic potential for a drain voltage of 2 Volts and
an applied external gate voltage of -1 Volt,
In Fig. 4.11, the carrier distribution is plotted and the corresponding SPR can
be compared in Fig. 4.12 for the same bias condition as the potential shown in Fig 4.10.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
118
*10'
0.9
0.8
JL£
0.7
0.6
0.5
0.4
0.3
0.2
0.5
1.5
2.5
3.5
4.5
Fig. 4.11Contour plot of the carrier density (the labels are normalized by 1017) for 2
Volts drain voltage and -I Volt applied gate voltage
We observe that the depletion region was created in the channel. The SPR
kept mesh points in the grid where the carrier density varies rapidly, namely, around the
depletion region and along the junction between the channel and the buffer layer. In
smooth regions, where the interpolation is more accurate, fewer points are needed. Thus,
few points are necessary to predict the carrier density deep down in the buffer layer. This
representation changes in time, which makes the grid dynamically adaptive. For different
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
119
bias points, the SPR is obviously different as the depletion region can be reduced or
enhanced.
• i •' HMN
• • • mSmm i <
3 S T
^
10
*
i
20
• •
301*
40;
5 0'
60
10
20
30
ru *703
40
50
60
Fig. 4.12Sparse Point representation of the carrier density at 2 Volts drain voltage
and -I Volt applied gate voltage.
The performance of the new scheme is investigated at a single bias point: 2
Volts drain voltage and -1 Volt applied gate voltage. Fig. 4.13 shows the drain, source
and gate currents computed at each time step.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
120
200
100
\S o u r c e Current
-100
a -300
-400
-500
-600
200
400
600
800
1000
1200
1400
1800
Time Iterations (dt = 8.1(T1S)
Fig. 4.13Time Domain drain, source and gate currents at 2 Volts drain voltage and / Volt applied gate voltage.
A conventional behavior is obtained as the source and drain currents converge
to the same value and the gate current tends to zero. Three absolute wavelet thresholds
were chosen: e =0.01, e =0.001 and e =0.0002. The relative threshold was defined as
e.Nd, where Nd is the doping of the channel layer. A set of simulations was performed for
these three threshold values and for two cases of scaling functions p=2 and p=4 . Fig.
4.14 represents the number of mesh points remaining in the SPR after thresholding of the
wavelet coefficients.(i.e after eliminating the mesh points that can be interpolated from
the remaining points with an error e ).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
121
1200
1000
e « 0.0001
<2 800
Z 600
e « 0.001
Z 400
200
e > 0.01
200
400
600
800
1000
1200
1400
1600
Time Iterations (dt * 8 .10~1S)
Fig. 4.14Mesh Adaptability for three different values of wavelet threshold. ( Solid
lines p=2 , dashed lines p=4 )
Different conclusions can be drawn from this numerical experiment. As the
threshold gets smaller the error that is tolerated in the representation of the carrier density
gets smaller, and more points are necessary to accurately represent the carrier density. At
the beginning of the simulation the carrier density is initialized to the doping profile thus
very few points are needed due to the smoothness of this initial condition. As time
evolves, the depletion region starts to appear and points are added in the SPR, this why
for the first 400 iterations the number of unknowns grows. When the depletion region is
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
122
created and the MESFET starts to reach its steady state, the carriers concentration stops
changing and the number of points in the SPR remains the same.
The remarks made above, illustrate the dynamic behavior of the mesh. More
involved is the behavior of the SPR for p=2 and p=4. Fig. 4.14 demonstrates that for e
=0.01 and e =0.001, the linear interpolation achieves a better compression ratio than the
cubic interpolation. For £ =0.0001 however, it is the opposite. When the carrier density is
almost constant (i.e in the buffer layer and in the channel far from the depletion region),
the linear interpolation performs better. Only two points can be used to describe the
correct solution: the coarser mesh of the linear interpolation requires fewer points than
the cubic interpolation. In the depletion region where the carrier density varies rapidly, a
significant error is acceptable. So, the linear interpolation still uses less unknowns than
the cubic one. At e =0.0001, it appears that the breaking point has been reached, and the
linear interpolation can no longer describe the depletion region with the required
accuracy without using the finest mesh. It is our understanding that this breaking point is
a function of the wavelet threshold, the bias point and the physical dimensions. The
buffer layer can be made smaller, which would result in a depletion region occupying
more space in the computational domain. The cubic interpolation could therefore
outperform the linear one for a much bigger error than it is the case in this experiment.
Finally, Fig. 4. IS presents the relative error defined in the previous section.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
123
E S 0 .0 0 1
e = 0.0001
8 10'
200
400
600
800
1000
1200
1400
1600
Time Iterations (dt = 8 .1<T1S)
Fig. 4.15Relative Error on the drain current for three different values of wavelet
threshold. ( Solid lines p = 2, dashed lines p=4)
This graph compares the solution o f a standard finite difference scheme, using
uniform mesh, with the interpolating scheme. Once again the results are plotted for three
wavelet thresholds and the two previous interpolating scaling functions. For the two cases
e =0.01 and e -0.001, the error on the drain current is respectively 50 % and S %, and
there is no major difference between the two scaling functions. This shows that for these
given wavelet thresholds, one would choose the linear interpolation to achieve better
CPU time, hi the case of e=0.0001, the linear interpolation gives a relative error of 0.8 %
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
124
while the cubic interpolation gives an error of 0.1 %. These results help to determine how
to choose the wavelet threshold.
In this study, we also examined the computational overhead created when the
SPR are generated. On one hand, this scheme shows good promises with a potential
reduction in the computation time between 74 to 95 %. On the other hand, the overhead
currently reduces the gain in computation time to a range of 25 to 55 % for the different
thresholds presented earlier. An appropriate data structure is expected to improve those
performances. It is also to be noted that more computationally expensive ODE solver
could be used which would reduce the share of CPU time used to deal with the SPR’s.
The initial size of the mesh that defines the finest grid also affects the CPU time
performances of the scheme because the size of the system of ODE changes.
4.4
Interpolating wavelets applied to a 2-D cavity
The interpolating scheme suits the semiconductor transport equations.
Therefore, it may as well be suitable for Maxwell’s equations. This would allow us to
create a unified technique for global modeling simulations without the need to use an
MRTD scheme. We propose to generate a non-uniform mesh using an interpolating
wavelet scheme. We will follow the same algorithm presented in the previous section.
We will use this scheme to solve a two-dimensional cavity (TE case) enclosing a cylinder
representing a discontinuity inside the cavity. This cylinder will allow us to see the
refinement process that occurs during the simulation.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
125
By removing the points that can be interpolated, we greatly compress the data.
The compression depends on the threshold value set as the minimum error for the
interpolation procedure. AH the computations are then done on the irregular mesh.
4.4.1
Sparse point representation
A 2-D cavity was simulated according to the example presented in [27]. The
cavity is discretized by a mesh of 65*33 with a space increment dx-3.0 mm. The space
increment is dt=5.0 ps. Fig. 4.16 shows the y component of the electric field at two
different time.
10
20
30
40
so
60
40
50
60
(a)
10
20
30
(b)
Fig. 4.16EIectric field in the y direction, (a) t=520 ps, (b) t=840 ps.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
126
We can see that at t=520 ps the wave did not reach the cylinder yet, whereas
at t=840 ps, the modes inside the cavity are developing and scattering due to cylinder is
occurring.
Fig. 4.17 shows the sparse point representation or in other terms, the nonuniform mesh generated by the scheme at the same times than Fig. 4.16. This
demonstrates the self-adaptibility of the mesh which gets finer in regions where the fields
are varying rapidly. At t=840 ps, scattering due to the cylinder needs to be modeled
accurately therefore more points are used around the cylinder.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
127
0
V •
•
9 9 9 9 9 9 9 • •
•
•
9 9 9 • •
• ••• •
99999 • 9 99999 •
•
9999
• • 999999
•
9999999999999999 •
99999
••
•
•••••••
5
- ••••••
19
999999
•••
9999999
10
•
• ••••• • ••• •
•••
999
• • • • • • • 9999999
999
•••
9999999 9999999
* 999999 99999 •
••
15
20
25
30
•
•
•
•
••••• • •
•
9999999 •
9999
9999999 •
999999999999999999999999999
••
1999
••
••••
999999
9999
99999 999999999999999999999 99999999999 • •
••
••
••••
999999
999999999999999999999999999
99999999
••••
99999999999 999999
9999999999999 • • 999999999999999 99999999999 • •
99999
••
••
9999999
• •
999999 • • • • • •
9999999999999999 •
••••
> • • • • • • • 9999999 • • • • • • • 9999999 • 99999 9999999 •
• 99999
• • • 999
• • • 9999
• • • 999
• • • • • • 9999999 9999999 9999999 • •
••• •
99 • • •
••• •
•
• • • • • • • • • • • • • • • • • 99999 • • • • • • • • •
• •
• • • • •A
w
•
•••••••
• • • 9 ,9 • •
20
30
40
•
•
•
•
•
•
•
■
•
•
•
•
•
I* .
10
•
•
•
.
- 9999
99
99 • •
••
••••
0
•
•
•
•
•>
•
•
•
•
•
•
•
50
60
(a)
0
• • • • • • • • * •• •• •• • •* ••••• •••* ••• ■ ••••* ••••• • ••* •••••••« •* •••••
•••••••
••••••■ • • • •
•••••• • ••••••• ••••••• ••••••• ••••• ••••••••••••••••••••• ••
5
10
• • • • : : : •••••••
• s : : : .................... . • : : : : : :
• • ■ ••••• ••••••• •••••••••••••••••••••••••••••
• •
••••••••• •••
■
•
• •
•••••••••••• •••••••••• • • • • • • • • • _
• • ••• • ••••••• ••••••••••••••••••••••••••• •
•
• ••••••••••••••••••••••• ••••••••••••••••• • •
• • ••••••••
15 - • ••••••••••••••••• ••••• ••••• •••••••••••
• • •833333.
.33383... . . . . . 88S88
3 8 ......... • • •
•••• •••••••
• . .• .• .• .• a• . . .•.•.•. •. •. . . . ••••• ••••••
..........
20
• .3 3 S 3 3 3 3 3 .S S 3 S 3 3 3 .........3 3 3 . . . . . . . . . . . a a
a
• ...... •
•••••
•
_••••
a a ....
8 8....................................... 888888 88888888888.* * * * 888888**
• "8888
••••••••••••..
.
• • • • • • • 9 IM IH • • • • • • • • • • • • • • • # • • • • • • • • • • • • • • • • • • • • • • •
30 •••••• • ••••••• ••••••• ••••••• ••••• ••••••••••••••••••••• •
•
25 * * 8 888***•••••• ••••*••888888 8888*88888888.8 • 8 8888888 * *
M M 999
0
I f M I I M f M IIM f f f M
10
20
30
40
f M
50
M M If • •
60
(b)
Fig. 4.17Sparse point representation of the y component of the electric field at (a)
t=520 ps and (b) t=840 ps
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
128
4.4.2
Compression characteristics
At every time step this non-uniform mesh is generated by the wavelet scheme.
Thus the number of unknowns varies in time. Fig. 4.18 shows this behavior. Three
wavelet thresholds are used. We can see that as the threshold gets smaller, more points
need to be used in the mesh. The error computing during the interpolation needs to be
smaller thus finer mesh is used.
2000
1800
1600
ep=0.05
ep=0.01
ep=0.001
5 2 1400
C 1200
1000
800
600
400
200
Time (s)
x10*'°
Fig. 4.18Number of unknowns remaining in the SPR for three different value of
wavelet threshold.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
129
4.5
Wavelet transform for efficient Poisson’s solver
In this section we present a technique that uses wavelet transform to solve a
linear system. This technique is based on the work by G. Beylkin and was introduced in
section 4.2.4. We apply this method to solve the potential of typical FET structure.
The discretized Poisson’s operator that results in a band diagonal matrix is illconditioned. Wavelet transform and diagonal-preconditioning are used to reduce the
condition number and therefore speed up the convergence.
The following section will briefly explain how to form the matrix representing
the Poisson’s operator and how to transform this matrix to make it symmetric. Symmetry
of the matrix is a necessary condition to obtain convergence using the conjugate gradient
method. Section 4.5.2 will present the main iterative solvers that can be used to solve the
system. In section 4.5.3, the wavelet transform will be introduced with the diagonal
preconditioning, while in the last section the convergence results will be presented.
45.1
Poisson’s operator for FET structure
Poisson’s operator is mainly a second order derivative. In this discussion we
will not worry about the right hand side term. We will mainly focus and the impact of the
boundary conditions regarding the symmetry of the matrix.
In Fig.4.I9, a typical FET structure is presented. Dirichlet boundary
conditions are used at the electrodes while Neumann boundary conditions are used at the
other walls.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
130
Gate
Source
Drain
Active layer
Buffer layer
Fig. 4.19 Conventional FET structure
The finite difference method is used to approximate the second order
derivative. We use an up and down ordering of the unknown. We set the following factor:
I
a-2.
1_
(4.24)
(^Ax2 + Ay 2]
The matrix resulting from this discretization is shown below. It can be noted
that in this case, only Neumann boundary conditions are used. Dirichlet boundary
conditions can be added easily, by deleting the rows and columns corresponding to the
grid points where the potential is enforced.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
131
2
a
‘ Ay2
1
Ax'
___1_
a
Ay2
Ay2
2
A2
0
2
0
Ax2
0
2
a
Ax2
1
a
Ax2
1
1
Ax2
Ay2
1
2
1
Ay2
Ax2
u/y
2
Ax2
Ay2
__ 1_
__ l_
Ay2
Ax2
__ 1_
a
Ax2
2
Ax2
0
0
0
2
a
Ay2
2
1
Ax2
Ay2
I
a
Ay2
2
2
Ax2
Ay2
a
We can see a 3*3 structure, corresponding to a 3*3 grid resulting in 9
unknowns. This matrix is non-symmetric and therefore the conjugate gradient algorithm
would not converge. A few modifications must be performed in order to make this matrix
symmetric.
Rows four and six need to be multiplied by one half. The same would happen
to any other middle block in the case of a large grid. The first and the last row of any
middle block would need to be multiplied by one half. The first and last block, which
discretize the left and right boundaries need the same transformation. The first and last
rows of these two blocks need to be multiplied by one fourth and the rest o f the rows
need to be multiplied by one half.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
132
These transformations make the discretized operator a symmetric matrix. The
size of this matrix is then reduced, by deleting the grid points corresponding to Dirichlet
boundary conditions. The right hand side term needs to be modified accordingly (only for
the Dirichlet points).
4.5.2
Iterative solvers
Once the problem has been discretized, an iterative solver needs to be used to
solve the remaining linear system. The basic preconditioning techniques are based on
splitting of the matrix A, it results in iterative schemes like the symmetric successive over
relaxation (SSOR). Another type of preconditioning, which can also be seen as an
acceleration technique is the preconditioning by projection. An example of such iterative
solvers is the conjugate gradient (CG) algorithm.
In the case of splitting matrices we have:
A =M -N
(4.25)
Ax=b<=> M x= b+ N x
(4.26)
So the system to solve becomes:
Which leads to the following iteration:
X (k) = X(M +M~l ( b - A X {k~l))
(4>27)
The basic algorithms are based on the splitting of A with respect to the diagonal matrix
D, the upper triangular matrix U and the lower triangular matrix L. Such a splitting is
shown in equation 4.28.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
133
A=L+D +U
(4.28)
If one considers the following M matrix:
M=D
(4.29)
M{=D+L
(4.30)
We obtain the standard Jacobi iteration.
If one considers the following Mi matrix:
We obtain the forward Gauss Seidel iteration. While if we consider M2 such as in
equation 4.31 we have the backward Gauss Seidel. By combining the two iterations, one
can implement the symmetric Gauss Seidel.
M Z=D+U
(4.31)
A more efficient iteration is obtained by using a relaxation factor in combination with
these splittings. We can set:
(4.32)
1
M, = - D + L
0)
Af, =— D+(J
'
(o
With Mi, we get theforward successive over relaxation (SOR) and with M2 we get the
backward SOR.By combining the two we get the well-knownSymmetric
Successive
Over Relaxation (SSOR). The final M matrix to be used in the iteration 4.27 is:
u =
^ h $ D + a L )D 'l ( D * M )
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
l4 J 3 )
134
These
couple
algorithms
represent
iterative
solvers
obtained
by
preconditioning with splitting. We can also obtain iterative solver by using projection
technique. The solution of the linear system is looked for within a subspace represented
by combining the columns of the matrix A.
X=VX
(4.34)
V represents the trial space. By combining the rows of A, we project the equations of the
linear system onto a subspace defined by the test space W. Equation 4.3S describes such
a projection.
W'AX=W'b
(4.35)
If we combine both ideas of the projection, we obtain an approximation of the solution in
equation 4.36.
X = V (W 'A vfW b
(4‘36)
The subspaces W and V can be chosen to be the same, in which case we have a Galerkin
type of algorithm. The subspaces can also be defined at each iteration in order to
converge to the solution faster. As we get closer to the solution, the subspaces where we
project the solution can be refined. Such a method is called steepest descent. By choosing
W and V differently we get different steepest descent schemes. The most well-known
algorithm being the conjugate gradient (CG).
Now that the main iterative solvers have been introduced, remains the
question of which one is best suited for one’s problem. Clearly, the accuracy needs to be
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
135
obtained, but a researcher needs also to worry about the efficiency. When dealing with
quasi-static device simulations, Poisson’s equation needs to be solved repetitively.
Therefore, an efficient solver has to be implemented. The CG algorithm is quite efficient
and can be improved by using SSOR preconditioning for the first iterations resulting in
what is known as the preconditioned conjugate gradient (PCG) algorithm. The
convergence speed can be also improved by using a sparser matrix. More zeros form the
matrix, so less operation need to be done. Finally, the condition number can be improved.
This number is the ratio of the largest and the smallest eigen-values of the matrix. The
closest to unity the condition number is, the less number of iterations are required to
converge to the solution.
4.5.3
Wavelet transform and diagonal preconditioning
The philosophy behind using wavelets in Moment’s method is to get a sparse
matrix instead of a dense matrix [33]. Researchers realized the difficulty of deriving the
MoM matrix using wavelet expansion. Wavelet transform can be directly applied to the
dense matrix and the linear system of equations can be solved in the wavelet domain. The
solution can be reconstructed afterwards by inverse wavelet transform.
However, wavelet transform is an orthogonal transformation so the condition
number of the matrix remains unchanged. Sparsity is obtained by thresholding, but it is
unclear what is the effect of such a process on the condition number o f the matrix and
what is the mathematical background behind it.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
136
Improvement in the condition number is only obtained by diagonal
preconditioning of the matrix in the wavelet domain. Under such treatment, the matrix is
better-conditioned and any sparse solver such as CG will converge faster than the original
system.
In section 4.S.1, we derived the matrix A, that represents the discretized
Poisson’s operator for a standard FET structure. We explained how to modify this matrix
to make it symmetric. Therefore, CG can be applied to solve the modified system Am. We
start with the standard linear system of equation 4.37:
Ax = b
(4.37)
After modification, A is symmetric and the Dirichlet boundary conditions are removed.
The resulting system is as follows:
(4.38)
That system is of reduced size compared to the original system. Wavelet transform is
applied to this system.
A
V ux i =b m
(4-39)
A diagonal preconditioning dependent on the number o f wavelet decomposition level is
used.
(4.40)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
137
One can link the wavelet transform to the types of iterative solver presented in
the previous section. We see that wavelet transform is nothing else than a projection of
the equations of the system on another subspace, and that the diagonal preconditioning is
a projection of the solution on another subspace.
Convergence results and condition number improvement will be presented in
the following section. We will know present the sparsity of the matrix and the effect of
the different wavelets that can be used.
Several level of decomposition can be used, and different wavelet function
can used. In this research we focused on the Daubechies wavelets, but it must be noted
that other wavelet could produce better or worse results. A wavelet can be characterized
by a filter of a certain length [44]. We used the first order Daubechie wavelet (Haar), the
eighth and the twentieth to show the impact on the sparsity of the matrix. In Fig. 4.20, we
present the sparsity obtained when using Haar wavelets.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
138
100
150
200
250
300
350
400
450
500
0
50
100
150
200
250 300
nz * 7452
350
400
450
500
Fig. 4.20 Sparsity of A using Haar wavelets.
The filter is short so the sparsity is relatively good. The number of non-zeros
elements is 7432. We see in Fig. 4.21 and Fig. 4.22 that the sparsity is degrading because
the filters become bigger. The number of non-zeros elements goes from 47646 to 96798.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
100
150
200
250
300
350
400
450
Fig. 4.21Sparsity of A using db8 wavelets.
n c= M 7 »
Fig. 4.22 Sparsity of A using db20 wavelets.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
140
The number of decomposition levels also affects the structure of the matrix.
Fig. 4.23 shows the typical finger like structure for only one resolution level.
o
50
100
150
200
250
300
350
400
450
500
0
50
100
150
200
250 300
nz*22944
350
400
450
500
Fig. 4.23Sparsity of A for one resolution level
This structure is more understood for two and more level of decomposition.
At this level the sparsity looks like the fingers of the hand. In the next section we present
the results regarding the condition number. Comparison will be made depending on the
number of decomposition level and the type of wavelet used.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
141
a
90
100
190
200
290
200
350
400
490
900
0
90
100
190
200290300390400490900
nz=30906
Fig. 4.24Sparsity of A for two resolution levels
100
200
2901
300
390
400
490
500
Fig. 4.25 Sparsity of A for three resolution levels
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
142
4.5.4
Results
The grid discretizing the FET was set so that the modified matrix Am would
be of size 2n where n is a positive integer. The condition number of the matrix A is 1183.
In the following table, we present the variation of the condition number according two
the wavelet used in the FWT and two the number of resolution levels.
Table 4.1 Condition number of A versus wavelet type and resolution
levels
Dbl
Db4
Db20
One level
641
669
602
Two levels
408
307
265
Three levels
355
301
270
It is clearly shown that the condition number decreases as the wavelet is of
higher order and as the number of decomposition levels increases.
There is a tread-off to make between good conditioning and good sparsity.
From the previous figures it was shown the sparsity is destroyed, as the filters get longer.
So one may want to use a reasonable order of wavelets. Regarding the number of levels
used, one may consider the time that the FWT takes to be computed. To keep on gaining
CPU time, we may not want to take too long to compute the FWT even if it is done only
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
143
once. The use of fourth order Daubechie wavelet with two levels seems to be a good
compromise.
In that case, we solved the system using SSOR, CG and PCG iterative
algorithm. Fig. 4.26 shows the convergence behavior of the original system when A is
not symmetric. It can be seen that CG does not converge. When A is modified to Am,
then CG converges.
• r: : u •; • • i n i n : i : : : n •
■-
n i.m; i
sso r
m m ijn iiiiim r iH il
i.i
•
100
r
!
200
:*L - i ! r l l ! i - i l * : i ! I : i :
300
immnn.nmnnm.:;’
.P ? . ■
T r . £ * *.;H i
400
=; • • inw.ii;• • • u:i: \ •
500
600
Fig. 4.26Convergence behavior of SSOR, CG and PCG on A and Am
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
144
When performing the FWT on Am, we showed that the condition number was
improved. It is clearly shown in Fig. 4.27 where every iterative solver tested, converges
faster than on the modified system.
V\
— CGw
— pCGw
— CGm
1
•
100
200
300
;;j r n
400
r;n n ; i
500
;
600
iteration number
Fig. 4.27 Convergence behavior of SSOR, CG and PCG on Araand \
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
145
4.6
Summary
The basic theory of wavelets numerical techniques has been introduced and a
through review of the research done in the last decade has been presented. This explains
why we focused on using interpolating wavelets to solve the semiconductor transport
equations. MRTD schemes are two computationally expensive to deal with those
equations. It was also discussed that the efficiency of such schemes can be argued. The
complexity of the implementation of absorbing boundary conditions is a real hurdle and
the integration needed to compute the wavelet coefficients can be computationally costly.
Interpolating wavelets are a very good tool to generate non-uniform mesh that
can be adaptive in time, but stability remains an issue. Also multilevels of resolution
create an computational overhead. Semicondutor and electromagnetic simulations were
performed using the interpolating wavelet scheme.
Finally, FWT was used to reduce the condition number of the discretized
Poisson’s operator matrix. It was shown that there is a tread off between the number of
resolution, the order of the wavelet used, the computation time to perform the FWT and
the condition number. We demonstrated that solution of the Poisson’s equation can be
obtained much faster using a wavelet representation of the operator.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 5
CONCLUSIONS
Global modeling simulation is the answer to a correct understanding and
modeling o f particle-wave interactions inside an active device. In order to perform such a
simulation, the scientist needs to be familiar with semiconductor transport,
electromagnetic theory and numerical techniques.
Semiconductor transport principles were explained in chapter one with an
emphasis on Boltzmann’s transport equation. This equation needs to be solved accurately
to simulate sub-micrometer gate length transistor. The derivation o f the balance equations
used to solve the BTE was shown. Then, it was explained how to solve them and some
simulation results were presented. These equations characterize the transport o f electrons
through the device and model the different semiconductor effects that can be
incorporated. Surface charge and deep electron traps can be added to the models. Several
approximations o f the full hydrodynamic model were presented, in order to study the
range o f validity o f every one.
In chapter three, Maxwell’s equations were solved using a time domain
numerical technique known as FDTD. The active device was incorporated within the
mesh using a modified FDTD method, the lumped-FDTD. The values o f the lumped
elements that characterized the device were obtained from the semiconductor simulations
performed in the second chapter. These simulations results were modeled by an artificial
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
147
neural-network that allowed efficient large signal simulations o f an amplifier. The basic
principles of the FDTD method have been introduced and some useful references have
been used for the reader who wishes to know more. A very comprehensive matrix
formulation o f ANN was introduced that will allow the reader to implement easily his
own model. Finally, global modeling simulations were performed using the new
developed technique. Large-signal simulations were performed in a very efficient
manner.
State o f the art numerical techniques using wavelets are presented in chapter
four. It begins with a brief overview o f multiresolution and wavelet theory.
Multiresolution analysis appears to be a good candidate for global modeling simulations
because it implements the large number o f scales that are present in a manufactured
transistor. Wavelet numerical techniques can be used in different ways. We focused on
the interpolating wavelet scheme that relates the wavelet coefficients to the real values o f
the unknowns. It represents a more efficient way (than MRTD schemes) to deal with the
highly non-linear equations that are the balance equations o f the BTE. The interpolating
wavelet scheme was used efficiently to create a non-uniform mesh to solve the
hydrodynamic model. It was also used to solve Maxwell’s equations. Finally, the fast
wavelet transform allowed us to speed up the solution o f Poisson’s equations.
The Neuro-FDTD method has been warmly accepted within the FDTD and
neural networks communities and has been referenced in conferences and books. The
numerical techniques developed using wavelets represent a one o f a kind effort to use
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
148
wavelets to solve the semiconductor equations within the framework o f global modeling.
It has also been recognized within the electromagnetic community.
Future leads in research, toward a more efficient, accurate and practical global
modeling technique include: Using a standard technique such as finite difference to
model a manufactured transistor and compare the model with measurements. This would
make a very strong statement about the validity o f the different models currently used for
high-frequency and high-speed applications. Another lead would be to use the NeuroFDTD method to its fullest potential using more simulations results for training and using
a better lumped-FDTD technique. Wavelets can be used efficiently to solve Maxwell’s
equations, but it would be recommended to develop an hybrid solver using finite
difference to deal with the semiconductor transport equations.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
REFERENCES
[1]S.M.S. Imtiaz and S. M. El-Ghazaly, “Global modeling o f millimeter-wave circuits:
Electromagnetic simulation o f amplifiers,” IEEE Trans. Microwave Theory tech ., vol.45,
pp.2208-2216, December 1997.
[2]M. A. Alsunaidi, S. M. S. Intiaz and S. M. El-Ghazaly, “Electromagnetic wave effects
on microwave transistors using a full-wave time-domain model,” IEEE Trans.
Microwave Theory tech ., vol.44, pp.799-808, June 1996.
[3] I. Erdin, R. Khazaka and M. S. Nakhla, “Simulation o f High speed Interconnects in a
Multilayered Medium in the Presence o f Incident Field,” IEEE Transactions on
Microwave Theory and Techniques, Vol. 46 No. 12, December 1998.
[4] K. A. Remley et al., “Characterization o f Near- and Far-field Radiation from Ultrafast
Electronic Systems,” IEEE Transactions on Microwave Theory and Techniques, Vol. 46
No. 12, December 1998.
[5] P. Ciampolini, L. Roselli and G. Stopponi, “Integrated FDTD and Solid-State Device
Simulation,” IEEE Microwave and Guided Wave Letters, Vol. 6 No. 11, November 1996.
[6] Kryloff, N. and Bogdliuboff, N.,"Introduction to Nonlinear Mechanics”, Princeton
University Press, Princeton, NJ, 1968.
[7] Kundert, K.S. et al., Applying Harmonic Balance to Almost-Periodic Circuits”, IEEE
Trans, on Microwave Theory & Tech., 1988, vol. 36, no. 2, pp. 366-377.
[8] K.S. Yee, “Numerical solution o f initial boundary value problems involving
Maxwell's equations in isotropic media,” IEEE Trans. Antennas Prop., vol.AP-14,
pp.302-307, May 1966.
[9] J. F. Mao, “Twofold Mur’s first order ABC in FDTD method,” IEEE Transactions on
Microwave Theory and Tech., vol. 46, pp. 299-301, March 1998.
[10] J.P Berenger, “A perfectly matched layer for the absoprtion o f electromagnetic
waves,” Journal o f Computational Physics Vol 114, pp. 185-200, October 1994.
[11] A. Taflove, “Adavances in computational electrodynamics: the finite difference time
domain method,” ARTECHHOUSE, 1998.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
150
[12] K.L Shlager and J.B. Schneider, “A Selective Survey o f The Finite-Difference TimeDomain Literature,” IEEE Antennas and Propagation Magazine, vol.37, no. 4, pp.39-56,
1995.
[13] C. Geuzaine et al., “An FETD approach for the modeling o f antennas” IEEE
Transactions on Magnetics, vol. 36, pp.892-896, July 2000.
[14] T. Itoh, “Numerical techniques for microwave and millimeter-wave passive
structures,” WILEY, 1989
[15] E. Larique et al., “Linear and non-Linear FET modeling applying an electromagnetic
and electrical hybrid Software” IEEE Transactions on Microwave Theory and Tech., vol.
47, pp.915-918, June 1999.
[16]W. Sui, D. A. Christensen and C. H. Dumey, “Extending the two-dimensional FDTD
method to hybrid electromagnetic systems with active and passive lumped elements,”
IEEE Trans. Microwave Theory tech., vol.40, pp.724-730, April 1992.
[17]M. Piket-May, A. Taflove and J. Baron, “FD-TD modeling o f digital signal
propagation in 3-D circuits with passive and active loads,” IEEE Trans. Microwave
Theory tech., vol.42, pp.1514-1523, August 1994.
[18]C.N. Kuo, B. Housmand and T. Itoh, ‘Tull wave analysis o f packaged microwave
circuits with active and nonlinear devices: An FDTD approach,” IEEE Trans. Microwave
Theory tech ., vol.45, pp.819-826, May 1997.
[19]C. N Kuo, B. J Housmand and T.Itoh, “Modeling o f Microwave Active Devices
Using the FDTD Analysis Based on the Voltage-Source Approach,” IEEE Microwave
and Guided Wave Letters, pp.199-201, May 1996.
[20]V. A. Thomas et al., “FDTD Analysis o f an Active Antenna,” IEEE microwave and
guided wave letters, Vol. 4, NO 9, September 1994.
[21] M. Chen, S. T. Chew and T. Itoh, “Nonlinear Analysis O f A Microwave Diode
Mixer Using The Extended FDTD,” IEEE Int. Microwave Symp. pp. 67-70 June 1997.
[22]V. A. Thomas et al., “The use o f SPICE Lumped Circuit as Sub-grid Models for
FDTD Analysis,” IEEE microwave and guided wave letters, VOL. 4, NO 5, May 1994.
[23]R. H. Voelker and R. J. Lomax, “A finite-difference transmission line matrix method
incorporating a nonlinear device model,” IEEE Trans. Microwave Theory tech. , vol.38,
pp.302-312, March 1990.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
151
[24] S.Goasguen and S. M. El Ghazaly, “Interpolating Wavelet Scheme Toward SelfConsistent Solution o f Semiconductor and Maxwell’s Equations,” Multi-resolution
Workshop Int. Microwave Symp. Boston, June 2000.
[25] S.Goasguen and S. M. El Ghazaly, “Interpolating Wavelet Scheme Toward Global
Modeling o f Microwave Circuits,” IEEE Int. Microwave Symp. Boston, June 2000.
[26] S. Goasguen and S. M. El-Ghazaly, “A Wavelet Based Self-Adaptive Mesh For
Semiconductor Devices Simulations,” European Microwave Conference, Paris, October
2000.
[27] S.Goasguen, M. Tomeh and S. M. El-Ghazaly, “ A Global Modeling Approach
Using Interpolating Wavelets,” IEEE Int. Microwave Symp. Phoenix, May 2002.
[28] T-S. Homg and C-C. Wang, “Microstrip circuit-design using neural-networks,”
IEEE Microwave Theory tech. Symposium Digest, pp.413-416, June 1993.
[29] I. Daubechies, ‘T en lectures on Wavelets,” Society fo r Industrial and applied
Mathematics, 1992
[30] Shirakawa et al., “A large signal characterization o f a HEMT using a multilayered
neural network,” IEEE Trans. Microwave Theory tech., Vol. 45, NO. 9, September 1997.
[31] D. L. Donoho et al., “Data compression and harmonic analysis,” IEEE Transactions
on Information Theory, vol.44, pp.2435-2476, October 1998.
[32]K. C. Gupta, “Emerging trends in millimeter-wave CAD,” IEEE Trans. Microwave
Theory tech ., vol.46, pp.747-755, June 1998.
[33] B. Z. Steinberg and Y. Leviatan, “On The Use o f Wavelet Expansions in the
Method o f Moments,” IEEE Transactions on Antennas and Propagation, Vol. 41 No. 5,
pp. 610-619 May 1993.
[34] S. Goasguen and S. M. El-Ghazaly, “A Coupled FDTD Artificial Neural Network
Technique For Large Signal Analysis o f Microwave Circuits,” International Journal o f
RF and CAD Engineering, special issue on Neural Networks applications to RF.
[35] S.Goasguen and S. M. El Ghazaly, “A Practical Large-Signal Global Modeling
Simulation o f a Microwave Amplifier Using Artificial Neural Network,” IEEE
Microwave and Guided Wave Letters, August 2000.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
152
[36] S. Goasguen and S.M. El-Ghazaly, "New Trends in Global Modeling Techniques:
Artificial Neural Network and Wavelet Based Approach,” National Radio Science
Meeting URSI, Boulder January 2000.
[37] S.Goasguen, S. M. Hammadi and S. M. El Ghazaly, “A Global Modeling Approach
Using Artificial Neural Network,” IEEE MTTsymposium, Anaheim CA June 1999.
[38]A. H. Zaabab, Q-J. Zhang and M. Nakhla, “A neural network modeling approach to
circuit optimization and statistical design,” IEEE Trans. Microwave Theory tech., vol.43,
pp.1349-1358, June 1995.
[39]P. M. Watson and K. C. Gupta, “EM-ANN models for microstrip vias and
interconnects in dataset circuits,” IEEE Trans. Microwave Theory tech ., vol.44, pp.24952503, December 1996.
[40] P. M. Watson and K. C. Gupta, “Design and optimization o f CPW circuits using
EM-ANN models for CPW components,” IEEE Trans. Microwave Theory tech ., vol.45,
NO. 12, December 1997.
[41]A. Veluswami, M. Nakhla and Q-J. Zhuang, “The application o f neural networks to
EM-based simulation and optimization o f interconnects in high speed VLSI circuits,”
IEEE Trans. Microwave Theory tech ., vol.45, pp.712-723, May 1997.
[42]G. L. Creech et al., “Artificial neural networks for fast and accurate EM-CAD o f
microwave circuits,” IEEE Trans. Microwave Theory tech. , vol.45, pp.794-802, May
1997.
[43] J. W. Bandler et al., “Neuromodeling o f microwave circuits exploiting the spacemapping technology,” IEEE Transactions on Microwave Theory and Tech. , vol. 47,
pp.2417-2427, December 1999.
[44] S. G. Mallat, “A Theory for Multiresolution Signal Decomposition: The Wavelet
Representation,” IEEE Transactions on Pattern Analysis And Machine Intelligence, Vol.
11 No. 7, pp. 674-693 July 1989.
[45] M. Krumpholz and P. Russer, “A Field Theoritical Derivation o f TLM,” IEEE
Transactions on Microwave Theory and Techniques, Vol. 42 No. 9, pp. 1660-1668
September 1994.
[46] M. Krumpholz and L. P. B. Katehi, “MRTD: New Time-Domain Schemes Based on
Multiresolution Analysis,” IEEE Transactions on Microwave Theory and Techniques,
Vol. 44 No. 4, pp. 555-571 April 1996.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
153
[47] G.Rubinacci et al., “Interpolating wavelets for the solution o f Maxwell equations in
the time domain,” IEEE Trans, on Magnetics, Vol. 34, NO. 5, September 1998.
[48] M. Fujii and W. J. R. Hoefer, “A three-dimensional Haar-Wavelet-Based
Multiresolution Analysis Similar to the FDTD Method-Derivation and Application,”
IEEE Transactions on Microwave Theory and Techniques, Vol. 46 No. 12, December
1998.
[49] K. Goverdhanam, L. P. B Katehi and A. Cangellaris, “Applications O f
Multiresolution Based FDTD Multigrid,” IEEE Microwave Theory and Techniques
Symposium digest, pp. 333-336 June 1997.
[50] M. Krumpholz and L. P. B. Katehi, “New Prospects For Time Domain Analysis,”
IEEE Microwave and Guided Wave letters, Vol. 5 No. 11, pp. 382-384 November 1995.
[51] K. Sabetfakhri and L. P. B. Katehi, “Analysis o f Integrated Millimeter-Wave and
Submillimeter-Wave Waveguides Using Orthonormal Wavelet Expansions,” IEEE
Transactions on Microwave Theory and Techniques, Vol. 42 No. 12, pp. 2412-2422
December 1994.
[52]M. Righi, M. Mongiardo, R. Sorrentino et al., “Efficient TLM diakoptics for
separable structures,” IEEE Microwave Theory tech. Symposium Digest, pp.425-428,
June 1993.
[53]T-W. Huang, B. Houshmand and T. Itoh, “The implementation o f time-domain
diakoptics in the FD-TD method,” IEEE Trans. Microwave Theory tech., vol.42,
pp.2149-2155, November 1994.
[54] Eswarappa, G. I. Costache and W. J. R. Hoefer, 'Transmission line matrix modeling
o f wide-band absorbing boundaries with time domain diakoptics for S-parameters
extraction,” IEEE Trans. Microwave Theory tech ., vol.38, nu. 4, April 1990.
[55] T.W. Huang, B. Housmand and T. Itoh, “Fast Sequential FDTD Diakoptics Using
the System Identification Technique,” IEEE Microwave and Guided wave letters, pp.
378-380, October 1993.
[56]P. B. Johns and K. Akhtarzad, “The use o f time domain diakoptics in discrete Models
o f fields,” International journal fo r Numerical Methods in Engineering, Vol. 17, pp. 114,1982.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
154
[57] G. Kron, “A set o f Principles to interconnect the solution o f Physical Systems,”
Journal O f Applied Physics, pp.965-980 August 1953.
[58]B. Houshmand, T. W. Huang and T. Itoh, “Microwave structure characterization by a
combination o f FDTD and system identification methods,” IEEE Microwave and Guided
wave letters., vol.3, pp.262-264, August 1993.
[59]W. Kuempel and I. Wolf, “System identification method for transient analysis o f
(M)MIC-components using time iterative methods,” 22nd European Microwave Conf,
vol.l, pp.345-349, August 1992.
[60]W. Kuempel and I. Wolf, “Digital signal processing o f time domain field simulation
results using the system identification method,” IEEE Microwave Theory tech.
Symposium D igest., pp.793-796, June 1992.
[61]D. Haesloop, “A Neural Network Structure for System Identification,” American
Control Conference, pp. 2460-2465, January 1990.
[62]W. L. Ko and R. Mittra, “A combination o f FD-TD and Prony’s methods for
analysing microwave integrated circuits,” IEEE Trans. Microwave Theory tech ., vol.39,
pp.2176-2181, December 1992.
[63] C. Jacoboni and P. Lugli, “The Monte-Carlo Method for Semiconductor Device
Simulation,” Berlin: Springer-Verlag, 1990.
[64] K. Hess, “Monte Carlo Device Simulation: Full Band and Beyond,” Norwell, MA:
Fluwer, 1992.
[65] M. Lundstrom, “Fundamentals o f Carriers Transport,” Addison Wesley 1990.
[66] YJK Feng and A. Hintz, “Simulation o f Submicrometer GaAs MESFET’s Using a
Full Dynamic Transport Model,” IEEE Transactions on Electron Devices, Vol. 35, No. 9,
pp. 1419-1431, September 1988.
[67] F. Heliodore et al., ‘Two-Dimensional Simulation o f Submicrometer GaAs
MESFET’s: Surface Effects and Optimization o f Recessed Gate Structures,” IEEE
Transactions on Electron Devices, Vol. 35, No. 7, pp. 824-830, July 1988.
[68] J.R Zhou and D.K. Ferry, “Simulation o f Ultra-Small GaAs MESFET Using
Quantum Moment Equations,” IEEE Transactions on Electron Devices, Vol. 39, No. 3,
pp. 473-478, March 1992.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
155
[69] JJR. Zhou and DJC. Ferry, “Simulation o f Ultra-Small GaAs MESFET’s Using
Quantum Moment Equations: Velocity Overshoot,” IEEE Transactions on Electron
Devices, Vol. 39, No. 8, pp. 1793-1796, August 1992.
[70] C. M. Snowden, and D Loret, ‘Two-Dimensional Hot Electron Models for Short
Gate-Lengthi GaAs MESFET’s,” IEEE Transactions on Electron Devices, Vol. 34, No.
2, pp. 212-223, August 1987.
[71] S.Goasguen, M. Tomeh and S. M. El-Ghazaly, “ Full Wave Analysis o f FET Fingers
Using Various Semiconductor Physical Models,” IEEE Int. Microwave Symp. Phoenix,
May 2002.
[72] C. M. Snowden, “Introduction to Semiconductor Device Modelling,” World
Scientific, 1986.
[73]K. Blotekjaer, ‘Transport equations for electrons in two-valley semiconductors"/£££
Transactions on Electron Devices, pp 38-47, Vol. 17, January 1970.
[74]C. Hilsum, “Simple empirical relationship between mobility and carrier
concentration,” Electronic Letters, pp 259-260, Vol. 10,1974.
[75] I. Son and T. W. Tang, “Modeling Deep-Level Trap Effects in GaAs MESFET’s,”
IEEE Transactions on Electron Devices, vol. 36, pp.632-640, April 1989.
[76] Q. Li and R. W. Dutton, “Numerical Small-Signal AC Modeling o f Deep-LevelTrap Related Frequency-Dependent Output Conductance and Capacitance for GaAs
MESFET’s on Semi-Insulating Substrates,” IEEE Transactions on Electron Devices, vol.
38, pp. 1285-1288, June 1992.
[77] C. L. Li, T. M. Barton and R. E. Miles, “Avalanche Breakdown and Surface DeepLevel Trap Effects in GaAs MESFET’s ,” IEEE Transactions on Electron Devices, vol.
40, pp.811-816, April 1993.
[78] C.Fiegna et al., “Modeling The Effects o f Traps on the I-V Characteristics o f GaAs
MESFETs,” IEEE IEDM, pp.773-776,1995.
[79] T. M. Barton, and C. M. Snowden, ‘Two-Dimensional Numerical Simulation o f
Trapping Phenomena in the Substrate o f GaAs MESFET’s,” IEEE Transactions on
Electron Devices, vol. 37, pp.1409-1415, June 1990.
[80] W. Batty, A. J. Panks and C. M. Snowden, “Fully Coupled Electro-Thermal
Simulation o f MMICs and MMIC Arrays Based on a Physical Model,” IEEE MTT-S, pp.
693-696,1999.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
156
[81] C. M. Snowden, “Characterization and Modelling o f Microwave and MillimeterWave Large-Signal Device-Circuit Interaction Based on Electro-Thermal Physical
Models,” IEEE MTT-S, pp. 155-158,1997.
[82] J. M. Golio, “Microwave MESFETs and HEMTs,” Artech House, 1982.
[83]X. Zhang and K. K. Mei, “Time-domain finite difference approach to the calculation
o f the frequency-dependent characteristics o f microstrip discontinuities,” IEEE Trans.
Microwave Theory tech ., vol.36, pp.1775-1787, December 1988.
[84]D. M. Sheen et al., “Application o f the three-dimensional finite-difference timedomain method to the analysis o f planar microstrip circuits,” IEEE Trans. Microwave
Theory tech., vol.38, pp.849-857, July 1990.
[85]W. P. Harokopus Jr. and L. P. B. Katehi, ‘‘Electromagnetic coupling and radiation
loss considerations in microstrip (M)MIC design,” IEEE Trans. Microwave Theory tech.,
vol.39, pp.413-420, March 1992.
[86] T. Shibata and H. Kimura, “Computer Aided Engineering for Microwave and
Millimeter-Wave Circuits Using the FD-TD Technique o f Field Simulations,”
International Journal o f Microwave and Millimeter-Wave Computer-Aided Engineering,
Vol. 2, pp. 238-250,1993.
[87] C. H. Dumey et al., “A General Formulation for Connecting Sources and Passive
Lumped-Circuit Elements Across Multiple 3-D FDTD cells,” IEEE Microwave and
Guided Wave Letters vol. 6, pp. 85-87, Feb. 1996.
[88] J. Xu, A. P. Zhao and A. V. Raisanen, “A Stable Algorithm for Modeling Lumped
Circuit Source Across Multiple FDTD cells,” IEEE Microwave and Guided Wave Letters
vol. 7, No. 9, pp. 308-310, September. 1997.
[89] C.L. Wagner and J. B. Schneider, “Divergent Fields, Charge, and Capacitance in
FDTD Simulations,” IEEE Trans. Microwave Theory and tech., vol.46, no. 12, pp.21312136, December 1998.
[90] R. Gillard, S. Dauguet and J. Citeme, “Correction Procedures for the Numerical
Parasitic Elements Associated with Lumped Elements in Global Electromagnetic
Simulators,” IEEE Trans. Microwave Theory and tech., vol.46, no. 9, pp.1298-1306,
September 1998.
[91] C. Christodoulou and M. Georgiopoulos, “Applications o f neural networks in
electromagnetics” Artech House, 2001
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
157
[92] L. P. B. Katehi, J. Harvey and E. Tentzeris, “Time-Domain Analysis Using
Multiresolution Expansions,” in ‘Advances in computational electrodynamics: the FDTD
method’ by Alen Taflove, Artech House, 1998.
[93] M. Tentzeris, J. Harvey and L. P. B. Katehi, “Time adaptive Time-Domain
Techniques for the Design o f Microwave Circuits,” IEEE Microwave And Guided Wave
Letters, Vol. 9 No. 3, pp. 96-99 March 1999.
[94] B. Jawerth and W. Sweldens, “An Overview o f Wavelet Based Multiresolution
Analyses,” SIAM Review, Vol. 36, No. 3, pp. 377-412, September 1994.
[95] W. Sweldens and P. Schroder, “Building Your Own Wavelets At Home,” Tech.
Report IM I 1995:5, Dept o f Mathematics, University o f South Carolina, 1995.
[96] G. Beylkin “On Multiresolution Methods In Numerical Analysis” Documenta
Mathematica, Extra Volume ICM 1998, pp. 481-490.
[97] G; Beylkin, J. Dunn and D. Gines “Order N Static And Quasi-Static Computations
in Electromagnetic Using Wavelets” April 26, 1994
[98] E. M. Tentzeris et al., “Stability and Dispersion Analysis o f Battle-Lemarie-Based
MRTD Schemes,” IEEE Trans. Microwave Theory Tech., vol. 47, pp.1004-1013, July
1999
[99] J. Keiser, “Wavelet Based Approach to Numerical Solution o f Nonlinear Partial
Differential Equations,” Ph-D Thesis, University o f Colorado, Boulder, 1995.
[100] D.L. Donoho, “Interpolating Wavelet Transforms,” Tech. Report 408, Dept, o f
Statistics, Stanford University, November 1992.
[101] W. Sweldens, “The Lifting Scheme: A Construction o f Second Generation
Wavelets,” Lucent Laboratories, Revised November 1997.
[102] S. Grivet-Talocia, “On the accuracy o f Haar-based multiresolution time-domain
schemes,” IEEE Microwave and Guided wave letters, vol. 10, pp.397-399, October 200
[103] M. Holmstrom, “Solving Hyperbolic PDEs Using Interpolating Wavelets,” Tech.
Report, Dept o f Scientific Computing, Uppsala University, Sweden, Dec 1996.
[104] M. Holmstrom, “An Adaptive Finite Difference Method for Time Dependent
PDEs,” Tech. Report, Dept o f Scientific Computing, Uppsala University, Sweden,
August 1997.
Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission.
158
[105] J. Walden, “Filter Banks Methods for Hyperbolic PDEs,” Technical Report, Dept
o f Scientific Computing, Uppsala University, Sweden, Oct 1996.
[106] J. Walden, “Orthonormal Compactly Supported Wavelets for Solving Hyperbolic
PDEs,” Tech. Report, Dept o f Scientific Computing, Uppsala University, Sweden, Dec
1995.
[107] M. Holmstrom and J. Walden, “Adaptive Wavelet Method for Hyperbolic PDEs,”
Tech. Report, Dept of Scientific Computing, Uppsala University, Sweden, Dec 1995.
[108] N. Saito and G. Beylkin “Multiresolution Representations Using The Auto­
correlation Functions O f Compactly Supported Wavelets” Revised January 9, 1992 for
Transactions on Signal Processing.
[109] O. Vasilyev and S. Paolucci, “A Dynamically Adaptive Multilevel Wavelet
Collocation Method for Solving Partial Differential Equations in a Finite Domain,”
Journal o f Computational Physics, 125, pp. 498-512, 1996.
[110] O. Vasilyev, S. Paolucci and M. Sen, “A Multilevel Wavelet Collocation Method
for Solving Partial Differential Equations in a Finite Domain,” Journal o f Computational
Physics, 120, pp. 33-47,1995.
[111] W. Cai and J.Z Wang, “Adaptive Wavelet Collocatioon Methods For Initial Value
Boundary Problems O f Nonlinear PDEs,” ICASE Report, No. 93-48, NASA Langley
Research Center, Hampton. 1993.
[112] R A . Lippert, T.A. Arias and A. Edelman, “Multiscale Computation with
Interpolating Wavelets,” Journal O f Computational Physics 140, pp.278-310,1998.
[113] S. Bertoluzza and G. Naldi “A Wavelet Collocation Method For The Numerical
Solution O f Partial Differential Equations” Applied and Computational Harmonic
Analysis 3,1-9,1996.
[114]M. Toupikov, G. Pan and S. M El-Ghazaly, “Global Modeling o f Microwave
Devices using Wavelets,” IEEE M TT Symposium digest, Baltimore June 1998.
[115] L. Jameson, “On The Differentiation Matrix For Daubechies-Based Wavelets On
An interval,” ICASE Report, No. 93-94, NASA Langley Research Center, Hampton.
1993.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
159
[116] L. Jameson, “A Wavelet-Optimized, Very High Order Adaptive Grid and Order
Numerical Method,” ICASE Report, No. 96-30, NASA Langley Research Center,
Hampton. 1996.
[117] L. Jameson, “Wavelet-Based Grid Generation,” ICASE Report, No. 96-31, NASA
Langley Research Center, Hampton. 1996.
[118] J.S. Hesthaven and L. Jameson, “A Wavelet Optimized Adaptive Multi-Domain
Method,” ICASE Report, No. 97-52, NASA Langley Research Center, Hampton. 1997.
[119] L. Jameson, “On The Daubechies-Based Wavelet Differentiation Matrix,” ICASE
Report, No. 93-95, NASA Langley Research Center, Hampton. 1993.
[120] L. Jameson, “On The Wavelet Optimized Finite Difference Method,” ICASE
Report, No. 94-9, NASA Langley Research Center, Hampton. 1994.
[121] A. Averbuch, G. Beylkin et al. “Multiscale Diversion O f Elliptic Operators” Signal
and Image Representations in Combined Spaces, pp. 1-16,1995.
[122] S. Jaffard “Wavelet Methods For Fast Resolution O f Elliptic Problems” SIAM. J.
Numerical Analysis, Vol. 29, No. 4, pp. 965-986, August 1992.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
BIOGRAPHICAL SKETCH
Sebastien Goasguen was bom on March 22nd 1974 in Rennes, France. He obtained his
Bachelor’s degree in electrical engineering from the Polytechnic Institute o f Toulouse,
France in 1997. Then in 1998, he graduated with honors from King’s College o f London,
England with a Master of Science in electronics research. He arrived at Arizona State
University in 1998. Since then he has been a graduate research associate worldng towards
his Ph-D degree. His area o f research is large and range from linearization techniques and
MMIC design to numerical techniques and semiconductor modeling, hi September 2001,
he will take a research position at Purdue University to work on nanotechnology and
molecular electronics. Sebastien has authored or co-authored 20 international
publications including 5 journal papers.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Документ
Категория
Без категории
Просмотров
0
Размер файла
5 387 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа