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.

1/--страниц